Magally Torres-Leguizamon

Je ´ro ˆme Mathieu

Thibaud Decae ¨ns

Lise Dupont
email: lise.dupont@u-pec.fr

Genetic Structure of Earthworm Populations at a Regional Scale: Inferences from Mitochondrial and Microsatellite Molecular Markers in Aporrectodea icterica (Savigny 1826)

Despite the fundamental role that soil invertebrates (e.g. earthworms) play in soil ecosystems, the magnitude of their spatial genetic variation is still largely unknown and only a few studies have investigated the population genetic structure of these organisms. Here, we investigated the genetic structure of seven populations of a common endogeic earthworm (Aporrectodea icterica) sampled in northern France to explore how historical species range changes, microevolutionary processes and human activities interact in shaping genetic variation at a regional scale. Because combining markers with distinct modes of inheritance can provide extra, complementary information on gene flow, we compared the patterns of genetic structure revealed using nuclear (7 microsatellite loci) and mitochondrial markers (COI). Both types of markers indicated low genetic polymorphism compared to other earthworm species, a result that can be attributed to ancient bottlenecks, for instance due to species isolation in southern refugia during the ice ages with subsequent expansion toward northern Europe. Historical events can also be responsible for the existence of two divergent, but randomly interbreeding mitochondrial lineages within all study populations. In addition, the comparison of observed heterozygosity among microsatellite loci and heterozygosity expected under mutation-drift equilibrium suggested a recent decrease in effective size in some populations that could be due to contemporary events such as habitat fragmentation. The absence of relationship between geographic and genetic distances estimated from microsatellite allele frequency data also suggested that dispersal is haphazard and that human activities favour passive dispersal among geographically distant populations.

Introduction

Earthworms represent one of the largest reservoirs of animal biomass and the main invertebrate group of soil ecosystem engineers in most terrestrial temperate ecosystems [START_REF] Lavelle | Soil Ecology[END_REF]. They play a key role in soil functioning: they relocate surface litter or organic matter throughout the soil profile [START_REF] Edwards | The importance of earthworms as key representatives of the soil fauna[END_REF][START_REF] Gorres | Soil micropore structure and carbon mineralization in burrows and casts of an anecic earthworm (Lumbricus terrestris)[END_REF], affect microbial activity [START_REF] Scheu | Effects of the presence and community composition of earthworms on microbial community functioning[END_REF], and have a significant effect on organic matter mineralisation and soil biogeochemical cycles [START_REF] Postma-Blaauw | Earthworm species composition affects the soil bacterial community and net nitrogen mineralization[END_REF]. They modify soil structure via the construction of biogenic aggregates and galleries [START_REF] Lavelle | Functional domains in soils[END_REF], resulting in differences in aeration and drainage [START_REF] Bouma | Influence of Soil Macroporosity on Environmental-Quality[END_REF]. Earthworms also influence plant growth [START_REF] Brown | How earthworms affect plant growth: burrowing into the mechanisms[END_REF] and plant community structure [START_REF] Eisenhauer | Impact of above-and below-ground invertebrates on temporal and spatial stability of grassland of different diversity[END_REF][START_REF] Laossi | Effects of plant diversity on plant biomass production and soil macrofauna in Amazonian pastures[END_REF], and can be used as indicators of habitat quality [START_REF] Mathieu | Habitat quality, conspecific density, and habitat pre-use affect the dispersal behaviour of two earthworm species, Aporrectodea icterica and Dendrobaena veneta, in a mesocosm experiment[END_REF][START_REF] Ruiz | IBQS: A synthetic index of soil quality based on soil macro-invertebrate communities[END_REF] and as biomarkers in toxicity tests [START_REF] Capowiez | Earthworm cast production as a new behavioural biomarker for toxicity testing[END_REF].

Despite their fundamental impact on soil ecosystems, the spatial population dynamics of earthworms is poorly understood. In particular, there is little information on the amount and spatial distribution of genetic variation in earthworm populations. Few studies have simultaneously investigated the influence of historical events, such as glacial periods, and contemporary processes, such habitat fragmentation, on the genetic diversity of these species. In their review on the genetic structure of soil invertebrate populations, Costa et al. [START_REF] Costa | Genetic structure of soil invertebrate populations: Collembolans, earthworms and isopods[END_REF] cite only seven studies of earthworm populations. They conclude that earthworm populations generally show a complicated pattern of gene flow, with a weak relationship between genetic and geographic distances. Population genetic structure of earthworms is therefore likely to be strongly influenced by human activities. In an agricultural landscape, the spatial distribution of earthworm species is expected to be fragmented, with patches of suitable habitat being separated by large areas of unsuitable habitat. Furthermore, it has been shown that earthworms are negatively affected by the intensity of agriculture [START_REF] Ponge | The impact of agricultural practices on soil biota: A regional study[END_REF] and in particular by the use of pesticides [START_REF] Pelosi | Reduction of pesticide use can increase earthworm populations in wheat crops in a European temperate region[END_REF]. The consequences of landscape spatial structure on genetic diversity depend on the rate at which individuals move between patches of suitable habitat. In particular, restricted dispersers, such as earthworms, are likely to be prone to local extinction due to stochastic processes [START_REF] Amos | Factors affecting levels of genetic diversity in natural populations[END_REF]. However, it has also been suggested that the rate of gene flow and the amount of genetic variation may actually increase as habitats become more fragmented [START_REF] Young | The effects of forest fragmentation on genetic variation in Acer Saccharum Marsh (sugar maple) populations[END_REF]. More earthworm population genetics studies are needed to determine (i) how earthworms move between patches, (ii) how spatial structure affects the stability and dynamics of spatially structured earthworm populations, and (iii) how the landscape affects genetic diversity.

