Evaluating the adaptive potential of the European eel: is the immunogenetic status recovering?

The recent increased integration of evolutionary theory into conservation programs has greatly improved our ability to protect endangered species. A common application of such theory links population dynamics and indices of genetic diversity, usually estimated from neutrally evolving markers. However, some studies have suggested that highly polymorphic adaptive genes, such as the immune genes of the Major Histocompatibility Complex (MHC), might be more sensitive to fluctuations in population dynamics. As such, the combination of neutrally- and adaptively-evolving genes may be informative in populations where reductions in abundance have been documented. The European eel (Anguilla anguilla) underwent a drastic and well-reported decline in abundance in the late 20th century and still displays low recruitment. Here we compared genetic diversity indices estimated from neutral (mitochondrial DNA and microsatellites) and adaptive markers (MHC) between two distinct generations of European eels. Our results revealed a clear discrepancy between signatures obtained for each class of markers. Although mtDNA and microsatellites showed no changes in diversity between the older and the younger generations, MHC diversity revealed a contemporary drop followed by a recent increase. Our results suggest ongoing gain of MHC genetic diversity resulting from the interplay between drift and selection and ultimately increasing the adaptive potential of the species.


INTRODUCTION
Preserving natural biodiversity while allowing species to maintain their adaptive potential is a major challenge in modern conservation biology (Frankham, Briscoe & Ballou, 2002). Anthropogenic activities impact global ecosystems and reduce population sizes of species, whether by shrinking or fragmenting available habitats, overexploitation, or disruption of population dynamics (Allendorf et al., 2008;England, Luikart & Waples, 2010;Thomas et al., 2004). Small populations are more vulnerable to environmental, demographic and genetic processes (Keith et al., 2008). Genetic factors are particularly important as they may not manifest immediately after population reduction but their effects persist in the population even if the population size recovers to otherwise sustainable levels, e.g., through bottleneck effects (Spielman, Brook & Frankham, 2004). Populations display complex evolutionary dynamics and evolutionary genetics offers ideal framework for conservation biologists to monitor population changes and viability (Hendry et al., 2011). Although genetic studies applied to wild populations of non-model species have largely focused on the analysis of neutrally evolving loci (see McMahon, Teeling & Höglund, 2014), conservation managers have expanded their toolbox to include screening of adaptive loci (Hendry et al., 2011). Such measures are required to fill knowledge gaps regarding species' evolutionary and adaptive potential (Eizaguirre & Baltazar-Soares, 2014).
The genes of the Major Histocompatibility Complex (MHC) have repeatedly been shown to be suitable candidates to evaluate the immune adaptive potential of endangered populations (Sommer, 2005;Stiebens et al., 2013b). This highly polymorphic (high heterozygosity and gene duplications), multigene family (Apanius et al., 1997;Klein, Sato & Nikolaidis, 2007) plays a decisive role in controlling the vertebrate adaptive immune system by presenting self-and pathogen-derived peptides to T-cells (Janeway et al., 2001). Pathogen-mediated selection is acknowledged to be one of the primary factors of balancing selection maintaining the extreme MHC polymorphism in a population (Eizaguirre et al., 2012a;Piertney & Oliver, 2006;Spurgin & Richardson, 2012). Therefore, investigating shifts in MHC allele frequencies (Eizaguirre et al., 2012b) may be a particularly informative tool as an indirect way to detect the emergence of diseases (Sommer, 2005). Similarly, examining how MHC genetic diversity fluctuates in parallel with the incidence of diseases or parasites can provide indirect evidence for the impact of those selective agents on the dynamics of the host population (McCallum, 2008). Lastly, MHC genes have also shown to be informative of demographic events, particularly when selection plays an important role in population reductions (Sutton et al., 2011).
The European eel is a highly migratory, semelparous fish whose spawning grounds are located in the Sargasso Sea and whose foraging grounds cover coastal, mixohaline and freshwater habitats (Harrod et al., 2005) across much of Europe and even extend to North Africa and the Levantine coast. The post-hatching early larval transport is first facilitated by local currents in the Sargasso sea, connecting the spawning area with the Gulf Stream (Baltazar-Soares et al., 2014), and by the North Atlantic gyre, that completes the larval migration to the continental shelf (Bonhommeau, Chassot & Rivot, 2008;Kettle, Bakker & Haines, 2008b;Munk et al., 2010). Once larval eels enter shelf waters, they undergo a series of substantial changes in morphology and physiology: individuals become glass eels, with a fusiform, transparent body that facilitates active swimming towards coastal waters (Miller, 2009). The duration of the continental life stage varies on the location of growth habitats, and may last from as little as two years (EIFAAC/ICES, 2013) to several decades prior to metamorphose into silver eel. At this stage, energetic reserves are collected to allow sexual maturation and the long spawning migration back to the Sargasso Sea (Tesch, 2003).
Although the number of glass eels arriving at continental coasts across Europe experienced a first drop in the 1960s, the major recruitment collapse occurred at the beginning of the 1980s (EIFAAC/ICES, 2011). For the subsequent three decades, the recruitment of glass eels has remained as low as 1 to 10% of the values prior to the 1980s (EIFAAC/ICES, 2011). The low recruitment regime is hypothesized to have resulted from multiple impairing factors including productivity changes in the Sargasso Sea (Friedland, Miller & Knights, 2007), habitat degradation including river regulation, pollution and reduced freshwater habitats (Prigge et al., 2013;Robinet & Feunteun, 2002), changes in oceanic currents (Baltazar-Soares et al., 2014;Kettle, Bakker & Haines, 2008), introduction and spread of diseases, such as the EVEX (Van Ginneken et al., 2005) and the swim bladder parasite Anguillicola crassus (Kirk, 2003), resulting in a severe lack of spawners (Dekker, 2003). Experimental studies focusing on this invasive nematode for instance revealed that the European eel is unable to mount an effective immune response, thus becoming particularly susceptible to infection (Knopf, 2006).
The European eel is currently considered a single panmictic population (Als et al., 2011;Pujolar et al., 2014), even though punctual deviations from panmixia have been reported (Baltazar-Soares et al., 2014;Dannewitz et al., 2005;Wirth & Bernatchez, 2001). While the vast majority of genetic studies have focused on solving the population structure of the species, few have tried to determine the impact of the recruitment decline on the species' genetic diversity (Pujolar et al., 2011;Pujolar et al., 2013;Wirth & Bernatchez, 2003). Studies that aimed at doing so have evaluated neutrally evolving genetic markers. In 2003, Wirth andBernatchez (Wirth &Bernatchez, 2003) analyzed seven microsatellite loci in 611 European eel individuals and reported no measurable signature of the 1980s recruitment decline. In 2011, Pujolar and coworkers (Pujolar et al., 2011) analyzed the diversity of 22 microsatellite markers on 346 individuals. Again, genetic signature of a population reduction was absent. It was therefore suggested that the 1980s decline, although marked in terms of population biology was not sufficiently extreme to affect the diversity of these polymorphic and neutrally evolving loci (Pujolar et al., 2011). Employing a genome-wide reduced representation sequencing technique (RAD)-that identified in its vast majority neutral markers- Pujolar et al. (2013) did not detect either evidence for a recent decline of genetic diversity to be associated with the drastic drop in recruitment.
Here, we expand on these studies detailing an extensive evaluation of the current genetic status of the European eel species, which includes (1) screening of both neutral and adaptive markers and (2) a temporal approach directly comparing two distinct generations of eels. A temporal approach is regarded as a key requirement when investigating the signature of demographic events on genetic diversity in a wild population (Sutton, Robertson & Jamieson, 2015). Specifically for the eel system, evaluating genetic diversity in two distinct, non-overlapping age cohorts is important since the major component of genetic distribution of nuclear markers in this species seems to relate to different age cohorts (Dannewitz et al., 2005).

