Marine bivalves voucher DNA barcoding from Eastern Mediterranean, with evidence for Ostrea stentina invasion

Abstract

Bivalve molluscs are a diverse group of animals with particular economic and ecological importance. Their morphological characteristics frequently confuse their identification leading to mislabelling of edible species. Genetic diversity is critical to the resilience of marine bivalve populations in the face of environmental stressors such as ocean acidification and warming. In this study, we characterized the phylogeny and defined the first DNA barcodes of six marine bivalves [Ostrea edulis (Linnaeus, 1758) Arca noae (Linnaeus, 1758), Pinctada radiata (Leach, 1814), Venus verrucosa (Linnaeus, 1758), Calllista chione (Linnaeus, 1758) and Ruditapes decussatus (Linnaeus, 1758)] sampled from different coastal areas of Aegean and Ionian Seas using the molecular markers cytochrome c oxidase subunit I (COI) and 18S ribosomal RNA (18S rRNA). Further, COI gene was employed to investigate the population genetic diversity since 18S rRNA exhibited no conspecific differences. The sequence of 18S rRNA successfully discriminated the bivalves at family or superfamily level but occasionally proved insufficient for species identification. Contrariwise, COI was highly informative and could reliably distinguish all species. Population haplotype diversity was moderate to high and was always accompanied by generally low nucleotide diversity, indicating genetically closely related haplotypes. The invasive Pinctada radiata was found to be panmictic even among distant sampling areas, while Ostrea edulis was the only species that exhibited moderate levels of population subdivision. Finally, here we report for the first time the presence of Ostrea stentina in Thermaikos Gulf sampled among Ostrea edulis specimens, demonstrating a new invasive bivalve species in Eastern Mediterranean.

Introduction

DNA barcoding is a molecular technique that is based on species-specific nucleotide variation in standardized highly conserved genomic regions. This method was initially proposed as a universal tool for the identification of taxonomic biodiversity (Hebert et al., 2003). Nowadays, DNA barcoding is also utilized for environmental monitoring, exhibiting multiple applications in food authentication and traceability (Galimberti et al., 2013) in pharmaceutical industry, as well as in forensics (Shadrin, 2021). DNA barcodes constitute an effective tool for studying aquatic biodiversity (Radulovici et al., 2010; Bucklin et al., 2011; Tsoupas et al., 2022) offering the great advantage of rapid acquisition of molecular data compared to morphology-based species identification (Frézal and Leblois, 2008). Mitochondrial DNA markers are usually targeted for population genetics and phylogenetics since mtDNA exhibits features such as absence of introns, relative lack of recombination, increased mutation rate compared to nuclear DNA, neutral selection pressure as well as maternal inheritance in the most marine animals (Ladoukakis and Zouros, 2017). The universal primers for the fragment of the mitochondrial gene encoding for cytochrome oxidase subunit I (COI) are very robust allowing recovery of its 5′ end for most animal phyla (Zhang and Hewitt, 1997). Furthermore, COI appears to possess a broader range of phylogenetic signals compared to any other mitochondrial gene. The third-position nucleotides of COI present a high incidence of base substitutions, resulting in a three-time faster rate of molecular evolution than 12S or 16S rRNA (Knowlton and Weigt, 1998).

Molluscs are a worldwide distributed and diverse group consisting of a greater number of species than mammals, birds, reptiles, amphibians and fishes combined (Lydeard and Lindberg, 2003). Some morphological characters, commonly used in taxonomy of bivalve Molluscs, such as shell structure exhibit great divergence owing to population density and environmental effects (Lydeard and Lindberg, 2003). Oysters are among the most varying animals concerning the shell shape (Gunter, 1950). The Ostreidae family is a challenging group for taxonomists due to the great level of phenotypic plasticity. Morphological identification of some arcoid species is also particularly puzzling, due to the phenotypic variation and the presence of closely related taxa (Feng et al., 2011). Thus, DNA barcoding might be a useful tool for identification of oyster species (Liu et al., 2011) and other taxa with high levels of morphological plasticity. Veneridae is another diverse, marine bivalve subfamily. Members of this family are worldwide distributed to benthic communities of littoral environments and many of them are commercially important. The utilization of shell characters alone is questionable for investigating the phylogenetic relationships of this group (Kappner and Bieler, 2006; Chen et al., 2011a). The shell morphology as well as the morphometrics of pearl oysters are also highly variable and consequently has led to taxonomy confusion especially involving the species Pinctada radiata (Leach, 1814), Pinctada albino (Lamarck, 1819), Pinctada fucata (Gould, 1850) and Pinctada imbricate (Röding, 1798) (Yu and Chu, 2006). Pinctada radiata is originally distributed in the Indo-Pacific, whereas when arriving in the Mediterranean Sea it has been considered as one of the worst invaders (Streftaris and Zenetos, 2006). Nowadays populations of this bivalve inhabit the entire Mediterranean (Gavrilović et al., 2017; Theodorou et al., 2019; Moutopoulos et al., 2021; Png-Gonzalez et al., 2021).

