Grégoire G Leroy
email: gregoire.leroy@agroparistech.fr

Lucille Callède
email: lcallede@gmail.com

Etienne Verrier
email: etienne.verrier@agroparistech.fr

J.-C Meriaux

Anne Ricard
email: anne.ricard@toulouse.inra.fr

Coralie Danchin-Burge

Xavier X Rognon
email: xavier.rognon@jouy.inra.fr

Jean- Claude Mériaux

Genetic diversity of a large set of horse breeds raised in France assessed by microsatellite polymorphism

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

During the twentieth century, horse breeding has undergone large changes in Europe. Previously considered as an agricultural, industrial and war tool, horse is now essentially bred for hobby riding. Draught horses, in particular, have been less and less used as utility horses, and many draught breeds have undergone a dramatic decrease in population size: according to the Haras Nationaux, out of the nine French draught breeds, six have annual births below 1000. Measures for in situ conservation have been applied in France for several years but such measures are in general expensive. Therefore, it would be useful to identify priorities among conservation purposes and this requires characterising diversity and genetic relations between breeds [START_REF] Toro | Scientific basis for policies in conservation of farm animal genetic resources[END_REF].

During the last fifteen years, microsatellite markers have frequently been used to evaluate genetic distances and to characterise local breeds, [START_REF] Aberle | Genetic diversity in German draught horse breeds compared with a group of primitive, riding and wild horses by means of microsatellite DNA markers[END_REF][START_REF] Bjørnstad | Genetic structure of Norwegian horse breeds[END_REF][START_REF] Cañon | The genetic structure of Spanish Celtic horse breeds inferred from microsatellite data[END_REF][START_REF] Druml | Individual-based assessment of population structure and admixture in Austrian, Croatian and German draught horses[END_REF][START_REF] Luis | Genetic diversity and relationships of Portuguese and other horse breeds based on protein and microsatellite loci variation[END_REF][START_REF] Marletta | Analysis of genetic diversity and the determination of relationships among western Mediterranean horse breeds using microsatellite markers[END_REF][START_REF] Moureaux | Caractérisation génétique des races françaises de sport et de course et analyse de leur variabilité génétique. 21ème journée de la recherche équine[END_REF][START_REF] Plante | Genetic diversity in a feral horse population from sable island, Canada[END_REF][START_REF] Solis | Genetic diversity within and among four south European native horse breeds based on microsatellite DNA analysis: implications for conservation[END_REF]. Some methods have recently been developed to evaluate the genetic contribution of populations to within-breed and between-breed diversities [START_REF] Caballero | Analysis of genetic diversity for the management of conserved subdivided populations[END_REF][START_REF] Ollivier | Aggregate diversity: New approach combining within-and between-breed genetic diversity[END_REF].

With about 800 000 animals belonging to 50 different breeds (source: Haras Nationaux), France shows a large diversity of horse populations. Among these breeds, 21 have a French origin or have been bred in France for at least a century. According to the FAO, at least 15 populations have disappeared during the last 50 years, and eight indigenous breeds are still considered as endangered or endangered-maintained. Among those breeds, the majority are draught breeds, namely the Ardennais, Auxois, Boulonnais, Poitevin and Trait du Nord breeds, the other ones being the Merens warm-blooded breed and the Landais and Pottock pony breeds. Information on the genetic diversity of French endangered breeds could help breeders and providers, decide where they should place more emphasis.

In the present study, we first analysed the genetic diversity of 39 horse populations reared in France: within-breed diversity, breed relationship and population structure were investigated, using microsatellite data. Then, we focussed on 19 breeds of French origin or having been raised in France for at least a century, and evaluated the conservation priorities between these populations, using different approaches to evaluate within, between and total diversity.

Methods

Populations sampled and microsatellite analysis

French nomenclature divides horse breeds into three groups: warm-blooded, draught horses and ponies. In this study, 39 populations were considered (Table 1). These 39 populations comprised 31 recognised breeds (including 13 warm-blooded breeds, nine draught breeds, and nine pony breeds), the primitive Przewalski horse (used as an outgroup), and seven populations originating from the splitting of two recognised breeds, namely the Anglo-Arab (AA) and Selle Français (SF) breeds (divided into four and three groups, respectively). The 2005 studbook rules define those groups according to the proportion of foreign genes that can be found from genealogical analysis: AA6 and AA9 are considered as pure AA, whereas AA5 and AA10 can have ancestors from another origin, the proportion of Arab origin being higher for AA5 and AA6 than the others. SF8 has a large proportion of PS origin and can therefore be used to produce AA, SFA97 constitutes a group closed to direct foreign influences, whereas SFB98 individuals can have a parent from another breed (under some conditions).

For each of the 39 populations, 23 to 50 animals born between 1996 and 2005, were sampled amounting to 1679 animals. Except for the Przewalski horse, where no pedigree data was available, the sampled animals were known to have no common parents. For the conservation approach, the study focussed on 19 populations, either of French origin, or having been bred in France for at least 100 years (PS, AA and AR breeds). In this approach, 50 animals were randomly sampled among the four and three AA and SF subpopulations, respectively, to constitute two populations.

Eleven microsatellite markers were used to perform the analysis (AHT4, AHT5, ASB2, HMS1, HMS3, HMS6, HMS7, HTG4, HTG6, HTG10, VHL20), with all but two (HMS1 and HTG6) being recommended by the International Society of Animal Genetics for parentage testing and used in similar studies (except HMS1) [START_REF] Marletta | Analysis of genetic diversity and the determination of relationships among western Mediterranean horse breeds using microsatellite markers[END_REF][START_REF] Plante | Genetic diversity in a feral horse population from sable island, Canada[END_REF][START_REF] Solis | Genetic diversity within and among four south European native horse breeds based on microsatellite DNA analysis: implications for conservation[END_REF]. For the entire sample, amplifications and analyses were performed by the same laboratory, using a capillary sequencer (ABI PRISM 3100 Genetic Analyzer, Applied Biosystems).

Statistical analysis

Allele frequencies, mean number of alleles (MNA), observed (Ho) and non-biased expected heterozygosity (He), were calculated using GENETIX [START_REF] Belkhir | GENETIX 4.04, Logiciel sous Windows TM pour la Génétique des Populations. Laboratoire Génome, Populations, Interactions[END_REF]. Wright Fis, Fit and Fst coefficients were also computed using the same software. GENEPOP [START_REF] Raymond | GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism[END_REF] was used to evaluate pairwise genetic differentiation between breeds [START_REF] Goudet | Testing differentiation in diploid populations[END_REF] and departure from Hardy-Weinberg equilibrium, using exact tests and sequential Bonferonni correction [START_REF] Rice | Analysing table of statistical tests[END_REF] on loci. Global tests on Hardy-Weinberg equilibrium were also performed using GENEPOP. Allelic richness was computed using FSTAT [START_REF] Goudet | FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3)[END_REF].

