An integrated genetic-demographic model to unravel the origin of genetic structure in European eel (Anguilla anguilla L.).

The evolutionary enlightened management of species with complex life cycles often requires the development of mathematical models integrating demographic and genetic data. The genetic structure of the endangered European eel (Anguilla anguilla L.) has been thoroughly analyzed in several studies in the past years. However, the interpretation of the key demographic and biologic processes that determine the observed spatio-temporal genetic structure has been very challenging owing to the complex life cycle of this catadromous species. Here, we present the first integrated demographic-genetic model applied to the European eel that explicitly accounts for different levels of larval and adult mixing during oceanic migrations and allows us to explore alternative hypotheses on genetic differentiation. Our analyses show that (i) very low levels of mixing occurring during larval dispersal or adult migration are sufficient to erase entirely any genetic differences among sub-populations; (ii) small-scale temporal differentiation in recruitment can arise if the spawning stock is subdivided in distinct reproductive groups; and (iii) the geographic differentiation component might be overestimated if a limited number of temporal recruits are analyzed. Our study can inspire the scientific debate on the interpretation of genetic structure in other species characterized by complex life cycle and long-range migrations.

Introduction

A sustainable and sound management of renewable marine resources requires a comprehensive knowledge of the spatial structure of target populations in time and space (Waples 1998). This is particularly true for threatened species as conservation strategies must account for the connectivity level among subpopulations (Hauser and Carvalho 2008). The integration of molecular techniques with behavioral, demographic, and environmental information has proven to be extremely powerful in identifying even subtle signals of genetic structure in time and/or space (Selkoe et al. 2008). In fact, heterogeneities in the genetic structure of populations may emerge as a consequence of the complex interactions between biologic, demographic, and environmental processes at different spatio-temporal scales and at different stages of the ontogenic development of a species. As a consequence, understanding which are the key processes shaping the observed patterns of genetic structure can be a very challenging task. Moreover, the implementation of these results into management strategies becomes even more complex when the genetic structure deviates from a straightforward geographic pattern (as evidenced by Moritz 1994 and Fraser and Bernatchez 2001). This is particularly true for those aquatic species characterized by long ontogenic migrations (such as tuna fish, marlin, turtles, and whales) and shifts between marine and freshwater habitats such as for diadromous species (Maes and Volckaert 2007; McDowall et al. 2009).

The catadromous European eel Anguilla anguilla (L.) represents an emblematic case of a fish species with a complex life history, a high conservation interest, and an unresolved and much debated genetic structure (Maes and Volckaert 2007; Palm et al. 2009; Pujolar et al. 2009). Adult eels are believed to reproduce in a common spawning ground (SG) in the Sargasso Sea where larval stages have been found (Schmidt 1923; McCleave 1993). Eel larvae (leptocephali) are then transported through ocean currents across the North Atlantic Ocean to continental feeding grounds in fresh and coastal waters, ranging from northern Africa to the Scandinavian Peninsula. Eel larvae metamorphose to ‘glass eel’ before approaching the continental waters and soon after to ‘yellow eels’, the pre-reproductive stage. Yellow eels spend a highly variable time period in continental waters, ranging between 2 and 20 years on average, depending on (i) sex, with males being smaller and taking substantially less than females to reach sexual maturity and (ii) latitudinal position, with eels in Northern Europe taking substantially longer than those from Mediterranean regions (Vøllestad 1992). When reaching partial sexual maturity, eels undergo a final metamorphosis into ‘silver eels’ and migrate back to the Sargasso Sea where they reproduce and die (Tesch 2003). Considerable mixing of adult eels coming from different continental locations can take place during the North-Atlantic migration, whose duration is still debated (van Ginneken and Maes 2005; Aarestrup et al. 2009).

The urgency for the conservation of the European eel arose, because in the last decades, eel stock and glass eel recruitment started a dramatic and still ongoing decline (Dekker 2000). The ultimate culprit of this collapse is still unknown and has to be found in a combination of causes ranging from overexploitation to habitat loss, contamination by persistent organic pollutants and, possibly, global climate change (ICES 2008). Eel status is currently so threatened that A. anguilla was included in the IUCN Red List of critically endangered species (Freyhof and Kottelat 2008). The EU Regulation 1100 issued in 2007 required all Member States to implement eel management plans. However, it is difficult to assess whether these plans will be successful, as the large-scale geographic structure of the population and the processes driving the demographic dynamics of whole stock still remain poorly understood (Maes and Volckaert 2007). As a consequence, it remains uncertain whether eel should be managed as one single panmictic stock or several partially independent sub-populations.

Over the last decade, several studies investigated the genetic structure of the European eel to explore whether the species truly consists of a single panmictic population or several genetic units [see van Ginneken and Maes (2005) and Maes and Volckaert (2007) for comprehensive reviews]. The view of a fully panmictic population was initially challenged by three different studies (Daemen et al. 2001; Wirth and Bernatchez 2001; Maes and Volckaert 2002) that reported evidence for large-scale geographic differentiation using various classes of genetic markers. Both Wirth and Bernatchez (2001) and Maes and Volckaert (2002) suggested the existence of a stock subdivided into roughly three subpopulations corresponding to the Mediterranean, Atlantic, and North-Baltic regions, interconnected by a significant amount of gene flow. Data from all three studies showed signals ranging from clinal variations in allele frequencies or haplotype diversity to an isolation-by-distance (IBD) pattern (Daemen et al. 2001; Wirth and Bernatchez 2001; Maes and Volckaert 2002). These results were congruent with known oceanic current patterns distributing eel larvae in Europe and a clinal variation in vertebrae numbers (Boëtius 1980). It should be noted that the maintenance of such a stable large-scale geographic structure would require not only different active migration routes and spawning areas for the three subpopulations but also a restriction in larval mixing during the drift through the Gulf Stream back to the continental sites. However, the level of spatial segregation among subpopulations necessary to produce such large-scale genetic structure and the compatibility of larval/adult segregation levels with oceanographic drift conditions (Kettle and Haines 2006; Bonhommeau et al. 2009) have never been investigated in a quantitative framework.