Pelagic larvae of marine bivalves come up with a potent dispersal among conspecific populations promoting extensive gene flow. This gene flow is expected to occur over extended marine geographical scale owing to the lack of obvious barriers to pelagic larval scattering (Sá-Pinto et al., 2012). Passive dispersal leads to a lack of genetic heterogeneity between distant populations (Larsson et al., 2016). Furthermore, short duration of larval stages often leads to strong subdivision of bivalve populations (Ye et al., 2015). However, temporal, environmental or oceanographic factors instead of geographic factors may hamper gene flow with a subsequent effect on population sub-structure (Banci et al., 2017; Gutierrez et al., 2017). Genetic variability of distant populations may be affected by historical evolutionary events such as bottleneck effects, reduction in the effective population size and formation of barriers to gene flow (Geraghty et al., 2013). Despite larval dispersal, the divergence of specimens belonging to a local population is very common for marine invertebrates having pelagic larvae. The intrapopulation genetic heterogeneity is commonly attributed to physical or biological barriers to larval dispersal or to selection after settlement which leads to varying survival or fecundity of the recruits (Hedgecock, 1986).

The annual global production of bivalve molluscs for human consumption is estimated to be above 15 million tonnes corresponding to a total of 14% percentage of marine aquaculture production, and aquaculture provides about 90% of the edible bivalves (Wijsman et al., 2019). Bivalve mollusc aquaculture as well as natural bivalve populations, appear to be especially vulnerable to the threats arising from the climate change. Several factors are accountable for mass mortalities of bivalves (Fleury and Huvet, 2012; Soon and Zheng, 2020; Lattos et al., 2022), however, incidents of most mass mortalities are stimulated by the synergistic effects of more than one factor (Callaway et al., 2013). Harmful impacts might immerge from the increasing sea temperature, decreasing pH and elevated frequency of extreme climatic events, while possible synergies with other kinds of stressors, including diseases and harmful algal blooms could intensify the threat to marine bivalve molluscs (Rodrigues et al., 2015) leading to repeated mortality outbreaks. The capability of aquatic organisms to adjust to environmental disturbances heavily depends on the intraspecific genetic variability of the populations. Studying the genetics of marine bivalve molluscs plays a crucial role in understanding their evolution and adaptation to changing environments and also provides insights into their reproductive biology, physiology, and disease resistance, which have important implications for the sustainability of coastal ecosystems and the aquaculture industry. Obtaining population genetic knowledge is therefore crucial for both species management strategies (Sanna et al., 2017) and also for management practices in aquaculture (Ye et al., 2012).

Previous genetic studies of bivalves from Greece, including mitochondrial regions-based surveys, were mainly focused on Mytilus galloprovincialis (Lamarck, 1819) populations both cultured and wild ones (Kravva et al., 2000; Giantsis et al., 2012), indicating the absence of genetic structure probably enhanced by anthropogenic transportations (Giantsis et al., 2014). Apart from M. galloprovincialis, the genetic composition of Modiolus barbatus (Linnaeus, 1758) has been also studied investigating COI gene sequences from different populations from the Aegean Sea, providing the first DNA barcodes (Giantsis et al., 2019). Nevertheless, despite the great economic and ecological value and the great biodiversity levels that bivalves represent in Greece, no genetic data have been obtained for the remaining bivalve species from the Aegean and Ionian Seas. The genetic composition of Callista chione populations has not been investigated in any area of the Mediterranean Basin, while regarding the rest of the selected species in this study no research was focused on the Eastern Mediterranean.

Hence, in this study six marine bivalve species were gathered from a total of eight sampling sites in Greece. Three of the studied species belong to the family Veneridae [Venus verrucosa (Linnaeus, 1758), Callista chione (Linnaeus, 1758), Ruditapes decussatus (Linnaeus, 1758)] and the others were Arca noae (Linnaeus, 1758), Ostrea edulis (Linnaeus, 1758) and Pinctada radiata (Leach, 1814). All of the above species are commercially important. Ruditapes species had a price of 19.99 EUR per kg in the Italian market in 2021 (FAO, 2022, accessed online 24 May 2023) while the European flat oyster Ostrea edulis has an average price 3–5 times greater than Crassostrea gigas (Thunberg, 1793) (FAO, 2024).Venus verrucosa had in 2023 a price of about 18 EUR per Kg in Greek markets (personal information) while the market price in Croatia for Arca noae reached 7 EUR per kg in 2005 (Peharda et al., 2006). Furthermore, apart from their high commercial value, all the aforementioned species are of high export orientation.