The matrix of Reynolds unweighted distances D R [START_REF] Reynolds | Estimation of the Coancestry Coefficient: Basis for a Short-Term Genetic Distance[END_REF] was computed using POPULATION (Olivier Langella; http:// bioinformatics.org/~tryphon/populations/). Regarding the D R distance, a NeighborNet tree was drawn using SPLITSTREE 4.8 [START_REF] Hudson | Application of phylogenetic networks in evolutionary studies[END_REF]. A factorial correspondence analysis (without the Przewalsky horse) was also performed using GENETIX. Finally, the genetic structure of the populations was assessed using Bayesian clustering methods developed by Pritchard (STRUCTURE, [START_REF] Pritchard | Documentation for STRUCTURE software[END_REF]): using a model with admixture and correlated allele frequencies, we made 20 independent runs for each value of the putative number of sub-populations (K) between 1 and 22, with a burn-in period of 20 000 followed by 100 000 MCMC repetitions. Pairwise similarities (G) between runs were computed using CLUMPP [START_REF] Jakobsson | CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure[END_REF].

To evaluate the conservation priorities in a set of populations, taking into account contributions to within-population and between-population genetic diversity, Ollivier and Foulley [START_REF] Ollivier | Aggregate diversity: New approach combining within-and between-breed genetic diversity[END_REF] have proposed the following method. First, the between-breed contribution (CB) is evaluated, based on the Weitzman [START_REF] Weitzman | On diversity[END_REF] loss Vk of diversity when the population k is removed from the whole set of breeds (in this study we used D R distance). Then, the within-breed contribution (CW) is defined as:

CW = 1 -H(S/k)/H(S) (1)
where H(S) is the average internal heterozygosity of the whole set S and H(S/k) the average internal heterozygosity of the set when k is removed. Finally, the aggregate diversity D of a population is defined as:

D = F st CB + (1 -F st )CW. ( 2 )
The cryopreservation potential (CP) could be computed as the product between the breed contribution (CB) and the probability of extinction (P ex ) of the breed, assumed to be directly proportional to the inbreeding rate (ΔF). Following Simianer et al. [START_REF] Simianer | An approach to the optimal allocation of conservation funds to minimize loss of genetic diversity between livestock breeds[END_REF], P ex can be approximated as

P ex = c ΔF = c/(2Ne) = c (M + F)/8 MF ( 3
)
where Ne is the effective population size, M and F are the numbers of breeding males and females, respectively, used inside the breed in 2005, and c is a constant, to be chosen. Considering that the effective population size of a breed should not be lower than 50 to avoid extinction in the short term [START_REF]FAO: Secondary guidelines for development of national farm animal genetic resources management plans: Management of small populations at risk[END_REF], we considered that P ex = 1 for Ne = 50. Therefore, c was set to 100 (see equation 3).

Caballero and Toro [START_REF] Caballero | Analysis of genetic diversity for the management of conserved subdivided populations[END_REF] have developed a parallel approach. The total diversity GD T can be considered as the exact sum of the gene diversity within population GD WS and the gene diversity between populations GD BS considering the following equations:

GD T = 1 -Σ i Σ j f ij /n 2 (4
)

GD WS = 1 -Σ i f ii /n ( 5
)
GD BS = Σ i Σ j D ij /n 2 (6)
where n is the number of populations, f ij is the average coancestry between populations I and j, and D ij is the Nei minimum distance between populations I and j. The contribution of a population to the diversity is evaluated by computing the loss or gain of diversity ΔGD when the population is removed.

The authors have also proposed to evaluate the contributions (c i ) of each population, which can maximise the total diversity at the next generation, using the following equation:

GD TN = 1 -Σ i c i [f ii -Σ j D ij c j ]. (7)
The contributions can be computed by maximising GD TN in equation ( 7), with the following restrictions: for each population i, c i ≥ 0 and Σ i c i = 1. a W = warm-blooded horse, D = draught horse, P = pony, Pr = primitive horse b France = breeds of French origin or raised in France for at least 100 years; other countries = country of origin for breeds raised in France for less than 100 years c In brackets, number of individuals of each AA and SF subpopulation used when aggregating the four and three subpopulations, respectively

Results

Genetic variations

One hundred and nine alleles were found over all populations and all markers. The average number of alleles per locus was 9.8 ranging from seven (locus HTG4 and HMS1) to 15 (locus ASB2). Some rare alleles in the whole data set were found with a high frequency in the PRW population: for instance, with the HTG6 loci, the two most frequent alleles in the PRW population (70%) were seldom found in other breeds (less than 1%). Heterozygosities, mean number of alleles (MNA) and allelic richness (AR) are presented in Table 2. MNA and AR were highly correlated, (r = 0.98, P < 0.0001). He ranged from 0.43 in the FRI breed to 0.79 in the PFS breed, while Fis per breed ranged from -0.08 (TDN breed) to 0.11 (PRE breed). Some significant heterozygote deficits after corrections were found, for different loci and populations (see Table 2). Only one test exhibited significant excess (AA5 with HMS1). Using global tests, five populations (AB, AR, AUX, CAM, PRE) and two markers (HMS3 and HTG10) showed significant deficit in heterozygotes (P < 0.01). Other studies have shown similar results for these two markers [START_REF] Cañon | The genetic structure of Spanish Celtic horse breeds inferred from microsatellite data[END_REF].

Testing population differentiation, 11 pairs of populations were found non significantly differentiated out of the 741 tests performed: AA5 with AA6, AA9 with AA10, The Fis, Fit, and Fst values were 0.019, 0.116 and 0.099, respectively. We found a gene differentiation coefficient G ST [START_REF] Nei | Analysis of gene diversity in subdivided populations[END_REF] of 0.0989.

Breed relationships and clustering

The NeighborNet network (Figure 1) clearly separated draught horses (also including MER, HAF breeds) and warm-blooded horses, whereas most pony breeds were placed between these two groups. Nordic (IS, SHE, FJ) breeds formed a separate group. FRI and PRW populations were isolated from the other breeds, the closest groups being draught horses and Nordic breeds, for the FRI breed and PRW population, respectively.

Neighbour-Net for the 39 horse populations, based on Reynolds D R distance Figure 1 Neighbour-Net for the 39 horse populations, based on Reynolds D R distance.

In Figure 2, the 38 populations (PRW being excluded) were placed according to the two main axes of the correspondence analysis (accounting for 27.4% and 11.5% of the inertia, respectively). Axis 1 clearly differentiates warm-blooded horses, ponies and draught horses, whereas axis 2 separates Nordic horses (IS, SHE, FJ) from the other ones. The FRI breed seems to be isolated from the other populations, the closest populations being the draught breeds.

Neighbornet and FCA approaches were also used on 34 and 33 breeds, respectively (the four samples of AA breed and three samples of the SF breeds being aggregated into two samples of 50 animals each), showing similar results to previous figures (see Additional files 1 and 2).

Breed assignment to clusters provides complementary information on genetic relationships between populations. As K increases from 2 to 7, mean similarity coefficients among runs are respectively equal to 0.997, 0.993, 0.993, 0.773, 0.562, and 0.658, respectively. Likelihood increased until K reached 15-18 values (see additional file 3), indicating that the most significant subdivisions were obtained for such values. Since mean similarity coefficients were slightly lower for K = 16 (0.78) or 17 (0.81) than for K = 15 (0.83), the results are shown for this last value. Figure 3 shows the assignment of populations to clusters for each K, using runs having the highest pair-wise similarity coefficients. LAND between ponies and AR groups, and APPAL among the three clusters. When K reached 6, depending on the runs, FRI or PRW populations were alternately isolated, which led to a decrease of similarity across runs and explains the low similarity coefficient (0.562) in comparison with other K. When K = 7, these two populations were isolated. The different runs highlight some differences among sub-populations of AA and SF breeds, underlining a more important proportion of AR genes in AA6, AA5 and respectively SFA97 and SF98 groups. Some warm-blooded (FRI until K = 6, MER) and pony breeds (HAF) were classified with draught horses, while the CAM warm-blooded breed was clustered with ponies. As K reached 15, most breeds were shared among different clusters. The ARD, AUX and TDN breed constituted a single cluster while FJ/IS and LUS/PRE constituted two others. In a few cases, a single cluster was essentially associated to a single breed (BOUL, FRI, SHE, PRW).

Partition of diversity

In the set of the 19 French breeds, we found a gene diversity within population GD WS of 0.685, a gene diversity between populations GD BS of 0.073, and a total gene diversity GD T of 0.758. Table 3 shows between-breed, within-breed, and total contribution/variation of diversity according to Ollivier and Foulley [START_REF] Ollivier | Aggregate diversity: New approach combining within-and between-breed genetic diversity[END_REF] and Caballero and Toro [START_REF] Caballero | Analysis of genetic diversity for the management of conserved subdivided populations[END_REF] approaches. For within-breed diversity, CW and ΔGD WS ranged from -0.48 to 0.50 and from -0.0055 to 0.0069 respectively. In both cases, the POIT breed showed a particularly low within-breed diversity. CW and ΔGD WS were negatively correlated (r = -0.715, P = 0.001). For between-breed diversity, CB and ΔGD BS ranged from 0.85 to 12.60 and from -0.0041 to 0.0024, respectively. Here, the POIT breed showed a particularly high contribution to the between-breed diversity. The correlation between CB and ΔGD BS was not significant. D and ΔGD T , accounting for total diversity, were negatively correlated (r = -0.53, P < 0.019). They ranged from -0.32 to 1.25 and from -0.0042 to 0.0039, respectively. In both cases, the ARD and PS breeds showed a particularly low and high diversity, respectively.

Considering contributions to the between-breed diversity and probabilities of extinction, the BOUL, LAND and POIT breeds showed the highest cryopreservation potentials (2.95, 2.95 and 4.83, respectively).

Cluster assignment of each of the 39 populations to the K cluster Figure 3 Cluster assignment of each of the 39 populations to the K cluster. Among 20 runs, solutions having the most similar pair-wise similarity coefficients are presented here. Breeds not classified in their group according to French nomenclature are in italic.

Contributions of each population for an optimal GD T are given in Table 3: the composite PFS breed should contribute to 70% of the pool, for a total GD T of 0.79. Besides, to maximise the total gene diversity, seven of the 19 breeds should be maintained, namely the BOUL, COBND, LAND, PFS, POT, PS and SF breeds.

Discussion

Gene diversity and genetic relations among breeds

Differences between breeds explained 10% of the total genetic variation, which is quite similar to other analyses, where values ranged from 8% to 15% [START_REF] Aberle | Genetic diversity in German draught horse breeds compared with a group of primitive, riding and wild horses by means of microsatellite DNA markers[END_REF][START_REF] Bjørnstad | Genetic structure of Norwegian horse breeds[END_REF][START_REF] Cañon | The genetic structure of Spanish Celtic horse breeds inferred from microsatellite data[END_REF][START_REF] Plante | Genetic diversity in a feral horse population from sable island, Canada[END_REF]. According to previous studies using microsatellites, expected heterozygosities ranged from 0.47 for the FRI breed [START_REF] Luis | Genetic diversity and relationships of Portuguese and other horse breeds based on protein and microsatellite loci variation[END_REF] to 0.80 for the Sicilian Indigenous breed [START_REF] Luis | Genetic diversity and relationships of Portuguese and other horse breeds based on protein and microsatellite loci variation[END_REF]. In our study, only one result was found outside this range of values: 0.43 for the FRI breed, i.e. close to the value found by Luis et al. [START_REF] Luis | Genetic diversity and relationships of Portuguese and other horse breeds based on protein and microsatellite loci variation[END_REF]. Plante et al. [START_REF] Plante | Genetic diversity in a feral horse population from sable island, Canada[END_REF] recently analysed 22 Canadian and Spanish populations. Our estimated values of He were slightly lower (0.71 on average vs. 0.75, P = 0.048) for the eight breeds shared between their study and the present one. Differences on the within-breed diversity among studies using microsatellites can be explained, on the one hand, by the loci used and, on the other hand, by the populations analysed, incidentally belonging to similar breeds but having different recent histories. In the AR breed, we found a He value of (0.72) with a significant deficit of heterozygotes, which can be explained by the fact that this is an international breed in which mating between close relatives is common [START_REF] Moureaux | Genetic variability within French race and riding horse breeds from genealogical data and blood marker polymorphisms[END_REF]. Plante et al. [START_REF] Plante | Genetic diversity in a feral horse population from sable island, Canada[END_REF] and Luis et al. [START_REF] Luis | Genetic diversity and relationships of Portuguese and other horse breeds based on protein and microsatellite loci variation[END_REF] have found similar results for the same breed, but not Aberle et al. [START_REF] Aberle | Genetic diversity in German draught horse breeds compared with a group of primitive, riding and wild horses by means of microsatellite DNA markers[END_REF] who observed a lower heterozygosity (0.57) without a heterozygote deficit. The PER population seemed to have a particularly high genetic diversity in the Plante study (He = 0.78), in comparison with the French PER population (He = 0.68). Because PER populations have been bred in America since the end of the 19 th century, such results should be interpreted bearing in mind that the French PER population has probably suffered from recent bottlenecks due to several modifications of the selection aims.

The three approaches based on genetic relationships (genetic distances, FCA and clustering methods) gave similar results. The populations considered in the present study can be classified into four more or less differentiated clusters: warm-blooded, draught, Nordic and pony breeds. Similar patterns of clustering have been found in other studies [START_REF] Aberle | Genetic diversity in German draught horse breeds compared with a group of primitive, riding and wild horses by means of microsatellite DNA markers[END_REF][START_REF] Bjørnstad | Genetic structure of Norwegian horse breeds[END_REF][START_REF] Plante | Genetic diversity in a feral horse population from sable island, Canada[END_REF][START_REF] Solis | Genetic diversity within and among four south European native horse breeds based on microsatellite DNA analysis: implications for conservation[END_REF]. The draught horses constitute a quite homogenous group, including the nine French CW = contribution to within-breed diversity; CB = contribution to between-breed diversity; D = aggregate diversity;CP = Cryopreservation potential; ΔGD WS = Loss or gain of gene diversity within populations when breed is removed; ΔGD BS = Loss or gain of gene diversity between populations when breed is removed; ΔGD T = Loss or gain of total diversity when the breed is removed; C i = contribution of the breed to optimise GD T draught horse breeds and three breeds presently classified as pony (HAF) or warm-blooded (MER and FRI in a lesser extent) breeds. These three breeds were historically used as draught horse breeds and could therefore have been subject to crossbreeding with other draught horse populations in their past history. Pony breeds formed a group in an intermediate position in comparison to the other clusters. It also included the CAM breed, today recognised as a warm-blooded breed, but morphologically considered as a pony [START_REF] Quittet | Races chevalines en France 2[END_REF]. According to our analysis, FRI and PRW populations were found to be genetically isolated, which can be, to some extent, linked to a low genetic variability [START_REF] Takezaki | Genetic distances and reconstruction of phylogenetic tree from microsatellite DNA[END_REF] due to historical bottlenecks within these breeds [START_REF] Aberle | Genetic diversity in German draught horse breeds compared with a group of primitive, riding and wild horses by means of microsatellite DNA markers[END_REF][START_REF] Ducro | Genetic diversity in the Dutch Friesian horse[END_REF]. Moreover, another parameter explaining isolation of the PRW breed is the presence of rare alleles, which was in agreement with other studies [START_REF] Aberle | Genetic diversity in German draught horse breeds compared with a group of primitive, riding and wild horses by means of microsatellite DNA markers[END_REF] and expected for a population considered as a primitive wild horse.

Population differentiation tests and Bayesian approaches indicate clear differences between sub-populations of AA and SF. Such results may be largely explained by differences in the proportion of thoroughbred (PS) origins in the gene pool of these sub-populations. Within the AA breed, AA5 and AA6 populations appeared distinct from AA9 and AA10 populations and close to the PS breed. This was in agreement with the studbook rules: on the basis of pedigree data, AA5, AA6, AA9 and AA10 populations were indeed found to have respectively 94%, 89%, 44% and 59% of genes from PS origin (Sophie Danvy, personal communication). Within the SF breed, the SF8 (not differentiated from the PS breed) was distinct from SFA97 and SFB98 populations. This result was in agreement with previous results from pedigree data [START_REF] Dubois | Etat des lieux de la sélection Selle-Français. Quelle marge de manoeuvre possible pour d'autres plans de sélection? In 31ème journée de la recherche équine: 2 March 2005 Les Haras Nationaux[END_REF]: the SF8 was found to have 98% of genes from PS origin. The three draught breeds ARD, AUX and TDN, were found to be quite similar, which is linked to a common historical and geographical origin (north of France) [START_REF] Quittet | Races chevalines en France 2[END_REF]. Iberic breeds (LUS and PRE) were also found to be genetically quite close. These results and the fact that according to Bayesian approaches, the likelihood became stable before K reached the number of breeds, indicate that the most relevant division is situated at a level superior to that of the breeds [START_REF] Cañon | Geographical partitioning of goat diversity in Europe and the Middle East[END_REF]. Such a subdivision of the whole set can be explained by the existing crossbreeding management system in several horse populations.

Conservation priorities

In the present study, an almost comprehensive sampling of French breeds was achieved. The different approaches used gave an estimation of the contribution of each breed to the whole French horse stock. Petit [START_REF] Petit | Identifying populations for conservation on the basis of genetic markers[END_REF] has proposed allelic richness as a good parameter to evaluate the genetic diversity of a population, useful as an indicator of past bottlenecks [START_REF] Foulley | Estimating allelic richness and its diversity[END_REF]. In our study, the POIT breed was found to have the lowest allelic richness and also one of the low-est within-breed contributions to diversity according to the two other methods used in the study. Because of the strong correlation with the mean number of alleles, the concept of allelic richness interest seemed to be of limited value in our study.

The results given by the aggregate diversity and gene diversity approaches were slightly correlated. By definition, breeds with low contributions to aggregate and total diversities should have related breeds in the data set. Thus, ARD, TDN, and AUX breeds, which were genetically highly related, illustrate quite well such a hypothesis.

According to the approaches of Ollivier and Foulley [START_REF] Foulley | Estimating allelic richness and its diversity[END_REF] and Cabalero and Toro [START_REF] Caballero | Analysis of genetic diversity for the management of conserved subdivided populations[END_REF], populations that contributed a lot to the total diversity were mostly non-endangered breeds (AR, PS, SF, TF). There were, however, some differences between the two methods when considering the eight breeds classified as endangered or endangered/ maintained by the FAO (ARD, AUX, BOUL, LAND, MER, POIT, POT, TDN). Using the approach of Ollivier and Foulley [START_REF] Foulley | Estimating allelic richness and its diversity[END_REF], contributions to aggregate diversity D of BOUL, MER and POIT breeds were quite high, and taking into account population size, CP was the highest for BOUL, LAND and POIT breeds. Using the approach of Caballero and Toro [START_REF] Caballero | Analysis of genetic diversity for the management of conserved subdivided populations[END_REF], GD T decreased only when LAND and POT breeds were removed, and those two breeds plus the BOUL breed should have been kept to optimise GD T . The differences can be explained by the methods used in the two approaches, particularly considering the evaluation of the contributions to between-diversity. Using the approach of Caballero and Toro [START_REF] Caballero | Analysis of genetic diversity for the management of conserved subdivided populations[END_REF], some Weitzman criteria, such as the twin property [START_REF] Weitzman | On diversity[END_REF], were not applied: for instance, assuming that two populations are genetically identical but very different from the whole set, removing one of them will largely decrease GD BS , which will not be the case when using the Weitzman approach. However, one advantage of the approach of Caballero and Toro [START_REF] Caballero | Analysis of genetic diversity for the management of conserved subdivided populations[END_REF] is the fact that there is no need to give weight to withinand between-diversities to compute total diversity, since by definition GD T is the sum of GD WS and GD BS . In fact, our results outline that both approaches should be considered as complementary to identify which breeds have to be taken into account in a context of genetic resource management. Therefore, conservation priorities should concern particularly BOUL, LAND, MER, POIT and POT breeds.

Another advantage of the method of Caballero and Toro [START_REF] Caballero | Analysis of genetic diversity for the management of conserved subdivided populations[END_REF] is the possibility of computing the contribution of each population to optimise total diversity. Such an approach was designed to conserve a large diversity of alleles. Therefore, it is not surprising to notice that the three breeds (PFS, SF, BOUL) that should have the highest contribution to optimise genetic diversity represent the three identified genetic differentiated groups. The importance of the PFS breed is due to the fact that this synthetic pony breed has the largest number of alleles. SF, another composite breed, has a smaller variability but carries alleles representative of the warm-blooded breed group, while the BOUL breed carries alleles seldom present in the two other breeds but frequent in draught horses.

Finally, several considerations have to be taken into account before taking final conservation decisions [START_REF] Reist-Marti | Weitzman's approach and Conservation of Breed Diversity: an Application to African Cattle Breeds[END_REF], such as the special range of performances for given traits, current production systems associated to the breed, sociocultural value, or dynamics of the group of breeders.

Between 1998 and 2003, births remained more or less stable for BOUL, LAND, POIT and POT breeds, but decreased for the MER breed [START_REF] Danvy | Le point sur les ressources génétiques équines en France[END_REF]. In the endangered breeds, specific uses should be supported to maintain a demand for such horses (production of mules for the POIT breed, ecotourism for local breeds, draught activities, meat production). Genetic variability should also be managed, especially since some of these breeds constitute a pool of original genes (BOUL, MER and POIT) (see Figure 3). For instance, sires with different origins should be used [START_REF] Verrier | La gestion génétique des races d'effectifs limités: principes et application aux cas du cheval de trait Boulonnais et de l'âne grand noir du Berry[END_REF]. When populations of the same breed are raised in other countries (such as the POT breed in Spain [START_REF] Cañon | Geographical partitioning of goat diversity in Europe and the Middle East[END_REF]), regular exchanges should be organised between both countries to maintain a relatively large variety of reproducers.

Conclusion

Based on this study, horse breeds raised in France can be clustered into four groups. These groups were found to be meaningful according to the use of breeds, morphological characteristics and/or geographical origins. The combined use of different methods allowed us to identify breeds for which conservation efforts should be a priority, in order to preserve the maximum genetic variability. Since several horse studies have used similar panels of markers [START_REF] Marletta | Analysis of genetic diversity and the determination of relationships among western Mediterranean horse breeds using microsatellite markers[END_REF][START_REF] Plante | Genetic diversity in a feral horse population from sable island, Canada[END_REF][START_REF] Solis | Genetic diversity within and among four south European native horse breeds based on microsatellite DNA analysis: implications for conservation[END_REF], it would be interesting to merge the corresponding data.

For K = 2 ,

 2 there was a clear separation between draught and warm-blooded horses, with other populations showing intermediate results. When K reached 3, Nordic/primitive breeds, ponies, and some warm-blooded horses segregated more or less clearly from the two other clusters. As K increases to 4 and 5, the five clusters were constituted of Nordic/primitive breeds, draught horses, ponies, warm-blooded populations close to the AR breed and warm-blooded populations close to the PS breed. Some breeds were shared among the last three clusters, such as Correspondence analysis of allele frequencies for 38 of the populations studied (PRW is not included) Figure2Correspondence analysis of allele frequencies for 38 of the populations studied (PRW is not included). The projection is shown on the first two axes.

Table 1 : Basic information on the 39 populations studied Population code Breed Type a Country b Nb of foals registered in 2005 Sample size c

 1

	AA10	Anglo-Arab	W	France	282	(13)
	AA5				781	(11)
	AA6				244	(15)
	AA9				252	(11)
	APPAL	Appaloosa	W	USA	84
	AB	Arab-Barb	W	Morocco	71
	AR	Arab	W	France	1267
	ARD	Ardennais	D	France	645
	AUX	Auxois	D	France	130
	BA	Barb	W	Morocco	99
	BOUL	Boulonnais	D	France	290
	BR	Breton	D	France	3548
	CAM	Camargue	W	France	468
	CO	Connemara Pony	P	Ireland	456
	COBND	Cob Normand	D	France	495
	COMT	Comtois	D	France	4173
	FJ	Fjord	P	Norway	237
	FRI	Friesian	W	The Netherlands	53
	HAF	Haflinger	P	Austria	344
	IS	Iceland Pony	W	Iceland	96
	LAND	Landais	P	France	31
	LUS	Lusitanian	W	Portugal	312
	MER	Merens	W	France	443
	NF	New Forrest Pony	P	UK	119
	PER	Percheron	D	France	1309
	PFS	Poney français de selle	P	France	1069
	POIT	Poitevin	D	France	90
	POT	Pottok	P	France	170

Table 1 : Basic information on the 39 populations studied

 1

(Continued)

Table 2 : Values for parameters of polymorphism within the 39 populations studied

 2

	Population code	He	Ho	F is	HWE deficiency	MNA	AR
	AA10	0.71	0.72	-0.01	0	5.45	5.0
	AA5	0.73	0.71	0.03	0	5.73	5.4
	AA6	0.73	0.71	0.03	0	5.91	5.3
	AA9	0.69	0.70	-0.01	0	4.91	4.6
	APPAL	0.77	0.72	0.06	0	7.55	6.9
	AB	0.76	0.74	0.03	1	7.00	6.7
	AR	0.72	0.66	0.08	1	6.09	5.4
	ARD	0.64	0.62	0.03	0	6.09	5.5
	AUX	0.65	0.62	0.05	1	6.00	5.5
	BA	0.74	0.74	0.00	0	7.00	6.8
	BOUL	0.62	0.60	0.03	1	5.09	4.7
	BR	0.66	0.67	-0.02	0	6.36	5.8
	CAM	0.73	0.68	0.07	1	6.36	6.0
	CO	0.75	0.73	0.03	1	6.64	6.1
	COBND	0.72	0.73	-0.01	0	6.64	6.1
	COMT	0.69	0.67	0.03	2	6.00	5.6
	FJ	0.67	0.69	-0.03	0	6.00	5.6
	FRI	0.43	0.43	0.00	0	3.45	3.2
	HAF	0.65	0.62	0.05	0	4.82	4.6
	IS	0.70	0.68	0.03	1	6.27	5.7
	LAND	0.75	0.71	0.05	1	6.82	6.6
	LUS	0.74	0.71	0.04	1	6.27	5.9
	MER	0.70	0.71	-0.01	0	5.91	5.6
	NF	0.76	0.74	0.03	1	7.64	6.9
	PER	0.68	0.69	-0.01	0	6.64	6.0
	PFS	0.79	0.79	0.00	0	8.09	7.2
	POIT	0.57	0.58	-0.02	0	4.82	4.4
	POT	0.77	0.79	-0.03	0	7.82	7.1
	PRE	0.70	0.62	0.11	1	6.55	5.7
	PRW	0.59	0.56	0.05	0	3.73	3.7
	PS	0.69	0.70	-0.01	0	5.00	4.6
	QH	0.73	0.72	0.01	0	7.00	6.2
	SF8	0.71	0.73	-0.03	0	5.55	4.9
	SFA97	0.74	0.73	0.01	0	6.27	5.7
	SFB98	0.75	0.75	0.00	0	7.00	6.1
	SHE	0.69	0.65	0.06	0	6.00	5.2
	TDN	0.64	0.69	-0.08	0	5.36	5.3
	TF	0.70	0.69	0.01	1	6.27	5.5
	WAB	0.76	0.74	0.03	0	7.55	7.0

He = non biased heterozygosity; Ho = observed heterozygosity; MNA = mean number of alleles; AR = allelic richness; HWE deficiency: number of loci deviating from Hardy-Weinberg equilibrium after Bonferroni correction SF8 and PS, PS and SF8, AA10 with SF8 and PS, AB with BA, APPAL with QH, AUX with TDN, SFA97 with SFB98.

Table 3 : Contributions of the different breeds to genetic diversity according to different approaches Breed code Nb of breeding animals in 2005 Pr. extinction Aggregate diversity and cryopreservation potential (Ollivier and Foulley, 2005) Loss or gain of diversity when a breed is removed and contributions to optimal diversity (Caballero and Toro, 2002)

 3

		Males	Females		CW	CB	D	CP	ΔGD WS	ΔGD BS	ΔGD T	C i
	AA	119	1443	0.11	0.35	0.85	0.39	0.10	-0.0013	-0.0018	-0.0031	0%
	AR	480	2130	0.03	0.29	10.90	1.25	0.35	-0.0015	-0.0010	-0.0026	0%
	ARD	187	1417	0.08	-0.48	1.33	-0.32	0.10	0.0031	0.0001	0.0032	0%
	AUX	24	248	0.57	-0.19	3.14	0.11	1.79	0.0023	-0.0005	0.0018	0%
	BOUL	58	540	0.24	-0.27	12.35	0.87	2.95	0.0040	-0.0023	0.0018	6%
	BR	621	6380	0.02	-0.38	5.57	0.16	0.12	0.0016	0.0009	0.0024	0%
	CAM	118	837	0.12	0.00	7.99	0.73	0.97	-0.0018	0.0013	-0.0006	0%
	COBND	63	760	0.21	-0.06	2.42	0.16	0.52	-0.0017	0.0019	0.0002	2%
	COMT	856	7073	0.02	-0.25	3.63	0.11	0.06	0.0000	0.0015	0.0015	0%
	LAND	22	73	0.74	0.06	3.99	0.41	2.95	-0.0029	0.0016	-0.0014	2%
	MER	93	1012	0.15	-0.04	10.41	0.91	1.53	0.0000	0.0001	0.0001	0%
	PER	183	2461	0.07	-0.32	4.60	0.12	0.34	0.0006	0.0014	0.0020	0%
	PFS	100	949	0.14	0.39	1.93	0.53	0.27	-0.0055	0.0024	-0.0031	70%
	POIT	39	199	0.38	-0.43	12.60	0.75	4.83	0.0069	-0.0030	0.0039	0%
	POT	94	910	0.15	0.19	1.33	0.29	0.20	-0.0040	0.0024	-0.0016	5%
	PS	369	8049	0.04	0.50	6.17	1.02	0.22	-0.0001	-0.0041	-0.0042	1%
	SF	474	11700	0.03	0.45	1.33	0.53	0.04	-0.0024	-0.0013	-0.0037	15%
	TDN	16	183	0.85	-0.17	1.93	0.02	1.64	0.0032	-0.0009	0.0022	0%
	TF	527	15950	0.02	0.36	7.51	1.01	0.18	-0.0002	-0.0029	-0.0032	0%
	Sum				0	100	9.054		0	-0.043	0.043	100%

(page number not for citation purposes)

Acknowledgements

The authors thank the Haras Nationaux for the data provided and Wendy Brand-Williams for linguistic revision.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

JCM carried out the genotyping. AR contributed to the description of the populations and carried out the sampling collection. LC performed the preliminary analysis. GL carried out the computational analysis and prepared the manuscript. XR participated in the computational analysis and preparation of the manuscript. CDB participated in the preparation and the revision of the manuscript. EV participated in the design of the study and the revision of the manuscript. All authors read and approved the final manuscript.

Additional material

Additional file 1

Supplementary Figure 1. Neighbour-Net for the [START_REF] Reist-Marti | Weitzman's approach and Conservation of Breed Diversity: an Application to African Cattle Breeds[END_REF]


Vincent Baillard
email: vincent.baillard@univ-lorraine.fr

Cécile Sulmon

Anne-Kristel Bittebiere

Cendrine Mony

Ivan Couee

Gwenola Gouesbet

Marie Laure Delignette-Muller

Simon Devin

Elise Billoir

Effect of interspecific competition on species sensitivity distribution models: analysis of plant responses to chemical stress

Keywords: Species Sensitivity Distribution, environmental risk assessment, isoproturon, herbicide, biotic interaction, multi-stress

Species Sensitivity Distributions (SSD) are widely used in environmental risk assessment to predict the concentration of a contaminant that is hazardous for 5% of species (HC 5 ). They are based on monospecific bioassays conducted in the laboratory and thus do not directly take into account ecological interactions. This point, among others, is accounted for in environmental risk assessment through an assessment factor (AF) that is applied to compensate for the lack of environmental representativity. In this study, we aimed to assess the effects of interspecific competition on the responses towards isoproturon of plant species representative of a vegetated filter strip community, and to assess its impact on the derived SSD and HC 5 values. To do so, we realized bioassays confronting six herbaceous species to a gradient of isoproturon exposure in presence and absence of a competitor. Several modelling approaches were applied to see how they affected the results, using different critical effect concentrations and investigating different ways to handle multiple endpoints in SSD. At the species level, there was a strong trend toward organisms being more sensitive to isoproturon in

Introduction

Long-term anthropogenic activities disrupt the environment. Due to the complexity of the interactions between living organisms, the accurate evaluation of the effects of these disturbances is delicate. Environmental Risk Assessment (ERA) is used to quantify the risk posed by contaminants and other environmental perturbations to living organisms. It mainly consists in the derivation of Predicted Environmental Concentration (PEC) and their comparison to Predicted No Effect Concentration (PNEC) to assess the risk of chemicals on the environment [START_REF] Amiard | Conventional Risk Assessment of Environmental Contaminants[END_REF]). Assessment factors (AF) are applied to secure the procedure. Such AF take into account the lack of realism of some conditions of classical ecotoxicological experiments and thus aim to extrapolate to realistic environmental conditions by including biotic interactions, intraspecific variability or intra and interlaboratory variability [START_REF] Amiard | Conventional Risk Assessment of Environmental Contaminants[END_REF]). The different ERA methods are divided in 4 tiers [START_REF] Aagaard | Guidance on tiered risk assessment for plant protection products for aquatic organisms in edge-of-field surface waters[END_REF]), the higher tiers using more ecologically relevant data.

