K Guizien
email: guizien@obs-banyuls.fr

M Belharet

C Moritz

J M Guarini

Vulnerability of marine benthic metapopulations: implications of spatially structured connectivity for conservation practice

Keywords: larval dispersal, metapopulation, connectivity, numerical simulation, vulnerability, resistance, natural and anthropogenic stresses

Aim: Connectivity is currently considered to essential to enhance biodiversity conservation efforts and benefit adjacent areas. We evaluate the spatial structure of vulnerable marine benthic invertebrate populations based on connectivity to improve the placement of marine protected areas.

Location: Gulf of Lions, NW Mediterranean Sea.

Methods: We used a spatially explicit metapopulation model to explore how larval dispersal affects marine benthic invertebrate dynamics at local and regional scales.

Minimum recruitment success (the proportion of larvae that settle in a site and survive between spawning and first reproduction, which is required to ensure species persistence) is proposed as a measure of vulnerability. Three contrasting simulations were constructed: closed versus connected populations, habitat loss versus recruitment failure disturbances, and varying patterns of oceanographic connectivity. The rescue factor (the ratio of minimum recruitment success in connected and closed populations) quantified the decrease in vulnerability due to oceanographic connectivity.

Results: In the Gulf of Lions, connectivity reduced population vulnerability by a rescue factor of 5.3. The minimum recruitment success for regional persistence decreased from 3 % to 0.2 % when species' life expectancies increased from 2 to 50 years. Connectivity enabled the spread of individuals over the entire region, but there were higher densities in sink sites, where low retention rates do not allow for local persistence (western tip of the gulf). Source sites, defined as those sites bringing metapopulation resistance to recruitment failure and stresses from habitat loss, were located in the center of the gulf.

Main conclusions: Connectivity drives the spatial structure of regional adult distribu-tions and population vulnerability in the absence of any other structuring factor. In the Gulf of Lions, marine protected areas are located in the western tip, where there are population sink; however, if protection measures were taken in the central part, populations would benefit throughout the region via the rescue effect.

INTRODUCTION

Landscape connectivity is considered an essential factor for enhancing the resilience of biodiversity and ecosystem services in the face of climate change [START_REF] Heller | Biodiversity management in the face of climate change: a review of 22 years of recommendation[END_REF][START_REF] Morecroft | Resilience to climate change: translating principles into practice[END_REF]. The role of connectivity in biodiversity resilience (defined as the long-term persistence of a species) has been extensively examined theoretically in conservation biology and population ecology. Connectivity alters biodiversity in many ways, with both negative and positive outcomes. On the one hand, connectivity linked to anthropogenic actions (e.g. the opening of commercial routes, species introductions) has enabled the spread of alien species [START_REF] Hulme | Trade, transport and trouble: managing invasive species pathways in an era of globalization[END_REF][START_REF] Katsanevakis | Invading European Seas: Assessing pathways of introduction of marine aliens Ocean and[END_REF], which has been identified as one of the major drivers of biodiversity loss (MEA, 2005). On the other hand, connectivity associated with propagule dispersal during reproduction has favored the wide geographical distribution of sessile species, which are dependent on the availability of suitable habitat downstream (Scheltema, 1986). Connectivity has also been recognized as a process that can reduce extinction threats due to invasive species by linking species refuges [START_REF] Gilbert | Plant invasions and extinction debts[END_REF]. Despite the adaptive disadvantage represented by dispersal to new sites, larval dispersal is also a way for sessile populations to recover following disturbances, such as extreme meteorological events or anthropogenic stresses [START_REF] Baker | Climate change and coral reef bleaching: An ecological assessment of long-term impacts, recovery trends and future outlook[END_REF][START_REF] Miller | Genetic diversity and gene flow in collapsed and healthy abalone fisheries[END_REF][START_REF] Kendrick | The Central Role of Dispersal in the Maintenance and Persistence of Seagrass Populations[END_REF][START_REF] Bernhardt | Resilience to Climate Change in Coastal Marine Ecosystems[END_REF].

