This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Introduction

We examined patterns of genetic divergence in 26 Mediterranean populations of the
semi-terrestrial beachflea Orchestia montagui using mitochondrial (cytochrome oxidase subunit I), microsatellite (eight loci) and
allozymic data. The species typically forms large populations within heaps of dead
seagrass leaves stranded on beaches at the waterfront. We adopted a hierarchical geographic
sampling to unravel population structure in a species living at the sea-land transition
and, hence, likely subjected to dramatically contrasting forces.

Results

Mitochondrial DNA showed historical phylogeographic breaks among Adriatic, Ionian
and the remaining basins (Tyrrhenian, Western and Eastern Mediterranean Sea) likely
caused by the geological and climatic changes of the Pleistocene. Microsatellites
(and to a lesser extent allozymes) detected a further subdivision between and within
the Western Mediterranean and the Tyrrhenian Sea due to present-day processes. A pattern
of isolation by distance was not detected in any of the analyzed data set.

Conclusions

We conclude that the population structure of O. montagui is the result of the interplay of two contrasting forces that act on the species
population genetic structure. On one hand, the species semi-terrestrial life style
would tend to determine the onset of local differences. On the other hand, these differences
are partially counter-balanced by passive movements of migrants via rafting on heaps
of dead seagrass leaves across sites by sea surface currents. Approximate Bayesian
Computations support dispersal at sea as prevalent over terrestrial regionalism.

Keywords:

Introduction

In many marine species genetic structure is often shallow as a consequence of large
population sizes coupled with high potential for long distance dispersal. A high number
of polymorphic markers, large sample sizes and ad-hoc designed sampling schemes are hence necessary to unveil population structure [1]. On the other hand, phylogeographic patterns with a clear geographic component are
relatively common and easy to detect in terrestrial and/or freshwater species due
to the intimately fragmented nature of their habitats [2]. Phylogeographic patterns in species living in the transition area between sea and
land (supralittoral zone) have received relatively little attention compared to the
wealth of work done on marine and terrestrial species [3-6]. However, a proper understanding of the evolutionary trajectories governing coastal
communities is fundamental, especially in light of the increasing pressure that human
beings have been exerting on them in the last decades [7].

Talitrid amphipods are often the dominant animal group in the supralittoral zone,
where they are confined to a narrow area avoiding both the arid environment inland
and submersion in seawater. Active dispersal is restricted to short, mostly nocturnal,
foraging movements along the beach [8,9]; maximal active distance recorded for talitrids in nature is 200 m [10]. Talitrids lack any planktonic larval stage and consequently they cannot rely on
this stage for long-distance passive dispersal. They are, however, considered as facultative
rafters [11,12]. Generally speaking, rafting involves the traveling of organisms over the sea surface
via floating material. Talitrids are often associated to wrack, wood or other floating
material that could be washed away by waves and transported by sea surface waters
from one beach to another. Rafters can be transported from short to medium-long distances
and arrive in otherwise unreachable areas [12]. This process has a special impact on those species that are not able to actively
disperse over long distances. In particular, it could potentially increase connectivity
among populations, thus reducing the onset of local genetic divergence. Peracarid
crustaceans are found in both high diversity and abundance on floating items, which
seemingly play an important role in their life histories [12].

The semi-enclosed Mediterranean Sea is surrounded by thousands of kilometers of coasts.
Previous genetic studies focused on multiple Mediterranean populations of seven different
species of talitrids [8,13-18] have revealed two main patterns of genetic structuring, corresponding to two groups
of species with contrasting life histories. Those species (i.e. Orchestia montagui, O. mediterranea, Platorchestia platensis) that live within beached decaying seagrass, wrack or gravel, showed a shallow genetic
structuring even on a wide geographical scale. Species (i.e. Talitrus saltator, O. stephenseni, Macarorchestia remyi) that are exclusively associated to sandy beaches or stranded rotting logs and have
the tendency to burrow into the substrate were invariably highly structured genetically
[17,18]. De Matthaeis et al. [8,15] interpreted these alternate patterns in terms of probabilities for talitrids to realize
passive dispersal via rafting by being dragged away by waves. These probabilities
would be high for the species belonging to the first ecological group because they
colonize ephemeral habitats close to the waterline. On the other hand, they would
be low for species of the second group bound to sand beaches or rotting logs; these
species colonize the upper supralittoral zone and are not well suited to survive on
rafts [19].

Here we focus on Orchestia montagui Audouin 1826 and we expand the data available on the species both in terms of geographic
sampling and kind of molecular markers used. The species is abundant along the coasts
of the Mediterranean Sea and reaches the Black Sea on the East ([8,20]; Pavesi pers. obs.); it is commonly found in large numbers within banks of beached
dead leaves of Posidonia oceanica, a seagrass endemic to the area. A number of previous studies, all based on allozymes,
considered mostly populations from islets surrounding the island of Sardinia (Western
Mediterranean and Tyrrhenian Sea), where O. montagui is the most common talitrid, and a few populations from the Eastern Mediterranean
and Aegean Sea [8,9,13,21]. At that geographic scale and for allozymes, the species resulted poorly structured.
We recently developed eight highly polymorphic microsatellite loci for O. montagui[22]. We used these loci, together with a fragment of the mitochondrial DNA (mtDNA) coding
for the cytochrome oxidase subunit I (COI), to coarsely assess relationships among
six O. montagui populations from five Mediterranean basins (Tyrrhenian, Ionian, Adriatic, Western
and Eastern) [23]. In that study both classes of markers congruently retrieved a major phylogeographic
break between Ionian and Adriatic populations on one side and Tyrrhenian, Western
and Eastern Mediterranean populations on the other side, thus capturing a pattern
of divergence that escaped allozymes. These results suggested that the lack of structure
initially found with allozymes was probably due to their relatively poor resolving
power. We planned the present study to overcome the sampling limitations of [23] and to place for the first time under the same analytical frame previous allozyme
data along with new mtDNA and microsatellite information collected for the largest
number of sampling localities we could access along the coasts of the Mediterranean
Sea. We were able to expand our samplings to a total of 26 populations, and thus attain
a more detailed picture of patterns of genetic connectivity in O. montagui. These populations were sampled in six different Mediterranean basins (see details
in Table 1 and Figure 1) and screened for variation at mtDNA and microsatellites. To contrast genetic divergence
at different geographical scales, we included in the study populations that were separated
from a very few to hundreds of kilometers. In order to make possible combining the
data acquired for this study with those previously published in De Matthaeis et al.
[8] we re-sampled in thirteen out of the sixteen locations included in that paper.