Results from methods of higher tiers (i.e. more environmentally relevant methods) are divided by lower AF than those obtained by lower tier methods.

The Species Sensitivity Distribution (SSD; [START_REF] Kooijman | A safety factor for LC50 values allowing for differences in sensitivity among species[END_REF], [START_REF] Posthuma | Species sensitivity distributions in ecotoxicology[END_REF]) is a tier 2 ERA method made to produce PNEC at the community level for a given contaminant using data obtained from monospecific tier 1 bioassays.

The SSD is a distribution of CECs (Critical Effect Concentrations such as Effective Concentrations (EC) or Benchmark Doses (BMD)) obtained for different species for a given contaminant. The species under study have to be chosen to be consistent with the mode of action of the studied contaminant and representative of the considered environment [START_REF] Van Den Brink | Predictive Value of Species Sensitivity Distributions for Effects of Herbicides in Freshwater Ecosystems[END_REF]). The SSD can lead to the derivation of the Hazardous Concentration for five percent of the species (HC 5 ). This value, after application of an AF (ranging from 1 to 5 for SSDs; [START_REF] Amiard | Conventional Risk Assessment of Environmental Contaminants[END_REF]), is commonly used in SSD-based ERA (for example in the United States of America: [START_REF] Stephan | Guidelines for deriving numerical national water quality criteria for the protection of aquatic organisms and their uses[END_REF] and the European Union: European Chemical Bureau [2003]). The SSD method is however based on several assumptions listed in [START_REF] Forbes | Species Sensitivity Distributions Revisited: A Critical Appraisal[END_REF] that still need to be tested.