Subsequent studies analyzing eel genetic structure using only eels belonging to the same cohort focused on genetic differentiation at smaller scales. These studies evidenced that temporal genetic differences between groups of glass eels recruited at the same site within the same year (also referred to as glass eel waves) can exceed between-site large-scale spatial genetic differences of the same recruitment cohort (Maes et al. 2006, 2009; Pujolar et al. 2006, 2007, 2009). In agreement with these observations, such a small-scale genetic patchiness of glass eel waves would result from independent mating events involving small groups of breeders separated by space or by time (Maes et al. 2006; Pujolar et al. 2009). According to this theory, the smaller the breeding group, the higher the observed between-wave genetic differences. A rigorous quantitative analysis of the number of eels per breeding group that are required to produce the observed within-site temporal genetic differences was never carried out.

Finally, a common problem in studies of broadly distributed marine organisms is the potential bias in the estimated genetic structure generated by limited sampling effort (Pujolar et al. 2006). Insufficient sampling can severely influence the outcome of genetic analyses either by swamping the existing temporal/geographic genetic structure or by over-estimating signals of geographic genetic differentiation, when the actual level is much lower (Waples 1998; Pujolar et al. 2007). In this specific case of A. anguilla, a quantitative analysis of how the actual number of glass eel waves sampled in a recruitment year might affect the estimation of large-scale genetic structure has not yet been carried out.

Only very recently, samples of larvae from the SGs have been analyzed with genetic markers, providing additional strong evidence for the panmictic stock hypothesis and reiterating the complex transient nature of patterns of isolation by distance and time (T.D. Als, M.M. Hansen, G.E. Maes, M. Castonguay, L. Riemann, K. Aarestrup, P. Munk, H. Sparholt, R. Hanel and L. Bernatchez, unpublished data). However, although these results provide crucial new insights into the genetic structure of eels, the lack of spawning adult samples, the restricted total sampling size and the small temporal sampling window of the study only provide a snapshot of the total breeding population structure in the Sargasso Sea. It is hence apparent that the complexity of the life cycle and, specifically, the difficulty of monitoring sufficiently broadly and for a longer period the oceanic phase of eel hindered so far a clear interpretation of genetic structure in terms of the ecological, demographic, and environmental processes shaping it. In this study, we introduce a new approach to tackle part of this problem by developing the first integrated demographic-genetic model for the European eel. The model was developed with the aim of investigating alternative hypotheses on the demographic and environmental processes driving the observed large-scale continental geographic structure and the small-scale temporal differentiation of the eel stock. We intended to model and partition genetic variation in a temporal and spatial component using the most realistic demographic data available for the eel to this date, while investigating the influence of different levels of gene flow during and after spawning. As such, we aimed at reconciling the most recent results of complete panmixia, geographic, and temporal differentiation during the oceanic and continental phase with a modeling approach, to open up a new research route for more advanced modeling exercises including newly appearing genetic or demographic results. Specifically, we aimed at answering the following questions:

what is the expected number and size of reproductive groups required to produce the temporal genetic differentiation typically observed between glass eels waves during a recruitment year (Maes et al. 2006);

what level of adult and larval segregation during ocean migration is required to result in a significant and plausible geographic genetic differentiation (Daemen et al. 2001; Wirth and Bernatchez 2001; Maes and Volckaert 2002) or in complete panmixia (Als et al., unpublished data);

what is the influence of sampling only a limited number of glass eel waves on the level of geographic and temporal genetic differentiation.

Materials and methods

Experimental studies on genetic structure carried on in the last 10 years have been used as a benchmark to compare the results produced by our mechanistic, i.e. process-based, model of eel demography, and population genetics under alternative and competing assumptions on large-scale geographic differentiation and small-scale temporal differentiation. The demographic component describing the eel life cycle during the pre-reproductive continental phase was built over the extensive modeling work carried out in the last 20 years by Vøllestad and Jonsson (1988) in Norway, De Leo and Gatto (1995) in Italy, Dekker (2000) in the Netherlands, Bevacqua et al. (2007) in France, etc. As for the migration of eel larvae from the Sargasso Sea to the continental waters, we relied on the results of Kettle and Haines (2006) and Bonhommeau et al. (2009) to estimate larval survival, while we used information on the estimated abundance of silver eels leaving the continental waters (Dekker 2000) and eel fertility (Boëtius and Boëtius 1980) to derive a realistic range of oceanic survival of spawners and their reproductive success. With the exception of an over-simplified version of eel whole life cycle by Åström and Dekker (2007), which did not account for the spatial structure of the stock in the continental phase, no models have been published yet in the literature describing neither the full life cycle nor the genetic structure of the European eel.

The model

Population dynamic features were simulated by a life-stage- and age-structured model from a continental perspective, as described hereafter and in Fig. 1. In order to explore the hypothesis proposed in earlier studies (Daemen et al. 2001; Wirth and Bernatchez 2001, Wirth and Bernatchez 2003 and Maes and Volckaert 2002) on the large-scale geographic structure of eel stock and to enable the joint assessment of spatial and temporal variation while keeping the population structure straightforward, we simplified the model for analytic reasons assuming that:

Figure 1

Life cycle of the European eel and model parameters: h (parameter controlling glass eel survival, Appendix A); δ, γ, ρ (parameters controlling sexual differentiation and maturation, Appendix A); M and F (natural and fishing mortality on the continent, Appendix A); L (body length, Appendix A); σL and σS (larval and silver eel survival rates, Appendix B); ASI and LSI (adult and larval segregation indices, eqns 1, 2); NB (total number of breeders); NG (number of reproductive groups); and NBGe (effective number of breeders per group, eqn 3).

