Genetic diversity and population differentiation in the cockle Cerastoderma edule estimated by microsatellite markers

Abstract

The edible cockle Cerastoderma edule is a marine bivalve commercially fished in several European countries that have lately suffered a significant decrease in production. Despite its commercial importance, genetic studies in this species are scarce. In this work, genetic diversity and population differentiation of C. edule has been assessed using 11 microsatellite markers in eight locations from the European Atlantic coast. All localities showed similar observed and expected heterozygosity values, but displayed differences in allelic richness, with lowest values obtained for localities situated farther north. Global Fst value revealed the existence of significant genetic structure; all but one locality from the Iberian Peninsula were genetically homogeneous, while more remote localities from France, The Netherlands, and Scotland were significantly different from all other localities. A combined effect of isolation by distance and the existence of barriers that limit gene flow may explain the differentiation observed.

Introduction

Many bivalve fisheries are currently being overexploited due to the increasing demand for seafood on a global scale. To avoid compromising the survival of these fisheries, management strategies comprising a better understanding of bivalve biology are necessary. In particular, the characterization of genetic diversity and a good understanding of population structure are important for delineating management units and maintaining the evolutionary potential. Even though low or no genetic structure is expected for most marine organisms, given their high dispersal potential and the absence of evident barriers in most of the marine environment, there are increasing reports of differentiation of bivalve populations over large and small scales (Reeb and Avise 1990; Ridgway 2001; Luttikhuizen et al. 2003; Arias et al. 2011). Exploitation without taking into account the genetic background can cause alterations of the genetic subdivision of populations and loss of genetic variation (Allendorf et al. 2008), decreasing their fitness and adaptive potential.

The cockle Cerastoderma edule (Bivalvia: Cardiidae) is a filter-feeding bivalve that inhabits the marine soft-bottom intertidal and shallow subtidal along the European coast. Its distribution comprises the eastern Atlantic coast, from Senegal in western Africa to the Barents Sea, reaching into marginal seas like the Baltic and the Mediterranean (Hayward and Ryland 1995). Cockles are relatively short lived, mature quickly, and have a high fecundity (Honkoop and van der Meer 1998) with very variable recruitment and population size. Cockle larvae are exposed to tidal currents and drift for approximately 30 days before metamorphosing and settling to the seabed as postlarvae (Dare et al. 2004). Cerastoderma edule plays a key role in coastal food webs of the northeast Atlantic and is commercially fished in several countries including United Kingdom, Spain, France, and Portugal. According to FAO statistics, European production of this species has suffered a sharp decrease, going from 107,794 tons in 1987 to 24,799 tons in 2008. Its overexploitation is especially evident in The Netherlands, being the main producing country during the 80 and 90s, with a maximum of 76,349 tons in 1989, but with no registered production nowadays.

Despite its commercial importance, few studies regarding population genetics in C. edule have been carried out. Most of them are based on allozyme markers. Beaumont et al. (1980) analyzed genetic variation at the Octopine dehydrogenase (Odh) locus in four localities coming from Wales, France, and The Netherlands, and they found significant differences between localities situated east and west of the English Channel. By contrast, Hummel et al. (1994) studied seven allozyme loci and found high gene flow among localities from Denmark, The Netherlands, and France. Beaumont and Pether (1996), by means of nine loci, concluded that there was little differentiation of cockle beds around southern UK, without any clear separation between western and eastern localities. In a recent study, Krakau et al. (2012) used mitochondrial cytochrome oxidase subunit I (COI) to assess the genetic variability in 19 sampling sites. They found two dominant central haplotypes separated by low divergence.

This study aims to assess the genetic diversity and population differentiation of eight European localities using microsatellites as molecular markers. Being more polymorphic than allozymes and mitochondrial markers, microsatellites should be more powerful to detect genetic structure and outcrossing rates.

Materials and methods

Sample collection and DNA extraction

Samples of C. edule were collected from eight locations (Fig. 1): Aveiro (Portugal), San Simón, Boiro, Corcubión and Ferrol (northwest Spain), Brest (France), Oosterschelde (The Netherlands), and Isle of Lewis (Scotland). Total genomic DNA was extracted from the adductor muscle according to Fernández-Tajes and Méndez (2007). For samples collected from Spain, DNA extraction was carried out from fresh tissue, while for the rest of localities, cockle tissue was preserved in ethanol until extraction.