In particular, very few attempts have been realized to test the effects of biotic interactions on organisms response against contaminants. In the field, organisms are however affected by different types of biotic interactions that can be intra or interspecific. Those interactions can have positive (comensalism, symbiosis or facilitation for example) or negative (competition for example) effects on the organism fitness. Biotic interactions can therefore lead to modifications of the effects of contaminants the organisms are exposed to. For example, [START_REF] Foit | Competition increases toxicant sensitivity and delays the recovery of two interacting populations[END_REF] used a dose-response design to test the effects of intra and interspecific competition on Daphnia magna and Culex pipiens molestus response to fenvalerate. They found that competition tended to increase the toxic effect of fenvalerate on the two tested organisms and modified community dynamics, meaning that toxic exposure might disturb natural community for a longer period than predicted by monospecific bioassays alone. [START_REF] Gust | Daphnia magna's sense of competition: intra-specific interactions (ISI) alter life history strategies and increase metals toxicity[END_REF] investigated the effects of intraspecific competition on Daphnia magna's tolerance to copper and lead using a dose-response design. They concluded that intraspecific competition tended to increase Daphnia magna mortality caused by Cu and Pb. [START_REF] Viaene | Species interactions and chemical stress: Combined effects of intraspecific and interspecific interactions and pyrene on Daphnia magna population dynamics: Combined effect of species interactions and pyrene[END_REF] modeled the effects of intraspecific competition, interspecific competition by B. calyciflorus and predation by Chaoborus sp. (at low and high densities plus control for each type of biotic interaction) on the response of Daphnia magna against pyrene. They found that predation had the highest negative effect on daphnia density under chemical stress and that competition and predation tended to have antagonistic effects on daphnia tolerance to pyrene. One of the hypothesis for this unexpected effect was that pyrene would have modified the community composition, leading to reduced negative biotic interactions. These studies all used macroinvertebrate species and mostly Daphnia magna, and their conclusions need to be extended to other taxonomic groups. For plants, which represent a key compartment of the ecosystems, very few studies have been conducted to test interaction effects between biotic interactions, such as competition, and chemical stress responses. Such a question is however essential for sessile organisms like plants, since they need to be tolerant to all local negative environmental conditions to survive and maintain their population. Results mainly highlighted the existence of interactions between chemical stress responses and competitive interactions between plant species [START_REF] Boutin | Effects of Sublethal Doses of Herbicides on the Competitive Interactions Between 2 Nontarget Plants, Centaurea cyanus L. and Silene noctiflora L[END_REF], [START_REF] Damgaard | The effect of glyphosate on the growth and competitive effect of perennial grass species in semi-natural grasslands[END_REF]), reinforcing the need to take into account such biotic interactions in ERA procedures.