the eel stock in the continental phase is divided into three major subpopulations (Fig. 2) corresponding to Northern (N, north of 50°N), Atlantic (A, between 35°N and 50°N), and Mediterranean (M, south of 35°N) regions;

Figure 2

Schematic representation of the population structure hypotheses tested in our model. Following larval dispersal, the groups of larvae are delivered to the subpopulation corresponding to their spawning area or to one of the other two subpopulations according to ψ (top scheme). Similarly, following adult migration, silver eels can end in the spawning ground corresponding to their subpopulation or in one of the other two spawning grounds according to ϕ (bottom scheme). N, Northern; A, Atlantic; M, Mediterranean.

the SG is also divided into three putative distinct areas (N-SG, A-SG, and M-SG) representing the primary spawning area of the corresponding hypothetical continental subpopulations.

the three subpopulations can be connected, either weakly (discrete populations) or strongly (panmixia), depending on how many individuals from one subpopulation end up in a different subpopulation during their journey either from the Sargasso Sea to continental waters in the larval phase or from the continental waters to the Sargasso Sea as reproductive adults, as described hereafter.

Let ψ be the probability of progeny produced in spawning area k (k = 1 for N-SG, 2 for A-SG, and 3 for M-SG) to be delivered to continental region j (j = 1 for N, 2 for A and 3 for M):

where j, k = 1, 2, 3 and 0 ≤ LSI ≤ 1 and LSI represents a larval segregation index (0 corresponds to complete random mixing and 1 to complete segregation).

Similarly, the fraction ϕ of silver eels migrating from the continental region j to the spawning area k is computed as

where j, k = 1, 2, 3 and 0 ≤ ASI ≤ 1 and ASI represents an adult segregation index (Fig. 2). As a consequence, ASI = LSI = 1 represents the extreme case of three completely distinct and independent populations; ASI = LSI = 0 the opposite case of perfectly overlapping SGs or events (i.e. of a single homogenous population); any combination of ASI and LSI between 0 and 1 implies some gene flow among subpopulations. ASI > 0 can be also interpreted as temporal separation of breeders: for instance, in the extreme case of ASI = 1, we can either assume that there exist three different SGs or, alternatively, that there exist three subpopulations breeding potentially in the same spawning area but at different times during the reproductive season. This may happen if the migration to the Sargasso Sea is scarcely or not completely synchronized among eels coming from distant geographic regions, concordant to departure time variation (Anthony Acou, pers. comm.). In summary, the model can mimic any combination of spatial and temporal segregation that would lead to an overall level of adult segregation measured by the parameter ASI.

In order to account for small-scale structure of the breeding stock, as suggested by Pujolar et al. (2009) and Maes et al. (2006), breeding in each spawning areas was assumed to occur in small spawning groups where reproduction effectively takes place. Reproduction was never observed in the wild; laboratory studies are contrasting, suggesting a batch spawning as well as a mass spawning behavior (Pedersen 2003; van Ginneken and Maes 2005). In the present model, we assumed spawning to take place in a single mass event among eels belonging to the same reproductive group (van Ginneken and Maes 2005).

For a given number of breeders NB, the size of each group (NBG, number of breeders per group) depends upon the number of groups NG in which breeders are randomly divided:

Offspring originating in each spawning group determines a glass eel recruitment wave, assigned to a specific continental region. The sex ratio of reproductive eels was taken into account by defining an effective number of breeders per group as:

where η is the fraction of females in each spawning event (Crow and Kimura 1970; see Appendix B for the estimation of η).

Each cohort in each subpopulation in the continental phase is characterized by its own genotype frequency distribution (GFD). For those silver eels surviving the migration back to the Sargasso Sea for reproduction, the genotype of each breeding individual is randomly drawn from the GFD of the cohort and subpopulation it belongs to. Number of eggs produced by each female in each reproductive group is estimated from its body size according to Boëtius and Boëtius (1980) and eggs are randomly fertilized by breeding males of the same reproductive group. Offspring genotype frequencies are calculated as the product of male and female gamete frequencies (Crow and Kimura 1970, p. 44). Mutation occurs before union of gametes according to a stepwise mutation model with mutation rate μ = 5 × 10−4 (Balloux and Lugon-Moulin 2002). Breeders die after reproduction (semelparity). Offspring produced in each reproductive group remain subdivided into larval waves characterized by their own GFD. When these larval waves approach the continental coasts, individuals that survived the migration from the Sargasso Sea metamorphose and represent new glass eel recruitment waves that are then recruited in one of three continental sub-populations according to eqn (1). This means that each wave comes from only a single spawning area. The GFD of glass eel waves are used to estimate genetic differentiation without taking into account the sampling error encountered in empirical studies (see ‘Indices of genetic differentiation’). The estimation of genetic differentiation indices occurs before merging the glass eel waves of the same cohort in each continental region so as to derive the GFD of the newly recruited cohort.

The process of drawing genotypes for each reproducing individual from the corresponding GFD of the cohort it belongs to as well as the assignment of larval groups to different continental regions according to eqn (1) add a stochastic component to the model, all the other parameters being constant. To account for these sources of variability, the indices of genetic differentiation were averaged on 100 simulation replicates for each parameter set (ASI, LSI, NG, etc.). We initially set genotype frequencies in Hardy–Weinberg proportions, based on allelic frequencies drawn from a uniform distribution. Simulations were run for 300 years, which was largely sufficient for genetic differentiation to converge to stable values (drift-migration equilibrium), except for extreme cases of segregation (i.e. ASI or LSI ≥ 0.99) for which convergence was slower, requiring running the model for 3000 years.

