Mitochondrial DNA Reveals Genetic Structuring of Pinna nobilis across the Mediterranean Sea

Pinna nobilis is the largest endemic Mediterranean marine bivalve. During past centuries, various human activities have promoted the regression of its populations. As a consequence of stringent standards of protection, demographic expansions are currently reported in many sites. The aim of this study was to provide the first large broad-scale insight into the genetic variability of P. nobilis in the area that encompasses the western Mediterranean, Ionian Sea, and Adriatic Sea marine ecoregions. To accomplish this objective twenty-five populations from this area were surveyed using two mitochondrial DNA markers (COI and 16S). Our dataset was then merged with those obtained in other studies for the Aegean and Tunisian populations (eastern Mediterranean), and statistical analyses (Bayesian model-based clustering, median-joining network, AMOVA, mismatch distribution, Tajima’s and Fu’s neutrality tests and Bayesian skyline plots) were performed. The results revealed genetic divergence among three distinguishable areas: (1) western Mediterranean and Ionian Sea; (2) Adriatic Sea; and (3) Aegean Sea and Tunisian coastal areas. From a conservational point of view, populations from the three genetically divergent groups found may be considered as different management units.


Introduction
The assessment of genetic variation in marine species and the evaluation of connectivity among populations are crucial components for conservation purposes and resource management [1], [2]. However, direct estimates of dispersal are difficult to obtain for marine organisms [3] because of the multiple factors (e.g., pelagic larval duration, larval phylopatry, regional and local water circulation, and habitat specificity) that influence dispersal capability [4]. Molecular tools are often employed to make indirect inferences about the levels of connectivity among natural populations [5], [6]. Molecular approaches can be particularly useful for determining the influence of such factors as hydrodynamic conditions [7], [8], habitat specificity [9], [10] and various ecological parameters (e.g., [11], [12]) on gene flow and spatial genetic structure. However, the spatial genetic structure of a species may reflect not only its current dispersal capability but also its phylogeographic history. For example, the occurrence of different genetic lineages in Mediterranean species has been explained as the result of either the present dispersal capability or past geomorphologic processes or both (see [13] and references therein). Furthermore, human activities (e.g., intense exploitation, industrialisation, and coastal tourism) often lead to significant reduction of population size, which may result in an increased magnitude of genetic drift.
The fan mussel Pinna nobilis (Linnaeus, 1758) (Mollusca: Bivalvia) is an endemic Mediterranean species of great conservation interest that can be regarded as a flag species sensu Walpole and Leader-Williams [14] and Heywood [15]. As one of the largest bivalves, with maximum lengths in excess of 1 m, this popular species is able to capture the public imagination [16]. Pinna nobilis occurs at depths between 0.5 and 60 m in soft bottoms primarily characterised by seagrass meadows [16]. Its veliger-stage larvae drift in the water column [17] before settling in the sediment and anchoring via a byssus [18]. Although De Gaulejac and Vicente [19] hypothesised a period of 5-10 days, little is known about the length of the larval life cycle of P. nobilis. Moreover, Peharda and Vilibić [20] postulated that P. nobilis veliger might exhibit daily vertical migration, spending the daylight hours in darker, deeper waters and occupying superficial layers at night, as observed for other bivalves [21].
As a result of the population declines associated with human activity (primarily the harvesting of byssus for the manufacture of so-called ''sea silk'' and shell collecting [22], [23], secondarily the use as food of the abductor muscle in some Mediterranean regions [24], and the collection of the little, with poor commercial value, pearls [25]), P. nobilis has been included in the list of Mediterranean endangered species (Annex IV of the Habitat Directive and Annex II of the Barcelona Convention) since 1995. Remarkably, it is likely that some of the known geographic and ecological barriers [26], [27], [28], which define different marine ecoregions and biogeographic sectors in the Mediterranean [29], [30], may affect the degree of connectivity among P. nobilis populations. In the Mediterranean region, several authors have documented population subdivision in relation to physical barriers in other marine species with a high dispersal potential (e.g., [13], [31], [32]).
Despite the conservation importance of P. nobilis, few studies have investigated the biology and ecology of this species (e.g., [20], [23], [33], [34]), and fewer still have studied its population genetics. To date, only two molecular surveys have been conducted, one in the Aegean Sea [35] and one along the Tunisian coasts [36]. The former [35] surveyed the genetic variability of four Aegean populations by means of the mitochondrial DNA (mtDNA) markers Cytochrome c Oxidase subunit I (COI) and 16S ribosomal DNA (16S). A high level of haplotypic diversity was found for the COI gene, whereas the 16S gene showed a lower level of variability. These results suggested the lack of genetic structuring among the Aegean populations. In the second study [36], the authors used COI sequences to describe the genetic variability of five populations from the northern, eastern, and southern Tunisian coasts. A North-East decreasing gradient of genetic variability was found among these populations, which was explained in light of the variance in the hydrodynamic regime of the areas analysed.
Knowledge of the amount of genetic variability and distribution in space and time is crucial for a correct diagnosis of the conservation status and viability of populations and the threats to them [37], [38]. Accordingly, it is evident that a molecular analysis performed on a larger number of Mediterranean P. nobilis populations could provide a deeper insight into the possible influence of environmental or anthropogenic stress on the population dynamics and evolutionary history of this endangered species, allowing the development of effective conservation measures. Based on this background, the present study was aimed to investigate the large-scale patterns of spatial genetic variation of P. nobilis in the Mediterranean. Our primary focus was the appraisal of the genetic structuring and population connectivity of P. nobilis in the western Mediterranean, Ionian Sea, and Adriatic Sea ecoregions. We surveyed populations from the Sardinian-Corsican region, Elba Island, Sicily, and the Venetian Lagoon and compared these populations with the available data from the Aegean Sea [35] and Tunisian coasts [36]. This large area is characterised by a complex and variable pattern of marine circulation (e.g., [39]) and encompasses several biogeographic sectors [28]. To accomplish the objective of the present study, we employed two mitochondrial markers (COI and16S) previously used for the investigation of Aegean [35] and Tunisian populations [36].

