Genetic Structuring across Marine Biogeographic Boundaries in Rocky Shore Invertebrates

Biogeography investigates spatial patterns of species distribution. Discontinuities in species distribution are identified as boundaries between biogeographic areas. Do these boundaries affect genetic connectivity? To address this question, a multifactorial hierarchical sampling design, across three of the major marine biogeographic boundaries in the central Mediterranean Sea (Ligurian-Tyrrhenian, Tyrrhenian-Ionian and Ionian-Adriatic) was carried out. Mitochondrial COI sequence polymorphism of seven species of Mediterranean benthic invertebrates was analysed. Two species showed significant genetic structure across the Tyrrhenian-Ionian boundary, as well as two other species across the Ionian Sea, a previously unknown phylogeographic barrier. The hypothesized barrier in the Ligurian-Tyrrhenian cannot be detected in the genetic structure of the investigated species. Connectivity patterns across species at distances up to 800 km apart confirmed that estimates of pelagic larval dispersal were poor predictors of the genetic structure. The detected genetic discontinuities seem more related to the effect of past historical events, though maintained by present day oceanographic processes. Multivariate statistical tools were used to test the consistency of the patterns across species, providing a conceptual framework for across-species barrier locations and strengths. Additional sequences retrieved from public databases supported our findings. Heterogeneity of phylogeographic patterns shown by the 7 investigated species is relevant to the understanding of the genetic diversity, and carry implications for conservation biology.


Introduction
Through the study of species distribution as well as patterns of genetic structuring, a link between biogeography and phylogeography can be explored, which aids in the identification of discontinuities in the spatial distribution and genetic diversity of species [1]. The level of connectivity between geographic areas is mainly determined by the interplay between the dispersal capacity of a species (including pre and post-settlement processes [2]), the hydrological regime [3] and the occurrence of geological/ topographical boundaries, which may reduce the dispersal potential and hence the gene flow [4].
Major barriers to gene flow between oceanic basins have previously been identified (for example, Point Conception in California, [5]; the Eastern Pacific barrier, [6]; the Indo-Pacific Barrier, [7]; the South-eastern Australian barrier, [8]). These barriers support the occurrence of genetic discontinuities between the world's biogeographical regions [9]. Significant genetic structuring has also been detected at smaller spatial scales in correspondence with the boundaries between biogeographic provinces [10] or sectors [11,12].
In the Central Mediterranean Sea nine biogeographic sectors have been proposed based on species' distributions [13], but their consistency with phylogeographic sectors has only been investigated in a few cases [14][15][16]. These studies investigated genetic differentiation between the Mediterranean basins (eastern, west-ern, Adriatic Sea) in a single species. However, data from a single species, with specific life history traits, ecology, ethology, etc, cannot be used to generalize phylogeographic barriers or limitations to gene flow across other species. In fact, studies on different target species conducted in the same geographical areas often reveal contrasting genetic structures and/or phylogeographical patterns [17][18][19]. Comparative phylogeography has recently been suggested as a powerful tool to discern general patterns over species-specific traits [20,21]. It has proven effective in identifying across taxa barriers to gene flow in the Almeria-Oran Front [22], and along the north-eastern Pacific coast [23]. Currently no multispecies studies identifying barriers to gene flow or genetic discontinuities in the Mediterranean Sea exist. Given the practical implications of defining genetic units for the conservation and management of natural resources, a study of this nature is sorely needed [24].
In the present study, the occurrence of phylogeographic boundaries in the Central Mediterranean Sea, and their consistency across species was investigated. By means of a multifactorial hierarchical sampling design, three putative long-term barriers to gene flow were tested. These putative barriers were established based on known biogeographic boundaries [13]. Two replicated areas were sampled across each barrier, and within each area the genetic structure of seven common shallow water invertebrate species was analysed using a portion of the COI mitochondrial gene.
In addition to genetic differentiation measures, two types of multivariate analysis were used to test the consistency of the phylogeographical patterns across species. The results of this study provide a first contribution towards a comparative phylogeography in the Central Mediterranean Sea, with the aim of fully understanding the main past and present drivers of genetic diversity and differentiation.