The model component describing eel demography in the continental phase explicitly accounts for key aspects of its life history, such as natural and fishing mortality (Dekker 2000), sex-specific body growth (as silver females are remarkably larger than silver males; Melià et al. 2006a), size- and sex-specific maturation rates (as females takes substantially longer than male to reach sexual maturity; Bevacqua et al. 2006), and size-dependent fertility of females (as larger females produce substantially more eggs than smaller females; Boëtius and Boëtius 1980). Moreover, the model accounts for geographic variations in life-history traits, as yellow eels in the Mediterranean area grow faster and take less to reach sexual maturity than eels from the central and northern Europe (Vøllestad 1992) but suffer higher mortality rates. Accordingly, a detailed demographic model was developed so as to mimic these geographic and between-sex differences in age and size at sexual maturity and to allow us to realistically simulate the gene flow occurring when spawners from different cohorts and different geographic areas mate in the Sargasso Sea. A detailed description of the full demographic model and its parameterization is reported in Appendix A. By setting larval survival equal to 1.5 × 10−3 as in Bonhommeau et al. (2009), survival of silver eels migrating to the Sargasso Sea and density-dependent mortality of glass eels have been derived under steady-state hypothesis (Appendix B), by assuming that silver eel abundance at equilibrium is about 2 × 107 and glass eel abundance 109 as estimated by Dekker (2000). All the other demographic parameters have been set according to published works.

Indices of genetic differentiation

Genetic differentiation was estimated using genotype and allele frequencies of glass eels measured through hierarchical F-statistics. The components of the overall differentiation between recruitment waves FST were partitioned into a within-subpopulation within-cohort FSC (temporal component) and a between-subpopulation within-cohort component FCT (geographic component). In empirical studies, estimates of genetic differentiation are derived from sample data; therefore, statistical techniques are used to correct for sampling bias [e.g. the θ parameter developed by Weir and Cockerham (1984) to estimate FST]. In our model, however, glass eel genotype frequencies are assumed to be known and can be used directly to define the true FST and related indices as ratio of heterozygosities. FSC is computed using all the recruitment waves at a continental region:

where HC is the expected within-subpopulation heterozygosity, computed as the weighted average of the expected heterozygosities in the three geographically distinct subpopulations; HS is the expected within-recruitment wave heterozygosity, computed as the weighted average of the expected heterozygosities in all the recruitment waves within a region. Between-subpopulation differentiation is computed as

where HT is the expected heterozygosity in the whole population. Finally, the overall index of genetic differentiation FST is:

These indices satisfy the well-known relationship (1 − FST) = (1 − FSC)(1 − FCT).

A first set of simulations was run by using 5–10–20 alleles, respectively, and a number of loci ranging from 1 to 10. As reported in Appendix D, a sensitivity analysis showed that the number of loci and alleles did not significantly affect model results for FST, FSC, and FCT. Hence, even though multiple loci would provide greater precision in a single model run, for reasons of computational tractability, genetic structure was modeled by using only one locus with 10 alleles.

Scenario analysis

The model was used to estimate values of FST, FCT and FSC under different assumptions about (i) the level of small-scale genetic differentiation, modulated through the number of reproductive groups (NG = 1500, 3000 and 6000); (ii) the number of breeders (NB = 0.42 × 106, 0.83 × 106 and 1.66 × 106); and (iii) the connectivity level among subpopulations (LSI and ASI = 0, 0.5, 0.99 and 1; the combination ASI = LSI = 1 was not considered because it corresponds to three independent populations). The F-statistics were computed by both using all the glass eel waves recruiting at each continental region and a finite number (i.e. 5 and 10) of glass eel waves to mimic the fact that only a subset of arrival waves is actually sampled in field studies.

Results were used to investigate the relationship between FSC and NBGe and whether a particular value of NBGe yields an FSC equal to that estimated from microsatellites in Maes et al. (2006), i.e. 0.0018 ± 0.0014 (hereafter this specific value is referred to as NBGe*). Similarly, simulation results were used to explore if the same level of genetic differentiation could be derived by assuming large-scale geographic segregation. Finally, the model was used to investigate how FSC, FCT, and FST estimates change when a small number of glass eel waves per continental region are randomly sampled. Median values of the indices of differentiation were compared by a Kruskal–Wallis nonparametric ANOVA (KW test).

Results

Small-scale genetic differentiation FSC

The index of within-subpopulation genetic differentiation FSC decreased linearly with the inverse effective size of breeding groups, NBGe:

where c = 0.5137 ± 0.00059 (mean ± standard error estimated on 100 replicates; linear regression over FSC and 1/NBGe, R2 = 0.99, F1,98 = 752.640, P < 0.0001; Fig. 3). If the level of within-subpopulation differentiation FSC derived from molecular analysis is equal to 0.0018 ± 0.0014 (Maes et al. 2006), then the effective number of breeding eels leading to the observed FSC, computed from eqn (8), is NBGe* = 295. By taking into account the uncertainty in FSC from molecular estimates, the lower and upper bounds for NBGe* are 160 (when FSC = 0.0032) and 1284 (when FSC = 0.0004).

Figure 3

The estimated relationship between NBGe and FSC (solid line) and model-derived FSC values under different demographic scenarios. Notice that markers overlap when their relative scenarios correspond to the same NBGe. The dotted line and the grey shaded area indicate the observed FSC = 0.0018 ± 0.0014 (Maes et al. 2006).

The value of the parameter c appearing in eqn (8) was derived assuming intermediate levels of adult and larval migration (ASI = LSI = 0.50); as shown in Table 1; however, the value of FSC was affected only minimally by assumptions on the connectivity level among the three subpopulations. The estimated values of c and NBGe* are therefore robust to large uncertainties in the levels of adult segregation and larval mixing. The number of breeding eels NBG was defined as the ratio between the total number of breeders NB and NG (eqn 3). The derived relationship between FSC and NBGe (eqn 8) held true for any combination of breeders abundance NB and number of groups NG (Fig. 3) and yielded the same effective number of breeders NBGe* corresponding to the observed FSC.