Ethics Statement
No field studies involving impacting manipulation, dislocation, or removal of Pinna nobilis individuals were performed. For each location under protection all necessary permits were obtained for the sampling activities by the authority responsible for each protected area: - No specific permissions were required for the locations placed out of protected areas, since those locations are not privatelyowned or protected in any way.

Sampling
Overall, we sampled 236 specimens of P. nobilis from 25 locations (Table 1 and Figure 1), sited in the following ecoregions according to Spalding et al. [30]: (1) the central part of the western Mediterranean (the Sardinian Sea, the North Tyrrhenian Sea, the South Tyrrhenian Sea, and the Strait of Sicily); (2) the Ionian Sea (the South-East of Sicily); and (3) the Adriatic Sea (the Venetian Lagoon). Two individuals from the Levantine Sea (Cyprus) were also analysed as the easternmost Mediterranean outliers. We developed a specific non-lethal sampling method, performed by SCUBA divers, which does not cause significant damages to the shell and soft tissues of P. nobilis. The valves of a given individual were held open with a wooden stick (diameter = 0.5 cm), put in proximity (4-5 cm) of the hinge ligament, and a 20-50 mg sample of mantle tissue was excised using Hartman Alligator Forceps 3.5. The stick was then removed, and the tissue sample stored in a 5 ml tube. The tissue sample was then transferred to the laboratory in a refrigerated box and there preserved in 75% ethanol. This method ensures the survival of the sampled individuals.
Each 25 ml PCR mixture contained approximately 100 ng of total genomic DNA, 0.32 mM of each primer, 2.5 U of EuroTaq DNA Polymerase (Euroclone), 16 reaction buffer and 200 mM of each dNTP. MgCl 2 concentration was set at 3 mM, and 12.5 mg of BSA were added to the reaction mixture. PCR amplifications were performed according to the following steps: 1 cycle of 2 min at 94uC, 35 cycles of 1 min at 94uC, 1 min at 46uC and 1 min 30 s at 72uC. A post-treatment of 5 min at 72uC was applied. After electrophoresis on 2% agarose gels, the PCR products were purified using ExoSAP-IT (USB Corporation) and sequenced using an external sequencing core service (Macrogen Europe). Dual peaks of similar height, which could be interpreted as evidence of mitochondrial pseudogenes in the nucleus (Numts) or heteroplasmy, were not observed in any of the electropherograms. The PCR products did not show any occurrence of aspecificity, excluding the possibility of multiple nuclear mtDNA-like sequences.