Since the earliest metapopulation model of [START_REF] Levins | Evolution in Changing Environments[END_REF], a wealth of stochastic patch occupancy models have been used to address the fundamental question of the consequences of habitat fragmentation on the persistence of individual (Quinn and Hastings, 1987;[START_REF] Hanski | Metapopulation dynamics: Does it help to have more of the same ?[END_REF][START_REF] Hanski | A practical model of metapopulation dynamics[END_REF] and multiple species (Mouquet & Loreau, 2002). Increased satellite imaging capabilities have enabled the production of high-resolution maps of species' habitats in terrestrial environments, which have improved the spatially explicit modeling of species distribution (Schnell et al., 2013). Such an approach is limited in the marine environment, where habitat mapping is still incomplete. Conversely, dispersal is often described simplistically in terrestrial environments and is based solely on the distance between patches of a species' habitat (Roy et al., 2008;[START_REF] Dornier | Colonization and extinction dynamics of an annual plant metapopulation in an urban environment[END_REF]Skelsey et al., 2013;[START_REF] Johansson | Estimates of connectivity reveal nonequilibrium epiphyte occurence patterns almost 180 years after habitat decline[END_REF]. Yet in the case of propagule dispersal (pollen in the air, larvae in the ocean), transport mediated by an unsteady and non-uniform flow will often be asymmetrical [START_REF] Aiken | Environmental fluctuations and asymmetrical dispersal: generalized stability theory for studying metapopulation persistence and marine protected areas[END_REF] and can cover long distances. In the marine environment, larval transfer between sites has been calculated from larval dispersal simulations and summarized by oceanographic connectivity matrices [START_REF] Melbourne-Thomas | A multi-scale biophysical model to inform regional management of coral reefs in the western Philippines and South China sea[END_REF][START_REF] Thomas | Bivalve larvae transport and connectivity within the Ahe atoll lagoon (Tuamotu Archipelago), with application to pearl oyster aquaculture management[END_REF][START_REF] Jacobi | Identification of subpopulations from connectivity matrices[END_REF][START_REF] Morecroft | Resilience to climate change: translating principles into practice[END_REF][START_REF] Moritz | Disentangling the role of connectivity, environmental filtering, and spatial structure on metacommunity dynamics[END_REF]. Such connectivity matrices have been analyzed to elucidate population persistence using either graph theory methods [START_REF] Treml | Modeling population connectivity by ocean currents, a graph-theoretic approach for marine conservation[END_REF] or eigenvalue decomposition [START_REF] Hastings | Persistence of spatial populations depends on returning home[END_REF]. However, local demography may mitigate oceanographic connectivity [START_REF] Moilanen | On the limitations of graph-theoretic connectivity in spatial ecology and conservation[END_REF]; thus, assessing population persistence requires the integration of oceanographic connectivity and local demography into models of metapopulation dynamics [START_REF] Watson | Changing seascapes, stochastic connectivity and marine metapopulation dynamics[END_REF][START_REF] Williams | Stochastic dispersal and population persistence in marine organisms[END_REF]. Of the local demographic pa-rameters, recruitment success (the proportion of settling larvae that survive from spawning to the adult stage) ultimately drives the persistence of a population, although flow does bring larvae to a settling area. Persistence requires a 'minimum recruitment success' to compensate for the loss of adults in the population during the same period. The order of magnitude of recruitment success is virtually unknown for most species due to limitations in observing the supply of settled larvae [START_REF] Grüss | Estimating local settler-recruit relationship parameters for complex spatially explicit models[END_REF], and this parameter is the missing link between dispersal and local demography. Integrating demographic connectivity in spatially explicit metapopulation models is a challenge for the practical design of conservation strategies for the marine environment [START_REF] Kaplan | Model-Based Assessment of Persistence in Proposed Marine Protected Area Designs[END_REF][START_REF] Moffitt | The utility and limitations of size and spacing guidelines for designing marine protected area (MPA) networks[END_REF]. It would be particularly interesting to identify sites, by considering the functional connections within self-persistent sites or between sites, where protection efforts aimed at minimizing adult loss would ensure persistence [START_REF] Berumen | Persistence of self-recruitment and patterns of larval connectivity in a marine protected area network[END_REF]Olds et al., 2014).

Furthermore, if persistence is achieved within the protected network, then dispersal export may further benefit populations in non-protected areas [START_REF] Harrison | Larval export from marine reserves and the recruitment benefit for fish and fisheries[END_REF].

In this paper, we present the results of numerical simulations of the spatio-temporal dynamics of soft-bottom worm populations in the Gulf of Lions driven by oceanographic connectivity. The generic metapopulation model of [START_REF] Hastings | Persistence of spatial populations depends on returning home[END_REF], which incorporates spatially structured parameters for survival, fecundity and dispersal was modified to account for density-dependent recruitment (Roughgarden et al., 1985) and time-varying connectivity. A set of twenty different connectivity matrices were estimated from larval dispersal simulations and were randomly picked in the simulations of the metapopulation dynamics to reflect the inter-annual variability in dispersal that arises from the conjunction of environmental (currents) and biological (spawning timing) variability [START_REF] Guizien | Using larval dispersal simulations for Marine Protected Area design: application to the Gulf of Lions (NW Mediterranean)[END_REF]. 'Minimum recruitment success', which is required to ensure populations persistence under each oceanographic connectivity scenario was used as a measure of population vulnerability. Simulations were designed to test three hypotheses:

(1) Connectivity decreases population vulnerability at the regional scale, (2) population vulnerability to stresses is spatially structured by connectivity and (3) adult spatial distribution allows for the identification of sink sites but not sites essential for regional persistence (source). The findings were also used to propose an extension of the current protection measures in the Gulf of Lions.

METHODS

The generic metapopulation model for benthic invertebrates

The spatio-temporal dynamics of benthic invertebrate populations with a dispersive pelagic larval stage were simulated after modifying the generic metapopulation model of [START_REF] Hastings | Persistence of spatial populations depends on returning home[END_REF]. The model describes explicitly, in discrete times and for a set of sites with different areas connected by larval transfer, the spatial density dynamics of only the sedentary adult stage (Fig. 1). Density dynamics at a given site result from spatially structured local survivorship and reproductive success input potentially contributed by all the other sites in the system. The original model of [START_REF] Hastings | Persistence of spatial populations depends on returning home[END_REF] was modified to account for (1) recruitment limitation due to space availability at the destination site (computed as the proportion of free space based on the saturating density of adults N max = 1/α A , where α A is the mean cross-sectional area of one adult), and (2) the variability in propagule transfer rate. The model in matricial form is as follows:

N(t + ∆t) = min(G(t)N(t), N max ) (1)
with a time step of ∆t and a transfer matrix G defined by:

G ji = b j C ji (t)a i + s jj δ ji (2)
where N(t) is a vector ∈ R p containing the spatial density of adults at each site Propagule transfer rates are the elements of a square connectivity matrix C and are defined as the proportion of larvae coming from an origin site (columns) and arriving at a settlement site (rows), multiplied by the ratio of the areas of the origin to the settlement sites.

Recruitment success is the proportion of larvae settling in one site that has survived dispersal and will finally reach sexual maturity. This parameter, with values of between 0 and 1 and accounting for all mortality losses since egg release until the first reproduction by the new recruits, is unknown for all species and is very difficult to measure. An alternative approach is to seek a minimum value of this parameter that would avoid either population (based only on local retention) or metapopulation (allowing for transfer within the network of populations) extinction. The minimum recruitment success of a population corresponds to the recruitment success that yields a population growth rate equal to 1. Thus, the minimum recruitment success co-varies with the fecundity value; when the fecundity value is increased by a factor of 10, the minimum recruitment success is decreased by a factor of 10. Such a minimum recruitment success provides a measure of population vulnerability: the lower the minimum recruitment success, the lower the population vulnerability. The minimum recruitment success is independent of the saturating capacity choice because the site carrying capacity is not reached. The rescue factor was defined as the ratio between the recruitment success required to maintain individuals in any connected site in the region (regional persistence) and the minimum recruitment success required to maintain individuals in at least one isolated site in the region (local persistence). The rescue factor is, by definition, independent of species fecundity.

Application to the Gulf of Lions soft-bottom polychaete populations

The Gulf of Lions is a wide, micro-tidal continental shelf, dominated by soft-bottom habitat and it is delimited along the shore at its two extreme tips by rocky beds extending to a depth of 40 m [START_REF] Aloisi | Carte gologique du prcontinent languedocien au 1/250000ime[END_REF]. Polychaetes dominate the soft-bottom benthic communities of the gulf and their community composition is related to the grain size distribution and organic matter content, which is strongly related to water depth in this region except for some restricted areas at the river outlets (Labrune et al., 2007a). In the simulations, demographic parameters were assumed to be spatially uniform for the 10-30 m depth range to more clearly reveal the effect of connectivity on spatial structure.

Species were differentiated by their life expectancy (ranging from 2 to 50 years) and their mean fecundity (averaged over the duration of their reproductive life). According to the review of demographic parameters of 200 polychaete species by [START_REF] Mchugh | Do life history traits account for diversity of polychaete annelids ?[END_REF], no trade-off can be found between fecundity and life-expectancy, but fecundity is related to larval development mode. Fecundity values for polychaete species with dispersive larval stages range from 10 3 to 10 6 [START_REF] Mchugh | Do life history traits account for diversity of polychaete annelids ?[END_REF]. As fecundity can be related to food availability, fecundity in the oligotrophic Mediterranean sea is expected to be in the lower part of the range. Fecundity was thus varied between 10 3 and 10 5 eggs per female in the simulations, resulting in a propagule production rate of 5 10 1 to 5 10 3 larvae per female (using a sex ratio of 1 : 1 and a fertilization rate of 10%), independent of species life expectancy. As site saturating capacity is an upper limit introduced to avoid unrealistic population explosion when population growth rate is greater than 1 and is not expected to be reached at equilibria, this parameter, by definition, should be larger than the natural density. Furthermore, to avoid the introduction of an artificial trade-off between life expectancy and site saturating density, site saturating density was overestimated and set to 25,000 ind m -2 , which corresponds to a value larger than any density observed in the area (Labrune et al., 2007a). Given these demographic parameters, recruitment success was also considered to be spatially uniform and varied from 10 -6 to 1 (10 regular steps within each log 10 -interval) to detect the minimum value R at which either local or regional persistence was reached (see above).

The continuous sand facies bed spanning the 10 to 30 m isobaths along the gulf was discretized into thirty-two contiguous sites to describe the spatial distribution of the different worm species living in this habitat (Fig. 2). Site boundaries were placed according to the natural coastline features forming bays or approximately every 20 km along the rectilinear shoreline, which yielded sites ranging in area from 6.2 to 121 km 2 . Connectivity matrices, with local retention rates on the diagonal, and transfer rates between distant sites off of the diagonal, were quantified from Lagrangian larval dispersal simulations.

Simulations of larval dispersal in the gulf used three-dimensional currents and turbulent kinetic energy computed at a resolution of 750 m × 750 m following interpolation between the hourly outputs [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]. Spawning was simulated by releasing thirty neutrally buoyant larvae in the center of each of the 32 sites, on the 20 m isobath, every hour during the reproductive season. As temperate polychaetes reproduce in late winter-early spring in the NW Mediterranean Sea [START_REF] Bhaud | Quelques données sur le déterminisme écologique de la reproduction des Annélides Polychètes[END_REF], a 100-day reproductive season (from January 5 at 0h until April 13 at 23h) was considered and divided into ten consecutive 10-day spawning events in 2004 and 2006 (see details in [START_REF] Guizien | Using larval dispersal simulations for Marine Protected Area design: application to the Gulf of Lions (NW Mediterranean)[END_REF]. In the absence of knowledge regarding sinking behavior that might differentiate species' dispersal capabilities [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], connectivity was set as the same for all species in this study because most benthic invertebrates have a similar competency period of between 3 and 5 weeks (Mc Edward, 1995). Three connectivity matrices corresponding to three different pelagic larval durations (3-, 4-and 5-week) were built for each spawning period, and the proportion of larvae coming from the center of any site and arriving in each of the thirty-two sites was multiplied by the ratio of the areas of origin to the destination sites. A composite connectivity matrix was then calculated as the weighted sum of the 3-(1/4), 4-(1/2) and 5-(1/4) week pelagic larval duration connectivity matrices. In total, twenty variant composite matrices were used to describe the interannual variability in connectivity arising from a conjunction between meteorological variability and the timing of spawning during a reproductive season, and this variability was assumed to affect all polychaete species. Then, the average and standard deviation of these twenty variant matrices were computed, and used as the average conditions matrix in some of the simulations described in the next section.

Simulation groups and statistical methods

Three groups of spatio-temporal simulations over a duration of 50 life-expectancy cycles were carried out, starting with a low uniform spatial density of 0.4% of the site's carrying capacity. Asymptotic stable equilibrium was defined as the average of the population density between the 45 and 50 life-expectancy cycles. Connectivity was either set deterministically (the same matrix repeated every year) or stochastically (a random matrix picked up every year among the twenty variants). When connectivity was set stochastically, ensembles of 500 simulations were performed for each set of the other parameters in order to evaluate the variability arising from the connectivity sequence. Redundant patterns among the ensemble of 500 stochastic simulations in the regional spatial distribution were sought by grouping spatial distributions that deviated by less than 20% (square root of the sum of normalized squared deviations in each site). The spatial distribution of population density was described by the mean and standard deviation of the 500-simulation asymptotic stable equilibrium. The spatial distribution of relative density was described by the ratio between the population density in each site and the regional average of population density. The average regional coverage was defined as the ratio between the regional average of population density and the saturating capacity.

Simulation settings are summarized in Table 1. Simulation groups 1 and 2 explored the effect of connectivity on population persistence and spatial distribution. Simulation group 3 explored the resistance of a short lived species to two types of scenarios (Fig. 1B and C). By gradually increasing the number of unsuitable sites starting from each port, group 3A aimed to quantify the resistance of the metapopulation to habitat loss around the four main ports: Port-Vendres, Port La Nouvelle, Sète and Marseille (Fig. 2). Group 3B aimed to quantify the resistance of the metapopulation to different return periods of recruitment failure affecting either the Roussillon or the Languedoc Coasts. Recruitment failure is typically caused by extreme meteorological events (e.g. swell or heat waves, Durrieu de Madron et al., 2011) that affect populations over larger spatial scales than the scales affected by anthropogenic stresses.

RESULTS

Minimum recruitment success for local and regional persistence: the connectivity effect (group 1)

Connectivity matrices for the Gulf of Lions display a strong but highly variable spatial structure with variation coefficients for each connection spanning from 40 to 450 % (Fig. 3A andB). Amongst the twenty variants, three main connectivity structures could be identified and ranked according to drift direction: 65% with a dominant westward drift, 30% without dominant drift (diffusive connectivity) and 5 % with a dominant eastward drift. On average, the local retention rates (matrix diagonal terms) were higher in the central gulf (ranging from 0.2 % to 1.41 %) than in the eastern and western tips (< 0.05%). Sites located in the central gulf were also major larvae exporters, and can be split into two groups. Sites 10 to 13 predominantly exported individuals to the west, whereas sites 20 to 25 exported to both the east and the west (Fig. 3A). Larval transfer within and between the two groups of sites was also strong.

The spatial structure of local retention rates along the gulf coastline is depicted by the minimum recruitment successes required for the local persistence of species with life-expectancies ranging from 2 to 50 years and a fecundity of 10 4 eggs per female (Fig. 4 A). The minimum recruitment success would be multiplied by a factor of 10 for species with a fecundity of 10 3 eggs per female or divided by a factor of 10 for species with a fecundity of 10 5 eggs per female, but their relative values along the Gulf of Lions would be the same. In the nine sites along the western tip of the gulf (Roussillon coast, region of Port-Vendres) and also for three sites at its very eastern tip (region of Marseille), local retention rates were too small to achieve persistence for a species with a life expectancy of up to 10 years (up to 50-year life span in the five sites at the western tip). In contrast, the larger local retention rates in the central part of the gulf ensure the local persistence of the shortest-lived species (2-year) and of the longest-lived species (50-year) if recruitment success is larger than 70 % and 6 %, respectively. In fact, the recruitment success required for regional persistence without (in one isolated population with a maximum retention rate) and with connectivity decayed exponentially as species life expectancy increased, as it was divided by a factor of 2 between 2-and 7-year life expectancies (Fig. 4B).

Connectivity decreased vulnerability independently of the species life expectancy (Fig. 4B).

The rescue factor was independent of species life expectancy and varied from 1 to 6.6

among the twenty connectivity matrix variants (3.8 on average). When a connectivity matrix was picked randomly every year from the twenty variants, the minimum recruitment success required for regional persistence was the same in the 500 simulations and was slightly higher or equal to the value when the averaged connectivity was repeated every year (data not shown). The rescue factor was 5.3 for the averaged connectivity and 5.0 when the connectivity matrix was variable every year.

Connectivity between distant sites (i.e. larval import) tended to decrease the sensitivity of the persistence thresholds to the variability of the hydrody-namics. The minimum recruitment success required for regional persistence was highly variable among the twenty different hydrodynamic periods when only local retention was included, which reflected the large variability in local retention rates (Fig. 4 B, without connectivity). In contrast, the minimum recruitment success required for regional persistence was less variable among the connectivity matrix variants when larval import from distant sites was included, with the exception of one matrix characterized by very low transfer rate values (Fig. 4 B, with connectivity).

Connectivity and spatial distribution of population density at equilibrium (groups 1 and 2)

The relationship between connectivity matrix structure and the spatial distribution of population density at equilibrium is illustrated in this section for a species with a 2-year lifespan but would be the same for any other life expectancy provided that recruitment success was set to the minimum value required to ensure regional persistence. When the same connectivity matrix was used every year and regional persistence was ensured, the spatial distribution of population density was significantly correlated with the spatial distribution of cumulated import rates in the settlement sites of the corresponding connectivity matrix (on average, R = 0.78, p < 0.05). When larval transfer drifted to the east, density was much lower in the western part of the gulf (Fig. 5 A). When larval transfer drifted to the west, individuals were mainly distributed in the western part of the gulf (Fig. 5B). As larval transfer rates were of the same order of magnitude regardless of the different connectivity structures, regional persistence was ensured by similar values of minimum recruitment success. Interestingly, in both cases, individuals were not only present in the sites with the highest retention rates but also in sites where local retention was far too low to ensure local population persistence. When the averaged connectivity matrix was repeated every year in a simulation (deterministic case) with recruitment success set to its minimum value for regional persistence, density distribution was uneven in the gulf, which reflected the average spatial distribution of the cumulative import rates along the gulf (R = 0.94, p < 0.05). Individuals, although present everywhere, accumulated in the western half of the gulf with the highest densities in sites 13 to 17 and in sites 1 and 2 at the very western tip (Fig. 5C).

When connectivity was variable each year (conserving the averaged connectivity during the 100-year simulation) and recruitment success was set to the minimum value for regional persistence, density spatial distribution displayed the same accumulation pattern as in the deterministic case, whereas the relative density disequilibrium between the eastern and western part of the gulf was smoothed (Fig. 5D). It was not possible to identify a dominant pattern for the spatial distribution of population density, as the redundancy of the asymptotic equilibria density spatial distribution reached 4% at the maximum among the 500 simulations. In fact, spatial distribution of population density displayed higher variability among the 500 simulations than among years in each simulation. This variability was also spatially structured, and was larger in the locations identified as the accumulation sites of the averaged connectivity matrix (Fig. 5D) as these locations also correspond to the areas of highest transfer rate variability (Fig. 3A andB). When recruitment sucess was increased, densities increased and their spatial distribution became more even (data not shown).

Connectivity and resistance to habitats loss or periodic recruitment failure of a 2-year lifespan species (groups 3A and 3B)

The spatial structure of vulnerability to stresses affecting habitat availability or recruitment occurrence is described in this section only for a 2-year lifespan species, because the vulnerability of such short-lived species was already highest in the absence of any stress.

Wherever habitat loss was occuring (western, central or eastern part of the Gulf of Lions), the average regional coverage decreased with habitat availability and reached metapopulation extinction before all of the suitable habitats were lost (Fig. 6A andB). However, resistance to habitat loss displayed a strong spatial structure. Applying a recruitment success value of 3% (which was the minimum value required for regional persistence in the absence of stress), resistance to habitat loss was much lower in Sète (central part of the Gulf of Lions, where extinction occurred as soon as 3 sites were lost) than in Marseille (western tip, where extinction occurred when 10 sites were lost) and Port La Nouvelle (central part, where extinction occurred when 12 sites were lost) or finally Port-Vendres (eastern tip, where extinction occurred when 16 sites were lost, Fig. 6A). Yet ultimately, the spatial structure of the resistance to habitat loss can be shaped by the recruitment success value (Fig. 6A andB). When recruitment success was doubled, resistance to habitat loss was increased globally but the spatial structure of vulnerability to habitat loss was modified. For example, Port La Nouvelle became as resistant as Port-Vendres (Fig. 6B).

The spatial structure of vulnerability to habitat destruction was not sensitive to the history of the connectivity sequence, as variability lead to less than a 5% standard deviation in species average regional coverage amongst the 500 simulations, irrespective of where habitat loss started.

When recruitment failure affected the western part of the Gulf of Lions (Roussillon coast), populations never went extinct and average regional coverage was insensitive to recruitment failure frequency (Fig. 7A andB). In contrast, when recruitment failure affected the northern part of the Gulf of Lions (Languedoc coast), the species with a 2-year lifespan would go extinct with a recruitment failure return period greater than 4 years if recruitment success was the minimum required to ensure metapopulation persistence in the absence of stress (Fig. 7A). However, the metapopulation could resist recruitment failure of up to every other year (species life expectancy) when compensated by a higher recruitment success (Fig. 7B). In contrast to habitat loss, metapopulation resistance to recruitment failure was strongly dependent on the history of the connectivity sequence, when the northern part of the Gulf of Lions was impacted (the standard deviation in average regional coverage amongst the 500 simulations reaching approxiamtely 50%).

DISCUSSION

The role of population connectivity in species persistence is currently emphasized in conservation biology and is translated into policy in the current Strategic Plan for Biodiversity 2011-2020 which formulated a conservation goal for the marine environment through 'ecologically representative and well connected systems of protected areas, and integrated into the wider seascapes' (https://www.cbd.int/decision/cop/?id=12268). At the same time, the openness of marine populations is being reconsidered using numerical simulations of larval dispersal that suggest that the scales of dispersal may not be as large as previously thought, with drastic reductions in realized long-distance connectivity due to larval mortality [START_REF] Cowen | Connectivity of marine populations: open or closed?[END_REF]. However, this concept does not imply that the closed-population hypothesis is valid [START_REF] Cowen | Larval dispersal and marine population connectivity[END_REF]. By subdividing the nearshore area of the Gulf of Lions into consecutive sites of tens of kilometers in width, we showed that local persistence would require different values of local recruitment success to compensate for different retention rates along an unfragmented uniform soft-bottom habitat. This spatial structure was already identified in an earlier study concerned with the definition of self-persistent no-take reserves in this region [START_REF] Guizien | Using larval dispersal simulations for Marine Protected Area design: application to the Gulf of Lions (NW Mediterranean)[END_REF].

The findings of the present study confirm the very low or non existent potential for the conservation of soft-bottom species with a dispersive larval phase and a life span of up to 50 years in existing isolated marine protected areas located in the extreme tips of the Gulf of Lions (e.g. Parc Marin de la Côte Bleue, Réserve Naturelle Marine de Cerbère-Banyuls). These MPAs mainly target the conservation of hard-bottom species, which may experience higher retention rates in the boundary layer that develops around coralligeneous three-dimensional structures (Moulin et al., 2007). Of greater concern is the location of the recently created Parc Naturel Marin du Golfe du Lion (PNM, Fig. 2) which explicitly targets soft-bottom species as part of its sustainability goals for fisheries, while, according to our simulations, only including sink populations. The location of these sink populations is confirmed by the recurrent but transient extirpations of the 2-year life span serpulid Ditrupa arietina in Banyuls Bay followed by rapid recovery over the last 15 years [START_REF] Grémare | On the dramatic increase of Ditrupa arietina O.F. Müller (Annelida: Polychaeta) along the French and Spanish Catalan coasts[END_REF][START_REF] Medernarch | Propulation dynamics, secondary production and calcification in a Mediterranean population of Ditrupa arietina (Annelida: Polychaeta)[END_REF]Labrune et al., 2007a, F. Charles, Pers. Comm.). Furthermore, the presence of many soft-bottom polychaete spawners in the western part of the Gulf of Lions (Labrune et al., 2007a), where local retention rates are too low to ensure species persistence as isolated populations, also indicate the existence of efficient larval transfer at the meso-scale. The rescue factor is proposed as a metric to rank the decrease in population vulnerability due to connectivity (regional persistence is not achieved by local persistence at any site, but implies a larval flux between sites) in different regions and for different species. In the Gulf of Lions, the rescue factor for species with a 4-week planktonic larval duration reached a value of 5.3 on average. Interestingly, this value was only 6% less when connectivity was variable among years, in contrast to findings reported by [START_REF] Watson | Changing seascapes, stochastic connectivity and marine metapopulation dynamics[END_REF].

In the present study, scenarios of spatially structured disturbances based on existing meteorological and anthropogenic threats to a short-lived species were used to identify source sites that promote resistance to disturbances and thus are essential for regional persistence. As an example, Ditrupa arietina, which displays high densities in the absence of stress, is subject to local extirpations following successive recruitment failures due to mechanical disturbances and resuspension during swell events in Banyuls Bay [START_REF] Charles | Reproduction, recruitment and larval metamorphosis in the serpulid polychaete Ditrupa arietina (O.F. Müller)[END_REF]Labrune et al., 2007b). In the Gulf of Lions, swells impact different parts of the coast according to their origin: southerly swells impact the Languedoc coast whereas easterly swells impact the Roussillon coast (Fig. 2, [START_REF] Guizien | Spatial variability of wave conditions in the Gulf of Lions (NW Mediterranean sea)[END_REF]. The recruitment failure scenario applied to Ditrupa arietina indicates that there is no risk of regional extirpation linked to the easterly swells impacting the Roussillon coast, which is consistent with the recolonization observed after two successive years of recruitment failure in Banyuls Bay (Roussillon coast). Another potential threat to species persistence in the Gulf of Lions is habitat loss due to boat traffic hazards spreading from ports (oil spill, invasive species).

The present study identified source sites essential for metapopulation persistence around Sète. Our findings suggest that a re-evaluation of the boundaries of the Parc Naturel Marin should include the center of the Gulf of Lions and more generally account for connectivity when evaluating the ecological cost of habitat destruction.

This study showed that the few essential sites for regional persistence could not be depicted from the spatial distribution of population density along the gulf because, in the absence of disturbances, connectivity will enable individuals to spread to all importing sites. Consequently, traditional methods based on species inventories may fail to guide the design of protection measures to target source sites. Indeed, for a short-lived species, connectivity patterns lead to a spatial distribution in the Gulf of Lions with a higher relative abundance in the western part than in the eastern part. Such disequilibrated spatial distribution is also present in long-lived species. These findings, which only account for the dispersal process in homogeneous habitat, are consistent with the community diversity throughout the muddy shelf sediment in the Gulf of Lions, which is separated into two clusters corresponding to the eastern and western parts of the gulf (Labrune et al., 2007a). These two clusters were originally differentiated by organic matter content, i.e., a habitat descriptor. In fact, the spatial distribution of organic matter in the gulf is likely to be driven by oceanographic dispersal as well. Such a result highlights the existence of co-variations among factors affecting the spatial distribution of individuals that may blur the assessment of connectivity using statistical inferences [START_REF] Moritz | Disentangling the role of connectivity, environmental filtering, and spatial structure on metacommunity dynamics[END_REF].

In light of our findings, consideration should be given, when setting local protection measures to the oceanographic connectivity at the meso-scale (from tens to hundreds of kilometers) which drives regional distribution of population density and population vulnerability. In the Gulf of Lions, our results indicate that MPAs in the western tip are located in population sink for the soft-bottom species and establishing local protection measures in the central part of the gulf would benefit soft-bottom species over the entire region through the rescue effect. We highlight that metapopulation modeling integrating oceanographic connectivity, local demography and habitat mapping is a powerful tool that enables the inclusion of scenarios in the a priori design of conservation measures. Table 1: Settings for the three simulation groups. Connectivity was set either deterministically using the same matrix every year (either one of the twenty variants or the averaged value D) or stochastically picking a different matrix among the twenty variants every year with a uniform probability in order to conserve the averaged connectivity during the simulation (S). 5: Spatial distribution of the relative population density (ratio between the site population density and the regional average of population densities) of a 2-year lifespan species computed between 91 and 100 years and repeating the same connectivity matrix every year (simulations group 1): (A) for a connectivity matrix displaying a westward drift, (B) for a connectivity matrix displaying an eastward drift, (C) for the average connectivity matrix, and (D) for an ensemble of 500 simulations in which connectivity was picked every year with uniform probability amongst the 20 connectivity variants (simulation group 2). Black circles display the standard deviation of the 500 simulations replicated when connectivity was stochastic. R is the minimum recruitment success value to ensure species persistence in each connectivity case. A B Fig. 7: Median of the average regional coverage of a 2-year lifespan species between 91 and 100 years for the ensemble of 500 simulations of group 3B. Connectivity was picked every year with uniform probability amongst the 20 connectivity variants and the return frequency of recruitment failure in the northern (Languedoc) and western (Roussillon) part of the Gulf of Lions varied from every ten years to every year. Recruitment success R was set to (A) 3%, and (B) 6%. The shaded area is delimited by the 5% and 95% quantiles of the average regional coverage between 91 and 100 years.

Gulf of

7

  Lions (NW Mediterranean Sea) bathymetric map with depth contours displayed in gray lines (depth contours are 20, 50, 100, 500 m and spacing is 500 m below 500 m deep). The nearshore sandy habitat (depth between 10 and 30 m) was meshed into 32 sites from the center of which (indicated by circles and numbered from west to east) larvae were released for dispersal simulations. Hollow circles indicate the release points off the Roussillon coast and filled circles indicate the release points off the Languedoc coast used in the simulations, group 3B. Solid boxes show the extension of the two nested computational domains (1500 m and 750 m resolution grids) and the thick-lined box displays the extent of the Parc Naturel Marin du Golfe du Lion (PNM). The two stars indicate the location of the two other marine protected areas (Parc Marin de la Côte Bleue in the West and Réserve Naturelle Marine de Cerbère-Banyuls in the East). The four trade ports of the Gulf of Lions for the simulations group 3A are indicated by a filled square and a letter. From east to west, M stands for Marseille, S for Sète, PLN for Port La Nouvelle and PV for Port-Vendres. Banyuls Bay (not depicted) is located 10 km south of Port-Vendres. . . . . . . . . . 35 3 (A) Mean and (B) standard deviation of the twenty connectivity matrices computed for twenty 10-day spawning periods and a 4-week pelagic larval duration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36Average regional coverage of a 2-year lifespan species between 91 and 100 years for the ensemble of 500 simulations of group 3A. Connectivity was picked every year with uniform probability amongst the 20 connectivity variants and the number of habitat sites destroyed was increased from 1 to 32.Recruitment success R was set to (A) 3% and (B) 6%. The errorbars displaying variability among the 500 simulations are smaller than the symbols. 38 Median of the average regional coverage of a 2-year lifespan species between 91 and 100 years for the ensemble of 500 simulations of group 3B. Connectivity was picked every year with uniform probability amongst the 20 connectivity variants and the return frequency of recruitment failure in the northern (Languedoc) and western (Roussillon) part of the Gulf of Lions varied from every ten years to every year. Recruitment success R was set to (A) 3%, and (B) 6%. The shaded area is delimited by the 5% and 95% quantiles of the average regional coverage between 91 and 100 years. .

Fig. 1 :Fig. 2 :Fig. 4 :

 124 Fig. 1: (A) Schematic diagram for a spatially explicit model of a metapopulation, where gray intensity indicates adult abundance per unit area in each population and solid arrows indicate larval fluxes in the absence of disturbance. (B) Habitat loss disturbance: adult abundance in dashed sites is set to zero. (C) Recruitment failure disturbance: dashed arrows indicate useless larval transfers due to recruitment failure in sites diplayed in black.

  Fig.5: Spatial distribution of the relative population density (ratio between the site population density and the regional average of population densities) of a 2-year lifespan species computed between 91 and 100 years and repeating the same connectivity matrix every year (simulations group 1): (A) for a connectivity matrix displaying a westward drift, (B) for a connectivity matrix displaying an eastward drift, (C) for the average connectivity matrix, and (D) for an ensemble of 500 simulations in which connectivity was picked every year with uniform probability amongst the 20 connectivity variants (simulation group 2). Black circles display the standard deviation of the 500 simulations replicated when connectivity was stochastic. R is the minimum recruitment success value to ensure species persistence in each connectivity case.

Fig. 6 :

 6 Fig.6: Average regional coverage of a 2-year lifespan species between 91 and 100 years for the ensemble of 500 simulations of group 3A. Connectivity was picked every year with uniform probability amongst the 20 connectivity variants and the number of habitat sites destroyed was increased from 1 to 32. Recruitment success R was set to (A) 3% and (B) 6%. The errorbars displaying variability among the 500 simulations are smaller than the symbols.

  multiplied by the ratio of the areas of origin A i to settling A j sites, b j [adult per larva] is the recruitment success in site j, s jj [no units] is the adult survivorship rate at settling site j, δ ij is the Kronecker symbol (1 when i = j, 0 otherwise) and N max is the site's carrying capacity. The propagule production rate is equal to the number of larvae produced by each adult female F S R f where F is the fecundity rate, S R is the sex ratio in the adult population, and f is the probability of an egg being fertilized. The adult survivorship rate can be related to species life expectancy L E as s jj = exp(ln(0.01) ∆t/L E ), where life expectancy is the age at which 99% of the individuals of the same generation have died.

	i ∈ [1, p] at time t, a i [larva per adult] is the propagule production rate in site i, C ji
	[no units] is the propagule transfer rate from site i (column) to site j (row)

Acknowledgments: We thank METEO FRANCE for supplying atmospheric forcing data and the Banque Hydro for supplying the rivers discharge data for the hydrodynamical

Impact of an invasive species Crepidula fornicata on the hydrodynamics and transport properties of the benthic layer. Aquatic Living Resources. 20(1), 15-31. Mouquet, N., Loreau, M. (2002) Coexistence in Metacommunities: The Regional Similarity Hypothesis. The American Naturalist 159(4), 420-426.

Olds, A.D., Connolly, R.M., Pitt, K.A., Maxwell, P.S., Aswani, S., Albert, S. ( 2014


Mariana Padrón

Federica Costantini

Sandra Baksay

Lorenzo Bramanti

Katell

Federica Costantini

Katell Guizien
email: guizien@obs-banyuls.fr

Passive larval transport explains recent gene flow in a Mediterranean gorgonian

Keywords: marine connectivity, Pelagic Larval Duration, genetic structure, Marine Protected Areas, hydrodynamic model, Gulf of Lion

à 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

The Strategic Plan for Biodiversity 2011-2020 aims to conserve 10% of the ocean through well-connected systems of marine protected areas (MPAs) by 2020 (Aichi Target 11, https://www.cbd.int/sp/targets/rationale/target-11/). Hence, understanding the degree and patterns of connectivity among existing MPAs and areas where protection could be extended becomes urgent. Connectivity is also known as the mechanism through which populations of sessile marine species, and particularly coral species, can recover after local and global disturbances [START_REF] Roberts | Connectivity and management of Caribbean coral reefs[END_REF][START_REF] Hanski | Metapopulation dynamics[END_REF][START_REF] Cowen | Connectivity of marine populations: open or closed?[END_REF]. In the context of global change, relating contemporary connectivity to the seascape processes that shape the spatial distribution of marine species is essential for the design of effective biodiversity conservation strategies. Seascape processes that drive population connectivity of sessile species with a larval dispersal stage, such as corals or gorgonians, should include both seabed geomorphology and dynamic hydrography [START_REF] Manderson | Seascapes are not landscapes: an analysis performed using Bernhard Riemann's rules[END_REF]). Correlations at the local scale between demographic population descriptors and the environment may, indeed, not be sufficient, and the successful exchange of individuals between distant locations (i.e., connectivity) should be included. Recent studies on seascape connectivity have focused on evaluating how larval transport drives genetic patterns across the seascape using a combination of biophysical and genetic approaches (see [START_REF] Liggins | Taking the plunge: an introduction to undertaking seascape genetic studies and using biophysical models[END_REF] for a review). Nonetheless, [START_REF] Selkoe | A decade of seascape genetics: contributions to basic and applied marine connectivity[END_REF] reported that only 31% of the analyzed studies found ocean currents to be a good predictor of gene flow patterns and ascribed this poor predictability to inappropriate scales in genetic sampling. These findings suggest that demographic and genetic connectivity of benthic sessile species, for which connectivity among populations only occurs during the larval dispersive stage, are affected by several processes acting as successive biological filters at different temporal and spatial scales along their life cycles from spawning to reproduction by the settled individuals [START_REF] Pineda | Larval transport and dispersal in the coastal ocean and consequences for population connectivity[END_REF].

Prior to dispersal, population size and fecundity are among the key factors that determine the number of emigrant offspring and have been shown to influence connectivity patterns in marine species [START_REF] Treml | Reproductive output and duration of the pelagic larval stage determine seascape-wide connectivity of marine populations[END_REF][START_REF] Dawson | Dispersal potential and population genetic structure in the marine intertidal of the eastern North Pacific[END_REF]. After dispersal, gene flow among areas that are well connected through larval transport can be limited by settlement regulation due to space limitations [START_REF] Roughgarden | Demographic theory for an open marine population with space-limited recruitment[END_REF]. However, the primary filter of connectivity is larval transport, which is notably driven by the amount of time that larvae spend being transported by ocean currents, otherwise known as the pelagic larval duration (PLD). Hence, the PLD has been the primary factor considered in biophysical modeling when attempting to differentiate dispersal distance capabilities among different species [START_REF] Shanks | Propagule dispersal distance and the size and spacing of marine reserves[END_REF].

Genetic studies have tested the relationship between realized connectivity (gene flow) and the PLD following the same reasoning, and the results have led to contrasting conclusions [START_REF] Kinlan | Propagule dispersal in marine and terrestrial environments: a community perspective[END_REF][START_REF] Weersing | Population genetics, larval dispersal, and connectivity in marine systems[END_REF][START_REF] Riginos | Effects of geography and life history traits on genetic differentiation in benthic marine fishes[END_REF][START_REF] D'aloia | Patterns, causes, and consequences of marine larval dispersal[END_REF].

Several arguments have been invoked to explain the absence of a relationship between PLD and realized dispersal distance and migration rates, with some arguments referring to the way in which realized connectivity is inferred from measures of genetic differentiation among populations. For instance, the inverse relationship between the most common metric of genetic differentiation, FST, and the migration rate (Wright's Island Model 1931), may be blurred by the uncertainty of the estimation of the effective size of the population [START_REF] Faurby | Theoretical limits to the correlation between pelagic larval duration and population genetic structure[END_REF].

However, the failure to identify a relationship between PLD and realized connectivity is more likely explained by the lack of unequivocal correlation between PLD and dispersal distance. In fact, the variability of ocean currents [START_REF] Cowen | Scaling of connectivity in marine populations[END_REF][START_REF] Guizien | Vulnerability of marine benthic metapopulations: implications of spatially structured connectivity for conservation practice in the Gulf of Lions (NW Mediterranean Sea)[END_REF], larval behavior [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][START_REF] Paris | Surfing, spinning, or diving from reef to reef: effects on population connectivity[END_REF]) and larval mortality [START_REF] Cowen | Connectivity of marine populations: open or closed?[END_REF] can lead to very different dispersal patterns given the same PLD. Moreover, for any given species, PLD does not have a single value, but it varies along a competency window [START_REF] Arnold | Setting into an increasingly hostile world: the rapidly closing "recruitment window" for corals[END_REF], within which successful recruitment can occur at any time once a suitable habitat is found along the dispersal pathway.

To determine which biological filters transform larval release into gene flow (potential into realized connectivity), we used one of the most abundant and conspicuous species of the rocky sublittoral zone in the northwestern Mediterranean Sea, the white gorgonian Eunicella singularis.

In the hydrodynamic province of the Gulf of Lion (NW Mediterranean, [START_REF] Rossi | Hydrodynamic provinces and oceanic connectivity from a transport network help designing marine reserves[END_REF], the spatial distribution of E. singularis is discontinuous [START_REF] Gori | Spatial distribution patterns of the gorgonians Eunicella singularis, Paramuricea clavata, and Leptogorgia sarmentosa (Cape of Creus, Northwestern Mediterranean Sea)[END_REF] due to fragmentation of the rocky habitat [START_REF] Aloisi | Carte géologique du précontinent languedocien au 1/250000ième[END_REF]). E. singularis is believed to play an important ecological role due to its erect structure, which provides habitat for the local epifauna and increases the biomass and diversity of the rocky community [START_REF] Ballesteros | Mediterranean coralligenous assemblages: a synthesis of present knowledge[END_REF][START_REF] Mitchell | Secondary production of gorgonian corals in the northern Gulf of Mexico[END_REF]. Furthermore, its indirect economic value due to its contributions to shaping the seascape and attracting SCUBA diving and fishing activities makes E. singularis a species of interest in the monitoring plans of all MPAs in the Gulf of Lion. E. singularis is a long-lived (up to 15 years) internal brooding sessile species. Lecitotrophic larvae (planulae) are released once a year in June-July [START_REF] Théodor | Contribution à l'étude des gorgones. VII. Ecologie et comportement de la planula[END_REF][START_REF] Weinberg | The life cycle of a gorgonian: Eunicella singularis[END_REF]. Larvae can survive up to 122 days in ex situ experiments in the absence of predation [START_REF] Théodor | Contribution à l'étude des gorgones. VII. Ecologie et comportement de la planula[END_REF], with 50% to 80% survival after 40 days (K. Guizien, Personal communication). Competency for recruitment was observed after a few days and lasted up to a few weeks (L. Bramanti, Personal communication). Such larval traits suggest a wide dispersive potential of the species.

In the present study, we propose a methodological approach that combines demographic, biophysical and genetic methods to test the relationship between the dispersive potential of a sessile benthic species and the effective gene flow in a fragmented habitat layout. The study aims to (1) assess the genetic structure of Eunicella singularis within the fragmented habitat of the Gulf of Lion, (2) test the role of demographic filters (i.e., reproductive output and intra-specific recruitment limitations) in shaping gene flow, and (3) test the significance of various PLD-driven larval transport modes in shaping the spatial pattern of gene flow.

MATERIAL AND METHODS

Study area and regional hydrodynamics

The Gulf of Lion (hereafter GoL) is a large (~100 km), microtidal continental shelf with a 100-km radius located in the northwestern Mediterranean Sea (Figure 1). The first 30 meters of the bathymetry is mainly covered by soft sandy sediment; however, fragmented patches of hard substrate can be found in the eastern (Côte Bleue), central (Plateau des Aresquiers and Cap d'Agde), and western (Cap Leucate and Côte Vermeille) parts of the gulf [START_REF] Aloisi | Carte géologique du précontinent languedocien au 1/250000ième[END_REF].

Circulation in the GoL is driven by three main forces: 1) strong continental winds from the north (Mistral) and the northwest (Tramontane) that frequently generate upwelling and downwelling along the coast, as well as transient currents and eddies; 2) the Northern Current (NC), which usually flows westward along the shelf break of the GoL from the Ligurian Sea to the Catalan Sea; 3) river inputs, particularly from the Rhône, that deliver freshwater and high nutrient and organic matter loads [START_REF] Millot | The Gulf of Lions' hydrodynamics[END_REF][START_REF] Hu | Numerical simulations of eddies in the Gulf of Lion[END_REF][START_REF] Campbell | Nutrients and plankton spatial distributions induced by coastal eddy in the Gulf of Lion. Insights from a numerical model[END_REF].

Numerical simulations at the basin scale have identified that the GoL is an isolated hydrodynamic province, with the Northern Current acting as a hydrodynamic barrier that limits the transport of particles across the shelf break [START_REF] Rossi | Hydrodynamic provinces and oceanic connectivity from a transport network help designing marine reserves[END_REF]. Large-scale flow limits entries into the Gulf of Lion up to a 60-day PLD from the east and up to a 30-day PLD from the west. Moreover, despite an overall westward drift leading to large sediment exports through western canyons [START_REF] Ulses | Suspended sediment transport in the Gulf of Lions (NW Mediterranean): Impact of extreme storms and floods[END_REF], persistent anticyclonic eddies are observed in the southwestern part of the GoL during the summer [START_REF] Hu | Numerical simulations of eddies in the Gulf of Lion[END_REF], and connectivity patterns within the Gulf of Lion stabilize after 3 weeks.

Sampling design

From July 2013 to October 2014, detectable colonies (taller than 1 cm) of Eunicella singularis were counted and sampled at 17 geo-referenced stations distributed along the GoL, from the Côte Bleue to the Côte Vermeille (Table 1) (Figure 1). The stations were distributed along four rocky habitat patches separated by large areas of sandy bottom (Côte Bleue, Plateau des Aresquiers, Cap d'Agde and Côte Vermeille) at depths between 15 and 35 meters. The distance between stations varied from less than 1 km to 180 km. The population size structures were similar in the four rocky habitat patches, with 25% of the colonies smaller than 5 cm, 55% of the colonies between 5 and 20 cm and 20% of the colonies taller than 20 cm (L. Bramanti,

Personal communication).

The number of stations sampled in each habitat patch varied according to the along-shore extension of E. singularis coverage of the rocky bottom (L. Bramanti,

Personal communication).

In Côte Bleue, where the smallest coverage of E. singularis was found, only three stations were sampled; in Plateau des Aresquiers and Cap d'Agde, four stations were sampled, while six stations were sampled in the Côte Vermeille.

The mean population density at each station was calculated from colony counts in 4 quadrats (1 m 2 each). For each station, the 4 quadrats were located 5 m from the geo-referenced station, one in each of the four cardinal directions. Small apical fragments (<3 cm) of all colonies of E. singularis that were counted at each station were collected, preserved in 95% ethanol and kept at 6 °C until DNA extraction. For the present study, 13 to 36 colonies per station were randomly selected and genotyped, for a total of 434 colonies.

Genetic diversity and population structure

All colonies were genotyped at eight microsatellite loci specifically developed for E. singularis: C21, C30 and C40 (Molecular Ecology Resources Primer Development Consortium 2010) and E. verrucosa: Ever003, Ever004, Ever009, Ever013 and Ever014 [START_REF] Holland | Isolation and characterization of fourteen microsatellite loci from the endangered octocoral Eunicella verrucosa (Pallas 1766)[END_REF] after DNA extraction (Appendix S1). The scoring of the alleles was carried out using GeneMapper v.3.7 software (Applied Biosystems). We calculated the summary statistics of the genetic diversity (allelic richness, observed heterozygosity, expected heterozygosity), and tested for departure from the Hardy-Weinberg equilibrium, linkage disequilibrium and the presence of null alleles (Appendix S1). The statistical power to detect the level of genetic variation among populations was estimated using Powsim v.4.1 [START_REF] Ryman | POWSIM: a computer program for assessing statistical power when testing for genetic differentiation[END_REF]. The spatial genetic structure of E. singularis was analyzed by testing for the presence of a pattern using the isolation by distance (IBD) and the Loiselle kinship coefficient (Appendix S1). The patterns of population genetic structure were analyzed using FST, a Bayesian clustering method that identifies clusters of individuals who share similar patterns of variation (Structure v.2.3.4;[START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF] and a hierarchical analysis of molecular variance (AMOVA). For details, see the Methods in Appendix S1.

Contemporary gene flow

The genotyped individuals from the 17 stations were clustered into 4 groups defined by habitat discontinuity, corresponding to the four largest rocky habitat patches (Côte Vermeille, Cap d'Agde, Plateau des Aresquiers and Côte Bleue). Contemporary gene flow over the last 2-3 generations was estimated in the form of matrices containing migration rates including a proxy of self-recruitment among these 4 clusters using BayesAss 3.0 applied on multilocus genotypes (Wilson and Rannala 2013). It is important to note that the BayesAss 3.0 output is the proportion of individuals in the destination population that is descendant from a source population.

The confidence intervals and the uncertainty values around the average migration rates were estimated (Appendix S1).

Recent larval transport

Using the 3D ocean model Symphonie (S-2010 release 26, at http://sirocco2.omp.obsmip.fr/outils/Symphonie/Accueil/SymphoAccueil.htm), high-resolution flow simulations were carried out and used for larval dispersal simulations (Appendix S1). The latter were upscaled to the larval transport matrices for single reproductive events. Larval transport matrices contain transfer rate values defined as the proportion of larvae (neutrally buoyant particles) released from a source habitat patch (row) during a given release period that reach a destination habitat patch (column) after a given PLD. Fifteen variations of larval transport matrices were established per PLD, with each variant corresponding to a non-overlapping weeklong release period in June (five variants per year) over three years (2010( , 2011( , 2012, Appendix S1), Appendix S1). Six weeks of dispersal were simulated and eight PLD values ranging from 3.5 to 42 days were tested. The PLD range that was tested was greater than the period during which Eunicella singularis larvae competency for recruitment as observed in laboratory experiments (L. Bramanti, Personal Communication).

Moreover, the maximum PLD tested in the present study (42 days) was longer than the time required to observe stable connectivity patterns within the Gulf of Lion (Briton et al., submitted).

Comparison between larval transport and gene flow

and the mean population density in that source (proxy of reproductive output given the same population size structure in all habitat patches), and 2) the number of contemporary immigrants to a destination habitat patch (sum of migrant numbers by column) and the mean population density in that destination patch (proxy of intra-specific competition). All statistical analyses were conducted in Matlab (R2012a).

RESULTS

Genetic diversity and population structure

All sampled individuals were included in the analysis (no identical multilocus genotypes). All loci were polymorphic, with more than 4 alleles among stations on average and from 5 to 13 alleles over all loci, which is comparable to the degree of polymorphism obtained by [START_REF] Holland | Isolation and characterization of fourteen microsatellite loci from the endangered octocoral Eunicella verrucosa (Pallas 1766)[END_REF]. Linkage disequilibrium was not observed (Appendix S2, Table A2). The observed and expected heterozygosities were similar among stations with no significant heterozygote excesses, but there were significant heterozygote deficits across two populations (CV1 and CA1) (Appendix S2, Table A2) and in two of the eight loci (Ever003 and C40) (Appendix S2, Table A3). The frequency of null alleles was negligible (Appendix S2, Table A3).

The overall FST value was low (0.04), and almost all pairwise FST values among stations were significantly different (P-value < 0.001) (Appendix S2, Table A4). Nei's genetic distance values were closer among the stations within habitat patches than between them (Figure 2 A5).

The correlation between FST / (1-FST) and the logarithm of the geographic distance (R 2 = 0.28; P-value = 0.01) was significant but low among the 17 stations in the GoL and was not significant among the 4 patches defined according to rocky habitat fragmentation (R 2 = 0.81; Pvalue = 0.100) (Appendix S2, Figure A2). The estimated maximum size of the genetic patch based on Loiselle's kinship coefficient was 30 km, which is larger than the maximum size of any of the habitat patches in the Gulf of Lion (Appendix S2, Figure A3).

Contemporary gene flow

Figure 4 shows the spatial patterns of the number of contemporary migrants among the 4 separated habitat patches in the GoL. The uncertainty of the migration pattern, following an error metric adapted from the Euclidean matrix norm (see details in Appendix S1), was below 25%.

Local retention dominated over emigration in all habitat patches. However, local retention was much lower in Côte Bleue (mean = 51 migrants, CI (95%) = 48 -52 migrants) compared to all other habitat patches (from 90 to 121 migrants) and was comparable to the value of emigration from Cap d'Agde (45 ± 14 migrants). Cap d'Agde acted as a hub, distributing migrants to both Côte Vermeille (westward) (mean = 16 migrants, CI (95%) = 6-19 migrants) and Plateau des Aresquiers (eastward) (mean = 28 migrants, CI (95%) = 17-34 migrants).

Relationship between population density and migration

Neither the number of contemporary emigrants coming from a habitat patch (Figure 5A) nor the number of contemporary immigrants arriving to a habitat patch (Figure 5B) showed a significant correlation with the population density of the source (R 2 = 0.203, P > 0.05) or the destination (R 2 = 0.787, P > 0.05) habitat patch, respectively.

Although the lowest number of emigrants corresponded to the source habitat patch with the lowest population density, the opposite was not true. The highest number of emigrants did not come from the habitat patch with the highest population density. The number of immigrants did not decay with increasing population density in the destination habitat patch but conversely increased. However, a similar number of immigrants were found in the destination habitat patches with quite different population densities.

Comparison between recent larval transport and contemporary gene flow

In the following section, both larval transport and gene flow have been normalized for each source population and are subsequently no longer comparable between different source populations. Recent larval transport varied greatly among the PLDs ranging from 3.5 to 42 days: from a scenario of high local retention within all habitat patches to a scenario of high levels of larval transport among most habitat patches, except Côte Vermeille (Figure 6A Vermeille. This pattern matched the one observed for contemporary gene flow particularly well (Figure 6F).

The Mantel test showed a significant correlation (Mantel statistics <5%) between contemporary gene flow and recent larval transport (excluding retention rates) only for the PLDs ranging from 7 to 14 days with a determination coefficient (including retention rates) larger than 50% (Figure 7). However, the determination coefficient between contemporary gene flow and recent larval transport, which includes retention rates, was also high for a PLD of 3.5 days (R 2 >0.95). Interestingly, the determination coefficient did not decrease monotonically with PLD but rather peaked at a PLD of 7 days (R 2 = 0.98).

DISCUSSION

Habitat fragmentation fails to explain regional genetic structure of Eunicella singularis in the

Gulf of Lion

The analysis of the genetic differentiation of groups of Eunicella singularis colonies among a wide range of spatial scales (from less than 1 km to 180 km) showed high levels of genetic variability at small scales with a deficit of heterozygotes at some stations, as reported previously for the same species [START_REF] Costantini | Limited genetic connectivity between gorgonian morphotypes along a depth gradient[END_REF] and in other Mediterranean gorgonians such as Paramuricea clavata [START_REF] Mokhtar-Jamaï | From global to local genetic structuring in the red gorgonian Paramuricea clavata: the interplay between oceanographic conditions and limited larval dispersal[END_REF] and Corallium rubrum [START_REF] Costantini | Genetic structuring of the temperate gorgonian coral (Corallium rubrum) across the western Mediterranean Sea revealed by microsatellites and nuclear sequences[END_REF]. In these studies, genetic patchiness was attributed to inbreeding, putatively due to the low dispersal ability of the larvae. However, the genetic clustering analysis of E. singularis populations in the GoL revealed strong similarity among the individuals sampled in Cap d'Agde and Plateau des Aresquiers, which are separated by 30 km of sandy bottom that is unsuitable for colonization by the species [START_REF] Weinberg | Mediterranean octocorallian communities and the abiotic environment[END_REF]. Moreover, one station in the Côte Vermeille cluster displayed higher similarity to the stations at Cap d'Agde and Plateau des Aresquiers than to the other stations in Côte Vermeille. These similarities to the genetic profiles of distant stations led to an estimated genetic patch size of 30 km, meaning that the individuals found within an area of 30 km would be more genetically related than those taken randomly from any further distance. These results, together with the absence of genetic structure among E. singularis populations found within 15 km [START_REF] Costantini | Limited genetic connectivity between gorgonian morphotypes along a depth gradient[END_REF] in the neighboring hydrodynamic province (Balearic Sea, [START_REF] Rossi | Hydrodynamic provinces and oceanic connectivity from a transport network help designing marine reserves[END_REF], suggest that the heterozygote deficit in the E. singularis populations in the GoL should be attributed to on-going migration among open populations rather than to inbreeding in closed populations. Nevertheless, migration was not explained by isolation by distance [START_REF] Sexton | Genetic isolation by environment or distance: which pattern of gene flow is most common?[END_REF][START_REF] Nanninga | Environmental gradients predict the genetic population structure of a coral reef fish in the Red Sea[END_REF][START_REF] Thomas | Isolation by resistance across a complex coral reef seascape[END_REF], which pinpoints the limitation of using the IBD model to explain the chaotic genetic patchiness resulting from larval dispersal via ocean flow [START_REF] Selkoe | Taking the chaos out of genetic patchiness: seascape genetics reveals ecological and oceanographic drivers of genetic patterns in three temperate reef species[END_REF]. Indeed, undirected pairwise FST is not a good measure of genetic distance when migration among populations is asymmetrical [START_REF] Beerli | Estimation of migration rates and population sizes in geographically structured populations[END_REF][START_REF] Broquet | Quantifying genotyping errors in noninvasive population genetics[END_REF][START_REF] Van Strien | Isolation-by-distance in landscapes: considerations for landscape genetics[END_REF]. In such a case, assignment tests or coalescent methods should be used to describe gene flow patterns and determine the factors shaping them [START_REF] Manel | Assignment methods: matching biological questions with appropriate techniques[END_REF].

Gene flow is determined by larval transport

We tested three major successive filters that potentially influence gene flow during larval dispersal among habitat patches. First, the number of emigrants of E. singularis within separated habitat patches in the GoL was not explained by reproductive output. Second, the number of immigrants was not limited by settlement regulation due to intra-specific competition, although this mechanism is often presented [START_REF] Roughgarden | Demographic theory for an open marine population with space-limited recruitment[END_REF]Padron and Guizien 2015). It is interesting to note that conversely, the population density in the destination habitat patch increased with the number of immigrants, suggesting that the population densities observed in the GoL were below habitat saturating capacity. Third, larval transport explained a large portion of the gene flow among populations. In the GoL, the determination coefficient between contemporary gene flow and recent larval flow was extremely high for a 7-day larval transport duration among four habitat patches.

The combination of the absence of a density-dependent regulation of E. singularis migrants in the GoL and the fact that gene flow is driven by unsteady larval transport support the alternate dominance of recruits vs. adult colony populations hypothesis that was observed for the same species along the Costa Brava (NW Mediterranean, Spain) [START_REF] Linares | Size distribution, density and disturbance in two Mediterranean gorgonians: Paramuricea clavata and Eunicella singularis[END_REF]. The strong correlation between recent larval transport and contemporary gene flow found in this study contrast with other studies based on the correlation with multi-generational gene flow (White et al. 2010;[START_REF] Alberto | Isolation by oceanographic distance explains genetic structure for Macrocystis pyrifera in the Santa Barbara Channel[END_REF]. This underpins the importance of testing the relationship between larval transport and gene flow estimated for a single dispersive event rather than for multiple generations. Moreover, estimations at the single dispersive event scale enabled the discrimination of the efficiency of different PLDs in connecting populations.

The environment filters the PLDs that connect populations

Ex situ observations of survival expectancy and competency duration of E. singularis larvae (up to 40 days, K. Guizien and L. Bramanti, Personal Communication) suggest that all tested PLDs were equally probable in the absence of predation. In such a case, one would expect gene flow to be a composite of larval transport for all PLDs. However, gene flow among populations was best explained by PLDs ranging from 7 to 14 days in the habitat configuration of the GoL, and led to a spatially structured genetic population despite high local chaotic genetic patchiness. Similar finding was obtained among lobster populations along the Californian coast and related to upwelling intensity [START_REF] Iacchei | Combined analyses of kinship and FST suggest potential drivers of chaotic genetic patchiness in high gene-flow populations[END_REF]). In the present case of a lecithotrophic larva with a wide competence window, we suggest that larval mortality significantly altered larval transport patterns for periods longer than 14 days [START_REF] Cowen | Connectivity of marine populations: open or closed?[END_REF]. Correlatively, it suggests that mortality, particularly by predation, on Eunicella singularis larvae (2-3 mm) can be neglected up to a 14-day PLD in the GoL. [START_REF] Treml | Reproductive output and duration of the pelagic larval stage determine seascape-wide connectivity of marine populations[END_REF] found that longer maximum PLDs lead to a more connected seascape on an ocean basin scale. The study claimed that maximum PLD is a good predictor of larval transport at large geographic scales (200-500 km) but a poor predictor at a mesoscale (<200 km).

Our study evaluated the larval transport patterns among populations at a regional scale (up to 180 km) and showed that the PLD-indexed larval transport could significantly determine gene flow at a mesoscale. Our findings suggest that several PLDs should be tested along the competency period of a species rather than the average or the maximum PLD to provide a more accurate assessment of the importance of larval transport in shaping gene flow among sessile benthic species [START_REF] D'aloia | Patterns, causes, and consequences of marine larval dispersal[END_REF].

Implications for conservation

Our results have important implications for regional conservation and the spatial management of benthic sessile species in the GoL. Empirical gene flow patterns allow for the identification of key habitat patches of E. singularis: the high levels of gene flow between Cap d'Agde and Plateau des Aresquiers, along with transfers to other habitat patches, suggest that these two areas play a pivotal role in the spread of E. singularis in the region. Furthermore, the significant correlation between contemporary gene flow and recent larval transport for a PLD ranging from 7 to 14 days suggests that Cap d'Agde and Plateau des Aresquiers could be important sites for any species with similar dispersive traits (PLD, larval behavior, or reproductive period) dwelling on the same rocky substrate, as they should follow the same potential connectivity patterns. Hence, we advocate for enforcing the protection of those sites already identified as important for the regional persistence of soft-bottom species [START_REF] Guizien | Vulnerability of marine benthic metapopulations: implications of spatially structured connectivity for conservation practice in the Gulf of Lions (NW Mediterranean Sea)[END_REF], which also appears to be functionally important for the life cycle of some rocky species.

In summary, the methodological approach presented in this study highlights the hidden determinism within the apparent genetic chaos by incorporating interdisciplinary data (demography and ocean currents) and applying accurate genetic tools (coalescent-based analyses) to evaluate the gene flow among populations. Information for each station is provided in Table 1. Information for each station is provided in Table 1.

  ). While the stations in Côte Bleue and Cap d'Agde-Plateau des Aresquiers are close to each other, the stations in the Côte Vermeille (CV) appear more scattered along the first axis, with one station (CV5) falling near the stations of Cap d'Agde and Plateau des Aresquiers (CA and PA, respectively). The Bayesian clustering analysis identified three genetic clusters (K = 3) (mean LnP(K) = -7532.99; Delta K = 25.09). The resulting clustering partially corresponded to the rocky habitat fragmentation in the region (Figure 3): Côte Vermeille (6 stations), Cap d'Agde and Plateau des Aresquiers (8 stations), and Côte Bleue (3 stations). Some individuals in Côte Vermeille depicted higher similarity to individuals from Cap d'Agde and Plateau des Aresquiers rather than to other individuals in Côte Vermeille. The results of the AMOVA confirmed the clustering structure (P < 0.001) (Appendix S2, Table

  to 6E). Larval transport was restricted to the exchange between Cap d'Agde and Plateau des Aresquiers for the shortest PLD tested (3.5 days, Figure 6A), and increased in a westward direction for the PLDs longer than 14 days (Figures 6D, 6E and data not shown). For the PLDs of 7 and 10.5 days (Figure 6B and 6C), larval transport was westward from Plateau des Aresquiers and Côte Bleue, but both westward and eastward from Cap d'Agde. Moreover, it is interesting to note that the local retention decayed more rapidly in Cap d'Agde compared to the other habitat patches due to efficient larval export towards the closest habitat patches of Plateau des Aresquiers and Côte

  ., L.B., and K.G. designed the study. S.B. performed all the laboratory work. M.P., F.C. and K.G. analyzed all the data. All authors contributed to the writing of the manuscript.

Figure 1 .-

 1 Figure 1.-Locations of sampling stations and habitat patches of Eunicella singularis in the Gulf of Lion, NW Mediterranean. A) Black rectangles highlight the sampled habitat patches. Red areas show the distribution of the rocky habitat in the region. B) Sampling stations of E. singularis and rocky habitat extent in the Côte Bleue. C) Sampling stations of E. singularis and rocky habitat extent in the Côte Vermeille. D) Sampling stations of E. singularis and rocky habitat extent in Cap d'Agde and Plateau des Aresquiers. Red squares correspond to the areas of larval release considered in the hydrological model. Sampling stations are represented by gray dots.