Fig. 1

Sampling sites location: Aveiro (Ave); San Simón (Sim); Boiro (Boi); Corcubión (Cor); Ferrol (Fer); Brest (Bre); Oosterschelde (Oos); Isle of Lewis (ILe). Continuous line indicates a potential barrier to gene flow as identified by Barrier software

Genetic variation

Individual genotypes were determined for 11 polymorphic microsatellites according to Martínez et al. (2009). Allele frequencies, observed number of alleles per locus (N _A), observed heterozygosity (H _o), and unbiased expected heterozygosity (H _e) of Nei (1978) were obtained with the Genetix v.4.03 software (Belkhir et al. 2004). Allelic richness (Rs) per locus and per locus-locality combination was computed with Fstat v. 2.9.3.2 (Goudet 2002). A Friedman test was carried out to compare allelic richness between localities with the statistical package SPSS 16.0 (SPSS Inc.). A post hoc analysis on the Friedman test was carried out as described in Bortz et al. (2000).

Tests for agreement with Hardy–Weinberg equilibrium (HWE) and linkage disequilibrium between pairs of loci within each locality were carried out with Genepop v.4.0 software (Rousset 2008). The significance was determined by a Markov chain method using 10,000 dememorizations, 1,000 batches, and 5,000 iterations per batch. MicroChecker software (Van Oosterhout et al. 2004) was used to examine microsatellite data for evidence of null alleles and their frequency was estimated following Brookfield (1996). The frequency of private alleles in each locality was tested for correlation with sample size. Should multiple tests be carried out, P values were adjusted using the sequential Bonferroni correction (Rice 1989).

Both the Ewens–Watterson homozygosity test (Watterson 1978; Watterson 1986) and the Ewens–Watterson–Slatkin exact test (Slatkin 1994; Slatkin 1996) were carried out to check for departures from selective neutrality using Arlequin 3.11 software (Excoffier et al. 2005).

Population differentiation

Tests for genic and genotypic differentiation for all pairs of localities were performed using Genepop v.4.0 software. Weir and Cockerham’s F statistics (1984) computed over all localities, over all loci, and on a pair-wise basis between localities were obtained using Genetix v.4.03 software, which was also used to determine probability of significance of Fst values by a nonparametric permutation approach (10,000 permutations). Since the Fst value is highly dependent on the level of genetic variation, making the interpretation and comparison of the level of genetic differentiation between loci and studies difficult, a standardized measure (Fst’) was calculated in an analogous way to Hedrick’s approach (2005). That is, the observed value of Fst was divided by the maximum level of genetic differentiation that can be obtained for the observed amount of genetic variation. This maximum value was calculated by recoding the data such that all populations only contain unique alleles. The standardized measure was interpreted as equal to zero when negative values of Fst were obtained. When multiple tests were performed, significance values were adjusted using the sequential Bonferroni correction.

A hierarchical analysis of molecular variance (AMOVA) locus by locus was conducted using Arlequin 3.11 software to assess the component of genetic diversity attributable to variance among groups (different regional groupings were tested), variance among localities within groups, and variance within localities.

A Mantel test with 10,000 randomizations, as implemented in software IBDWS v.3.16 (Jensen et al. 2005), was also carried out to check for linear correlation between genetic differentiation (Fst/(1−Fst)) and the natural log of geographical distance between sample pairs (measured as the coastline distance (km) between sampling locations). Latitude and longitude coordinates were obtained for each sampling location and used along with the genetic data collected from the microsatellite analysis to generate a connectivity network of genetic distances based on Delaunay triangulation using Barrier v.2.2 software (Manni et al. 2004). Monmonier’s maximum difference algorithm was then used to identify the main putative genetic boundary across the oceanographic landscapes. Significance was established by bootstrap (100 replicates used).

Results

