Null alleles are ubiquitous at microsatellite loci in the Wedge Clam (Donax trunculus)

Introduction

Recent studies have reported an unusually high frequency of non-amplifying alleles, also known as null alleles, at microsatellite loci in bivalves (Hargrove et al., 2015; Chiesa et al., 2016) including the wedge clam (Nantón et al., 2014). Null alleles occur when mutations in the binding site of the targeted DNA sequence prevent the efficient annealing of at least one primer resulting in failure of amplification during the PCR reaction. Null alleles can also occur via segmental aneuploidy, where one chromosome has a deletion containing the primer binding site. For instance, in a study of allelic BAC sequences in the Pacific oyster obtained by the Oyster Genome Consortium, 42 of 101 microsatellite loci were found to occur in the hemizygous state, owing to indels of various sizes (P. Gaffney, pers. com., 2016). Null alleles can occur in homozygous state when the samples produce no amplification at all, and in heterozygous state, when the sample appears as a homozygous individual for a particular locus. Nantón et al. (2014), characterised the first nineteen polymorphic microsatellite markers in D. trunculus to assist in the management of the species and allow the delineation of conservation units essential in fisheries. Of the 19 loci characterised in this study, at least 10 showed the presence of null alleles with frequencies ranging from 0.109 to 0.277. Another study aiming to examine the genetic variability and relationships between two Manila clam (Ruditapes philippinarum) populations in Korea also suggested the occurrence of null alleles in 7 loci out of 10 (Kim et al., 2014). More recently Chiesa et al. (2016) demonstrated the occurrence of null alleles through rigorous analytical tests in this species. The presence of null alleles has also been reported in a study that assessed the genetic diversity of four wild and six hatchery stocks of the hard clam (Mercenaria mercenaria) in Florida. Null alleles were present in all loci with frequencies ranging from 0.015 to 0.296, but for four loci out of seven, the frequency of null alleles was particularly high (>0.146) (Hargrove et al., 2015).

Simulation studies suggest that null alleles with frequencies between 5% and 8% should have only minor effects on classical estimates of population differentiation, but that higher frequencies would bias such parameters (Chapuis & Estoup, 2007). Null alleles have been associated with heterozygous deficits in many studies and several have tested their presence and effects in population structure estimates using various analytical and simulation tools (e.g., Dąbrowski et al., 2015 and references therein). However, to the best of our knowledge no study has empirically tested for its consequences in estimating differentiation using population samples separated by well characterised biogeographic barriers.

The wedge clam constitutes an important fishing resource in southern Spain and Portugal due to its high economic value (Gaspar, Ferreira & Monteiro, 1999). In Europe alone, the recorded landings in the last 12 years total 11,202 tons with a maximum yield of 1,355 tons in 2005 (FAO-FIGIS, 2016). Overall, the species has experienced a steady decline in recent years reaching only 757 tons in 2014 (FAO-FIGIS, 2016). Despite its overexploitation threat, D. trunculus has only recently been studied from a population genetics (Marie et al., 2016). Understanding connectivity among populations through genetic structure provides tools to determinate the appropriate units and spatial scale for fisheries conservation and management (Waples & Gaggiotti, 2006; Funk et al., 2012). D. trunculus is an Atlantic–Mediterranean warm-temperate species found in the Black Sea, in the Mediterranean Sea (Bayed & Guillou, 1985) and from Senegal to the northern Atlantic coast of France (Tebble, 1966). It inhabits sandy beaches exposed to tidal rhythms, characterised by intense wave action and sediment instability (Brown & McLachlan, 1990). In these environments, populations are capable of reaching sufficiently high densities to support large commercial fleets (Gaspar, Ferreira & Monteiro, 1999). Because, wedge clams are filter feeders, they play important roles in the trophic structure of beaches, but can also accumulate xenobiotic compounds making them, ideal model organisms for environmental monitoring (Saavedra & Bachere, 2006). For this reason, they may also constitute potential risks for human health when they are eaten and thus have been extensively studied from an ecotoxicology perspective (Tlili et al., 2010; Yawetz et al., 2010; Bouzas et al., 2011; Company et al., 2011; Hamdani & Soltani-Mazouni, 2011; Tlili et al., 2011).