Our model earthworm species, Aporrectodea icterica is an abundant species commonly found in agricultural soils [START_REF] Capowiez | Lethal and sublethal effects of imidacloprid on two earthworm species (Aporrectodea nocturna and Allolobophora icterica)[END_REF]. It belongs to the endogeic ecological type (i.e. species living and foraging in or immediately below the rhizosphere making horizontal burrows through the soil to move around and to feed) [START_REF] Bouche | Lombriciens de France[END_REF], although A. icterica is also believed to feed at least partly on leaf litter [START_REF] Bernier | Earthworm feeding activity and development of the humus profile[END_REF]. This diploid, obligatory bi-parental species [START_REF] Casellato | On polyploidy in oligochaetes with particular reference to Lumbricids[END_REF] is native to the temperate zones of Europe [START_REF] Bouche | Lombriciens de France[END_REF], but is an invasive species in North America [START_REF] Schwert | The first North American record of Aporrectodea icterica (Savigny, 1826) (Oligochaeta, Lumbricidae), with observations on the colonization of exotic earthworm species in Canada[END_REF]. Its taxonomic status is firmly grounded and the species has distinct morphology making it easy to recognise. Its dispersal behaviour has been studied in laboratory [START_REF] Mathieu | Habitat quality, conspecific density, and habitat pre-use affect the dispersal behaviour of two earthworm species, Aporrectodea icterica and Dendrobaena veneta, in a mesocosm experiment[END_REF][START_REF] Caro | Are dispersal behaviours of earthworms related to their functional group?[END_REF] and it has been used in ecotoxicological studies [START_REF] Capowiez | Lethal and sublethal effects of imidacloprid on two earthworm species (Aporrectodea nocturna and Allolobophora icterica)[END_REF][START_REF] Capowiez | Sublethal effects of imidacloprid on the burrowing behaviour of two earthworm species: Modifications of the 3D burrow systems in artificial cores and consequences on gas diffusion in soil[END_REF][START_REF] Capowiez | Assessment of the effects of imidacloprid on the behavior of two earthworm species (Aporrectodea nocturna and Allolobophora icterica) using 2D terraria[END_REF]. At the genetic level, a recent study of two A. icterica populations revealed the existence of two mitochondrial lineages with divergence values ranging from 10% to 11% [START_REF] Torres-Leguizamon | Isolation of polymorphic microsatellite markers in Aporrectodea icterica (Savigny 1826)[END_REF]. Such highly divergent mitochondrial lineages have been reported in several other earthworm morphospecies (e.g. [START_REF] James | DNA Barcoding reveals cryptic diversity in Lumbricus terrestris L. 1758 (Clitellata): resurrection of L. herculeus[END_REF][START_REF] King | Opening a can of worms: unprecedented sympatric cryptic diversity within British lumbricid earthworms[END_REF][START_REF] Novo | High genetic divergence of hormogastrid earthworms (Annelida, Oligochaeta) in the central Iberian Peninsula: evolutionary and demographic implications[END_REF]). In A. icterica, nuclear analysis indicates that the two lineages interbreed [START_REF] Torres-Leguizamon | Isolation of polymorphic microsatellite markers in Aporrectodea icterica (Savigny 1826)[END_REF], demonstrating that they belong to the same species.

Deep mitochondrial divergences within morphospecies can be attributed to population isolation within distinct periglacial refugia [START_REF] King | Opening a can of worms: unprecedented sympatric cryptic diversity within British lumbricid earthworms[END_REF]. When the divergent lineages were found in sympatry, such as in A. icterica, it was suggested that lineages came into contact and mixed during recolonisation, during the warmer interglacial periods [START_REF] King | Opening a can of worms: unprecedented sympatric cryptic diversity within British lumbricid earthworms[END_REF]. Given the low vagility of earthworms, we hypothesise that this mixing is in large part due to human activities which have accelerated the rate of organism dispersal, and brought previously allopatric species into contact [START_REF] Wolf | Predicting the risk of extinction through hybridization[END_REF]. During this secondary contact, weak reproductive barriers between lineages and fertile hybrids with little or no reduction in fitness can lead to genetic assimilation and loss of genetic distinctness between the hybridizing lineages, and the possible extinction of one or both parental lineages [START_REF] Wolf | Predicting the risk of extinction through hybridization[END_REF]. For recent or in progress hybridization events, introgressed mitochondrial and nuclear genes are predicted to display cytonuclear disequilibrium [START_REF] Ballard | The incomplete natural history of mitochondria[END_REF] (i.e. non-random association of alleles or genotypes at a nuclear locus with haplotypes of cytoplasmically inherited organellar DNA [START_REF] Asmussen | Sampling Theory for Cytonuclear Disequilibria[END_REF]).

Here, we focus on the genetic structure of seven earthworm populations sampled at a regional scale (,100 km 2 ), comparing the spatial regional patterns obtained using mitochondrial (mtDNA) and nuclear (nDNA) molecular markers. We discuss the role of evolutionary forces including genetic drift and contemporary gene flow, large-scale landscape changes (e.g. glacial periods) and anthropogenic effects in structuring the genetic diversity and in the differentiation of populations.

Materials and Methods

Sampling and DNA extraction

In April 2010, 218 A. icterica individuals were collected from seven populations in Normandy (northern France). Six populations were located in .6 year-old pastures (on average, clay = 16%, silt = 64% and sand = 20%, mean pH: 6.1, C: 23 g.kg 21 , N: 2.3 g.kg 21 , C/N: 10), within a distance of 3 to 10 km from the city of Yvetot (I03, I07, I19, I20, I25 and I27). Each proprietor gave his agreement for sampling to J. Mathieu who should be contacted for future permissions. These pastures are grazed by dairy cattle from mid-March to mid-September with a stocking rate of 2-6 animal units ha 21 depending on the season, and spread with cow manure each year. Plant cover consisted mainly in Festuca elator L., Phleum pratense, Trifolium repens L., and Lolium. The seventh population was located 35 km away from Yvetot (IR), near the University of Rouen in a location for which no specific permission was required (Fig. 1 and Table 1 for GPS coordinates). This field study did not involve endangered or protected species.