To our knowledge the effects of biotic interactions on SSD results were investigated in only one study by De Laender et al. [2008]. In this study, the effects of ecological interactions on SSD results were assessed by a mechanistic dynamic ecosystem model involving two phytoplancton, three zooplancton and one fish compartments and interspecific interactions between them (predation, competition...). A toxic effect submodel, embedded in the ecosytem model, described the effects of the toxicants on the parameters of the ecosystem model (potentially affected parameters were mortality rate for zooplancton and fish and photosynthesis rate for phytoplancton). One thousand theoretical contaminants were used with randomized 10% effect concentration (EC 10 ) values drawn from a log-normal distribution (SSD without taking ecological interactions into account). Those EC 10 values were used in the mechanistic dynamic ecosystem model, thus simulating the subsequent modifications of the ecosystem fate.

New EC 10 values derived from those simulations were calculated to obtain SSD taking ecological interactions into account. For about 25% of those 1000 contaminants (254) taking into account biotic interactions led to a change in mean (190), standard deviation (94) or both (30) of the SSD compared to the one derived from EC 10 produced without taking into account those interactions.

Moreover, results showed that this trend was higher for herbicides, implying, that the derived HC 5 would have great chances to be different. However, this work needs to be confirmed by data-based studies and to be extended to other types of communities.

Large amount of pesticides are used worldwide to treat a wide diversity of human-made systems. This is the case of croplands to ensure high agricultural productivity, but also of vegetated sport fields and green spaces, to control specific plant composition [START_REF] Alavanja | Pesticides Use and Exposure Extensive Worldwide[END_REF]). Pesticides have effects on the targeted organisms within the treated areas but also on non-targeted organisms outside those areas [START_REF] Crone | Nontarget effects of broadleaf herbicide on a native perennial forb: a demographic framework for assessing and minimizing impacts[END_REF]). Indeed, due to drift, leaching, and runoff, those contaminants are found in terrestrial [START_REF] Geiger | Persistent negative effects of pesticides on biodiversity and biological control potential on European farmland[END_REF]) but also largely in aquatic ecosystems [START_REF] Anderson | Neonicotinoids in the Canadian aquatic environment: A literature review on current use products with a focus on fate, exposure, and biological effects[END_REF], [START_REF] Annett | Impact of glyphosate and glyphosate-based herbicides on the freshwater environment: Impact of glyphosate-based herbicides[END_REF], [START_REF] Mcmahon | Fungicide-induced declines of freshwater biodiversity modify ecosystem functions and services[END_REF]). Isoproturon, for example, is a substituted urea photosynthesis inhibitor, which is often used as a pre-and post-emergence herbicide in cereal crops. Before its banning in the EU in 2017 because of potential groundwater contamination and risks to aquatic life (European Commission [2016]), isoproturon showed a median concentration of 0.02 mg/kg of soil in investigated European agricultural soils [START_REF] Silva | Pesticide residues in European agricultural soils -A hidden reality unfolded[END_REF]). The degradation and fate of isoproturon has also been studied. [START_REF] Johnson | Water movement and isoproturon behaviour in a drained heavy clay soil: 2. Persistence and transport[END_REF] showed that 7 days after isoproturon treatment, concentrations were about 9 mg/kg of soil, after 78 days this concentration was about 1 mg/kg of soil and after 155 days this concentration was above 0.5 mg/kg of soil. More generally, isoproturon was also found to contaminate rivers from one year to another [START_REF] Dragon | The migration of pesticide residues in groundwater at a bank filtration site (Krajkowo well field, Poland)[END_REF]), and to persist in the soil for at least 3 years after application allowing it to be mobilised by increasing rainwater [START_REF] Johnson | Penetration of herbicides to groundwater in an unconfined chalk aquifer following normal soil applications[END_REF]). It is thus clear that a wide diversity of communities, within or outside pesticide-treated areas, are exposed to numerous pesticides, with concentrations ranging from residual to acute levels. But these communities, associating several species, are also the site of numerous biotic interactions, such as competition. This is particularly the case for plant communities, which, as sessile organisms, may thus be submitted simultaneously to both pesticide exposure and potentially strong biotic interactions. Plant communities are therefore a relevant community model to test the assumption: "Interactions between species do not influence the sensitivity distribution" highlighted by [START_REF] Forbes | Species Sensitivity Distributions Revisited: A Critical Appraisal[END_REF] about the SSD method. Considering the studies listed above, we expected that 1) interspecific competition would have a negative effect on organism responses to toxicants in our study and that 2) the protective concentration for 95% of species would be lowered by interspecific competition. To test these assumptions, a dataset of monospecific bioassays results on 6 herbaceous species exposed to isoproturon in presence and absence of competition was produced. We then built SSDs with and without competition and compared them. For this purpose, we tried to characterize the effects of interspecific competition on SSD results using a variety of CECs as well as endpoints and ways to handle them.