Motivated by our interest in the conservation genetics of this important shellfish resource, we used sequences generated by one NGS platform (Roche-454) for microsatellite characterisation. Evidence of successful PCR amplification and allelic size variation across a subset of population samples allowed us to optimise 16 loci into four multiplex panels, permitting the PCR amplification of several markers simultaneously. We then genotyped two adjacent locations in the southern Atlantic Coast of Spain. One of them is situated within the boundaries of the Doñana National Park, exploited, due to historic rights, only by artisanal hand-operated dredges with restricted quotas. The second is approximately 85 km apart and is exploited by a commercial boat-operated dredging fleet. We also genotyped a sample from the south Mediterranean Coast of Spain. These samples, located on both sides of the Gibraltar strait, which is a well-known major barrier to gene flow (Perez-Losada et al., 2007; Sá-Pinto et al., 2012), allowed us to test the strength of the analyses using loci that have high frequencies of null alleles in this species. Finally, we also show through empirical and analytical tests that null alleles are ubiquitous in this species as previously suggested (Nantón et al., 2014).

Materials & Methods

Results and Discussion

A total of 18,795 reads with an average length of 232 bp were obtained from 454 pyrosequencing of a microsatellite enriched library and assembled into 595 contigs. We identified 267 contigs (44.9%) that contained microsatellite loci with at least 10 di- or seven tri-nucleotide repetitions. Of these, there were a total of 85 contigs (31.8%) we identified as PAL of di- or tri-nucleotide microsatellites and 182 contigs (68.2%) as nonPAL (i.e., contigs containing di- or trinucleotide microsatellites without suitable flanking PCR-primer sites). Comparing across the two classes of repeats in PAL sequences, di-nucleotide repeats represented 66.8% and tri-nucleotides 33.2%.

From the microsatellite containing sequences, we selected 52 loci for testing their usefulness in multiplex genotyping and 37 produced PCR amplification when visualised on agarose gels and ethidium bromide staining. Of these 37 loci, nine gave multiple bands, faint bands, or fragments clearly defined but of a different size than the one expected from the original sequence. These loci were discarded from subsequent analyses. Of the remaining 28, 21 could be genotyped and displayed polymorphism in 48 individuals. These 21 loci were deposited in GenBank under accession numbers HG792255.1–HG792275.1. Five loci, D.tru28, 30, 31, 39 and 46, amplified two or three loci simultaneously and were discarded from further analyses. For the remaining loci, the number of alleles per locus ranged from eight to 28 ($\bar{\chi }=18.25$, σ = 5.86). Expected and observed heterozygosities, and polymorphic information content ranged from 0.15 to 0.82 ($\bar{\chi }=0.53$, σ = 0.21), 0.54 to 0.94 ($\bar{\chi }=0.81$, σ = 0.11), and 0.51 to 0.93 ($\bar{\chi }=0.79$, σ = 0.12), respectively. F_IS estimates ranged from −0.07 to 0.74 ($\bar{\chi }=0.35$, σ = 0.04) (Table 1). Private alleles were present in all loci for the south Atlantic Coast of Spain and in 11 loci for the South Mediterranean Coast of Spain (Table S1 ). There was no evidence of large allele dropout or genotyping errors due to stutter peaks for the 16 loci. However, all loci but two (D.tru4 and D.tru16), showed a significant excess of homozygotes (deviation from Hardy–Weinberg equilibrium, P < 0.001). There was no evidence of linkage disequilibrium between any pair of loci. Genotyping was consistent across the 48 samples, which were each tested twice yielding the same genotypes. The multiplex optimisation yielded a total of four panels for the 16 primer pairs. The numbers of possible multiplex combinations was limited due to the allele size ranges of the markers and differences in the annealing temperatures of the primers. We were thus able to design two 3-plex, and two 5-plex multiplex reactions after excluding the duplicated loci.