Individuals were preserved in pure ethanol for DNA analysis. Total genomic DNA was extracted from a segment of the anterior end of the earthworm using the CTAB extraction protocol: digestion using proteinase K, followed by protein precipitation with CTAB, a chloroform:isoamyl alcohol (24:1) wash and DNA precipitation with sodium acetate (3 M) and ethanol.

Mitochondrial DNA amplification and sequencing

A fragment of the cytochrome c oxidase subunit I mitochondrial gene (COI) was amplified and sequenced using the universal primers LCO1490 and HCO2198 [START_REF] Folmer | DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates[END_REF]. For the I20 and IR populations, sequences were taken from Torres-Leguizamon et al. [START_REF] Torres-Leguizamon | Isolation of polymorphic microsatellite markers in Aporrectodea icterica (Savigny 1826)[END_REF] (GenBank accession numbers JN381881-JN381930). Each amplification mixture (25 ml) contained 10 ng DNA, 12.5 ml of Taq PCR Master Mix (Qiagen, Hilden, Germany) and 0.25 mM of each primer. Polymerase chain reactions were performed using an initial denaturation step at 94uC for 3 min, followed by 40 cycles of the three following steps: denaturation at 94uC for 30 s, annealing at 49uC for 1 min and extension at 72uC for 1 min 30 s. The final extension was done at 72uC for 10 min. PCR products were purified using Microclean (Microzone Limited, Haywards Heath, UK). Both strands of amplicons were sequenced using Big-Dye Terminator Cycle sequencing kit version 1.1 (Applied Biosystems, Foster City, CA, USA) according to the standard protocol used in the genomic platform at the Mondor Institute of Biomedical Research (Cre ´teil, France). Sequences were deposited in GenBank (Accession numbers: KF856627-KF856710).

Mitochondrial DNA sequence analysis COI sequences were aligned manually using BioEdit v. 7.0.5.3 [START_REF] Hall | BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT[END_REF]. For each population, haplotype diversity (h) and nucleotide diversity (p) were estimated using DNAsp v. 5.10 software [START_REF] Librado | DnaSP v5: a software for comprehensive analysis of DNA polymorphism data[END_REF]. Haplotypic richness after rarefaction to a population size of 11 individuals was estimated using Contrib software [START_REF] Petit | Identifying populations for conservation on the basis of genetic markers[END_REF]. Departure from neutrality was tested using Fu's Fs [START_REF] Fu | Statistical tests of neutrality against population growth, hitchhiking and background selection[END_REF] and Ramos-Onsins and Rozas' R2 statistic [START_REF] Ramos-Onsins | Statistical properties of new neutrality tests against population growth[END_REF], which are powerful tests for detecting recent population expansion under assumptions of neutrality. R2 is appropriate for small sample sizes [START_REF] Ramos-Onsins | Statistical properties of new neutrality tests against population growth[END_REF]. The significance of R2 and Fs were evaluated by comparing their observed values with their null distribution, generated by 10 000 random replicates using the empirical population sample size and observed number of segregating sites implemented by DnaSP v. 5.10. [START_REF] Librado | DnaSP v5: a software for comprehensive analysis of DNA polymorphism data[END_REF].

To describe the phylogenetic relationships between haplotypes, a statistical parsimony network was constructed using TCS v. 1.21 [START_REF] Clement | TCS: a computer program to estimate gene genealogies[END_REF]. The divergence among haplotypes was calculated in MEGA 5 using the mean uncorrected p-distance [START_REF] Tamura | MEGA5: Molecular Evolutionary Genetics Analysis Using Maximum Likelihood, Evolutionary Distance, and Maximum Parsimony Methods[END_REF].

To investigate the regional structure of A. icterica populations at the mitochondrial level, the overall genetic differentiation between populations was first estimated by calculating the global W st in Arlequin v. 3.5.1.2 software [START_REF] Excoffier | Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows[END_REF]. We then performed a spatial analysis of molecular variance using SAMOVA software [START_REF] Dupanloup | A simulated annealing approach to define the genetic structure of populations[END_REF]. Clusters are identified based on geographic proximity and genetic homogeneity [START_REF] Dupanloup | A simulated annealing approach to define the genetic structure of populations[END_REF]. The geographic coordinates of each sampled locality were used as spatial information. Simulations were conducted with 'K' (number of groups) ranging from two to seven and each simulation annealing process was repeated 100 times.

The clustering giving the highest W ct value, corresponding to the optimal number of groups, was selected.

Microsatellite genotyping

Individuals were genotyped using seven microsatellite loci: Ai45, Ai56, Ai68, PB10D, 2PE40, 2PE70, C4 [START_REF] Torres-Leguizamon | Isolation of polymorphic microsatellite markers in Aporrectodea icterica (Savigny 1826)[END_REF]. For individuals from the I20 and IR populations, genotypes were taken from Torres-Leguizamon et al. [START_REF] Torres-Leguizamon | Isolation of polymorphic microsatellite markers in Aporrectodea icterica (Savigny 1826)[END_REF]. The seven loci were amplified by a touchdown PCR procedure that included an initial denaturation step of 3 min at 94uC, followed by 35 s at 94uC, 45 s at the initial temperature T a +8uC, 10 cycles in which the temperature was decreased by 1uC per cycle, 1 min at 72uC, 25 cycles of 35 s at 94uC (except for PB10D and C4 for which 30 cycles were done), 45 s at Ta 22uC, 1 min at 72uC, with a final elongation step of 10 min 72uC. Each amplification mixture (15 ml) contained 10 ng/ml DNA, 1X reaction buffer (GoTaq Flexi buffer 5X), 2.5 mM of MgCl 2 (except for Ai56 and Ai68 for which 1.5 mM was used), 0.5 mM dNTPs, 0.25 mM of each primer and 0.5 units of GoTaq Flexi DNA polymerase (Promega, Madison, WI, USA). PCR products were loaded on an ABI 310 sequencer along with the LIZ500 size standard; alleles were scored using Genescan software (all from Applied Biosystems, Foster City, CA, USA).

Statistical analysis of microsatellite variation