Statistical Analysis
Sequences were aligned using the program Clustal W [40] and deposited in GenBank (Accession Nos.: COI, JX854788-JX855023; 16S, JX854562-JX854787). Estimates of the number of polymorphic sites (S), the number of haplotypes (H), haplotype diversity (h), nucleotide diversity (p), and the mean number of pairwise differences (d) were obtained using the software package DnaSP 5.10 [41].
After having added to our COI and 16S sequences the COI [35], [36] and 16S [35] sequences available in the literature, we constructed two different datasets. One dataset included 236 COI sequences obtained in the present study along with 26 and 49 sequences from the Aegean Sea [35] and Tunisian coasts [36], respectively (see Table 1 for details). The second dataset included the concatenated COI and 16S regions and consisted of the 219 sequences of the present study along with 25 from the Aegean Sea [35]. Concatenating these two mitochondrial genes was appropriate because COI and 16S did not show significant heterogeneity by partition-homogeneity test (P = 0.06) performed with PAUP* 4.0b10 [42]. a) Analysis of the COI dataset. The presence of population genetic structure was assayed by the Bayesian model-based clustering algorithm implemented in BAPS 5.2 [43]. Clustering was performed using the module for linked molecular data and by applying the codon linkage model, which is appropriate for sequence data. Each analysis was run ten times with a vector of Kvalues = 1 to 22, each with six replicates. Haplotypes were organised into haplogroups following the partition of sequences into the distinct genetic clusters evidenced by Bayesian clustering. The Bayesian clustering algorithms probabilistically assign individuals to groups based on nucleotide frequencies of DNA sequence data without presuming pre-defined populations [44]. For each sampling location, we computed the proportion of a given haplogroup to build maps of haplogroup frequency in the Mediterranean.
Genetic relationships among haplotypes were investigated by a median-joining network [45] using the software package Network 4.5.0.1 (www.fluxus-engineering.com). Mutations were inversely weighted according to the number of times they originated. Therefore, different weights were assigned to the point mutations: sites involved in a single mutational event = 90; sites involved in repeated mutational events = 90/number of repetitions.
Patterns of genetic differentiation at the population level were assessed by Arlequin 3.5.1.3 [46] using a matrix of Tamura and Nei's [47] genetic distances with a gamma correction according to the best-fitting model of sequence evolution obtained with jModeltest 2.1.1 [48]. First, we calculated pairwise W ST values between sampling localities. Samples with less than five individuals were excluded from population level analysis due to a lack of statistical power. We next conducted a hierarchical analysis of molecular variance (AMOVA) [49]. Alternative groupings of populations were taken into account. Initially, we ran analyses in which populations were grouped a posteriori according to the pattern of genetic structure suggested by the pairwise W ST values. Significance of both pairwise W ST values and AMOVA was assessed by a permutation test (with 10,000 random replicates). Then, we defined a priori regional groups that corresponded to the main biogeographic subdivisions of the Mediterranean [29]: (1) western vs. eastern Mediterranean; (2) western vs. eastern Mediterranean vs. Adriatic. Sicilian populations located close to the border separating the two Mediterranean basins were twofold tested to be either entirely part of western Mediterranean, or part of western (San Vito lo Capo-SVC, Mondello-MON, Milazzo-MLZ) and eastern (Pachino-PAC, Ognina-OGN) Mediterranean. Where necessary, probability values were corrected for multiple testing by applying the False Discovery Rate method [50].
Aspects of historical demography were inferred by examining the fit of the observed mismatch distribution of DNA pairwise differences to a model of demographic population expansion [51], [52]. Such distributions are unimodal when populations have experienced recent expansion and multimodal at demographic equilibrium or when populations are significantly subdivided. Mismatch distributions and relative Rogers' parameters [53] (t, h 0 ,  h 1 ) were assessed using Monte-Carlo simulations of 10,000 random samples as implemented in Arlequin 3.5.1.3. The fit of the observed mismatch distribution to the expected distribution under a model of demographic expansion was assessed by Monte-Carlo simulations of 10,000 random samples. The sum of squared deviations (SSD) between observed and expected mismatch distributions was used to test the probability of obtaining a simulated SSD that was larger than or equal to the observed one. Tajima's D [54] and Fu's Fs [55] neutrality tests were used to infer departures from population equilibrium models. Significant negative values of Tajima's D are expected to occur in cases of recent population expansion or after a selective sweep [55]. Positive values are expected in the case of balancing selection, population subdivision or recent bottlenecks [56]. Similarly, significant negative Fu's Fs values indicate an excess of rare haplotypes, which might be caused by recent population expansion [55], whereas positive values indicate balancing selection, population structure or moderate bottlenecks [57]. Combining different neutrality tests can help to distinguish among the different evolutionary processes responsible for departures from equilibrium; Fu's Fs can better detect demographic expansions, whereas Tajima's D can better detect bottlenecks and populations contractions [58].