Before the advent of next-generation sequencing technologies, no genetic marker has found such widespread use as microsatellites in the last two decades. This study joins an ever increasing number of others that have successfully employed NGS to detect a virtually unlimited number of polymorphic microsatellites across a wide variety of taxa (e.g., (Wang et al., 2012b; Wasimuddin et al., 2012; Whitney & Karl, 2012; Miller et al., 2013). Our study takes advantage of Roche-454 titanium chemistry. Microsatellite-containing sequences from shotgun and enriched libraries have also been reported across a broad range of taxa in recent years (Angeloni et al., 2011; Santure et al., 2011; Wang et al., 2012a; Wang et al., 2012b; Zhang et al., 2012). Together, these studies demonstrate that the main obstacle of costly de novo isolation for using microsatellite markers for population studies and characterisation has been rapidly overcome by the advent of NGS technologies. More importantly, however, is the fact that a growing proportion of research projects, working on modest financial resources, are now capable of developing a suitable number of microsatellite markers at minimal laboratory costs for virtually any taxa (Schoebel et al., 2013). The number of microsatellite loci detected in this study illustrates how abundant these sequences are in eukaryotic genomes. Comparing the number of loci detected across studies and taxa is inherently difficult because genomes vary substantially in their frequency of microsatellites and genome size (Schoebel et al., 2013). However, the number of loci detected here favourably compares with other studies (e.g., Wang et al., 2012a; Wang et al., 2012b; Zhang et al., 2012). For example, in a recent study that used this technology, Schoebel et al. (2013) carried an extensive survey of microsatellite abundance in 17 non-model species including plants, fungi, invertebrates, birds and a mammals and determined if flanking regions were suitable for primer development. They concluded that that depending on the species, a different amount of 454 pyrosequencing data might be required for successful identification of a sufficient number of microsatellite markers for ecological genetic studies. Irrespective of the difficulty to compare microsatellite abundance in different taxa, an important message of our study is that a virtually unlimited number of microsatellite loci can be identified from the large amount of sequence data generated with NGS technologies. It is worth noticing as well that 454 pyrosequencing has been replaced by other technologies which generally yield many more sequences, but with shorter reads. The implications of these shorter reads may have an important impact on marker development as microsatellites require the availability of long flanking sequences for primer design. However, we have to consider that, despite the many advantages of microsatellite markers, admittedly, the development of SNP markers by NGS or SNP genotyping by sequencing have increasingly become attractive alternatives for many reasons such as availability of high numbers of annotated markers, low-scoring error rates, improved genotyping results for poor quality samples, highly informative in terms of segregation among populations, and the ability to examine neutral variation and regions under selection, among others (Brumfield et al., 2003; Rosenberg et al., 2003; Morin et al,, 2004; Liu et al., 2005; Grewe et al., 2015).

Once the multiplex panels were optimised, we genotyped a total of 195 individuals from the three sampled locations. Table 1 summarises the genetic diversity found in each locus for all the samples analysed. The initial uncorrected typing error rate was 0.0023 per reaction or 0.0015 per allele. The frequency of single-locus genotypes missing in the data set was <0.05. Adjacent allele scoring error analysis of the final corrected data set yielded no significant deviations (X² < 2, P > 0.15), hinting that there were very few remaining scoring errors. From individuals genotyped as controls in each plate, we calculated a mistyping rate of 0.0019 per allele genotyped. Type A and type B error rates were 0.015 and 0.021 for multiplex 1, 0.013 and 0.008 for multiplex 2, 0.026 and 0.051 for multiplex 3, and 0.009 and 0.014 for multiplex 4, and 0.0152 and 0.023 across all loci and were mainly caused by excessive stuttering of some long alleles. Once again, all loci but D.tru4 and D.tru16 showed a significant excess of homozygotes (deviation from Hardy–Weinberg equilibrium, P < 0.001) (Table 1). Exact tests for genotypic linkage disequilibrium confirmed the absence of physical linkage among most loci (P < 0.05 after Bonferroni correction).

The results of the re-amplification of a subsample of individuals under lower stringency mostly confirmed the above findings. However, only 10 loci produced codominant fragments (i.e., one allele of maternal and one of paternal origin) under relaxed annealing temperatures (Table S1 : Original Annealing Temperature, New Annealing Temperature). Locus D.tru16, which initially had not deviated from HWE, amplified multiple loci and thus could not be used for this comparison. Furthermore, two additional loci showed null allele frequencies not different from zero (D.tru6 and D.tru23). This result suggests that non-amplifying alleles can be detected by reducing the stringency of the annealing temperature. For the remaining eight loci, no differences were observed in the heterozygous deficit or the null allele frequencies supporting the presence of null alleles in the wedge clam genome (Table 2).

The fact that most loci had large heterozygous deficits was an intriguing result which deserves a thorough discussion. The first consideration that needs to be evaluated in addition of the occurrence of null alleles, is the presence of large allele dropout. Another consideration is the level of stuttering as it may hamper the reliable scoring of alleles differing within one or two repeated units in some loci. Another possible explanation could be the differences in allele size as reported by previous studies that have found that loci with longer alleles tend to have higher dropout rates than those with shorter alleles (Sefc et al., 2003; Buchan et al., 2005; Broquet, Ménard & Petit, 2007). In our study, allele size ranges ranged from 60 to 24 ($\bar{\chi }=46.2$, σ = 12.12) and from 110 to 30 ($\bar{\chi }=54.7$, σ = 26) for the original and relaxed annealing temperatures respectively (Table 2). It is thus expected that the largest allele at a locus would amplify with less efficiency than the shortest one as the latter would outcompete the former during amplifications (Sefc et al., 2003; Buchan et al., 2005; Broquet, Ménard & Petit, 2007). This will artificially increase the frequency of null alleles even in the absence of true null alleles.

The analysis of null allele frequencies per population using the maximum likelihood estimates implemented in ML-NullFreq confirmed the presence of null alleles in all loci but D.tru4 and D.tru16. For these loci, null-allele frequencies estimates ranged between 0% and 31.2% ($\bar{\chi }=0.160$, σ = 0.098) and all loci but D.tru4, D.tru16 and D.tru23 had a null-allele frequencies above 5% (Table 3). Results from FreeNa showed lower null allele frequencies for each locus (except for the locus D.tru19) than results obtained using ML-NullFreq. Null allele frequencies estimates ranged between 0% and 29.3% ($\bar{\chi }=0.149$, σ = 0.090; Table 3). This high frequency of null alleles found in this study is consistent with many reports in other bivalve species (Launey et al., 2002; Nantón et al., 2014; Chiesa et al., 2016). Clams are filter feeders and thus, are prone to accumulate xenobiotic compounds (Tlili et al., 2010; Yawetz et al., 2010; Bouzas et al., 2011; Company et al., 2011; Hamdani & Soltani-Mazouni, 2011; Tlili et al., 2011). It is plausible that the mutagenic action of some pollutants such as heavy metals (e.g., mercury or cadmium) (Wong, 1988) can be responsible for an increased mutation rate in clams. This could be further investigated by comparing samples representing genetically distinct populations where organisms are either heavily exposed to pollutants in beaches adjacent to major urban areas or those occupying habitats in remote locations where the level of pollution is likely to be low and far from industrial pollution.

As expected due to the absence of biogeographic barriers and the proximity of the sampling sites in the south Atlantic coast of Spain, the fixation index between samples from Isla Canela and Doñana National Park was not significantly different from zero (Arlequin: F_ST = 0.0028–P > 0.05; FreeNa: F_ST including ENA correction =0.0025) indicating uninterrupted migration of larvae across the area. Furthermore, as predicted, the fixation index between these two localities and the Mediterranean sample was significantly different from zero (Arlequin: F_ST = 0.033–P < 0.001; FreeNa: F_ST including ENA correction =0.032). The STRUCTURE analysis found that the two populations from the south Atlantic coast of Spain were genetically homogeneous while they were significantly different from the Mediterranean population with the highest ΔK value equal to the number of biogeographic regions (K = 2) and the programme showed that the 10 independent runs from K = 1 to K = 3 produced consistent results whether or not the RECESSIVEALLELES option was used in the analysis (Fig. 1). Calculation of the statistic ΔK from the STRUCTURE runs indicated that the two biogeographic regions analysed have clearly differentiated populations.

When running STRUCTURE with the complete data set without a different locus at a time for each run, our results also showed that three microsatellite loci were very informative (D.tru4, D.tru29 and D.tru32; Fig. S1, Table S2 ). Indeed, compared to the reference average level of individuals’ assignment (84.7% ± 2.5 and 90.6% ± 1.7 for the South Atlantic Coast of Spain and the South Mediterranean Coast of Spain, respectively), without these three loci, the average level of individuals’ assignment dropped to 65.5% ± 6.5, 68.2% ± 3.7 and 77.9% ± 2.6 in the South Atlantic Coast of Spain and to 38.2% ± 5.4, 89.5% ± 2.2 and 89.9% ± 2.8 in South Mediterranean Coast of Spain, respectively, when excluding either locus D.tru4, D.tru29 or D.tru32. Surprisingly, two of these loci (D.tru29 and D.tru32) have large frequencies of null alleles (0.311 ± 0.030 and 0.228 ± 0.027 for the locus D.tru29 and D.tru32 respectively; Table 3) confirming that even when microsatellite loci display a large frequencies of null alleles, they can, nevertheless, be useful to identify the genetic structure of populations. This conclusion was also corroborated with the F_ST results using FreeNa and Arlequin, which showed that the use of the null allele’s correction option had no effect in the F_ST estimates.

To the best of our knowledge, it is unusual that 87.5% of newly developed loci are affected by any form of genotyping errors or allele dropout during PCR due to the presence of long alleles. Although, simulation studies suggest that null alleles with frequencies between 5% and 8% should have only minor effects on classical estimates of population differentiation (Chapuis & Estoup, 2007), in our case null allele frequencies will render all these loci, but three, completely useless in population genetic analyses if the heterozygous deficits were due to null alleles. The fact that we found panmixia between samples from adjacent localities with no apparent barriers to gene flow and a well-defined structure for populations separated by a well-characterised biogeographic barrier suggests that the large heterozygous deficits we found, cannot be entirely due to null alleles or genotyping artefacts. If this was the case, one would have expected significant, albeit artificial, genetic differences between the two adjacent localities. Alternatively, it suggests that despite the presence of null alleles, clearly differentiated populations can be detected, although levels of genetic diversity are most likely underestimated. Admittedly, however, results pertaining to genetic structure presented here should be interpreted cautiously until they are corroborated by other approaches, namely SNP genotyping.

Conclusions

In conclusion, we have shown that null alleles at microsatellite loci are unusually abundant in the wedge clam and that genotyping artefacts are unlikely to completely explain the large heterozygous deficits observed in all population samples. We have also shown that the high frequencies of null alleles observed at these loci do not appear to have a significant effect in the population genetic parameters commonly assessed for microsatellite loci. The expected population structure between Atlantic and Mediterranean samples was confirmed for the tests carried out. Furthermore, the unusually high frequency of null alleles reported in the wedge clam and in other bivalve species is worth investigating further. As mentioned above, it is reasonable to hypothesise that pollutants may increase mutation rates in bivalves, resulting in high frequencies of null allele. Considering the economic value of D. trunculus, additional investigations are essential to validate this assumption.

Supplemental Information