Study scheme
A total of 683 eels were analyzed in this study, 202 of which corresponded to mature silver eels caught in freshwater while undertaking their spawning migration. The other 481 individuals were glass eels collected from 2009 to 2012 immediately upon their arrival in coastal waters (Fig. S1). Number of individuals per sampling site, year of capture, developmental stage, and geographical locations can be found in Table 1. Note that we mostly sampled in one catchment per country and used countries' acronyms to refer to the locations where the specimens were collected. All samples were included in the analyses of mtDNA and microsatellites. Although the main objective of this work was to evaluate the genetic status of the eel species on a temporal scale, we also analyzed spatial patterns of genetic differentiation for comparison with previously published studies (Als et al., 2011;Maes & Volckaert, 2002;Pujolar et al., 2014;Wirth & Bernatchez, 2001).
Amongst the 683 individuals, 327 were sequenced at the exon 2 of the MHC class II B gene. This number was achieved as a compromise between costs of 454 sequencing and the need to provide sample sizes sufficient to provide a robust screen for spatial and temporal patterns of local adaptation. Locations and respective sample sizes of fish screened for the MHC gene are highlighted in Fig. S2. Temporal analyses were performed after dividing the dataset in two distinct age groups: the older ''silver eels'' age group and the younger ''glass eels'' group. The first group included individuals born soon after the drop in recruitment (late 1990s and early 2000s) while the second group consisted of very recently recruited individuals (2009 onwards) (ICES, 2015). DNA was extracted from fin clips (''silver eels'') or tail clips (''glass eels'') with Qiagen DNeasy Kit c Blood and Tissue kit (Hilden, Germany) following the manufacturer's protocol.

