INTRODUCTION
Surf clams (also known as duck clams or trough shells), belonging to the genus Mactra Linnaeus 1767, live in the surf zone of exposed beaches and are widely distributed along mud–sandy coasts of the Pacific Ocean, Atlantic Ocean, Black Sea and Mediterranean Sea (Conroy et al., Reference Conroy, Smith, Michael and Stotter1993). They represent commercially important bivalves in many countries and are extensively utilized as seafood, raw materials for manufacturing flavouring materials and live feed at various aquaculture farms (Hou et al., Reference Hou, Lü, Zou, Bi, Yan and He2006).
The rayed trough-shell Mactra corallina (=M. stultorum) Linnaeus 1758 inhabits sandy bottoms at depths between 5 and 30 m, and it is distributed along coasts of the Black Sea, Mediterranean Sea and the eastern Atlantic Ocean from Norway to Senegal. It is a medium sized marine bivalve with a very thin and delicate shell with concentric growth lines. This species is present with two different and well-defined morphotypes, which, although they live sympatrically, are generally classified as two different sub-species. These morphotypes are easily distinguishable by the colour of the shell: the white variant, named M. corallina corallina Linnaeus 1758, has a shell of a hyaline white with weak ivory radial bands, whereas M. corallina lignaria Monterosato 1878 shows brownish radiating bands (D'Angelo & Gargiulo, Reference D'Angelo, Gargiulo and Fabbri1987; Fischer et al., Reference Fischer, Schneider and Bauchot1987).
The correct specific name for the rayed trough-shell M. corallina is a longstanding issue for zoologists and malacologists. As reported in the Mediterranean marine molluscs checklist (Chiarelli, Reference Chiarelli1999), three species belonging to the genus Mactra are present: M. stultorum (=M. corallina) Linnaeus 1758, M. glauca Von Born 1778 and M. olorina Philippi 1846. Within M. corallina, two taxa, M. c. corallina and M. c. lignaria, are recognized.
Nevertheless, based on analyses of partial region of 18S rDNA by PCR-SSCP, Livi et al. (Reference Livi, Cordisco, Damiani, Romanelli and Crosetti2006) found preliminary genetic evidences that the traditional classification of M. c. corallina and M. c. lignaria as subspecies was in contrast with the high genetic distance observed between the two taxa. Besides, M. c. corallina formed a highly supported cluster with a further unknown genetic profile, giving evidence of a third taxon belonging to the M. corallina complex (Livi et al., Reference Livi, Cordisco, Damiani, Romanelli and Crosetti2006).
In his handbook Carta d'Identità delle Conchiglie del Mediterraneo Parenzan (Reference Parenzan1976) describes five distinct phenotypes ascribable to the genus Mactra. But actually the most plausible hypothesis is that M. corallina is a complex formed by two or more species (Livi et al., Reference Livi, Cordisco, Damiani, Romanelli and Crosetti2006).
The official Italian checklists of marine fauna (compiled in their latest version in 2006 and available at http://www.sibm.it/CHECKLIST/principalechecklistfauna.htm) refer to these clams as belonging to the single species M. stultorum whereas the FAO identification handbook of Mediterranean species (Fischer et al., Reference Fischer, Schneider and Bauchot1987) and Riedl (Reference Riedl and Muzzio1991) indicate M. corallina as the valid name for this species and M. stultorum as a synonym.
We decided to adopt the FAO specific designation and thus we refer to the white variant as M. c. corallina and to the brown habitus as M. c. lignaria as described in D'Angelo & Gargiulo (Reference D'Angelo, Gargiulo and Fabbri1987).
This work represents a first attempt to resolve the confused and contradictory systematics of bivalves belonging to M. corallina putative species by analysing molecular and morphological characters of the two morphotypes observed. Analysed samples were collected along the north Adriatic coasts of Cesenatico (Italy). In the present study we analysed molecular data obtained by four DNA markers: a nuclear ribosomal DNA subunit (18S) and the mitochondrial genes cytochrome oxidase I (COI), small (12S) and large (16S) ribosomal subunits, in order to provide a stable and robust phylogenetic estimate of the target. In addition, a morphological analysis was carried out on the basis of five parameters of the shell.
MATERIALS AND METHODS
Sampling and DNA extraction
Samples were collected in the north Adriatic Sea in front of Cesenatico (Italy) during a single diving in September 2006 and stored at –80°C. To avoid the problem of collecting paralogous mtDNAs, as found in doubly uniparental inheritance (DUI) bivalve species (see Passamonti & Ghiselli, Reference Passamonti and Ghiselli2009, and references therein, for a review on the issue), foot muscle tissue was dissected from each individual using a sterile cutter and stored in 80% ethanol at 4°C for the following DNA extraction. DUI has not been searched for in Mactra, because of the lack of specimens with fully developed gonads, but even if it would be present, foot muscle is expected to mostly carry mtDNA of maternal origin (Garrido-Ramos et al., Reference Garrido-Ramos, Stewart, Sutherland and Zouros1998). Total genomic DNA was prepared from 25 mg of muscle tissue according to the DNeasy Tissue Kit (Quiagen) protocol.
DNA amplification, cloning and sequencing
Sequences from partial 12S, 16S, 18S and COI were obtained. PCR amplifications were carried out in a 50 µl volume, as follows: 5 µl reaction buffer, 150 nmol MgCl2, 10 nmol each dNTP, 25 pmol each primer, 20 ng genomic DNA, 1.25 units of DNA polymerase (Invitrogen, Carlsbad, CA, USA), water up to 50 µl. Thermal cycling consisted of 35 cycles at 94°C for 60 seconds, the specific annealing temperature (48°C for 12S and 16S; 50°C for 18S and COI) for 60 seconds, and 72°C for 60 seconds. An initial denaturation step (94°C for 5 minutes) and a final extension holding (72°C for 7 minutes) were added to the first and last cycle, respectively. Primer pairs were SR-J14197 ÷ SR-N14745 for 12S (Simon et al., Reference Simon, Buckley, Frati, Stewart and Beckenbach2006), 16SbrH(32) ÷ 16Sar(34) (5′–CGCCTGTTTAACAAAAACAT–3′) for 16S (modified from Palumbi et al., Reference Palumbi, Martin, Romano, McMillan, Stice and Grabowski1996), 18SF ÷ 18SR for 18S (Livi et al., Reference Livi, Cordisco, Damiani, Romanelli and Crosetti2006), and LCO1490 ÷ HCO2198 (Folmer et al., Reference Folmer, Black, Hoeh, Lutz and Vrijenhoek1994) for COI. Amplified DNAs were treated with Wizard® SV Gel and PCR Clean-Up System (Promega). For a single Mactra corallina lignaria individual it was necessary to clone the 18S rDNA gene fragment with Ultramax DH5α–Competent Cells (Invitrogen) following the manufacturer's instructions.
Purified amplifications were either cycle sequenced using the ABIPrism BigDye Terminator Cycle Sequencing kit (Applied Biosystems) and run on an ABI310 Genetic Analyser (Applied Biosystems) or sent to Macrogen (Seoul, EE Korea) for sequencing. Polymorphisms were confirmed by sequencing both strands.
Sequence analysis
Haplotypes (GenBank Accession Numbers FJ830395–FJ830446; Appendix 1) were aligned using the MAFFT multiple sequence alignment tool (Katoh et al., Reference Katoh, Misawa, Kuma and Miyata2002) available online at http://align.bmr.kyushu-u.ac.jp/mafft/online/server. Q-INS-i (Katoh & Toh, Reference Katoh and Toh2008) and G-INS-i (Katoh et al., Reference Katoh, Kuma, Toh and Miyata2005) algorithms were chosen for ribosomal- and protein-coding genes, respectively. Sequences of species belonging to different families of heterodont bivalves were downloaded from the NCBI databank and added to alignment as reference data. In order to compare orthologous characters, only female mtDNA sequences from GenBank were used for DUI species. Gaps were coded as presence/absence data following the simple indel coding method of Simmons & Ochoterena (Reference Simmons and Ochoterena2000) with the software GapCoder (Young & Healy, Reference Young and Healy2003).
The analysis of molecular variance (AMOVA) framework (Excoffier et al., Reference Excoffier, Smouse and Quattro1992) implemented in Arlequin v3.11 software (Excoffier et al., Reference Excoffier, Laval and Schneider2005) was used to test the overall genetic heterogeneity of surf clam samples. In this statistical method, a hierarchical AMOVA was carried out on the partitioning of molecular variability at arbitrarily chosen levels (i.e. from the individual to the group of samples level). In the present analysis, groups were obtained by pooling bivalve samples in two groups corresponding to the two subspecies Mactra corallina corallina and M. c. lignaria. Kimura 2-parameters distances (K-2-P; Kimura, Reference Kimura1980) were computed with MEGA4 software (Tamura et al., Reference Tamura, Dudley, Nei and Kumar2007) with pairwise deletion of gaps/missing data and with a uniform mutation rate. Φst and Fst fixation indices (mitochondrial and nuclear genome respectively) as implemented in Arlequin were calculated to assess the genetic divergence. Statistical significance was estimated by comparing the observed distribution with a null distribution generated by 1000 permutations, in which individuals were randomly re-distributed into samples.
A barcoding-like approach was used to analyse genetic distances computed as formerly described. Frequencies of intra- and inter- specific distances were separately plotted in histograms to provide a visual output of genetic differentiation between the two morphs.
Phylogenetic relationships were inferred through Bayesian analyses implemented in MrBayes 3.1.2 (Huelsenbeck & Ronquist, Reference Huelsenbeck and Ronquist2001; Ronquist & Huelsenbeck, Reference Ronquist and Huelsenbeck2003). All analyses employed a cold chain and three incrementally heated chains. Starting trees for each chain were randomly chosen and the default values of MrBayes were used for all settings (including prior distributions). Each metropolis coupled Markov Chain Monte Carlo (MCMC) was run for ten million generations, with trees sampled every 100 generations. Burn-in was visually determined for each gene fragment by plotting average standard deviation of split frequencies over generation seeking for apparent convergence. Chains had always converged to a stable average standard deviation of split frequencies values <0.01.
Posterior probabilities (PP) were used to assess clade support. Analyses were performed using the evolutionary models selected for each gene fragment by the Bayesian information criterion of Modeltest (Posada & Crandall, Reference Posada and Crandall1998). Selected models were K81uf+Γ (Kimura, Reference Kimura1981) for 12S and 16S, K80+Γ (Kimura, Reference Kimura1980) for 18S, and TVM+Γ for COI. They were implemented into MrBayes with the more similar and more complex model available in the program. Mytilus galloprovincialis (female) was used as outgroup to root phylogenetic trees. Nodes with PP < 0.95 were collapsed with the exception of 12S gene fragment data (PP < 0.85). Trees were graphically edited by MrEnt v2.0 (Zuccon & Zuccon, Reference Zuccon and Zuccon2006).
Morphological analysis
Five morphological variables were measured: (i) shell length (antero-posterior, L); (ii) height of specimens (ventro-dorsal, H); (iii) maximum width of shell (left–right, W); (iv) distance between the points of intersection of the adductor muscles impressions and the pallial line (AP); and (v) distances between the points of intersection of the adductor muscles impressions and the apex of the umbo (UA and UP). Parameters were measured to 0.01 cm with a caliper. On the basis of such measures, the ratios H/L, W/L and W/H were obtained. Plots were graphically edited by R (Ihaka & Gentleman, Reference Ihaka and Gentleman1996). Morphological data were statistically treated with Pearson's coefficient (r) to assess correlation between different sizes; ratios were examined by analysis of F test and the Welch two samples t-test to assess mean differences. The F test is a statistic used to test the hypothesis that two parameters have the same variance against the alternative hypothesis that the variances are different. Degrees of freedom were calculated taking into account number of groups (i.e. gl1 = 2−1 = 1) and number of specimens (i.e. gl2 = [35– 1]+[28 – 1] = 61). The critical values of F with P = 0.975 were calculated with the function qf(p, gl1, gl2) as implemented in R statistical computing software (Ihaka & Gentleman, Reference Ihaka and Gentleman1996; R Development Core Team, 2009). Welch's t-test is an adaptation of the Student's t-test intended for use with two samples having possibly unequal variances. Values of t-test were calculated using the function t.test(x1, x2) implemented in R software.
RESULTS
Genetic data
Twenty individuals for each morphotype were collected. A total of 34 specimens, 15 ascribed to Mactra corallina corallina and 19 to M. c. lignaria, were amplified and sequenced for the 12S, 16S, 18S and COI genes (partial sequences). Fragments of 397, 513, 516 and 571 bp respectively were obtained. Variable sites (including maximum parsimony informative sites), haplotype frequencies, specimen numbers and GenBank IDs are given in Appendix 1.
Data obtained by aligning the 12S partial sequence appeared quite soon less powerful than other gene fragments probably because of sampling artefacts. Actually, technical problems occurred during amplification and sequencing of the 12S and only four individuals of each group gave suitable PCR amplicons and electropherograms. Twenty-six repeated null amplifications were observed (11 in M. c. corallina and 15 in M. c. lignaria), accounting for the presence of point mutations in the annealing site of either primer. Further analyses will be required to unravel this latest issue.
In any case, examining sequence alignments for all the analysed gene fragments, high genetic divergences were observed between specimens of the two different morphs here considered (i.e. var. corallina and var. lignaria). Diagnostic sites were 7 out of 397 for 12S, 8 out of 513 for 16S, 25 out of 516 for 18S, and 43 out of 571 for COI (Appendix 1).
No mutation was observed at the amino acid level for the COI gene. Most point-mutations occurred at the third position of the codon. Six out of 60, however, were found at 2nd position (343, 358, 370, 412, 475 and 478).
Levels of genetic variability within the same morphotype were remarkably low and some shared haplotypes were observed (Appendix 1). A weak polymorphism was observed in the 18S fragment within both morphotypes, in the proportion of one different haplotype out of eleven in M. c. lignaria (sample n. 14 C2; C/T transition in position 467) and one out of six in M. c. corallina (sample n. 32; C/A, A/G, C/A transversion/transition in position 198, 200 and 202 respectively). Incidentally, the M. c. lignaria observed single 18S variant was found in a cloned sequence (see Appendix 1).
The higher proportion of overall molecular variance was always found at ‘between morphotypes’ hierarchical level (from 77.78%, P < 0.05; to 99.23%, P < 0.01; Table 1). All fixation indices values were high and significant or even highly significant. With the only exception of the 12S fragment (Fst = 0.778, P = 0.025), fixation indices values were higher than 0.90 and ranged from 0.902 (COI) to 0.992 (18S; Table 1).
Table 1. Analysis of partition of molecular variance (AMOVA) and fixation indices values (Fst for diploid data, Φst for haploid data). *, P = 0.05; **, P = 0.01; ***, P = 0.001.
Figure 1 shows histograms obtained by plotting intra- and inter- specific K-2-P distances for the four analysed gene fragments. Intra- and inter- morphotype distances are well separated and the gap between these distances ranges from about 0.005 (16S) to about 0.064 (COI), respectively.
Fig. 1. Histogram illustrating K-2-P distances distribution among Mactra corallina/M. lignaria group, as resulting from the four characterized genes. K-2-P distance values are reported on x-axis, whereas their frequencies are reported on y-axis. A, 12S; B, 16S; C, 18S; D, COI; light grey: intra-specific distances; dark grey: inter-specific distances.
The Bayesian analysis performed with different combinations of data yielded differently resolved but comparable and well supported tree topologies (Figures 2–5). In all trees, the two morphotypes clustered separately from all other sequence data with 0.95 < PP < 1.00. Mactra c. corallina was resolved as a monophyletic group for 12S (PP = 0.88), 18S (PP = 1.00) and COI (PP = 1.00). Similarly, M. c. lignaria was resolved as monophyletic for 16S (PP = 0.96), 18S (PP = 1.00) and COI (PP = 1.00). Both morphotypes were paraphyletic in other cases (i.e. 16S and 12S respectively). At a higher taxonomic level, the superfamily Mactroidea (= Mactracea) Lamarck 1809 (Mactridae Lamarck 1809+Anatinellidae Gray, 1853+Cardiliidae Fischer, 1887+Mesodesmatidae Gray 1840) appear to be monophyletic in all obtained trees, with PP values ranging from 0.97 to 1.00, while the superfamily Veneroidea Rafinesque 1815 showed a complex situation that would require further investigations.
Fig. 2. Bayesian phylogeny of Mactra corallina/M. lignaria samples inferred by 12S sequence data. Individuals belonging to the corallina morphotype are marked with a square whereas individuals belonging to the lignaria morphotype are marked with a triangle. For correspondences to the GenBank accession number, see Appendix 1.
Fig. 3. Bayesian phylogeny of Mactra corallina/M. lignaria samples inferred by 16S sequence data. Taxon symbols as in Figure 2.
Fig. 4. (A) Bayesian phylogeny of Mactra corallina/M. lignaria samples inferred by 18S sequence data. Taxon symbols as in Figure 2. Grey arrow heads point to Mesodesmatidae species; (B) histogram illustrating intergeneric K-2-P distances distribution among Mactridae: K-2-P distance values are reported on x-axis, whereas their frequencies are reported on y-axis; data from established genera of Mactridae are shown in white, whereas data from inter-specific comparisons among Mactra corallina/M. lignaria group are shown in dark grey, as in Figure 1C.
Fig. 5. Bayesian phylogeny of Mactra corallina/M. lignaria samples inferred by COI sequence data. Taxon symbols as in Figure 2.
Morphological data
Morphological analyses showed that only three parameters (i.e. L, H and W) were statistically significant, while AP, UA and UP did not present any element of significance on discriminating the two morphotypes (data not shown). As a consequence, the last three parameters were not considered and here we will take into account ratios that only involve the former three parameters.
The analysis of Pearson's correlation reflects the degree to which two variables are related. The correlation between the considered sizes gives the following r values: in M. c. corallina rH/L = 0.915, rW/L = 0.741 and rW/H = 0.749; in M. c. lignaria rH/L = 0.941, rW/L = 0.781 and rW/H = 0.777.
Both in M. c. corallina and M. c. lignaria, all morphological features considered were positively correlated. In particular, high values of r were found for correlation between H and L. Morphometric ratios found are given in Figure 6.
Fig. 6. Morphometric ratios in Mactra corallina and M. lignaria.
The F test applied to W/L and W/H ratios showed statistically significant values, while for H/L the null hypothesis cannot be rejected (Table 2). Similarly, the t-test assessed a significant difference in W/H and W/L ratios. No significant difference was found in H/L ratio (Table 2).
Table 2. Analysis of F test with P = 0.975 calculated with the function qf(p, gl1, gl2) (degrees of freedom: gl1 = 1 and gl2 = 61) and of the Welch two samples t-test calculated using the function t.test(x1, x2) applied to H/L, W/L and W/H ratios.
DISCUSSION
The development of molecular tools for species identification scored an increased importance because of difficulties of discriminating them on the basis of morphological characters only. This is mostly true for organisms at early developmental stages and in cases of morphological stasis of adults or presence of sibling species (Øines & Heuch, Reference Øines and Heuch2005; Livi et al., Reference Livi, Cordisco, Damiani, Romanelli and Crosetti2006).
Molecular assays presented in this paper brought to light a stable genetic divergence between M. c. corallina and M. c. lignaria. The clams analysed in this work were caught during a single dive in the very same area. The sympatric occurrence of the two morphotypes, coupled with the genetic divergence detected, is strong evidence of separate gene pools, thus supporting a reproductive isolation between the two morphs. Therefore, the taxon previously described as M. corallina should be rather considered as two different biological species, M. corallina and M. lignaria. A very similar experimental procedure, although based on allozyme analysis, was reported in Backeljau et al. (Reference Backeljau, Bouchet, Gofas and de Bruyn1994), who identify Chamelea gallina and C. striatula, previously considered as two subspecies of C. gallina, as two distinct and reproductively isolated biological species; actually, despite the probable overlap in breeding season between the two Chamelea morphotypes, they maintained a large genetic distance in sympatric conditions, giving evidence of two different biological species (Backeljau et al., Reference Backeljau, Bouchet, Gofas and de Bruyn1994).
For our Mactra, more genetic data obtained are consistent with two different species: the magnitude of genetic distances observed between M. c. corallina and M. c. lignaria were comparable to, if not greater than, distances detected among different genera belonging to the family Mactridae (K-2-P distance values, Figures 1 & 4B). The intra-specific pairwise K-2-P genetic distances were an order of magnitude lower than inter-specific comparisons (Figure 1). This divergence is also clearly shown by the high and statistically supported values of fixation indices, which were close to one and indicated the presence of a sharp dichotomy between genotypes, and the unbalanced partition of molecular variance with the majority of percentage detected at the higher hierarchical level, i.e. ‘among morphotypes’. In the phylogenetic trees, albeit in two cases a soft paraphyly was observed (Figures 2 & 3) we observed a separation of M. c. corallina clusters from M. c. lignaria clusters, supported by robust node values.
Finally, the observed variability in the 18S gene well falls within the range of expected variability for this locus. This gene, generally highly conserved within species, shows variability higher in bivalves than in other taxa (Adamkewicz et al., Reference Adamkewicz, Harasewych, Blake, Saudek and Bult1997; Passamaneck et al., Reference Passamaneck, Schandler and Halanych2004). Moreover, the unique different haplotype found in M. c. lignaria was collected from a clone, which might have brought to light a rare variant (i.e. intra-individual variability among 18S repeats within the nuclear genome).
Preliminary morphological analyses seem also concordant with genetic data, although only one shell character (other than the colour) was significantly different; in fact, the main morphological character discriminating the two morphs seems to be the W value (maximum width of shell, i.e. the convexity) which differentiates morphometrical ratios in specimens with the same length or height. According to the data, the ratios W/L and W/H assume a clear (and classic) diagnostic value and allows us to take the following value to discriminate the two groups: in M. c. corallina W/L > 0.50 and W/H > 0.60, while in M. c. lignaria W/L < 0.50 and W/H < 0.60.
The effective reproductive isolation between M. c. corallina and M. c. lignaria (and/or sterility of hybrids) has still to be directly demonstrated, but obtained data are sound enough to support the species level for both morphs. Nevertheless, an additional sampling along the Adriatic coasts has already been planned to better describe the genetic landscape of Mactra, which seems to represent a complex of at least two (but probably more) different species (Livi et al., Reference Livi, Cordisco, Damiani, Romanelli and Crosetti2006).
Finally, the phylogenetic position of Mactra was addressed in this study. On the basis of 18S and 28S rRNA genes, it was previously found that the superfamily MACTROIDEA, traditionally classified near to the superfamily CARDIOIDEA (=CARDIACEA) Lamarck 1809 with an implicit sister-group relationship, showed grater affinity to UNGULINIDAE H. & A. Adams 1857 and the group of VENERIDAE Rafinesque 1815—CORBICULARIDAE Gray 1847—ARCTIDAE Newton 1891—CHAMIDAE Blainville 1825, but no connection with CARDIOIDEA (Taylor et al., Reference Taylor, Williams, Glover and Dyal2007). In our preliminary phylogenetic analysis, the genus Mactra was always monophyletic, although the 16S sequence of Coelomactra antiquata obtained from GenBank generates a polyphyly in the clade of Mactra (polyphyly supported by a significant PP nodal value of 0.98). Moreover, the superfamily MACTROIDEA clustered separately in all trees and was statistically well supported. Finally, in the 18S tree, individuals belonging to families MACTRIDAE and MESODESMATIDAE were intermingled (Figure 4A). This situation suggests further investigation focused on these species to assess the monophyly of the genus Mactra and to validate the family status of MESODESMATIDAE.
ACKNOWLEDGEMENTS
Two anonymous referees provided valuable comments on the manuscript. Dr Federica Rongai of the Laboratory of Applied Biotechnology for Aquaculture and Fishery of the Aquaculture Institute of Cesenatico (Italy) and Dr Davide Gambarotto of the Laboratory of Molecular Zoology of the Department of Experimental Evolutionary Biology of Bologna (Italy) supported the laboratory work.
Appendix 1. Alignment of the two variants of Mactra corallin a analysed (lig: lignaria, cor: corallina), related frequencies (f), specimen numbers as in figures 2 to 5 and GenBank accession number. Only variable sites are reported.