Due to the lack of knowledge, the main objective of the present study was the evaluation of the population genetic structure of the aforementioned six bivalve species in different regions of Greece's territorial waters due to their commercial and ecological importance. In addition phylogenetic analyses were performed so as to investigate the taxonomic relationships of these bivalves with their closely related species and the relatedness of the derived haplotypes with haplotypes from other areas. Therefore, COI and 18S rRNA were sequenced to investigate the phylogenetic relationships and subsequently mtDNA (COI) sequences were used to assess the levels and the patterns of genetic diversity of the populations.

Material and Methods

Sampling areas and sampling procedure

Adult marine bivalve molluscs were collected from eight areas of Greece (Figure 1, Table 1) by professional scuba divers. Species identification of the sampled specimens was based on morphology. Samplings were performed from April 2020 to September 2022. After capturing, animals were transferred on ice to the laboratory, where a small piece (20–30 mg) of muscle tissue was dissected with sterilized tools from every animal and kept at −20°C until DNA extraction.

Figure 1. Geographic locations of the sampling sites and the species collected from each site.

Table 1. Genetic diversity indices of the bivalve populations from the different sampling areas

n, Sample sizes; H _n, haplotype number; H _d, haplotype diversity; π, nucleotide diversity.

DNA extraction, PCR amplification and sequencing

DNA extractions of all specimens were performed using the NucleoSpin Tissue^® DNA Extraction Kit (Macherey-Nagel, Düren, Germany) and according to the manufacturer's protocol. Genomic DNA was extracted from a piece of adductor muscle to avoid co-extraction of the male-type mitochondrial DNA (Plazzi and Passamonti, 2010; Kitada et al., 2013). Quality (260_nm/280_nm and 260_nm/230_nm) and concentration of extracted DNA were evaluated using the Quawell UV-Vis 5000 spectrophotometer (Quawell Technology, San Jose, CA, USA) and DNA was kept in the freezer (−20°C) until PCR amplifications. Primers used for the amplification of 18S rRNA and cytochrome c oxidase subunit I (COI) are listed in Table 2. Amplifications were performed in 0.2 ml PCR tubes on a Fast Gene ULTRA Cycler (Nippon Genetics, Japan). The reaction mixture (total volume 20 μl) contained 0.5 μl of each primer (10 pmol μl⁻¹), 10 μl of Fast Gene Taq 2x Ready mix (Nippon Genetics, Europe), 30 ng of extracted DNA as template and nuclease-free water up to 20 μl. Thermal cycling regime was an initial denaturation step of 94°C for 3 min followed by 38 cycles of 95°C for 30 s, 48–54°C (Table 2) for 35 s and 72°C for 40 s or 1 min (Table 2) followed by a final cycle of 72°C for 8 min. The size and quality of each PCR product were checked on 1.5% agarose gel stained with ethidium bromide. PCR products were purified using the NucleoSpin^® Gel and PCR Clean-up kit (Macherey-Nagel, Düren, Germany) and samples were subjected to Sanger sequencing.

Table 2. Primers, annealing temperature and extension time used for the PCR amplification

Plylogenetic analyses and population genetics

Chromatograms obtained from Sanger sequencing were visually checked and analyzed using BioEdit (Hall, 1999) and Finch TV 1.4.0 (Geospiza, Inc.). All sequences derived both from 18S rRNA and from COI genes were aligned in MEGA X (Kumar et al., 2018) using the clustal-W (Thompson, 1994) and MUSCLE algorithm respectively, separately for each species. The obtained haplotype sequences were then deposited in GenBank under the accession numbers OR127258-OR127275 (18S rRNA) and OR126934-OR126983 (COI) (Supplementary Tables 1 and 2). Sequences from GenBank database corresponding to closely related species were also downloaded (Supplementary Tables 1 and 2) and included in the dataset for the construction of the phylogenetic trees. Total aligned sequences in any case were subjected to test for selection of the best fit to the data substitution model before tree construction through the Maximum Likelihood (ML) statistical method. Best-fit model was chosen through the Akaike Information Criterion (AIC) which was automatically calculated by MEGAX software. ML trees were created in MEGA X. Support for nodes in the trees was determined with the bootstrap confidence levels making use of 1000 replicates.