Table 1

Mean (and standard deviation) of genetic differentiation indexes under different scenarios of adult and larval segregation

		ASI
	LSI	0	0.5	0.99	1
FSC (× 10−3)	0	1.76 (0.01)	1.76 (0.02)	1.75 (0.01)	1.76 (0.02)
	0.50	1.76 (0.02)	1.75 (0.01)	1.76 (0.02)	1.76 (0.01)
	0.99	1.78 (0.01)	1.75 (0.01)	1.81 (0.02)	1.81 (0.02)
	1	1.82 (0.01)	1.79 (0.02)	1.81 (0.01)	–
FCT (× 10−6)	0	1.48 (0.44)	1.06 (0.39)	1.5 (0.57)	1.59 (0.36)
	0.50	1.69 (0.69)	1.47 (0.43)	1.67 (0.83)	1.65 (0.50)
	0.99	1.50 (0.32)	1.65 (0.48)	6.88 (2.30)	11.85 (5.95)
	1	1.61 (0.33)	1.75 (0.56)	25.47 (17.07)	–
FST (× 10−3)	0	1.76 (0.01)	1.76 (0.02)	1.75 (0.01)	1.77 (0.02)
	0.50	1.76 (0.02)	1.75 (0.01)	1.76 (0.02)	1.76 (0.01)
	0.99	1.78 (0.01)	1.75 (0.01)	1.81 (0.03)	1.82 (0.02)
	1	1.82 (0.01)	1.79 (0.02)	1.84 (0.03)	–

These figures have been derived by setting NBGe* = 295 and all the other parameters as described in ‘Methods’ and in Appendix A.

Large-scale and overall genetic differentiation, FCT, and FST

Spatial differentiation among the three continental subpopulations, FCT, increased with the levels of adult and larval segregation, ASI, and LSI (Table 1). However, FCT was always very small (in the order of 10−6–10−5) and never reached values that would be detectable using molecular markers (10−4–10−3).

The overall level of genetic differentiation, FST, was minimally affected by the levels of adult and larval segregation. For most of the range of variation of ASI and LSI, FST was not significantly different from within-subpopulation differentiation FSC. FST was affected by the level of population connectivity only when ASI and LSI were set to extreme values (ASI = 0.99 and LSI = 1), rising from 0.00176 to 0.00184.

Sampling bias

The number of glass eel waves actually sampled to derive the F-statistics significantly affected the estimates of FSC, FCT and FST (Fig. 4). When only 5 or 10 glass eel waves per continental region were analyzed (LSI = 0.5 and ASI = 0.5), both FSC and FST were smaller than the corresponding value computed using all the arrival waves (KW test: , P < 0.001 and , P < 0.05 respectively), while FCT was larger (KW test: , P < 0.001).

Figure 4

Influence of sampling bias on hierarchical F-statistics. Box and whisker plot of FSC, FCT, and FST estimated using different number of glass eel arrival waves and setting ASI = LSI = 0.50 and NBGe = 295. The box delimits the lower quartile and the upper quartile values; the horizontal lines within the box indicate the median. The whiskers extend upon a range equal to 1.5 times the inter-quartile range.

Discussion

Small-scale differentiation

Our analysis confirmed the hypothesis that genetic differentiation in glass eels can result from the subdivision of breeders in a high number of isolated spawning events and predicted that genetic differentiation is inversely proportional to the number of breeders in each reproductive group. The model predicted the observed value of FSC by assuming that each reproductive event involves as few as ca. 300 individuals. Assuming that the overall reproductive stock NB is of the order of 106, our model suggests that there might be as many as a couple of thousands reproductive isolated events; most likely, these groups are partitioned by a combination of spatial and temporal mechanisms (Wirth and Bernatchez 2001; Maes et al. 2006).

Our findings do not rule out that other mechanisms not explicitly simulated in our model might contribute to produce the observed genetic structure. For instance, the assumption of random mating within reproductive groups implies that all breeding silver eels contribute equally to offspring, apart from differences in female fecundity that are taken into account (eqn C3). However, differences in reproductive success may arise as a result of random processes, especially in the marine environment. Extreme heterogeneity in physical and dynamical conditions can lead to unequal parental contribution to the gamete pools. High variance in reproductive success among adults can considerably increase genetic drift and thus genetic differentiation between larval groups (Hedgecock 1994; Hedrick 2005). Environmental variation can also differentially affect larval groups, leading to family correlated survival; this further increases variability in parental contribution and leads to smaller effective number of breeders (Waples 2002), stronger genetic drift and higher potential for genetic differentiation between larvae (Maes et al. 2006; Pujolar et al. 2006).

The interpretation of our results based on genetic differentiation values as low as FSC = 0.0018 should be done with caution. Such low values, estimated from molecular markers, might be subject to considerable sources of inter-locus variance leading to wide confidence intervals (FSC = 0.0018 ± 0.0014; Waples 1998; Maes et al. 2006) or can be influenced by genotyping errors (e.g. allelic dropouts). Indeed, inter-locus variance introduced large uncertainty in the estimation of the number of breeding eels per reproductive group NBGe*, which ranged from 160 to 1294 effective breeders.

It is not possible to derive the actual number of male and female spawning eels from FSC. The effective number of breeding eels NBGe* derived through the model from molecular-based measures of genetic differentiation (FSC) corresponds to virtually infinite combinations of male and female silver eels in varying sex ratios, according to eqn (4). For instance, the same NBGe* = 295 could correspond to 111 females and 217 males, giving a sex-ratio of roughly two males for each female (like the one assumed in the model; Appendix C); or to 369 females and 92 males, giving a sex ratio of four females for each male; or to any combination satisfying eqn (4) and the constraints of stable population size given in Appendix B.