Table 1.List of populations, genetic indices and different localities/samplings included in
the study

Figure 1.Geographic distribution of sampling sites for the O. montagui populations included in the study. Details on each sampling site are given in Table 1. The names of the six sampled Mediterranean sub-basins are indicated along with the
pattern of the main surface sea currents.

The aim of this study is thus to determine the genetic structure of the beachflea
O. montagui in the central portion of the Mediterranean Sea. The species is tightly linked to
stranded heaps of the seagrass P. oceanica and hence passive dispersal via rafting is a plausible hypothesis. Here we specifically
ask whether such passive dispersal in O. montagui is a strong enough factor to promote genetic homogeneity over long distances. P. oceanica wrack may provide food and shelter; given the wide distribution of this seagrass
in the Mediterranean Sea we suspect events of transports to be quite frequent. On
the other hand, being strictly associated to a potential vector of dispersal does
not necessarily imply that gene flow is realized. We indeed found a remarkably deep
geographic partition in the population genetic structuring of the talitrid M. remyi, a species strictly associated to stranded rotting logs that should theoretically
aid passive dispersal [18]. Here we used multiple genetic data sets (mtDNA, microsatellites and allozymes),
which allowed testing the following two alternative scenarios (marine and terrestrial,
respectively). Under the marine scenario if passive dispersal at sea were predominant,
we would expect patterns of genetic structuring to be relatively shallow within basins
and largely congruent with those of marine Mediterranean taxa with planktonic stages
in the life cycle. Multiple authors recurrently placed a major breakpoint to gene
flow in marine species at the Strait of Sicily (STS): this separates the Western from
the Eastern Mediterranean basin (the so-called STS boundary; [24-27]). We also expect sea surface currents - the main carriers of floating wracks - to
shape patterns of relationships among populations. Alternatively, if the bound to
terrestrial environment prevailed over dispersal at sea, we would expect more pronounced
genetic divergence also at the within-basin level with evidence of genetic regionalism.
Support for both scenarios was evaluated with Approximate Bayesian Computation [28]. Finally, the study represents an opportunity to compare three marker systems characterized
by different evolutionary properties regarding inheritance (maternal for mtDNA vs.
bi-parental for microsatellites and allozymes) and pace of molecular evolution (from
relatively slow for the protein-coding allozymic loci to very fast for the non-coding
microsatellites).

Materials & methods

The study complies with the current laws of Italy and Germany regarding proper treatment
of animals for research and was approved by the Ethical Committee of University of
Rome "Sapienza", Italy.

Sampling

O. montagui specimens were collected by hand or using an aspirator on P. oceanica decaying leaves. A total of 26 populations were sampled for this study on rocky and
sandy beaches along shores of the following Mediterranean basins: Tyrrhenian Sea (10
populations), Western Mediterranean basin (12 populations), Eastern Mediterranean
basin (one population), Ionian Sea (one population), Adriatic Sea (two populations).
All 26 populations were analyzed at the mtDNA level and 21 (out of 26) were genotyped
for microsatellites.

Additionally, we re-analyzed the dataset based on 21 allozymic loci of De Matthaeis
et al. [8]. That study included 16 populations and 363 individuals. One population was from
the Aegean Sea (an area that we could not cover with our new sampling) while 13 have
been re-sampled for the present study. The different data sets combined hence include
29 populations: information about populations, sampling localities and markers used
is listed in Table 1; all populations analyzed with mtDNA and microsatellites were collected anew for
the study while De Matthaeis et al. [8] was the only source of the allozymic data. Figure 1 shows the geographic location of sampling sites.

Genetic analyses

Including the 153 samples analyzed in Pavesi et al. [23], we extracted total genomic DNA from 544 individuals using the “Qiagen DNeasy Blood
& Tissue Kit” (Hilden, Germany) according to the manufacturer’s protocol.

Mitochondrial DNA