Nucleotide variation at the region amplified from the nuclear 18S rRNA gene (about 535 bp) was not informative being identical to all specimens belonging to the same species (except for one polymorphism in Callista chione), thus these sequences were not further considered for population genetics. Only mtDNA (COI) sequences were used to define levels and patterns of genetic diversity of the bivalves within and among the sampling sites. Genetic heterogeneity within populations was estimated as haplotype diversity (H_d) and nucleotide diversity (π) using DNAsp v5 (Librado and Rozas, 2009). Molecular variation and corresponding PhiPT values among and within populations and also pairwise PhiPT values were calculated through GenAlEx 6.5 (Peakall and Smouse, 2006, 2012).

Results

Plylogenetic analyses through 18S rRNA

A 535 bp fragment from the 18S rRNA gene was sequenced and aligned from eight specimens for each species and each collection site. The derived data were used for the Maximum Likelihood tree construction together with sequences from closely related species which were obtained from GenBank. The phylogenetic tree of 18S rRNA is illustrated in Figure 2. Based on the depicted phylogeny, different families of bivalves are clustered separately as does the superfamily of Arcoidea. Arca noae sequence from Greece presented a large similarity to two Arca navicularis (Bruguière, 1789) 18S rRNA sequences (Table 3). One SNP (99.81% similarity) was detected in our sequence compared to two Arca noae sequences from Croatia and the USA (Figure 2). Pinctada radiata 18S rRNA demonstrated over 99% sequence similarity to many other congeneric species. Interestingly, the 18S rRNA sequences of Pinctada fucata and Pinctada maculata (Gould, 1850) were almost identical to our sequence, as was Pinctada radiata from Venezuela, while slightly lower similarity was found relative toPinctada radiata from USA (Table 3). The 18S rRNA partial gene sequence of Ostrea edulis was identical to an Ostrea chilensis (Küster, 1844) deposited sequence and was also very alike to many other oyster species (Table 3). Interestingly, the invasive Ostrea stentina (Payraudeau, 1826) collected from Thermaikos was not clustered in the same clade but there was small node support (Figure 2). Ruditapes decussatus18S rRNA-derived sequence led to its clustering with conspecific specimens from other areas but Polititapes rhomboides (Pennant, 1777) and Venerupis corrugata (Gmelin, 1791) were also grouped within the same branch (Figure 2). Callista chione is found within the same branch with Callista chione from Spain and Italy, while Venus verrucosa was grouped in the same branch both with other Venus verrucosa specimens from Europe and also with Dosinia exoleta (Linnaeus, 1758) and Venus casina (Linnaeus, 1758). The 18S rRNA sequences retrieved from the three venerid bivalves in this study had greater than 98.6% similarity and also all of them had large sequence affinity to many relative species (Table 3). All the aforementioned inferences indicate the invalidity of 18S rRNA gene as a DNA barcode for marine bivalves.

Figure 2. Phylogeny of 18S rRNA haplotypes made using the maximum likelihood method. Superfamily of Arcoida and families of Margaritidae, Ostreidae and Veneridae are depicted with different colours. Bold sequences are sequences obtained in this study. All sequences from related species which were included in the analyses exceeded 98% similarity to the examined species. Parentheses indicate unaccepted species names as deposited in GenBank. Bootstrap values above 50 are demonstrated on the tree.

Table 3. 18S rRNA sequence similarity of the studied bivalves with other deposited sequences in GenBank

Plylogenetic analyses through COI

Contrariwise to 18S rRNA, COI haplotypes in the present study can clearly distinguish and characterize the three species belonging to Veneridae. The different subfamilies where the Venerids belong were clustered separately, while closely related species were placed into different clades (Figure 3). The haplotypes retrieved from V. verrucosa in this study were more genetically related to Mediterranean haplotypes from Spain, Portugal and Croatia than more distant Atlantic haplotypes such as ones from west France and Ireland. Similarly, the retrieved haplotypes of R. decussatus were more similar to the haplotypes from Turkey and specifically from Umurbey and Izmir which are located closer to Vistonikos Bay than Fetiye which is on the western side of Turkey (Figure 3). In the phylogenetic tree (Figure 3) there are two well-supported sister groups of Ruditapes decussatus whose genetic divergence is in accordance with their geographic locality. One group is representative of eastern Mediterranean and includes sequences from the present study, and the second group is representative of western Mediterranean and European Atlantic coastlines. The three haplotypes of COI sequence affiliated with A. noae were 98.96–99.85% similar to two A. noae deposited sequences (Figure 4). O. edulis COI sequence can discriminate this species from other Ostreidae since only one species was found to possess highly similar sequence. This species was Ostrea angasi (Sowerby, 1871) which had 98.65% similarity to a haplotype from Maliakos (Figure 5). Calculating the genetic distances of the first 15 most similar sequences to our O. edulis sequence and also the first 15 most similar to O.angasi (AF540598) sequences from GenBank (these haplotypes were not included in the tree due to visualization purposes), the genetic distances within species was 0.0029 for O. angasi and 0.0068 for O. edulis (average 0.00485) while between species it was 0.0167. The conspecific genetic distances did not exceed the genetic distances between different species based on COI gene of these oysters within Ostrea genus. Consequently a barcoding gap exists between intraspecific and interspecific divergences. Hebert et al. (2004) suggested a standard threshold 10x the mean intraspecific variability for the studied group. The two oyster species exhibited a 0.485% average intraspecific variability, thus a 4.85% threshold could identify the two species.

