Next Article in Journal
The Increase of Simple Sequence Repeats during Diversification of Marchantiidae, An Early Land Plant Lineage, Leads to the First Known Expansion of Inverted Repeats in the Evolutionarily-Stable Structure of Liverwort Plastomes
Next Article in Special Issue
Within-Generation Polygenic Selection Shapes Fitness-Related Traits across Environments in Juvenile Sea Bream
Previous Article in Journal
CRISPR/Cas9-mediated Disruption of Fibroblast Growth Factor 5 in Rabbits Results in a Systemic Long Hair Phenotype by Prolonging Anagen
Previous Article in Special Issue
Seascape Genetics and the Spatial Ecology of Juvenile Green Turtles
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Genes
Volume 11
Issue 3
10.3390/genes11030298
genes-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Population Genetic Structure Is Unrelated to Shell Shape, Thickness and Organic Content in European Populations of the Soft-Shell Clam Mya Arenaria

by
Michele De Noia
1,2,3,4,†,
Luca Telesca
2,3,†,
David L. J. Vendrami
1,†,
Hatice K. Gokalp
1,
Grégory Charrier
5,
Elizabeth M. Harper
2,3 and
Joseph I. Hoffman
1,3,*
1
Department of Animal Behavior, University of Bielefeld, Postfach 100131, 33615 Bielefeld, Germany
2
Department of Earth Sciences, University of Cambridge, Downing Street, Cambridge CB2 3EQ, UK
3
British Antarctic Survey, High Cross, Madingley Road, Cambridge CB3 OET, UK
4
Institute of Biodiversity, Animal Health & Comparative Medicine, College of Medical, Veterinary & Life Sciences, University of Glasgow, Glasgow G12 8QQ, UK
5
University Brest, CNRS, IRD, Ifremer, LEMAR, F-29280 Plouzané, France
*
Author to whom correspondence should be addressed.
†
These authors contributed equally to this work.
Genes 2020, 11(3), 298; https://doi.org/10.3390/genes11030298
Submission received: 27 January 2020 / Revised: 6 March 2020 / Accepted: 6 March 2020 / Published: 11 March 2020
(This article belongs to the Special Issue Genetic Diversity of Marine Populations)
Browse Figures
Review Reports Versions Notes

Abstract

:
The soft-shell clam Mya arenaria is one of the most ancient invaders of European coasts and is present in many coastal ecosystems, yet little is known about its genetic structure in Europe. We collected 266 samples spanning a latitudinal cline from the Mediterranean to the North Sea and genotyped them at 12 microsatellite loci. In parallel, geometric morphometric analysis of shell outlines was used to test for associations between shell shape, latitude and genotype, and for a selection of shells we measured the thickness and organic content of the granular prismatic (PR), the crossed-lamellar (CL) and the complex crossed-lamellar (CCL) layers. Strong population structure was detected, with Bayesian cluster analysis identifying four groups located in the Mediterranean, Celtic Sea, along the continental coast of the North Sea and in Scotland. Multivariate analysis of shell shape uncovered a significant effect of collection site but no associations with any other variables. Shell thickness did not vary significantly with either latitude or genotype, although PR thickness and calcification were positively associated with latitude, while CCL thickness showed a negative association. Our study provides new insights into the population structure of this species and sheds light on factors influencing shell shape, thickness and microstructure.

1. Introduction

The soft-shell clam Mya arenaria Linnaeus, 1758, is a marine bivalve that occurs in numerous intertidal infaunal communities across Europe and North America. While this species has a relatively high dispersal potential during the planktonic larval stage and as juveniles [1], the contemporary geographic distribution of M. arenaria appears to have been largely influenced by human-mediated translocations. Although deliberate introductions have been documented in the literature [2], the majority of introductions were probably unintentional, occurring either as a byproduct of oyster transplants [3] or via shipping, as this species can survive for extended periods in ballast water [4,5].
Present in Europe during the Pliocene [6], the soft-shell clam is believed to have been extirpated from this continent during the Pleistocene glaciations [7]. Radiocarbon dating of shells found along the coast of the North Sea suggests that M. arenaria was subsequently reintroduced into Europe by the Vikings during the 13th to 15th centuries [8,9], making it one of the oldest marine invaders of European coasts [10]. Following recolonization from the North Sea and probable further man-mediated introductions, the soft-shell clam can nowadays be found around most of the continent, including in the Mediterranean Sea [11,12], on the Iberian Peninsula [13], along the Atlantic coast of France [14], around the British Isles [15,16] and in the North Sea [1,17,18].
Despite M. arenaria being an important species in many European coastal ecosystems, little is known about its population genetic structure, especially when compared to other European shellfish species [19,20,21,22,23]. In addition to two studies of M. arenaria from North America, both of which included a single European population [24,25], only three studies to our knowledge have focused on the population genetics of the soft-shell clam in Europe [5,26,27]. In all cases, little in the way of population genetic structure could be detected and genetic diversity appeared to be low. However, these studies were exclusively based on either allozymes or mitochondrial DNA (mtDNA), which provide limited power to detect genetic differentiation, especially when population structure is shallow.
More recently, microsatellite markers have been developed for the soft-shell clam, which, in contrast to the studies described above, have uncovered relatively strong population genetic structure and higher levels of genetic diversity in North America [25,27,28]. One of these studies [27] also generated microsatellite data for four European M. arenaria populations, although these were locally restricted to the Southern coast of Ireland, North Wales and the Netherlands. This study found clear evidence for genetic differentiation at microsatellites but not mtDNA [27], implying that population genetic structure may be present over a broader scale in Europe.
Characterizing population genetic structure not only provides information on patterns of gene flow and genetic drift, but it may also shed light on the extent to which morphological or other phenotypic traits are plastic versus under genetic control [29,30]. The underlying logic of this approach is to relate phenotypic differences among populations to the underlying genetic structure. On the one hand, the presence of pronounced phenotypic differences among genetically indistinguishable populations may be suggestive of a strong influence of the environment on the traits in question. On the other hand, coincident patterns of genetic and phenotypic variation point towards a possible role of genetics. Shell morphological traits of marine invertebrates provide an ideal test bed for investigating the relative contributions of genes and the environment towards phenotypic trait variation, as protocols have already been established for quantifying shell morphological and microstructural variation [31,32,33] and it has been suggested that phenotypic plasticity may play an important role in determining individual variation in shell features of many marine organisms including limpets [34], sea snails [35] and scallops [36].
One obvious shell morphological trait where great variation among individuals can be observed is shape. Shell shape variation has been studied in a variety of molluscan species and it appears to be significantly influenced by both environmental and genetic effects, as well as their interaction [37,38,39]. Abiotic factors associated with shape variation include water temperature and acidity [40], while an effect of latitude has also been reported [41], probably due to changes in environmental conditions linked to latitude. Moreover, biotic variables such as predation have also been demonstrated to be important in determining individual variation of shell shape [42]. However, little is known about the specific causes of shell shape variation in M. arenaria. While Swan [43] and Emerson [44] have provided evidence for a potential role of substrate typology and water flow, to our knowledge no study has investigated the potential effect of genetics on shell shape variation.
The shell of M. arenaria is composed of four main layers. These are the outermost periostracum, the granular prismatic (PR), the crossed-lamellar (CL) and the complex crossed-lamellar (CCL) layers [45] (Figure 1 and Figure S1). While the periostracum consists of sclerotized proteins, the other layers are mainly composed of aragonite (CaCO3) crystals and inter-crystalline biomineral organic matrix. Differences in the energetic cost of producing and maintaining different shell structures and components [46,47] as well as geographical variation in physical and biotic stressors that are likely to exercise a selective pressure on shell morphology, are expected to influence variation in shell microstructural composition and thickness [48,49]. The fact that M. arenaria is widespread and locally abundant, combined with the availability of recently developed microsatellites [25,28], makes this species a good candidate for investigating the relative contributions of genetic and geographical factors towards variation in shell morphological traits.
Here, we collected a total of 266 M. arenaria samples from nine locations around the coastline of Europe, spanning a latitudinal cline from Lisbon to St. Andrews (Figure 2). Genetic data for 12 microsatellite loci were generated and analyzed in combination with data on shell shape, total thickness and the thickness of the PR, CL and CCL layers. For a subset of samples originating from the extreme northern and southern sampling sites, the organic content of these layers was also measured. This allowed us to characterize the pan-European population genetic structure of this species and to test for effects of genotype and latitude on key shell characteristics.

2. Materials and Methods

2.1. Sample Collection and DNA Extraction

Between 15 and 40 M. arenaria specimens of wild origin were collected from the eulittoral zone of each of nine locations within Europe (Figure 2, Table 1). The length of the shell of every individual was measured with digital calipers (0.01 mm precision) and used as a within-population proxy for age [50]. The shells were retained and processed to determine shell shape as well as to quantify shell thickness and organic content as described subsequently. Approximately 1 cm3 of tissue, either adductor muscle or mantle, was taken from each individual and stored in 95% ethanol at −20 °C for genetic analysis. Whole genomic DNA was extracted following an adapted phenol-chloroform protocol [51].

2.2. Microsatellite Genotyping