Direct observations of the actual sex-ratio of breeding eels in the field are not available and breeder sex ratio likely depend upon several uncertain factors and processes that are scarcely known or difficult to track, such as environmental sex determination occurring in continental waters (De Leo and Gatto 1996; Laffaille et al. 2006); sex-dependent growth and maturation rates (Bevacqua et al. 2006; Melià et al. 2006a,b); sex-related differences in fishing mortality, as the majority of fishing gears for eels are usually size selective (Bevacqua et al. 2009) and the larger females remain in coastal or inland waters for much longer periods than males; sex-related differences in natural mortality during the oceanic migration, as mortality is related to differences in body size and in relevant energy stores (van Ginneken and van den Thillart 2000), and the effect of endocrine disruptors in both sexual differentiation and in the accumulation of lipid energy stores which is notoriously different between female and male eels (Palstra et al. 2006; Belpaire et al. 2009). As a consequence, we did not have other option that restricts the results of the present model to be expressed in units of ‘effective’ breeders rather than actual individuals.

Large-scale geographic differentiation

Our analysis showed that even in the case of the largest feasible values of adult and larval spatial segregation (i.e. LSI = 0.99 and ASI = 1 or LSI = 1 and ASI = 0.99), the index of between-subpopulation genetic differentiation FCT remained below 10−5 (Table 1). Such low values for FCT estimated through molecular studies would hardly be significant neither in a statistical nor in a biologic sense (Waples 1998) and could be because of technical issues during genotyping. Hence, it is evident that, under the model assumptions, the three subpopulations would hardly be genetically different whatever the level of spatio-temporal segregation of the breeding stock. In particular, even strong levels of adult segregation (ASI) are compatible with the observed lack of genetic differentiation in glass eels, so that we cannot rule out the hypothesis that reproduction may take place in isolated areas of the SG or in distinct time periods. This hypothesis could be hopefully corroborated (or rejected) by the results of the genetic analysis of eel larvae sampled close to the SGs to test for spatio-temporal structure (Als et al., unpublished data). However, as shown in the previous section, the lack of large-scale geographic differentiation in this study does not preclude the existence of a fine-grade genetic differentiation caused by a small-scale temporal patchiness or cryptic spatial structure in the spawning stock (Dannewitz et al. 2005).

The absence of large-scale geographic genetic differentiation is compatible with very high levels of adult and larval segregation. In other words, genetic differentiation would be practically 0 even if the three subpopulations were connected by very small migration rates. Even if the level of geographic genetic differentiation observed in glass eels is not significantly different from 0 (Dannewitz et al. 2005), this means that we cannot rule out the possibility that the three subpopulations exist as independent demographical units. This may happen because the migration rate required to ensure demographic connectivity among subpopulations and synchronization of population dynamics is generally much higher than the rate sufficient to maintain genetic homogeneity (Waples and Gaggiotti 2006; Koizumi et al. 2008). The practical consequences of the existence of three distinct demographic units for the conservation of the stock would be crucial, as each subpopulation would require a distinct management plan. However, the generalized drop in recruitment observed in the overall stock (ICES 2008) would suggest that European eel truly consists of one single demographic unit. The actual attempt of an international coordination of local and national management plans under the EU directive 1100/2007 seems thus justified.

Finally, the absence of a spatial pattern of neutral genetic differentiation does not imply the absence of crucial genetic variation at the more relevant adaptive functional genetic level (linked to life-history traits). Although a recent eel population structure assessment using randomly spread EST-linked (Type I) markers did not evidence any spatial genetic differentiation (Pujolar et al. 2009), the detection of genetic variation underlying important phenotypic traits requires additional in depth analyses of the eel genome (ongoing study) (Maes and Volckaert 2007; Nielsen et al. 2009).

The role of sampling effort

Our results evidence that the observed signal of genetic structure can be overestimated at the geographic (large-scale) level, while underestimated at the temporal (small-scale) level when the number of arrival waves sampled to compute the F-statistics is too small. As shown in Fig. 4, even if as many as 10 arrival waves are sampled at each continental location, FSC can be severely underestimated and, consequently, FCT overestimated. This issue remains a problem in empirical studies, as 10 arrival waves per continental region substantially exceeds current sampling effort. As a practical recommendation, we advise population geneticists to stratify samples by cohort (as previously suggested by Dannewitz et al. 2005) and to collect genetic material over the whole recruitment season from different glass eel waves (as pointed out by Pujolar et al. 2006, 2007) for several locations (Maes et al. 2006). It would be interesting to see if empirical data sets support our findings that FCT decreases with the number of waves analyzed. This could be done, for instance, by using different number of glass eel wave samples from the same data set to estimate the F indices.

The importance of using many glass eel waves is to provide a correct representation of the genetic structure of the entire European eel stock partially derived from the model assumption that larval groups do not mix during the oceanic journey (cohesive dispersal). Nonetheless, a significant level of larval mixing is likely to occur because of oceanic conditions (Siegel et al. 2003) or behavioral mechanisms (Fraser et al. 2005; Bonhommeau et al. 2009). Additionally, glass eel waves can mix in front of continental coasts while waiting for favorable conditions to migrate into inland waters (Boëtius and Boëtius 1989; Tsukamoto et al. 1998). Finally, the potential for adult batch spawning suggested by some laboratory experiments (Pedersen 2003) may result in additional gene flow between larval groups, deriving from adults reproducing more than once in distinct groups, with consequences similar to mixing of larvae coming from different groups. Hence, the bias introduced by a limited number of samples on FCT might be less severe than that predicted in our model if the hypothesis of cohesive larval dispersal were relaxed.

Reconciling earlier and current knowledge on eel genetics