Figure 3. Phylogeny of COI haplotypes from Veneridae made using the maximum likelihood method. All different subfamilies are depicted with different colour and bold sequences correspond to sequences obtained in this study. Since all the examined species exhibited lower than 90% similarity from the closest relative, 2–3 closely related species were included from GenBank along with conspecific haplotypes. Bootstrap values above 50 are demonstrated on the tree.

Figure 4. Phylogeny of COI haplotypes from Arca noae made using the maximum likelihood method. Bold sequences with different colour correspond to sequences obtained in this study from the different sampling areas. Since Arca noae exhibited lower than 75% similarity from the closest relative (except of Tetrarca tetragona), two additional relative species were included from GenBank along with conspecific haplotypes. Bootstrap values above 50 are demonstrated on the tree.

Figure 5. Phylogeny of COI haplotypes from Ostrea edulis made using the maximum likelihood method. Bold sequences with different colour correspond to sequences obtained in this study from the different sampling areas. Since Ostrea edulis exhibited lower than 89% similarity from the closest relative (except of Ostrea angasi), two additional relative species were included from GenBank along with conspecific haplotypes. Phylogenetic relationship with Ostrea stentina is also shown. Bootstrap values above 50 are demonstrated on the tree.

Invasive species genetic analyses

Interestingly, one oyster haplotype from Thermaikos Gulf, which was morphologically characterized as Ostrea edulis by the divers, considerably differed from the rest of O. edulis haplotypes. The 18S rRNA sequence of this specimen exhibited three different nucleotides while the aligned COI sequence revealed great differences compared to O. edulis. This sequence was subjected to BLAST and found to possess greater than 99% sequence similarity to many Ostrea stentina deposited sequences in GenBank. This COI sequence was clustered within the same cluster as O. stentina from Spain and Japan and was slightly different from USA and Mauritus-originated haplotypes (Figure 5). With respect to the invasive Pinctada radiata, the unique COI sequence was at least 8.5% different from the closest congeneric species (Figure 6). The 18S rRNA sequence of Pinctada radiata was identical to two congeneric species and extremely related to several other species (Figure 2).

Figure 6. Phylogeny of COI haplotype from Pinctada radiata made using the maximum likelihood method. Bold sequences with different colour correspond to sequences obtained in this study from the different sampling areas. Since Pinctada radiata exhibited lower than 92% similarity from the closest relative, four additional relative species were included from GenBank along with a conspecific haplotype from United Arab Emirates. Bootstrap values above 50 are demonstrated on the tree.

Population genetics

Veneridae

Ruditapes decussatus specimens were gathered only from Vistonikos Bay. Two polymorphic sites were detected in COI fragment of 685 bp, resulting in a nucleotide diversity of 0.00075 and three haplotypes (Table 1). The haplotype diversity was found 0.6 among Callista chione individuals that were collected from three areas (Table 1). A total of three polymorphic sites were also detected in 685 bp of COI with one polymorphism being present in all three populations and one polymorphism only in Saronikos Gulf, yielding three different haplotypes in this area, while two different haplotypes were detected in Maliakos Gulf and Chalkidiki. Site nucleotide diversity and haplotype diversity varied from 0.00049 to 0.00111 and 0.33 to 0.67, respectively (Table 1). PhiPT values were negative and statistically not significant among all areas and all molecular variances found distributed within the populations (Table 4). Venus verrucosa was also sampled from three areas. The greatest number of polymorphic sites (overall 18) in COI fragments was found in this species resulting in an overall number of 20 different haplotypes (Table 1). Haplotype diversity of the populations was consequently great varying from 0.89 to 0.98. The lowest nucleotide diversity was estimated in Lesvos Island and the greatest in Saronikos Gulf (Table 1). PhiPT values were low or negative and statistically not significant between the populations while 3% of the observed molecular variance was found among the three populations and 97% within the populations (Table 4).