Figure 2 .

 2 Figure 2.-Plot of the Principal Coordinate Analysis (PCoA) from the genetic distance matrix assessed with pairwise Nei's genetic distance estimates of E. singularis in the Gulf of Lion. The first two axes explain 61.51% of the variation.

Figure 3 .-

 3 Figure 3.-Population structure of E. singularis in the Gulf of Lion, as revealed by Structure software. Each individual is represented by a vertical color line, with each color representing the proportion of membership to each cluster. Black lines represent the division among the identified clusters (K = 3). Values between parentheses correspond to the overall proportion of membership of the individuals in each cluster.

Figure 4 .-

 4 Figure 4.-Patterns of contemporary migration among habitat patches of Eunicella singularis in the Gulf of Lion. The matrix shows the number of migrants coming from any potential source and reaching any destination habitat patch. CV = Côte Vermeille, CA = Cap d'Agde, PA = Plateau des Aresquiers, and CB = Côte Bleue. The number of migrants between any two habitat patches is represented by the grayscale intensity.

Figure 5 .

 5 Figure 5.-A) Relationship between the number of emigrants (sum of number of migrants by row) from a source habitat patch and the population density in that source habitat patch. B) Relationship between the number of immigrants (sum of number of migrants by column) to a destination habitat patch and the population density in that destination habitat patch. X error bars correspond to the standard deviation around the mean population density values. Y error bars correspond to the confidence intervals around the number of migrants estimated by BayesAss.