b) Analysis of the concatenated COI and 16S
dataset. Statistical treatments were conducted as with the COI dataset. Note that for this dataset, redefined population groupings were required for some analyses (e.g., AMOVA) due to the lack of 16S sequences from Tunisian coasts. c) Estimates of divergence time. The modal value (t) of the mismatch distribution of pairwise genetic differences was also used to estimate time since expansion using the formula t = t/2u [51]. Here, t is the time since expansion in generations and u is the mutation rate of the entire DNA fragment, equal to u = ms, where m is the mutation rate per nucleotide per generation, and s is the sequence length. Uncertainty in the estimates of expansion time since expansion was accounted for computing the 95% confidence interval (CI) of t from the parametric bootstrap approach with 10,000 pseudo-replicates [52]. Time since expansion expressed in years was estimated assuming a constant molecular clock, a 2years generation time [17], and a divergence rate per nucleotide site per million years of 0.52% as estimated in Luttikhuizen et al. [59] for bivalves.
We also used the Bayesian coalescent approach implemented in BEAST 1.7.4 [60] to infer historical demography, as it allows improved recovery of the historical signal within DNA sequences [61]. As there was a lack of 16S sequences for Tunisian populations and of an appropriate divergence rate for 16S region for the species, analysis was carried out on the COI dataset only.
We used the Bayesian skyline plot model [62] assuming a piecewise constant model with ten coalescent intervals; Markov-Chain-Monte-Carlo simulations were run under the GTR+G+I model with four gamma categories using a relaxed uncorrelated molecular clock [62]. Divergence rate of the COI region was modeled using a uniform distribution that ranged from 0.14% to 0.52% (divergence rate per site per million years) according to Luttikhuizen et al. [59]. For each simulation three independent runs were carried out; depending on the size (number of sequences in the dataset) 100 million up to 400 million iterations were used (thinning was adjusted to obtain a final sample of 10,000 records). Tracer 1.5 [63] was used to assess convergence of runs and the Effective Sample Size (ESS) of each parameter. In order to obtain an adequate effective sample size (ESS $200), the three independent runs performed for each simulation were combined using the Logcombiner utility of the BEAST 1.7.4 package. The resulting file was used to estimate population size change through time, that was visualized by the Bayesian skyline plot computed with Tracer 1.5.