Genetic variation statistics by locus, locality, and overall are shown in Table 1. All loci analyzed were polymorphic in all localities; number of alleles per locus ranging from eight (CeATC2-44) to 48 (CeATC2-11). Allelic richness across localities per locus varied between 5.052 (CeATC1-36) and 19.627 (CeATC1-5). Expected heterozygosity per locus ranged from 0.384 (CeATC1-36) to 0.929 (CeATC1-5), and observed heterozygosity from 0.229 (CeATC2-12) to 0.821 (CeATC1-22). For each locus-locality pair, allelic richness ranged from 4.235 (CeATC1-36 in Aveiro) to 21.348 (CeATC1-5 in Ferrol).

Table 1 Summary statistics of genetic diversity within C. edule

Regarding localities, Friedman test (P = 0.016) detected significant differences in allelic richness, and post hoc analysis showed that Isle of Lewis’ allelic richness was significantly different from those of Ferrol and Brest (P = 0.004 and P = 0.037, respectively). Sixty-three out of the 248 alleles detected across loci were private, although their frequency was always lower than 0.05. The number of private alleles was strongly correlated with sample size (F _[1,6] = 22.992, P = 0.003, R ² _adj  = 0.793). All localities showed similar values of expected and observed heterozygosity, ranging from 0.729 (Corcubión) to 0.756 (Oosterschelde) and from 0.535 (Brest) to 0.602 (Boiro), respectively, and for all of them, values of expected heterozygosity were higher than those of observed heterozygosity.

Linkage disequilibrium tests performed for all pair of loci across localities gave no significant values (P > 0.05) for the 440 comparisons analyzed. Fifty-five out of the 88 locality-locus combinations showed no significant deviations from HWE after sequential Bonferroni correction. Four out of the 11 loci conformed to HWE in all localities (CeATC1-22, CeATC1-36, CeATC1-52, and CeATC2-34), while three (CeATC1-5, CeATC2-12, and CeATC2-51) showed significant departures from HWE in almost all localities. All combinations departing from HWE expectations but one showed positive Fis values, which indicates the existence of a heterozygote deficit. The exception comprises locus CeATC2-44 in Aveiro, which shows significant departure from HWE, although Fis value for the combination is very close to zero. This is due to the fact that even though global homozygote and heterozygote frequencies do not differ from expected values, distribution of heterozygote genotypes is not consistent with predicted values. Deviations from HWE were also detected for all localities across loci, all of which display positive Fis values ranging from 0.200 (Boiro) to 0.272 (Brest).

According to Microchecker software, eight loci showed evidence of the presence of null alleles, three (CeATC1-5, CeATC2-12, and CeATC2-51) at high frequencies (>0.2) for several localities. Allele frequencies were corrected for putative null alleles using FreeNa software (Chapuis and Estoup 2007). Nevertheless, the estimated allele frequencies and the outcome of the statistical tests of differentiation were very similar to the original ones and only these are presented. Moreover, taking into account that the presence of null alleles at frequencies above 0.2 leads to a considerable overestimation of Fst estimators and genetic distance (Chapuis and Estoup 2007), these three loci were removed from subsequent analyses.

Ewens–Watterson homozygosity test reported deviations from neutrality at locus CeATC1-5 in Aveiro, San Simón, Boiro, and Brest, while Ewens–Watterson–Slatkin exact test did so for the same locus-locality combinations and also for Oosterschelde at locus CeATC1-54 and for Isle of Lewis at loci CeATC2-34 and CeATC1-54 (data not shown).

Global multilocus Fst value was 0.015, significantly different from zero (P < 0.001), and the standardized value was Fst’ = 0.056. Estimates of Fst per locus ranged from 0.001 (CeATC2-11) to 0.059 (CeATC2-4), three of them (CeATC1-52, CeATC2-4, and CeATC2-44) being significant at the 5 % level (Table 2). These results fit in with those obtained when analyzing genic and genotypic differentiation using Genepop software, which detected significant differences (P < 0.05) for these three loci (data not shown). Of the 28 pair-wise Fst tests performed between localities, 20 were significant after Bonferroni correction (Table 3). Maximum differentiation was observed between Isle of Lewis and Ferrol (Fst = 0.049). Aveiro, San Simón, Boiro, and Corcubion localities were not significantly different from one another, while Ferrol was significantly different from all other localities except Corcubión. Brest, Oosterschelde, and Isle of Lewis were significantly differentiated from each other and from all other localities analyzed.