All of the samples were genotyped at 12 previously characterized microsatellite loci [25,28]. These were polymerase chain reaction (PCR) amplified in four separate multiplexed reactions using a Type It Kit (Qiagen) with the following PCR profile: one cycle of 5 min at 95 °C; 30 cycles of 30 s at 94 °C, 90 s at the specified annealing temperature (Ta) and 30 s at 72 °C; and final elongation step of 15 min at 60 °C (see Table S1 for Ta). Fluorescently labeled PCR products were resolved by electrophoresis on an ABI 3730xl capillary sequencer and allele sizes were scored by three independent observers using GeneMarker v. 2.6.2 (Softgenetics®). Samples that failed to genotype at four or more loci were excluded from subsequent analyses.

2.3. Genetic Summary Statistics

The R package Pegas v. 0.12 [52] was used to test for deviations from Hardy-Weinberg equilibrium at each locus using 1000 Monte Carlo replicates, while Genepop on the Web [53,54] was used to check for deviations from linkage equilibrium using default parameters. The resulting p-values were corrected for the table-wide false discovery rate (FDR) according to the procedure described by Benjamini & Hochberg [55]. Next, the R package diveRsity v. 1.9.90 [56] was used to calculate the number of alleles (Na), allelic richness (Ar), observed heterozygosity (Ho) and expected heterozygosity (He). Finally, we calculated standardized multilocus heterozygosity (sMLH) for each individual using the R package inbreedR v. 0.3.2 [57].

2.4. Analysis of Population Structure

To test for the presence of population structure, we implemented a number of different approaches. First, we calculated pairwise Fst values using Arlequin v. 3.5.2.2 [58], where statistical significance was determined based on 1000 permutations of the dataset. The resulting values were then used together with measures of shortest coastline distances among populations to test for the presence of isolation-by-distance by implementing a Mantel test. Second, we conducted a principal component analysis (PCA) of the microsatellite dataset using the R package adegenet v. 2.1.1 [59,60]. Third, we used Structure v. 2.3.3 [61] to carry out a Bayesian cluster analysis. Structure uses an iterative approach to group individuals into K groups by dividing the dataset in such a way that maximizes Hardy-Weinberg and linkage equilibria within the resulting groups. Each individual is then attributed a group membership value (Q) that varies from 0 to 1, with the latter indicating full group membership. We ran five simulations for each value of K between one and 10. We set the burn-in period and Markov chain Monte Carlo repetitions to 105 and 106, respectively. The most likely number of genetic groups was evaluated using the maximal average value of Ln P (D), a model choice criterion that estimates the posterior probability of the data, and the ΔK procedure described by Evanno et al. [62].

2.5. Elliptic Fourier Analysis of Shell Outlines

A geometric morphometric approach [63] based on elliptic Fourier analysis (EFA) of the shell outlines [64] was used to describe shape variation both within and among populations. Geometric information was extracted from the shell outlines and described as periodic function [32] through decomposition into the harmonic sum of progressively simplified trigonometric functions [65]. Low-frequency harmonics were used to approximate coarse-scale variation in the outlines, whereas higher-frequency harmonics captured fine-scale variation [32].
Outlines of orthogonal lateral and anterior views of the left shell valves (total n = 262; Figure 1b) were digitized and used as input data. The outlines for both views were processed independently, geometrically aligned and later combined for analysis following the protocol of Telesca et al. [31]. We then implemented an EFA on the resulting coordinates from shapes invariant to outline size, translation and rotation. After calibration, we chose seven harmonics to retain 99% of the cumulative harmonic power [66] (Figure S2). Four coefficients per harmonic (28 descriptors per view) were then extracted for each shell outline and used as variables quantifying geometric information [31,67].
A PCA was performed on the matrix of harmonic coefficients to characterize shell shape variation among individuals. The first three principal components (PCs), capturing 90.1% of the total shape variance, were used as new shape variables and analyzed using multivariate analysis of variance (MANOVA) to test for significant effects of location of origin and shell length (size) on shape variances. To visualize shell outline differences at the extremes of the morphospace, we generated deformation grids and iso-deformation lines through mathematical formalization of thin plate spline (TPS) analysis [68]. Shell morphometric analyses were carried out using Momocs v1.2.9 [32].
Generalized linear mixed models (GLMMs) were then used to explain variation in shell shape with respect to a number of predictors. Specifically, the scores of the first three shape PCs were modelled as a function of standardized latitude, the first two genetic PCs, sMLH, shape (a categorical variable with three levels: shape PC1, shape PC2, shape PC3) and their two-way interactions. Shell length was also included in the model to control for possible effects of within-population size variation, and collection site was fitted as random effect. Scores from the first three shape PCs were used as response variables within the same model to account for the interdependence of multiple shape variables that simultaneously describe variation in shell outlines as a whole [31]. This model was then optimized by rejecting non-significant interaction terms and factors that minimized the AICc value (Table S2). The final model was then of the form:
S h e l l S h a p e i j k = L a t i t u d e i k + L e n g t h i k + g P C 1 i k + g P C 2 i k + s M L H i k + S h a p e P C j + L e n g t h i k × L a y e r j + S i t e i j + ε i j k   S i t e i j   ~   N 0 ;   σ S i t e 2   ε i j k ~   N 0 ; σ 2
where ShellShapesijk is the kth thickness observation from PCj(j = PC1, PC2, PC3) and site i (i = 1, …, 9), Siteij is the random intercept and εijk is the error, which are assumed to be normally distributed with an expectation of zero and variances σ S i t e 2 and σ 2 respectively. Mean effect sizes of the predictor variables were estimated from the optimal model fitted on standardized variables [69]. Due to the difficulty of reliably estimating p-values in mixed models [70], we considered as significant any effect whose 95% confidence interval (CI) did not overlap zero. CIs were generated using a bias-corrected parametric bootstrap approach with 10,000 iterations of the data.

2.6. Analysis of Shell Thickness

A total of 167 shells were characterized in terms of both total thickness and the thickness of the three individual layers (PR, CL and CCL). The left shell valves were set in polyester resin (Kleer-Set FF, MetPrep, Coventry, UK) blocks and sliced longitudinally along their axis of maximum growth (Figure 1c) using a diamond saw. Shell cross-sections were progressively polished with silicon carbide paper (grit size: from P800 to P2500) and diamond paste (grading: from 9 µm to 3 µm). Polished cross-sections were treated with Mutvei’s solution [33] to highlight different growth marks and microstructures within the shells (Figure 1c). The thickness of each shell layer was measured on photographs of polished sections that were acquired using a stereo-microscope (Leica M165 C equipped with a Leica DFC295 HD camera, Leica, Wetzlar, Germany). Since larger individuals had undergone environmental abrasion, which removed the granular prismatic layer near the umbo, we measured the thickness of the whole shell and individual layers at the midpoint along the shell cross-section.
Next, we constructed two separate GLMMs to investigate how various predictor variables were associated with whole shell thickness and the thickness of the three individual layers. First, whole shell thickness was modelled as a function of standardized latitude, genetic PC1 and PC2, sMLH, shell length (which was included to control for possible effects of within-population size variation on layer thickness) and collection site, which was included as random effect. In the second model, the thicknesses of the three shell layers were combined into a single response variable and the following predictor variables were analyzed: standardized latitude, the first two genetic PCs, sMLH, shell layer (a categorical variable with three levels: PR, CL, CCL) and their two-way interactions. Shell length and collection site were also included, as in the model of whole shell thickness. The thickness of the three layers was analyzed within the same GLMM to simultaneously model common and divergent effects on each layer as well as to reduce the probability of type I error. Non-significant interaction terms were then rejected and only factors that minimized the AICc value (Table S2) were retained in the final model. This was of the form:
T h i c k n e s s i j k = L a t i t u d e i k + L e n g t h i k + g P C 1 i k + g P C 2 i k + s M L H i k + L a y e r j + L a t i t u d e i k × L a y e r j + L e n g t h i k × L a y e r j + S i t e i j + ε i j k   S i t e i j   ~   N 0 ;   σ S i t e 2   ε i j k ~   N 0 ; σ j 2 × e 2 δ × g P C 2 i k
where Thicknessijk is the kth thickness observation from layer j (j = PR, CL, CCL) and site i (i = 1, …, 9), Siteij is the random intercept for layer j, which is assumed to be normally distributed with an expectation of zero and variance σ S i t e 2 . εijk is the normally distributed error with a mean of zero and different variances per layer j changing exponentially with the variance covariate gPC2ik. This variance structure was chosen over others because it minimized the AICc value (Table S2). Mean effect sizes of predictor variables were estimated following the same procedure described for the GLMM on shell shape (see Section 2.7).

2.7. Analysis of Organic Shell Content

We performed thermogravimetric analyses (TGA) to estimate the weight proportion (wt%) of organic matrix within the two dominant shell layers, PR and CCL. A random subsample of five M. arenaria specimens were selected from two populations Lisbon (LIS) and Saint Andrews (SAN) to test for differences in shell organic content between low and high latitudes. We removed the periostracum by sanding before isolating pieces of PR and CCL. Shell pieces were cleaned, air-dried and then finely ground. We tested 10 milligrams of this powdered shell with a thermogravimetric analyzer (TGA Q500, TA Instruments, New Castle, DE, USA). The samples were subjected to constant heating from ~25 °C to 700 °C at a linear rate of 10 °C min−1 under a dynamic nitrogen atmosphere and weight changes were recorded. The proportion of organic matrix (wt%) was then estimated as the proportion of weight loss during the thermal treatment between 150 °C and 550 °C (Figure S3). The wt% of organic matrix in the PR and CCL layers (n = 20) was then modeled as a function of collection site and shell layer to test for latitudinal differences in shell organic content. Pairwise contrasts with a standard Bonferroni correction were then used to test for differences in wt% between the layers both within and between sampling locations.