Materials and methods

Dataset

Tested species and competitor species

Six herbaceous grass species were chosen as representative of a model vegetated filter strip community. They were selected to represent the natural plant diversity existing in terms of isoproturon tolerance and of competitive ability: Dactylis glomerata, Lolium multiflorum, Arrhenatherum elatius, Trisetum flavescens, Poa pratensis, Poa trivialis. Bromus erectus was chosen as a competitor species to ensure a competitor pressure as constant as possible along the experiment duration. All the seeds were obtained from the Phytosem seed company (Gap, Hautes-Alpes, France).

Experimental design

The pesticide used for exposure was isoproturon, a phenyl-urea photosystem II inhibitor. Despite its recent banning in EU, isoproturon environmental persistence, use in some no-EU countries, and common mode of action, similar to numerous worldwide used herbicides, still make this molecule a relevant pesticide model to study [START_REF] Alberto | Low doses of triazine xenobiotics mobilize ABA and cytokinin regulations in a stress-and low-energy-dependent manner[END_REF], [START_REF] Eker | Determination of The Inhibitory Effect of Isoproturon on Unacclimated Microbial Culture[END_REF], [START_REF] Johnson | Penetration of herbicides to groundwater in an unconfined chalk aquifer following normal soil applications[END_REF]).

Ecotoxicological bioassays were realized for the six tested species and the two competition modalities (absence or presence of the competitor). Experiments were realized in microcosms under controlled conditions (20°C ; 16 hours day of light at 120µmol photons/m 2 /s). Three-liter-round microcosms (20 cm diameter, 12.9cm height) containing inert sterilized vermiculite as substrate were vegetalised by transferring seedlings at 2-leaves phenological stage previously grown from seeds in the absence of chemical stress. For each microcosm, one seedling of the tested species was placed at the center of the enclosure, and for the competition treatment, 37 seedlings of Bromus erectus were planted according to a standardized hexagonal pattern corresponding to 3 circles of, respectively, 6, 12, and 19 seedlings at equal distance [START_REF] Birch | Rectangular and hexagonal grids used for observation, experiment and simulation in ecology[END_REF]). Microcosm contamination with isoproturon began after a 4-days acclimatization.

Plant exposure to isoproturon was carried out by substrate watering in order to induce root chemical exposure. Continuous chemical stress was performed by watering the device with 150 mL of contaminated nutritive solution [Hoagland basal salt mix (No2, Caisson Laboratories, Smithfield, UT, USA) at 0.82 g/L, pH6] twice a week during the 25 days experiment. Five isoproturon (Cluzeau Indo Labo, Sainte Foy la Grande, France) concentrations (0.25, 0.5, 1, 1.5 and 1.75 µM, corresponding to, respectively, 51.5, 103, 206, 309, and 360 µg/L) plus a control were used. For each concentration (6 concentrations including the control) and each competition modality (absence or presence of competitor), eight replicates were realized, leading to 96 microcosms per tested species.

Studied endpoints and metrics

At the end of the experiment, eleven endpoints were measured on the tested species, and an additional metric of dry mass (DM) ratio was calculated. These endpoints were chosen in order to detect the effects of the different treatments (herbicide exposure, competition) applied to the tested species in light of the literature. In addition to global traits (DM traits, shoot height, root length), some traits were preferentially chosen for their responses to herbicide (meximum efficiency of photosystem II determined from Fv/Fm, pigment contents), and to competition (root/shoot DM ratio, Specific Leaf Area (SLA), Leaf Dry Matter Content (LDMC), and ligula height corresponding to stem height). Fv/Fm chlorophyll fluorescence and pigment contents were measured as previously described in [START_REF] Serra | Low environmentally relevant levels of bioactive xenobiotics and associated degradation products cause cryptic perturbations of metabolism and molecular stress responses in Arabidopsis thaliana[END_REF]. SLA and LDMC were quantified as described in [START_REF] Cornelissen | A handbook of protocols for standardised and easy measurement of plant functional traits worldwide[END_REF]. Bromus shoot dry biomass was also weighted, for each microcosm, at the end of the experiment. Analysis of these DM data showed that mean Bromus shoot DM was similar between species and isoproturon treatments, thus allowing to compare all treatment modalities carried out in the experiment.

Data analysis

Critical Effect Concentration values are the elementary components necessary to build SSD. Those CECs were obtained following a stepwise procedure: 1) transformation of data for some endpoints, 2) selection of responsive endpoints for a given species, 3) fitting of concentration-response curves, 4) derivation of CEC values from the fits. The whole modelling process was implemented under the R environment (version 3.5.2; R Core Team [2019]).

Data transformations

Concentration-response curves fitting using non-linear regression assumes a Gaussian error model. Accordingly, the first step was to apply a transformation on response for some endpoints as an attempt to improve homoscedasticity and residuals normality. For each endpoint, we built two ANOVA models where measured values were explained by species, competition modality and isoproturon concentration combinations so that the remaining variability is only inter-replicate variability: one model was built with raw data and the other after transformation of the data.

We then visually inspected the residuals of the models to see if the transformation improved homoscedasticity and residuals normality and in that case, we kept it. Different transformations were tested, depending on the considered endpoint. For the Fv/Fm chlorophyll fluorescence endpoint, which is a propor-tion, we applied a logit transformation (log (p/ (1 -p))) to change its scale, after having previously changed the 0 values (3 out of 471 Fv/Fm measurements) to 0,01 (corresponding to about one fourth of the lowest non-zero value of 0,0395).

For every other endpoint, a log transformation was tested.

Selection of responsive endpoints

Some of the measured endpoints did not present any variation against the isoproturon treatment. Concentration-response modelling was meaningless in such cases and led to numerical issues. A selection step was therefore applied to keep "responsive" endpoints for each species i.e. those that have exhibited variations according to isoproturon concentration. With six species, eleven endpoints and two competition modalities, we had 132 subdatasets. It would have been complicated to inspect visually every subdataset to find the responsive ones. We thus decided to test for their responsivity numerically using the same procedure as in [START_REF] Larras | DRomics: A Turnkey Tool to Support the Use of the Dose-Response Framework for Omics Data in Ecological Risk Assessment[END_REF]. We used a linear trend test to assess the significance of a regression line linking the endpoint values to the isoproturon concentrations. A Benjamini-Hochberg correction was applied on the p-values to reduce false positive selection, due to the high number of realized tests. A 0.05 default threshold was used as false discovery rate. This procedure led to select different endpoints for each species. However, for a given species, an endpoint was selected only if it was responsive with and without competition to enable a comparison of SSDs in both situations.

Concentration-Response Curves modelling

A concentration-response relationship was fitted on data for each selected endpoint, species and competition modality combination. Non-linear regressions were realized with the drc R package (version 3.0-1; [START_REF] Ritz | Dose-Response Analysis Using R[END_REF]) using a log-probit model (see equation 1).

y ij = c + (d -c) * φ (b * (log (x i ) -log (e))) (1)
where i refers to the i th isoproturon concentration and j refers to the j th replicate, x i is the isoproturon concentration, y ij is the response level of the endpoint,

φ

Critical Effect Concentration derivation

For each dose-response curve, two types of CEC were calculated. First, the Effective Concentration which leads to x% of maximum effect (i.e. between the parameters c and d) (EC x ), was calculated. It is the most commonly used CEC in ecotoxicology as it begins to be widely accepted that NOEC (No Observed Effect Concentration) and LOEC (Low Observed Effect Concentration) suffer from some important weaknesses [START_REF] Jager | Bad habits die hard: The NOEC's persistence reflects poorly on ecotoxicology[END_REF]). Although any x value could be used, we have here studied EC 50 , as it is the most widely used value in ecotoxicological studies and it is a direct parameter of the log-probit model we used, and EC 10 , as this value is often used as a no-effect concentration proxy in risk assessment [START_REF] Iwasaki | Does the Choice of NOEC or EC10 Affect the Hazardous Concentration for 5% of the Species?[END_REF]).

Secondly, we calculated the Benchmark Doses (BMD Zsd ) as an alternative to EC x . This CEC, described in the EFSA guidance (EFSA Scientific Committee ), but we also calculated BMD 2sd . This latter would correspond to a change up to one of the bounds of the 95% confidence interval around the predicted value at the control [START_REF] Larras | DRomics: A Turnkey Tool to Support the Use of the Dose-Response Framework for Omics Data in Ecological Risk Assessment[END_REF]).

SSD modelling

SSDs were built via three different scenarios that are commonly used in scientific publications and regulatory texts. A fourth approach, adapted to our multiple endpoint dataset but not used in regulatory texts has also been tested.