In each population, the genetic diversity was analysed by computing allelic frequencies, number of alleles (N all ) and unbiased expected heterozygosity (H e ) using Genepop v. 4 software [START_REF] Raymond | GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism[END_REF]. To take into account differences in sample size, allelic richness (A) after rarefaction to a population size of 22 individuals was estimated using FSTAT v. 2.9.3. software [START_REF] Goudet | FSTAT, a program to estimate and test gene diversities and fixation indices[END_REF]. Exact tests for genotypic linkage disequilibria and deviations from Hardy-Weinberg equilibrium (HWE) were computed using Genepop v. 4 [START_REF] Raymond | GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism[END_REF]. The sequential Bonferroni method was applied to adjust for multiple comparisons. Weir and Cockerham's (1984) estimator of the inbreeding coefficient F is was calculated using Genepop v. 4 [START_REF] Raymond | GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism[END_REF]. The presence of null alleles was tested using Micro-Checker v. 2.2.3 software, in which the Oosterhout method [START_REF] Van Oosterhout | MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data[END_REF] was implemented and potential frequency of null alleles was estimated.

We tested for deviation from mutation-drift equilibrium in the study populations using the approach detailed in Cornuet & Luikart [START_REF] Cornuet | Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data[END_REF] and implemented in their software BOTTLENECK v. 1.2.02. Using a Wilcoxon test, observed heterozygosity was compared with the heterozygosity expected under equilibrium considering a two-phase mutation model (TPM) recommended for microsatellite data [START_REF] Dirienzo | Mutational Processes of Simple-Sequence Repeat Loci in Human-Populations[END_REF]. Recently founded populations are expected to show a transient excess of expected gene diversity, whereas expanding populations (e.g. recovering from a bottleneck) or populations resulting from immigration from differentiated sources should show a deficit in expected gene diversity [START_REF] Cornuet | Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data[END_REF].