COI Dataset
Specific primers amplified an internal portion of COI 338 base pair (bp) long. Therefore, after the addition to our data of those from Katsares et al. [35] and Rabaoui et al. [36], we obtained a 338 bp sequence alignment from 311 Pinna nobilis individuals. Sixty-two different haplotypes (Table S1), defined by 45 polymorphic sites, were found. In our dataset (236 individuals), many haplotypes were found only at a single locality (63%). This percentage increased to 79% when the samples from the Aegean Sea [35] and Tunisian coasts [36] were included. Total mean haplotype and nucleotide diversity, calculated on 311 individuals, were h = 0.910 and p = 0.007, respectively. The lowest values of haplotype and nucleotide diversity were found among Aegean and Tunisian populations (h = 0.720 and p = 0.004 for Aegean Sea and h = 0.605 and p = 0.002 for Tunisian coasts); overall higher values of genetic diversity were found in the other sites. Estimates of COI genetic diversity for each site are reported in Table S2.
a) Genetic structuring. The Bayesian assignment analysis identified four distinct haplotype groups, hereafter denominated P1, P2, P3 and P4 (Table S3). The groups P1, P3 and P4 includeed samples from Sardinia, Corsica, Elba Island, Sicily, the Venetian Lagoon and Cyprus and lack internal geographic structuring, whereas the Aegean and Tunisian samples all exhibited P2 haplotypes, except for one Aegean individual from Epanomi-EP and one Tunisian specimen from Monastir-M, both belonging to group P3. The geographic distributions of haplogroups over the Mediterranean map are provided in Figure 2. The group P1 was the most widespread (41.8%), being present in Sardinia (with the exception of all samples from Isola di Mal di Ventre-IMV), Corsica, Elba Island, Sicily and the Venetian Lagoon (two individuals only). The group P2 was found in 23.5% of individuals and characterised populations from the Aegean Sea and Tunisian coasts. Although P3 was the group spread most widely across the Mediterranean (occurring in all populations except for the majority of the Aegean and Tunisian samples), it occurred at low frequency (26.7%), reaching higher proportion in areas as distant as the northern Adriatic and western Sardinia. The group P4 was characterised by the lowest average frequency (8%) and occurred in north-eastern Sardinia, southern Corsica, at Elba Island, northern Sicily (plus one individual from south-eastern Sicily, Pachino-PAC) and in the Venetian Lagoon (two individuals only).
The median-joining network analysis identified a low level of divergence among haplotypes, as, with a few exception, connections were characterised by one or two point mutations, (Figure S1, and Figure 3). Three haplotypes (PN1, PN6, PN7; Table S1), present in 134 of 311 (43.1%) individuals, were the most common among all samples, with the exception of the Aegean and Tunisian samples, of which 85.3% (64 of 75 individuals) showed two main haplotypes (PN52, PN54; Table S1). No evidence of geographic structuring among regions (all characterised by many instances of haplotype sharing) was detected for the entire sample, with the exception of populations from the Aegean Sea and Tunisian coasts that did not share haplotypes with the remaining populations ( Figure S1). When haplotypes were grouped on the network according to the results of Bayesian assignment, each cluster (P1, P2, P3 or P4) was generally characterised by one or two widelydispersed haplotypes, with radiating star-like branching patterns ( Figure 3).
The pairwise W ST estimates indicated significant genetic differentiation between the Aegean and Tunisian samples and the other sites (Table S4). The Venetian Lagoon was also significantly differentiated in most comparisons, with the only exceptions being three Sardinian populations (Le Saline-SAL, Porto Ottiolu-OTT, Marceddì-MAR) and one Sicilian (Milazzo-MLZ). The AMOVA (Table 2) maximized the largest differences among groups defined a posteriori when the samples from the Aegean Sea and Tunisian coasts were grouped together and the samples from Sardinia, Corsica, Elba Island, Sicily and the Venetian Lagoon were considered a unique separate group (W CT = 0.369, P,0.001) ( Table 2, A). However, a similar value of molecular variance (W CT = 0.367, P,0.001) was obtained when treating the Venetian Lagoon as a separate third group (Table 2, B). When alternative a priori groupings of samples were tested, based on biogeographic criteria, the AMOVA showed a decrease in the proportion of W CT variance ( Table 2, C and D).
b) Historical demography. The mismatch distribution carried out on the entire dataset showed a unimodal distribution of pairwise DNA differences ( Figure 4A), fitting the Rogers and Harpending's model [51] of demographic expansion (SSD = 0.002, P.0.05). When assuming the occurrence of discrete groups of populations that, according to AMOVA, represent three geographic areas of (1) Sardinia, Corsica, Elba Island and Sicily, (2) the Venetian Lagoon (Adriatic Sea), and (3) Aegean Sea and Tunisian coasts, we found that the two former groups displayed a unimodal mismatch distribution fitting the demographic expansion model (Figures 4B and 4C), whereas the mismatch distribution of the Aegean and Tunisian samples ( Figure 4D) did not correspond to the expected distribution under a demographic expansion model (SSD = 0.014, P,0.05). Neutrality tests showed a significant departure from equilibrium, with the exception of Tajima's D for the Venetian Lagoon (Adriatic) population ( Figure 4C). The negative values resulting from the neutrality tests were consistent with population expansions in all three groups. The non-significant value of D encountered in the Venetian Lagoon may reflect the reduced power of this test compared to the Fu's Fs to detect demographic expansions [58].
The time since expansion calculated according to Rogers and Harpending [51] for P. nobilis in Mediterranean was estimated to be the Pleistocene. In particular, for populations from Sardinia, Corsica, Elba Island and Sicily, we found that the demographic expansion started around 1.6 Mya (95% CI: 0.512-2.643 Mya). A similar estimate was calculated for the Venetian Lagoon population, with time since expansion dating back to about 1.3 Mya (95% CI: 0.643-1.956 Mya). A more recent demographic expansion was inferred for the Aegean-Tunisian group starting around 0.5 Mya (95% CI: 0.319-0.848 Mya).
The demographic history inferred by the Bayesian skyline plot revealed a signature of demographic expansion starting approximately 0.270 Mya and ending around 0.100 Mya ( Figure 5A). However, only samples from Sardinia, Corsica, Elba Island and Sicily showed a demographic expansion that corresponded to this time interval ( Figure 5B). Conversely, the Venetian Lagoon ( Figure 5C) and Aegean-Tunisian samples ( Figure 5D) had a virtually constant size throughout their demographic history.