Study areas and species
Three major potential biogeographic/phylogeographic boundaries, located in the Central Mediterranean Sea (Figure 1), and hereafter referred to as putative barriers (PB), were investigated. The first putative barrier (PB1, in Location 1) separates the Ligurian Sea and the North Tyrrhenian Sea, corresponding to a virtual line between the coast of Tuscany, and the islands of Elba and Corsica. The two sides of PB1 are characterized by divergent currents [25]. The second putative barrier (PB2, in Location 2) separates the Tyrrhenian and Ionian Seas, and corresponds to a line across the Strait of Messina and Sicily. Very strong tidal currents characterize the narrow Strait of Messina, flowing mainly from north to south [25]. The third putative barrier (PB3, in Location 3) separates the Ionian and Adriatic Seas, around the Apulian Cape of Santa Maria di Leuca. Strong currents with different patterns in opposite directions characterize the water circulation along this putative barrier [26].
Seven invertebrate species were selected based on their prevalence in the central Mediterranean as well as to cover a variety of life history traits. The selected species belong to different taxa: Gastropoda (3), Polyplacophora (1), Ascidiacea (1); Crustacea (1), and Porifera (1). Data retrieved on life history traits and larval ecology are summarized in Figure 2. The longest Pelagic Larval Duration (PLD) estimates correspond to Balanus perforatus, Patella caerulea, and Halocynthia papillosa. Hexaplex trunculus is a brooding species and has no theoretical PLD. PLD estimates for Chondrosia reniformis, Chiton olivaceus, and Osilinus turbinatus are of about one week.
In each area 30 individuals per species were collected, fixed in 80% ethanol and maintained at 4uC prior to processing. Only the sponge was sampled as pieces rather than as whole organism. To minimize the possibility of clonality, samples were collected at least 5 m apart [27]. No specific permits were required for the areas and the species studied.