Genetic estimates of diversity, population differentiation and demography
All 683 glass eels were sequenced using Sanger sequencing for the mitochondrial NADH dehydrogenase 5 (ND5) exactly replicating (Baltazar-Soares et al., 2014). Haplotype diversity (Hd) and nucleotide diversity (π) were calculated for each sampling location in DnaSP v5 (Librado & Rozas, 2009). Genetic structure was estimated using Arlequin v3.5 with 10.000 permutations (Excoffier & Lischer, 2009). Moment-based demographic parameters that test for changes in effective population size were calculated for each sampling location in DnaSP v5 under the assumption of mutation-drift equilibrium. Tajima's D (Tajima, 1989) and raggedness ' r (Rogers & Harpending, 1992) were also calculated in DnaSP v5. Ninety-five percent confidence intervals were estimated through coalescence simulations using 1.000 permutations. We evaluated the nucleotide mismatch pairwise distributions within each geographical location (Rogers & Harpending, 1992). These distributions were compared to expected distributions under a constant population size and sudden population expansion (Librado & Rozas, 2009).
Genetic signatures of a bottleneck were tested for each location using the tests available in BOTTLENECK (Cornuet & Luikart, 1996). These methods are sensitive to recent and severe reductions on effective population size (N e ) (Cornuet & Luikart, 1996). A two-phase mutation model was assumed with 10% of the loci allowed to evolve through stepwise mutation (Kimura & Ohta, 1978). Allele frequency distributions were also calculated for each location.

Genetic estimates of diversity, differentiation and demographyinter-generation level
In order to compare genetic diversity and demographic histories between the two cohorts and to avoid sampling bias from disproportionate number of samples in the two eel age groups (''glass eels'' n = 481 and ''silver eels'' n = 202), we performed 10 rounds of resampling of the data without replacement using PopTools (Hood, 2010), hereafter referred to as ''replicates.'' Replicates were performed based on 50 individuals. This standardization is critical to validate future comparisons, as it has been long acknowledged that sample size affects the detection of genetic signatures of recent bottlenecks (Luikart et al., 1998) and the estimation of effective population size (Waples & Do, 2010).
Deviations from Hardy-Weinberg equilibrium (HWE) were calculated for each replicate in Arlequin v3.5 (10,000 permutations). Nei's unbiased heterozygosity (He), observed heterozygosity (Ho), allelic richness (Ar) and F IS were calculated and compared between groups of replicates, i.e., ''glass eels'' and ''silver eels,'' with two-sided t -tests in FSTAT (1,000 permutations) (Goudet, 1995). The distribution of genetic variance between ''glass eels'' and ''silver eels'' was assessed with an analysis of molecular variance (AMOVA, Arlequin v3.5) amongst groups of replicates. Demographic history was inferred using two approaches. First, we evaluated the possible genetic signature of the recent population decline using BOTTLENECK (1,000 iterations) for each age group as previously described. Second, we estimated the effective population size (N e ) of each replicate of ''silver eels'' and ''glass eels'' using the linkage-disequilibrium method implemented in NeEstimator V2.01 (Do et al., 2014). We utilized P crit = 0.05, since lower P crit can overestimate Ne (Waples & Do, 2008). All estimates were obtained with the composite Burrows method (Weir, 1990). The unweighted harmonic mean was calculated for each group according to the following where j is the number of replicates, i is a given replicate and