Figure 6 .-

 6 Figure 6.-Patterns of recent larval transport and contemporary gene flow among habitat patches of E. singularis in the Gulf of Lion. A) Recent larval transport for a PLD of 3.5 days; B) recent larval transport for a PLD of 7 days; C) recent larval transport for a PLD of 10.5 days; D) recent larval transport for a PLD of 14 days; E) recent larval transport for a PLD of 21 days; F) contemporary gene flow. The Grayscale indicates the probability values.

Figure 7 .-

 7 Figure 7.-Determination coefficient between contemporary gene flow and recent larval transport (including retention rates) for a PLD ranging from 3.5 to 42 days. Filled circles depict the PLDs for which the Mantel test showed significant correlations between contemporary gene flow and recent larval transport matrices, excluding retention rates.

Figure 1 .

 1 Figure 1.-Locations of sampling stations and habitat patches of Eunicella singularis in the Gulf of Lion, NW Mediterranean. A) Black rectangles highlight the sampled habitat patches. Red areas show the distribution of the rocky habitat in the region. B) Sampling stations of E. singularis and rocky habitat extent in the Côte Bleue. C) Sampling stations of E. singularis and rocky habitat extent in the Côte Vermeille. D) Sampling stations of E. singularis and rocky habitat extent in Cap d'Agde and Plateau des Aresquiers. Red squares correspond to the areas of larval release considered in the hydrological model. Sampling stations are represented by gray dots.