DNA extraction and COI amplification
Total DNA was extracted from a piece of circa 5 g of tissue using a REDExtract-N-Amp kit (Sigma-Aldrich). The COI gene was chosen for analysis due to its intermediate rate of variation that provides information on genetic structure across a wide range of spatial and temporal scales [28]. Universal primers described in [29] were initially used for the amplification of a fragment of the COI gene. After the first sequences were retrieved, six pairs of species-specific primers were designed with the software PRIMER v. 3.0 (http://www.fokker.wi.mit.edu/primer3/input.htm) in order to increase amplification success and the quality of sequences obtained (Table S1).
PCR amplification reactions were carried out in a 25 ml total volume with 5 ml of 5x GoTaq Flexi Buffer (Promega), 4 ml of MgCl 2 25 mM, 0.5 ml of dNTP mix 10 mM, 0.5 ml of each primer, 0.15 ml of GoTaq DNA polymerase (Promega) and 2.5 ml of template DNA (from a 1:50 dilution of the extracted DNA). A hold at 94uC for 3 min was followed by 30 cycles of denaturation at 94uC for 45 s, annealing at a primer specific temperature for 45 s and extension at 72uC for 90 s, finishing with a final extension at 72uC for 5 min on a thermal cycler (Applied Biosystems, GeneAmp 2700) (Table S1). PCR products were purified and sequenced by MACROGEN INC.

Data analysis
Preliminary analysis showed no differentiation in any of the investigated species between the two replicated areas of each site. As a result, individuals from the two areas were pooled and hereafter data have been analysed considering only the sites (A, B, C, D, E, F; Table 1). We detected no genetic differentiation between Locations 1 and 3 of Balanus perforatus, 1500 km apart, so decided not to analyse Location 2 samples.

a. Genetic diversity, differentiation and gene flow estimates
Number of haplotypes (h), haplotype and nucleotide diversity (Hd and p, respectively), were calculated for each species by site, location, and for the whole study area, using DnaSP v. 4.50.3 [31]. Diversity values and relationships between the haplotypes were   [32]. Loops were solved when possible following criteria derived from coalescent theory [33]. The fixation index (Fst) between sites across each PB was used to quantify the differentiation due to genetic structure. Large scale genetic structuring across the whole study area was also assessed, as the genetic differentiation between all pairs of sites. Pairwise Fst values were calculated in Arlequin v 3.1 [34] with 10000 permutations, and their significance corrected for multiple testing following [35].
Analysis of Molecular Variance (AMOVA) was performed to investigate sources of variability within and between the three locations studied. Significance was tested with 16000 permutations in Arlequin v 3.1 [34].
Isolation by distance across the whole study area was tested with a Mantel test as implemented in Genepop, following the Isolde procedure which compares a semi matrix of Fst/(1-Fst) with a semi matrix of the natural logarithm of the shortest oceanic distance in km.
Migration rates between areas were estimated with the software MIGRATE-n v. 3.3.1 [36], using Markov chain parameters.
COI sequences of the target species available in Genbank (Table S2) were aligned and represented in haplotype networks with the haplotypes obtained in this study, to provide some insight on large-scale genetic diversity of the species and patterns across other acknowledged barriers to gene flow.

b. Estimation of effective population size variations
The history of effective population size (Ne) was assessed by testing the fit of the mismatch distributions to the theoretical distribution in an expansion scenario [37] using Arlequin. Fu's and Tajima's neutrality tests were also carried out, as they are more sensitive to population expansions than the mismatch distributions [38]. These tests were performed for each site, for the whole data set, and for undifferentiated groups of sites identified by the previous population differentiation analysis.
For those sites that were not significantly different, and did not deviate from the expansion model, the coalescence time t (i.e. time since the start of a population expansion) was estimated following the formula T = 2ut [37], where T is the date of growth or decline in units of mutational time and u = 2 mk, where k is the number of generations per year and m is the mutation rate per nucleotide.
For the same undifferentiated sites, past changes in effective population sizes and estimates of expansion time were also inferred by generating Bayesian Skyline Plots (BSPs) using the software BEAST v. 1.7.5 [39]. The best-fit models of nucleotide substitution for each species data set were selected with jModeltest v. 0.1.1 [40] and fed into BEAST. Mutation rate estimates used were: 2.4% MY for molluscs [41], 3.1% MY for barnacles [42], and 2.86% MY for ascidians [43]. Given the absence of molecular clock estimates for sponges, the BSP for this species was not calculated. Confidence intervals for Ne through time were calculated from the posterior probability distributions using TRACER [44].

c. Multispecies analysis
To test for consistency of the phylogeographic patterns across species, two types of multivariate analysis were used.
The multivariate Analysis of Similarity (ANOSIM), as implemented in Primer v.6 [45], was used to test the differences across each of the three putative barriers. With this analysis we aimed to provide an across-species analysis of haplotype frequencies across the putative barriers, reassuming most of the information on the genetic structure of all species, and minimizing the effect of between species sample size and/or genetic diversity differences. The frequencies of each haplotype of all seven species were the variables, and the Areas were the samples. In order to accomplish with the independence of the variables, the haplotype frequencies were calculated for each haplotype across all sites (instead of for each site across all haplotypes). We used Site and Location as crossed factors in order to test for the effect of each of the PB. A Bray Curtis similarity matrix between all sites was represented graphically using a Multidimensional Scaling plot (MDS). Apparent groupings across the whole study area were also statistically tested with ANOSIM. Discriminant analysis of principal components (DAPC), available in the Adegenet package [46] for R (R core team 2012), was also performed. This technique is designed for multivariate genetic data (multi-locus genotypes). It maximizes the variation between groups by first performing a principal component analysis (PCA) on pre-defined groups or populations and then uses the PCA factors as variables for a discriminant analysis (DA), thus ensuring their independence. To convert our sequences into a multivariate dataset, each polymorphic position of each sequence was codified with a four-digit code, corresponding to the four nucleotides. Only the polymorphic sites were considered, creating a ''pseudogenotype'' of each specimen. Polymorphic sites across species were counted as null alleles in order to test all species simultaneously. In this case, the number of variables corresponds to the total number of specimens and the samples correspond to the six sites. This analysis allowed us to take into account not only the haplotype frequencies differences, but also how different these haplotypes are, minimizing the effect of the species (although unavoidable due to the different levels of polymorphism, and hence the number of ''pseudo-loci'' included per species). Patella caerulea and Hexaplex trunculus showed two haplogroups in their networks. P. caerulea presented two common haplotypes separated by a single mutation and surrounded by several low frequency haplotypes, 26 of which were private (Figure 3a). In H. trunculus, nine mutation steps separated the haplogroups, which could therefore be considered as two divergent lineages. The most common lineage occurred across all sites, whereas the second lineage was found in sites D, E, and F (except for H_27, Figure 3b). Two other highly divergent haplotypes were found in site A (H_16) and site E (H_23).

Genetic diversity and connectivity patterns
Three species presented star-like shaped haplotype networks with different degrees of ramification. The Balanus perforatus network (Figure 3g) showed a high number of low frequency haplotypes. Osilinus turbinatus also showed a profusely branched network, with several middle frequency haplotypes (H_1, H_11, H_13, H_31; Figure 3c). The Chiton olivaceous network had a simple star shape: one very common haplotype surrounded by several   private and three low-frequency haplotypes, with a maximum distance of three mutation steps (Figure 3e). The five haplotypes of Chondrosia reniformis (one common haplotype, three privates, and a middle frequency haplotype in sites E and F, Figure 3d) were separated by a maximum of 2 mutation steps. Halocynthia papillosa presented several closely related common haplotypes, and few private haplotypes (11 out of 21, Figure 3f).
Fst values failed to detect significant genetic differentiation across the putative barriers PB1 and PB3 for any of the seven species studied (Table 2). Significant genetic differentiation (p, 0.015 after FDR correction) was only detected across the PB2 for Patella caerulea and Hexaplex trunculus (Table 2a- Chiton olivaceus showed significant genetic differentiation between sites A, B, C and site F. Site E was also significantly different from site A (Table 2e). A significant Mantel test (p = 0.014) suggests that this pattern of genetic differentiation could be attributed to isolation by distance. Migration estimates were higher between nearby than between distant locations (M 1-. 2  Halocynthia papillosa showed significant genetic differentiation between several pairs of sites (Table 2f,) with no relation to distance (Mantel test p = 0.158). AMOVA analysis on the three locations wasn't significant (p = 0.064; Table S3), however it attributed a high percentage of variation to the differences between them (Va% = 13.56). Migration rate estimates showed wide confidence intervals (data not shown). The acorn barnacle Balanus perforatus showed no genetic differentiation between sites across putative barriers (PB1 and PB3), or even across sites located more than 1000 km apart (Table 2g).
Among the 32 sequences of Patella caerulea retrieved from Genbank (Table S2), 12 corresponded to the H_4 haplotype (7 from the Western Mediterranean, 2 from Lampedusa, and 3 from west of the Almeria-Oran Front) and six sequences corresponded to the H_1 haplotype (from Greece, Lampedusa, Tuscany, Valencia, and 2 from west AOF). The haplotypes H_11, H_26 and H_17 (Figure 3a, northwestern haplogroup) were found in 6 individuals from the western Mediterranean region, and 4 haplotypes were not found in our dataset. The 11 COI sequences of Hexaplex trunculus retrieved belonged to individuals from the South of Portugal (Table S2). Nine of them corresponded to the H_4 haplotype and the other two haplotypes were one nucleotide different from H_4.
Among the 10 Osilinus turbinatus COI sequences retrieved (Table  S2), 3 (2 from Sardinia and one from SE Spain) presented the H_4 haplotype (Fig. 3c); 3 corresponded to H_9 (Turkey, Cyprus and Croatia); and 2 belonged to H_24 (SE Spain). Two individuals from Croatia and one from Cyprus presented private haplotypes.
The single sequence from Chondrosia reniformis published in Genbank corresponded to the H_2 haplotype and belonged to an individual from Israel (Fig. 3d). The Chiton olivaceus sequence came from an individual from north-eastern Spain and presented the H_2 haplotype (Fig 3e). The six sequences of Halocynthia papillosa (Table S2) belonged to individuals from the north-western region of the Mediterranean. Two of these presented the H_12 haplotype, while the other two presented the H_13 haplotype and the last two showed the H_18 haplotype (Fig. 3f).

Estimation of population size variations
In four species, mismatch distributions and neutrality tests indicated expansion in most sites, as well as in the whole study area ( Figure S1; Table S4). The brooding gastropod Hexaplex trunculus did not show signs of expansion in most sites (Table S4b), probably due to the admixture of two genetic pools, visible in the mismatch distribution ( Figure S1b) and previously detected as the divergent lineages in the haplotype network. The solitary ascidian Halocynthia papillosa and the sponge Chondrosia reniformis did not show signs of expansion (Table S4).
Bayesian skyline plots (BSP) showed large confidence intervals in most cases ( Figure S2), yet they roughly reflected the estimates of time since expansion calculated with Rogers & Harpending's formula. The northern haplogroup of Patella caerulea and the southern haplogroup of Hexaplex trunculus demonstrated stable

Multispecies analysis
A matrix of 159 different variables corresponding to the frequencies of the haplotypes of all species (excluding Balanus perforatus) across sampling areas was generated. The analysis of similarity (ANOSIM) indicated no significant effect of the factor site crossed by location, i.e. no multispecies significant differences across any of the three putative barriers. The MDS of the Bray-Curtis matrix of distances between all six sites (Table 3) showed higher similarity of sites A, B (Location 1) and C (Location 2). Site D (Location 2) was isolated, and sites E and F (Location 3) clustered together (Figure 4a). ANOSIM on these three groups showed significant differences between them (R = 1; p = 0.017).
After transforming haplotypes into pseudo-genotypes, we obtained a matrix of 802 individuals with genotypes of 152 loci each (most of them zeros), corresponding to all polymorphic sites across all species (33 of Patella caerulea, 50 of Hexaplex trunculus, 8 of Halocynthia papillosa, 37 of Osilinus turbinatus, 19 of Chiton olivaceus, and 4 of Chondrosia reniformis). The resulting DAPC showed a distribution of the samples (sites) similar to the MDS plot: sites A, B and C grouped together. Site D was separated from all the others, and sites E and F clustered together (Figure 4b).

Discussion
In the study area two zones of genetic discontinuity have been identified: across the biogeographic boundary between the Tyrrhenian and Ionian Seas, and between the western and eastern coasts of the Ionian Sea.
The strongest genetic structure detected was across the Tyrrhenian-Ionian biogeographic boundary. The presence of this long-term boundary was reflected in the genetic structure of two species out of the seven studied. It probably corresponds to the boundary between biogeographic sectors 3 and 6 as defined by [13] (see Figure 1). The genetic structure across these sectors has been described by [47][48][49], but the location of the boundary remained vague. In this area, due to the complex topography of straits and islands, Pleistocene sea level fluctuations have left a signature in the genetic structure of many species [48], subjected to periodic divergences and reconnections of the populations on the two sides of the Messina Strait. Past history events might explain the origin of haplogroups or divergent lineages, whereas present day events, such as the current regimes, aid in maintaining and enhancing the differences due to limited or asymmetric connectivity between sites.
Patella caerulea showed significant genetic structure across this Tyrrhenian-Ionian PB, reflected in the two haplogroups with different frequencies across each side of the boundary. In spite of the significant differentiation, a bidirectional but asymmetric migration across the barrier was estimated. Other limpet species showed comparable genetic structure in this area [50]. In addition, the sequences retrieved in Genbank gave further support to this structure: western and eastern Mediterranean samples clustered more often with haplotypes of the ''northern haplogroup'', and ''southern haplogroup'', respectively. It is noteworthy that individuals from the Atlantic side of the Almeria-Oran Front (AOF), the best-known barrier to gene flow in the Mediterranean [20], presented common haplotypes, suggesting this barrier had no large effect on P. caerulea genetic structure [51]. Expansion estimates indicate long-term persistence of this species and, in accordance with the older expansion estimates, an older origin of the northwestern haplogroup. This is in agreement with the widespread distribution and abundance of this species and with the more probable Atlantic origin of this genus [52].
The subtidal gastropod Hexaplex trunculus, a brooding species, showed the sharpest genetic structure across the Tyrrhenian-Ionian biogeographic boundary. Despite this fact, several haplotypes were common to sites located up to 1500 km apart, on both sides of the putative barrier, probably indicating high migration rates from the Ligurian-Tyrrhenian towards the Ionian-Adriatic populations, but almost none in the opposite direction. This sharp asymmetry has been described in other regions and species, and is suggested to reflect persistent current regimes affecting larval dispersal during the reproductive period [53]. Estimates of expansion of both lineages are consistent with the persistence of this species during the LGM, in more than one glacial refuge, followed by a unidirectional secondary contact [54]. Expansion estimates are older in the southern lineage, which, together with the higher diversity values in this area, indicates a major contribution of the Eastern Mediterranean genetic pool to the total genetic diversity of this species [41]. All sequences retrieved from Genbank belonged to individuals from South Portugal, so no further details on the structure of other Mediterranean sites could be inferred. Most individuals showed haplotypes of the northern haplogroup, reinforcing their Mediterranean origin, but also suggesting no effect of the AOF on the long-term genetic structure of this species, as was the case for Patella caerulea. An unexpected area of genetic discontinuity was disclosed between the eastern and western Ionian Sea. Two species, presumably not affected by the Tyrrhenian-Ionian boundary previously described, showed significant genetic structure across the Ionian Sea, i.e. between the Adriatic-Ionian sites and the rest of the study area. The Ionian Sea has never been considered to host phylogeographic barriers, either related to historical events nor according to present day oceanographic processes. No records referring to genetic structure across this region were found in the literature, despite the fact that the genetic structures of some species support this finding [48], although no explicit tests to identify the barrier had been carried out. The location of the biogeographic boundary separating the eastern and western Mediterranean is controversial [55]. This newly detected phylogeographic boundary provides new insight to understanding the biogeography of the area, although further data on the genetic structure of other species in the area would be needed.
Osilinus turbinatus presented a high diversity of haplotypes and significant structure across the Ionian Sea, with bidirectional and symmetric migration estimates. Published sequences supported the detected structure and individuals from the eastern Mediterranean and the Adriatic Sea presented haplotypes related to those found exclusively in the Adriatic-Ionian sites. The estimates of expansion of this species coincide with the Last Glacial Maximum (LGM between 24,000 and 27,000 ya, [20]).
The sponge Chondrosia reniformis, in spite of evincing low diversity values, also showed significant genetic structure across the Ionian Sea. A similar differentiation has been described in other widespread sponge species in the same geographical area [56]. The only sequence retrieved from Genbank belongs to an eastern Mediterranean individual and shares the more common haplotype. Although further data from the Ionian Sea and the eastern Mediterranean would be needed to confirm our finding, this study suggests that for Osilinus turbinatus and C. reniformis the phylogeographic barrier between eastern and western Mediterranean is located in the Mid Ionian Sea.
Significant genetic structure was detected across the northern and southern Tyrrhenian Sea in Hexaplex trunculus and Osilinus turbinatus. In both cases significant Fst values were likely due to the occurrence of several private haplotypes, south Tyrrhenian haplotypes in the case of H. trunculus, and Ligurian-North Tyrrhenian private haplotypes in the case of O. turbinatus.
Halocynthia papillosa showed several significant pairwise differences between sites across the three locations. The positive values of the neutrality tests, the presence of several common haplotypes, the smaller number of private ones, and the flat shape of the BSP, strongly indicate a long-term stability of the population size. This species may have persisted in two or more glacial refugia, followed by a secondary contact between them [57]. Although this ascidian has large dispersal capacity estimates [43] our results support larval retention, as has been previously suggested ( [58], Figure 2), allowing genetic structuring in the range of 300 and 1500 km.
Balanus perforatus, the species with the longest PLD estimates, lacked significant genetic structure across distances of up to 1500 km. Neutrality tests indicated a very old expansion, starting around 60,000 ya. This pattern of persistence is similar to that shown by other widespread and abundant species inhabiting comparable habitats along the northeastern Pacific [59,60].
The genetic structure found among populations of Chiton olivaceus across the whole study area is congruent with a pattern of isolation by distance. This species showed the only estimate of expansion time clearly posterior to the LGM, indicating reduced persistence and recent colonization from glacial refugia [60]. This would also be in agreement with the low abundances and low molecular diversities detected with respect to the other species studied, and to other polyplacophora species [3].
As previously shown [61], our results support the weak correlation between PLD and genetic structure. The scale of larval dispersal is related to the complex interaction of larval life history traits, oceanographic regimes, habitat quality and distribution, and the variability of each of these factors through time. However, it is noteworthy that contrasting dispersal capacity patterns, as those of Hexaplex trunculus and Balanus perforatus, influence connectivity across the boundaries of biogeographic areas, as already detected in fishes [62]. The genetic structure observed is most likely due to the different effects of glacial and interglacial isolations and reconnections on their populations. Levels of gene flow are not the only factors determining present day genetic structure among populations; divergence times and population size can greatly affect the observed genetic structure [63], however the lack of basic biological data for most of our focus species, does not allow to disentangle the role these different drivers play in the observed genetic structures.
The two multivariate analyses allowed a general picture to be provided by merging in a common framework the results obtained for individual species. The statistical tests as well as the graphical representation provided a general picture of the phylogeographic structure in the Central Mediterranean, as well as the relative importance of the discontinuity areas detected, according to the number of species affected. This approach enables to summarize the available genetic information in a simple, straightforward result, which might be useful for decision-making or management purposes. Species with high haplotype and/or nucleotide diversities could bias these analyses; however the multispecies approach counteracts this effect, so the results are not obscured by the structure of the most variable species.
According to the results of the present study, the Tyrrhenian-Ionian biogeographical boundary acts as a phylogeographical barrier, as suggested by the genetic structure of two out of seven invertebrate species. The hierarchical sampling design adopted allowed the identification of an area of genetic discontinuity in the mid-Ionian Sea, reflected in the genetic structure of two species. The data retrieved from public databases reinforced the observed haplotype distributions, and gave further evidence on the variability of the effect of known barriers to gene flow on the genetic structure of different species. In the Mediterranean Sea, Pleistocene glaciations had a profound influence on the present day distribution of species and genetic diversity. The estimates of expansion and divergence allowed further understanding of the genetic structure found at present. Due to the ecological importance of the Central Mediterranean Sea as a hot spot of marine biodiversity [64], these results are of importance for the better understanding of the ecological functioning of its coastal ecosystem. General patterns pertinent to the management and conservation of this important ecosystem have been shown thanks to the comparative phylogeography approach used. Further research analyzing the genetic structure of other species, inclusive of additional Putative Boundaries, while also utilizing a larger set of molecular markers, would provide more insight on a multispecies, multi-scale picture of connectivity in the marine realm. Table S1 Specific primers. Names, sequence, amplified fragment length, and optimal annealing temperature for each of the seven pairs of primers used. (DOCX) Table S2 Previously published sequences. Accession numbers, geographical origin and identification with the haplotypes found in this study, as named in Figure 3, of sequences retrieved from Genbank. (DOCX)