To investigate the genetic structure among populations, a G-test of allelic differentiation was carried out using Genepop v. 4 [START_REF] Raymond | GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism[END_REF]. In addition, we performed a principal component analysis (PCA) on gene frequency data using PCAGen v. 1.2.1 software (available at http://www2.unil.ch/popgen/softwares/pcagen.htm). The significance of each principal component was assessed from 1000 permutations. We also used the software STRUCTURE v. 2.3.1 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF] to estimate the number of genetic clusters (K) present among all populations. This software generates clusters of individuals based on the Hardy-Weinberg model of genotypic distribution. Simulations were run using the admixture model without prior population information. We modelled cluster assignments for K ranging from 1 to 10. We performed 25 independent runs for each K value to confirm consistency across runs. In all simulations, we applied a burn-in period of 10 000 iterations and 100 000 Markov chain Monte Carlo iterations. To determine the most likely value of K, we used the DK method [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF].

To estimate recent migration rates and test for significant cases of assignment to populations other than the population of origin (i.e. first-generation dispersers) we used the Bayesian method [START_REF] Rannala | Detecting immigration by using multilocus genotypes[END_REF] implemented in GeneClass2 v. 2.0. [START_REF] Piry | GENECLASS2: A software for genetic assignment and first-generation migrant detection[END_REF] paired with a Monte Carlo resampling method for computation of assignment probabilities for each population [START_REF] Paetkau | Genetic assignment methods for the direct, real-time estimation of migration rate: a simulation-based exploration of accuracy and power[END_REF] using 10 000 simulated individuals.

To test for the null hypothesis of independence between genetic and geographic distances, the logarithm of Euclidian geographic distances were plotted against F st /(12F st ) to compute a linear relationship following the recommendations of Rousset [START_REF] Rousset | Genetic differentiation and estimation of gene flow from Fstatistics under isolation by distance[END_REF] and Mantel tests [START_REF] Mantel | The detection of disease clustering and a generalized regression approach[END_REF] were performed using Genepop v. 3.4 [START_REF] Raymond | GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism[END_REF] across 100 000 permutations.

Cytonuclear disequilibrium analysis

Departures from random cytonuclear associations were tested using the CNDm programme [START_REF] Basten | The exact test for cytonuclear disequilibria[END_REF]. The analyses were carried out by encoding mitochondrial haplotypes as two synthetic lineages (L1 and L2). Normalised cytonuclear disequilibria (CND) were calculated following Asmussen & Basten [START_REF] Asmussen | Sampling Theory for Cytonuclear Disequilibria[END_REF] for allelic associations, and significance levels were tested using Fisher's exact test.

Results

mtDNA genetic variation

The amplified COI fragment contained a homopolymer poly-C. In most of the reactions, the sequence became mixed after the poly-C, most probably because of polymerase stutter. Sequences were thus truncated (fragment length ,200 bp). Consensus sequences shorter than 374 bp were excluded of the analysis. Over the whole mtDNA data set (134 sequences), we detected 15 haplotypes defined by 34 parsimony informative sites (9%) among 44 variable sites (12%). Within populations, haplotype diversity (h) ranged from 0.135 to 0.742 and nucleotide diversity (p) ranged from 0.00864 to 0.05973 (Table 1). Populations displaying low haplotypic richness were I03, I07, I19 and I20 (r (11) #1) while I25, I27 and IR showed relatively high values of haplotypic richness (2.5 to 3). None of the Fu's Fs and R2 values were significant.

The haplotype distribution at the population level is shown in Figure 2. This haplotype network illustrates the relationships between the 15 haplotypes and shows a clear separation of A. icterica haplotypes into two divergent lineages L1 and L2 (Figs. 1 and2). Both lineages showed a high percentage of divergence (8.7%). L1 consisted of 47 individuals and 8 haplotypes, two of which were abundant (H1 and H2). Within L1, haplotypes were Table 1. Genetic diversity in A. icterica populations. more divergent (i.e. separated by several mutational steps) than within L2. L2 included sequences from 87 individuals and the most abundant haplotype (H5) was found in all populations, representing over half the L2 individuals (60.5%). The remaining seven L2 haplotypes were relatively infrequent but all closely related. They differed from the most common haplotype H5 by at most only two mutation steps. W st analysis showed significant genetic structure at the level of the whole study (7 populations, Wst = 0.324, p,0.001). In terms of regional mitochondrial structure, the SAMOVA showed that the highest significant value (W ct = 0.409) was obtained when the populations were split into two groups (Fig. 1): group A corresponded to the I20 and I25 populations and group B was composed of the remaining populations in the Yvetot area (I03, I07, I27, I19) and the IR population (not shown in Fig. 1).

nDNA genetic variation

Among the seven microsatellite loci, the number of alleles per locus ranged from one to seven (Table 2). None of the loci showed significant linkage disequilibrium. Genetic diversity indices varied among populations (Table 1), with the I03 and I20 populations showing lower values (A = 2.54, H e = 0.339 and A = 2.48, H e = 0.395, respectively) than the other populations (3.14,A, 3.84 and 0.497,H e ,0.515). Depending on the population, the Hardy-Weinberg expectations (HWE) test showed a significant deviation for some of the loci (Ai68, Ai56, C4, PB10D and 2PE70). However, null alleles were suspected for several loci (Table 2). The estimated frequencies of null alleles ranged from 12.1% (locus 2PE70, population I19) to 40.9% (locus C4, population I19). The data set was thus corrected for null alleles and both data sets (original and corrected) were used for analyses based on allelic frequencies. When testing for mutation-drift equilibrium, a significant gene diversity excess was detected only in the I27 population using the original dataset, but also in the I07, I20, I25 and IR populations using the corrected dataset without null alleles, suggesting that these populations are good candidates for recent demographic disequilibrium arising from a population bottleneck.

Significant genetic structure was revealed at the level of the whole study (G-test, p,0.001). No pattern of isolation by distance was observed among the six Yvetot populations (p = 0.342 and p = 402 using the original and the corrected data set respectively). Clustering analysis (Figs. 3A and3B) clearly indicated genetic similarities among the I03, I27 and IR populations (Cluster 1) and among I07, I20, I25 populations (Cluster 2). The case of I19 was more ambiguous. The highest DK value was obtained for K = 2 (DK = 60.90, Fig. 3A), although the DK value for K = 3 was comparable (DK = 47.87, Fig. 3B). For K = 2, the I19 population was grouped with Cluster 1, whereas for K = 3, it formed a third group. The results of the PCA on allelic frequencies were in agreement with the results of the clustering analysis (Fig. 3C). The populations were separated into two major groups along the second axis of the PCA, which was highly significant (original dataset, p = 0.002 and corrected data set, p = 0.001). Population I19 was genetically isolated from the two clusters. Within each group, I03 and I20 populations were highly differentiated from the other two populations in their respective group (Fig. 3B).

Contemporary gene flow was detected at this regional scale with 64 individuals (29%) identified as first-generation dispersers. Among these 64 individuals, 30 belonged to Cluster 1, 25 belonged to Cluster 2 and the 9 remaining migrants belonged to the I19 population (Table 3). Most of the first-generation dispersers were assigned to populations belonging to the same cluster (19 first-generation dispersers from Cluster 1 out of 30 and 18 first-generation dispersers from Cluster 2 out of 25). Populations I03, I20 and I19 showed the lowest number of dispersers (5, 6, and 3 migrants, respectively).

Relationships between mitochondrial lineages and microsatellite alleles

The test for overall non-random association between microsatellite alleles and mitochondrial lineages revealed significant cytonuclear disequilibrium after Bonferroni correction for four of the most polymorphic microsatellite loci (Table 4). Three alleles showed significant association with L1, whereas 5 alleles were significantly associated with L2.

The test was also carried out within each population that showed both lineages, with the rarest representing at least 30% (IR, I20 and I27). The association between allele 129 at the Ai56 locus and L1 was suggested in the IR population, but was not significant after Bonferroni correction (p = 0.037). In addition, there was a trend for an association between allele 178 of the PB10 locus and L2 (p = 0.092) in the I20 population.

Discussion

Low genetic diversity within A. icterica populations

The level of polymorphism detected in A. icterica populations using microsatellite markers and COI sequences was surprisingly low (Table 5 andTable 6). The seven microsatellite markers showed low genetic variability with only 3 to 11 alleles over all loci. This polymorphism was lower than that reported in all other microsatellite datasets on earthworm populations (Table 5). For instance, the mean number of alleles per locus (N A ) ranged from 2.57 to 4.14, but values of 5 to 17 alleles have been reported in other earthworm species (Eisenia fetida [START_REF] Somers | Novel microsatellite loci for the compost earthworm Eisenia fetida: A genetic comparison of three North American vermiculture stocks[END_REF] and Lumbricus terrestris [START_REF] Velavan | High genetic diversity and heterogeneous parasite load in the earthworm Lumbricus terrestris on a German meadow[END_REF], respectively). Similarly, the sequenced fragment of the COI gene (374 bp) displayed low genetic variability in comparison to other earthworm species, despite the relatively restricted geographical scale and the short length of sequenced fragment tested in this study. For instance, only 12% of sites were polymorphic, but 33% (Hormogaster elisae [START_REF] Novo | High genetic divergence of hormogastrid earthworms (Annelida, Oligochaeta) in the central Iberian Peninsula: evolutionary and demographic implications[END_REF]) to 36.5% (Metaphire sieboldi [START_REF] Minamiya | A phylogeographic study of the Japanese earthworm, Metaphire sieboldi (Horst, 1883) (Oligochaeta: Megascolecidae): Inferences from mitochondrial DNA sequences[END_REF]) of sites are polymorphic in other earthworm species.

Ancient bottleneck events due to population isolation in periglacial refugia may be partly responsible for the current low genetic variation in this earthworm species. Among contemporary events, there are two major explanations for the low level of polymorphism in A. icterica: the occurrence of recent population bottlenecks and/or recurrent inbreeding due to reproduction between relatives. High inbreeding due to preferential mating among relatives (see for instance [START_REF] Dupont | Combined effects of bottlenecks and selfing in populations of Corella eumyota, a recently introduced sea squirt in the English Channel[END_REF]) is unlikely since deviation from HWE was inconsistent across loci and populations and could be attributed to null alleles. In some A. icterica populations, inbreeding may nevertheless occur due to a decline in effective population size. Our results indeed suggested that some populations were recovering from a recent population bottleneck. Bottlenecks can occur following colonisation events because the number of initial colonists is often small and genetic drift may Table 3. Inference of gene flow between populations belonging to each group defined using Structure: the percentage of individuals assigned to each locality as estimated by GeneClass2 is presented. result in reduced genetic variation in a newly established population [START_REF] Amos | Factors affecting levels of genetic diversity in natural populations[END_REF]. However, in an outcrossing species such as A. icterica, the likelihood of a genetic bottleneck is low because even only a few immigrants can introduce a large increase in genetic variation [START_REF] Novak | Genetic bottleneck in alien plant species: influence of mating systems and introduction dynamics[END_REF]. Only successive and drastic bottlenecks could have severely affected the genetic variation of A. icterica. Agricultural practices such as crop rotation can contribute to a fragmentation of the species habitat and thereby cause successive genetic bottlenecks [START_REF] Amos | Factors affecting levels of genetic diversity in natural populations[END_REF]. In addition, geographic isolation of populations due to natural and artificial barriers to gene flow can accentuate the loss of genetic variability through genetic drift. Low levels of gene flow in fragmented habitats can even lead to extinction of local populations due to stochastic processes. Extinction-recolonisation is a classic metapopulation scenario, with periodic extinction of all individuals in a particular patch and subsequent recolonisation of this patch from surrounding areas [START_REF] Amos | Factors affecting levels of genetic diversity in natural populations[END_REF].

Relationship between genetic and geographic distances at a regional scale

A relatively high level of genetic differentiation was revealed among localities, regardless of the marker used. Interestingly, no relationship between genetic and linear geographic distances was observed at this regional scale (populations separated by less than 13 km in the Yvetot area), corroborating other earthworm population genetics studies (reviewed in [START_REF] Costa | Genetic structure of soil invertebrate populations: Collembolans, earthworms and isopods[END_REF], but see [START_REF] Novo | High genetic divergence of hormogastrid earthworms (Annelida, Oligochaeta) in the central Iberian Peninsula: evolutionary and demographic implications[END_REF]). The lack of relationship between genetic differentiation and geographic distances was confirmed by the cluster analyses. At the nuclear level, populations were clustered into two major groups (Cluster 1 = I03, I27 and IR and Cluster 2 = I07, I20, I25), within which most of the first-generation dispersers were detected. It is noteworthy that two geographically close populations (I03 and I07) did not belong to the same cluster. There are two hypotheses that can explain the lack of correlation between the genetic and geographic distances. First, stochastic events, such as environmental changes, demographic factors (i.e. chance differences among individuals in survival or fecundity) and genetic drift may be more important than active dispersal in partitioning genetic variation at this scale (i.e. 1 to 15 km). Among earthworms, which are believed to be able to actively disperse at distances ranging from 4 to 14 m year 21 (review in [START_REF] Mathieu | Habitat quality, conspecific density, and habitat pre-use affect the dispersal behaviour of two earthworm species, Aporrectodea icterica and Dendrobaena veneta, in a mesocosm experiment[END_REF][START_REF] Eijsackers | Earthworms as colonizers of natural and cultivated soil environments[END_REF]), A. icterica is considered to be relatively vagile, being able to travel up to 500 m year 21 under conditions that trigger dispersal [START_REF] Mathieu | Habitat quality, conspecific density, and habitat pre-use affect the dispersal behaviour of two earthworm species, Aporrectodea icterica and Dendrobaena veneta, in a mesocosm experiment[END_REF][START_REF] Caro | Are dispersal behaviours of earthworms related to their functional group?[END_REF]. Tracing active dispersal events requires a study at a finer scale (,500 m 2 ).

Second, passive dispersal due to rain, floods, streams, birds, cattle or various human activities [START_REF] Eijsackers | Earthworms as colonizers of natural and cultivated soil environments[END_REF] may promote gene flow between geographically distant populations. In agricultural regions, such as in the Yvetot area, earthworms or cocoons are likely to be passively dispersed via various human activities that involve transporting soil or plant material, for instance (see [START_REF] Eijsackers | Earthworms as colonizers of natural and cultivated soil environments[END_REF] for review).

Discordant patterns of mitochondrial and microsatellite genetic structure

Two divergent (8.7%) mitochondrial lineages were observed within studied populations of A. icterica. In the Yvetot area, most of the individuals from the I20 and I25 populations belonged to L1 whereas the majority of samples from I03, I07, I19 and I27 belonged to L2 (Fig. 1). In Rouen, the population was predominantly composed of individuals belonging to L2. Genetic differentiation was confirmed in the SAMOVA analysis, with a grouping along the same lines.

Divergent sympatric mitochondrial lineages often reveal the existence of cryptic species, particularly when divergence is confirmed in the nuclear compartment of the genome and/or when reproductive isolation between lineages has been demonstrated [START_REF] King | Opening a can of worms: unprecedented sympatric cryptic diversity within British lumbricid earthworms[END_REF][START_REF] Dupont | New insight into the genetic structure of the Allolobophora chlorotica aggregate in Europe using microsatellite and mitochondrial data[END_REF]. In A. icterica, our results indicate that the two divergent lineages were randomly interbreeding with respect to mtDNA haplotypes over a relatively restricted geographical area. Deep mtDNA divergence despite clear interbreeding can reflect long periods of geographical isolation followed by secondary contact favouring gene flow, homogenising the nuclear genome over time.

Glaciation, which became increasingly severe throughout the Pleistocene, is known to have drastically modified species distributions [START_REF] Hewitt | Post-glacial re-colonization of European biota[END_REF]. Most organisms presently distributed across Europe retreated to refugia during glaciation ca. 18 000 years BP, mostly in the peninsulas of Iberia, southern France, Italy and the Balkans, and, in some cases, near the Caucasus and the Caspian Sea [START_REF] Hewitt | Post-glacial re-colonization of European biota[END_REF]. Although no common phylogeographic histories across Europe have been proposed, Taberlet et al. [START_REF] Taberlet | Comparative phylogeography and postglacial colonization routes in Europe[END_REF] highlighted some concordance in two postglacial colonisation routes: [START_REF] Lavelle | Soil Ecology[END_REF] refugia in southern France. However, there has been also recolonisation from eastern Europe and north-eastern France, and two micro-refugia were probably conserved in very specific locations in the Vosges (north-eastern France) and in Brittany (north-west of France). However, these recolonisation sources are difficult to assess [START_REF] Mathieu | Glaciation as an historical filter of below-ground biodiversity[END_REF]. Nevertheless, we can assume that the divergent A. icterica lineages originated from distinct refugia and that they merged during post-glaciation recolonisation. MtDNA divergence may thus be the result of neutral differences within the species, representing a historical mark of divergent lineages that have remerged [START_REF] Kvie | Deep sympatric mtDNA divergence in the autumnal moth (Epirrita autumnata)[END_REF][START_REF] Webb | Random interbreeding between cryptic lineages of the Common Raven: evidence for speciation in reverse[END_REF]. Over time, haplotypes are lost due to genetic drift (i.e. lineage sorting), resulting in populations monophyletic for a single gene lineage [START_REF] Harrison | Molecular-Changes at Speciation[END_REF]. Deep mtDNA divergence can only be maintained in a panmictic population with large effective population size, which effectively slows lineage sorting [START_REF] Webb | Random interbreeding between cryptic lineages of the Common Raven: evidence for speciation in reverse[END_REF]. However, we argue that A. icterica has small effective population size and has undergone serial population bottlenecks during the process of post-Pleistocene recolonisation in northern Europe, further accentuated by recent bottlenecks due to habitat fragmentation. We therefore assume that A. icterica lineages have come into contact too recently for lineage sorting to be completed. Furthermore, human activities could be, at least in part, responsible for the merging of two divergent lineages. For instance, the contribution of historical human activities to the current pattern of spatial genetic structure was documented for numerous plant species (e.g. [START_REF] Petit | Identification of refugia and post-glacial colonisation routes of European white oaks based on chloroplast DNA and fossil pollen evidence[END_REF][START_REF] Thomas | Present Spatial Diversity Patterns of Theobroma cacao L. in the Neotropics Reflect Genetic Differentiation in Pleistocene Refugia Followed by Human-Influenced Dispersal[END_REF][START_REF] Tomimatsu | Phylogeography of Camassia quamash in western North America: postglacial colonization and transport by indigenous peoples[END_REF]). Overall, our results suggest both past and ongoing anthropogenic impacts on A. icterica population genetic structure.

Here, we investigated whether the process of remerging can be traced back by studying the cytonuclear disequilibrium within contemporary populations of A. icterica. In the global dataset, some nuclear alleles were non-randomly associated with one of the two mitochondrial lineages. Because of the low number of populations displaying enough copies of both lineages and because of the relatively low number of individuals for which both COI haplotypes and microsatellites genotypes were scored, this nonrandom association could not be confirmed at the population level. Thus we cannot completely exclude the possibility that the observed cytonuclear disequilibrium is due to genetic structuring at the scale of the study.

Conclusions

Overall, this study confirmed general patterns of genetic variation observed in earthworms, such as (i) the importance of historical events for explaining their current genetic variation and (ii) the weak relationship between genetic and geographic distances suggesting the importance of passive dispersal in structuring earthworm populations. Nevertheless, some uncertainties persist such as the underlying cause of the mito-nuclear discordance and the respective roles of active versus passive dispersal in partitioning population genetic structure. In particular, it is critical to investigate how individual dispersal interacts with landscape structure. Further study is now needed to examine the extent to which barriers to movement and corridors that facilitate dispersal determine earthworm population connectivity in heterogeneous landscapes.

Figure 1 .

 1 Figure 1. Geographical distribution of Yvetot populations and repartition of cytochrome oxidase subunit I gene (COI) lineages. Lineage 1 (L1) is shown in yellow and Lineage 2 (L2) is shown in purple. Groups revealed by the SAMOVA analysis of mitochondrial data and STRUCTURE analysis of microsatellite data are shown. Land use is also indicated. doi:10.1371/journal.pone.0101597.g001

  Sample locations and codes are given, together with mitochondrial data: sample size (N mt ), number of haplotypes (N h ), haplotypic richness after rarefaction to a population size of 11 (r[START_REF] Mathieu | Habitat quality, conspecific density, and habitat pre-use affect the dispersal behaviour of two earthworm species, Aporrectodea icterica and Dendrobaena veneta, in a mesocosm experiment[END_REF] ), haplotype diversity (h) and nucleotide diversity (p), and microsatellite data: sample size (N ms ), mean number of alleles (N all ) allelic richness after rarefaction to a population size of 22 (A r ), expected heterozygosity (H e ), estimator of the inbreeding coefficient F is (significant values are in bold). doi:10.1371/journal.pone.0101597.t001

Figure 2 .

 2 Figure 2. Cytochrome oxidase subunit I gene (COI) haplotype 95% statistical parsimony network for Yvetot and Rouen A. icterica samples. H1 and H5 represent presumed ancestral sequence. Circle size is relative to the proportion of each haplotype in the sample. Mutational steps are indicated by small black circles. doi:10.1371/journal.pone.0101597.g002

Figure 3 .

 3 Figure 3. Genetic clustering of A. icterica populations based on analysis of microsatellite data. A and B. STRUCTURE Q plots representing the number of genetic nuclear groups for K = 2 and K = 3 respectively in A. icterica for I03, I07, I20, I25, I27 and IR populations. Each individual is represented by a vertical bar showing degree of admixture. C. Principal components analysis (PCA) of microsatellites allele frequencies for the whole dataset. Level of significance was derived from 1000 permutations and significant P-value is shown. doi:10.1371/journal.pone.0101597.g003

  Source localities are given in rows, recipient localities in columns. doi:10.1371/journal.pone.0101597.t003

Table 2 .

 2 Hardy-Weinberg equilibrium P-value (significant values after sequential Bonferroni correction are in bold) together with the estimation of null allele frequency (in parentheses) for A. icterica microsatellite markers in each study populations.

		C4 (11)	0.093 (0.155)	0.000 (0.385)	0.000 (0.409)	0.000 (0.252)	0.000 (0.314)	0.000 (0.346)	0.000 (0.306)
		2PE70 (3)	0.442	0.001	0.013 (0.121)	0.001 (0.274)	0.002 (0.257)	0.163	0.002 (0.284)
		2PE40 (3)	ML	ML	0.144	ML	ML	ML	1.000
		PB10D (4)	0.077	0.001 (0.256)	0.330	0.935	0.110	0.000 (0.210)	1.000
		Ai68 (8)	1.000	0.067 (0.146)	0.002 (0.216)	0.755	0.463	0.000 (0.182)	0.002 (0.230)
	Locus ID (Number of alleles)	Ai45 (4) Ai56 (7)	0.296 1.000	0.500 0.609	0.116 0.459	0.105 ML	0.031 0.405	0.613 0.217	0.315 0.000 (0.318)
		Population	I03	I07	I19	I20	I25	I27	IR

ML: monomorphic locus. doi:10.1371/journal.pone.0101597.t002

Table 4 .

 4 from Iberia and southern France towards Scandinavia and (2) from a Balkan refugium towards south-eastern France. Recent analyses of earthworm communities have shown that past climate changes have left a deep footprint on present-day earthworm diversity patterns, from community to macroecological scales[START_REF] Mathieu | Glaciation as an historical filter of below-ground biodiversity[END_REF]. It appears that earthworms recolonised France from two large Cytonuclear linkage disequilibrium between A. icterica mitochondrial lineages and microsatellite alleles, estimated using CNDm software[START_REF] Basten | The exact test for cytonuclear disequilibria[END_REF]).

	Microsatellite locus (Number of alleles*)	mtDNA lineage
		L1	L2
	Ai45 (3)	-	-
	Ai56 (5)	129	133
	Ai68 (8)	-	118
	PB10D (4)	220	178, 182
	2PE40 (3)	-	-
	2PE70 (3)	-	-
	C4 (9)	176	178
	Alleles significantly associated with mitochondrial lineage, after Bonferroni correction, are indicated.
	* the analysis was only executed for samples for which both COI haplotypes and multilocus microsatellite genotypes were scored.
	doi:10.1371/journal.pone.0101597.t004

Table 5 .

 5 Polymorphism of microsatellite loci in earthworm species.

											e ) are indicated.
	N loci N A H e Reference	8 6.63-9.63 0.725-0.774 [63]	7 2.57-4.14 0.339-0.550 This study	11 7.18 0.654 [73]	16 5.00-5.75 0.630-0.660 [57]	10 12.5 0.821 [74]	4 7.32 0.890 [75]	8 9.75 0.669 [76]	10 12.8 0.853 [77]	3 17 0.817 [58]	loci ), mean number of alleles per locus (N A ), and multilocus gene diversity (H
	Morphospecies Sampling design Geographical range N ind	Allolobophora chlorotica 2 populations in 2 countries ,500 km 62	Aporrectodea icterica 7 populations in northern France ,100 km 2 218	Aporrectodea longa 1 site ,0.50 km 2 31	Eisenia fetida 3 vermiculture stocks NA 70	Hormogaster elisae 1 site ,0.50 km 2 26	1 site 0.064 km 2 75	Lumbricus rubellus 1 site ,0.50 km 2 34	Lumbricus terrestris 1 site ,0.50 km 2 32	1 site ,0.50 km 2 281	For each earthworm species for which microsatellites have been used for population genetics study, the number of polymorphic loci (N	NA: information not provided in the study.	doi:10.1371/journal.pone.0101597.t005

Table 6 .

 6 COI sequence polymorphism in earthworm morphospecies.

										S ) and
	Number of	Morphospecies locations Sampling area Geographical range N ind L seq N h P S (%) N L Reference	Allolobophora chlorotica 38 5 European countries ,430 000 km 2 153 582 54 NA 7 [63]	Aporrectodea icterica 7 Northern France ,100 km 2 134 374 15 12 2 This study	Aporrectodea trapezoides 47 11 countries ,5 510 000 km 2 178 456 37 34 2 [78]	Dendrobaena octaedra 6 Southern Finland ,52 000 km 2 118 441 24 NA 1 [79]	Hormogaster elisae 7 Central Iberian Peninsula ,100 km 2 82 658 38 33 6 [30]	Metaphire sieboldi 64 Southern Japan ,300 000 km 2 71 690 NA 36.5 1 [59]	Rhinodrilus alatus 21 Southeastern Brazil ,16 000 km 2 69 593 59 34 6 [80]	h ), proportion of variable sites (P seq ), the number of haplotypes (N For each earthworm morphospecies for which COI has been recently used (since 2009) in a population genetics study, the length of sequence alignment (L	the number of divergent lineages (N L ) are indicated.	NA: information not provided in the study.	doi:10.1371/journal.pone.0101597.t006

PLOS ONE | www.plosone.org

July 2014 | Volume 9 | Issue 7 | e101597

Acknowledgments

The authors are grateful to V. Torres-Ramirez, S. Barot, F. Dubs, N. Bottinelli, G. Caro, A. Livet and C. Lecarpentier, for their help in sampling. We thank Naoise Nunan and Carolyn Engel-Gautier for English language editing as well as two anonymous reviewers for their helpful comments.

This work was funded by the French National Research Agency (ANR) as a part of the project EDISP No ANR-08-JCJC-0023. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Author Contributions