Concatenated COI and 16S Dataset
The concatenated dataset comprised 338 bp for COI and 450 bp for 16S. The sequences of 219 specimens were determined here. Twenty-five COI and 16S sequences from Katsares et al. [35] were added. Among these 244 individuals, 103 haplotypes (Table S5), defined by 78 polymorphic sites (S) were found. Total mean haplotype diversity and nucleotide diversity were h = 0.961 and p = 0.005, respectively. A large number of haplotypes (77%) was found only at a single locality. The Aegean populations showed the lowest average values of haplotype diversity and nucleotide diversity (h = 0.333 and p = 0.001), whereas high levels of genetic diversity were found for almost all other sites. Estimates of genetic diversity for the concatenated COI and 16S dataset are provided in Table S2. Furthermore, estimates of genetic diversity and haplotype frequencies for the 16S region are given in Table  S2, and Table S6, respectively. a) Genetic structuring. The Bayesian analysis identified four haplotype groups, hereafter denominated N1, N2, N3 and N4 (Table S7). The most frequent group (N2), representing half of the individuals analysed, included the two most common haplotypes (PNCS 1, PNCS 6; Table S5), found in all populations, with a few exceptions in Sardinia (Isola di Mal di Ventre-IMV), Cyprus and the Aegean Sea (Epanomi-EP, Aggeloxori-AG, and Korinthiakos Gulf-KO). The second group (N3) (32.8% of individuals) included samples from Sardinia, Corsica, Elba Island, Sicily, the Venetian Lagoon and Cyprus. The groups N1 and N4 were less frequent, occurring in 9% and 8.2% of individuals only, respectively. N1 occurred in individuals from Sardinia, Corsica, Elba Island, Sicily and the Venetian Lagoon, while N4 included individuals from all Aegean populations with the exception of those from Xios-XI.
The haplogroups frequency distribution largely resembled that obtained by the COI analysis ( Figure 6). In particular, groups N1 and N3 approximately mirrored the distributions of P4 and P3, respectively. A striking exception is represented by haplogroup N2, whose geographic distribution overlapped that of the COI haplogroup P1, spreading as far as the Aegean Sea. At Xios-XI, the N2 group reached frequencies as high as those found in western Mediterranean and Ionian Sea populations. In contrast, the geographic distribution of haplogroup N4 was restricted to the Aegean Sea.
The median-joining network of haplotypes ( Figure S2 and Figure 7) showed a pattern of high genetic variability distributed across the Mediterranean. The two most frequent haplotypes (PNCS 1, PNCS 6; Table S5) were found in 23.8% of individuals and were shared only among individuals from Sardinia, Elba Island and Sicily. Different star-like branching patterns radiating out from the inferred root haplotypes occurred in the network ( Figure S2). A few point mutations (1 to 6) distinguished haplotypes from each other, and a diffused unresolved reticulation of branches was found among individuals. Haplotypes belonging to the Venetian Lagoon and Aegean Sea populations occupied a peripheral position in the network. The two Cyprian haplotypes were very divergent from one another ( Figure S2). No evidence of spatial structuring was apparent at either the local or large geographic scales. However, an exception was represented by Aegean individuals, which did not share haplotypes with other samples. Considering groupings detected by Bayesian assignment analysis (Figure 7), the network showed that the most diffused group (N2) was characterised by two main haplotypes and many derived ones, that diverged for a single point mutation. The group  Table 1. The number of mutations on the network branches are reported in the Figure S1. doi:10.1371/journal.pone.0067372.g003 N3 was also characterised by a high level of haplotype variation, with unique haplotypes radiating from the central ones. The group N1 exhibited a similar pattern, although a reduced level of polymorphism was found. The group N4 showed the lowest level of diversity.
The pairwise W ST estimates obtained for the concatenated dataset (Table S8) were similar to those reported for COI. A significant genetic differentiation occurred between the Aegean samples of Epanomi-EP and Aggeloyesori-AG, and the other samples. As previously noticed, the Venetian Lagoon sample was significantly differentiated from the others in most comparisons with the only exceptions being two Sardinian populations (Le Saline-SAL, Porto Ottiolu-OTT). Nonetheless, combining data from different genes resulted in a larger number of significant comparisons (Table S8). In particular, samples from Corsica appeared to be genetically differentiated not only from samples outside the Sardinian-Corsican region, but also from their northern Sardinian counterparts, and samples from Xios-XI were significantly different from those from the Aegean Sea.
Consistently with results of the COI analysis, AMOVA (Table 3) detected maximum genetic differentiation among groups defined a posteriori on the pairwise W ST values, when samples from the Aegean Sea were grouped separately from the others (W CT = 0.396, P,0.001) (Table 3, A). A similarly high level of molecular variance was found when the Venetian Lagoon was considered as a third group (W CT = 0.356, P,0.001) (Table 3, B). AMOVA tested on a priori groupings of samples also showed a decrease in the proportion of variance (Table 3, C and D); particularly no significant differences were found between groups corresponding to the western and eastern Mediterranean (W CT = 0.045, P.0.05) ( Table 3, Figure 8C, 8D) were likely to reflect the reduced power of the tests to detect departures from equilibrium, due to small sample size. Indeed, when we used only samples from Aegean Sea, Fu's Fs statistic did not displayed a significant departure from equilibrium for the COI dataset too (data not shown).