There are several components to take into account when analyzing marine organism populations to avoid drawing conclusions based on incomplete or biased sampling of populations or the genome (Waples 1998; Waples et al. 2008). The main issues to consider are lack of samples from the spawning stock (adults); samples comprising different cohorts and life stages; lack of temporal stability of the genetic differentiation patterns; temporal restricted sampling of batch spawning groups during the protracted spawning season, a low signal-to-noise ratio because of inter-locus variance.

For instance, most studies showing a clinal pattern in eel genetic structure used samples from different cohorts or life stages (Daemen et al. 2001; Wirth and Bernatchez 2001; Maes and Volckaert 2002), potentially leading to a geographic genetic differentiation if life-stage difference is related to geographic location and life stages differ genetically (Dannewitz et al. 2005). On the other hand, not all studies exhibited this age–geography relation and the IBD pattern remained visible when analyses were repeated only within a single glass eel cohort (Thierry Wirth, pers. comm.). Moreover, if temporal inter-cohort variation could entirely explain the observed geographic variation, one would not expect to observe a clinal variation, but rather an unordered differentiation pattern.

Alternatively, the reasons for clinal genetic variation could lie in a real biologic pattern emanating from various processes: (i) a fluctuating rate of introgression between the European and American eel (Anguilla rostrata), affecting the apparent but weak genetic differences between A. anguilla samples from North to South, because of nonrandom dispersal. The very fact that introgression seems concentrated in the Northern distribution area, as seen in various meristic and genetic studies (Tesch 2003; Albert et al. 2006; Gagnaire et al. 2009), points to a nonrandom dispersal of hybrid offspring ultimately leading to such a clinal pattern (only visible if the studied markers exhibit sufficient genetic differentiation between species or are linked to such markers). However, because of the fluctuation in introgression rates every generation and during every reproductive event, such patterns are not stable in time, (ii) Although larval mixing during dispersal can be extensive, the most recent results analyzing larvae collected in the Sargasso Sea showed the probability to detect an IBD pattern if siblings are included in a recruitment wave (Als et al., unpublished data). This could also have been the case in several studies considering strong larval retention until the continent, and (iii) Finally, the combination of drift effects and differential oceanic dispersal history could lead to a transient pattern of isolation by distance or time only in certain years.

Hence, although many earlier eel studies clearly suffered from one or several from these shortcomings, the emergence of a subtle but transient genetic structure in European eel points to the possibility of real biologic differentiation only picked up during certain years and under specific conditions. The further development of more complex geographic models or Lagrangian oceanographic models of larval dispersal might shed light on this aspect and will be object of future work. Additionally, the genetic model presented here will form the basis for further modeling including various new components, such as additional markers (type I genetic markers, SSRs or SNPs; Maes and Volckaert 2007), the evaluation of a continuous geographic distribution area, the inclusion of introgression rates, and the inclusion of markers subjected to selection as better tracers of geographic selection mechanisms in the genome.

Model limitations and improvements

Despite our efforts to provide a realistic representation of the complete eel life cycle and its complex demographic and genetic components, our modeling approach is obviously not exempt from limitations. In addition to the issues already discussed (mass spawning; variation in reproductive success; cohesive larval dispersal; representation of the geographic distribution of eel in terms of three continental sub-populations; the impossibility to derive the actual number of male and female parents), there are three further assumptions on which our simulation exercise was built that deserve attention, that is: all larval groups have the same survival rate; survival and fecundity rates are constant in time; demographic parameters do not vary within the three large subpopulations. The implications of these simplifications and possible future developments are discussed in the following.

Same survival rate for all larval groups

The model assumes that all larval groups experience the same survival rate. In reality, survival rate can be extremely variable among groups because of extensive variability in marine environmental conditions (Hedgecock 1994). As a consequence, some larval groups may undergo mass mortality and never reach the continental shelves. Such variability in larval production can lead to prominent reduction of the effective number of breeders (Hedrick 2005) and a consequent reinforcement of genetic differentiation. The implications for interpreting the model result are not easy to predict. One possible implication is that the total number of breeders NB corresponding to the observed FSC would be higher; but this would be of little practical importance for the results, as the focus is on the effective NBGe rather than the total and NB has a large confidence interval (Appendix B). This particular aspect may be improved by coupling the present model with a Lagrangian circulation model (e.g. Kettle and Haines 2006) capable of describing larval dispersal and the probability of successful recruitment to continental waters.

Temporal invariability of vital rates

Following Bonhommeau et al. (2009), we assumed that the eel stock was at a steady state, while in reality eel recruitment was characterized by a 90% drop in the last 30 years. This assumption was required to derive otherwise unknown demographic parameters (i.e. fraction of silver eels successfully migrating to the Sargasso Sea for reproduction). F-statistics were thus computed assuming a constant recruitment whose abundance corresponds to the more recent estimates (Dekker 2000). This implies that our model expresses the expected relationship between demography and genetic structure at a steady state and does not represent recent population dynamics. Our modeling approach could be used in the future to explore scenarios with variable population size, in particular to analyze the effects of the recent population collapse (ICES 2008). This would allow us to describe temporally unstable patterns of genetic differentiation that may be because of the inter-annual variability of oceanic conditions (see previous section).

Within-subpopulation variation in demographic parameters

Finally, we have not explicitly considered the effect of small-scale habitat variations on growth and survival during the continental phase of eel life cycle. In fact, in salt marsh and coastal lagoons, eels grow faster but suffer higher mortality and have smaller size at sexual maturity than in rivers and freshwater lakes. However, the most important coastal lagoons (in terms of eel production) are located in the Mediterranean Sea and our model parameterizations for growth and survival already reflects this large-scale pattern.

While aware of these and other limitations, we stress that our model structure mostly reflected the state of the art of knowledge of life-history traits and ecological processes available from literature. More detailed mechanistic models of specific demographic, biologic, and environmental processes could have been certainly developed without the possibly of realistically calibrating them on data. The simplifications we listed earlier were thus deliberately applied to calibrate the model according to the available literature.

