C Moritz
email: m.moritz@laposte.net

Christine N Meynard

V Devictor

K Guizien

J.-M M Guarini

Nicolas Mouquet

C Labrune

N Meynard

Disentangling the

  L'archive ouverte pluridisciplinaire

The concept of metacommunity, defined as a set of distinct communities connected by dispersal, has marked a turning point in community ecology in shifting the focus from biological interactions to a more complex vision including environmental filtering, dispersal or random effects [START_REF] Leibold | The metacommunity concept: a framework for a multi-scale community ecology[END_REF][START_REF] Holyoak | Metacommunities: a framework for large-scale community ecology[END_REF][START_REF] Logue | Empirical approaches to metacommunities -a review and comparison to theory[END_REF]. While the initial metacommunity frameworks (i.e. neutrality, patch dynamics, species sorting, and mass effect) were intentionally oversimplified, they helped to build a more process-based large-scale ecology [START_REF] Leibold | The metacommunity concept: a framework for a multi-scale community ecology[END_REF][START_REF] Mcgill | Towards a unification of unified theories of biodiversity[END_REF][START_REF] Massol | Linking community and ecosystem dynamics through spatial ecology[END_REF]. For example, the neutral theory, which assumes ecological equivalence between species, has helped defining the role of dispersal and stochastic processes in structuring ecological communities [START_REF] Bell | Neutral macroecology[END_REF][START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF][START_REF] Chave | Neutral theory and community ecology[END_REF]). Similarly, niche-based studies recognising ecological differences between species have emphasized the influence of environmental heterogeneity and patch extinction-colonisation dynamics (Tilman 1994, Mouquet and[START_REF] Mouquet | Coexistence in metacommunities: the regional similarity hypothesis[END_REF]. It is now accepted that communities are structured by a combination of these processes to different degrees [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF][START_REF] Holyoak | Metacommunities: a framework for large-scale community ecology[END_REF][START_REF] Cavender-Bares | The merging of community ecology and phylogenetic biology[END_REF]). The potential of adopting a more integrated approach in which the relative influences of environmental gradients, dispersal, and biological interactions on different types of metacommunities are all quantified is now widely recognised [START_REF] Logue | Empirical approaches to metacommunities -a review and comparison to theory[END_REF]. In other words, under the current metacommunity framework, these mechanisms are not exclusive, but rather complementary [START_REF] Tuomisto | Dispersal, environment and floristic variation of western Amazonian forests[END_REF][START_REF] Cottenie | Integrating environmental and spatial processes in ecological community dynamics[END_REF][START_REF] Ellis | Evaluating the long-term metacommunity dynamics of tree hole mosquitoes[END_REF][START_REF] Meynard | Bird metacommunities in temperate South American forest: vegetation structure, area and climate effects[END_REF][START_REF] Drakou | Freshwater fish community structured more by dispersal limitation than by environmental heterogeneity[END_REF][START_REF] Zhang | Partitioning temperate plant community structure at different scales[END_REF][START_REF] Meynard | Beyond taxonomic diversity patterns: how do a, b and g components of bird functional and phylogenetic diversity respond to environmental gradients across France?[END_REF]. Therefore, rather than trying to define which is the single most important structuring force, the focus is now on understanding their relative importance, or what proportion of the overall variance in community structure can be explained by each process.

Among the patterns that have been used to explain community composition, community similarity ('beta diversity') has been widely documented (Cottenie 2005, Tuomisto and[START_REF] Tuomisto | Analyzing or explaining beta diversity? Understanding the targets of different methods of analysis[END_REF]. In particular, the decrease in community similarity with geographic distance, the so-called 'distance decay', has proven to be common to different groups of organisms [START_REF] Soininen | The distance decay of similarity in ecological communities[END_REF]): closely-located communities are generally more similar in terms of species composition than those located further apart [START_REF] Tuomisto | Dispersal, environment and floristic variation of western Amazonian forests[END_REF]. This pattern is expected if EV-2 dispersal is an important limiting factor in structuring ecological communities. However, spatial autocorrelation in biological data can also arise from environmental filtering: sites that are located at short distances are more likely to be environmentally similar and thus suitable for the same species [START_REF] Legendre | Numerical ecology[END_REF]Legendre 1998, Lichstein et al. 2002). These two different mechanisms (dispersal and environmental filtering) will therefore produce similar predictions regarding distance decay, and are therefore difficult to separate in empirical studies [START_REF] Meynard | Beyond taxonomic diversity patterns: how do a, b and g components of bird functional and phylogenetic diversity respond to environmental gradients across France?[END_REF].

It may nevertheless be possible to disentangle the relative roles of dispersal and environmental filtering in particular natural systems where each one of these produces different spatial structures. In many ecological systems, dispersal is likely to be influenced not only by geographic distance but also by its direction. For example, dispersal of seeds in some terrestrial plants [START_REF] Minami | Various flying modes of wind-dispersal seeds[END_REF] and dispersal of freshwater invertebrates in some aquatic systems [START_REF] Vanschoenwinkel | Wind mediated dispersal of freshwater invertebrates in a rock pool metacommunity: differences in dispersal capacities and modes[END_REF]) are driven by winds, and larval dispersal of marine organisms is influenced by oceanic currents and, to some lesser extent, by larval behaviour [START_REF] Cowen | Population connectivity in marine systems: an overview[END_REF]. In these systems, the directionality of dispersal influences connectivity between communities, which prevents the latter to be fully overlapping with environmental gradients. In such cases, comparing changes in community composition with geographic distance does not account for directional dispersal. Instead, a connectivity measure based on landscape properties (e.g. hydrodynamic, aerodynamic or vector-mediated processes) and the biology of the dispersing organisms should be more relevant to test the relative effects of dispersal and environmental filtering in shaping community structure and composition.

Marine environments are ideal systems in which to implement such analyses. In these systems, connectivity is an important determinant of marine community structure [START_REF] Kinlan | Propagule dispersal and the scales of marine community process[END_REF][START_REF] Watson | Realized and potential larval connectivity in the southern California Bight[END_REF] and can limit population abundances [START_REF] Bode | Larval dispersal reveals regional sources and sinks in the Great Barrier Reef[END_REF][START_REF] Guizien | Dispersal of Owenia fusiformis larvae by wind-driven currents: turbulence, swimming behaviour and mortality in a three-dimensional stochastic model[END_REF]. Moreover, marine organisms are strongly influenced by asymmetric dispersal (i.e. where connectivity between communities A to B differs from connectivity between B to A, [START_REF] Salomon | Effects of asymmetric dispersal on the coexistence of competing species[END_REF], and water flows in the ocean create highly complex structures, generating gyres favouring retention in some areas [START_REF] Hill | The western Irish Sea gyre: a retention system for Norway lobster (Nephrops norvegicus)? -Oceanol[END_REF], or streams favouring long-distance routes [START_REF] Botsford | Principles for the design of marine reserves[END_REF]. In such complex systems, connectivity between communities will thus not be a simple linear, Gaussian or exponential function of geographic distance but, rather, will be influenced by the characteristics (e.g. speed, direction) of the medium transporting organisms [START_REF] Sammarco | Localized dispersal and recruitment in Great Barrier Reef corals: the Helix experiment[END_REF]Andrews 1988, Cowen et al. 2000). Although the importance of local biotic and abiotic conditions in the maintenance of marine organisms has been demonstrated [START_REF] Galeron | Variation in structure and biomass of the benthic communities at three contrasting sites in the tropical northeast Atlantic[END_REF][START_REF] Gogina | Distribution of benthic macrofaunal communities in the western Baltic Sea with regard to near-bottom environmental parameters. 2. Modelling and predictions[END_REF], quantifying the relative effects of realised flow connectivity (effective links between communities and their strengths) versus environmental filtering in the recruitment success of benthic organisms has rarely been achieved [START_REF] Scheltema | On dispersal and planktonic larvae of benthic invertebrates: an eclectic overview and summary of problems[END_REF][START_REF] Roughgarden | Recruitment dynamics in complex life cycles[END_REF][START_REF] Caley | Recruitment and the local dynamics of open marine populations[END_REF][START_REF] Connolly | Theory of marine communities: competition, predation, and recruitment-dependant interaction strength[END_REF].