Adaptive marker: diversity and demography of the MHC
We amplified the exon 2 of the MHC class II gene that encodes for the peptide-binding groove of the molecule following protocols optimized for the European eel (Bracamonte, Baltazar-Soares & Eizaguirre, 2015). We used the forward primer AaMHCIIBE2F3 (5 -AGTGYCGTTTCAGYTCCAGMGAYCTG-3 ) and reverse primer AaMHCIIBE2R2 (5 -CTCACYTGRMTWATCCAGTATGG-3 ) which allow the amplification of different allelic lineages of the MHC class II β genes (Bracamonte, Baltazar-Soares & Eizaguirre, 2015). Sequencing was performed on a 454 c platform at LGC genomics (Belgium) following (Stiebens et al., 2013a;Stiebens et al., 2013b). Briefly, two independent reactions were prepared for each individual. After a first PCR of 20 cycles, a reconditioning step (dilution 1:5) was performed, and the template was used for a second PCR of 20 cycles. The reconditioning step combined with independent reactions was shown to significantly decrease the number of PCR artifacts (Lenz & Becker, 2008) and facilitate allele call (Stiebens et al., 2013a). The second set of PCR was performed using the specific MHC primers extended by the 454 adaptors and a 10 bp individual tag. Allele calling and respective assignment to individuals followed (Stiebens et al., 2013a;Stiebens et al., 2013b) and primarily relied on matching alleles present in both independent reactions (Sommer, Courtiol & Mazzoni, 2013). Genotyping using this method has previously been compared to Sanger sequencing and showed its high accuracy (Bracamonte, Baltazar-Soares & Eizaguirre, 2015). Even though variants may stem from different loci, we will refer to them as alleles hereafter.
Individual MHC allele numbers detected from the different paralogs, nucleotide diversity, and individual average nucleotide p-distance (Eizaguirre et al., 2012a) were calculated for each sampling location and for each group, i.e., ''silver eels'' and ''glass eels,'' in DnaSP v5 and using custom Perl scripts. MHC allele pools were compared amongst sampling locations and between ''silver eels'' and ''glass eels'' with analyses of similarity (ANOSIM) using Primer v6 (Clarke, 1993) following (Eizaguirre et al., 2011) 1,000 permutations. Correlation between MHC divergence and neutral structure was calculated using a Mantel test between pairwise F ST matrices (mtDNA and microsatellites) and pairwise Bray-Curtis similarity matrices (MHC).
Minimum number of recombination events (Rm) and estimates of recombination rate (R) were calculated in DnaSP v5 (Hudson & Kaplan, 1985), as well as the relative (R/θ) contribution of recombination (R) and point mutations (θ ) in the generation of genetic diversity (Reusch & Langefors, 2005). Gene conversion was investigated using ψ that measures the probability of a site to be informative for a conversion event (ψ > 0, (Betran et al., 1997)) between ''glass eels'' and ''silver eels,'' using a sliding window method (window length = 2, step size = 1) implemented in DnaSP v5.
In order to test for the mode of evolution of the MHC in A. anguilla, overall positive selection was estimated with a Z -test implemented in MEGA v5 (Tamura et al., 2011). We tested for signs of codon-specific positive selection using maximum likelihood site models with CODEML implemented in PAML v4.4 (Yang, 2007) and the mixed effects model of evolution (MEME) (Murrell et al., 2012) implemented in the Datamonkey web server (Delport et al., 2010;Pond & Frost, 2005). The maximum likelihood procedures evaluate heterogeneous ratios (ω) among sites by applying different models of codon evolution. Three likelihood-ratio tests of positive selection were performed comparing the models M1a (nearly neutral) vs M2a (positive selection), M7 (ß) vs M8 (ß + ω), and M8a (ß + ω = 1) vs M8 (Yang, 2007). In the models M2a and M8, positively selected sites are inferred from posterior probabilities calculated by the Bayes inference method (Yang, Wong & Nielsen, 2005). We further tested for sites that experienced episodic events of positive selection by using MEME. This model considers that ω varies between sites (fixed effect), and between branches at a site (random effect) (Murrell et al., 2012). The null expectation is that all branches have ω < 1. In short, this model allows each site to have its own selection history, contrary to fixed effect models (as the ones implemented in CODEML) that assume constant selective pressures within a branch (Murrell et al., 2012).
To assume selection as the main evolutionary mechanism responsible for changes in genetic diversity, sites under positive selection were concatenated for downstream analyses (Positively Selected Sites, PSS). Theoretically, sites interacting with parasite-derived antigens are expected to be under positive selection while other sites, within the same exon, may evolve differently but still maintain the integrity of the MHC molecules. Therefore, we also concatenated the remaining sites (nPSS) and performed identical analyses as for the PSS. Nucleotide mismatch pairwise distributions were calculated for PSS and nPSS of both ''silver eels'' and ''glass eels'' under the assumptions of a constant population size and sudden expansion using DnaSP v5. The historical profile of MHC genetic diversity was investigated with Bayesian skyline plots (BSP) (Drummond et al., 2005) in BEAST v1.8 (Drummond & Rambaut, 2007). Although these statistical procedures are often used to infer demographic events based on fluctuations of neutral genetic diversities, they have also been used to estimate the strength of adaptive evolution at the organism level (Bedford, Cobey & Pascual, 2011), and of functional genomic regions (Padhi & Verghese, 2008).
Because of several expected characteristics of the MHC, including a deviation from a neutral mode of evolution, recombination or gene conversion events (Spurgin et al., 2011) and trans-species polymorphism (Lenz et al., 2013) that also occurs in this species (Bracamonte, Baltazar-Soares & Eizaguirre, 2015) as well as the existence at least three loci (Bracamonte, Baltazar-Soares & Eizaguirre, 2015), we did not attempt to associate the substitution rate to a clock-calibrated evolution. As such, we fixed a molecular clock and assumed three different mutation rates: 0.2, 1 (the default parameter) and 5 substitutions per time unit respectively. The substitution model was chosen in jModeltest (Tamura-Nei: Tn93) (Darriba et al., 2012;Guindon & Gascuel, 2003) and also inserted as parameter in BEAST's runs. Markov chain run was set to a length of 1 x 10 8 . The historical profile of MHC diversity was reconstructed for PSS and nPSS of both ''silver eels'' and ''glass eels.'' Piecewise constant skyline model allowing for five skyline groups were used in three independent MCMC runs to verify consistence in parameter space. We then compared the marginal probability distributions of several parameters amongst the runs. Lastly, we constructed lineage-through-time plots. These plots reflect accumulation of lineages through time translated for a given dated phylogeny (Nee, Mooers & Harvey, 1992).

Molecular indices, population structure and demography amongst sampling locations
Analyses of 355 bp of the mtDNA ND5 in 683 European eels revealed 102 haplotypes including 73 singletons (Data S1). Forty-eight randomly picked singletons were verified by independent extraction, amplification and re-sequencing to eliminate possible risks of sequencing errors. After we eliminated the possibility of sequencing errors (48 out of 48 singletons were verified by independent sequencing) we included all 73 singletons in the subsequent analyses. Amongst sampling locations, haplotype diversity ranged between 0.575 (BU) and 0.934 (GL). Nucleotide diversity ranged between 0.003 (G_SPA) and 0.008 (GL), with an average of 0.005 (±0.001) amongst sampling locations (Table 1). Pairwise F ST comparisons computed from haplotype frequencies amongst the 26 geographically confined groups revealed 17 significant pairwise differentiations however none passed corrections for multiple tests following the false discovery rate (threshold p = 0.013 for 26 tests (Narum, 2006)). All results are shown in Table S1.
Tajima's D values were negative amongst almost all sampled locations, suggestive of population expansion or of population subdivision (Tajima, 1989). The three exceptions were GER (Germany), which belongs to a closed system, the Schwentine river, where recruitment is solely mediated by stocking (Prigge et al., 2013), as well as G_TITA (Italy) and G_WENG (England), where the values might reflect artificial stochasticity due to low sample sizes. Mismatch distribution analyses performed at the population level showed the typical pattern of a historical population expansion where the peak differs from zero (Rogers & Harpending, 1992) (Fig. 1A).