Folmer et al.’s [29] universal primers were used to amplify via Polymerase Chain Reaction (PCR) a 571
base pair long fragment of the mtDNA coding for the COI gene in 1 to 15 individuals
per population for a total of 295 individuals (details are given in Table 1). COI double-stranded amplifications, purification of PCR products and sequencing
reactions were carried out as described in [18,23,30]. Editing and alignment of sequences were performed using Sequencher 4.1 (Gene Codes,
Ann Arbor, MI, USA) and checked by eye. We used Paup* v.4.0b10 [31] to calculate the number of parsimony-informative characters. We used Arlequin v.3.11
[32] to calculate genetic diversity (H), nucleotide diversity (πn) and the mean number of pairwise differences between all pairs of haplotypes (π). The same software was used to assess the population genetic structure by hierarchical
analysis of molecular variance (AMOVA; [33]). Populations were combined according to their geographic origin and to the results
of phylogeographic analyses (see Results). AMOVA was initially run on populations
combined in three groups. These pooled all Tyrrhenian samples (populations 1–2 and
15–24 in Table 1), North Western Sardinian + Ionian and Adriatic samples (populations 3–7 and 26–28)
and South Western Sardinian + Gozo Is. populations (8–14 and 25). An additional AMOVA
was run on five groups of populations. These included Corsican + South Eastern Sardinian
+ Tyrrhenian peninsular sites (populations 1–2; 15 and 23–24), North Western Sardinian
localities + Gozo Is. (populations 3–7 and 25), South Western Sardinian localities
(populations 8–14), North Eastern Sardinian localities (populations 16–21) and Adriatic
+ Ionian samples (populations 26–28). We also calculated pairwise FST values on all pairs of populations as a measure of their differentiation. Statistical
significance of these values was assessed by 1,000 permutations after sequential Bonferroni
correction for multiple tests [34]. We also performed a Mantel Test [35] with Arlequin, to test for the presence of Isolation-by-distance (IBD) in the dataset;
variables were pairwise FST values and corresponding geographic coordinates distances Degrees-Minutes-Seconds
(DMS). Geographic distances were measured as the shortest marine distances between
each pair of populations, according to the circulation of surface currents [36,37]. We used the spatial analysis of molecular variance (SAMOVA v.1.0) to detect genetic
barriers among groups of populations that were geographically homogeneous and maximally
differentiated from each other [38]. Arlequin was also used to assess the fit of demographic [39] and geographic [40] expansion models to O. montagui genetic variation (populations grouped according to the SAMOVA results). Mismatch
distributions were used to fit the models of sudden demographic and geographic expansion
and to estimate the τ parameter. In both models τ is a mutation-scaled measure of
time since expansion (demographic or geographic) and can be used to estimate time
in generations (T) since expansion, using T = τ/2 μ, where μ is the crustacean COI
mutation rate per site per year of Knowlton et al. [41] already adopted in Pavesi et al. [23]. Goodness of fit was assessed by the sum of square deviations (SSD) between the observed
and the expected mismatch, and its significance was determined by a parametric bootstrap
with 10,000 replicates. Demographic equilibrium was also tested on each SAMOVA group
by Tajima’s D[42] and Fu’s Fs[43] statistics in Arlequin. These parameters have been shown to be sensitive to signals
of population expansion; their statistical significance was tested by simulating random
samples (10,000 replicates). P-values were obtained as the proportion of simulated values smaller than or equal
to the observed values (α = 0.05).

In order to estimate genealogical relationships among all haplotypes, we used Tcs
v.1.13 [44] to carry out a statistical parsimony analysis with a 95% connection limit [45]. To better visualize the haplotype network layout, we used HapStar [46] combined with Inkscape, a free open-source Svg graphics editor [47]. HapStar directly combines the network connection output data generated from Arlequin
with a force-directed algorithm to automatically lay out the network for easy visualization.

Microsatellites

We genotyped all 544 individuals (7 to 30 per population) at eight microsatellite
loci (msats, Table 1). PCR amplifications were performed as described in [22,23]. Fragment sizes were determined on an ABI 3100 automatic sequencer using the GeneMapper
v.3.5 software and an internal size standard (GS500LIZ from Applied Biosystems). We
used Microchecker v.2.2.3 [48] to identify genotyping errors (null alleles) and the scoring of stutter peaks in
all microsatellites. Null alleles can cause deviations from Hardy–Weinberg proportions.
In order to test for the impact of null alleles on our data set, we calculated global
FST values with FreeNA software. These values were computed, as described in Chapuis
and Estoup [49], with 10,000 bootstrap iterations, alternatively using and not using the Excluding
Null Alleles (ENA) method. We used Arlequin v.3.11 [32] to calculate expected (HE) and observed (HO) heterozygosity, deviations from Hardy–Weinberg equilibrium and tests of linkage
disequilibrium among loci. We corrected P-values using the sequential Bonferroni correction for multiple comparisons [34]. We used Fstat v.2.9.3 [50,51] to calculate allelic richness (AR) for each population, and GenAlex v.6.3 [52] to determine the number of alleles per locus, the number of private alleles (PA) and their relative frequencies in each population. With Arlequin v.3.11 we calculated
pairwise FST values and performed two AMOVA analyses (on three and five groups, respectively)
and Mantel test as described for the mtDNA dataset in the previous paragraph. For
the AMOVA we grouped populations based on the results of Structure v.2.3.2 [53] (three and five groups; see below).

We run Structure to evaluate genetic subdivision and to assign individuals to the
most probable number of clusters (K). The software uses a Monte Carlo Markov Chain
(MCMC) Bayesian clustering model that assumes that different loci are at Hardy-Weinberg
and linkage equilibrium with one another within populations. We performed 10 independent
runs for each value of K (K = 1–21) assuming an admixture model and using 100,000
replicates of the MCMC after a burn-in of 10,000 iterations. The number of clusters
was calculated from the ΔK estimates using Evanno et al.’s method [54]. Assignments of individuals to the inferred clusters were estimated according to
the highest Q-values (probability of membership).