Here, we account for both flow connectivity (hydrodynamically-driven exchange of individuals between different locations) and geographic distance (linear representation of connectivity) to better understand the spatial structure of a marine invertebrate metacommunity along the coastline of the northwest Mediterranean Sea. The data set used consists of 21 transects, perpendicular to the shore, along the coastline of the Gulf of Lions, where geostrophic currents display a typical directional windinduced regime [START_REF] Millot | The Gulf of Lions' hydrodynamics[END_REF]). We focused on the polychaetes, which are a highly diverse group, both taxonomically and functionally [START_REF] Mchugh | Do life history traits account for diversity of polychaete annelids?[END_REF]. In the Gulf of Lions, polychaetes make up a large part of the abundance and biomass of the benthos and are strongly structured along environmental gradients [START_REF] Grémare | Long-term comparison of macrobenthos within the soft bottoms of the Bay of Banyuls-sur-mer (northwestern Mediterranean Sea)[END_REF][START_REF] Labrune | Assessment of soft-bottom polychaete assemblages in the Gulf of Lions (NW Mediterranean) based on a mesoscale survey[END_REF]. A measure of connectivity that estimates larval dispersal based on ocean current simulations was used as an alternative and more realistic measure of dispersal than that predicted by geographic distance. Note that as we explicitly account for space and the environment, our statistical approach can be related to the variation partitioning strategy outlined and discussed against other approaches in [START_REF] Logue | Empirical approaches to metacommunities -a review and comparison to theory[END_REF]. We tested how much the environmental gradient was structuring community composition and how much of the remaining variation was explained by geographic distance and our measure of flow connectivity. We used the correlation between species turnover (beta diversity) and flow connectivity as an estimation of dispersal limitation in structuring the polychaete metacommunity at different spatial scales. If dispersal is dominant, and spatial autocorrelation is only reflecting dispersal limitations, we would expect beta diversity to be strongly correlated with flow connectivity. On the contrary, if environmental filtering is dominant, and spatial autocorrelation is reflecting environmental structure rather than dispersal limitations, we would expect communities to be strongly correlated with the environment.

Methods

Study area

The Gulf of Lions (Fig. 1) is located in the northwestern part of the Mediterranean Sea, between Cape Creus (Spain) and Cape Croisette (Marseille, France). It is characterised by a large continental shelf, up to 80 km wide at Cape Agde (France), with a mean depth of approximately 76 m. At the boundary of the continental shelf, the shelf break is crossed by submarine canyons. Towards the continent, the local orography plays a key role on the main wind regimes: the Tramontane wind blows from north-northwest between the Pyrenees and the Massif Central (western side of the Gulf ), and the Mistral wind blows from the north in the Rhône valley between the Massif Central and the Alps (eastern side of the Gulf ). Since the Mediterranean Sea has a microtidal regime (tidal amplitudes lower than 0.5 m), the coastal circulation mainly results from the interaction of these winds at regional scale, river inputs, and a large scale thermo-haline northern current flowing from east to west along the continental shelf break [START_REF] Millot | The Gulf of Lions' hydrodynamics[END_REF][START_REF] Estournel | Observation and modeling of the winter coastal oceanic circulation in the Gulf of Lions under wind conditions influenced by continental orography (FETCH experiment)[END_REF].

EV-3

Polychaete survey and environmental data

Between 19 and 29 September 1998, sampling was conducted in the Gulf of Lions along the portion of coast between the French-Spanish border and the mouth of the Rhône River to collect benthic samples from soft bottoms. This section of coast is approximately 110 km from south to north and 140 km from west to east (Fig. 1). Twenty-one cross-shore transects were sampled at two to five stations (five stations on 15 transects, four stations on one transect, three stations on three transects, two stations on two transects, 92 stations in total), with similar sampling protocol. At each station, four sediment samples were collected using a 0.1 m 2 van Veen grab, from which three were taken for benthic macrofauna identification to the species level, and one for granulometric and organic content analysis. Macrofauna samples were immediately sieved with a 1 mm mesh and the fauna retained was fixed in 5% formalin. Samples were then sorted at the laboratory to separate polychaetes from the remaining macrofauna. Identification of polychaete species was done as described in [START_REF] Labrune | Assessment of soft-bottom polychaete assemblages in the Gulf of Lions (NW Mediterranean) based on a mesoscale survey[END_REF]. Polychaetes were identified to the lowest practical taxonomic level and counted. Unidentified taxa were only counted when they could not be mistaken for other identified taxa and counted as separate species. Sediment granulometry analysis showed that the marine sea floor is composed of several sediment types stretching along the coast, ranging cross-shore from sandy to muddy sediments [START_REF] Labrune | Assessment of soft-bottom polychaete assemblages in the Gulf of Lions (NW Mediterranean) based on a mesoscale survey[END_REF]).

Spatial and environmental data analysis

Geographic distance, flow connectivity, and environmental distance between each pair of stations were calculated. Geographic distance between each pair of stations was calculated using latitudinal and longitudinal WGS84 coordinates of each station. Transforming this linear geographic distance matrix to a Gaussian or exponential connectivity matrix was not performed here, contrary to other studies because it would intensify connectivity between close sites and reduce connectivity between distant sites, patterns that are already captured, more accurately, in the flow connectivity matrix. Flow connectivity was assessed using numerical simulations of larval dispersal in the Gulf of Lions including realistic hydrodynamic forcing of 2004. Hydrodynamic simulations (using the SYMPHONIE computer code: [START_REF] Johns | On the wind-driven coastal upwelling in the Gulf of Lions[END_REF][START_REF] Marsaleix | Modélisation tridimensionnelle de la circulation océanique dans le Golfe du Lion[END_REF][START_REF] Estournel | Observation and modeling of the winter coastal oceanic circulation in the Gulf of Lions under wind conditions influenced by continental orography (FETCH experiment)[END_REF] have been validated in a number of studies in this area using in situ and satellite observations (e.g. Rhône river plume, [START_REF] Estournel | The plume of the Rhône: numerical simulation and remote sensing[END_REF][START_REF] Marsaleix | A numerical study of the formation of the Rhône River plume[END_REF][START_REF] Estournel | The Rhône River plume in unsteady conditions: numerical and experimental results[END_REF]northern current inflow, Auclair et al. 2001;wind-induced currents, Estournel et al. 2003). Dispersal of virtual active larvae released from spawning points situated along the coasts of the Gulf of Lions was calculated using the Lagrangian procedure described in [START_REF] Guizien | Dispersal of Owenia fusiformis larvae by wind-driven currents: turbulence, swimming behaviour and mortality in a three-dimensional stochastic model[END_REF]. Flow connectivity among the sampling stations was quantified by counting the number of virtual larvae coming from the spawning points and reaching areas delimited around the stations, between 5 and 50 m depth, after four weeks of dispersal (Supplementary material Appendix 1 and Fig. A1 for details).

An environmental Euclidian distance matrix was calculated between each pair of stations using field environmental measurements directly related to the sediment: depth, grain size (median particle diameter: d50), total carbon and nitrogen content in sediment, organic carbon in sediment, organic carbon to nitrogen ratio in the sediment. We also used other environmental variables, namely chlorophyll a (chl a) and sea surface temperature (SST), that are known to affect water column properties and thus indirectly sediment properties (e.g. through organic matter enrichment). Estimates of these two variables were derived from MODIS imagery ( http://modis.gsfc.nasa.gov/ ) as seasonal averages for 2006. All these environmental variables were scaled to have a mean equal to 0 and a variance equal to 1 prior to distance calculation. EV-4 bands of stations are relatively homogeneous environmentally, particularly with respect to sediment and depth, allowing us to decouple the effects of environment and both connectivity and geographic distance while controlling for the influence of depth.

Results

Both environment and flow connectivity were spatially structured (Fig. 2). When all environmental properties were considered, environmental differences increased with geographic distance, i.e. closer sites were more similar between them than sites located far apart (Fig. 2A). This pattern was reinforced when stations of the same depth were considered (30 m Fig. 2B, and see Fig. A2 for other depths), i.e. environmental similarity among closely located sites was stronger when looking within a depth band. Mantel tests showed that connectivity was not significantly correlated with geographic distance (Table 1, Mantel test). We noted nonetheless a triangular relationship between connectivity and geographic distance (Fig. 2C-D). This reveals that 1) sites far away from each other have weaker chances to be connected, even if they are on the same pathway in terms of hydrodynamics, and that 2) exchange rates of closely located sites can vary from small to large values, likely because of hydrodynamics directionality. Retention (within-community connectivity) reached 4% and occurred for nearly every site during the elapsed time of four week-dispersal, which means that each local community is partly supported by self-recruitment (Fig. A1). Inter-community connectivity rates varied from 0 (absolutely no connectivity) to 3%, and connectivity was greater below the retention diagonal than above, meaning that the flow was predominantly oriented along the coast from east to west (Supplementary material Appendix 1 Fig. A1). This pattern indicates that connectivity is unevenly distributed over the study area. Overall, flow connectivity and geographic distance provided complementary information about spatial structure of the metacommunity.