Molecular indices, population structure and demography between generations
Haplotype diversity and genetic diversity between ''silver eels'' and ''glass eels'' were very similar: Hd silver eels = 0.821, Hd glass eels = 0.842; π silver eels = 0.0048, π glass eels = 0.0049. No Figure 1 Location-specific demography accessed with neutral evolving markers. Two types of analyses specific to each marker are presented in this figure: mismatch pairwise distribution of mtDNA's single nucleotide polymorphisms within locations (A), and frequency distribution of allelic states across microsatellite loci (B). The locations here exposed are G_AD2010 (cohort of 2010 from Adour), BU (Burrishole), SLC (LoughComber), Q (Quoile) and were chosen among the others in order to show the variety of demographic signatures produced by mtDNA data. Therefore in mismatch pairwise distribution graphs, full blue lines represent expected distribution under sudden population expansion, full red lines represent expected distribution under constant population size and dotted lines the observed distribution. The xaxis shows the number of mismatches and y-axis its frequency. It is possible to observe the signature of an expanding population in G_AD2010, stable population or recovery from bottleneck BU and SLC, and stable population in Q. Allele frequency distribution plots obtained from same locations show the signature of non-bottlenecked population (B). In these graphs, the bars correspond to allele frequencies, x-axis corresponds to allele frequency classes and y-plots to number of alleles. (continued on next page. . . )

Figure 1 (...continued)
evidence for temporal genetic structure based on haplotype frequency distributions was detected between these groups, F ST = 0.000, p = 0.138. Both groups also had negative and significant Tajima indices: D silver eels = −2.053, D glass eels = −2.357, both p < 0.05 (Table 2). Mismatch distribution analyses revealed that both ''silver eels'' and ''glass eels'' display the distribution of expanding populations, suggesting that the overall pattern is not driven by a single generation and has a true biological origin visible in both cohorts (Fig. 2).

Molecular indices, population structure and demography amongst locations
Across populations, He ranged between 0.687 (G_NIRL) and 0.758 (PT), Ho between 0.557 (G_TITA) and 0.667 (DK) and the average number of alleles per locus varied between 3.546 (G_VFRA) and 15.773 (G_AD2012). Allelic frequencies are reported in Data S2. F IS varied between 0.0380 (G_NIRL) and 0.254 (Q) ( Table S1). Even though F ST estimates are very low, pairwise comparisons revealed 7 statistically significant pairwise comparisons after correcting for multiple testing (Table S1). Since six of those included comparisons between locations with low sample size, i .e. G_NIRL, G_WENG, G_TITA or G_BNIRL, the significance could be attributed to stochasticity due to low sample sizes. STRUCTURE analyses did not show any signs of population clustering as expected under the weak observed differentiation (Fig. S3). None of the sampled locations showed either heterozygote excess or a mode shift in allele frequencies, genetic signatures characteristic of population bottlenecks (Fig. 1B).
None of the replicates showed evidence of heterozygote excess or deficiency. Averaged allele frequencies of neither ''silver eels'' nor ''glass eels'' deviated from an expected L-shape distribution (Fig. S3). However, we found that the averaged allele frequencies distribution observed in the ''silver eels'' group showed the signature of a 5-generations-old bottleneck identified from computer simulations (Luikart et al., 1998). This is particularly evident in the distribution of the two most common allele classes, 0.8-0.9 and 0.9-1.0. This signature was not visible anymore in the ''glass eel'' group (Fig. S4).
Estimates of effective population size (N e ) amongst replicates ranged between 0-625.2 for ''silver eels'' and 0-2708.9 for ''glass eels''. The harmonic mean of effective population size estimates amongst replicates,N e , resulted in 480.9 <N e silver eels < 2941.7 and 1380.2 <N e glass eels < 3506.0 (Table S2).

Molecular indices and population structure
We sequenced a 247 bp fragment of the exon 2 of the MHC class II region (91% of the total size of the exon) in 327 individuals using 454 sequencing technology. We detected 229 different amino acid coding variants. Among those, 226 (98%) were found to be unique in the dataset but present in both independent replicated reactions and therefore kept as true variants (Data S3). A total of 116.276 sequence reads were used in this study. Amongst locations, MHC nucleotide diversity ranged between 0.102 (LL) and 0.138 (BT). The mean number of alleles per individual ranged between 2 (Q, SE = 0.298) and 4 (G_BU, SE = 0.392) (Table S3) and revealed to overall significantly differ amongst sampled locations (F 17 = 1.674, p = 0.046). However, post-hoc pairwise comparisons showed no significant differences between pairs of populations after correction for multiple testing (all p > 0.05). The mean nucleotide divergence (p-distance) ranged between 0.078 (BL) and 0.141 (FI) ( Table S4) and was significantly different among sampled locations (F 17 = 1.860, p = 0.021). Post-hoc pairwise comparisons revealed two significant comparisons after correction for multiple testing (GER vs FI, t = −3.551, p = 0.045, GER vs G_AD2011, t = −3.961, p = 0.010), suggesting a reduced MHC diversity the stocked freshwater system of the Schwentine river, in Germany.