Table 2 Fst values per locus and overall

Table 3 Pair-wise Fst values (Weir and Cockerham 1984)

Global analysis of molecular variance (AMOVA) detected a 1.55 % variation among populations, while the percentage of variation within populations amounted to 98.44 %. Furthermore, based on the pair-wise Fst analysis results, a hierarchical analysis of molecular variance was carried out to test for the component of genetic diversity attributable to variance among different regional groups. When Aveiro, San Simón, Boiro, and Corcubión localities were grouped, AMOVA analysis identified an among groups significant component, explaining 1.73 % of the total variance. Most of the variance was distributed within localities (98.09 %) and variance component among localities within groups was nonsignificant (Table 4). When Corcubión was excluded from this group and moved to form a group with Ferrol, variance component among groups kept significant (P < 0.001), explaining 1.59 % of the total variance, but the variance component among localities within groups, although not significant, became slightly higher (data not shown).

Table 4 Analysis of molecular variance (AMOVA) performed for five groups inferred from pair-wise analysis; group 1 = Ferrol; group 2 = Corcubión, Boiro, San Simón, Aveiro; group 3 = Oosterschelde; group 4 = Isle of Lewis; group 5 = Brest

Using Fst pair-wise estimates and based on combined data from the eight loci, Barrier software identified a barrier to gene flow in the English Channel (Fig. 1). Mantel test revealed a significant correlation between Fst and geographic distance (R ² = 0.441, P = 0.005) when all eight localities were analyzed (Fig. 2).

Fig. 2

Bivariate plot of log-transformed geographic distance against Fst/(1−Fst) values

Discussion

This study uses, for the first time in C. edule, microsatellite markers to analyze genetic variation and differentiation in samples collected along the European Atlantic coast. High levels of polymorphism and heterozygosity values were detected, with allele number per locus ranging from eight to 48 and global expected and observed heterozygosity values of 0.752 and 0.572, respectively. This contrasts with values obtained for the same species using allozyme loci that showed lower number of alleles (8, Beaumont et al. 1980; 3–11, Beaumont and Pether 1996), and lower expected (0.28, Hummel et al. 1994; 0.322, Beaumont and Pether 1996) and observed (0.436, Beaumont et al. 1980; 0.24, Hummel et al. 1994; 0.278, Beaumont and Pether 1996) heterozygosity values. Thus, microsatellite markers were, as expected, more variable than allozyme markers.

All localities showed similar values of observed and expected heterozygosity, but significant differences between them were detected for allelic richness, with the lowest values obtained for localities sampled farther north; Isle of Lewis (7.958) and Oosterschelde (9.915), the latter being no significantly different from any other locality analyzed. This contrasts with what was detected by Krakau et al. (2012), who analyzed a fragment of the mitochondrial cytochrome c oxidase I gene and found higher genetic diversity values in samples of C. edule from the north of Europe than in those collected in the south. Indeed, within the southwestern group of their study, where almost all our localities are included, diversity values tend to increase as they move north. Even though our results do not agree with the existence of a northern refuge as postulated by the authors, differences might be caused by the use of nuclear markers against mitochondrial ones. Still, more sampling would be required to settle those differences.

However, our results are in line with the expectation that genetic diversity has an inverse relation with latitude due to the effect caused by Pleistocene glaciations, when marine fauna took shelter in the south (Hewitt 2004).

Compared to other bivalves, the mean number of alleles (22.545) in C. edule is in the range obtained with the same type of markers for oysters (18.533, Launey et al. 2002), clams (29.044, Vadopalas et al. 2004), and scallops (35.167, Kenchington et al. 2006). Overall expected and observed heterozygosity values obtained for C. edule were similar to those obtained for the other species of the genus (0.742 and 0.653 for C. glaucum, Tarnowska et al. 2010) but lower than those observed for other bivalve species (0.839 and 0.914, Launey et al. 2002; 0.937 and 0.596, Vadopalas et al. 2004; 0.797 and 0.731, Kenchington et al. 2006).