We compared the results obtained from Structure with those from Tess v.2.3 [55]. The software is based on a hierarchical mixture model whose neighborhood system
is obtained from the Voronoi tessellation [55]. We performed 10 runs (K = 2–21) with an admixture model with a total of 50,000 sweeps
with a burn-in of 10,000. The number of clusters was calculated plotting the Deviance
Information Criterion (DIC) against values of K and choosing the value of K that corresponded
to a plateau of the curve [56].

Allozymes

We re-analyzed the dataset of De Matthaeis et al. [8] consisting of 21 allozymic loci, 16 populations and 363 individual (2 to 105 per
population; see Table 1). We calculated pairwise FST values, AMOVA and IBD using Arlequin. We performed the AMOVA on three groups combining
populations according to their geographical position and patterns of relationships
as detailed in De Matthaeis et al. [8] (Tyrrhenian + Aegean Sea; North Western Sardinia; South Western Sardinia). We analyzed
the genetic structure of populations by running Structure (K = 1–16) and Tess (K = 2–16)
with the settings described in the microsatellite section.

Approximate Bayesian computation

We used Approximate Bayesian Computation in DIYABC 1.0.4.39 [57] to evaluate the two competing scenarios (marine vs. terrestrial) detailed in the
Introduction. To reduce computation time, we selected nine populations (SCP, SNC,
STA, SAR, CCL, FRA, GOZ, PMP, PSF) as representative of each sampled area. We chose
populations with complete data sets for both mtDNA and microsatellites and with the
lowest divergence from the other populations from the same area. The marine scenario
was designed following the subdivision of the Mediterranean Sea in biogeographic sectors
proposed by Bianchi & Morri [58]; these would group the West Mediterranean basin and the Tyrrhenian Sea together (hence
SCP, SNC, STA, SAR, CCL, FRA were forced to merge into a common ancestor), the Maltese
Archipelago (GOZ) would be a transition area between the former sector and the Ionian
one (PMP) while the South Adriatic Sea (PSF) would form a separated biogeographic
area. The terrestrial scenario was designed on the basis of the geographic origin
of populations. We considered the following areas; North-Eastern (STA), North-Western
(SAR), South-Western (SNC), South-Eastern Sardinia (SCP), Tyrrhenian Sea (CCL, FRA;
one insular and one peninsular population), Maltese Archipelago (GOZ), Ionian (PMP)
and Adriatic (PSF) Sea. Three data sets were considered in the analysis (mtDNA, microsatellites,
both markers combined); for each scenario we simulated 2 × 106 data sets. For mtDNA we assumed the default settings in DIYABC for mitochondrial
markers; for microsatellites the loci were assumed to follow a generalized stepwise
mutation model. We calculated the posterior probabilities of the two competing scenarios
for each data set by i) evaluating the relative proportions of each scenario found in the selected closest
data sets (Direct estimate; 500 simulated data sets) and through ii) a polychotomous logistic regression on the 1% of the simulated data sets closest
to the observed data sets (default settings in DIYABC). Before estimating the posterior
probabilities of scenarios, we checked whether the combination of scenarios and prior
distributions of their parameters have produced data sets that are similar enough
to the observed ones by performing a Principal Component Analysis (PCA) of the first
100,000 simulated data sets to evaluate the position of the observed data compared
to the simulated.

Values of H, πn and π are listed in Table 1. H ranged between 0.362 (SNG) and 0.971 (TTA); π ranged between 0.381 (SNG) and 9.571 (PSF) whereas πn ranged between 0.001 (SNA, SNG, SAR) and 0.017 (PSF). The AMOVA detected significant
population structure when we tested for three and five groups (Table 2): in both cases the major source of variation was due to differences within populations.
A small percentage of the detected variation was apportioned among groups; this increased
when five groups were considered. Pairwise FST values are shown in Table 3: values ranged between 0.049 (SAF-TTA) and 0.596 (SSS-SNC). Most comparisons were
significant after the Bonferroni correction with some exceptions. In particular, CSR
(from Corsica Is.) diverged significantly from only two Sardinian populations (SSS
and SNC) whereas the very close CCL (Corsica Is.) was significantly differentiated
from it (and many others) but not from the far away PSC (Adriatic Sea). Significant
was also the comparison between PMP and PSF (Ionian and Adriatic Sea, respectively).
A Mantel test performed on all populations revealed no statistically significant relationship
between genetic and geographic distances (R = 0.000005; P = 0.426).

SAMOVA suggested the existence of three groups (FCT = 0.45; P < 0.05), finding two barriers isolating the Adriatic from the Ionian populations
and then these two basins from the rest. The analysis of mismatch distributions supported
demographic expansion in the Adriatic group (PSSD = 0.283) only; geographic expansion was supported in all groups (PSSD ≥ 0.184). Tajima’s D and Fu’s Fs were not significant for the Ionic (D = −0.792, PD = 0.217; FS = 3.854, PFS = 0.95) and Adriatic (D = −1.468, PD = 0.06; FS = 0.932, PFS = 0.706) groups while they were negative and highly significant for the remaining
populations lumped in a single group (D = −2.157, PD = 0.000; FS = −25.189, PFS = 0.000). In light of the relatively poor performance of mismatch distribution statistics
in detecting population demographic expansion compared to D and FS statistics [59], we tend to favor the latter two, which support demographic expansion in the Western
Mediterranean and Tyrrhenian populations. In the Adriatic group the timing of demographic
expansion (τ = 34.63; T = 1.731 million years) preceded the spatial expansion (τ = 7.52;
T = 376,000 years). In the Ionic basin the geographic expansion dates back to 500,000 years
ago (τ = 10). In the Western Mediterranean and Tyrrhenian Sea the two expansions were
more recent and coeval (τ = 1.01; T = 50,500 years), justifying the star-like shape
of the haplotype network (see below). The relatively small τ values suggest that expansions
always started from small effective population sizes.