Screening for novel genetic diversity through events of gene conversion
Gene conversion segments with an average nucleotide length of 4 bp were detected within the ''glass eels'' group but not in the ''silver eel'' group. The average ψ of the whole segment was found to be 0.0002 (Fig. 3). This value is high enough to ascertain the occurrence of conversion events, but not robust enough to determine the exact length of the observed tracts (Betran et al., 1997).

Testing for positive selection
Model-based tests using CODEML revealed 11 sites under positive selection while MEME identified 27 sites that have experienced episodic events of positive selection (Table S4). The discrepancies between the two methods reflect the different assumptions underlying the fixed effect models implemented in CODEML and the mixed effect models of MEME. Positively selected sites detected by both methods matched 10 out of 19 antigen binding sites identified in humans by X-ray crystallography (Reche & Reinherz, 2003) (Fig. 3). Due to the functional role of the MHC, all amino acid sites that have experienced at least episodic events of selection were selected for further analyses (Fig. 3)

Demography and historical profile of MHC genetic diversity
All mismatch distributions indicated a clear deviation from a constant population size, fitting a scenario where a major demographic event occurred (Figs. 4 and 5). The frequency distribution of pairwise differences showed different peaks for both PSS (Fig. 4)   (PSS pairwise differences = 20; nPSS pairwise differences = 10) (Fig. 5). Those peaks reflect old lineages that are maintained in genes exhibiting trans-species polymorphism as is expected of the MHC (Klein, Sato & Nikolaidis, 2007) and recently reported for the European eel (Bracamonte, Baltazar-Soares & Eizaguirre, 2015). It was also possible to observe peaks in PSS and nPSS in the frequency of pairwise differences equaling 1. Those peaks are suggestive of increases in genetic diversity. No differences in the demographic profiles were detected between ''silver eels'' and ''glass eels' ' (Figs. 4 and 5).  Bayesian demographic reconstructions revealed a steep decline in genetic diversity that occurred close to the present time. This pattern is characteristic of a genetic bottleneck and is shared by all reconstructions independently of the substitution rates (Fig. 4). However, lineage-through-time plots reveal a very recent burst of lineage diversification, also at t = 0 and common to both generations (Fig. 6). The three independent MCMC runs clearly overlap for the distributions of the posterior, likelihood and skyline estimates (Fig. S5), assuring that the profiles observed were not a product of the Bayesian stochasticity, but rather a real and reproducible pattern.

DISCUSSION
We investigated how the steep decline in European eel recruitment observed in the 1980s and subsequent population reduction may have affected this species' genetic diversity. We expanded on previous studies (Pujolar et al., 2011;Pujolar et al., 2013;Wirth & Bernatchez, 2003) and searched for contemporary signature of a population reduction by considering neutral markers (mtDNA or microsatellites) as well as a highly polymorphic region of the immune genes of the MHC. Although MHC is an excellent marker to evaluate genetic diversity, as a proxy for adaptive potential, in endangered populations (Sommer, 2005;Sutton, Robertson & Jamieson, 2015), no formal study of its variation and evolution exist for the European eel.