Figure 2 .-

 2 Figure 2.-Plot of the Principal Coordinate Analysis (PCoA) from the genetic distance matrix assessed with pairwise Nei's genetic distance estimates of E. singularis in the Gulf of Lion. The first two axes explain 61.51% of the variation.

Figure 3 .

 3 Figure 3.-Population structure of E. singularis in the Gulf of Lion, as revealed by Structure software. Each individual is represented by a vertical color line, with each color representing the proportion of membership to each cluster. Black lines represent the division among the identified clusters (K = 3). Values between parentheses correspond to the overall proportion of membership of the individuals in each cluster.

Figure 4 .-

 4 Figure 4.-Patterns of contemporary migration among habitat patches of Eunicella singularis in the Gulf of Lion. The matrix shows the number of migrants coming from any potential source and reaching any destination habitat patch. CV = Côte Vermeille, CA = Cap d'Agde, PA = Plateau des Aresquiers, and CB = Côte Bleue. The number of migrants between any two habitat patches is represented by the grayscale intensity.

Figure 5 .

 5 Figure 5.-A) Relationship between the number of emigrants (sum of number of migrants by row) from a source habitat patch and the population density in that source habitat patch. B) Relationship between the number of immigrants (sum of number of migrants by column) to a destination habitat patch and the population density in that destination habitat patch. X error bars correspond to the standard deviation around the mean population density values. Y error bars correspond to the confidence intervals around the number of migrants estimated by BayesAss.