Discussion
This study provides the first analysis on mtDNA genetic variation of Pinna nobilis populations from a wide study area, roughly corresponding to the central part of the western Mediterranean, Ionian Sea and Adriatic Sea ecoregions (sensu Spalding et al. [30]). Supplementing the sequences obtained in the present study with those from previous investigations conducted in the Aegean Sea and along the Tunisian coasts provided a deep insight into (1) the large-scale patterns of spatial genetic variation of P. nobilis and (2) the role of Mediterranean biogeographic boundaries in shaping the species' genetic structuring. The high variability of the mtDNA COI region of P. nobilis further supports the effectiveness of this marker in investigating the genetic architecture of marine species (see [64] and references therein). In addition, the sample sizes adopted in the present survey allowed us to identify 16S haplotypes undetected in the previous study by Katsares et al. [35]. As the results obtained from the analysis of the COI sequences and of the merged COI and 16S sequences were largely consistent, we discuss them jointly in the following sections.

Genetic Structure of Pinna nobilis in the Mediterranean
Our study showed the occurrence of a moderate level of genetic structuring in P. nobilis populations in the study area. The low genetic divergence among the haplogroups and haplotypes, as We detected a primary, significant pattern of genetic structuring along the West-East direction, with samples from the Aegean Sea (excluding the Xios-XI sample for the concatenated COI and 16S dataset) and Tunisian coasts, which were genetically differentiated from the remaining samples. This pattern could be explained as a possible consequence of a reduced gene flow between the western and eastern Mediterranean. Several studies have shown that a hindrance to gene flow exists between the western and eastern Mediterranean, even for species with high dispersal capabilities and large population sizes (e.g., [65], [66], [67]). This result was partially consistent with findings reported for other bivalves with a high potential for dispersal, along with significant genetic differentiation and a pattern of isolation by distance (e.g., [68], [69], [70]).
Other studies have also shown that the genetic break between the western and eastern Mediterranean corresponds to the Siculo-Tunisian Strait (e.g., [71]); nonetheless, for different marine organisms the geographic position of the genetic break may shift eastward with respect to the strait (e.g., [29], [72] and references therein). In our case, the pattern of genetic structuring between the Aegean with Tunisian samples and other populations (Sardinian-Corsican region, Elba Island and Sicily) did not mirror the West-East separation of the Mediterranean along the Siculo-Tunisian Strait (see [73]). Moreover, the sample from the southern coast of the strait (Tunisian populations), which was homogeneous with the Aegean populations, was genetically divergent from the samples of the northern coast of the strait (Ionian Sicily), suggesting the occurrence of a further genetic break in a North-South direction. This finding could be due to the hydrodynamic characteristics of this area: the Algerian Current, which flows eastward from the Strait of Gibraltar, along the northern coasts of Africa and passes over the Siculo-Tunisian Strait [39], represents a hydrodynamic barrier between the facing coasts of Tunisia and Sicily. Such a barrier, which reduces the connectivity among populations via marine currents, has also been invoked for the identification of distinct mitochondrial lineages in the vermetid gastropod Dendropoma petraeum in adjacent locations in Tunisia and Sicily [74]. Furthermore, the genetic affinity of the Sicilian samples with the western Mediterranean ones, led us to hypothesise a shift of the boundary of the genetic break eastward from the Siculo-Tunisian   Strait, as proposed by Pérès and Picard [75] and Bianchi [29]. Although few studies on population genetics involve this area, a similar picture has been reported for Patella rustica [76], whereby the western mtDNA lineage reaches as far as the southern coast of the Italian Peninsula, well beyond the Siculo-Tunisian Strait. A genetic break located at the level of the southern tip of Calabria has also been described for Posidonia oceanica [77]. However, at this stage of knowledge, we cannot rule out that the genetic pattern observed around the western-eastern edge of the Mediterranean may reflect a primary genetic break in a North-South direction along the Siculo-Tunisian Strait. Furthermore, as has already been reported for other species [77], the genetic homogeneity found between the samples from the Tunisian coasts and the Aegean Sea could be ascribed to the persistence of a relict population in the eastern, semi-enclosed basin, rather than to an effective gene flow between the two Mediterranean basins. More individuals from unsampled regions of both western and eastern Mediterranean are needed to disentangle this point.
Our results also suggest the occurrence of a weak but significant genetic divergence of the Venetian Lagoon population from the above-mentioned two groups of populations. This finding is further evidence of the semi-enclosed nature of the Adriatic Sea, which represents a well-defined phylogeographic region within the Mediterranean [13]. A genetic discontinuity between the Adriatic Sea and the rest of Mediterranean has also been found in other species, such as the common cuttlefish Sepia officinalis [78], and the sea urchin Paracentrotus lividus [79]. Furthermore, several studies have revealed the occurrence of a genetic break among the Adriatic and Aegean seas, that was explained by the hydrographical isolation of these two basins [80], [81], [82], [83].
Interestingly, samples from Cyprus, at the easternmost part of the Mediterranean, exhibited haplotypes typical of the westernmost populations. The human-mediated introduction of planktonic larvae passively transported in the ballast water of commercial ships is a common finding in marine molluscs invasion [84], [85], [86], [87] and may account for this outcome. However, the small number of samples from Cyprus prevented any statistically supported inference.