Conclusions

The present work represents the first attempt to model the full life cycle of the European eel and its genetic structure, and we are confident that our findings have brought new insights into the complex life-history traits governing population genetic and demographic processes in this diadromous species. The main findings of our work are threefold: first, we showed that even very low larval mixing during dispersal or adult mixing occurring during migration to the Sargasso Sea are sufficient to overshoot geographic genetic differences among sub-populations. Second, we showed that the observed within-cohort genetic differentiation can be generated by a spatio-temporal subdivision of breeders in distinct reproductive groups constituted by a limited effective number of spawning eels (of the order of few hundred individuals). Therefore, complete panmixia seems unlikely and fine scale subdivision of the breeding stock likely occurs in the spawning areas. In addition, a potential latitudinal geographic variation would be only visible during periods of introgression of American eel material, followed by nonrandom mixing. Third, we further showed that limited sampling effort of few arrival waves may lead to an overestimation of large-scale genetic differentiation.

The practical implications of these findings are manifold. First, given the limited evidence of a strong large-scale differentiation, it seems more than reasonable to protect eels as a single tock unit, as implemented in the European policy, and to stock eel recruits from areas of greater abundance to those of lower abundance. Second, great attention should be given to guarantee a stratified sample of a sufficient numbers of individuals and/or glass eel waves in order to provide a robust representation of eel genetic structure. Third, on the research side, more effort should be put on casting light on the numerous aspects of eel life cycle that are still unknown or highly debated, such as the mechanisms driving (e.g. density) or affecting (e.g. endocrine disruptors) sexual differentiation and survival both in the continental and in the oceanic phase (De Leo and Gatto 1996; Laffaille et al. 2006; Palstra et al. 2006; Belpaire et al. 2009) and migration routes and time to and from the spawning areas (Aarestrup et al. 2009; Bonhommeau et al. 2010).

Our approach can be applied in the future to test various scenarios of the genetic structure of other eel species from different parts of the globe, as long as basic demographic and population genetic data is available. This will eventually improve our general understanding of eel biology and ecology. In particular, the model can be extended to take into account the association between adult spawning and thermal fronts (van Ginneken and Maes 2005) and the consequences of temporal instability of thermal fronts on genetic structure. Our mechanistic, i.e. process-based, approach can inspire further work to develop a new generation of integrated seascape models such as in Galindo et al. (2006, 2010) based on a careful combination of oceanographic, genetic and demographic data for species with complex life cycle and/or long migration ranges. Such approach can ultimately be implemented into new management strategies (De Leo and Gatto 1995), as to include potential evolutionary effects on adaptive variation of geographically focused (Northern vs southern Europe), sex-biased (mainly females) or life-stage-specific overfishing, leading to irreversible fishery-induced evolutionary response of the harvested stocks (Heino et al. 2002; Jørgensen et al. 2007). Integrated modeling approaches as presented here will foster initiatives for evolutionary enlightened fishery management and conservation of species with complex life cycle and long migration ranges.

Table A1

Model parameters

N*	663 mm	Length at maturity for females†
A*	664
M*	572
N	406 mm	Length at maturity for males†
A	385
M	405
αN	2.77	Correction factor for age at maturity
αA	1.59
αM	1.26
βYM,N	1.18	Correction factor for length at maturity (males)
βYM,A	1.12
βYM,M	1.18
βYF,N	1.23	Correction factor for length at maturity (females)
βYF,A	1.19
βYF,M	1.06
γMAX,F	1	Maximum rate of metamorphosis for females‡
γMAX,M	1	Maximum rate of metamorphosis for males‡
ε	0.777 × 106	Number of eggs per reproductive female
η	0.34	Fraction of females among breeders
ηF	26.2 mm	Slope parameter of the metamorphosis function for females‡
ηM	15.4 mm	Slope parameter of the metamorphosis function for males‡
λF	541 mm	Semi-saturation constant of the maturation function for females‡
λM	344 mm	Semi-saturation constant of the maturation function for males‡
μ	5 × 10−4	Mutation rate§
ξdiff,N	5 years	Duration of the undifferentiated stage
ξdiff,A	3
ξdiff,M	3
ξYM,N	15 years	Duration of the yellow eel stage for females†
ξYM,A	10
ξYM,M	6
ξYF,N	10 years	Duration of the yellow eel stage for males†
ξYF,A	5
ξYF,M	5
ρ	0.5	Fraction of undifferentiated eels becoming males
σL	0.0015	Larval survival fraction¶
σYU,N	0.84	Undifferentiated survival fraction
σYU,A	0.69
σYU,M	0.65
σYM,N	0.69	Yellow male survival fraction
σYM,A	0.44
σYM,M	0.44
σYF,N	0.79	Yellow female survival fraction
σYF,A	0.71
σYF,M	0.53
σS	0.028	Silver eel survival fraction
F	0.63	Cumulative fishing mortality rate**
h	1.38 × 10−8	Strength of density-dependent glass eel survival
kYF	1.66 × 10−3	Brody coefficient for yellow females††
kYM	3.05 × 10−3	Brody coefficient for yellow males††
kYU	1.00 × 10−3	Brody coefficient for undifferentiated††
Ldiff	221 mm	Length at sexual differentiation††
L0	60 mm	Length at age recruitment††
L∞,YF	573 mm	Asymptotic length for yellow females††
L∞,YM	386 mm	Asymptotic length for yellow males††
L∞,YU	399 mm	Asymptotic length for undifferentiated††
M	2.52	Cumulative natural mortality rate**

The three values refer to the Northern, Atlantic, and Mediterranean sub-populations, respectively.

Tesch (1977), Vøllestad (1992) and Durif (2003);

Bevacqua et al. (2006);

Balloux and Lugon-Moulin 2002;

Bonhommeau et al. (2009);

Dekker (2000);

Melià et al. (2006a).