Figure 6 .

 6 Figure 6.-Patterns of recent larval transport and contemporary gene flow among habitat patches of E. singularis in the Gulf of Lion. A) Recent larval transport for a PLD of 3.5 days; B) recent larval transport for a PLD of 7 days; C) recent larval transport for a PLD of 10.5 days; D) recent larval transport for a PLD of 14 days; E) recent larval transport for a PLD of 21 days; F) contemporary gene flow. The Grayscale indicates the probability values.

Figure 7 .-

 7 Figure 7.-Determination coefficient between contemporary gene flow and recent larval transport (including retention rates) for a PLD ranging from 3.5 to 42 days. Filled circles depict the PLDs for which the Mantel test showed significant correlations between contemporary gene flow and recent larval transport matrices, excluding retention rates.

Acknowledgements

This work was (co-) funded through a MARES Grant. MARES is a Joint Doctorate program selected under Erasmus Mundus coordinated by Ghent University (FPA 2011-0016). Check www.-mares-eu.org for additional information. This work was also partially funded by the French National Program LITEAU IV of the Ministère de l'Ecologie et de l'Environnement Durable under project RocConnect-Connectivité des habitats rocheux fragmentés du Golfe du Lion (PI, K. Guizien, Project Number 12-MUTS-LITEAU-1-CDS-013. The authors particularly thank the scientific managers of the Gulf of Lion MPA: S. Blouet, E. Charbonnel, B. Ferrari and J. Payrot for valuable interactions during the study design. The authors gratefully acknowledge the helpful assistance during sampling of R. Bricout, F. Cornette, S. Fanfard, B. Hesse, C. Labrune, L. Lescure, J.-C. Roca, P. Romans, and the staff of the Réserve Naturelle Marine de Cerbère-Banyuls, Aire Marine Protégée Agatoise, Parc Naturel Marin du Golfe du Lion and Parc Marin de la Côte

Bleue. We also acknowledge American Journal Experts for English editing service.

Data Accessibility

Microsatellite genotypes will be made available in Dryad.

With gene flow being estimated as the proportion of individuals in the destination population that is a descendant from a source population, and larval transport being the proportion of the larvae released from a source population that arrives in a destination population, comparing gene flow and larval transport requires a series of data transformations (see details in Appendix S1). Briefly, as the actual number of larvae released from a source population is unknown, the transformation of larval transport into the number of larvae reaching a destination population is not possible. Hence, gene flow was transformed into the proportion of migrants that leaves a source population and arrives at each destination population. For each PLD, recent larval transport was estimated as the average of the 15 larval transport matrix variants (one for each single reproductive event). Recent larval transport matrices were normalized so that the sum of the destination probabilities from the same source (i.e., all columns within each row) was equal to one, which removed the proportion of larvae dispersed out of the 4 habitat patches sampled for genetics as this proportion is not included in the gene flow assessment. It is worth noting that after these normalizations, for both the gene flow and larval transport matrices, only rows, not columns, are statistically independent. Hence, the significance of the Pearson product-moment correlations (excluding diagonal terms, i.e., retention rates) between the contemporary gene flow and recent larval transport matrices (4x4) was tested using a Mantel test [START_REF] Mantel | The Detection of Disease Clustering and a Generalized Regression Approach[END_REF]) for all 24 possible permutations, after modifying the random QAP.m Matlab routine by Puck Rombach (2011). In addition, the determination coefficient between contemporary gene flow and recent larval transport was computed including diagonal terms, i.e., retention rates.

Relationship between population density and migration

Finally, the effect of demographic filters on the number of migrants (see calculation in Appendix S1) was tested, namely, the modulation of emigrants by reproductive output from the source populations and the modulation of immigrants by settlement limitation in the destination populations. Two correlation analyses were performed to test the relationship between 1) the number of contemporary emigrants from a source habitat patch (sum of migrants number by row)