Table 4. Analysis of molecular variance (AMOVA)

Pinctada radiata and Ruditapes decussatus were excluded, as sequences of the former were identical in all sampling areas and the later was sampled only from one area. Asterisk (*) depicts statistically significant differences at P < 0.05.

Ostrea edulis

Ostrea edulis was gathered from two Gulfs, namely Thermaikos and Maliakos. A total of five polymorphic sites were observed in 693 bp of COI yielding a total number of six different haplotypes. Haplotype diversity and nucleotide diversity were 0.7 and 0.00116 in Thermaikos Gulf while in Maliakos Gulf these values were 1 and 0.0077, respectively (Table 1). PhiPT value between these populations was 0.35 but barely statistically insignificant (P = 0.05, Table 5). Molecular variance was shared among and between populations, being 35% and 65%, respectively.

Table 5. Analysis of pairwise PhiPT values

Pinctada radiata and Ruditapes decussatus were excluded, as sequences of the former were identical in all sampling areas and the later was sampled only from one area. Asterisk (*) depicts statistically significant differences at P < 0.05.

Arca noae

Arca noae were sampled from Thermaikos Gulf and Lesvos Island. Three SNPs were found in COI fragment of 703 bp, yielding an overall number of four haplotypes, from which two were found in both sampling areas, while the other two haplotypes were observed only in Thermaikos Gulf. As a result, animals originating from Thermaikos Gulf exhibited increased haplotype and nucleotide diversity (Table 1). PhiPT value of the two areas was negative and not statistically significant and all estimated molecular variance was attributed within the two populations (Table 4).

Pinctada radiata

Pinctada radiata specimens were collected from Kefalonia Island, Maliakos Gulf and Evoikos Gulf. Sequences obtained from a 679 bp fragment of COI were identical to every single specimen of the 17 sampled individuals.

Discussion

Evaluation of applicability of the different barcoding markers

The determination of species by DNA barcoding relies on the sequence variation within and between the taxa (Mikkelsen et al., 2007). The efficacy of barcoding is eliminated proportionally from overlaps between genetic variation within species and overlap in divergence distinguishing sister species (Meyer and Paulay, 2005). The exclusive utilization of 18S rRNA sequencing was not able to identify the presence of Ostrea stentina as it had almost identical sequence with Ostrea edulis, exhibiting only three nucleotide polymorphisms. DNA barcoding using the COI fragment, on the other hand, proved proper for distinguishing two different oyster species in Thermaikos Gulf. Hamaguchi et al. (2017) recorded for the first time the presence of O. stentina in two distant areas of Japan employing also mitochondrial DNA markers. The high similarity of COI sequence between O. angasi and O. edulis triggered our calculations of genetic distances and revealed that DNA barcoding depending on barcoding gap and the threshold approach (Hebert et al., 2004) can effectively identify even these closely related oyster species. Liu et al. (2011) working on mitochondrial genes (COI and 16S ribosomal RNA) highlighted the efficacy of species identification in Ostreidae family via DNA barcodes and phylogenetic analysis. Hamaguchi et al. (2017) also declared a clear distinction among native Japan oyster species by DNA barcoding using mitochondrial markers including COI.

As observed with the oysters in the present study, the similar morphology as well as the sympatric distribution of some bivalve species may frequently prevent their accurate identification. Consequently, incorrect labelling of edible species is repeatedly observed along the production chain (Harris et al., 2016; Giusti et al., 2020; Parrondo et al., 2021) while proper population management is hindered (Vierna et al., 2014). Misidentification can lead to incorrect estimations of population size, growth, and reproductive potential, leading to flawed management strategies and also to inappropriate regulations, causing overexploitation or underutilization of certain species. Additionally, undetected invasive bivalve species can have detrimental effects on native ecosystems and their false identification, or their overlook may result in inadequate control practices, allowing invasive species to proliferate and outcompete native species.

Food fraud often occurs in markets, where cheaper morphologically similar bivalve species are sold at the prices of more expensive species (e.g. Costa et al., 2008). COI efficacy in discriminating the different bivalve species highlights this marker as a decisive tool for detecting mislabelled or substituted species. Since sometimes COI alone may not be informative, a second genetic marker could be applied along with COI. As for the provenance, COI could be combined with the fatty acid profile and/or isotope analysis since they vary among different species and also among geographic locations and have been already validated as useful tracers in bivalves (Zhang et al., 2019; Fonseca et al., 2022). In bivalves that are generally characterized by similar morphology combined with low genetic differentiation, the combination of at least two of the aforementioned methodologies is generally proposed. In addition, stringent regulations and penalties for food fraud should be implemented by local authorities aiming to minimize the deliberate food fraud.