The statistical parsimony network is shown in Figure 2. Altogether, 53 out of 82 haplotypes were singletons. The analysis yielded one network
linking all haplotypes. Three main haplogroups could be recognized: two haplogroups
were star-like shaped and most of the sampled Mediterranean basins were represented
in there with the exception of the Ionian and Adriatic ones. These were instead loosely
gathered in a third sub-group, together with some Tyrrhenian and Western Mediterranean
haplotypes and many intervening missing haplotypes. In four circumstances rare haplotypes
from different basins were very divergent from all others, being separated by a high
number of substitutions (bracketed in Figure 2).

Microsatellites

All microsatellite loci were polymorphic in the 21 studied populations with the number
of alleles per locus ranging from 7 to 16. Number of individuals analyzed, HE and HO, AR and PA are provided in Table 1. No scoring errors and no evidence of linkage disequilibrium were detected between
any pairs of loci after Bonferroni correction. For some loci Microchecker indicated
possible presence of null alleles; global FST values calculated alternatively including and excluding null alleles returned very
similar values (Table 4) suggesting very little effect of null alleles on FST estimates.

AMOVA indicated a significant genetic structure and most of the genetic variability,
when populations were grouped either in three (83.01%) or five groups (83.15%), was
apportioned within populations (Table 2). Population differentiation (pairwise FST, Table 3) was significant for all pairwise comparisons except two comparisons (TTA-SCP and
CCL-SCP). The Mantel test did not detect a pattern of IBD (R = 0.0001; P = 0.465).

Results from Structure and Tess indicated respectively K = 3 (LnP(K) = −13309.7) and
K = 5 (DIC = 18998.5) as the best fitting our data (Figure 3). We also analyzed results for K = 5 in Structure (LnP(K) = −13310.1) and K = 3 (DIC = 20242.1)
in Tess for a better visualization of data and further comparison. For K = 3 in Structure,
the three clusters corresponded to 1) Western Mediterranean basin (Northern Sardinian
populations: SAR-SAP-SAB-SAF-SSS) plus Ionian and Adriatic populations (PMP and PSF);
2) Tyrrhenian basin (STA-SMM-SMS-SMT-TTA-FRA-CCL); 3) Western Mediterranean basin
(Southern Sardinian populations: SPC-SNS-SNG-SNA-SNC) plus a single population from
the Tyrrhenian basin (SCP). The Eastern Mediterranean population (GOZ; Maltese Archipelago)
was shared between cluster 1 and 3. Tess with K = 3 also retrieved the Tyrrhenian
cluster, but the Ionian and Adriatic populations clustered on their own while the
third cluster grouped all populations from the Western Mediterranean basin plus GOZ
(Eastern Mediterranean basin).

Figure 3.Patterns of genetic divergence at the eight microsatellite loci as revealed by Structure
(above) and Tess (below) for K = 3 and K = 5. Each vertical line represents one individual. The length of each line reflects the
percentage probability of each individual’s membership to each cluster. For K = 3
clusters are shown in green, blue and red. For K = 5 clusters are indicated in light
green, red, blue, orange and dark green. Codes for populations are as in Table 1.

Intriguingly, individuals carrying the four deeply divergent COI haplotypes (Figure 2) do not possess a similarly deviant microsatellite genotype. Membership of the Tyrrhenian/Ionian,
S-W Sardinian, and Adriatic individuals to the respective clusters ranges between
0.84 and 0.98 (Structure and Tess). Tess assigned the N-W Sardinian individual to
the Northern Sardinian group with a membership of 0.75 while Structure suggested for
this individual a membership shared between cluster 1 and 2 (0.36, 0.43).

Allozymes

We re-analyzed a dataset of 16 populations (13 of them in common with our data set;
Table 1) genotyped at 21 allozyme loci by De Matthaeis et al. [8]. Pairwise FST values ranged between 0.03 (FRA-SNA) and 0.49 (STA-SAP) (Table 5). STA was significantly differentiated from all populations, whereas the Greek population
(GRE) was significantly differentiated only from STA (Tyrrhenian basin). Northern
and Southern Sardinian populations (Western Mediterranean basin) were consistently
differentiated from one another with the exceptions of SPC (Southern Sardinia). Interestingly,
significant differentiation emerged also at a scale as low as that of the Asinara
Is., where SAR was found to be significantly divergent from both SAP and SAB. The
AMOVA analysis revealed that 87.3% (P = 0.0001) of the genetic variance was apportioned within populations. A Mantel test
revealed a lack of significant IBD when the Greek population (GRE) was excluded from
the analysis (R = 0.0001; P = 0.006), whereas IBD was detected when GRE was included in it (R = 0.0001; P = 0.001).

Table 5.Pairwise FST values for all comparisons among populations with allozymes