1) In a first scenario, the geometric mean of the obtained CECs for the different endpoints has been used for each species. This approach has the advantage to use all the information available for the different endpoints. This is also the most commonly used approach in SSDs in the literature (Xu et al. [2015]). 2)

In a second scenario, we used the lowest CEC value obtained for each species.

This method is very protective, but it is only using a single value, thus making it sensitive to potential outliers. 3) In a third scenario, only the total dry mass endpoint was considered for each species. This was a responsive endpoint for all of the species under study and is usually measured in ecotoxicological studies on herbaceous grass species [START_REF] Del Signore | Development and application of the SSD approach in scientific case studies for ecological risk assessment: Development and application of the SSD[END_REF]). 4) The fourth approach did not directly use CEC values. We first built sensitivity distributions of the different endpoints for each species [START_REF] Arts | Sensitivity of submersed freshwater macrophytes and endpoints in laboratory toxicity tests[END_REF], [START_REF] Hanson | New Technique for Estimating Thresholds of Toxicity in Ecological Risk Assessment[END_REF]) and called them "Endpoint Sensitivity Distribution" (ESD) by analogy to SSDs. We then calculated the fifth percentile of these distributions for each species as a protective concentration for 95% of the endpoints and used these values to build the SSDs themselves. This ESD-using approach was considered to be a good compromise between the two first approaches (i.e. the geometric mean approach and the lowest value approach) as it takes into account all of the available datapoints (like the geometric mean approach) but gives a more protective result (like the lowest value approach) without being too sensitive to outliers.

SSDs have been modeled with log-logistic distributions. The fits were done using the fitdistcens function from the fitdistrplus R package (version 1.0-14;

Delignette-Muller and Dutang [2015]) to integrate censored data and HC 5 values were derived.

Results

Data transformation

After applying the procedure to chose which data to transform, the logit transformation was applied to Fv/Fm chlorophyll fluorescence and the log transformation to every other endpoints except chlorophyll and carotenoid contents.

Selection of responsive endpoints

Table 1 presents the results of the plant species screening and endpoint responsiveness. Five to nine endpoints were selected depending on the species. Some endpoints were selected for all of the species. It is the case for Fv/Fm chlorophyll fluorescence and the dry masses of root, shoot and total. In contrast, some endpoints were not responsive to isoproturon exposure, such as LDMC and pigment contents (Table 1). We noticed that some endpoints for some species were not selected because of the very high variability between replicates (for chlorophyll and carotenoid contents for example). In total, 42 couples of (species, endpoint) were selected, leading to the construction of 84 concentration-response curves (with and without competitor).

Species P. trivialis P. pratensis T. flavescens A. elatius L. multiflorum D. glomerata

Ligula height X X X 8 9 7 5 5 8

Table 1: Responsive endpoints for the different tested species. An "X" indicates that the endpoint is considered responsive for the given species. DM stands for dry mass.

Concentration-response curves modelling

Figure 1 shows an example of curve fitting for P. pratensis for the different selected endpoints. An example of fit for another species can be found in SI (A). No numerical errors were encountered with this fitting procedure and the models have been fitted correctly in every case. The different fits have been visually checked and were consistent with the data. For every endpoint, the concentration-response curve is decreasing with increasing isoproturon concentrations, the only exception being SLA whose value is increasing with the isoproturon concentration.

Critical Effect Concentration derivation

Figure 1 also shows an example of CEC calculation for P. pratensis. The values of the different CECs calculated for the selected endpoints for each species can be found in SI (B). For BMD 2sd , a value could not always be calculated as the BMR Zsd was sometimes beyond the asymptote for high concentrations, meaning that the amplitude of the response was in that case lower than 2 (for BMD 2sd ) times the model's residual standard deviation. This happened 12 times on 84 BMD 2sd calculations and never happened for BMD 1sd calculations even if it was conceptually possible. As the EC x effect level is a percentage of maximum effect (between the two asymptotes), it was successfully calculated in every case.

There were however cases where the calculated EC x or BMD Zsd were above the maximum concentrations. We considered those values as censored values in the interval [maximum tested dose; +∞ [. This happened 3 times among the 84 calculated EC 10 , 20 times among the 84 calculated EC 50 and 10 times among the 84 calculated BMD 2sd . This did not happen for BMD 1sd . Figure 2 shows the values of BMD 1sd calculated for the different species and endpoints (an example for another CEC given in SI (C)). In most cases, the calculated CEC values were lower in presence of a competitor, thus showing that interspecific competition had a negative impact on organisms tolerance on most of the studied species and endpoints.

We calculated for each BMD [32.4, 64.3]).

Figure 3 summarizes the CEC values calculated for the different CEC types, species and endpoints and displays the value with competition against the value without competition. All CEC being pooled, about 10% of points are between the two dashed lines, thus showing that, in this case, competition had here minor or no effect (less than a 5% variation in one way or another). Many more points are below the dashed lines (70%) than above (20%), thus showing that interspecific competition had mostly a negative effect on organism tolerance to chemical stress in our experiment. This is consistent with the hypothesis stated in introduction. This trend toward negative effect can be seen for every calculated CECs and every species and endpoints tested even if it is not systematic.

This trend toward negative effect however seems less visible for EC 50 (44%) than for other CECs (69% for EC 10 , 73% for BMD 1sd and 74% for BMD 2sd ).

ESD and SSD modelling

An example of the fitted ESD for the different species for BMD 1sd is shown in Figure 4. An exemple of figure for another CEC can be found in SI (D).

Figure 5 shows examples of SSDs calculated from BMD 1sd and using the different methods for handling multiple endpoints (minimum, mean, total dry mass and "ESD"). An example of SSD obtained for another CEC is shown in SI (E). Table 2 summarizes the shifts between HC 5 values with and without competition through the ratios between HC 5 values without competition and HC 5 values with competition. We can see in this table that every shift ratio was above one, thus showing that the interspecific competition tended to lower plant tolerance at the community level. We can also see that the competition effect can lower the HC 5 up to 11-fold, which strongly surpasses the Assessment Factor of 5 typically used in SSD-based ERA. This however happened only once and for the least robust method. handling methods. The value before the slash is the HC 5 without competition, the value after the slash the HC 5 with competition and the value between brackets the ratio between the two HC 5 .

EC
Figure 2: BMD 1sd values for the different species and endpoints. A black "X" indicates that the endpoint was not selected for the considered species and that no CEC values was calculated. The grey marks are the CEC values without competition and the black ones are the values with competition. The lines linking those two marks for a given species and endpoint combination are solid lines when competition had a negative effect on organisms tolerance (CEC with competition lower than CEC without competition) and are dashed lines when competition had a positive effect on organisms tolerance (CEC with competition higher than CEC without competition), thus leading to rather facilitation for the species. DM stands for dry mass. graphic displays a way to handle multiple endpoints: "Dry mass" is the SSD produced using the total dry mass for each species, "ESD" is the SSD produced using the 5 th percentile from Endpoint Sensitivity Distribution of each species, "Mean" is the SSD produced using the mean value for each species and "Minimum" is the SSD produced using the minimum value for each species.

Discussion

Ecotoxicological relevance of bioassays

Usual ecotoxicological designs for plants only monitor very few endpoints as root length (ISO/TC 190/SC 4 Caractérisation biologique [2012]) or mass and seedling emergence (OCDE [2006]) and consider they are either representative of the organism health (total dry mass for example) or they are primary maker in response to contamination compared to the other possible endpoints and are thus involved in the plant health [START_REF] Krewski | New Directions in Toxicity Testing[END_REF]). Among the eleven (or twelve if including the DM ratio, cf Table 1) endpoints measured, dry mass endpoints were found responsive to chemical stress for every species, thus reflecting plant health in accordance with the literature. In this study, CECs related to dry masses were often the lowest among the different endpoints for each species (Figure 2), advocating for their use in environmental risk assessment. The Fv/Fm photosystem II efficiency endpoint was also responsive to herbicide exposure. This endpoint, with regards to the mode of action of isoproturon which specifically inhibits photosystem II reaction center [START_REF] Grouselle | Herbicide isoproturon-specific binding in the freshwater macrophyte Elodea densa-a single-cell fluorescence study[END_REF]), is thus a relevant early endpoint to monitor such specific contaminations in plants. In our study, pigment content endpoints exhibited low responsiveness compared to previous works using photosystem II inhibitor herbicides [START_REF] Ramel | Natural variation reveals relationships between pre-stress carbohydrate nutritional status and subsequent responses to xenobiotic and oxidative stress in Arabidopsis thaliana[END_REF], [START_REF] Sulmon | Sugar-induced tolerance to atrazine in Arabidopsis seedlings: interacting effects of atrazine and soluble sugars on psbA mRNA and D1 protein levels[END_REF]). LDMC was also found weakly responsive to chemical exposure, confirming that this parameter was not relevant for plant ecotoxicity tests. Finally, the responsiveness of Root/Shoot DM ratio, SLA, and of length endpoints was found to be species specific underlying the involvement of species specific pattern of organ allocation in response to stress [START_REF] Eziz | Drought effect on plant biomass allocation: A meta-analysis[END_REF], Xiong et al. [2018]).

A similar pattern of responses was observed for some of the endpoints. These different groups of endpoints could be discriminated based on their responsiveness in the different species as well as the ranking of species in terms of tolerance.

A first group contained the ligula height, the maximum shoot height and the root length. They were indeed selected for the same species: D. glomerata, P.

pratensis, P. trivialis (root length was also selected for A. elatius) (Table 1).

These endpoints also shared the same tolerance ranking for concerned species:

D. glomerata > P. pratensis > P. trivialis (Figure 2). A second group gathered the dry masses (root, shoot and total) and concerned all species (Table 1). Again, a common tolerance ranking was observed: L. multiflorum > D.