Location-specific patterns of genetic diversity and demographic estimates
Although it was not the primary focus of this study, spatial comparison of mtDNA genetic diversity provided information on the genetic differentiation amongst continental locations. Under neutrality, haplotype and nucleotide diversities are predicted to be a function of population size (Frankham, Briscoe & Ballou, 2002). Therefore, in the drastically declined European eel population, we expected to find an overall low genetic diversity. Furthermore, and due to reports of panmixia (Als et al., 2011), we also expected coherent patterns amongst sampled locations. Instead, we found variation in nucleotide diversity (0.003-0.008), haplotype diversity (0.575-0.934) and Tajima's D estimates amongst locations. The high variation in genetic indices amongst geographical areas suggests that processes act differently across the continental distribution of the European eel. Our observations may result from the post-hatching transatlantic migration, since simulations showed variability in spawning grounds would leave genetic signatures across continental locations under low recruitment (Baltazar-Soares et al., 2014). Conversely, it could be explained by a scenario where mtDNA haplotypes are linked to genes under single-generation local selection (Pujolar et al., 2014). In this scenario, selection would act mainly in the local foraging environment, and not in the spawning ground, with specific pressures sorting out genotypes in given locations. Expanding the study towards a more genomic approach with adult fish sampled from the spawning ground would reveal further insights into the most prominent scenario.
Due to high levels of polymorphism, neutrally evolving microsatellites are thought to be sensitive enough to detect subtle shifts in population dynamics (England, Luikart & Waples, 2010). Here as well, we hypothesized that the chronically low recruitment in European eel observed since the 1980s had negatively affected estimates of genetic diversity in a coherent spatial pattern. Estimates of allelic richness (Ar = 2.710-2.960) and heterozygosity estimates (He = 0.687-0.758) were very similar, and neither mode shifts nor heterozygote excesses were observed (Fig. 1, Fig. S1). These results are in line with a previous study (Pujolar et al., 2011), which was conducted focusing on 12 locations (3 locations sampled across a temporal range) and that employed 22 EST-linked microsatellites. The apparent homogeneity of the allelic indices amongst locations matches the expectations based on neutral nuclear markers for a panmictic population where successfully recruited mature fish would mate randomly in the spawning ground (Als et al., 2011). Note, after correction for multiple testing, we detected significantly different pairwise F ST estimates amongst some locations. Although this could suggest deviations from panmixia, it is likely this pattern is linked to stochasticity due to low sample sizes in those locations (G_NIRL, G_TITA and G_WENG).
Regarding the MHC, we found differences in the mean number of alleles amongst locations although none of the post-hoc pairwise comparisons revealed to be significant after correction for multiple testing. MHC allele pool composition did not vary amongst sampled locations. We found that the mean nucleotide distances vary amongst locations, but only comparisons between the German populations and two other sampled locations revealed to be significant after correction for multiple testing (Germany vs Finland, Germany vs France 2011). Whether it relates to the fact that those samples correspond to a freshwater system, the Schwentine in Germany, where all eels are stocked remains to be investigated. Stocking however might be an issue in other sampling areas as well and therefore understanding why this population demonstrates a reduction in diversity should be the focus of further studies.

Recruitment decline of the 1980s did not affected genetic estimates of neutral evolving markers
The major objective of this study was to compare patterns of temporal variation of genetic diversity post-recruitment collapse observed in the 1980s. Here as well, our mitochondrial DNA results showed (i) no evidence for a genetic bottleneck, (ii) no differences in haplotype and nucleotide diversities between ''silver eels '' and ''glass eels,'' and (iii) no signature of bottleneck in the frequency distribution of pairwise mismatches. The later rather points towards an historical population expansion (Rogers & Harpending, 1992). Given the low variability of the mtDNA marker-in comparison with the set of highly polymorphic microsatellites and MHC gene-we suggest that such an event extends back in time to a scenario of expansion related to ice-sheet retreat after the last glacial maxima, as previously proposed for such a pattern (Jacobsen et al., 2014).
Investigating the distribution of the genetic variance observed at microsatellites, we found significant differentiation between ''silver eels'' and ''glass eels'' replicates. This pattern of genetic variance distribution is in line with previous reports that also attributed higher genetic variance amongst temporal, rather than spatial, partitions of A. anguilla along the European coasts (Dannewitz et al., 2005). It also confirmed our a-priori assumption of each group representing a distinct generation, and therefore excludes possible confounding factors associated with overlapping generations from the interpretation of demographic estimates (Cornuet & Luikart, 1996;Waples & Do, 2010). Note, even though less likely, we cannot exclude that this observed structure also relates to a spatially structured spawning area of this species (Baltazar-Soares et al., 2014;Dannewitz et al., 2005).
Using microsatellites, we detected no evidence of heterozygote excess in any of the replicates, nor any differences between allelic richness of ''silver eels'' and ''glass eels.'' This supports previous reports that the recruitment collapse and subsequent low abundance of the eel population did not leave the expected genetic signatures of reduced genetic diversity (Pujolar et al., 2011) suggestive of a system which replenishes genetic diversity rapidly.
However, in-depth demographic analyses suggest that the eel effective population size might actually not be stable. Several lines of evidence support this interpretation: firstly, we estimated ∼20% higher effective population size in ''glass eel'' replicates (harmonic mean N e = 3506.0) compared to ''silver eels'' (N e = 2941.7). These contemporary estimates are near the lower confidence intervals of historic effective population sizes previously reported (5,000 < N e < 10,000; (Wirth & Bernatchez, 2003), but within contemporary estimates of 3,000 < N e < 12,000 (Pujolar et al., 2011)). Noteworthy, the confidence intervals calculated in this study for each generation overlap, raising the need to cautiously interpret those results. Secondly, we observed fewer alleles in the most frequent class of allele frequencies in ''silver eels''. The apparent reduction of the most frequent allele class may suggest that the species demography is experiencing a transitory stage from a severe bottleneck, partly detected in ''silver eels.'' It is important to mention that such a signature would only be detected under a severe population reduction a few generations in the past. Hence, we could speculate that the hypothetical bottleneck detected only in ''silver eels'' may relate to the drops in European eel recruitment that occurred in the beginning of the 1960s (EIFAAC/ICES, 2011). It is possible that the 1960s low recruitment had a major impact on the overall genetic diversity of the species. By the crash in the 1980s, the population would have already been depleted from its original genetic diversity, at least for neutrally evolving markers. Such a scenario requires further studies to be confirmed and would rely on historical samples to be analyzed.