Historical Demography
The results of historical demographic analysis and the estimates of time since expansion lead to the hypothesis of a Pleistocene scenario for P. nobilis, with expansions/contractions of its populations due to sea-level fluctuations of this period, repeatedly leading to the isolation, or partial isolation, of the Mediterranean basins [88]. Conversely, vicariant events, leading to ancient divergence of populations, are unlikely to have occurred in P. nobilis because we would expect a bimodal mismatch distribution of pairwise DNA differences for the entire dataset, which may occur when different genetic lineages are present [11], [89]. It should be noted that the estimates of time since expansion may be biased by the lack of a mutation rate specific for P. nobilis and the time dependency of molecular rate estimates [90]. In the present investigation, we used a range of inter-specific divergence rates estimated for bivalves [59]. However, the estimated molecular rates were much higher over short timescales at the intra-specific level than between species ( [90], [91] and references therein but see also [92]). Despite these caveats, our results, nonetheless, provide initial insights into the dynamics of the colonisation of P. nobilis populations in the Mediterranean.
All populations analysed in this study are likely to have shared a common origin. Our data suggest an eastward expansion across the Mediterranean, likely mediated by marine currents (e.g., the Algerian Current), and followed by one or more founder events, that lead to lower diversity and private haplotypes in the Aegean and Tunisian populations. A similar trend was also reported for populations of the sea cucumber Holothuria mammata from the

Final Remarks
Rabaoui et al. [36] suggested that the conservation and management of P. nobilis should benefit from improvements in the sampling plan (with the inclusion of samples from the western Mediterranean). Although part of the range of the species remains unsampled, our study represents an effort in this direction, as it triples the existing genetic information on P. nobilis with new data from a wide area of the western Mediterranean.
The general picture revealed by our study is the presence of a large, homogeneous group of populations of P. nobilis spanning two Mediterranean marine ecoregions (the western Mediterranean and Ionian Sea), which genetically diverges from the Adriatic population and those from the Aegean Sea and Tunisian coasts. The semi-enclosed nature of the Adriatic Sea may explain the genetic divergence of the Venetian Lagoon population. Hindrances to gene flow, that are related to biogeographic boundaries among these Mediterranean sectors, and the Pleistocene changes in sea levels, can be invoked to explain the pattern of genetic structuring shown in this study. Our results further suggest the occurrence of a West-East genetic break located eastward of the Siculo-Tunisian Strait, and also provide evidence of a further genetic break in this area which is in the North-South direction. In this context, it would be worth noting that other samples from the western (e.g. Spain and Balearic Islands) and eastern Mediterranean (e.g. from Egypt and Israel) would be very informative for the understanding of the species' evolutionary history.
Interestingly, although the data for the local or regional levels of genetic erosion due either to human disturbance or to bioecological changes are unavailable, the higher levels of mitochondrial variability found in the Sardinian-Corsican region, Elba Island, Sicily and the Venetian Lagoon relative to the populations from the Aegean Sea and Tunisian coasts support the occurrence of an eastward demographic expansion. The lack of strong genetic divergence among these three geographic areas (e.g., only one point mutation separates the Aegean and Tunisian samples from those from Sardinia, Corsica, Elba Island, Sicily and the Venetian Lagoon) is evidence of a common origin of P. nobilis populations in the Mediterranean.
From a conservation point of view, the three genetically divergent groups, (1) the samples from Sardinia, Corsica, Elba Island and Sicily, (2) the samples from the Aegean Sea and Tunisian coasts and (3) the Venetian Lagoon sample, should be considered as different management units. Moreover, the samples from the Venetian Lagoon may be considered a ''peripheral isolate'' sensu Frey [93], with particular conservation relevance, given its peculiar geographic position.