All but one locus-locality combinations that departed from HW expectations and all localities analyzed reported heterozygote deficits as showed by positive Fis values. Heterozygote deficits relative to HWE are common in marine bivalves and have been well documented not only for allozyme loci (Zouros and Foltz 1984; Gaffney 1990; Fairbrother and Beaumont 1993; Beaumont and Pether 1996) but for microsatellite markers as well (Launey et al. 2002; Hedgecock et al. 2004; Astanei et al. 2005; Carlsson and Reece 2007). Several biological and technical factors could account for the observed heterozygote deficit, such as inbreeding, Wahlund effect, selection and the presence of null alleles. The latter are the most likely cause for the heterozygote deficiencies detected here, as eight of the loci analyzed showed evidence for the presence of null alleles, three of them at frequencies above 0.2. Null alleles are a frequent feature in bivalves, a taxa often beset by amplification problems (Selkoe and Toonen 2006), and have been reported as the cause of heterozygosity deficit in many bivalve species such as geoduck clams (Vadopalas et al. 2004), oysters (Galindo-Sánchez et al. 2008), and mussels (Gardeström et al. 2008).

Although marine species’ high potential for dispersal, along with the continuity of the marine habitat, tends to oppose isolation and divergence of populations, results obtained in this survey support the existence of genetic differentiation among C. edule localities along the European Atlantic coast. A significant genetic structure was revealed at the level of the whole study (Fst = 0.015, P < 0.001). This Fst value is similar to those obtained in other studies that used microsatellite markers to assess genetic differentiation in other bivalve species such as oysters (Launey et al. 2002; Galindo-Sánchez et al. 2008). Morton et al. (1993) suggested that measures of differentiation for highly variable loci may be as much as an order of magnitude lower than that of traditional markers, such as allozymes. Therefore, estimates of differentiation were expected to be low because of the high variability of microsatellite loci, as the practical upper limit for Fst is actually the level of expected homozygosity, which is low for this type of markers, implying a reduction in their maximum Fst value (Charlesworth 1998; Hedrick 1999). Nevertheless, despite the effect of polymorphism on deflating Fst expectations, highly variable markers are valuable in exact tests of differentiation not only because of their sensitivity to gene flow (Ross et al. 1999) but also because their high variability gives them greater statistical power (Estoup et al. 1998). This has indeed been shown for several studies on aquatic species where microsatellites revealed a finer resolution of genetic structuring than allozymes (Ruzzante et al. 1999; Lundy et al. 2000; Lage et al. 2001; Mattiangeli et al. 2002; Knutsen et al. 2003; Taylor et al. 2003; Carlsson et al. 2004; Jorgensen et al. 2005). Considering the differences relative to Fst values depend on the marker used, a standardized value of Fst was computed (Fst’ = 0.056), indicating moderate genetic differentiation according to Wright’s (1978) guidelines. This value is higher than those obtained with allozyme loci in the same species; Hummel et al. (1994) obtained an overall Fst value of 0.04, detecting no significant differentiation between locations analyzed, while Beaumont and Pether (1996), with an average Fst of 0.015, detected little differentiation between samples coming from the UK southern coast.

Fst pair-wise values showed no differentiation between four localities situated northwest of the Iberian Peninsula (Aveiro, San Simón, Boiro, and Corcubión), which is supported by the AMOVA analysis, with a nonsignificant variance component among localities within groups that comprise 0.19 % of the total variance (Table 4). This indicates the existence of high gene flow rates between them. Despite the distance between Aveiro and the Spanish localities, which doubles that between San Simón, Boiro, and Corcubión, all four localities were genetically similar. Aveiro is a shallow lagoon system with high tidal flow that allows a large water renovation. This may promote larval drift to the coastal waters (Joaquim et al. 2010) and be responsible for larvae dispersal, preventing genetic differentiation. Although Ferrol is also located northwest of the Iberian Peninsula, it showed significant differences to all other four Iberian localities analyzed but Corcubión. This may be due to their geographic position with respect to Cape Finisterre, where a major oceanographic boundary has been described (López-Jamar et al. 1992). Another possible explanation calls on the hypothesis of “Sweepstakes Reproductive Success”, which suggests that the extremely large variance in individual reproductive success could produce, among other effects, the presence of genetic heterogeneity at reduced spatial scales (see review in Hedgecock and Pudovkin 2011).