2.8. Data accessibility

The microsatellite data used to generate this study are publibly available from the Figshare repository (https://figshare.com) (DOI: 10.6084/m9.figshare.11949546).

3. Results

In order to investigate patterns of genetic and morphological variation in M. arenaria, we collected a total of 266 samples from nine locations spanning a European latitudinal cline. A total of 247 individuals were successfully genotyped at 12 microsatellite loci. Genetic variability was relatively high, with each locus carrying on average 17 alleles and observed heterozygosity averaging 0.69 (Table 1 and Table S1). A number of significant deviations from Hardy-Weinberg equilibrium (HWE) were found after table-wide correction for the false discovery rate (Table S1), but the vast majority of loci deviated from HWE in fewer than three populations so null alleles or genotyping errors are unlikely to be responsible. Nevertheless, locus Ma26 showed significant deviations from HWE in seven populations and we therefore took the conservative measure of excluding this marker from subsequent analyses.

3.1. Population Structure

Pairwise Fst values varied between 0.02 and 0.14 (Table S3), with comparisons involving the Italian population (ITA) producing the highest overall values (mean = 0.12). The majority of Fst values were statistically significant after FDR correction (Table S3). A significant pattern of isolation-by-distance was also obtained for the full dataset (Mantel’s r = 0.67, p < 0.05). Significance was lost after excluding ITA (Mantel’s r = 0.15, p = 0.3) but the overall relationship became significant again after both ITA and LIS were excluded (Mantel’s r = 0.52, p < 0.05). Results of the PCA confirmed the greater magnitude of differentiation of the Italian population by clearly resolving individuals from ITA as a separate cluster (Figure 3a), while the other eight populations were distributed more or less as a continuum along the second principal component axis.
Arguably the most powerful tests of population structure need not rely on knowledge of the populations from which individuals were sampled. We therefore used the program Structure to identify the most likely number of genetic groups (K) within our dataset and to derive group membership coefficients (Q) for each individual. The number of groups is often identified using the maximal value of Ln P(D), a model-choice criterion that estimates the posterior probability of the data. However, Ln P(D) often plateaus or continues to increase slightly after the true value of K has been reached [62]. Our data yielded just such a pattern, with Ln P(D) rising steeply until around K = 4 and then increasing gradually towards a peak at K = 9 (Figure 3b). We therefore calculated ΔK, an ad hoc statistic based on the second order rate of change of the likelihood function with respect to K that has been shown by Evanno et al. [62] to be effective at capturing the uppermost hierarchical level of population structure under most circumstances. ΔK peaked at three (Figure 3b), indicating support for three main genetic groups corresponding to (i) Italy; (ii) Brest, Plymouth and St. Andrews; and (iii) the remaining North Sea populations, Kiel and Lisbon (Figure 3c). Increasing K to four additionally resolved St. Andrews as a separate group (Figure 3d), while further increases in K did not appreciably change the overall pattern.

3.2. Shell Shape Variation

PCA performed on harmonic coefficients revealed marked variation of shell outlines among individuals in both lateral and anterior views. The first three PCs accounted for 90.1% of the shape variability and a scatterplot of PC1 versus PC2 revealed appreciable variation among the nine populations across the morphospace (Figure 4). Multivariate analysis of variance (MANOVA) indicated a significant influence of collection site on shell shape (Wilk’s λ = 0.51, approx. F8, 252 = 7.85, p < 0.01) but there was no effect of shell size (Wilk’s λ = 0.99, approx. F1, 252 = 0.60, p = 0.62).
Subsequently, we attempted to identify specific shell features making the greatest contributions towards shell shape variation through comparisons of outlines at the extremes of the morphospace (Figure S4). We furthermore decomposed shell shape variation according to the contributions of the first three PCs through visual inspection of shell outlines constructed for increasing values of each PC (Figure S5). Finally, we used the mean shape and TPS analyses to illustrate the main outline deformation required to pass from one extreme of the morphospace to the other (i.e., from population Plymouth (PLY) to SAN, Figure S6). We found that PLY was characterized by more elongated and deeper shells than SAN, which exhibited rounder shells and flatter anterior view profiles.
Finally, we constructed a generalized linear mixed model (GLMM) to explore the effects of latitude, body size and genotype (expressed as genetic PC1 and PC2 scores and individual standardized multilocus heterozygosity, sMLH) on shell shape variation, as described in the Materials and methods. Shell shape was not significantly associated with latitude or genotype (Figure 5a, Table 2), but a significant association was found between shell length (a proxy for age) and the third morphological PC (p < 0.01, effect size = 0.27, 95% CI = 0.07 to 0.45). Consequently, shell shape appears not to be influenced by any of the predictor variables with the possible exception of a weak effect of age.

3.3. Variation in Shell Thickness

In order to investigate how shell thickness may be related to latitude, shell length and genotype, we constructed two GLMMs, the first of whole shell thickness and the second of the thickness of the individual layers (see Materials and methods for details). Variation in whole shell thickness was not associated with any of the predictor variables, with the exception of shell length (Figure 5b, Table S4), indicating an apparent absence of any latitudinal or genetic effects but a likely positive effect of age. In the GLMM of the thickness of the individual shell layers, we again detected a significant influence of shell length but no effect of genotype (Figure 5b, Table 2). However, this time, significant effects of latitude on the individual shell layers were detected, both in the form of a main effect of latitude and an interaction between latitude and shell layer. Specifically, the PR layer was proportionately thicker at higher latitudes and the CCL layer was proportionately thicker at lower latitudes, while the thickness of the CL layer did not appear to vary with latitude (Figure 5c).

3.4. Variation in Organic Content

To investigate differential patterns of organic content deposition in the two main shell layers at different latitudes, we quantified the proportion of organic matrix in the PR and CCL layers for each of five individual shells taken from the two extremes of the latitudinal range (LIS and SAN). The PR layer was characterized by a significantly higher wt% of organic content in shells from warm temperate regions in comparison to cold temperate regions (mean difference = 0.53 wt%, z = 3.69, p < 0.01), indicative of increasingly calcified prismatic layers at high latitudes (Figure 5d). By contrast, no significance difference was found in the organic content of the CCL layer (mean difference = −0.11 wt%, z = −0.86, p = 0.83). In addition, significant differences in the organic content of the two layers were detected in the low-latitude population (mean difference = 0.84 wt%, z = 6.01, p < 0.01), but not in the high-latitude population (mean difference = 0.20 wt%, z = 1.49, p = 0.44).

4. Discussion

We used microsatellites to characterize the population genetic structure of M. arenaria across Europe and to test for associations between genetic variables and shell morphological traits. We uncovered evidence for strong population genetic structure, which was unrelated to variation in shell shape, thickness, microstructure and organic content. Instead, although none of our predictor variables explained variation in shell shape, latitude appeared to be the best predictor of variation in shell thickness and organic content. We therefore conclude that most if not all of the observed variation in shell shape and thickness is probably due to environmentally driven phenotypic plasticity.

4.1. Population Genetic Structure

Our data are suggestive of the presence of four main genetic groups situated in the Adriatic, the Celtic Sea, on the east coast of Scotland and along the continental coast of the North Sea and the entrance of the Baltic Sea. This observation lends support to a previous study by Cross et al. [27], that documented significant genetic differences between three populations from the British Isles and a single population from the Netherlands. Unfortunately, a direct comparison of the two studies is not possible as we were unable to source samples from the same localities as Cross et al. [27]. This leaves open the question of whether additional genetic groups might be detected around the British Isles as well as in other localities that were not included in the current study, such as the Baltic and White Sea.
Although our sampling design was far from exhaustive, the broad geographical coverage of our study allowed us to capture a number of interesting patterns. First and most obviously, samples from the Adriatic Sea (ITA) were strongly differentiated from the Atlantic populations. Deep genetic divergence between the Mediterranean and the Atlantic is a common pattern found in native European marine invertebrates [71,72], which typically results from a combination of historical vicariance and the cessation of contemporary gene flow due to the presence of the Almería-Oran oceanographic front. However, this is unlikely to provide an explanation in the case of M. arenaria for two reasons. First, the soft-shell clam is believed to have gone extinct in Europe during the last glacial maximum [7], which would mean that events pre-dating this period could not have left a genetic trace. Second, after this species was reintroduced into Europe, anthropogenic translocations were commonplace and we are aware of at least two documented instances in which M. arenaria was introduced into the Mediterranean [11,12]. Additionally, Lasota et al. [5] suggested that one potential origin of Adriatic M. arenaria populations could be the Atlantic coast of North America, which would be consistent with the large magnitude of genetic differentiation observed in the current study. To investigate this further, more extensive geographical sampling would need to be combined with genetic assignment testing in order to evaluate the most likely source of origin(s) of the Mediterranean M. arenaria populations. Genomic data would also be desirable as these might allow divergence times to be estimated and the strength of support for alternative recolonization pathways to be evaluated sensu [20].
The pattern of genetic structure we uncovered among the Atlantic M. arenaria populations may superficially resemble that observed in other marine species native to Europe, where Atlantic populations are divided into a northern and a southern genetic cluster [20,73,74]. However, this apparent similarity cannot have originated from the same processes. In the case of marine species native to Europe, this pattern emerged as a consequence of natural postglacial recolonization, either because of the presence of separate glacial refugia in northern and southern Europe or as a result of a founder effect event that occurred during northward recolonization from a single refugium located in southern Europe [74]. By contrast, the recolonization of Europe by M. arenaria was at least partially anthropogenic, following its introduction to Denmark between the 13th and 15th centuries [8,9]. This recolonization is known to have proceeded both northward and southward, providing a possible explanation for the overall pattern of isolation by distance in our data.
Finally, we could show that M. arenaria samples from Portugal showed high genetic affinity with populations from along the continental coast of the North Sea and the entrance of the Baltic Sea (Le Crotoy (LCT), Balgzand (TEX), Sylt (SYL) and Kiel (KIE)). This finding is again consistent with the notion that man-mediated translocation played an important role in shaping the genetic structure of M. arenaria across Europe. In this particular case, the most parsimonious explanation for the observed pattern would be that soft-shell clams from the North Sea coast were introduced into Portugal, either deliberately or inadvertently due to the fact that M. arenaria can be easily transported with ballast water [4,5]. In line with the first of these explanations, Conde et al. [13] have already argued that soft-shell clams from Lisbon, as well as from other two Portuguese populations, may have originated as a consequence of intentional introductions.

4.2. Shell Shape Variation

An EFA of shell outlines uncovered differences in shell shape among the collection sites, but none of the predictors fitted in our GLMM could explain a significant proportion of the total variation. Although it is not necessarily surprising that the genetic variables were unrelated to variation in shell shape, studies of other bivalve species have reported strong effects of local environmental conditions, substrate type and predation pressure [40,41,44,75]. Consequently, we were somewhat surprised not to have detected any influence of latitude on shell shape in M. arenaria. One possible explanation could be that our predictor variables were too crude to capture meaningful variation in key abiotic or biotic variables. In particular, although our broad geographic coverage maximized variation in a suite of climatic and other variables, this may have hindered the detection of relatively subtle effects. This issue could potentially be circumvented by sampling populations over a much finer geographic scale, as this would effectively control for large-scale sources of variability while facilitating more quantitative investigation of specific factors such as substrate typology and sediment size.

4.3. Variation in Shell Thickness, Microstructure and Organic Content

Geographic patterns of molluscan skeletal production have typically been explained by two paradigms: poleward reduction of predation pressure [76,77] and increased calcification costs at high latitudes, which results from a combination of decreased CaCO3 saturation and reduced metabolic rates [47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79]. For M. arenaria, predation acts most heavily upon young individuals that are present only in the upper layer of the substrate, while predation risk for older individuals buried deeper in the substrate is negligible [80,81]. Given the size classes analyzed in this study, it could therefore be argued that predation pressure is unlikely to have influenced shell morphology, although we cannot discount the possibility that variation in the risk of predation during early life could have played a role.
Surprisingly, our results are also not consistent with poleward skeletal reductions due to increased calcification costs at high latitudes, as we did not observe any change in total shell thickness along a latitudinal cline. However, shell organic content was reduced in samples from the more northerly, colder environment, which is in line with suggestions that shell organics have higher production costs than CaCO3 deposition [46,47]. Similar divergent patterns have been documented for laternulid clams [82], where thicker shells have been suggested to have a selective advantage when individuals are moved through sediment by ice disturbance [83]. Likewise, the calcification pattern we observe might represent a trade-off to preserve a constant shell thickness across latitude. Specifically, the production of less energetically-expensive shell layers may be favored at high latitudes as a means of enhancing protection from physical disturbance.

4.4. Genetic Versus Plastic Contributions

Previous studies of a variety of marine and freshwater invertebrates have provided support for a plastic nature of shell morphology. This conclusion has been reached either due to the absence of a link between population genetic structure and morphological variation [30,34,36] or based on the results of reciprocal transplant experiments, which have uncovered a prominent role of environmental variation in shaping shell morphology [84,85]. We built upon the approaches used in previous population genetic studies by integrating genetic variation in the form of principal component scores, which capture multiple aspects of genetic variation and facilitate morphological-genetic comparisons at the level of individuals rather than populations. Furthermore, as heterozygosity is associated with morphological variation in several bivalve species [86,87,88], we included individual sMLH as a predictor in our models. We found no main effects of any of these genetic variables on either shell shape, thickness or organic content, providing evidence for a primary role of phenotypic plasticity, although genomic approaches capable of generating tens of thousands of genetic markers should offer greater power for testing for gene-phenotype associations [36].

5. Conclusions

We used microsatellites to characterize the population genetic structure of M. arenaria across Europe and to relate this to latitudinal variation in shell shape, thickness, microstructure and organic content. We uncovered strong population structure, consistent with the known involvement of humans in reintroducing this species to Europe as well as a long history of both deliberate and unintentional long-distance translocations. Additionally, we were unable to find any genetic effects on individual variability in shell shape and thickness, consistent with previous studies of other shellfish species reaching the conclusion that shell morphology is largely plastic [34,35,36]. Specifically, the best predictor of the thickness of individual shell layers in M. arenaria was latitude, which is associated with variation in numerous variables, both biotic and abiotic.

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4425/11/3/298/s1. Figure S1: Scanning Electron Microscope (SEM) images of M. arenaria shell microstructure. Shell cross-sections were progressively polished (up to 1μm), ultrasonically cleaned and air-dried prior to mounting and sputter coating (Emitech K550X). Images were taken using a Quanta-650F SEM (Department of Earth Sciences, University of Cambridge, UK). The microstructural nomenclature of Bieler et al. [45] has been used. Scale bars are shown on the bottom right of each image. (a) The Periostracum, indicated by white arrows (b) the outer granular prismatic layer, (c) the middle crossed-lamellar layer and (d) the inner complex crossed-lamellar layer (comprising a sequence of crossed-lamellar and prismatic layers shown by white arrows). Figure S2: Results of the calibration procedure used for elliptic Fourier analysis (EFA) of shell outlines. Cumulative harmonic Fourier power is shown separately for (a) lateral and (b) anterior shell views. The power is proportional to the harmonic amplitude and can be considered a measure of shape information [63]. We evaluated the appropriate number of harmonics to retain so that their cumulative power captured 99% of the total Fourier power. Average shell shapes reconstructed for different numbers of harmonics (1, 3, 5, 7 and 9) are shown. Figure S3: Example outcomes of Thermal Gravimetric Analysis (TGA, green line) and Derivative Thermogravimetry (DTG, blue line). The TGA curve represents weight changes with increasing treatment temperature for the granular prismatic layer of M. arenaria. Four regions of weight loss with increasing temperature are highlighted: (i) the evaporation of physically adsorbed water at 30‒150 °C; (ii) the degradation of extra-crystalline organic matrix at 150‒400 °C; (iii) the release of intra-crystalline organics at 400–550 °C; and (iv) the rapid decomposition of calcium carbonate (CaCO3) into calcium oxide (CaO) and carbon dioxide (CO2) starting at ~550 °C. The DTG line represents the derivative of the thermal curve and shows the rate of weight loss during heating. The peak indicates the temperature at which the organic mass loss was fastest. Figure S4: Contributions of PCs toward variation in shell shape for increasing PC values: mean −3 SD (blue), mean (black) and mean +3 SD (red). PC1 contributed mainly towards variation in shell roundness and depth. By contrast, PC2 described variation in shell roundness and in the symmetry of the anterior view profile. PC3 contributed towards minor variation in the lateral view profile. Figure S5: Contribution of the first four shape variables (PCs) to shape variation. Average shell shapes for the lateral and 3anterior view are shown for increasing values along each PC (Mean −3 SD, Mean, Mean +3 SD) and shapes at the extremes of each variable are compared (Mean ±3 SD). Figure S6: Patterns of shell shape variations in European M. arenaria samples. (a) Mean shell shapes for each collection site. (b) Differences between mean shapes at the extremes of the morphospace represented through (i) iso-deformation lines (bottom), representing the outline regions subjected to different degrees of change (blue: low deformation; red: strong deformation), and (ii) deformation grids (top), depicting the bindings required to pass from an extreme (PLY) to another (SAN). Table S1: Details of the 12 microsatellite loci used in this study together with their annealing temperatures and polymorphism characteristics in 247 soft-shell clams sampled from nine different populations. Table S2: Alternative models of shell shape and shell layer thickness for optimal fixed and variance structures (shown in bold). Degrees of freedom (k), log likelihood estimation, corrected AICc values and likelihood estimation methods are reported for each model. Table S3: Pairwise FST values (below diagonal) and corresponding p-values (above diagonal) calculated from 11 microsatellites. Bold p-values were significant only before table-wide FDR correction, while bold and underlined p-values also remained significant after FDR correction. Table S4: Summary of the results of the GLMM of whole shell thickness. Estimated statistics bootstrapped 95% CIs for regression parameters are reported for the modelled relationships between shell thickness, latitude, shell length and genetic variables. Regression parameters were considered statistically significant and highlighted in bold when the bootstrapped 95% CI did not overlap zero.

Author Contributions

Conceptualization, M.D.N, L.T., D.L.J.V., E.M.H. and J.I.H.; methodology, M.D.N., L.T., D.L.J.V. and H.K.G.; formal analysis, M.D.N., L.T., D.L.J.V. and H.K.G.; data curation, M.D.N., L.T., D.L.J.V., H.K.G., G.C., E.M.H. and J.I.H.; writing—original draft preparation, L.T., D.L.J.V. and J.I.H.; writing—reviewing and editing, M.D.N., L.T., D.L.J.V., H.K.G., G.C., E.M.H. and J.I.H.; supervision, E.M.H. and J.I.H.; project administration, E.M.H. and J.I.H.; funding acquisition, E.M.H. and J.I.H. Original artworks: L.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the European Union Marie Curie Seventh Framework Programme under grant agreement no. 605051.

Acknowledgments

We are grateful to Kirti Ramesh (University of Gothenburg, Gothenburg, Sweden), Birgit Hussel (Alfred Wegener Institute, List/Sylt, Germany), Iain Johnston (Scottish Oceans Institute, St. Andrews, UK), Sarah Dashfield (Plymouth Marine Laboratory, Plymouth, UK), Henk van der Veer and Rob Dekker (Royal Netherlands Institute for Sea Research, Texel, Netherlands) and Frederico Batista (Centre for Environment, Fisheries and Aquaculture Science, UK) for providing M. arenaria samples. We acknowledge support for the Article Processing Charge by the Deutsche Forschungsgemeinschaft and the Open Access Publication Fund of Bielefeld University.

Conflicts of Interest

The authors declare no conflict of interest

References

  1. Strasser, M. Mya arenaria—An ancient invader of the North Sea coast. Helgoländer Meeresunters. 1999, 52, 309. [Google Scholar] [CrossRef] [Green Version]
  2. Hanna, G.D. Introduced mollusks of western North America. Occ. Pap. Calif. Acad. Sci. 1966, 48, 1–108. [Google Scholar]
  3. Carlton, J.T. Man’s role in changing the face of the ocean: Biological invasions and implications for conservation of near-shore environments. Conserv. Biol. 1989, 3, 265–273. [Google Scholar] [CrossRef]
  4. Leppäkoski, E. The Baltic and the Black Sea-Seriously contaminated by nonindigenous species? In Proceedings of the Conference and Workshop Nonindigenous Estuarine and Marine Organisms (NEMO), Seattle, WA, USA, April 1993; pp. 37–44. [Google Scholar]
  5. Lasota, R.; Pierscieniak, K.; Garcia, P.; Simon-Bouhet, B.; Wolowicz, M. Large-scale mitochondrial COI gene sequence variability reflects the complex colonization history of the invasive soft-shell clam, Mya arenaria (L.)(Bivalvia). Estuar. Coast. Shelf Sci. 2016, 181, 256–265. [Google Scholar] [CrossRef]
  6. Strauch, F. Phylogenese, Adaptation und Migration einiger nordischer-mariner Molluskengenera (Neptunea, Panomya,-Cyrtodaria und Mya). Abh. Senckenberg. Naturforsch. Ges. 1972, 531, 1–211. [Google Scholar]
  7. Hessland, I. On the quaternary Mya period in Europe. Ark. Zool. 1946, 37A, 1–51. [Google Scholar]
  8. Petersen, K.S.; Rasmussen, K.L.; Heinemeier, J.; Rud, N. Clams before Columbus? Nature 1992, 359, 679. [Google Scholar] [CrossRef]
  9. Essink, K.; Oost, A.P. How did Mya arenaria (Mollusca; Bivalvia) repopulate European waters in mediaeval times? Mar. Biodivers. 2019, 49, 1–10. [Google Scholar] [CrossRef]
  10. Reise, K.; Gollasch, S.; Wolff, W.J. Introduced marine species of the North Sea coasts. Helgoländer Meeresunters. 1998, 52, 219. [Google Scholar] [CrossRef] [Green Version]
  11. Crocetta, F.; Turolla, E. Mya arenaria Linné, 1758 (Mollusca: Bivalvia) in the Mediterranean Sea: Its distribution revisited. J. Biol. Res. 2011, 16, 188–193. [Google Scholar]
  12. Zenetos, A.; Koutsoubas, D.; Vardala-Theodorou, E. Origin and vectors of introduction of exotic molluscs in Greek waters. Belg. J. Zool. 2004, 134, 161–168. [Google Scholar]
  13. Conde, A.; Novais, J.M.; Domínguez, J. The presence of Mya arenaria in the Ria de Aveiro is the third confirmed record of this invasive clam on the Portuguese coast. Mar. Biodivers. Rec. 2012, 5. [Google Scholar] [CrossRef] [Green Version]
  14. Nicolas, D.; Le Loc’h, F.; Désaunay, Y.; Hamon, D.; Blanchet, A.; Le Pape, O. Relationships between benthic macrofauna and habitat suitability for juvenile common sole (Solea solea, L.) in the Vilaine estuary (Bay of Biscay, France) nursery ground. Estuar. Coast. Shelf Sci. 2007, 73, 639–650. [Google Scholar] [CrossRef] [Green Version]
  15. Hughes, R.G.; Lloyd, D.; Ball, L.; Emson, D. The effects of the polychaete Nereis diversicolor on the distribution and transplanting success of Zostera noltii. Helgol. Mar. Res. 2000, 54, 129. [Google Scholar] [CrossRef] [Green Version]
  16. Seaward, D.R. Distribution of the Marine Molluscs of North West Europe (No. 165); Joint Nature Conservation Committee: Peterborough, UK, 1990. [Google Scholar]
  17. Ysebaert, T.; Meire, P.; Coosen, J.; Essink, K. Zonation of intertidal macrobenthos in the estuaries of Schelde and Ems. Aquat. Ecol. 1998, 32, 53–71. [Google Scholar] [CrossRef]
  18. Maximovich, N.V.; Guerassimova, A.V. Life history characteristics of the clam Mya arenaria in the White Sea. Helgol. Mar. Res. 2003, 57, 91. [Google Scholar] [CrossRef]
  19. Vendrami, D.L.; Houston, R.D.; Gharbi, K.; Telesca, L.; Gutierrez, A.P.; Gurney-Smith, H.; Hoffman, J.I. Detailed insights into pan-European population structure and inbreeding in wild and hatchery Pacific oysters (Crassostrea gigas) revealed by genome-wide SNP data. Evol. Appl. 2019, 12, 519–534. [Google Scholar] [CrossRef] [Green Version]
  20. Vendrami, D.L.; De Noia, M.; Telesca, L.; Handal, W.; Charrier, G.; Boudry, P.; Hoffman, J.I. RAD sequencing sheds new light on the genetic structure and local adaptation of European scallops and resolves their demographic histories. Sci. Rep. 2019, 9, 7455. [Google Scholar] [CrossRef]
  21. Zbawicka, M.; Drywa, A.; Śmietanka, B.; Wenne, R. Identification and validation of novel SNP markers in European populations of marine Mytilus mussels. Mar. Biol. 2012, 159, 1347–1362. [Google Scholar] [CrossRef]
  22. Launey, S.; Ledu, C.; Boudry, P.; Bonhomme, F.; Naciri-Graven, Y. Geographic structure in the European flat oyster (Ostrea edulis L.) as revealed by microsatellite polymorphism. J. Hered. 2002, 93, 331–351. [Google Scholar] [CrossRef]
  23. Daguin, C.; Bonhomme, F.; Borsa, P. The zone of sympatry and hybridization of Mytilus edulis and M. galloprovincialis, as described by intron length polymorphism at locus mac-1. Heredity 2001, 86, 342. [Google Scholar] [CrossRef] [PubMed]
  24. Strasser, C.A.; Barber, P.H. Limited genetic variation and structure in softshell clams (Mya arenaria) across their native and introduced range. Conserv. Genet. 2009, 10, 803. [Google Scholar] [CrossRef] [Green Version]
  25. St-Onge, P.; Sévigny, J.M.; Strasser, C.; Tremblay, R. Strong population differentiation of softshell clams (Mya arenaria) sampled across seven biogeographic marine ecoregions: Possible selection and isolation by distance. Mari. Biol. 2013, 160, 1065–1081. [Google Scholar] [CrossRef]
  26. Lasota, R.; Hummel, H.; Wolowicz, M. Genetic diversity of European populations of the invasive soft-shell clam Mya arenaria (Bivalvia). J. Mar. Biol. Assoc. U.K 2004, 84, 1051–1056. [Google Scholar] [CrossRef] [Green Version]
  27. Cross, M.E.; Bradley, C.R.; Cross, T.F.; Culloty, S.; Lynch, S.; McGinnity, P.; Prodöhl, P.A. Genetic evidence supports recolonisation by Mya arenaria of western Europe from North America. Mar. Ecol. Prog. Ser. 2016, 549, 99–112. [Google Scholar] [CrossRef] [Green Version]
  28. Krapal, A.M.; Popa, O.P.; Iorgu, E.I.; Costache, M.; Popa, L.O. Isolation and characterization of new microsatellite markers for the invasive softshell clam, Mya arenaria (L.)(Bivalvia: Myidae). Int. J. Mol. Sci. 2012, 13, 2515–2520. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Chen, C.; Zhang, H.; Wang, A.; Lu, M.; Shen, Z.; Lian, C. Phenotypic plasticity accounts for most of the variation in leaf manganese concentrations in Phytolacca americana growing in manganese-contaminated environments. Plant Soil 2015, 396, 215–227. [Google Scholar] [CrossRef]
  30. Zieritz, A.; Hoffman, J.I.; Amos, W.; Aldridge, D.C. Phenotypic plasticity and genetic isolation-by-distance in the freshwater mussel Unio pictorum (Mollusca: Unionoida). Evol. Ecol. 2010, 24, 923–938. [Google Scholar] [CrossRef]
  31. Telesca, L.; Michalek, K.; Sanders, T.; Peck, L.S.; Thyrring, J.; Harper, E.M. Blue mussel shell shape plasticity and natural environments: A quantitative approach. Sci. Rep. 2018, 8, 2865. [Google Scholar] [CrossRef] [Green Version]
  32. Bonhomme, V.; Picq, S.; Gaucherel, C.; Claude, J. Momocs: Outline analysis using R. J. Stat. Softw. 2014, 56, 1–24. [Google Scholar] [CrossRef] [Green Version]
  33. Schone, B.R.; Dunca, E.; Fiebig, J.; Pfeiffer, M. Mutvei’s solution: An ideal agent for resolving microgrowth structures of biogenic carbonates. Palaeogeogr. Palaeoclimatol. Palaeoecol. 2005, 228, 149–166. [Google Scholar] [CrossRef]
  34. Hoffman, J.I.; Peck, L.S.; Hillyard, G.; Zieritz, A.; Clark, M.S. No evidence for genetic differentiation between Antarctic limpet Nacella concinna morphotypes. Mar. Biol. 2010, 157, 765–778. [Google Scholar] [CrossRef]
  35. Trussell, G.C.; Etter, R.J. Integrating genetic and environmental forces that shape the evolution of geographic variation in a marine snail. In Microevolution rate, Pattern, Process; Springer: Dordrecht, The Netherlands, 2011; pp. 321–337. [Google Scholar]
  36. Vendrami, D.L.; Telesca, L.; Weigand, H.; Weiss, M.; Fawcett, K.; Lehman, K.; Hoffman, J.I. RAD sequencing resolves fine-scale population structure in a benthic invertebrate: Implications for understanding phenotypic plasticity. R. Soc. Open Sci. 2017, 4, 160548. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Shields, J.L.; Barnes, P.; Heath, D.D. Growth and survival differences among native, introduced and hybrid blue mussels (Mytilus spp.): Genotype, environment and interaction effects. Mar. Biol. 2008, 154, 919–928. [Google Scholar] [CrossRef]
  38. Fuentes, J.; Reyero, I.; Zapata, C.; Alvarez, G. Influence of stock and culture site on growth rate and mortality of mussels (Mytilus galloprovincialis Lmk.) in Galicia, Spain. Aquaculture 1992, 105, 131–142. [Google Scholar] [CrossRef]
  39. Dickie, L.M.; Boudreau, P.R.; Freeman, K.R. Influences of stock and site on growth and mortality in the blue mussel (Mytilus edulis). Can. J. Fish. Aquat. Sci. 1984, 41, 134–140. [Google Scholar] [CrossRef]
  40. Fitzer, S.C.; Vittert, L.; Bowman, A.; Kamenos, N.A.; Phoenix, V.R.; Cusack, M. Ocean acidification and temperature increase impact mussel shell shape and thickness: Problematic for protection? Ecol. Evol. 2015, 5, 4875–4884. [Google Scholar] [CrossRef] [Green Version]
  41. Illesca, A.; Oyarzún, P.A.; Toro, J.E.; Gardner, J.P.A. Morphometric variability of smooth-shelled blue mussels from the Pacific coast of South America. Biol. J. Linn. Soc. 2018, 20, 1–16. [Google Scholar] [CrossRef]
  42. Brönmark, C.; Lakowitz, T.; Hollander, J. Predator-induced morphological plasticity across local populations of a freshwater snail. PLoS ONE 2011, 6, e21773. [Google Scholar] [CrossRef]
  43. Swan, E.F. The growth of the clam Mya arenaria as affected by the substratum. Ecology 1852, 33, 530–534. [Google Scholar] [CrossRef]
  44. Emerson, C.W. Influence of sediment disturbance and water flow on the growth of the soft-shell clam, Mya arenaria L. Can. J. Fish. Aquat. Sci. 1990, 47, 1655–1663. [Google Scholar] [CrossRef]
  45. Bieler, R.; Mikkelsen, P.M.; Collins, T.M.; Glover, E.A.; González, V.L.; Graf, D.L.; Harper, E.M.; Healy, J.; Kawauchi, G.Y.; Sharma, P.P.; et al. Investigating the Bivalve Tree of Life—An exemplar-based approach combining molecular and novel morphological characters. Invertebr. Syst. 2014, 28, 32. [Google Scholar] [CrossRef] [Green Version]
  46. Palmer, A.R. Calcification in marine molluscs: How costly is it? Proc. Natl. Acad. Sci. USA 1992, 89, 1379–1382. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  47. Watson, S.-A.; Morley, S.A.; Peck, L.S. Latitudinal trends in shell production cost from the tropics to the poles. Sci. Adv. 2017, 3, e1701362. [Google Scholar] [CrossRef] [Green Version]
  48. Leung, J.Y.S.; Russell, B.D.; Connell, S.D. Mineralogical plasticity acts as a compensatory mechanism to the impacts of ocean acidification. Environ. Sci. Technol. 2017, 51, 2652–2659. [Google Scholar] [CrossRef]
  49. Telesca, L.; Peck, L.S.; Sanders, T.; Thyrring, J.; Sejr, M.K.; Harper, E.M. Biomineralization plasticity and environmental heterogeneity predict geographical resilience patterns of foundation species to future change. Glob. Chang. Biol. 2019, 25, 4179–4193. [Google Scholar] [CrossRef] [PubMed]
  50. Filippenko, D.; Naumenko, E. Patterns of the growth of soft-shell clam Mya arenaria L. (Bivalvia) in shallow water estuaries of the southern Baltic Sea. Ecohydrol. Hydrobiol. 2014, 14, 157–165. [Google Scholar] [CrossRef]
  51. Sambrook, J.; Fritsch, E.F.; Maniatis, T. Molecular Cloning: A Laboratory Manual, 2nd ed.; Cold Spring Harbor Laboratory Press: Suffolk County, NY, USA, 1989. [Google Scholar]
  52. Paradis, E. pegas: An R package for population genetics with an integrated–modular approach. Bioinformatics 2010, 26, 419–420. [Google Scholar] [CrossRef] [Green Version]
  53. Raymond, M. GENEPOP (version 1.2): Population genetics software for exact tests and ecumenicism. J. Hered. 1995, 86, 248–249. [Google Scholar] [CrossRef]
  54. Rousset, F. genepop’007: A complete re-implementation of the genepop software for Windows and Linux. Mol. Ecol. Resour. 2008, 8, 103–106. [Google Scholar] [CrossRef]
  55. Benjamini, Y.; Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B (Methodol.) 1995, 57, 289–300. [Google Scholar] [CrossRef]
  56. Keenan, K.; McGinnity, P.; Cross, T.F.; Crozier, W.W.; Prodöhl, P.A. diveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors. Methods Ecol. Evol. 2013, 4, 782–788. [Google Scholar] [CrossRef] [Green Version]
  57. Stoffel, M.A.; Esser, M.; Kardos, M.; Humble, E.; Nichols, H.; David, P.; Hoffman, J.I. inbreedR: An R package for the analysis of inbreeding based on genetic markers. Methods Ecol. Evol. 2016, 7, 1331–1339. [Google Scholar] [CrossRef]
  58. Schneider, S.; Roessli, D.; Excoffier, L. Arlequin version 2.000; A Software for Population Genetics Data Analysis; Genetics and Biometry Laboratory, University of Geneva: Geneva, Switzerland, 2000. [Google Scholar]
  59. Jombart, T. adegenet: A R package for the multivariate analysis of genetic markers. Bioinformatics 2008, 24, 1403–1405. [Google Scholar] [CrossRef] [Green Version]
  60. Jombart, T.; Ahmed, I. adegenet 1.3-1: New tools for the analysis of genome-wide SNP data. Bioinformatics 2011, 27, 3070–3071. [Google Scholar] [CrossRef] [Green Version]
  61. Pritchard, J.K.; Stephens, M.; Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 2000, 155, 945–959. [Google Scholar] [PubMed]
  62. Evanno, G.; Regnaut, S.; Goudet, J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol. Ecol. 2005, 14, 2611–2620. [Google Scholar] [CrossRef] [Green Version]
  63. Adams, D.C.; Rohlf, F.J.; Slice, D.E. Geometric morphometrics: Ten years of progress following the “revolution”. Ital. J. Zool. 2004, 71, 5–16. [Google Scholar] [CrossRef] [Green Version]
  64. Kuhl, F.P.; Giardina, C.R. Elliptic Fourier features of a closed contour. Comput. Graph. Image Process. 1982, 18, 236–258. [Google Scholar] [CrossRef]
  65. Rohlf, F.J.; Archie, J.W. A comparison of Fourier methods for the description of wing shape in mosquitoes (Diptera: Culicidae). Syst. Zool. 1984, 33, 302. [Google Scholar] [CrossRef]
  66. Crampton, J.S. Elliptic Fourier shape analysis of fossil bivalves: Some practical considerations. Lethaia 1995, 28, 179–186. [Google Scholar] [CrossRef]
  67. Claude, J. Morphometrics with R; Springer: New York, NY, USA, 2008. [Google Scholar]
  68. Bookstein, F.L. Morphometric Tools for Landmark Data: Geometry and Biology; Cambridge University Press: Cambridge, UK, 1991. [Google Scholar]
  69. Schielzeth, H. Simple means to improve the interpretability of regression coefficients. Methods Ecol. Evol. 2010, 1, 103–113. [Google Scholar] [CrossRef]
  70. Bolker, B.M. Linear and generalized linear mixed models. In Ecological Statistics; Fox, G.A., Negrete-Yankelevich, S., Sosa, V.J., Eds.; Oxford University Press: Oxford, UK, 2015; pp. 309–333. [Google Scholar]
  71. Côrte-Real, H.B.S.M.; Hawkins, S.J.; Thorpe, J.P. Population differentiation and taxonomic status of the exploited limpet Patella candei in the Macaronesian islands (Azores, Madeira, Canaries). Mar. Biol. 1996, 125, 141–152. [Google Scholar] [CrossRef]
  72. Almaça, C. Evolutionary and zoogeographical remarks on the Mediterranean fauna of brachyuran crabs. In Mediterranean Marine Ecosystems; Springer: Boston, MA, USA, 1985; pp. 347–366. [Google Scholar]
  73. Krakau, M.; Jacobsen, S.; Jensen, K.T.; Reise, K. The cockle Cerastoderma edule at Northeast Atlantic shores: Genetic signatures of glacial refugia. Mar. Biol. 2012, 159, 221–230. [Google Scholar] [CrossRef] [Green Version]
  74. Maggs, C.A.; Castilho, R.; Foltz, D.; Henzler, C.; Jolly, M.T.; Kelly, J.; Viard, F. Evaluating signatures of glacial refugia for North Atlantic benthic marine taxa. Ecology 2008, 89, S108–S122. [Google Scholar] [CrossRef] [Green Version]
  75. Gardner, J.P.A.; Thompson, R.J. Influence of genotype and geography on shell shape and morphometric trait variation among North Atlantic blue mussel (Mytilus spp.) populations. Biol. J. Linn. Soc. 2009, 96, 875–897. [Google Scholar] [CrossRef]
  76. Aronson, R.B.; Thatje, S.; Clarke, A.; Peck, L.S.; Blake, D.B.; Wilga, C.D.; Seibel, B.A. Climate change and invasibility of the Antarctic benthos. Annu. Rev. Ecol. Evol. Syst. 2007, 38, 129–154. [Google Scholar] [CrossRef] [Green Version]
  77. Harper, E.M.; Peck, L.S. Latitudinal and depth gradients in marine predation pressure. Glob. Ecol. Biogeogr. 2016, 25, 670–678. [Google Scholar] [CrossRef]
  78. Clarke, A. Temperature and extinction in the sea: A physiologist’s view. Paleobiology 1993, 19, 499–518. [Google Scholar] [CrossRef]
  79. Watson, S.-A.; Peck, L.S.; Tyler, P.A.; Southgate, P.C.; Tan, K.S.; Day, R.W.; Morley, S.A. Marine invertebrate skeleton size varies with latitude, temperature and carbonate saturation: Implications for global change and ocean acidification. Glob. Change Biol. 2012, 18, 3026–3038. [Google Scholar] [CrossRef]
  80. Zwarts, L.; Wanink, J. Siphon size and burying depth in deposit- and suspension-feeding benthic bivalves. Mar. Biol. 1989, 100, 227–240. [Google Scholar] [CrossRef]
  81. Zwarts, L.; Wanink, J.H. How oystercatchers and curlews successively deplete clams. In Coastal Waders and Wildfowl in Winter; Evans, P., Goss-Custard, J., Hale, W., Eds.; Cambridge University Press: Cambridge, UK, 1984; pp. 69–83. [Google Scholar]
  82. Watson, S.-A. Latitudinal Gradients in Marine Invertebrate Shell Morphology: Production Costs and Predation Pressure; University of Southampton: Southampton, UK, 2009. [Google Scholar]
  83. Harper, E.M.; Clark, M.S.; Hoffman, J.I.; Philipp, E.R.E.; Peck, L.S.; Morley, S.A. Iceberg scour and shell damage in the Antarctic bivalve Laternula elliptica. PLoS ONE 2012, 7, e46341. [Google Scholar] [CrossRef] [Green Version]
  84. Hinch, S.G.; Bailey, R.C.; Green, R.H. Growth of Lampsilis radiata (Bivalvia: Unionidae) in sand and mud: A reciprocal transplant experiment. Can. J. Fish. Aquat. Sci. 1986, 43, 548–552. [Google Scholar] [CrossRef]
  85. Moore, H.B. The relation of shell growth to environment in Patella vulgata. J. Molluscan Stud. 1934, 21, 217–222. [Google Scholar]
  86. Borrell, Y.J.; Pineda, H.; McCarthy, I.; Vazquez, E.; Sanchez, J.A.; Lizana, G.B. Correlations between fitness and heterozygosity at allozyme and microsatellite loci in the Atlantic salmon, Salmo salar L. Heredity 2004, 92, 585. [Google Scholar] [CrossRef] [Green Version]
  87. Yezerinac, S.M.; Lougheed, S.C.; Handford, P. Morphological variability and enzyme heterozygosity: Individual and population level correlations. Evolution 1992, 46, 1959–1964. [Google Scholar] [CrossRef] [PubMed]
  88. Mitton, J.B. Relationship between heterozygosity for enzyme loci and variation of morphological characters in natural populations. Nature 1978, 273, 661–662. [Google Scholar] [CrossRef]
Genes 11 00298 g001 550
Figure 1. (a) Depiction of Mya arenaria in its natural habitat; (b) lateral and anterior shell views; (c) a dorsoventral cross-section of the left shell valve along the axis of maximum growth showing the internal shell structure. Shown are the granular prismatic (PR) layer, the crossed-lamellar (CL) layer and the complex crossed-lamellar (CCL) layer.
Figure 1. (a) Depiction of Mya arenaria in its natural habitat; (b) lateral and anterior shell views; (c) a dorsoventral cross-section of the left shell valve along the axis of maximum growth showing the internal shell structure. Shown are the granular prismatic (PR) layer, the crossed-lamellar (CL) layer and the complex crossed-lamellar (CCL) layer.
Genes 11 00298 g001
Genes 11 00298 g002 550
Figure 2. Map showing nine M. arenaria sampling locations across Europe: Comacchio, Italy (ITA); Lisbon, Portugal (LIS); Brest, France (BRE); Plymouth, UK (PLY); Saint Andrews, UK (SAN); Le Crotoy, France (LCT); Balgzand, Netherlands (TEX); Sylt, Germany (SYL); Kiel, Germany (KIE)
Figure 2. Map showing nine M. arenaria sampling locations across Europe: Comacchio, Italy (ITA); Lisbon, Portugal (LIS); Brest, France (BRE); Plymouth, UK (PLY); Saint Andrews, UK (SAN); Le Crotoy, France (LCT); Balgzand, Netherlands (TEX); Sylt, Germany (SYL); Kiel, Germany (KIE)
Genes 11 00298 g002
Genes 11 00298 g003 550
Figure 3. Results of genetic analysis of 247 soft-shell clams genotyped at 11 microsatellites. (a) A scatterplot of individual variation in principal component (PC) scores derived from principal component analysis (PCA) of the microsatellite dataset. The amounts of variation explained by each PC are given as percentages and the colored ellipses represent 95% confidence intervals for each population; (b) Results of the Structure analysis showing mean and standard deviations of estimated Ln probabilities of the data [P(D)] (dark grey) and ΔK (light gray) for each value of the number of groups (K); Panels (c) and (d) show estimated group membership coefficients obtained from Structure analyses with the number of groups (K) set to three and four respectively. Each individual is represented by a vertical line partitioned into segments of different color, the lengths of which indicate the posterior probability of membership to each group. The populations in panel (a) have been color-coded according to the four group solution shown in panel (d).
Figure 3. Results of genetic analysis of 247 soft-shell clams genotyped at 11 microsatellites. (a) A scatterplot of individual variation in principal component (PC) scores derived from principal component analysis (PCA) of the microsatellite dataset. The amounts of variation explained by each PC are given as percentages and the colored ellipses represent 95% confidence intervals for each population; (b) Results of the Structure analysis showing mean and standard deviations of estimated Ln probabilities of the data [P(D)] (dark grey) and ΔK (light gray) for each value of the number of groups (K); Panels (c) and (d) show estimated group membership coefficients obtained from Structure analyses with the number of groups (K) set to three and four respectively. Each individual is represented by a vertical line partitioned into segments of different color, the lengths of which indicate the posterior probability of membership to each group. The populations in panel (a) have been color-coded according to the four group solution shown in panel (d).
Genes 11 00298 g003
Genes 11 00298 g004 550
Figure 4. Scatterplot of individual variation in the first two principal components (PCs) from a PCA performed on elliptic Fourier analysis coefficients of lateral and anterior left shell views. The amounts of variation explained by each PC are given as percentages and the ellipses represent 95% confidence intervals for each population. The ellipses are color coded according to the main genetic groups discovered by Structure (shown in Figure 3d). Extreme and average reconstructed shell outlines are shown in grey.
Figure 4. Scatterplot of individual variation in the first two principal components (PCs) from a PCA performed on elliptic Fourier analysis coefficients of lateral and anterior left shell views. The amounts of variation explained by each PC are given as percentages and the ellipses represent 95% confidence intervals for each population. The ellipses are color coded according to the main genetic groups discovered by Structure (shown in Figure 3d). Extreme and average reconstructed shell outlines are shown in grey.
Genes 11 00298 g004
Genes 11 00298 g005 550
Figure 5. Summary of models of shell shape, thickness and organic content. Panels (a) and (b) show the mean effect sizes and bootstrapped 95% confidence intervals (CIs) of predictor variables estimated from generalized linear mixed models (GLMMs) of shell shape and shell thickness respectively. Results are summarized in panel (a) separately for each of three shape principal components (PCs), which are respectively color-coded in red, yellow and blue respectively. Panel (b) summarizes the results of the whole shell thickness model (in black) and the shell layers model, with the granular prismatic (PR) layer shown in red, the crossed-lamellar (CL) layer shown in yellow and the complex crossed-lamellar (CCL) layer shown in blue. Regression parameters were considered statistically significant when the bootstrapped 95% CI (error bars) did not overlap zero (asterisks denote significant differences from zero). (c) Relationship between latitude and the thickness of the PR, CL and CCL layers. Mean values (solid lines) and confidence intervals (shaded areas) were predicted while controlling for shell length. (d) Latitudinal differences in shell organic content of the PR (red bars) and CCL (blue bars) layers between representative warm temperate (LIS) and cold temperate (SAN) populations. Error bars represent 95% CIs, asterisks represent statistically significant comparisons (p < 0.05) and NS denotes non-significant comparisons.
Figure 5. Summary of models of shell shape, thickness and organic content. Panels (a) and (b) show the mean effect sizes and bootstrapped 95% confidence intervals (CIs) of predictor variables estimated from generalized linear mixed models (GLMMs) of shell shape and shell thickness respectively. Results are summarized in panel (a) separately for each of three shape principal components (PCs), which are respectively color-coded in red, yellow and blue respectively. Panel (b) summarizes the results of the whole shell thickness model (in black) and the shell layers model, with the granular prismatic (PR) layer shown in red, the crossed-lamellar (CL) layer shown in yellow and the complex crossed-lamellar (CCL) layer shown in blue. Regression parameters were considered statistically significant when the bootstrapped 95% CI (error bars) did not overlap zero (asterisks denote significant differences from zero). (c) Relationship between latitude and the thickness of the PR, CL and CCL layers. Mean values (solid lines) and confidence intervals (shaded areas) were predicted while controlling for shell length. (d) Latitudinal differences in shell organic content of the PR (red bars) and CCL (blue bars) layers between representative warm temperate (LIS) and cold temperate (SAN) populations. Error bars represent 95% CIs, asterisks represent statistically significant comparisons (p < 0.05) and NS denotes non-significant comparisons.
Genes 11 00298 g005
Table
Table 1. Table of sampling locations, including the number of samples used for the genetic and morphometric analyses. Four genetic diversity statistics are also shown. Observed heterozygosity (Ho), expected heterozygosity (He), the number of alleles (Na) and allelic richness (Ar) are given as values averaged across loci, with standard deviations reported in parentheses.
Table 1. Table of sampling locations, including the number of samples used for the genetic and morphometric analyses. Four genetic diversity statistics are also shown. Observed heterozygosity (Ho), expected heterozygosity (He), the number of alleles (Na) and allelic richness (Ar) are given as values averaged across loci, with standard deviations reported in parentheses.
Population IDLocationSamples Used for Genetic AnalysisSamples Used for Shape AnalysisSamples Used for Thickness AnalysisHoHeNaAr
ITAComacchio (Italy)3028180.75 (0.13)0.76 (0.07)8.23 (2.37)6.24 (1.6)
LISLisbon (Portugal)2829210.67 (0.19)0.74 (0.18)8.36 (2.24)6.56 (1.73)
BREBrest (France)3235190.76 (0.19)0.77 (0.12)9.54 (3.32)7.08 (2.56)
PLYPlymouth (UK)2730190.66 (0.13)0.75 (0.09)8.63 (3.44)6.6 (2.24)
SANSaint Andrews (UK)2322200.74 (0.12)0.72 (0.1)5.91 (1.64)5.16 (1.45)
LCRLe Crotoy (France)3140210.71 (0.14)0.75 (0.08)8.18 (2.44)6.36 (1.54)
TXLBalgzand (Netherlands)3030210.71 (0.15)0.76 (0.1)8.72 (2.76)6.36 (1.65)
SYLSylt (Germany)1213100.62 (0.14)0.73 (0.13)6.36 (2.65)6.63 (2.65)
KIEKiel (Germany)3435180.59 (0.15)0.74 (0.11)8.91 (2.38)6.29 (1.52)
Total 2472621670.69 (0.06)0.75 (0.02)8.09 (1.19)6.36 (0.52)
Table
Table 2. Summary of the results of GLMMs of shell shape and shell layer thickness. Estimated statistics and bootstrapped 95% confidence intervals (CIs) for regression parameters are reported for the modelled relationships described in Equations (1) and (2) in the Materials and methods. Because both shell shape and layer were analyzed as categorical variables, shape PC1 and the PR layer were used as the reference levels in the respective models. Regression parameters were considered statistically significant and highlighted in bold when the bootstrapped 95% CI did not overlap zero.
Table 2. Summary of the results of GLMMs of shell shape and shell layer thickness. Estimated statistics and bootstrapped 95% confidence intervals (CIs) for regression parameters are reported for the modelled relationships described in Equations (1) and (2) in the Materials and methods. Because both shell shape and layer were analyzed as categorical variables, shape PC1 and the PR layer were used as the reference levels in the respective models. Regression parameters were considered statistically significant and highlighted in bold when the bootstrapped 95% CI did not overlap zero.
CoefficientEstimateSE95% CItp-Value
Shell shape GLMM †
(Intercept)0.0020.11−0.37; 0.370.010.99
Shape (PC2)−0.010.17−0.52; 0.49−0.070.94
Shape (PC3)0.020.19−0.51; 0.530.120.91
Latitude−0.0040.10−0.23; 0.22−0.030.97
Length × Shape (PC1)−0.110.08−0.35; 0.13−1.270.20
Length × Shape (PC2)−0.070.10−0.27; 0.13−0.700.49
Length × Shape (PC3)0.270.100.07; 0.452.700.0072
gPC10.030.06−0.09; 0.150.470.64
gPC2−0.010.06−0.13; 0.12−0.150.88
sMLH0.020.04−0.06; 0.110.550.58
Shell layers thickness GLMM *
(Intercept)274.607.15257.38; 291.4438.41<0.0001
Layer (CL)−181.297.30−204.95; −158.30−24.82<0.0001
Layer (CCL)−19.4514.63−45.03; 5.82−1.330.18
Latitude × Layer (PR)24.307.369.02; 39.273.300.0010
Latitude × Layer (CL)−2.683.40−17.10; 11.27−0.790.43
Latitude × Layer (CCL)−23.7213.36−37.43; −9.69−1.780.076
Length × Layer (PR)25.366.9212.29; 38.933.660.0003
Length × Layer (CL)6.742.67−7.63; 21.442.530.012
Length × Layer (CCL)30.719.7914.88; 46.793.140.0018
gPC1−0.982.65−9.80; 7.82−0.370.71
gPC2−1.622.45−10.62; 7.44−0.660.51
sMLH0.822.14−6.80; 8.640.380.70
† The random intercepts for the shell shape PCs were normally distributed with mean 0 and variances 0.23, 0.41 and 0.46 (for PC1, PC2 and PC3), respectively. * The random intercepts for the PR, CL and CCL layers were normally distributed with mean 0 and variances 5.78, 3.42 and 30.41, respectively. The variance structure indicates different standard deviations per layer (PR: 1.00; CL: 0.32; CCL: 1.15) and an exponential of the variance covariate gPC2 structure.

Share and Cite

MDPI and ACS Style

De Noia, M.; Telesca, L.; Vendrami, D.L.J.; Gokalp, H.K.; Charrier, G.; Harper, E.M.; Hoffman, J.I. Population Genetic Structure Is Unrelated to Shell Shape, Thickness and Organic Content in European Populations of the Soft-Shell Clam Mya Arenaria. Genes 2020, 11, 298. https://doi.org/10.3390/genes11030298

AMA Style

De Noia M, Telesca L, Vendrami DLJ, Gokalp HK, Charrier G, Harper EM, Hoffman JI. Population Genetic Structure Is Unrelated to Shell Shape, Thickness and Organic Content in European Populations of the Soft-Shell Clam Mya Arenaria. Genes. 2020; 11(3):298. https://doi.org/10.3390/genes11030298

Chicago/Turabian Style

De Noia, Michele, Luca Telesca, David L. J. Vendrami, Hatice K. Gokalp, Grégory Charrier, Elizabeth M. Harper, and Joseph I. Hoffman. 2020. "Population Genetic Structure Is Unrelated to Shell Shape, Thickness and Organic Content in European Populations of the Soft-Shell Clam Mya Arenaria" Genes 11, no. 3: 298. https://doi.org/10.3390/genes11030298

APA Style

De Noia, M., Telesca, L., Vendrami, D. L. J., Gokalp, H. K., Charrier, G., Harper, E. M., & Hoffman, J. I. (2020). Population Genetic Structure Is Unrelated to Shell Shape, Thickness and Organic Content in European Populations of the Soft-Shell Clam Mya Arenaria. Genes, 11(3), 298. https://doi.org/10.3390/genes11030298

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop