Abstract

Late Quaternary (Eemian) deposits of the Netherlands contain shells that resemble those of living Mytilus galloprovincialis. Similar broad-shelled mytilids also occur in estuaries of the southwestern Netherlands together with slender individuals typical of M. edulis. We sampled living mussels along a depth gradient in the Oosterschelde to a) investigate whether a relation exists between shell shape and depth, b) test if the broad-shelled specimens might represent M. galloprovincialis (or a hybrid with M. edulis) and c) assess by inference if the Quaternary specimens might be attributed to M. galloprovincialis as well. In order to do so, we compared genetic (length polymorphism of Me 15/16, COIII sequences and AFLPs) and shell-morphological characteristics (juvenile L/W ratios and so-called Verduin parameters) of the same specimens. The obtained dataset indicates that all studied mussels from the Oosterschelde should be attributed to M. edulis, including those with broad shell outlines. No correlation of shell-morphology and depth-distribution was found. The worn and generally damaged state of the Eemian specimens precluded measurement of the Verduin parameters, while juvenile L/W ratios turned out not to be diagnostic. Therefore the shell characters examined in this study are insufficient to demonstrate the possible presence of M. galloprovincialis shells in Quaternary deposits of the Netherlands.

Introduction

The taxonomy of common European Mytilus species is complex. The wide variety of shell shapes as well as the common occurrence of genetic interchange between populations shows that the discrimination between species is not straightforward (Filipowicz et al., 2008; Śmietanka et al., 2009). Conventionally, three Mytilus species are distinguished within the European seas, viz. Mytilus edulisLinné, 1758, M. galloprovincialisLamarck, 1819 and M. trossulus Gould, 1850.

Mussels living along the coast of the southwestern Netherlands show a wide variety of shell shapes. On exposed North Sea localities, populations are entirely dominated by slender shells typically identified as Mytilus edulis. Broad specimens with a relatively sharp dorsal ridge and a pointed umbo are very rare on these exposed shores, although flat and broad specimens with a characteristic light brown periostracum do occur on floating material (Entrop, 1965) and are also known from wrecks on the North Sea floor (pers. obs. FW). The identity of these light brown mussels with a flat and broad shell shape remains to be established. Bluish broad-shelled Mytilus specimens are not uncommon in harbour basins and estuaries. Intermediate forms between broad and slender-shelled bluish mussels also occur (Fig. 1). Some of the broad forms resemble shells of Mytilus galloprovincialis, a species well known from the Mediterranean but also occurring in southwestern Europe (SanJuan et al., 1994, 1997). Similar broad-shelled mussels with a pointed umbo are also known from Eemian interglacial faunas of the Netherlands (Fig. 1). The closest populations of M. galloprovincialis, as confirmed by genetic analyses, are found in the western Channel region (Hilbish et al., 2002). In Dutch Mytilus populations, M. galloprovincialis-type alleles have been found at very low frequencies (Luttikhuizen et al., 2002; Kijewski et al., 2009). The centre of the modern distribution of M. galloprovincialis is to the south of The Netherlands. Since the Eemian interglacial fauna contains various taxa with a more southerly centre of distribution than the present Dutch fauna (e.g. Nucula nucleus (Linné, 1758), Lucinella divaricata (Linné, 1758), Solen marginatusPulteney, 1799, Pholas dactylusLinné, 1758, Gibbula magus (Linné, 1758)) it would be realistic to assume that the Eemian broad-shelled mussels might belong to M. galloprovincialis. However, the exact identity of the Eemian as well as the modern broad-shelled mussels in The Netherlands remains to be investigated. The generally damaged or worn preservation state of Eemian mussels precludes assessment of shell morphological characters including relative adductor scar lengths, as shells are incomplete (suppleted material) or adductor scars are worn (fossils washed ashore). Their identification can only be based on indirect comparison with living broad-shelled mussels.

[ *image not found: m8002a01:FIG1 ]

In this paper we investigate the identity of extant broad-shelled specimens of Mytilus from the southwestern Netherlands by comparing their shell morphological and genetic data and we indirectly evaluate whether broad-shelled Eemian mussels could be attributed to M. galloprovincialis.

Distribution of modern European mussel species

The three European Mytilus species have different, partially overlapping distribution ranges. Mytilus trossulus has a subarctic high-boreal distribution (also outside Europe). It is well adapted to cool temperatures and wide salinity variations. Mytilus edulis is recorded to be common along the entire European Atlantic coasts, southwards to Morocco (Jaziri et al., 2003), and also occurs in the western part of the Mediterranean (Quesada et al., 1995). Mytilus galloprovincialis is regarded as the dominating mussel species in the Mediterranean, but also lives along the Atlantic coast of Europe, as far north as northwestern Ireland (e.g. Gardner, 1992). The species appears the least adapted to salinity variations, and in general thrives in higher temperature settings than the other two species.

Hybridization occurs in virtually every known case where Mytilus taxa occur sympatrically, or where they are in geographic contact (Koehn, 1991). This resulted in M. edulis-galloprovincialis and M. edulis-trossulus hybrids in Europe (Hilbish et al., 2002; Rawson and Hilbish, 1998; Riginos et al., 2002; Kijewski et al., 2009) and M. trossulus-galloprovincialis hybrids outside Europe (e.g. Inoue et al., 1997; Rawson et al., 1999). The common occurrence of hybrids questions the validity of the taxonomic status of Mytilus species. In convention with common terminology, we use the before mentioned names, as if they were well-established species, despite these taxonomic uncertainties (see discussion).

Identification

The discrimination of Mytilus species based on shell characters is problematic. Historically, identification was based on the combination of geographic provenance and the overall shape of the shell. Not a single shell-character (Koehn, 1991; Gosling, 1992) nor a combination of shell characters (McDonald et al., 1991) discriminates completely between M. edulis and M. galloprovincialis. Transplantation experiments as well as field surveys have shown that shell outlines in Mytilus species are strongly influenced by environmental settings, such as the degree of exposure to waves (Seed, 1968; Akester and Martel, 2000), or population density and type of predators (Innes and Bates, 1999 and references therein). Multivariate approaches based on large character data sets, independently verified by molecular data, have yielded incomplete species discrimination (McDonald et al., 1991; Innes and Bates, 1999). However, a set of complex shell characters (including shapes and sizes of various adductor and retractor muscle scars) proposed by Verduin (1979) might discriminate between M. edulis and M. galloprovincialis, although this is in need of confirmation by molecular data. Furthermore, the Verduin study used few populations located far apart, and lacked areas in between from where we now know that hybrid populations exist.

The identity of broad-shelled mussels found in Late Quaternary deposits of the Netherlands and extant ones living in the southwestern Dutch coastal region still needs to be established. In this paper, we examine shell-characters (L/W ratios and parameters after Verduin, 1979) and genetic markers (Me15/16, COIII and AFLP) of slender- and broad-shelled specimens of Mytilus, from the Oosterschelde estuary in the southwestern Netherlands, in order to establish whether the broad-shelled mussels may include M. galloprovincialis or hybrids.

Material and methods

Sample collection

In total 229 specimens of Mytilus sp. were collected by SCUBA diving at nine different depths at the second northernmost pillar at the north side of the Zeeland Bridge (Oosterschelde, the Netherlands, 51°37’42”N, 3°54’49”E, Table 1). We sampled at different depths in order to establish the distribution of morphological variation or potential distribution of species. Using regular scissors, individual animals were removed from the hard substrate by cutting their byssus threads. Under water, the mussels were collected in separate, labelled bags. Subsequently, the mussels were transported to the Netherlands Institute of Ecology (NIOO, Yerseke) for registration and dissection. Using a scalpel and pincers, mantle and anterior muscle tissue were separated from the rest of the animal. Half of each animal was stored in absolute ethanol (from here on denoted by suffix A) and the other half (denoted by suffix B) was stored in a saturated CTAB – 20% DMSO solution. The corresponding shells were labelled, cleaned and stored for shape analysis. Tissue samples were numbered as to enable recognition of the various parts of each individual. Shells and tissue samples in ethanol (A) were transported to the Netherlands Centre for Biodiversity, Naturalis, The Netherlands, for the study of shell morphometrics (shells were kept as voucher specimens) as well as for genetic analysis of markers Me 15/16 and COIII. Tissues of the same specimens stored in CTAB (B) were sent to the Radboud University Nijmegen for AFLP analysis.

[ *image not found: m8002a01:TAB1 ]

Morphometrics

Two sets of shell parameters were measured, a length/width (L/W) set and a ‘Verduin’ parameter set. Juvenile shell outlines (shell L circa 15 mm) were drawn along growth lines preserved in adult mussel shells using a stereomicroscope equipped with a drawing attachment. From these outlines the maximum length and width (the latter perpendicular to length axis) were measured (Fig. 2A). These subadult outlines were used in order to minimise eco-phenotypic effects occurring in later growth stages (Seed, 1968). Based on the total distribution of L/W ratios, twenty-one specimens that represented slender and broad morphologies were selected for molecular analyses (described below). For these specimens, the length of the anterior adductor muscle scar and of the internal radius of the hinge plate (parameters utilized by Verduin, 1979, illustrated in Fig. 2B) were measured and standardized for the height (H) of the shell (in this case of the subadult outline).

[ *image not found: m8002a01:FIG2 ]

Sample treatment and DNA extraction

Total genomic DNA from tissue stored in ethanol (A) was extracted with a general CTAB protocol, using prolonged incubation and precipitation times. Approximately 3 mm3 of tissue was digested with 20 µl proteinase K (20 mg/ ml) overnight at 60°C in 500 µl CTAB-buffer. Extractions were done with 500 µl of chloroform : isoamylalcohol (24:1) and DNA was precipitated overnight at -20°C in 350 µl of isopropanol. The precipitate was washed with 500 µl ethanol/ammonium acetate, air-dried and resuspended in 50 µl MQ. DNA extractions on tissues stored in CTAB (B) were done with a DNeasy Tissue Kit (Qiagen), following the manufacturer’s protocol. Different extraction protocols were used (CTAB in Leiden, DNeasy Tissue Kit in Nijmegen), according to the expertise within the institutes at that time (2003). If the extraction method used would in any way affect the laboratory outcome, this would be equal for all the specimens analysed. Therefore, we deduce that the data from the separate analyses trajectories are compatible and that a comparison is legitimate.

Cytochrome c oxidase subunit III (COIII: mitochondrial marker)

In order to assess the recent history within species or populations of Mytilus, we sequenced a 416 bp fragment of Cytochrome c oxidase subunit III(COIII), a mitochondrial gene that shows high substitution rates. Mytilus, like other members of the Mytilidae and Unionidae (Curole and Kocher, 2002) transmits its mitochondria by means of doubly uniparental inheritance (DUI: Zouros et al., 1994a, 1994b, also known as gender associated inheritance: Skibinski et al., 1994a, 1994b). This results in two independently evolving mitochondrial lineages (termed F and M). Males possess both lineages, whereas females only possess the F-lineage. Our shell morphometric and molecular analyses were carried out without prior anatomical gender assessment. Although all the DNA extracts had to contain F-lineage mitochondria, M-lineage mitochondria were potentially present as well. At the time the lab work for this study was carried out (2003) this was a technical difficulty, because only generic or M-specific COIII primers (M. edulis) were available (nowadays both F- and M-specific primers are available; Śmietanka et al., 2009).

To remove potentially present M-lineage COIII products, we slightly modified the nested-PCR procedure described by Stewart et al. (1995). Initial amplification was done with primer For1 5’-TATGTACCAGGTCCAAGTCCGTG-3’ and Rev1 5’-ATGCTCTTCTTGAATATAAGCGTACC-3’ using the same reaction conditions as described for Me15/Me16, but with an annealing temp. of 54°C and a final extension at 70°C. Amplification products were split and half of each product was digested with Mbo I, while the other half was digested with Ssp I (New England Biolabs, following manufacturer’s protocol). Recognition sites for these restriction enzymes are located on the M-lineage COIII sequences of M. edulis and M. trossulus, respectively. In 2003, it was not known if M-mitotype COIII sequences of M. galloprovincialis possessed a recognition site for these enzymes. Digestion products were separated on a 1% agarose gel and the largest fragment of the MboI digest was extracted from the gel and purified with a Qiaquick Gel Extraction Kit (Qiagen), following manufacturer’s protocol. A 1:100 dilution (MQ) of the purified product was used as a template (1 µl) for a half nested PCR with primers For2 5’-GTAACTCAAGCCCATAAGAG-3’ and Rev1. Products of this second PCR were cleaned with a Qiaquick PCR purification kit (Qiagen) and sequenced (BigDye Terminator, Applied Biosystems) in both directions, using the same primers. Sequence products were purified using AutoSeq G-50 columns (Amersham) and run on an ABI 377 sequencer (Applied Biosystems). Forward and reverse sequences were assembled and edited with Sequencher version 4.2 (Gene Codes corp.) and the resulting contig sequences were aligned manually in MacClade version 4.08 (Maddison and Maddison, 2005). A minimum spanning network was created using SplitsTree4 (Huson and Bryant, 2006).

Statistical analysis

In order to see if the L/W ratios were distributed normally a Shapiro-Wilk test was done. To test the association between L/W ratios and depth, a univariate ANOVA (sum of squares type III) was performed in GLM (General Linear Model), in which the L/W ratios (continuous) were used as the dependent variable and depth as a covariate. A Fisher exact test was used to test for an association between L/W ratios (categorized as either slender, < 1.65 or broad, > 1.65) and COIII groups (this covariate was categorized as group A or B; see the results section). These statistical tests were performed in SPSS release 11.0.4 (SPSS, 2005).

To test for a relation between AFLP-data and L/W-ratios a Principal Component Analyis (PCA) and Correspondence Analysis (CA, Nenadic and Greenacre, 2007) were performed in R (R foundation for statistical computing, Ihaka and Gentleman, 1996). Using scatter plots, we inspected whether the two L/W-ratio categories had very different combinations of scores of the new independent variables obtained in either PCA or CA. Using logistic regression, we tested whether the scores of new variables significantly explained the probability of being either broad or not.

Results

Shell morphology

The L/W ratios of the 229 measured specimens (Fig. 3) range between 1.42 and 2.00 (mean 1.65, standard deviation 0.10). The L/W ratios are continuous with the exception of two specimens that have a slender subadult shell outline. Slender- and broad-shelled forms are connected through intermediates and the entire L/W ratio range describes an almost a perfect normal distribution (Fig. 3). After exclusion of the two just mentioned outliers, the Shapiro-Wilk test statistic becomes 0.994 and p = 0.430. No correlation between L/W ratios and sampling depth was observed (GLM univariate ANOVA, F = 1.026, p = 0.324); distribution of L/W ratios over different depths is illustrated by Fig. 4. The combined Verduin parameters of the studied Oosterschelde samples fall within the range of Mytilus edulis (Fig. 5, Table 2). The range of morphological variation of the 21 Oosterschelde shells covers a large part of the total morphological variation found by Verduin (1979) for various Western European Mytilus edulis populations. Mytilus shells from the Zeeland Bridge show considerable morphological variation in L/W ratios as well as Verduin parameters. The latter fall within the range of M. edulis (Fig. 5). Hence the Verduin parameters imply that the broad morphs should be attributed to Mytilus edulis.

[ *image not found: m8002a01:FIG3 ]

[ *image not found: m8002a01:FIG4 ]

[ *image not found: m8002a01:FIG5 ]

[ *image not found: m8002a01:TAB2 ]

Adhesive protein gene

The Me15/16-PCR products for the 21 selected specimens all had a length of 180 bp and belong to M. edulis, following Inoue et al. (1997).

COIII

Digestion of the initial PCR product (primers For1 and Rev1) with MboI resulted in one fragment if only product from the F-lineage was present (or if the MboI recognition site would have been lost in the M-lineage). If next to F-lineage product, M-lineage product was present as well, three fragments were observed. Digestion of initial PCR product with SspI resulted in two fragments if COIII was amplified from the F-lineage, and in four fragments if it was amplified from the M-lineage. Hence products from both lineages possess an SspI recognition site, whereas only the M-lineage possesses a restriction site for MboI. Subsequent nested PCR on the largest fragment of the MboI-digest with primers For2 and Rev1 resulted in a 462 bp product for all specimens. A minimum spanning haplotype network (Fig. 6) was created based on the nucleotide sequences from this product (primer sites were excluded, length: 416 bp). The average uncorrected P-distance between the obtained sequences was 0.012 (standard deviation 0.011). This relatively high average can mostly be ascribed to the large distances added by the specimens 203, 247 and 248 (indicated as group B). These three specimens were separated from the most common haplotype (largest circle, Fig. 6) by 14 mutational steps or more, whereas all other specimens (group A) were separated from this haplotype by 6 mutational steps or less. If the sequences from group B are excluded, the average uncorrected P-distance decreases to 0.007 (standard deviation 0.006). All substitutions were synonymous and the largest uncorrected P-distance was 0.036 (between 203 and 81, 162 and 222). No association between L/W ratio and COIII groups A and B could be shown (Fisher’s exact test, p (2-sided) = 1.000). Sequences were deposited in GenBank (accession numbers DQ353870-DQ353890).

[ *image not found: m8002a01:FIG6 ]

AFLP

Of the 248 putative gene loci that were scored, only 13 were non-monomorphic, indicative of a homogeneous population. Only the first four principal components each explained more than 10% of the total variation in the AFLP sample. Cumulatively these four principal components explained 72% of the total variance. In the correspondence analysis, only the first two principal inertias explain more than 10% of the summed inertias, cumulatively they contribute 61 %. Both the PCA and CA scatter plots (results not shown) suggested some separation in L/W groups based on the AFLP-data, but always with an intricate dependence on the scores. Backward model selection using likelihood ratio tests, showed that none of the principal components or none of the coordinates of the principal inertias significantly explained the variation in L/W ratios. Therefore we decided to depict the AFLP results by means of a minimal spanning network (Fig. 7, SplitsTree4, Huson and Bryant, 2006), which more clearly illustrates the lack of separation of L/W groups. Also Fig. 7 shows a single population in which slender-shelled specimens (L/W ratio > 1.65; represented by grey circles) and broad-shelled specimens (L/W ratio < 1.65; represented by white circles) are connected repeatedly by only a single step.

[ *image not found: m8002a01:FIG7 ]

Discussion

The combined Verduin parameters and molecular data (Table 3) show that the analysed Mytilus specimens from the Zeeland Bridge should be attributed to Mytilus edulis. Three specimens (38, 52 and 61) have Verduin parameter values that lie within the M. edulis range, but are close to the M. edulis - M. galloprovincialis border (Fig. 5). None of the other datasets (Me15/16 and the minimum spanning networks based on COIII sequences and AFLP-data) indicate that these specimens are any different from the remainder eighteen.

[ *image not found: m8002a01:TAB3 ]

As for the molecular datasets, Me15/16 unanimously shows the presence of M. edulis and does not indicate potential M. edulis – M. galloprovincialis hybrids, let alone ‘pure’ M. galloprovincialis. The AFLP data show a random distribution of morphotypes over the network (Fig. 7). Only the network based on COIII sequences (Fig. 6) points to the presence of more than a single group; group B contains the haplotypes obtained from specimens 203, 247 and 248 and group A contains nine other haplotypes. All haplotypes in group A are within six mutational steps from the main haplotype (Fig. 6). At least fourteen steps separate the haplotypes of group B from the main haplotype. A BLAST-search (http://blast.ncbi.nlm.nih.gov/Blast.cgi) showed that sequence DQ353885 (specimen 203, group B) was most similar to F-lineage COIII sequences of specimens identified as either M. galloprovincialis (from Greece and the Black Sea: accession numbers DQ403170, AY130171, AY130172, AY130174, DQ445474 and DQ445468) or M. edulis (from Ireland and France, accession numbers AY130154 and AY130161). Śmietanka et al. (2009) also obtained partial COIII sequences for Mytilus sp. with a reasonable coverage of the entire European coast, including the Westerschelde. Hence we expected that sequences from this study would be most similar to those of Śmietanka et al. This was not the case, due to the fact that the forward primer used by Śmietanka et al. (5’-TCTTGGTACAACTGCGGGAA-3’, Skibinski et al., 1994b) is located ca. 115 bp downstream from For2. A BLAST-search with only the overlapping part (301 bp) of sequences from this study and that of Śmietanka et al. (2009), indeed shows that sequence DQ353885 (specimen 203) is most similar to sequence FJ549913 from the Westerschelde. A minimum spanning network based on the 21 COIII sequences obtained with this study and the additional 26 COIII sequences of Śmietanka et al. (network not depicted) reveals seven additional haplotypes. Nevertheless, all of these haplotypes are within six mutational steps from either the main haplotype of group A or B. Based on these results we have no indication for the presence of more than two groups of haplotypes in the Dutch delta region. Since in our Oosterschelde material Me15/16 shows no indication for hybridisation and the AFLP-results do not indicate more than a single population, we assume incomplete lineage sorting (ancestral polymorphism) is a more plausible explanation for these COIII results, than introgression.

Despite the presence of two distinct groups (A and B, Fig. 6) of COIII haplotypes in our material, all molecular results indicate that the Oosterschelde mussels belong to a single population of Mytilus edulis. This also corresponds to our identifications based on the Verduin parameters. The existence of broad-shelled M. edulis and the lack of M. galloprovincialis in the Oosterschelde, makes it for now impossible to establish the identity of similar broad-shelled Mytilus shells from Eemian (Late Quaternary) deposits of the North Sea Basin. Moreover it shows that merely the relative width (L/W ratio) of broad-shelled mytilids, is not sufficient to discriminate between M. edulis and M. galloprovincialis. We call for further studies using molecular analyses to investigate the ability of Verduin parameters to discriminate between M. edulis and M. galloprovincialis shells.

Whether the taxonomic units, jointly known as the Mytilus species complex, are distinct enough to merit full specific status has been a matter of debate since the 1860s (Gosling, 1992). The choice of species concept is of importance in such discussions. Nowadays it becomes common practice to refer to ‘pure’ M. edulis and ‘pure’ M. galloprovincialis solely based on genetic characters (Kijewski et al., 2009) and species are ‘taxonomically identified’ by the length of PCR products (Śmietanka et al., 2009), regardless of the morphology (shell shape, colouration of periostracum and mantle, etc.) that was classically used to distinguish Mytilus species.

In the fossil shells, the Verduin parameters could not be established due to wear of critical characters (adductor scars). Furthermore, the attribution of extant mussels to M. galloprovincialis as proposed by Verduin (1979) still needs to be validated in populations by molecular data. All of the specimens from the Oosterschelde (this study), including broad-shelled ones (low L/W ratio), were attributed to M. edulis based on the Verduin parameters.

Conclusion

Based on an assessment and comparison between shell-morphological and molecular characteristics from single mussel specimens (Mytilus sp.) collected along a depth profile in the Oosterschelde we conclude that:

(1) broad-shelled specimens from this location (identified by low L/W ratios) do not belong to Mytilus galloprovincialis, but instead all of the analysed Oosterschelde specimens (both slender and broad-shelled) must be attributed to M. edulis.

(2) based on COIII sequence data, two groups of haplotypes can be discerned, whereas Me15/16 and AFLP data do not show any structure within this mussel population.

(3) there is no association between L/W ratio (broad- or slender-shelled specimens) and the depth at which specimens were collected.

(4) there is no association between the two COIII haplotype groups and shell morphological variation (COIII haplotype groups and L/W ratios are independent variables, whereas the Verduin parameters indicate only a single group: M. edulis).

(5) due to the existence of broad-shelled M. edulis (low L/W ratio), we cannot use the L/W ratio to classify similar broad-shelled specimens from Late Quaternary (Eemian) deposits of the North Sea Basin as either M. edulis or M. galloprovincialis.

Acknowledgements

We want to thank Dr. T. van Dooren, Prof. Dr. G. Vermeij and an anonymous reviewer for their suggestions to improve this manuscript. We are especially grateful to Dr. T. van Dooren for his help with the analyses in R. Furthermore we thank Dr. K. Vrieling for comments on an earlier version of this manuscript. Part of this study was supported by the European Committee (Research Directorate General, Environmental Program Marine Ecosystems) through BIOCOMBE project (contract EVK3-2001-00146).

Zouros E, Ball AO, Saavedra C, Freeman KR. 1994b. An unusual type of mitochondrial DNA inheritance in the blue mussel Mytilus. Proceedings of the National Academy of Sciences of the United States of America 91: 7463-7467.