We found a hump-shaped relationship between species richness and depth (p(depth 2 )  0.001, Fig. 3A; see also [START_REF] Labrune | Assessment of soft-bottom polychaete assemblages in the Gulf of Lions (NW Mediterranean) based on a mesoscale survey[END_REF]), with species richness peaking at 30 m (Fig. 3A). Coarse sediments found in most 10 to 30 m-depth stations also tended to host more individuals (R 2  0.17, p  0.001, Fig. 3B). This abundance effect was driven by one dominant species, Ditrupa arietina, at these stations, but not by a change in the number of species present (Fig. 3C). As expected, depth and grain size were correlated (R 2  0.37, p  0.001, Fig. 3D): coarse (sandy) sediment dominated in shallow water, whereas the mud content of sediments increased with depth. The environmental gradient was therefore stronger when all stations (thus all depths) were included in the analysis, and weaker when each depth was considered separately, except around the mouth of the Rhône River, due to river particle deposits during flooding events [START_REF] Marion | In situ record of sedimentary processes near the Rhône River mouth during winter events (Gulf of Lions, Mediterranean Sea)[END_REF] allowing muddy sediments to be mixed with sandy sediments at shallow depths [START_REF] Labrune | Assessment of soft-bottom polychaete assemblages in the Gulf of Lions (NW Mediterranean) based on a mesoscale survey[END_REF]).

Mantel and partial Mantel tests of the relationships between environmental properties and connectivity and

Estimating beta diversity

Similarity between each pair of polychaete communities was assessed using the Simpson beta-diversity index calculated as b sim  min(b, c)/(min(b, c)  a), where a represents the total number of species occurring in both communities, b represents the total number of species that occur in the neighbouring community but not in the focal one, and c represents the total number of species that occur in the focal community but not in the neighbouring one [START_REF] Koleff | Measuring beta diversity for presence-absence data[END_REF]. This index focuses on species turnover between sites while standardizing the effects of changes in species richness between these sites [START_REF] Koleff | Measuring beta diversity for presence-absence data[END_REF]. It is related to the Jaccard and Sorensen indices, but is less sensitive to sample size, and generally performs better than other beta diversity indices [START_REF] Koleff | Measuring beta diversity for presence-absence data[END_REF][START_REF] Baselga | Partitioning the turnover and nestedness components of beta diversity[END_REF]). We used this index because of large differences in polychaete species richness among the communities sampled, with a few species (e.g. Ditrupa arietina) that largely dominated in some places. This index reaches a maximum when the proportion of species shared in two communities (a) is small, and when the proportion of species gained and lost when moving from the focal (b) to the neighbouring (c) community is similar.

Statistical analysis

We assessed whether and how similarity in species composition among stations was related to environmental gradients, flow connectivity, and geographic distance. Significance of the Pearson product-moment correlations was computed using Mantel tests [START_REF] Mantel | Detection of disease clustering and a generalized regression approach[END_REF]) and partial Mantel tests [START_REF] Smouse | Multiple regression and correlation extensions of the Mantel test of matrix correspondence[END_REF]) using the vegan library [START_REF] Oksanen | vegan: community ecology package[END_REF] in the statistical package R ver. 2.13.0 (R Development Core Team). These tests used a modified regression analysis between two (Mantel test) or three (partial Mantel tests) matrices. The rows and the columns of the matrices were randomly permuted 10 000 times to calculate the correlation coefficient r. Finally, a Monte Carlo simulation method was used to estimate the statistical significance of each test. The Mantel tests were used to assess which component, among the environment, geographic distance and connectivity, was best correlated with beta diversity. The partial Mantel tests were complementary of the Mantel tests to further determine if the effects of one of the components were hidden behind the effects of another. For instance, if the environment is spatially autocorrelated, partialing out geographic distance in a Mantel test using the environment should remove significance of the environment.

Scale of the analysis

Depth can be an important factor determining benthic community structure [START_REF] Garrabou | Structure and dynamics of north-western Mediterranean rocky benthic communities along a depth gradient[END_REF][START_REF] Wei | Global patterns and predictions of seafloor biomass using random forests[END_REF]. All Mantel and partial Mantel tests were thus performed for all stations, and separately for three groups of stations sampled at depths of 1) 10 and 20 m, 2) 30 m, and 3) 40 and 50 m. This division corresponds roughly to sediment layout (bands of similar sediment types stretching along the coastline) and to the four clusters of polychaete assemblages previously defined in [START_REF] Labrune | Assessment of soft-bottom polychaete assemblages in the Gulf of Lions (NW Mediterranean) based on a mesoscale survey[END_REF]. These EV-5 A,B,C) and between the main environmental characteristics (grain size and depth, D). Hump-shaped correlation was found between depth and number of species (A): polynomial model with richness  20.023  depth 2  1.3  depth  7.29, p(depth 2 )  0.001; significant linear correlations were found between grain size and number of individuals (B): R 2  0.17, p  0.001 and depth (D): R 2  0.37, p  0.001; no significant correlation was found between grain size and the number of species (C).

EV-6 geographic distance revealed the complexity of the spatial organisation of these communities. When all stations were considered, beta diversity was significantly correlated with environmental variations (Table 1, Mantel test). This was confirmed by partial Mantel tests which revealed that, when geographic distance or flow connectivity were taken into account, beta diversity was still significantly correlated to the environment (Table 1, partial Mantel tests). Conversely, when the environment was taken into account, beta diversity was neither significantly correlated with geographic distance nor with flow connectivity (Table 1). This suggests that, at least at this scale, the environment was the most important driver of community composition. When the same analysis was performed within depth bands, results showed that the environment, geographic distance and flow connectivity were all significantly correlated with beta diversity (Table 1 and Fig. 4 for 30 m depth; other depth analyses can be found in the Supplementary material Appendix 1 Table A1 and Fig. A3). When the environment was controlled for in the partial Mantel tests (Table 1), both geographic distance and flow connectivity remained significantly correlated with beta diversity. When geographic distance was partialed out, the environment was no longer correlated with beta diversity. On the contrary, when flow connectivity was partialed out, the environment remained significantly correlated to beta diversity. This suggests that the effect of the environmental filtering on species turnover depends on its spatial structure (e.g. the environment is strongly interrelated to geographic distance), but is independent of connectivity.

Discussion

This study provides an example of how the spatial structure of environmental variability can be partially decoupled from connectivity when studying their effects on community composition in dispersive environments. The patterns of spatial structure and connectivity, both when all sites are considered together and when depth bands are separated, reveal an important role of environmental gradients in shaping polychaete communities. Geographic distance per se and flow connectivity appear to be important drivers determining the composition of these communities when patterns of diversity are considered within depth bands, i.e. within more homogeneous environments.

A first important result is that the spatial scale considered influences the most relevant factors determining community structure. We found that alongshore flow connectivity was important when environmental heterogeneity was low, i.e. within depth bands, but that cross-shore environmental gradients were preponderant over connectivity when including all depth bands. This predominant role of the environmental filtering on polychaetes, previously highlighted in empirical studies from [START_REF] Labrune | Assessment of soft-bottom polychaete assemblages in the Gulf of Lions (NW Mediterranean) based on a mesoscale survey[END_REF][START_REF] Labrune | Structure and diversity of shallow soft-bottom benthic macrofauna in the Gulf of Lions (NW Mediterranean)[END_REF], suggests that community composition arises partly from species sorting [START_REF] Holyoak | Metacommunities: a framework for large-scale community ecology[END_REF]. In marine environments, depth and associated characteristics such as sediment layout and grain size are believed to strongly influence the composition of benthic communities [START_REF] Garrabou | Structure and dynamics of north-western Mediterranean rocky benthic communities along a depth gradient[END_REF][START_REF] Labrune | Assessment of soft-bottom polychaete assemblages in the Gulf of Lions (NW Mediterranean) based on a mesoscale survey[END_REF], 2008), setting distinct niches for the organisms, hence driving dissimilarities between communities.