Localities of Brest, Oosterschelde, and Isle of Lewis resulted in differences not only between them, but from all other localities analyzed. This could be a product of their geographical location, as gene flow modulates over greater distances and eventually reaches a level where the exchange of migrants is not strong enough to counter the effects of random genetic drift (Lind et al. 2007). The Mantel test revealed the existence of a significant correlation between pair-wise Fst estimates and geographical distance when all localities were analyzed. However, even though the analysis was statistically significant, graphic representation (Fig. 2) did not adjust well to a linear model; from a certain geographic distance onwards, genetic distance tends to increase sharply, more consistent with an exponential shape than a linear one. Therefore, samples at intermediate distances should be analyzed to confirm that genetic divergence between localities is a result of isolation by distance. It is important to note that geographical distance measures may not reflect real dispersal pathways, which are likely to be more complex. Isolation by distance may explain just part of the genetic differentiation observed. Nearshore circulation patterns may be affected by shoreline irregularities such as estuaries or bays, resulting in eddies and fronts that may retain larvae over time scales comparable with dispersal duration, and thus may account for enhanced local recruitment at specific locations (McShane et al. 1988; Dupont et al. 2007). This could be the case for the locality of Brest, where genetic isolation could also be influenced by the characteristics of the bay: a semi-enclosed marine ecosystem (Richard et al. 2006) favorable to the isolation of the larval pool and, in the long term, genetic isolation. Genetic isolation of this bay has already been observed for the molluscan species Crepidula fornicata (Dupont et al. 2007).

Barrier software suggested the presence of a barrier to gene flow situated in the English Channel. The precise situation of this barrier cannot be determined since no locations along the English Channel were analyzed. Other studies have found gene flow restrictions in this region, although the exact boundary location varied. Thus, Roman and Palumbi (2004) in their analysis of population structure in Carcinus maenas by means of mitochondrial cytochrome c oxidase I DNA found a significant break between western and northern Europe located between The Netherlands and Germany. In C. edule, the existence of such genetic discontinuity has also been described by Krakau et al. (2012) using a fragment of the same mitochondrial sequence. Other authors situate the genetic break further south; Jolly et al. (2005) found a phylogeographic break for the polychaete Pectinaria koreni that separates populations coming from Brittany and the English Channel. Also, Dupont et al. (2007) found a strong hierarchical structure for the molluscan species C. fornicata with two groups located on both sides of the Cotentin Peninsula. In spite of the still undefined location of the genetic barrier, it is clear that the English Channel currents play an important role in gene flow, with restrictions that may promote population differentiation.

Although microsatellite markers detected genetic differentiation among localities surveyed, the level of differentiation detected in C. edule is lower than in the other species of the genus, C. glaucum. In the latter, microsatellite markers detected a high level of genetic structure between localities along the coast of Europe (Tarnowska et al. 2010). Allozyme markers also showed a higher level of differentiation in C. glaucum with respect to C. edule (Hummel et al. 1994; Mariani et al. 2002; Nikula and Väinölä 2003). Hummel et al. (1994) pointed out that these differences could be explained by the different character of the species habitat, with C. glaucum individuals found in smaller, more or less limited water bodies, thus resulting in a more limited gene flow.

In conclusion, this work makes novel use of microsatellite markers to provide estimates of genetic diversity and population differentiation in C. edule from the Iberian Peninsula and other European locations. All localities displayed similar heterozygosity levels, but significant differences in allelic richness were observed, with the lowest values found in localities situated farther north. Localities from the Iberian Peninsula were genetically homogeneous, except for the one located farther north, and were different from those of other European sites, which may be partially explained by an isolation by distance process. Scottish, French, and Dutch localities displayed contrasting results, likely due to a combined effect of isolation by distance and the existence of barriers limiting gene flow between them. More localities should be analyzed by means of microsatellite markers in order to achieve complete knowledge of population structure and to establish the mechanisms that determine genetic differentiation in the cockle C. edule.