glomerata > A. elatius > P. pratensis > P. trivialis. The tolerance values of T.

flavescens was however more variable between those different endpoints. Shoot and total dry masses were also the only endpoints where a clear positive effect of competition on tolerance was visible on Figure 2 for L. multiflorum.

Our experimental design allowed to test the effects of competition on the responses of tested species to chemical stress. As Bromus erectus was preidentified as the most tolerant to isoproturon and the most competitive species, it was used as competitor to ensure a constant competition pressure. The competition tested was thus interspecific competition. One could possibly argue that the competition effects identified could not be differentiated from intraspecific competition. Indeed, an interesting experimental design would have been to also test intraspecific competition by exposing the tested species in presence of other individuals of the same species. However, aside from the fact that it would have considerably overburdened the experimental design, it would then have been impossible to interpret correctly the results as the different species display different competitive abilities and isoproturon tolerances. At high isoproturon concentrations, species with low tolerance would have been highly

affected, thus resulting in a reduced competition pressure.

Relevance of SSD approaches in the context of ERA

The number of species used in this study is under the classical SSD standards. [2002] have also proposed that the quality of data is at least as important as its quantity and the fact that the six species studied here come from a consistent community and that several endpoints have been measured for each may be arguments in favor of the quality of our dataset. Finally, the built SSDs are not intended to have a direct regulatory use but to address the effect of interspecific competition on SSD. Even with only six species, our study was able to provide some insight about this question.

Our study used a great number of endpoints. Such approach provides new insights to assess the importance of measured endpoints in ERA through our innovative proposal for handling multiple endpoints. Indeed, by mixing several endpoints at once in our SSD building, we challenged their relevance and their relative individual sensitivity. It seems important to note that the effect of using different endpoints and ways to combine them into SSDs have not been studied yet [START_REF] Del Signore | Development and application of the SSD approach in scientific case studies for ecological risk assessment: Development and application of the SSD[END_REF]). Our results showed that HC 5 calculated with the ESD 5 th percentile was the lowest for 7 of the 8 CEC types and competition modalities (Table 2). Its constructions should also ensure that it is robust towards potential outliers. Regulatory guidances usually propose to use the minimum values or the geometric mean of the different values for each species to build SSDs, the latter being the most commonly used in practice. In light of our work, the 5 th percentile from ESD method we investigated here could be proposed as a decent alternative. This method is however not applicable if the number of endpoints is low (≤ 5). Such a method would also allow to construct more realistic SSD, adjustable to target communities, and suitable to the diversity of contaminants and related mixtures found in ecosystems. Indeed, using a larger set of endpoints enables to include general, species specific, and also contaminant specific endpoints, regarding both community diversity and contaminant mode of action.

Another issue is the question of the "non-responsive" endpoints. We made the choice to exclude them from the analysis, considering that they would not present a response to the isoproturon concentration, even above the maximum concentration of our experimental design. Excluding those data from our dataset could lead to possible biases. Indeed, this exclusion implies that they are not relevant for our study whereas they may be only sparsely or not at all sensitive, thus meaning we may underestimate the total tolerance of tested species and related community.

An interesting continuation of this work would be to set up this experimentation again using other pesticides, especially not banned ones as isoproturon is. This would be of great help to assess if our conclusions on isoproturon could be extended to other compounds or if the responses we observed were specific, notably concerning the selection of "responsive" endpoints. It would also sup-

Conclusion

This study showed that competitive interactions affected, both at the species level (bioassays results) and the community level (SSD results), the responses of plants to isoproturon. Such results thus highlight the relevance of accounting for biotic interactions to construct SSD models in a context of community dynamics prediction and of ERA. SSDs could thus be used more efficiently to design and predict the evolution of key plant communities, such as vegetated filter strips.

They are indeed natural or sown plant communities whose function is to protect aquatic ecosystems by preventing pesticide leaching from crops to surrounding rivers. Installation of these grass strips between croplands and rivers is now regulatory (European Council [1991]). In this context, SSD models could be a relevant tool to predict plant community dynamics under conditions of recurrent pesticide exposures, in order to design and maintain functional buffer grass strips. More data are however necessary to assess these points and research must be carried on those topics.

Associated content

Supporting information. SI A: Example of concentration-response fits for L. multiflorum; SI B:

  is the cumulative probability density of the normal law, e is the concentration at which the maximum slope occurs (equal to the EC 50 ; > 0) b is a shape parameter. if b < 0, c is the response level at high concentrations and d is the response level at low concentrations. if b > 0, c is the response level at low concentrations and d is the response level at high concentrations. If b = 0, the whole model becomes a constant model set at the arithmetical mean between c and d.

  et al.[2017]), has the advantage to take into account data variability. Indeed, the BMD Zsd is the concentration at which the Benchmark Response (BMR Zsd ) is reached, the latter being equal to a change of z times the model residual standard deviation from the control mean. The z of BMD Zsd is therefore theoretically speaking close to the x of EC x as it defines the level considered to have a critical effect. The EFSA guidance proposes to use a z value of 1[START_REF] Efsa Scientific Committee | Update: use of the benchmark dose approach in risk assessment[END_REF]

Figure 1 :

 1 Figure 1: Examples of concentration-response fits for P. pratensis. The points represent the mean of data for each concentration and the curves the fitted models. Data without competition and with competition are highlighted in grey and black respectively. The dashed vertical lines represent the EC 50 levels and the plain vertical lines the BMD 1sd levels. EC 10 and BMD 2sd are not displayed here for reasons of clarity. DM stands for dry mass.

Figure 3 :

 3 Figure 3: Comparison of calculated CEC values in presence of a competitor (y-axis) against the CEC values in absence of a competitor (x-axis) for the four types of calculated CECs.The different species and the different endpoints are respectively described by different point shapes and by their shade of grey. The black solid line represents the first bisector and the two black dashed lines give a 5 percent variation above and below the first bisector. We considered that the points situated between those two dashed lines showed no variation with regard to the competition modality. DM stands for dry mass.

Figure 4 :

 4 Figure 4: Example of Endpoint Sensitivity Distribution fits for BMD 1sd and the six tested species. The stairs represents the empirical cumulative distribution function of the data used to model the ESDs and the curves gives the ESDs themselves. In grey are the data without competition and in black the data with competition. The vertical dashed lines represent the fifth percentile.

Figure 5 :

 5 Figure 5: Example of Species Sensitivity Distribution fits produced with BMD 1sd . The stairs represents the empirical cumulative distribution function of the data used to model the SSDs and the curves give the SSDs themselves. In grey are the data without competition and in black the data with competition. The vertical dashed lines represent the HC 5 levels. Each

  port the potential need to measure more endpoints to assess the plant health and give new insights on the ESD method we propose.Interspecific competition had a negative effect on organism responses to chemical stress, as shown by the different CEC values calculated. It can also be seen in Table2that this negative effect was propagated from single CEC values to HC 5 produced from SSDs. The effect of interspecific competition we evidenced is robust, since HC 5 values always exhibited the same trend, even if its intensity depended on the CEC that was used and on the multiple endpoints handling method. Those results are consistent with the fact that, in SSDs, biotic interactions should have important effects on organism responses to herbicides and other pesticides targeting low trophic levels[START_REF] De Laender | Do we have to incorporate ecological interactions in the sensitivity assessment of ecosystems? An examination of a theoretical assumption underlying species sensitivity distribution models[END_REF]). Competition effect is not taken into account in the first two tiers of ERA as the first consists in using the lowest CEC of classical bioassays and the second in incorporating those bioassay results in predictive models. It is however included in the other tiers as those tiers are not based on monospecific bioassays. The third ERA tier uses results from experiments in mesocosms where different species coexist. Those species are therefore interacting with each other and biotic interactions are occurring. The fourth ERA tier uses data from real ecosystems and is therefore considered as the most environmentally relevant by integrating real interactions and sources of variability. Those two tiers however consider data from specific environments and conclusions can be difficult to expand to other situations.

  Zsd value the corresponding percentage of response (i.e. the x in EC x corresponding to each BMD Zsd value). BMD 1sd corresponded in average to EC 29.7 values (interquartile range of x equivalent for BMD 1sd : [16.5, 38.3]) whereas BMD 2sd corresponded in average to EC 49.0 values (interquartile range of x equivalent for BMD 2sd :

Table 2 :

 2 Summary of HC 5 values (in µM) for the different CEC and multiple endpoints

  Table of all the CEC values calculated; SI C: Graphical summary of EC 50 values; SI D: ESD produced with EC 50 values; SI E: SSD produced with EC 50 values. Figures equivalent to Fig 1 but for other species and figures equivalent to Fig 2, 4 and 5 but for other CEC are available upon request to the authors.

Max shoot height X X X Root length X X X X LDMC X SLA X X X X Fv/Fm X X X X X X Chlorophyll X Carotenoid Root DM X X X X X X Shoot DM X X X X X X Total DM X X X X X XRoot DM/shoot DM X X

Acknowledgements

This work has been funded by the "Ecosphère continentale et côtière" (EC2CO) interdisciplinary program from the Centre National de la Recherche Scientifique (CNRS, France). Vincent BAILLARD received financial support for salary from French Research ministry. We are grateful to Valerie Gouesbet (Ecobio Research Unit, Experimental Ecology platform), to Lucie Wronka and Malo Bourget for help with experimentations, and to Sandra Rigaud (Ecobio Research Unit) for administrative and accounting work. We thank Jean-Pierre Caudal (Ecobio Research Unit, Experimental Ecology platform) and the Experimental Ecology platform of Ecobio research unit for providing all the infrastructure needed to carry out the experiments.

We would like to thank the two anonymous reviewers whose suggestions and comments helped improve this manuscript.

Xiong, G., Zhang, A., Fan, D., Ge, J., Yang, D., Xie, Z., and Zhang, W. (2018).