Nevertheless, and consistent with our results, recent theoretical and empirical studies have demonstrated that dispersal is also a significant process whose strength can strongly impact community structure [START_REF] Mouquet | Immigration and local competition in herbaceous plant communities: a three-year seed-sowing experiment[END_REF][START_REF] Cadotte | Dispersal, spatial scale, and species diversity in a hierarchically structured experimental landscape[END_REF][START_REF] Venail | Diversity and productivity peak at intermediate dispersal rate in evolving metacommunities[END_REF]. For instance, the laboratory experimental study carried out by [START_REF] Venail | Diversity and productivity peak at intermediate dispersal rate in evolving metacommunities[END_REF] highlighted that different levels of dispersal could induce changes of diversity in a bacterial metacommunity, and [START_REF] Cadotte | Dispersal, spatial scale, and species diversity in a hierarchically structured experimental landscape[END_REF] found that community structure changed over time and space when dispersal occurred at different spatial scales. In our study area, the decoupled effect of environment and flow connectivity on polychaete community structure may be due to the specific ocean current regime occurring in the area, which is partly uncorrelated with environmental patterns (sediment layout especially). Two nearby sites (i.e. in the same transect perpendicular to the coast, at different depths) may not be strongly connected when the water flow is parallel to the coast [START_REF] Guizien | Dispersal of Owenia fusiformis larvae by wind-driven currents: turbulence, swimming behaviour and mortality in a three-dimensional stochastic model[END_REF], although some of them may be similar in terms of their physical properties (e.g. granulometry, [START_REF] Labrune | Assessment of soft-bottom polychaete assemblages in the Gulf of Lions (NW Mediterranean) based on a mesoscale survey[END_REF]). Exceptions may be encountered in specific areas such as towards the middle of the Gulf, where gyres and cross-shore connectivity may occur, allowing stronger connectivity between sites at different depths in this region.

This decoupling between environmental variability and connectivity is also consistent with empirical results in other systems. For example, [START_REF] Laliberté | Assessing the scale-specific importance of niches and other spatial processes on beta-diversity: a case study from a temperate forest[END_REF] found that the abundance of tree seedlings was first controlled by a Table 1. Correlations of species distance with environment, geographic distance, and connectivity, and correlation of geographic distance with connectivity, for all sites (92) and 30 m-deep sites measured using Mantel and partial Mantel tests. The tests were based on Pearson's product-moment correlation (10 000 permutations) for species distance matrix against matrices of environmental distance, geographic distance, and connectivity. EV-7 also within specific environments to account for structural factors of community composition. To refine the results, analysis including more specific sediment types within depth bands could be conducted, especially in the region of the mouth of the Rhône River which displays a complex sediment layout and where mud is found at shallower depths than in the rest of the Gulf.

Beyond the particular marine system studied, dispersal of propagules also occurs for terrestrial organisms, either through vectors (e.g. vertebrates transporting seeds: Howe and Smallwood 1982) or trough atmospheric (seeds transported by wind: [START_REF] Nathan | Mechanistic models of seed dispersal by wind[END_REF]) and oceanographic processes (seed water-dispersal in coastal areas: [START_REF] Brunbjerg | Species sorting dominates plant metacommunity structure in coastal dunes[END_REF]). On land, some of these processes might be directional, but may probably be more variable due to smaller inertia in air than in water [START_REF] Gaylord | Macroalgal spore dispersal in coastal environments: mechanistic insights revealed by theory and experiment[END_REF] and to the vectors' behaviour (e.g. plants: [START_REF] Ozinga | Dispersal potential in plant communities depends on environmental conditions[END_REF]). In marine ecosystems, dispersal resulting in connectivity between regions can play a major role on community dynamics because it enables the colonisation of distant areas by marine species and favours their mixing at large scales [START_REF] Caley | Recruitment and the local dynamics of open marine populations[END_REF][START_REF] Bode | Larval dispersal reveals regional sources and sinks in the Great Barrier Reef[END_REF][START_REF] Guizien | Dispersal of Owenia fusiformis larvae by wind-driven currents: turbulence, swimming behaviour and mortality in a three-dimensional stochastic model[END_REF]. However, the exact role of connectivity vs environmental gradients or local environments leading to recruitment success remains difficult to assess because of a temporal decoupling of these two factors: larvae are transported by currents and settle when competent, and only after that their chances of developing into adults will depend on the environmental conditions found on their settling site [START_REF] Butman | Larval settlement of soft-sediment invertebrates: the spatial scale of pattern explained by active habitat selection and the emerging role of hydrodynamical processes[END_REF]. Conversely, currents, if too strong, may prevent some organisms from being transported to a suitable habitat or settlement site [START_REF] Guizien | Dispersal of Owenia fusiformis larvae by wind-driven currents: turbulence, swimming behaviour and mortality in a three-dimensional stochastic model[END_REF].

Our study illustrates the importance of another potential caveat concerning connectivity issues in marine environments, which is the lack of accurate estimate of flow connectivity itself. Indeed, empirical studies trying to measure connectivity have only been conducted at small or large spatial scales, for particular species and with specific methods [START_REF] Incze | Neustonic postlarval American lobsters, Homarus americanus, in the western Gulf of Maine: spatial and interannual variations. -Can[END_REF][START_REF] Carson | Evaluating the importance of demographic connectivity in a marine metapopulation[END_REF][START_REF] Saenz-Agudelo | Connectivity dominates larval replenishment in a coastal reef fish metapopulation[END_REF]), but no synthetic connectivity index exists to date. Numerical simulations applied to different areas of the globe have allowed approximations of connectivity patterns for local populations [START_REF] James | The structure of reef fish metapopulations: modelling larval dispersal and retention patterns[END_REF][START_REF] Bode | Larval dispersal reveals regional sources and sinks in the Great Barrier Reef[END_REF][START_REF] Marinone | Connectivity in the northern Gulf of California from particle tracking in a three-dimensional numerical model[END_REF][START_REF] Ayata | How does the connectivity between populations mediate range limits of marine invertebrates? A case study of larval dispersal between the Bay of Biscay and the English Channel (northeast Atlantic)[END_REF]), but constructing accurate hydrodynamic models coupled to particle dispersal is demanding in terms of their design, broad environmental gradient, and then, within this gradient, by other spatial processes such as dispersal occurring at smaller scales. Another theoretical study pointed out that different levels of dispersal and regional environmental heterogeneity combined could affect community structure [START_REF] Mouquet | Consequences of varying regional heterogeneity in source-sink metacommunities: a mechanistic model[END_REF]. Overall, our results provide empirical support for considering that metacommunity structure results from a combination of several mechanisms interacting at different spatial scales [START_REF] Cadotte | Dispersal, spatial scale, and species diversity in a hierarchically structured experimental landscape[END_REF]Fukami 2005, Logue et al. 2011): environmental filtering predominating large scales, and dispersal being more important within similar environments. The possible existence of decoupling between spatial autocorrelation, environmental heterogeneity and connectivity needs to be further integrated using linear (e.g. this study) and non-linear methods to assess whether community structure, spatial structure, and connectivity interact non-linearly [START_REF] Bode | Different dispersal abilities allow reef fish to coexist[END_REF].

Patterns of similarity decay with geographic distance are not sufficient to study the influence of dispersal on metacommunity structure or dynamics since they can be linked to each asset of connectivity, environmental similarity and geographic distance. Decay of similarity is often examined using geographic distance [START_REF] Soininen | The distance decay of similarity in ecological communities[END_REF], as in our analysis, but it has also been found with other distance measures (e.g. environmental distance, [START_REF] Tuomisto | Dispersal, environment and floristic variation of western Amazonian forests[END_REF]stream distance, Maloney and Munguia 2011). In the Gulf of Lions, we observed that environmental distance increased with geographic distance, so we used this environmental distance to perform tests on community similarity matrices, and found a decay of polychaete community similarity with increasing environmental distance. Other studies support this idea [START_REF] Tuomisto | Dispersal, environment and floristic variation of western Amazonian forests[END_REF][START_REF] Duque | Distance decay of tree species similarity in protected areas on Terra Firme forests in Colombian Amazonia[END_REF][START_REF] Maloney | Distance decay of similarity in temperate aquatic communities: effects of environmental transition zones, distance measures and life histories[END_REF]. For instance, [START_REF] Tuomisto | Dispersal, environment and floristic variation of western Amazonian forests[END_REF] found an increase in environmental distance with geographic distance and a decay of similarity both with geographic and environmental distances. A recent study of [START_REF] Maloney | Distance decay of similarity in temperate aquatic communities: effects of environmental transition zones, distance measures and life histories[END_REF] also reported different rates of decay with both environmental and geographic distances, suggesting that global decay patterns can emerge from a combination of different factors. These authors also found different rates of decay for the different regions they considered. Our analyses provide similar results, i.e. decay rates varied when different depth bands were considered. This decay within specific subscales demonstrates that complex community dynamics need to be investigated not only at large regional scale but EV-8 (e.g. seeds for plants). Our results may also be useful for diversity conservation since the design of marine protected areas and terrestrial natural reserves requires a sound understanding of community dynamics, its dispersal patterns and its interaction with habitats to protect, at both local and regional scales [START_REF] Economo | Biodiversity conservation in metacommunity networks: linking pattern and persistence[END_REF]. New analyses will be needed to meet these requirements, including elements of mechanistic landscape ecology and statistical analysis of biogeographic patterns within the same framework.

parameter isation, and computer resources needed to run them, especially when designed at fine resolution or close to the coast where fine-scale physical processes occur. Here, connectivity was determined using a realistic 3D hydrodynamical model with a simulation resolution and coverage adapted to study dispersal at large scale [START_REF] Guizien | Dispersal of Owenia fusiformis larvae by wind-driven currents: turbulence, swimming behaviour and mortality in a three-dimensional stochastic model[END_REF]. Because the coupling with a species-specific larval behaviour may not account for the different behaviours and pelagic larval duration (PLD) of all the polychaete species present in the Gulf, we calculated exchange rates using a four-week PLD, which is a good estimate of mean PLD for benthic invertebrates [START_REF] Mcedward | Ecology of marine invertebrate larvae[END_REF]. Subsequently incorporating different types of larvae in the model (active, passive, and with different PLD) was not possible in our simulations, but attempts to do so were achieved in [START_REF] Guizien | Using realistic Lagrangian larval dispersal simulations for marine protected area design: application to the Gulf of Lions (northwest Mediterranean)[END_REF]. The difficulty to account for these discrepancies lies in calculating a final connectivity matrix that incorporates appropriate proportions of all larval species, which first implies running separate simulations, species by species. Only this could lead to more realistic estimations of flow connectivity, and consequently, of connectivity effects on metacommunity structure. This could at the same time allow assessing the importance of connectivity on particular species of interest [START_REF] Treml | Reproductive output and duration of the pelagic larval stage determine seascape-wide connectivity of marine populations[END_REF]. Another limitation might come from having used a flow connectivity matrix estimated for 2004 and a one-year regional diversity survey in 1998: coastal circulation and diversity temporal variability may partly blur the correlation between the average connectivity matrix and the diversity survey, even if inter-annual variability of the average connectivity matrix can be ignored. Moreover, the single year diversity survey may have been partly biased by peculiar connectivity in the years before diversity assessment for short-lived species. It would therefore be beneficial to gather connectivity matrices for different years, especially around the year of diversity surveys.

Metacommunities are structured by a number of processes, in particular by environmental and dispersal filters [START_REF] Leibold | The metacommunity concept: a framework for a multi-scale community ecology[END_REF][START_REF] Holyoak | Metacommunities: a framework for large-scale community ecology[END_REF], which ecologists have mostly studied as separate paradigms [START_REF] Logue | Empirical approaches to metacommunities -a review and comparison to theory[END_REF]. Here, by using a measure of connectivity decoupled from geographic distance, we have shown than these paradigms may be spatially interlocked, leading to a more comprehensive understanding of metacommunity dynamics in this region and allowing to sort out the importance and scale of action of environmental and dispersal processes. Incorporation of connectivity estimates in the analysis, a characteristic that few studies have done to date, rendered the study more realistic concerning the processes occurring in dispersive environments. Environmental properties that were used in this study are known to influence polychaete communities [START_REF] Arvanitidis | Seascape biodiversity patterns along the Mediterranean and the Black Sea: lessons from the biogeography of benthic polychaetes[END_REF][START_REF] Labrune | Assessment of soft-bottom polychaete assemblages in the Gulf of Lions (NW Mediterranean) based on a mesoscale survey[END_REF]), but other variables may be taken into account in further analysis (e.g. oxygen, phosphate and silicate). Further work and other methods should be implemented to capture more precisely the effects of dispersal only on community structure, and assess in which type of marine ecosystem dispersal is a critical process. The approach presented here could also be applied to other freshwater and terrestrial ecosystems in which connectivity is known, and displaying complex or unusual dispersal patterns (e.g. asymmetric, oriented, gyres), and to other organisms dispersing particles

Figure 1 .

 1 Figure 1. Location of the 92 sites sampled during the oceanographic campaign in 1998. Transects were located along the Gulf of Lions, between the French-Spanish border and the eastern part of the Rhône River. Stations were located between 10 and 50 m water depth. Adapted from Labrune et al. (2007).

Figure 2 .

 2 Figure 2. Relationships between geographic distance and environmental distance (A): R 2  0.11, p  0.001; (B): R 2  0.30, p  0.001 and flow connectivity measured as exchange rates (C, D) for all stations and for stations located at 30 m-depth. The shading of each square is proportional to the density of dots in that square.

Figure 3 .

 3 Figure 3. Relationships between environmental and community characteristics (number of individuals and number of species,A, B, C) and between the main environmental characteristics (grain size and depth, D). Hump-shaped correlation was found between depth and number of species (A): polynomial model with richness  20.023  depth 2  1.3  depth  7.29, p(depth 2 )  0.001; significant linear correlations were found between grain size and number of individuals (B): R 2  0.17, p  0.001 and depth (D): R 2  0.37, p  0.001; no significant correlation was found between grain size and the number of species (C).

Figure 4 .

 4 Figure 4. Relationships between beta diversity and flow connectivity, geographic and environmental distances for 30 m-depth sites. Logarithm added of one for exchange rate (in percent), geographic distance (in km) and Euclidian environmental distance is given.

Acknowledgements -This work was partly supported by the French Ministère de l'Enseignement supérieur et de la Recherche through a doctoral grant. VD received funding from the Fondation pour la Recherche sur la Biodiversité (FRB, research projects FABIO and PHYBIO). We thank Antoine Grémare (EPOC-OASU, Arcachon, France), François Charles and Jean-Michel Amouroux (LECOB, Banyuls-sur-mer, France), Joao Gil and Rafael Sarda (CEAB, Blanes, Spain) who have organised and realised data collection, identification and database design. We are grateful to members of POC (Toulouse, France) for providing the ocean model SYMPHONIE. Finally, we thank Julian Caley and Camille Mellin for helpful advice and comments on the manuscript.


Camille Meslin

Sylvie Mugnier

Isabelle Callebaut

Michel Laurin

Géraldine Pascal

Anne Poupon

Ghylène Goudet-Guitton

Philippe Monget
email: philippe.monget@tours.inra.fr

Ge ´raldine Pascal

Ghyle ` Ne Goudet

Evolution of genes involved in gamete interaction: Evidence for positive selection, duplications and losses in vertebrates

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

A lot of barriers are able to prevent a successful reproduction. Dobzhansky proposed a classification of species isolation from prezygotic to postzygotic isolation [5]. In postzygotic isolation, a zygote is formed but the hybrid offspring is not viable or sterile. In prezygotic isolation, species could not mate because of a difference in sexual behavior or for mechanical issues. Another step of prezygotic isolation is the lack of interaction between male and female gametes, the sperm-egg interaction being a several steps process that leads to a successful fertilization. Although details vary between species, a successful fertilization requires three main steps in sperm-egg interaction [START_REF] Yanagimachi | Mammalian fertilization[END_REF]. The first step is the recognition and the binding between sperm and the extracellular matrix surrounding eggs (ECM), called the chorion in teleosts, the vitelline envelope in birds, and the zona pellucida (ZP) in mammals. The second step, required for the sperm to penetrate through the egg coat, is the acrosome reaction. The third step is the binding and the subsequent fusion between sperm and the oolemma, the plasma membrane of the oocyte. Many experiments have been performed in order to determine the importance of candidate proteins during fertilization. Surprisingly, knock-out experiment of B4GT1 or targeted mutation that abolished the protease activity of ACR in mouse show that these proteins are important but not essential [START_REF] Baba | Sperm from mice carrying a targeted mutation of the acrosin gene can penetrate the oocyte zona pellucida and effect fertilization[END_REF][START_REF] Lu | Sperm from beta 1,4-galactosyltransferase-null mice are refractory to ZP3-induced acrosome reactions and penetrate the zona pellucida poorly[END_REF]. Several studies showed that genes involved in sperm-egg interaction exhibit a particular evolution. For example, a recent study shows that several ZP genes have been lost during evolution of vertebrates [START_REF] Goudet | Phylogenetic analysis and identification of pseudogenes reveal a progressive loss of zona pellucida genes during evolution of vertebrates[END_REF]. ZP1 gene, one of the glycoproteins that compose the ZP, has been lost during evolution in cattle, pig, cat and dog, whereas ZP4 gene has been lost in mouse. Moreover, ZPAX and ZPD genes, which are present in chicken, have been lost in mammals. In addition, several genes involved in fertilization such as ZP3, ADAM1, ADAM2, ACR and CD9 have been shown to evolve under positive selection [START_REF] Swanson | Pervasive adaptive evolution in mammalian fertilization proteins[END_REF][START_REF] Swanson | Positive Darwinian selection drives the evolution of several female reproductive proteins in mammals[END_REF]. The authors hypothesized that these positively selected amino acids might play a role in the specificity of sperm-egg interaction between species because they are located in the putative sperm-egg binding domain. Despite the fact that several studies suggested the involvement of positive Darwinian selection on sperm-egg interaction genes, none has been done exhaustively. In our study, we have listed all genes which have been experimentally shown to be involved in one vertebrate and in at least one step of the sperm-egg interaction, on both sperm and/ or oocyte, and on 19 species with fully-sequenced genomes. We have studied all these genes with an evolutionary perspective by searching for gene gains and losses, as well as positive selection, a divergent evolution of these genes among vertebrate species potentially leading to reproductive isolation.

Results and Discussion

Rates of Gene Gain, Loss and Pseudogenization are Clade-dependent Variables

We identified 69 genes from the literature (Figure 1 and Table S1) that encode proteins experimentally known to be involved in at least one of the three steps described above. These genes are well studied in mouse and cattle that are common models for studies on physiology and reproduction. We classified all genes into three main groups according to which step they are involved in, excepted CRISP1 that is involved in all three steps of fertilization [START_REF] Busso | A novel function for CRISP1 in rodent fertilization: involvement in sperm-zona pellucida interaction-Evidence for the involvement of testicular protein CRISP2 in mouse sperm-egg fusion[END_REF][START_REF] Cohen | Participation of cysteine-rich secretory proteins (CRISP) in mammalian spermegg interaction[END_REF][START_REF] Da Ros | Impaired sperm fertilizing ability in mice lacking Cysteine-RIch Secretory Protein 1 (CRISP1)[END_REF]. Thanks to phylogenetic analyses, we found several gene duplications, as well as pseudogenizations, among genes involved in the three main steps of fertilization. In particular, we found evidence for (i) gain of FUT5 and Zp3r genes in primates and rodents, respectively, (ii) pseudogenization of SPA17 and ADAM15 genes in chicken, ZAN, ADAM1a and ADAM3 in human, ADAM3 in chimpanzee, spermadhesins in dog, IGSF8 and SPESP1 in horse, as well as ZPAX in dog, and the known ZP pseudogenes published elsewhere [START_REF] Goudet | Phylogenetic analysis and identification of pseudogenes reveal a progressive loss of zona pellucida genes during evolution of vertebrates[END_REF] (ZPB/ZP4 in mouse, ZP1 in cow and dog, ZPAX in primates and cow); and (iii) specific duplication events of ADAM genes in therian mammals, spermadhesins in ungulates, as well as the ACR gene in Xenopus (14 copies) and HSPE1 in macaque (11 copies). For each gene, we have systematically assessed, by optimizing gene presence, number, and type (functional or pseudogene) on a time-calibrated reference phylogeny, the absolute rate of the following four types of events: gene gain, duplication, inactivation (resulting in a pseudogene) and complete loss (presumably through deletion, although incomplete genome annotation may have led to overestimation of loss rate). The overall rate of gene appearance (genes that appear only in some monophyletic groups), pseudogenization and deletion events differs strongly between clades and type of event (Figure 2a). By taking into consideration the number and the geological age of sampled lineages, we found that these three types of events are far more frequent in eutherians than in birds and teleosts. A phylogeny-informed binomial test shows that these differences between eutherians and teleosts are statistically significant, differences between eutherians and birds are also statistically significant for deletions and gene appearance. However, a possible bias stems from the fact that our list of gene is composed almost exclusively of genes identified in eutherians and thus, the rates of the different events may be overestimated in this clade. No pseudogene has been identified in teleosts. The duplications occur at more comparable rates across these three taxa but tend to occur slightly more frequently in eutherians than in teleosts and aves. The consensual conclusion is that duplication rates for most of the genes are higher in teleosts than in other taxa due to the third whole genome duplication, specific to teleosts [START_REF] Jaillon | Genome duplication in the teleost fish Tetraodon nigroviridis reveals the early vertebrate proto-karyotype[END_REF][START_REF] Sato | Teleost fish with specific genome duplication as unique models of vertebrate evolution[END_REF]. Consequently, our results are surprising. They reflect the greater geological age of lineages of teleosts, resulting in a phylogenetic diversity twice as great as eutherians, despite a lower number of included species (Figure 2a). Thus, a greater number of events per lineage is required in teleosts to yield comparable rates.

Genes of Teleosts are More Subjected to Positive Selection

We have also submitted the protein sequence to the Phylea-sProg web server [START_REF] Busset | PhyleasProg: a user-oriented web server for wide evolutionary analyses[END_REF] (http://phyleasprog.inra.fr/) to determine selective pressures that shape the evolution of these genes and used a new phylogenetically-informed test to compare the amount of positive selection in various clades. The calculation of positive selection is globally based on rate of non-synonymous (dN) and synonymous (dS) substitutions [START_REF] Yang | PAML 4: phylogenetic analysis by maximum likelihood[END_REF]. The accumulation of synonymous substitutions (saturation of dS) may create a problem to detect positive selection with sequences from highly divergent species such as human and zebrafish with some methods. However, a previous study showed that the method we use to detect positive selection, i.e. maximum likelihood estimate, is robust to dS saturation [START_REF] Studer | Pervasive positive selection on duplicated and nonduplicated vertebrate protein coding genes[END_REF]. The calculation was possible for 48 genes and we found amino acids under positive selection for 33 of them (Figure 1). Among teleosts, birds and eutherians taxa, we found various patterns of positive selection. The FN1 gene has sites under positive selection in 7 species for which data are available. The number of genes under positive selection in each species ranges from 0 (horse and human) to 10 (opossum). Positive selection occurs in the three groups of genes, suggesting the presence of an evolutionary pressure in the three steps of fertilization.

Our comparison of the intensity of positive selection shows that genes involved in sperm-oolemma fusion, including genes of the integrin family (IGSF8, ITGa9, ITa v , ITGa6), are mostly subjected to positive selection in teleosts, but much less so in mammals (Figure 2b). The difference between teleosts and birds (p,0.05) and between teleosts and eutherian mammals (p,0.001) are statistically significant whereas the difference between birds and eutherian mammals is not. This result is surprising for teleosts because we expected to find more amino acids under positive selection for genes involved in the sperm-ECM binding, which would contribute to the first barrier in free-spawning conditions. In particular, positive Darwinian selection has been shown to play a role in the rapid evolution of genes involved in sperm-egg binding in three other free-spawning taxa, sea urchin, mussel and abalone [START_REF] Aagaard | ZP domain proteins in the abalone egg coat include a paralog of VERL under positive selection that binds lysin and 18-kDa sperm proteins[END_REF][START_REF] Mccartney | Adaptive Evolution of Sperm Bindin Tracks Egg Incompatibility in Neotropical Sea Urchins of the Genus Echinometra[END_REF][START_REF] Riginos | Positive selection on an acrosomal sperm protein, M7 lysin, in three species of the mussel genus Mytilus[END_REF][START_REF] Vacquier | Positive Darwinian selection on two homologous fertilization proteins: what is the selective pressure driving their divergence?[END_REF]. The comparisons of the intensity of positive selection between taxa for sperm-ECM binding and between taxa for acrosome reaction are not significant. When analyzed globally, genes involved in the three steps of gamete interaction show a similar, but less marked, pattern in which the only statistically significant difference is between teleosts and eutherian mammals (p,0.05) (Figure 2c). Besides, in contrast to previous studies [START_REF] Swanson | Positive Darwinian selection drives the evolution of several female reproductive proteins in mammals[END_REF], we did not find positive selection for ZPC/ZP3, in agreement with more recent data obtained on a larger number of species [START_REF] Berlin | Testing for adaptive evolution of the female reproductive protein ZPC in mammals, birds and fishes reveals problems with the M7-M8 likelihood ratio test[END_REF]. We also show that several genes are positively selected in more than one species, such as ACE that exhibits positive selection in chimpanzee, opossum, Xenopus, and zebrafish.

Amino Acids under Positive Selection are not Located in known Functional Domains of Proteins

In order to project our previous results of positive selection on 3D structures, we have built 23 tridimensional models of 6 proteins (ACE, ARSA, C1QBP, ITa V , PDIA3 and UCHL3) from species where positive selection has been detected (Figure 3 and Figure S1 to S5). These 6 proteins were chosen because PDB template structures were available [START_REF] Berman | The Protein Data Bank[END_REF]. Firstly, we show that amino acids under positive selection are never located in the domain involved in known biological role (for ACE, ARSA, PDIA3 and UCHL3). This suggests that these residues are involved in other secondary functions such as interaction with partners or in the folding of the protein. Secondly, we show that amino acids under positive Rate of events per million years (Myr) per lineage for the set of genes studied. This includes up to 69 genes, but in most cases, the number of relevant genes is lower because some types of events cannot happen in some taxa (losses or duplications cannot happen if a gene is primitively absent, and gains cannot occur in a clade if the gene appeared before the base of a given clade), and we have incorporated that factor into our calculations to obtain comparable rates across taxa. The geological age of lineages is used in the rate calculation through the phylogenetic diversity of the clades. b and c. Proportion of genes under positive selection in Teleostei, Aves and Eutheria for genes belonging to the third step, the sperm-oolemma fusion (B), and for all genes (C). Three types of confidence intervals are shown: 78%, 90% and 97% representing testing at 0.05, 0.01 and 0.001 probability thresholds, respectively. Statistically significant comparisons, for which confidence intervals do not overlap, are indicated. doi:10.1371/journal.pone.0044548.g002 selection, for the same protein, have different location between species. Thus, we suppose that the function of these proteins in the process of fertilization may differ between species. This hypothesis is supported by a recent study that has emphasized the importance of structural architecture of the egg coat and its importance in fertilization [START_REF] Monne | A Structural View of Egg-Coat Architecture and Function in Fertilization[END_REF]. Thus, we postulate that positive Darwinian selection could contribute to slight modifications of the architecture of egg or sperm coat by changing key amino acids in some taxa. These architectural modifications may also create or reinforce reproductive barriers. Experimental investigations must be done (i) to identify amino acids involved in sperm-egg interaction in these proteins, (ii) to study the expression and the role of these genes in gametes and (iii) to identify their function in sperm-egg interaction in a greater number of species.

In conclusion, we study the evolution of 69 genes involved in the sperm-egg interaction. This work should help to improve our understanding of the evolution and the mechanism of fertilization in vertebrates, which still remains puzzling. In particular, the evolution of proteins present on male and female gametes can lead to reproductive barriers at the level of sperm-egg interaction in vertebrates. The variation between species regarding positive selection patterns as well as gain or loss of certain genes involved in fertilization could actually lead to the appearance or the reinforcement of reproductive barriers. The relationship between positive selection and molecular adaptation has been experimentally shown previously [START_REF] Sawyer | Positive selection of primate TRIM5alpha identifies a critical species-specific retroviral restriction domain[END_REF][START_REF] Weinberger | Positions under positive selection -key for selectivity and potency of scorpion {alpha}-toxins[END_REF]. Our work provides a list of potential targets, i.e. amino acids under positive selection, for targeted mutagenesis experiments to see if these amino acids are indeed involved in reproductive barriers between species. However, this list of genes should be completed after further experiments on other species. All this should improve significantly our understanding of the appearance of prezygotic isolation mechanisms, and hence, of the speciation process, which is of paramount importance in generating biodiversity.

Materials and Methods

Selection of Data

This study has sampled the genome of 19 bony vertebrate species that have been fully sequenced (Bos taurus, Sus scrofa, Equus caballus, Canis familiaris, Rattus norvegicus, Mus musculus, Homo sapiens, Pan troglodytes, Macaca mulatta, Monodelphis domestica, Ornithorhynchus anatinus, Gallus gallus, Taeniopygia guttata, Xenopus tropicalis, Tetraodon nigroviridis, Takifugu rubripes, Gasterosteus aculeatus, Oryzias latipes and Danio rerio). For all identified genes, the corresponding Ensembl protein ID was retrieved from Ensembl database and submitted to for four species based on human PDB structure (PDB: 1O8A) [START_REF] Natesh | Crystal structure of the human angiotensin-converting enzyme-lisinopril complex[END_REF]. This enzyme is composed of two homolog N and C domains, each containing a catalytic site composed by two histidines in a conserved -HEXXH-zinc-binding motif [START_REF] Soubrier | Two putative active centers in human angiotensin I-converting enzyme revealed by molecular cloning[END_REF] and a glutamate as the third ligand 24 residues downstream this motif [START_REF] Williams | Identification of two active site residues in human angiotensin I-converting enzyme[END_REF]. These three binding sites are indicated in yellow, whereas amino acids under positive selection are in pink. ACE protein is expressed by sperm. Gene knockout experiments of ACE in mouse results in an impaired uterotubular sperm migration and a reduced ability to bind the zona pellucida of the egg [START_REF] Esther | Mice lacking angiotensin-converting enzyme have low blood pressure, renal pathology, and reduced male fertility[END_REF][START_REF] Hagaman | Angiotensin-converting enzyme and male fertility[END_REF][START_REF] Krege | Male-female differences in fertility and blood pressure in ACE-deficient mice[END_REF]. doi:10.1371/journal.pone.0044548.g003

the PhyleasProg web server [START_REF] Busset | PhyleasProg: a user-oriented web server for wide evolutionary analyses[END_REF]. For each submitted Ensembl protein ID, the PhyleasProg web server interrogates the Ensembl database to retrieve orthologs and paralogs, the orthology/ paralogy relationships assignment being based on a phylogenetic approach in Ensembl. All phylogenetic trees were reconstruct by PhyleasProg using the TreeBeST program (http://treesoft. sourceforge.net/treebest.shtml) and carefully examined before interpreting selective pressure results, also calculated by Phylea-sProg. Selective pressure results are systematically verified and eventually corrected by synteny analysis as previously described by Tian, et al. [START_REF] Tian | Gene birth, death, and divergence: the different scenarios of reproduction-related gene evolution[END_REF]. Thus, calculations were performed with supposed correct orthologs. To score the presence, number and type of genes in the data matrix, a search for pseudogenes was systematically performed in the concerned genomes for genes for which no ortholog has been identified in at least one of the species of interest. The pseudogene status was inferred for a gene in a genome only if we have found a stop codon or an indel in the sequence identified by the similarity search in the syntenic locus in comparison with the other species of interest. Moreover, a search for duplicates was realized by using phylogenetic trees available in Ensembl and accessible via the Ensembl ID used for PhyleasProg and summarized in the Table S1.

Inference of Positive Selection

The inference of selection was performed by PhyleasProg with branch-site models of codeml of the PAML package, with each ortholog branch tested for positive selection. Multiple alignments were systematically and carefully examined to avoid all false positive results. In particular, amino acids predicted to be under positive selection that were at the boundary of the alignments were not considered because doubtful. We also eliminated genes for which positive selection was due to sequence errors in Ensembl according to a comparison with other available sequences from other database such as RefSeq in NCBI (Figure S6 to S8). Each branch of each phylogenetic tree was tested for positive selection. So we performed multiple test corrections by controlling for the false discovery rate (FDR) using the R package QVALUE [START_REF] Storey | Statistical significance for genomewide studies[END_REF]. Results are considered significant with a threshold of q = 10% of false positives. Sites with posterior probabilities of Bayes Empirical Bayes analyses superior to 95 or 99% were considered as positively selected. Datasets with less than 10 sequences, the minimum threshold required to obtain significant results, with excessively divergent sequences, or with sequences of genes for which annotations are not reliable were not retained for subsequent analysis (Figure 1). MODELLER [START_REF] Eswar | Protein structure modeling with MODELLER[END_REF] was used to build homology models of the 3D structure of proteins when possible and amino acids under positive selection and essential functional residues were localized on them.

Comparison of Evolutionary Rates and the Intensity of Positive Selection Across Clades

Character history was inferred using parsimony optimization onto our rooted tree. In most cases, this yielded an unambiguous polarity. In a few cases, such as for the character ZPB/ZP4, the ancestral condition could not be unambiguously resolved using parsimony; it was then resolved using additional criteria. For instance, in that case, one of the two main clades, Teleostei, had the state "absent", whereas the other main clade, Tetrapoda, had the state "single copy". In such cases, we assumed that absence was primitive, and scored a gain for Tetrapoda. All cases of losses located on other parts of the tree (not on one of the two daughterbranches of the basalmost node) were inferred from parsimony optimizations whose polarity, in that case, is unproblematic. Similarly, when one of the two basalmost branches can be inferred to have had a single copy whereas the other could most parsimoniously have two or more, we inferred that the presence of a single copy was primitive, based on the principle that the most common condition is likely to be primitive. This situation was encountered in the character ZPAX, for which Tetrapoda have a single copy, whereas Teleostei can most parsimoniously be inferred to have had two or three. To assess evolutionary rates, we optimized the 69 discrete characters representing gene presence, number and type (functional gene or pseudogene) onto the reference tree (Figure S9) and scored the occurrence of events of four types (gene appearance, duplication, transformation into a pseudogene, or loss) in three clades (Teleostei, Aves, Eutheria). We considered the basal branch subtending each of these three taxa as being part of it, and this is consistent with the fact that 19 sampled taxa in this study do not include representatives of both basalmost branches of these clades (e.g. paleognaths are not among the two birds sampled). We then divided the number of events of each type in each clade by that clade's sampled phylogenetic diversity index [START_REF] Faith | Conservation evaluation and phylogenetic diversity[END_REF] (the sum of branch lengths, including of the basal branch). We then divided the resulting number by the proportion of the 69 genes in which such changes could occur in each clade. The last operation compensates for the fact that several genes present in eutherian mammals are absent in teleosts, and this avoids underestimating the rate of gene duplication in teleosts. Thus, origin of these genes could occur in teleosts, but not duplication, inactivation, or deletion. To statistically compare the rates between teleosts and eutherians, we took into consideration the sampled phylogenetic diversity of each clade and the number of genes for which each event was possible to assess the relative probability of occurrence of each event under the null hypothesis that these events are clade-independent. We then assessed the probability that a distribution at least as extreme as the observed one arises under the null hypothesis using a binomial test implemented in GraphPad (http://www.graphpad. com/quickcalcs/binomial2.cfm).

To compare the intensity of positive selection in the three clades of interest (Teleostei, Aves, Eutheria), we computed the inferred ancestral value of each clade using squared-change parsimony [START_REF] Maddison | Squared-change parsimony reconstructions of ancestral states for continuous-valued characters on a phylogenetic tree[END_REF] in Mesquite [START_REF] Maddison | Mesquite: a modular system for evolutionary analysis[END_REF]. Confidence intervals of these values were computed using phylogenetic independent contrasts [START_REF] Felsenstein | Phylogenies and the comparative method[END_REF] using the PDAP:PDTREE module [START_REF] Midford | PDAP Package for Mesquite[END_REF] of Mesquite. These were used to assess statistical significance of the differences in ancestral values of each clade. Our method is derived from the continuous analysis developed for heterochrony analysis of continuous characters [START_REF] Germain | Evolution of ossification sequences in salamanders and urodele origins assessed through event-pairing and new methods[END_REF], a method whose statistical properties (especially power and type I error rate) have already been assessed using simulations. However, we improved a bit the method to more fully incorporate uncertainty about nodal estimates; instead of comparing a nodal value with the confidence interval (CI) of another clade, we compare the confidence intervals of both clades. This allows using narrower CIs because the probability that both nodal values actually lie outside the computed CI at a given error threshold is equal to the square of that threshold. As a result, instead of comparing a nodal value with the 95% CI (0.05 probability threshold), we compared the CIs at 78% of both clades because the square of 0.22 is just under 0.05 (we rounded off to the nearest higher percentage, which should make the test slightly conservative). Thus, if the two compared CIs do not overlap, we conclude that the difference is statistically significant. This method thus incorporates uncertainty in the estimates of all compared nodal values. The fact that data from all taxa are used to compute ancestral values should also make the test conservative because this should tend to make ancestral estimates of sister-taxa more similar to each other than if the ancestral values were estimated for each clade separately.

Figure 1 .

 1 Figure 1. Signature of evolution of genes involved in sperm-oocyte interaction in Vertebrates. Genes (names are indicated on the top) are classified into three main functional groups (from left to right): sperm-zona pellucida binding, acrosome reaction, and sperm-oolemma fusion. Species are indicated in the tree (left). Blue square: duplication; purple square: positive selection; pink square: duplication + positive selection. Black square: pseudogene; Black dot: no event. Black circle: no calculation of positive selection. White square: no gene found. doi:10.1371/journal.pone.0044548.g001

Figure 2 .

 2 Figure 2. Rate of gene appearance, gene duplication, pseudogenization, loss and intensity of positive selection in Teleostei, Aves and Eutheria. a.Rate of events per million years (Myr) per lineage for the set of genes studied. This includes up to 69 genes, but in most cases, the number of relevant genes is lower because some types of events cannot happen in some taxa (losses or duplications cannot happen if a gene is primitively absent, and gains cannot occur in a clade if the gene appeared before the base of a given clade), and we have incorporated that factor into our calculations to obtain comparable rates across taxa. The geological age of lineages is used in the rate calculation through the phylogenetic diversity of the clades. b and c. Proportion of genes under positive selection in Teleostei, Aves and Eutheria for genes belonging to the third step, the sperm-oolemma fusion (B), and for all genes (C). Three types of confidence intervals are shown: 78%, 90% and 97% representing testing at 0.05, 0.01 and 0.001 probability thresholds, respectively. Statistically significant comparisons, for which confidence intervals do not overlap, are indicated. doi:10.1371/journal.pone.0044548.g002

Figure 3 .

 3 Figure 3. 3D structures of N and C domain of ACE in platypus, Xenopus, chimpanzee and zebrafish. ACE structures have been modeled for four species based on human PDB structure (PDB: 1O8A)[START_REF] Natesh | Crystal structure of the human angiotensin-converting enzyme-lisinopril complex[END_REF]. This enzyme is composed of two homolog N and C domains, each containing a catalytic site composed by two histidines in a conserved -HEXXH-zinc-binding motif[START_REF] Soubrier | Two putative active centers in human angiotensin I-converting enzyme revealed by molecular cloning[END_REF] and a glutamate as the third ligand 24 residues downstream this motif[START_REF] Williams | Identification of two active site residues in human angiotensin I-converting enzyme[END_REF]. These three binding sites are indicated in yellow, whereas amino acids under positive selection are in pink. ACE protein is expressed by sperm. Gene knockout experiments of ACE in mouse results in an impaired uterotubular sperm migration and a reduced ability to bind the zona pellucida of the egg[START_REF] Esther | Mice lacking angiotensin-converting enzyme have low blood pressure, renal pathology, and reduced male fertility[END_REF][START_REF] Hagaman | Angiotensin-converting enzyme and male fertility[END_REF][START_REF] Krege | Male-female differences in fertility and blood pressure in ACE-deficient mice[END_REF]. doi:10.1371/journal.pone.0044548.g003

PLOS ONE | www.plosone.org

September 2012 | Volume 7 | Issue 9 | e44548

Acknowledgments

We acknowledge Pierre Pontarotti, Hugues Roest Crollius and Marc Robinson-Rechavi for helpful discussion and useful comments on the manuscript.

CM is funded by a MENRT PhD fellowship. This work was supported by INRA. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Supporting Information

Figure S1 3D structure of ARSA in chicken, chimpanzee and medaka. The three structures were modeled with human ARSA structure as template (PDB: 1N2K) [105]. The active site of the enzyme is represented in yellow, amino acids under positive selection in pink. ARSA is expressed by the spermatozoa; its involvement in the fertilization process was demonstrated with antibodies directed against ARSA. When the sperm is pretreated with these antibodies, their ability to bind the zona pellucida is reduced [START_REF] Rattanachaiyanont | Anti-SLIP1-reactive proteins exist on human spermatozoa and are involved in zona pellucida binding[END_REF]. (TIF)

Figure S2 3D structures of C1QBP in opossum and zebrafish. Both structures were modeled with human P32 as template (PDB: 1P32) [106]. Amino acids under positive selection are in pink. The C1QBP glycoprotein is localized on sperm; its participation in the fertilization process was demonstrated using anti-C1QBP antibodies. The interaction between sperm and zona pellucida is suppressed when the sperm is pretreated with these antibodies UCHL3 protein is present on both sperm acrosome and oocyte cortex. A study showed that UCHL3 is an ubiquitin C-terminal hydrolase and is involved in antipolyspermy defense during porcine fertilization [START_REF] Yi | Ubiquitin Cterminal hydrolase-activity is involved in sperm acrosomal function and antipolyspermy defense during porcine fertilization[END_REF].

(TIF)

Figure S6 An example of a sequence error in Ensembl.

The sequence presented in (a) is the protein sequence of the gene PDIA3 from Xenopus. Amino acids in red and orange are those which have been detected as positively selected by PhyleasProg. A tBLASTn has been performed with this sequence against EST database of NCBI. The result of the BLAST is presented in (b), the domain which contains amino acids under positive selection in (a) is not retrieved. A new sequence for Xenopus was searched in NCBI RefSeq database. The sequence of the best BLAST hit is submitted to a tBLASTn against the EST database of NCBI for verification (c). The new identified sequence from NCBI replaced in this case the sequence from Ensembl. (TIF)

Figure S7 An example of a doubtful amino acid under positive selection. The amino acid in red in the sequence of the rat (Rattus norvegicus) was found as positively selected in the alignment of protein sequences of the gene CD46. Because it is located at the boundary of the alignment, this amino acid is not considered, because doubtful. (TIF)