MHC reveals signatures of selection
Extending the evaluation of genetic diversity to the evolutionary analysis of the adaptive immune genes of the MHC was motivated by two main reasons. The first relates to studies suggesting MHC diversity to be more sensitive than neutrally evolving markers in the detection of demographic shifts (Sommer, 2005;Sutton et al., 2011). The second relates to the invasion of European freshwater systems by the nematode parasite, Anguillicola crassus, for which the MHC was found to respond to in the paratenic host, the three-spined stickleback (Eizaguirre et al., 2012b). Using the exon 2 of the MHC class II β gene, we evaluated (1) genetic diversity, which might have been affected by the recruitment collapse and subsequent population reduction and (2) allele frequency shifts between generations which would be a signature consistent with directional parasite-mediated selection.
Using next-generation sequencing, we identified a total of 229 MHC alleles amongst 327 individuals. This indicates that the diversity within this species is not low and directly compares to observations made in wild populations of other fishes that are not qualified as endangered, such as for instances, the half-smooth tongue sole (88 MHC class II alleles amongst 160 individuals) (Du et al., 2011). Of note, the characterization of the MHC class II genes in this species revealed that up to six different alleles may exist per individual, suggesting the presence of at least three loci (Bracamonte, Baltazar-Soares & Eizaguirre, 2015;Bracamonte, 2013). Because next generation sequencing is thought to generally overestimate the number of MHC alleles detected (Babik et al., 2009;Lighten, Oosterhout & Bentzen, 2014;Sommer, Courtiol & Mazzoni, 2013), we took multiple precautions to avoid artifacts (reconditioning steps, reduced numbers of PCR cycles, duplicates). Despite such precautions and sequence confirmation using cloning (performed in Bracamonte, Baltazar-Soares & Eizaguirre, 2015), we cannot exclude that some variants were called alleles even though artifactual (Sommer, Courtiol & Mazzoni, 2013). Again, as mtDNA sequencing revealed a large number of haplotypes (N = 102), the large diversity at the MHC displayed in this species may not be surprising.
Generally, our results are suggestive of a pattern of selective sweep at the MHC between the two generations examined here. We found ''silver eels'' to exhibit lower mean number of alleles and lower mean nucleotide distance than ''glass eels,'' suggesting that the ''silver eel'' generation was under a selective pressure that reduced its pool of MHC alleles to fewer and more similar alleles. This observation suggests that either the selective pressure was widespread amongst continental locations or that it acted when all eels experienced similar conditions; as for instance, during the fastening spawning migration. A hypothetical selective pressure imposed by A. crassus meets both criteria. Not only is this parasite ubiquitous in European freshwater systems (Wielgoss et al., 2008) but it also impairs the swimming performance of infected eels (Palstra et al., 2007).
Overall, the higher genetic diversity of ''glass eels'' together with the identification of two short gene conversions in this study suggest an ongoing regeneration of the species immune adaptive potential. Indeed, gene conversation is a mechanism capable of generating novel diversity within the MHC region, and seems to be a predominant mechanism in genetically depauperate populations (Spurgin et al., 2011).
From our results we can hypothesize two scenarios. A first scenario involves a link between a sudden reduction in population size, a loss of genetic diversity and a constant selective pressure extending after the bottleneck. In this scenario, genetic drift would affect overall genetic diversity but since selection would continue to act, genetic diversity of positively selected regions would remain low (Eimes et al., 2011). A second scenario relates to multiple MHC loci carrying similar alleles due to recent duplications and to the hypothesis that a population would experience a size reduction and an event of selection within the same time frame. In this scenario, the selective pressure through the bottleneck would lead to a faster fixation of resistant alleles (Eimes et al., 2011).

Limits of the study
Evaluating demography in the European eel is complex particularly due to the difficulty of sampling them at their mating ground and the reliance on indirect genetic evidence. Even though patterns of selection and recent recovery of the MHC diversity seem robust, due to the high allelic variation reported in this gene, we acknowledge that other coalescence models could revealed refined patterns (Árnason & Halldórsdóttir, 2015;Wakeley, 2013). It is indeed possible that multiple loci analysed together as one locus-although with same functional basis-could create a pattern similar to multiple coalescence. To the best of our knowledge, multiple merger models such as those described by (Árnason & Halldórsdóttir (2015) have not yet been applied to investigate the evolutionary signature that linked copies of the same (functional) gene produce on the evaluation of genetic diversity and demography. This thus limits working hypotheses and would be difficult to interpret. Nonetheless, while further investigations are obviously needed to clarify fluctuations in genetic diversity in such a complex but evolutionary relevant immune gene, we argue that our work represents a first empirical step along this line of research.

Conclusions
In summary, our work reveals signatures of recent reduction in MHC genetic diversity and suggests signs of ongoing recovery of this gene's diversity contributing to the immunogenetic adaptive potential of this endangered species (Radwan, Biedrzycka & Babik, 2010;Sommer, 2005;Stiebens et al., 2013b). Future research will be needed to provide conclusive evidence as to which scenario holds, accommodating also newer theories to further verify the validity of our findings. A future perspective would be to extend the time-series analyses by incorporating screening of MHC diversity in ongoing monitoring practices. This would be a valuable approach to access the evolution of the species adaptive potential.