The 18S rRNA sequences were proved inappropriate for unveiling taxonomic issues to all examined bivalves. P. fucata and P. maculata were identical to P. radiata GenBank deposited sequences, while P. radiata sequence including our haplotype appeared greater than 99% similarity to many congeneric species. On the other hand, COI sequence clearly discriminated P. radiata from congeneric species and could be a useful marker to discriminate the rest of Pinctada species. 18S rRNA sequence was also confusing for the investigation of phylogenetics of Venerids, Ostrea edulis and Arca noae as well. The three species of the Veneridae family presented greater than 99% 18S rRNA sequence similarity to many other species even from different genera, while the similarities between them were greater than 98.6%. COI proved an exceptional marker for delimiting the venerid species studied herein from their closely related congeneric samples. Chen et al. (2011b) have also proved the efficacy of COI in distinguishing 11 species of the Tapetinae subfamily. The efficacy of cytochrome oxidase subunit I (COI) revealed in this study is not unusual as its strong signal has proved to be appropriate for resolving the phylogenetic relations of bivalves (David and Savini, 2011). Mitochondrial genes can be particularly informative markers, especially concerning shallow nodes (Plazzi et al., 2011; Sharma et al., 2012).

It should be noted that despite the fact that in most animals mtDNA is inherited strictly by maternal lineage, there is a unique mode of mitochondrial DNA (mtDNA) transmission within bivalve molluscs (Ladoukakis and Zouros, 2017) that is called doubly uniparental inheritance (DUI).

DUI is defined by the existence of a female (‘F’) mitochondrial genome transmitted from mothers to both daughters and sons, and a male (‘M’) mitochondrial genome exclusively transmitted from fathers to sons. This intriguing inheritance pattern has been observed in several bivalve species, including the Veneridae family (Chacón et al., 2020). The potential presence of DUI sex-linked mitochondrial DNA should be carefully considered, as it may lead to inaccurate clustering and false calculation in population genetics indices. To avoid extracting the male-type mitochondrial DNA, DNA extractions in this study were performed from adductor muscle which carries very little M-type mtDNA in DUI species (Plazzi and Passamonti, 2010; Kitada et al., 2013). Additionally, chromatograms and obtained sequences of COI were carefully checked before the analysis to exclude the case of potential double peaks.

Genetic diversity and differentiation

Haplotype diversity (H _d) of Ruditapes decussatus was found in this study 0.6. Habtemariam et al. (2015) found lower levels of H _d and similar to us nucleotide diversity in a 503 bp fragment of COI studying two wild populations from northern Spain. Sanna et al. (2017) also using COI gene observed the absence of genetic structure within the western Mediterranean and population substructuring among the Aegean and Marmara seas and other populations from the Mediterranean basin and the South European Atlantic coast as well. Another study concerning R. decussatus also revealed low intrapopulation differentiation (Borsa et al., 1994). The divergence of Ruditapes decussatus between the two observed groups in Figure 3 (Eastern Mediterranean and Western Mediterranean-European Atlantic coasts) could be suggestive of the occurrence of sister species within R. decussatus and should be further investigated in the future. Callista chione exhibited a total of three different haplotypes among the three sampling locations and just lower mean haplotype diversity than R. decussatus. All molecular divergence detected was distributed within the populations of C. chione. The third Venerid, Venus verrucosa exhibited high values of nucleotide and haplotype diversity at all sampling areas, while the majority of the observed molecular variance (97%) was found distributed within the different populations.

P. radiata owns a relatively long-lived pelagic larval phase. Al-Saadi (2013) found very low genetic differentiation in Pinctada radiata among different locations across the Qatari coast, while all variation in COI sequence was detected within the sampling areas. In our study this species was found panmictic at three sampling geographic areas, presenting zero differences in the both sequences of COI and 18S rRNA. Other authors also detected low haplotype and nucleotide diversity in P. radiata specimens among different sampling areas (Meyer et al., 2013; Barbieri et al., 2016) while at Adriatic Sea (Sobra Bay) and at United Arab Emirates high haplotype along with low nucleotide diversity was found by Gavrilović et al. (2017) and Cunha et al. (2011) respectively. Although these sequences were not analysed in comparison to our haplotypes owing to their much shorter available sequences, our results are in total agreement with these inferences concerning haplotype and nucleotide diversity. BLAST results of the P. radiata COI sequence obtained from the three locations of Greece were 99.22–100% similar to the ten haplotypes reported by Barbieri et al. (2016) from the Eastern Mediterranean, while the haplotype H2 from these authors which was the prevalent haplotype among seven geographic areas of the Eastern Mediterranean was 100% similar to our haplotype. The low level of population subdivision seems to be mainly assigned to the long pelagic larval stage of this bivalve.