Both Structure and Tess analyses performed on this dataset revealed K = 3 ((LnP(K) = −1767.52
and DIC = 2887.23 respectively) as the best describing the species population structuring
(Figure 4a, b). Nevertheless, membership of individuals to the inferred clusters was always
very low in Structure and no real structuring was evident (Figure 4a). Results from Tess (Figure 4b) gave instead a clear subdivision in three clusters. A first one included populations
located in the northern (respect to our samplings) part of the Tyrrhenian (populations
FRA, SMP, SMT, SMM, SCA, STA) and Western Mediterranean Sea (populations SAP, SAB,
SAR); a second cluster grouped populations from the southern part of the Tyrrhenian
(SCP) and the Western Mediterranean Seas (SNC, SNL, SNA, SPC). A population from the
Maddalena Is. (SMS, Tyrrhenian Sea) was included in the second cluster in spite of
the northern position of the island respect to Sardinia (see Figure 1). Memberships of these populations to this cluster were low. The third cluster included
GRE only (Aegean Sea). Finally, to further compare microsatellite and allozyme data
we ran the latter in Structure and Tess enforcing K = 5 (results not shown graphically).
In Structure no pattern emerged (membership Q < 0.4 in all cases except GRE); Tess
couldn’t detect any additional cluster, as individuals were partitioned again among
the three above-mentioned clusters.

Figure 4.Patterns of genetic divergence at twenty-one allozymic loci. a,b Allozymes results for Structure (a) and Tess (b) (K = 3). Histograms show the average membership to one of the three clusters of
each population (see Table 1 for details). Clusters are indicated by light green, blue and dark green and correspond
to Tyrrhenian Sea and Northern Sardinia (Western Mediterranean Sea) (light green),
Southern Sardinia (Western Mediterranean Sea) (blue) and Aegean Sea (dark green) (see
text for details).

Approximate Bayesian Computations

The PCA analyses conducted to check whether the combination of scenarios and prior
distributions of their parameters have produced data sets similar enough to the observed
ones revealed in all circumstances that observed data were surrounded by many simulated
data (not shown). Figure 5a, b, c shows the probabilities calculated with the two approaches (direct estimate
and logistic regression) for both scenarios and for each marker system separately
and together. Results are remarkably congruent across data sets and approaches; ABC
consistently identified the marine scenario as largely supported over the terrestrial
one.

Figure 5.Results of the Approximate Computation Analyses. a,b,c ABC results for mtDNA (a), microsatellites (b) and both marker systems combined (c). On the left are shown the posterior probabilities based on direct estimates, on
the right posterior probabilities are obtained through a logistic regression. In all
circumstances the marine scenario is better supported than the terrestrial one.

Discussion

In the last two decades we have been accumulating genetic data on Mediterranean talitrids
with the aim to comparatively describe the genetic structure of different species
at the scale of the whole Mediterranean Basin (see [60] and references therein for a review on the issue). So far, the present contribution
is the largest one both in terms of individuals screened and number of markers used
focused on a single talitrid species.

Present genetic data suggest that the beachflea O. montagui has a complex phylogeographic structure. By using ABC, we found that a scenario with
marine forces assumed as predominant in shaping the species genetic structure is clearly
better supported than a scenario of regionalism based on the geographic origin of
populations (terrestrial scenario). The marine scenario is favored regardless the
marker system analyzed (mtDNA and microsatellites; either separately or combined).
It is worth noting that, in spite of this robust evidence, patterns of relationships
yielded by mtDNA and microsatellites (and also allozymes but to a lesser extent) are
not mirroring one another but major differences emerged that can be explained considering
the intrinsically different evolutionary properties of these markers.

Populations east and west of the theoretical STS boundary are significantly differentiated
at the mtDNA level (SAMOVA results) while microsatellites are ambivalent on the issue
(compare the most likely groupings identified by Structure and Tess; three and five
respectively). The population sampled on Gozo Is. (Maltese Archipelago; Eastern Mediterranean
basin) shares both mtDNA haplotypes and microsatellite alleles with Western Mediterranean,
Tyrrhenian, Ionian and Adriatic populations. The lack of isolation by distance in
the mtDNA and both nuclear datasets excludes geographic distance as the main determinant
of genetic divergence. On the other hand, the vast majority of pairwise FST values and the AMOVA results are statistically significant regardless of the marker
considered. This evidence points to an overall scenario of complex genetic fragmentation
that needs to be carefully evaluated to correctly resolve population structure in
the species.

As pointed out by Thiel and Haye [61], population connectivity of rafters depends on a variety of factors (distances between
habitats, sea current patterns, dispersal capability and colonization potential of
organisms) and unraveling the relative importance of each one is a hard task that
would require coupling genetic and field data. O. montagui is a small enough organism to hypothesize that it could cling, find shelter and feed
on P. oceanica wrack while at sea for as long as the wrack would persist and/or strand in a new
location. The reproductive biology of amphipods may enhance dispersal and recruitment;
by incubating the embryos and developing juveniles in the brood pouch females provide
extended parental care [62]. Migrating gravid females would favor persistence of rafters over generations. The
four deeply divergent COI haplotypes that we found in our dataset (Figure 2) could indeed represent immigrants from geographical areas that we did not cover
with our sampling. The mtDNA legacy of these alleged immigrants has persisted in five
locations due to the mtDNA maternal inheritance while their ancestry has vanished
at the microsatellite level due to the bi-parental inheritance of these markers. The
sporadic occurrence of haplotypes profoundly divergent from an otherwise locally homogeneous
set of lineages has been already reported in the talitrid T. saltator[17].

In spite of this convincing evidence, we suggest caution in interpreting the previous
and following findings in light of our incomplete sampling coverage in the Eastern
Mediterranean Sea and the small sample size for one Adriatic population.

Divergence among and within basins: historical vs. contemporary processes

The SAMOVA analysis conducted on mtDNA data detected barriers to gene flow between
Adriatic and Ionian populations and between these two and the remaining basins. Patterns
of mtDNA relationships within each of the above groups are not similar. Adriatic and
Ionian haplotypes are scattered in the network of Figure 2 and many intervening missing haplotypes need to be postulated to connect them. Conversely,
the vast majority of haplotypes from the Tyrrhenian, Western and Eastern Mediterranean
basins clusters around the two most frequent mtDNA variants and differs from them
by one single mutational step. This indicates that haplotypes coalesce earlier into
their last common ancestor in the Tyrrhenian, Western and Eastern Mediterranean basins
than they do in the Adriatic and Ionian ones. A variety of reasons either intrinsic
to the species ecology (different rates of population crashes) or basin-specific (differences
in rates of formation of favorable conditions, i.e. heaps of stranded seagrasses,
intermittent rafting routes, marine genetic barriers) could be responsible. Field
studies aimed to gain a detailed understanding of the factors regulating the species’
niche are needed to properly understand the causes of these alternate mtDNA patterns.

Western Mediterranean/Tyrrhenian, Ionian and Adriatic mtDNA clusters (i.e. the three
groups recognized by SAMOVA but see also ABC results) perfectly overlap with the corresponding
biogeographic provinces identified within the Mediterranean Sea [58]; this evidence suggests historical factors to be responsible for the mtDNA pattern.
Barriers to mtDNA gene flow can be traced back to the Pleistocene [23], when the sea level repeatedly dropped at each ice age peak and the different Mediterranean
basins shrunk and became largely isolated from one another ([58] and references therein). The demographic and geographic expansions detected in our
mtDNA data sets are indeed recurrently placed after major drops in sea level [63].

Microsatellites (and to a lesser extent allozymes) cluster individuals preferentially
concordant to their geographic origin, while many mtDNA haplotypes from geographically
far away locations are recurrently gathered together in the haplotype network and
differ by one or few substitutions. The geographically close Ionian and Adriatic populations
have quite similar microsatellite profiles (Figure 3) whereas divergence emerges between Northern and Southern Western Sardinian populations.
Divergence was detected also among the Tyrrhenian, Western and Eastern Mediterranean
basins. A closer look at the haplotype network of Figure 2 reveals that Northern and Southern Western Mediterranean populations have inverted
frequencies of the two most common haplotypes, giving partial support (yet not detected
by SAMOVA) to the microsatellite evidence. Altogether, patterns of relationships at
the mtDNA level have relatively less geographical structure than those yielded by
microsatellites. How could we then reconcile these discrepancies in light of the inherently
different evolutionary properties of mtDNA and microsatellites [64]? The finer degree of microsatellite structuring would argue in favor of a female-biased
gene flow. In our recent study on O. montagui based on a subset of six populations, we invoked male-biased dispersal to explain
why the mtDNA fragmentation was not apparent at the microsatellite level [23]. In the current study, which relies on a better coverage of the species geographic
distribution, the pattern is reversed. In light of this, we tend here to favor the
female-biased gene flow hypothesis. This would also be in agreement with field studies
showing that populations of Mediterranean talitrids, including O. montagui, generally have a strongly female-biased sex ratio ([65] and references therein).

O. montagui is able to disperse actively only over short distances; individuals could nonetheless
be dragged away by waves within floating banks of P. oceanica and – theoretically –carried over ample ranges. The pattern of genetic structuring
yielded by microsatellites is, in our opinion, reflecting these present-day processes
and is influenced by the surface circulation of water masses, which are the principal
carriers of floating material at sea. Genetic connectivity mirrors quite remarkably
patterns of Mediterranean surface currents. Water masses flow in a south–north direction
along the Tyrrhenian coast of the Italian peninsula and bifurcate along the Tuscan
coasts into two branches (see Figure 1 modified from [37]). The first branch flows northward along the eastern coast of Corsica while a second
one moves southward following the Sardinian eastern coasts. This pattern would explain
why microsatellites consistently cluster Tyrrhenian populations together. Along the
Sardinian western coasts water masses flow in a northward direction and a major gyre
exists in the area, which is likely the cause of the microsatellite divergence found
between South- and North-Western Sardinian populations (Western Mediterranean basin).
The currents flowing in and out of the Adriatic Sea would promote gene flow between
the Adriatic and Ionian populations. Finally, the highly admixed genetic profile of
the sole Eastern Mediterranean population included in the study (GOZ) might be due
to the position of the Maltese Archipelago at the crossroads between the main current
flowing from the Atlantic Ocean to the Aegean Sea and that descending along the Ionian
and South-eastern Sicilian coasts. The former would carry immigrants of Western Mediterranean
and Tyrrhenian origin while the latter would be the source of Ionian and Adriatic
alleles.

Comparisons to other Mediterranean species

We could not detect a decrease in genetic variability from the Western to the Eastern
Mediterranean Sea, as previously shown by Viñas et al. [66] in swordfish. Multi-taxa evidence supports a marked decrease in variability between
the Mediterranean and the Black Sea but not necessarily in a west–east direction within
the Mediterranean Sea [67].

When we compare our data within those available on a similar geographic scale for
other Mediterranean species an interesting result emerges. Several marine species
are genetically more homogeneous in the Western Mediterranean basin than they are
in the Eastern Mediterranean, Ionic or Adriatic Sea. This was recurrently observed
in fishes (i.e. gobids; with Adriatic, Aegean and Tunisian groups, [25]; anchovies, with a clear distinction between Adriatic, Aegean and the rest of Mediterranean
Sea, [68]). O. montagui partially follows this scheme; the uniqueness of the Adriatic and Ionic populations
is evident at the mtDNA level only. Likewise, Western and Tyrrhenian basins are shallowly
fragmented at the mtDNA level while microsatellites revealed within-groups patchiness.
Interestingly enough, allozymes detected a fixed alternative allele between the Aegean
population (GRE) and all the others analyzed in [8]. Unfortunately, since we did not have access to any samples from the Aegean Sea for
mtDNA or microsatellites, we cannot be conclusive on the issue.

Previous studies on Mediterranean talitrid amphipods coarsely found that species with
a profound population genetic structure are generally confined to sandy beaches where
they burrow into sand (i.e. T. saltator, D. deshayesii) or dig into rotting logs (M. remyi) [8,15,18,60]. Various species of Orchestia, including O. montagui, associated with heaps of decaying wrack were thought to be genetically homogeneous,
even across vast geographical ranges [8-15]. The spatial scale at which genetic divergence would become apparent was found to
vary considerably in talitrids, from few kilometers with no sharing of mtDNA haplotypes
in M. remyi[18] to hundreds of kilometers in O. montagui ([8]; allozyme data). In T. saltator, allozymes clearly distinguished reciprocally isolated Tyrrhenian, Adriatic and Aegean
groups [8].

The intrinsic active mobility of organisms is only one determinant of the amount of
the realized dispersal [69]. Active dispersal is always limited in talitrids, independently from the single species’
ecology. Hence, differences in their degree of genetic structuring must have arisen
because of differences in the amount of passive dispersal the various species are
able to achieve. M. remyi lives in strong association with rotting logs stranded on sandy beaches; T. saltator burrows deeply into the sand to escape dehydration. In both circumstances, the chances
for individuals to be dragged away by waves are far lower than those of O. montagui, which colonizes a temporally ephemeral habitat close to the water line. This would
explain why our focus species was thought to be genetically poorly structured [8-15]. On the other hand, we believe that the conclusions to which previous studies on
Mediterranean talitrids have arrived should be critically evaluated, especially for
those species identified as able to sustain high levels of gene flow at the scale
of the whole Mediterranean Sea. We have here shown that O. montagui can no longer be considered a species with a shallow genetic structure, as De Matthaeis
et al. [8-15] concluded on allozyme basis. We have demonstrated that long distance dispersal through
surface currents might happen in this species but these episodes are not frequent
enough to ensure high levels of gene flow that would counteract the onset of appreciable
local differences. A possible limiting factor could be identified in the relatively
short persistence in time of the vector of passive dispersal (the P. oceanica wrack) at sea surface due to its negative buoyancy [61].

By comparing patterns yielded by microsatellites and allozymes for the same set of
populations, we obtained an identical number of clusters (3) but membership of individuals
was consistently higher in the microsatellite analyses, pointing – as expected – to
a much better resolving power of these markers. The molecular footprint of high gene
flow under low-resolution molecular assays that O. montagui was believed to bear has in fact converted to an appearance of intermediate gene
flow under more stringent assays.

Conclusions

We were able to unveil the genetic structuring of O. montagui within the Mediterranean Sea at a fine geographical scale. By combining molecular
markers with different evolutionary properties we distinguished historical processes
from present-day forces. Pleistocene climatic changes caused the isolation and divergence
of mtDNA lineages. Current dynamics – namely surface marine currents jointly with
rafting via seagrass wrack – produced new and different patterns of structuring, which
were revealed by microsatellites. ABC gave strong support to a scenario of genetic
fragmentation molded predominantly by marine forces.

The lack of a statistically significant pattern of IBD at both classes of markers
under a scenario of genetic fragmentation (see FST and AMOVA results) suggests that the species has not yet approached equilibrium between
gene flow and drift [70], probably because expansions (both demographic and geographic) started from small
effective population sizes and are relatively young, as revealed by the mismatch analyses.
On the other hand, this conclusion should be considered with caution in light of the
fact that we could cover only a small fraction of the eastern part of the species’
ranges.

Our results also pinpoint the difficulties in describing the population genetic structure
in an ecotonal species, because of its simultaneous exposure to forces typical of
two ecosystems (land and sea). The semi-terrestrial life style acquired by O. montagui prevents the species from realizing long-distance dispersal at sea at a rate comparable
to that of truly marine organisms, especially those with a planktonic stage in their
life cycle. The sea, on the other hand, functions as a means for passive dispersal
by carrying individuals of O. montagui with floating materials across sites. The species’ genetic structure can indeed be
largely explained in light of the pattern of circulation of surface sea currents.
This passive process, being intrinsically unpredictable and variable, cannot counteract
the onset of local genetic differences. Passive dispersal via floating material is
extremely difficult to observe in nature and no information exists either on the rate
of survival of these amphipods during the process or on their chances to settle on
new sites following dispersal. Future avenues of research should foresee the development
of ecological experiments capable of tracking single biological particles (talitrids
in this case) in floating wracks to ultimately link genetic and oceanographic data
in a unified framework [71].

Competing interests

The authors declare that they have no competing interests.

Authors’ contribution

LP collected the samples, designed the project, carried out the experiments and the
statistical analyses, interpreted the data, drafted and critically revised the manuscript.
VK collected the samples, designed the project, carried out part of the statistical
analyses, interpreted the data and critically revised the manuscript. RT provided
funding for all experiments and critically revised the manuscript. EDM provided the
allozymic data and critically revised the manuscript. All authors read and approved
the final manuscript.