Ostrea edulis exhibited high haplotype diversity at the two sampling sites and was the only species whose observed molecular variance was statistically significantly distributed both among (35%) and within (65%) of the populations. The PhiPT value between the sampling areas was barely statistically insignificant (P = 0.052) which may be attributed to the small sampling size since only nine oysters were examined in total. Studies on allozyme differentiation of the flat oyster (Jaziri, 1990; Saavedra et al., 1995) found that the overall population differentiation was low with a notable divergence among Mediterranean and Atlantic populations. Sobolewska and Beaumont (2005) examining neutral microsatellite loci of the European flat oyster from Scotland, Ireland, France, Netherlands and Norway noticed a generally low genetic subdivision of the populations. Šegvić-Bubić et al. (2020) also working on microsatellites found relatively low to moderate differentiation between 28 individual populations of the Adriatic Sea. The short pelagic larval duration of this species or a differential selection after larvae settlement at the two sampling sites could explain the observed increased molecular variance among the studied populations. Moreover, a blockage to bivalve larval dispersal may be present among Maliakos and Thermaikos Gulf but this suggestion should be further investigated since we did not sample other bivalve species from these two sites.

Intraspecific genetic diversity enhances the potential of animals to overcome environmental or ecological challenges (Lind et al., 2007). All bivalve populations examined in this study except for Pinctada radiata exhibited medium to high haplotype diversity and also a relative low nucleotide diversity in COI. These results suggest a high number of closely related haplotypes. Genetic knowledge of marine bivalve molluscs has the potential to strongly support conservation efforts, enhance aquaculture production and assist the long-term viability of these economically and ecologically important species. The total molecular variance of Callista chione, Arca noae and Venus verrucosa was found distributed within the populations. The different populations of the three aforementioned species along with the panmictic P. radiata appear to constitute a single management unit in the examined areas. C. chione and A. noae had the lowest values of intrapopulation molecular variance. On the other hand, Ostrea edulis exhibited a population subdivision among Thermaikos and Maliakos Gulf. Moreover, O. edulis had the second greater mean haplotype and nucleotide diversity among the examined six bivalve species since only V. verrucosa presented greater values. Employing a more informative nuclear gene for bivalve population studies such as internal transcribed spacer 2 (David and Savini, 2011) or microsatellite loci, along with COI gene would lead to more reliable conclusions concerning the population genetics of these bivalve species across the Greek coastline.

Distribution and abundance of invasive bivalves

Pinctada radiata is ranked as one of the most successful and widely distributed invasive marine organisms in the Mediterranean Sea (CABI, 2016). It has spread throughout the Mediterranean basin and has been recorded in several parts of Greece (Theodorou et al., 2019; Moutopoulos et al., 2021). Its simultaneous portraiture as a biofouling agent and a potential fisheries product puzzles the situation of the impact of this invasion. Although it is believed that it presents high levels of adaptation in new environments (Mohamed et al., 2006; Katsanevakis et al., 2008), the genetic profile from Greek Seas, demonstrates extremely low levels of genetic diversity as well as a panmictic homogeneous genetic composition. These data indicate that, yet, Pinctada radiata has not developed a distinct genetic identity in this marine area. Similarly, no genetic diversity has been observed in other places where this invaded. It has been assumed that invasions in different countries have different origins (Zenetos et al., 2004). This absence of genetic diversity provides evidence for a generally recent migration event in parallel with the widespread potential.

On the other hand, Ostrea stentina is an oyster species clearly differentiated from Ostrea edulis, which inhabits several marine areas where other oyster species such as O. edulis are distributed, in sympatric populations. Here we report for the first time its presence in Eastern Mediterranean, corresponding to a new invasion with unknown effects. Genetic analysis showed a complete absence of genetic distance in comparison with haplotypes from Spain and Japan, indicating a recent migration pattern along with a fast geographic spread as well. Although the two molecular markers applied in the present study are not sufficient to characterize the exact route of this invasion and examine if it corresponds to a random via shipping invasion enhanced by the pelagic larvae, the case of human-mediated expansion through trade, that has been also proposed for other marine bivalves in the area (Giantsis et al., 2019), cannot be excluded.