Skip to main content

Parapatric genetic divergence among deep evolutionary lineages in the Mediterranean green crab, Carcinus aestuarii (Brachyura, Portunoidea, Carcinidae), accounts for a sharp phylogeographic break in the Eastern Mediterranean



Recently, population genetic studies of Mediterranean marine species highlighted patterns of genetic divergence and phylogeographic breaks, due to the interplay between impacts of Pleistocene climate shifts and contemporary hydrographical barriers. These factors markedly shaped the distribution of marine organisms and their genetic makeup. The present study is part of an ongoing effort to understand the phylogeography and evolutionary history of the highly dispersive Mediterranean green crab, Carcinus aestuarii (Nardo, 1847), across the Mediterranean Sea. Recently, marked divergence between two highly separated haplogroups (genetic types I and II) of C. aestuarii was discerned across the Siculo-Tunisian Strait, suggesting an Early Pleistocene vicariant event. In order to better identify phylogeographic patterns in this species, a total of 263 individuals from 22 Mediterranean locations were analysed by comparing a 587 basepair region of the mitochondrial gene Cox1 (cytochrome oxidase subunit 1). The examined dataset is composed of both newly generated sequences (76) and previously investigated ones (187).


Our results unveiled the occurrence of a highly divergent haplogroup (genetic type III) in the most north-eastern part of the Mediterranean Sea. Divergence between the most distinct type III and the common ancestor of both types I and II corresponds to the Early Pleistocene and coincides with the historical episode of separation between types I and II. Our results also revealed strong genetic divergence among adjacent regions (separating the Aegean and Marmara seas from the remaining distribution zone) and confirmed a sharp phylogeographic break across the Eastern Mediterranean. The recorded parapatric genetic divergence, with the potential existence of a contact zone between both groups in the Ionian Sea and notable differences in the demographic history, suggest the likely impact of paleoclimatic events, as well as past and contemporary oceanographic processes, in shaping genetic variability of this species.


Our findings not only provide further evidence for the complex evolutionary history of the green crab in the Mediterranean Sea, but also stress the importance of investigating peripheral areas in the species’ distribution zone in order to fully understand the distribution of genetic diversity and unravel hidden genetic units and local patterns of endemism.


In the marine environment, genetic variability and population genetic structure of a species are shaped by both contemporary and historical marine barriers, such as narrow and shallow water passages between landmasses [1], salinity gradients [2], different types of currents [3], as well as palaeoecological history [4, 5]. In particular, historical processes have left marked imprints on the genetic structure of extant marine populations. For example, the extent of genetic divergence between genomes varies among marine species and can be used to estimate the time of their separation. Marine populations inhabiting specific regions (i.e., those confined to historical refugial zones) might have harboured persisting genetic lineages through time and therefore accumulated genetic differences, leading to their possible speciation through several climatic cycles (i.e., the Quaternary glaciations periods) [6].

Palaeoclimatic and palaeogeographical evolution of the Mediterranean, primed by sea-level low stands during glacial periods of the Pleistocene and the consequent shift in abiotic (i.e., temperature and salinity) and biotic factors (i.e., productivity), prompted population isolation and divergence, following local events of extinction and recolonization. This in turn enhanced the formation of intraspecific genetic units within various marine taxa, due to historical or contemporary hydrographic barriers to gene flow [4, 5]. The best known oceanographic discontinuities in the Mediterranean are: (1) the Almería-Oran Front as the main genetic breakpoint between the Atlantic Ocean and the Mediterranean Sea [4, 7], (2) the transition between the Western and the Eastern Mediterranean, due to unidirectional water flow at the Siculo-Tunisian Strait [5], and (3) the hydrographic isolation of the Aegean, Ionian and Adriatic seas [5]. In addition to these oceanographic barriers, the geological history of the Eastern Atlantic and Mediterranean Sea, including the breaking up of the Tethys Sea, the Messinian Salinity Crisis [8], and the Pleistocene glaciations [9], might have left marked footprints on the genetic structure of species and made the Mediterranean region a notably dynamic hotspot of diversity [10,11,12]. For example, intensified modifications to the coastline (sea level regressions) during repeated Pleistocene glaciations could have limited the biotic exchange across physical barriers, such as the Gibraltar Strait [13] and the Siculo-Tunisian Strait [14], and strongly influenced the formation of distinct phylogenetic lineages within Mediterranean marine species.

In this context, the phylogeography of Mediterranean species underwent a common set of processes resulting from fragmentation within glacial refugia, range expansions via postglacial colonisation, and secondary contact zones among historically divergent lineages [4, 12, 15,16,17,18].

The Mediterranean green crab, Carcinus aestuarii (Nardo, 1847) (Brachyura, Portunoidea, Carcinidae), represents a good model to test the impact of current and historical marine barriers on population subdivision and genetic structuring across various marine biogeographic boundaries in the Mediterranean Sea. This species is a very common inhabitant of estuaries and lagoons of the Mediterranean Sea and has also been reported from the southern Black Sea [19, 20]. It inhabits various types of environments, ranging from sheltered and often brackish habitats (including subtidal and intertidal mud, as well as sand flats at open coastal sites or lagoons and estuaries) to saltmarshes, and seagrass beds (authors’ personal observations). It is a voracious omnivore and aggressive competitor. Furthermore, it has a wide tolerance toward variations of salinity, temperature, and dissolved oxygen. Hence, it has a high ability to adapt to a wide variety of habitats [19, 21]. The Mediterranean green crab also exhibits high fecundity and a relatively long dispersal stages, with a larval planktonic phase of approximately six weeks [22, 23].

Owing to maritime commerce and ballast transport, several reports have pointed out to the accidental introduction of C. aestuarii specimens into several regions outside their native range, such as to the Canary Islands [24], Tokyo Bay, Japan [25] and South Africa [22, 26].

Among the above mentioned eco-biological characteristics of C. aestuarii, the particular ecological specialization of this species to estuarine and brackish-water habitats in the Mediterranean could be considered as important factor susceptible of generating significant phylogeographic patterns, owing to limited connectivity among estuarine populations. Restricted gene flow is highly expected in estuarine and brackish-water taxa owing to larvae retention within their natal estuary (e.g. [27,28,29]). In addition to this mechanism, geographically extended distances separating estuaries (e.g. [29]), the differential in physical characteristics of estuarine and their adjacent coastal environments [30], as well as the potential different physiological challenges between estuarine and coastal waters [31], could also impede inter-estuarine larval dispersal. Numerous investigations have shown that genetic connectivity among populations of estuarine taxa is more restricted than that occurring in taxa inhabiting the open coast [27, 32, 33].

Recent population genetic surveys on C. aestuarii, across native and invaded regions [22, 23, 34, 35], unveiled extensive genetic variability and marked population differentiation associated with the main oceanographic discontinuities that characterize these areas. In particular, Marino et al. [23] and Ragionieri and Schubart [35] found significant genetic differentiation among populations from western and eastern Mediterranean coastlines in Europe. This noticeable pattern of genetic structure has been later confirmed across the central African Mediterranean coast, when Deli et al. [36] revealed a sharp haplotypic discontinuity among eastern and western sites in Tunisia. Both retrieved genetic groups were found to be genetically and morphologically differentiated across the Siculo-Tunisian Strait [36, 37]. A more recent investigation by Deli et al. [12], detailing phylogeography and population genetic structure of the Mediterranean green crab across the Siculo-Tunisian Strait, showed concordant patterns of mitochondrial and nuclear divergence among western and eastern Mediterranean populations from northern Africa. The study also revealed a marked divergence between two highly separated haplogroups, suggesting an Early Pleistocene vicariant event in the genetic differentiation of this highly dispersive decapod species [12].

These insights, inferred so far from population genetic surveys of C. aestuarii throughout the western and central Mediterranean coasts [12, 23, 34,35,36], trigger the necessity of detailed phylogeographic and population genetic analyses of this species across the poorly investigated Eastern Mediterranean Basin. Sea-level Pleistocene oscillations potentially caused the isolation or partial isolation of the Black Sea, Aegean Sea and Eastern Mediterranean Basin [38]. These historical isolation processes have been maintained by the impact of the contemporary hydrographic isolation of the Adriatic, Ionian and Aegean seas [4, 5]. In addition to these historical and contemporary hydrographic isolation patterns, different selective forces related to the environmental conditions of the Ionian and Aegean seas could account for the phylogeographical break observed in the Eastern Mediterranean. Previous population genetic studies on many vertebrate and invertebrate species confirmed this tendency and identified a major genetic break related to the hydrographic isolation of the Aegean Sea [5, 16, 39,40,41,42,43,44,45].

In light of these considerations, it will be interesting to re-examine population genetic structure in a highly dispersive decapod species, like Carcinus aestuarii, across these potential barriers to gene flow. This may contribute to a better identification of the phylogeographic patterns in this species and unveil other potential evolutionary lineages within C. aestuarii across the poorly investigated easternmost part of its distribution range. It may also allow depicting the historical events that might have shaped the genetic structure of the Mediterranean green crab. In order to achieve these goals, new mitochondrial Cox1 (cytochrome oxidase subunit 1) sequences were obtained from 76 specimens collected from seven sites across the Eastern Mediterranean Basin, stretching across the Greek and Turkish coasts. This newly generated dataset of sequences was merged with previously examined data by Ragionieri and Schubart [35] and Deli et al. [12] for phylogeographic re-analyses.


Sampling strategy and genomic DNA extraction

A total of 263 samples of Carcinus aestuarii were included in the present phylogeographic study, of which 76 were newly analyzed, and the remaining 187 had been previously investigated [12, 35] (Table 1 and Fig. 1). Newly examined specimens of C. aestuarii were collected from seven locations stretching across the Greek (Lefkada and Alexandroupolis) and Turkish coasts (Izmir Bay and Enez Dalyan Lagoon in the Aegean Sea; Dardanelles Strait, Prince’s Islands and Bosphorus Strait in the Sea of Marmara) (Fig. 1). Previously examined data by Ragionieri and Schubart [35] include the populations of Pomer (Croatia), Amvrakikos, and Navarino (both Greece) corresponding to the Adriatic and Ionian seas respectively; while those retrieved from Deli et al. [12] comprise all surveyed populations from North Africa and Venice Lagoon. In the present study, we included all previously investigated eastern Mediterranean populations in order to optimize the population genetic structure analysis across the poorly surveyed Eastern Mediterranean Basin, whereas from the available western Mediterranean locations we only included those from Deli et al. [12]. This strategy was adopted, as the primary aim of this study is the examination of the phylogeography and evolutionary history of eastern Mediterranean populations of C. aestuarii. Western Mediterranean locations are included as reference populations to re-analyze population genetic structure and test for genetic subdivision across the Siculo-Tunisian Strait. Those from the North-African coast were included in this study, because they allowed Deli et al. [12] to confirm and delineate the geographic break across the Siculo-Tunisan Strait based on the concordant patterns of mitochondrial and nuclear phylogeographic structure. This choice has been validated after confirming that incorporation of additional western Mediterranean sequences from Ragionieri and Schubart [35] did not affect the outcome of phylogeographic analyses and population structure. From each crab, muscle tissue was dissected from a removed pereiopod (after releasing the animal back into its original environment) and stored in absolute ethanol at − 20 °C until genetic analysis. Total genomic DNA was isolated from muscle tissue using the Wizard® genomic DNA purification kit (Promega), the Puregene kit (Gentra Systems: Minneapolis, MN55447, USA) or the Roche High Pure PCR Template Preparation Kit (Indianapolis, USA) following the instructions of the suppliers.

Table 1 Sampling information on the green crab Carcinus aestuarii including collection sites, countries, Mediterranean basins, regions, geographic coordinates, and number of examined specimens (N) per each location. Genetic diversity measures (including number of haplotype (Nh), number of polymorphic sites (Nps), haplotype (h) and nucleotide (π) diversities, and mean number of nucleotide differences (K)) and historical demographic results (inferred from: Tajima’s D test (D), Fu’s F S test (F S ), Ramos-Onsins and Rozas’s R2 test (R2), and mismatch distribution raggedness index (rg)) were also provided for each investigated location and the total dataset
Fig. 1

Sampling locations of the green crab Carcinus aestuarii across the Mediterranean Sea. Distribution patterns and proportions of Cox1 types (I, II, and III) along the examined locations are shown in couloured circles. S-T S: Siculo-Tunisian Strait; PHB: Peloponnese Hydrographic Break (represented by the quasi-circular anti-cyclonic feature southwest of Peloponnese). The base map was constructed with the software DIVA-GIS 7.5.0 (

Cox1 gene amplification and sequencing

The decapod primers COL6a and COH6 (specifically designed for decapod crustaceans; see [46]) or the universal primers LCO1490 and HCO2198 [47] were used to amplify parts of the mitochondrial Cox1 gene. The adopted PCR mixture and thermocycling conditions were detailed in Deli et al. [12]. After being loaded on a 1.5% agarose gel and visualized under UV light, strong products were outsourced for sequencing with primer COL6a or HCO2198 to LGC Genomics (Berlin) or Macrogen Europe (Netherlands). The obtained sequences were visually inspected with Chromas Lite 2.1.1 [48], aligned with Clustal W as implemented in BioEdit [49], and trimmed to a 587 bp fragment for subsequent analyses. Sequences of Cox1 haplotypes, corresponding to the retrieved three haplogroups, were submitted to GenBank (accession numbers: MG798798-MG798903).

Statistical analyses

Intra-population genetic diversity and detection of selection signatures

We assessed the nucleotide composition of the analysed Cox1 fragment with MEGA version 7.0.18 [50]. In order to assess genetic diversity for each population as well as for the total dataset, we computed the number of haplotypes (Nh), number of polymorphic sites (Nps), haplotype diversity (h; [51]), nucleotide diversity (π; [51, 52]), and mean number of nucleotide differences (K) using DnaSP version 5.10 [53].

To test in how far natural selection is operating on the C. aestuarii Cox1 gene, we used the codon-based Z-test of selection for analysis, averaging over all sequence pairs, as implemented in MEGA version 7.0.18 [50]. Codon-based models of molecular evolution are able to infer signatures of selection from alignments of homologous sequences by estimating the relative rates of synonymous (dS) and nonsynonymous substitutions (dN) [54]. The relative abundance of synonymous and nonsynonymous substitutions within the examined gene sequences were computed, averaged, and compared for each population as well as for the total dataset. The rejection of the null hypothesis of strict neutrality (H0: dN = dS) depends on the tested alternative hypothesis (HA): (a) lack of neutral evolution (dN is different from dS: dN ≠ dS), (b) positive selection (dN > dS) or (c) purifying selection (dN < dS). A two-tailed test was implemented to reject neutral evolution, while one-tailed tests were used to check for positive and purifying selection, respectively. The variance of the difference between synonymous and nonsynonymous substitutions was computed using the bootstrap method (1000 pseudoreplicates). Analyses were conducted using the Nei and Gojobori [55] procedure.

Intraspecific evolutionary relationships among Cox1 haplotypes

A statistical parsimony network, constructed with the software TCS version 1.21 [56] under the 95% probability criterion for a parsimonious connection [57, 58], allowed inferring intraspecific evolutionary relationships among the Cox1 haplotypes of C. aestuarii. Based on the outcome of the TCS analysis, we also examined the distribution pattern of the discerned divergent haplotypes (or Cox1 haplogroups) across the surveyed geographic region.

Relationship between phylogeny and the geographical distribution of haplotypes

Based on the particular pattern of evolutionary relationships among Cox1 haplotypes (yielding various differentiated haplogroups that were restricted to specific geographic regions), we assessed the relationship between phylogeny and the geographical distribution of the recorded haplotypes by measuring levels of population subdivision, using both unordered (GST) and ordered haplotypes (NST). Estimation and comparison of these parameters was based on the methods described by Pons and Petit [59, 60] using PERMUT & CPSRR version 2.0 [60]. If NST is significantly higher than GST, it usually indicates the presence of phylogeographic structure [60, 61].

Divergence estimation among resulting mitochondrial Cox1 haplogroups

In order to elucidate the evolutionary history of C. aestuarii across the surveyed geographic region, we estimated divergence time between the discerned mitochondrial haplogroups. For this purpose, we applied different models and calibration strategies in order to obtain a comprehensive estimate, minimizing uncertainties due to different model assumptions. First, we applied a known biogeographical event calibration using interspecific sequences (considering the speciation event between Mediterranean and Atlantic green crab species). Subsequently, the entire intraspecific dataset (all examined sequences of C. aestuarii in this study) was used for the divergence estimate, applying a specifically determined clock rate of the examined genetic marker for the genus Carcinus [23].

In a first analysis, we considered the closure of the Strait of Gibraltar at the beginning of the Messinian Salinity Crisis (5.59 million years ago; [8]) as calibration point for divergence estimation. Indeed, the completely interrupted contact between the Mediterranean Sea and the Atlantic Ocean, during the Messinian Salinity Crisis, is thought to provide the responsible geographic barrier for the speciation of the Mediterranean green crab C. aestuarii and its Atlantic sister species C. maenas (see [26, 62]). Divergence estimations were carried out in BEAST version 1.7.5 [63]. Prior to the analysis, the most appropriate model of sequence evolution for the dataset was selected using MODELTEST version 3.7 [64] based on Akaike Information Criterion scores. We included only unique haplotype sequences, corresponding to the encountered C. aestuarii haplogroups and the Atlantic sister species C. maenas respectively. Notably, for such kind of analysis, the simplest tree priors are the one parameter Yule model [65] and the two-parameter Birth-Death model [66, 67]. We used the latter model, since it has been suggested as an appropriate null model for species diversification [66]. In order to test for the right clock model, analyses were carried out first with a strict clock and repeated with an uncorrelated lognormal relaxed clock [68]. Since the parameter of the standard deviation of the uncorrelated lognormal relaxed clock was significantly different from zero (ucld.stdev = 0.28, 95% HPD: 1.13 10− 4 -0.54), highlighting variation in rates among branches, final analyses were run enforcing a relaxed molecular clock model. Uncertainty on the divergence time was modelled using a normal prior with a standard deviation of 55,000 years [23]. The normal distribution is considered a useful calibration prior when applying a biogeographical date.

In order to check the consistency of dating results, additional analyses were performed involving only the examined intraspecific data of C. aestuarii from this study and implementing the coalescent tree prior that is typically used when all the samples are from the same species [69]. Both strict molecular clock and lognormal relaxed molecular clock models were compared using Bayes factors (BF) to test which of these two clock models best fitted our intraspecific data. We used TRACER version 1.5 [70] to compare twice the difference in the marginal model posterior likelihoods (MMPLs) as estimated from the harmonic mean of the sample of posterior trees for each scenario. Values of 2ln (BF) > 10 were considered as very strong evidence for a given model to be more likely than another [71]. As the Bayes factors indicated a much better fit for the strict clock model (MMPL = − 9484.4) than for a relaxed clock model (MMPL = − 9984.3) (2 ln(BF) =12.428), the final analyses were carried out with a strict molecular clock, and assuming the generalised time reversible (GTR) model of sequence evolution [72], as calculated by MODELTEST version 3.7 [64]. The specifically estimated mutation rate for Carcinus of 3.86% per Myr (see [23]) was used to calibrate the genealogy and date tMRCA of Cox1 lineages. This complementary strategy was implemented to minimise errors of divergence times estimation inferred from the use of deep calibration points [73].

For all kinds of Bayesian analyses, the Markov chain Monte Carlo (MCMC) simulations were run for 100 million steps and sampled every 1000 steps. The corresponding outputs were reviewed in TRACER version 1.5 [70] for robustness, and the resultant trees were summarized in TreeAnnotator (implemented in BEAST). The final results are presented with FigTree version 1.4.0 [74].

Population genetic structure and phylogeographic examination

Overall population genetic differentiation (assessed by one-level Analysis of molecular variance [75]) as well as detailed pairwise comparison of genetic differentiation were estimated in ARLEQUIN version 3.1 [76], using the two fixation indices: ΦST (implementing the Tajima-Nei model, appropriate for unequal nucleotide frequencies [77]) and FST (based on haplotypic frequency). For both kinds of analyses, the resulting significant values were calculated from 10,000 permutations. B-Y FDR correction [78] was then applied to yield the exact level of significance (critical value = 0.00830 with 231 hypothesis tests and alpha = 0.05). In order to characterize patterns of genetic structure, i.e., identify differentiated genetic groups of populations, we performed a non-metric Multidimensional Scaling (MDS), through PAST version 2.17 [79], based on Tajima-Nei genetic distances. In order to test the hypothesis that patterns of genetic differentiation are caused by isolation by distance (IBD), we ran the Mantel test [80] for pairwise matrices between geographical and genetic distances. The Mantel test was performed with the software AIS (Alleles in Space) version 1.0 [81]. The statistical significance of the test was assessed by running 10,000 random permutations. Population genetic structure of C. aestuarii was also examined (by means of two-level AMOVAs) under various biogeographic hypotheses, testing the significance of population structure among Mediterranean basins (Western Mediterranean vs. Eastern Mediterranean), or among defined regions within basins (Algerian Basin vs. Afro-Sicilian Basin vs. Adriatic Sea vs. Ionian Sea vs. Aegean Sea vs. Sea of Marmara). Partitioning of C. aestuarii genetic variation was also assessed based on the outcome of haplotype network (among the main groups of populations defining each haplogroup) and pairwise comparisons of genetic differentiation, as well as the MDS plot.

The spatial analysis of molecular variance (SAMOVA) approach, implemented in SAMOVA version 1.0 [82], without a prior population structure definition, was also used to infer the likely number of hierarchical groups explaining most of the retrieved genetic structure within C. aestuarii. The software was run with 100 random initial conditions for 10,000 iterations, with number of tested groups (K) ranging from 2 to 8.

Demographic history

Three neutrality tests (Tajima’s D [83], Fu’s Fs [84], and Ramos-Onsins and Rozas’s R2 [85]) were used to assess deviation from neutrality, and examine the demographic history of the Mediterranean green crab. Both D and Fs indices were estimated in ARLEQUIN, while R2 statistic was calculated in DnaSP. The examination of deviation from neutrality by all three indices was based on 1000 coalescent simulations. A scenario of population expansion is likely supported by significantly negative D and Fs values as well as significant R2 (in small population sizes). The Harpending’s raggedness index rg [86] was also used to examine demographic changes in C. aestuarii according to the population expansion model implemented in ARLEQUIN. A total of 10,000 replicates allowed testing the significance of the rg index. The four above mentioned parameters (D, Fs, R2, and rg) were applied to each examined population, the overall dataset, as well as the genetically differentiated geographic groups (as inferred mainly by SAMOVA).

Since deviations from neutrality are usually caused by changes in effective population size, we also applied Bayesian Skyline plots (BSP) [87] to explore the magnitude of historical demographic events. Notably, this Bayesian approach could allow the inference of detailed and realistic population size function [88], and also yield accurate estimation of expansion events [89]. BSP plots were generated for the geographic groups retrieved by SAMOVA. Analyses were carried out in BEAST version 1.7.5 considering a GTR model (as already calculated by MODELTEST version 3.7) and a strict molecular clock (confirmed as the best model fitting the examined intraspecific data when compared with lognormal relaxed molecular clock using Bayes factors (BF) method). The specific mutation rate of 3.86% per Myr, as estimated for Carcinus by Marino et al. [23], was implemented in the analysis in order to date expansion event. Pattern of effective population size evolution through time was assessed taking into account a generation time of approximately two years in the green crab [90]. Two independent MCMC (each with 50,000,000 iterations) were carried out. Following the removal of the first 10% iterations (5,000,000) as burn-in, the remaining replicates were combined in LogCombiner [63] and summarized as BSPs after checking their convergence (Effective Sample Sizes (ESS) of all parameters > 200 for each group) in TRACER version 1.5.


Genetic diversity and detection of selection signatures

Sequences corresponding to the mtDNA Cox1 gene from 263 individuals and 22 locations of C. aestuarii were included in the analyses. Of these, 76 were newly obtained, proofread and aligned. The resulting alignment had to be trimmed to a length of 587 basepairs. A total of 97 nucleotide sites were variable, of which 58 were parsimony-informative. Nearly 40% of the examined sequences were unique and allowed the identification of 106 haplotypes (Fig. 2 and Table 1). The nucleotide composition of the analyzed fragment showed an A-T bias (C = 18.81%; T = 36.12%; A = 26.62%; G = 18.45%), which is typical for invertebrate mitochondrial DNA [91].

Fig. 2

TCS haplotype network of Carcinus aestuarii, based on the alignment of 587 bp of the mitochondrial gene Cox1, showing the relationships among the recorded haplotypes. Haplotype 2 corresponds to the ancestral haplotype. Small black circles correspond to missing (or hypothetical) haplotypes. Each line between two points represents one mutational step. Circle sizes depict proportions of haplotypes; the smallest corresponds to 1 and the largest to 67 individuals

Genetic diversity analyses of this mitochondrial dataset revealed high total haplotype (h = 0.912 ± 0.012) and nucleotide (π = 0.0211 ± 0.0009) diversities (Table 1). A high level of the mean number of nucleotide differences (K) was also inferred (K = 12.389). The values of this parameter ranged from 0.166 in the population of Benikhiar to 12.911 in Navarino (Table 1).

The codon-based Z-test of selection allowed the rejection of the null hypothesis of strict-neutrality (dN = dS) for all examined populations (except Benikhiar), as well as for the total dataset (Table 2). Notably, the null hypothesis of strict-neutrality was rejected because of purifying selection (dS-dN = 6.256, P = 0.000), as no significant positive selection was detected for all investigated populations of C. aestuarii as well as for the whole dataset (dN-dS = − 6.013, P = 1.000) (Table 2).

Table 2 Codon-based Z-test of selection for analysis averaging over all sequence pairs. Alternative hypotheses of neutrality, positive selection and purifying selection were tested for each population as well as for the total dataset of Carcinus aestuarii. For each population, examined specimens were assigned to the corresponding Cox1 type

Intraspecific evolutionary relationships among Cox1 haplotypes

The phylogeographic relationships among the 106 recorded haplotypes, as inferred by the TCS statistical parsimony procedure, revealed a remarkable divergence among three haplogroups. These haplogroups will from now on be referred to as types I, II (see also Deli et al. [12]), and III. They are all characterized by a star-like shape centred around three main haplotypes: 2, 29 and 71 respectively (Fig. 2). Type III was found to be the most divergent, separated by at least 28 mutations from type I, and 17 mutations from type II. The most frequent haplotype 2 was found in 67 individuals and in all populations except those from northern Greece (Alexandroupolis) and Turkey (Izmir Bay, Enez Dalyan Lagoon, Dardanelles Strait, Prince’s Islands and Bosphorus Strait) (Additional file: Table S1). Haplotype 2 was distinguished from haplotypes 29 and 71 by 21 and 35 mutational steps, corresponding to sequence divergence rates of 3.40% and 5.96%, respectively. The main haplotypes 29 and 71, representative of types II and III respectively, were separated by 24 mutational steps, corresponding to a rate of sequence divergence of 4.08%.

Relationship between phylogeny and the geographical distribution of haplotypes

While type I sequences proved to be present in almost all examined populations, the genetic distinctiveness of the other two haplogroups hints at the existence of a remarkable regional geographic structure within C. aestuarii across the Eastern Mediterranean Basin. Indeed, type II sequences were only found in Tunisian (Monastir, Chebba, Sfax and Djerba) and Libyan (Tripoli and Mosrata) populations. On the other hand, type III sequences were mostly found in specimens from Greece (with increasing frequency in Lefkada, Navarino, and Alexandroupolis) and Turkey (Figs. 2 and 3a) to the extent that Alexandroupolis and the Turkish populations do not contain Cox1 type I, present in all other examined populations (Fig. 3a). Hence, the exclusive pattern of distribution of Cox1 type III in the Aegean Sea and the Sea of Marmara (Figs. 2 and 3b) and its lacking in the remaining investigated distribution range supports the existence of a phylogeographic break across the Eastern Mediterranean. The likely existence of geographic separation was supported and confirmed by the outcome of PERMUT analyses. Calculations of NST (0.211) and GST (0.121) revealed that the NST value is significantly higher than the GST value (P < 0.05), inferring significant relationship between phylogeny and the geographical distribution of haplotypes, and indicating the existence of marked phylogeographic structure within the examined material.

Fig. 3

Distribution patterns of types I, II and III of Carcinus aestuarii Cox1 gene across western and eastern Mediterranean basins (a) as well as geographic regions (b) (as defined in Table 1)

Distribution patterns of the divergent Cox1 types across the surveyed geographic region and clinal variation analysis

All three recorded Cox1 types occur in the Eastern Mediterranean (Fig. 3a). This clearly contrasts with distribution patterns in the Western Mediterranean, where almost exclusively one genetic type (type I) prevails (as evidenced by the present study (Fig. 3a) and the former study by Ragionieri et al. [35]). It should be noted that no population harbours all three types. Both genetic types I and II have been shown to co-exist only in the Tunisian and Libyan populations; while types I and III are recorded together in the two Greek populations of Lefkada and Navarino (Fig. 3a). Of particular interest to this study is the pattern of geographic distribution and transition of Cox1 types I and III along the Greek coastline, which might reflect a marked longitudinal cline. The proportion of type I, for example, decreased notably from the westernmost population of Amvrakikos (exclusively present) to the easternmost location of Alexandroupolis (where it was lacking) (Fig. 3a). This can also be clearly observed at the regional level. Indeed, a gradual decrease in proportion of type I was shown from the Adriatic Sea to the Aegean Sea (Fig. 3b).

Divergence estimation among Cox1 haplogroups

The outcome of Bayesian phylogenetic analysis, as implemented in BEAST, confirmed the results inferred from the TCS parsimony procedure, yielding a noticeable separation between the highly differentiated type III and the monophyletic group composed of both types I and II. Assuming i) a 5.59 Myr split (see [8]) between C. aestuarii and C. maenas as calibration point and ii) a Birth-Death prior and uncorrelated lognormal relaxed clock, divergence among type III and the monophyletic group composed of types I and II was estimated to occur approximately around 1.54 Mya (95% HPD - high posterior density interval: 0.88–2.67 Mya). Using a strict molecular clock with a species-specific mutation rate of 3.86% per Myr (calculated and used by Marino et al. [23] for Carcinus), and assuming a GTR model of sequence evolution and a coalescent tree prior, involving all the intraspecific data of C. aestuarii, the estimated time of the split between Cox1 type III and both types I and II was relatively younger (0.8 Mya; 95% HPD: 0.54–1.05 Mya).

Population genetic structure and phylogeographic examination

The one-level AMOVA gives evidence for strong and highly significant genetic differentiation among the examined populations of C. aestuarii based on Tajima-Nei distances (ΦST = 0.615, P < 0.001) and haplotype frequencies (FST = 0.109, P < 0.001). This differentiation was more pronounced based on nucleotide divergence (more than 61% of the variation among populations). Pairwise comparisons of genetic differentiation, estimated from nucleotide divergence and haplotype frequencies, also yielded significant differences for most comparisons and revealed, in particular, a clear genetic distinctiveness of the populations of Alexandroupolis, Izmir Bay, Enez Dalyan Lagoon, Dardanelles Strait, Prince’s Islands and Bosphorus Strait, after B-Y FDR correction (Table 3). The outcome of pairwise comparisons between these latter populations did not reveal significant differences, highlighting the existence of pronounced genetic divergence between populations from the Aegean Sea and the Sea of Marmara, and those assigned to the remaining distribution area. This trend of separation was also confirmed by the outcome of a Multidimensional Scaling (MDS) analysis, based on Tajima-Nei distances (Fig. 4), suggesting a kind of parapatric genetic divergence among two adjacent geographic groups of C. aestuarii. A significant relationship was found between genetic and geographic distances (r = 0.189, P = 0.002) by means of a Mantel Test, supporting an isolation by distance hypothesis to better explain the population separation.

Table 3 Pairwise comparisons of genetic differentiation estimated from nucleotide divergence (ΦST, below the diagonal) and haplotype frequency (FST, above the diagonal). Asterisks indicate significant values (P < 0.05) calculated from 10.000 permutations. These significance values were subjected to a B-Y FDR correction [78], rendering a critical value of P < 0.0083; values that remain significant after the corection are shown in boldface
Fig. 4

Multidimensional scaling plot based on ΦST (Tajima-Nei distances) values between examined populations of Carcinus aestuarii

Population genetic structure was examined by means of two-level AMOVA, testing for partitioning of genetic variation under alternative biogeographic hypotheses (based on the geographic origin of examined populations, the outcome of haplotype network and the outcome of pairwise comparisons of genetic differentiation). Our results showed significant genetic structure under these various grouping schemes (Table 4). Indeed, besides the significant genetic subdivision across the Siculo-Tunisian Strait (ΦCT = 0.164, P < 0.05; FCT = 0.082, P < 0.01), a significant and more pronounced genetic separation within C. aestuarii (ΦCT = 0.598, P < 0.001; FCT = 0.079, P < 0.001) was also revealed when testing differentiation among Mediterranean sub-basins (biogeographic hypothesis 2; Table 4), defined according to Spalding’s marine ecoregions [92]. Notably, even if all AMOVA results for most of the tested alternative biogeographic hypotheses highlighted significant but relatively similar results, partitioning of the genetic variance among the two adjacent genetic groups (separating the Aegean and Marmara Seas from the remaining distribution area) yielded the highest ΦCT and FCT levels (ΦCT = 0.754, P < 0.001; FCT = 0.159, P < 0.001). It also explained most of the population genetic structure of C. aestuarii (more than 75% of genetic variance explained between groups, based on nucleotide divergence; and nearly 16% of inter-group variance discerned according to haplotype frequencies; Table 4). The outcome of spatial analysis of molecular variance (SAMOVA), proposing the number of population groups based on geographical and genetic distances without a prior assumption of group composition, showed that partitioning of variance among groups (ΦCT) was highest with two hierarchical groups (K = 2: ΦCT = 0.750, P < 0.001; Table 5). It should be highlighted that even if slight differences exist between the generated ΦCT values (inferred from the different predefined numbers of group (K)), the selected significant grouping of K = 2 (corresponding to the highest partitioning of variance among groups) matches perfectly the results inferred from the outcomes of pairwise comparisons of genetic differentiation and MDS plots. Besides, the concordance in the outcome of the two different approaches defining the population structure pattern (with (AMOVA) or without (SAMOVA) a prior structure parameter) consolidate the hypothesis that most of population genetic structure within C. aestuarii is explained assuming two hierarchical groups. Accordingly, this perfect agreement between results of different analyses assures the correct delineation of population structure and verifies and confirms the choice of the SAMOVA grouping of K = 2. Such pattern of genetic differentiation highlights the existence of a barrier to gene flow between two delineated (and adjacent) geographic groups within C. aestuarii in the Eastern Mediterranean Basin. The first group, harbouring mainly types I and II, covered mostly the North African coast as well as the Adriatic and Ionian seas; while the second group, including exclusively type III specimens, characterized the Aegean and Marmara seas.

Table 4 Analysis of molecular variance assessing population genetic structure and testing for partitioning of Carcinus aestuarii genetic variation under alternative biogeographic hypotheses
Table 5 Results of spatial analysis of molecular variance (SAMOVA), depicting patterns of population structure of Carcinus aestuarii for each predefined number of group (K)

Demographic history

Overall, the applied neutrality tests revealed significant deviations from mutation-drift equilibrium for a total of twelve populations (11 by the Tajima’s D, 4 by the Fu’s Fs, and 5 by the Ramos-Onsins and Rozas’s R2 test; Table 1) suggesting that these populations seem to have experienced a recent expansion event. However, considering the fact that both D and Fs tests have little power to detect departure from neutrality, when sample size is small (less than 15–20 individuals, which is the case for most local populations in the present study), more weight should be put on the outcome of the R2 test, which is more powerful than either D or Fs tests in case of small sample sizes. Hence, population expansion can be roughly inferred for five populations (Bizerte, Pomer, Amvrakikos, Alexandroupolis, and Enez Dalyan Lagoon); whereas it is uncertain in all the remaining cases. The negative and significant value of Fu’s F S test, together with the small and non-significant value of Harpending’s raggedness index rg, also revealed significant deviations from neutrality for the whole dataset (Table 1), consistent with a scenario of a sudden demographic expansion. Evidence of departure from mutation-drift equilibrium was also recorded for both delineated geographic groups (as identified by SAMOVA). Indeed, all examined neutrality tests resulted in significant values (with marked negative output for both Tajima’s D and Fu’s Fs), associated with small and non-significant values of the Harpending’s raggedness index, for group 1 (D = − 1.496, P = 0.025; Fu’s F S  = − 24.740, P = 0.000; R2 = 0.041, P = 0.035; rg = 0.021, P = 0.999) and group 2 (D = − 1.963, P = 0.008; Fu’s F S  = − 20.118, P = 0.000; R2 = 0.038, P = 0.001; rg = 0.018, P = 0.955).

Demographic history of the two geographic groups of C. aestuarii, delineated by SAMOVA, was also inferred and detailed from the coalescent approach of the BSP analysis. SAMOVA-group 1 (North African coast + Adriatic and Ionian seas) showed a relatively sudden and recent increase of effective population size over time, following a short phase of size decrease, preceded by quite a stationary period (Fig. 5a). This clearly contrasts with the outcome of BSP plot for SAMOVA-group 2 (Aegean and Marmara seas), yielding remarkably progressive increase in the effective population size (Fig. 5b). Overall, BSP results were well concordant with those inferred from neutrality tests and revealed that the expansion time occurred approximately at about 35,000 years ago (CI: 25,000–42,000 years ago) for SAMOVA-group 1 and about 51,000 years ago (CI: 42,500–69,000 years ago) for SAMOVA-group 2. Notably, the increase of effective size was much more pronounced in the SAMOVA-group 1 (Fig. 5a).

Fig. 5

Bayesian skyline plot for the two genetically defined groups of Carcinus aestuarii, as identified by SAMOVA; a: group 1; b: group 2. Populations, defining both geographic groups, are reported in the results and in Table 5 (K=2). BSP plots showing changes in effective population size (Ne) over time (measured in years before present). The thick solid line depicts the median estimate, and the margins of the blue area represent the highest 95% posterior density intervals


The results of the present study provide new and interesting insights into the genetic composition of C. aestuarii across its distribution range, unveiling the existence of a new haplogroup (type III) in the most north-eastern part of the Mediterranean Sea, in addition to two haplogroups (types I and II) previously described by Deli et al. [12]. Type III was previously recorded in four specimens of C. aestuarii from the Greek population of Navarino (Peloponnesus) and briefly discussed in Ragionieri and Schubart [35] but the corresponding Cox1 sequences are here analysed for the first time. Overall, concordant results inferred from phylogenetic relationships among the recorded haplogroups, the patterns of their distribution, pairwise comparisons of genetic differentiation (illustrated with MDS plots), and two-levels AMOVA revealed strong genetic divergence among two adjacent groups (type III versus types I-II), separating the Aegean Sea and the Sea of Marmara from the remaining examined populations. Spatial genetic structure, evidenced by SAMOVA, highlighted the existence of a barrier to gene flow between these two delineated geographic groups within C. aestuarii and stressing a sharp phylogeographic break across the Eastern Mediterranean.

The existence of these two genetically and geographically differentiated groups confirmed the assumption and previous results that C. aestuarii exhibits complex population structure across its distribution range. Earlier investigations by Ragionieri and Schubart [35] and Deli et al. [12] showed marked genetic differentiation between populations from the Eastern and Western Mediterranean basins. In particular, Deli et al. [12] provided evidence for a Pleistocene vicariant event in this species across the Siculo-Tunisian Strait, inferred from marked concordant patterns of mitochondrial and nuclear phylogeographic structure.

The remarkable coexistence of all three Cox1 types (I, II, and III) in the Eastern Mediterranean, being regionally differentiated, allowed us to advance the hypothesis that such particular pattern of mitochondrial diversity within C. aestuarii could be the footprint of a complex evolutionary history of the species within the Mediterranean Sea. Type III represents the most distinct haplogroup with high levels of nucleotide divergence. We suggest it to be a genetic isolate that might have survived different episodes of Pleistocene climate changes. During the Quaternary glacial periods, sea level regressions limited the biotic exchange through the Strait of Gibraltar [13] and the Siculo-Tunisian Strait [14]. These historic complete isolation events must have shaped the genetic structure of the marine Mediterranean fauna. In particular, the Mediterranean crustacean fauna has been postulated to originate by repeated or continuous colonization events with adaptation to specialized habitats and adaptive radiation [93]. This could have led to marked genetic differences and a high proportion of endemism in different parts of the Mediterranean Sea [4, 94]. Earlier results by Deli et al. [12] allowed recognition of type II Cox1 haplotypes in the south-western part of the Eastern Mediterranean Basin, postulated to correspond to an eastern Mediterranean endemic isolate, originating from climate oscillations during the Early Pleistocene. Now we can ascertain that this type II seems to be regionally restricted to the African coast of the Eastern Basin.

The particular geographic concentration of C. aestuarii with type III mtDNA in the Sea of Marmara and the adjacent Aegean Sea potentially suggests that this isolated haplogroup could have originated in the Black Sea and subsequently dispersed into the Aegean. In this context, we hypothesize that the highly divergent Cox1 type III could be the result of historical isolation mediated by the closure of the Bosphorus Strait during Pleistocene climate shifts [95]. These processes probably caused a total restriction of gene flow and led to genetic divergence between Black Sea populations and their Eastern Mediterranean counterparts. Later on, resumed biotic exchange, following the opening of the Bosphorus Strait 10,000 years ago [95], could have restored the gene flow between the Black Sea and the Eastern Mediterranean, which may partly explain the recorded contemporary genetic structure of the green crab. This advanced scenario for C. aestuarii echoes a similar explanation of the discerned phylogeographic pattern in other marine invertebrate from the Black Sea and Aegean Sea, i.e., the bivalve Mytilus galloprovincialis (see [96]). Based on the outcome of haplogroup distribution, Kalkan et al. [96] hypothesized that one of the refugial regions of M. galloprovincialis could be the Black Sea, with subsequent dispersal of one of the clades (essentially counterparts of type III in the green crab) into the Aegean (via the Black Sea water current), following the onset of connection between both basins. A similar scenario was also proposed for the littoral prawn Palaemon elegans genetic type III, with a potential Black Sea origin [11]. Recently, Fratini et al. [97] discerned a specific haplogroup in the marbled crab Pachygrapsus marmoratus confined to the Black Sea. The authors stressed on the importance of the biogeographic barrier between the Aegean and Black seas, susceptible of disrupting dispersal and gene flow for many marine species during past and present times and triggering intraspecific diversification in the Mediterranean.

This likely scenario is reinforced by the fact that the Black Sea was repeatedly isolated from the Mediterranean Sea during Pleistocene glaciations and diluted with fresh water during those periods [98]. Svitoch et al. [99] reported that the isolation event resulting in the highest freshwater condition occurred around 18,000–20,000 years ago when the salinity dropped to around 2–4%. Accordingly, we might think that it was impossible for marine invertebrates, i.e., decapods, to survive these conditions. Nevertheless, the Quaternary paleogeography of the Black Sea indicates that it was never an exclusively freshwater basin [100]. Notably, during the Early Pleistocene (around one million years ago) when type III had started evolving, the Black Sea environment was still brackish (during the Pontian and Chaudian epochs), and the green crab could have survived as a ‘Pontic relict’ species [100]. Taking into account these insights, and backing the possibility that C. aestuarii could also survive in hypo-saline conditions (given its wide ecological valence [19]), it can therefore make much sense to argue for a potential Black Sea origin of the type III mtDNA in the Mediterranean green crab. Reuschel et al. [11] previously suspected that the genetic type III of P. elegans could be a lineage that survived for extended periods in the Black Sea (for example during the Messinian Salinity Crisis). This lineage must be relatively tolerant to brackish waters, considering that it recently invaded the Baltic Sea (characterized by increasingly low salinity from west to east). Similarly, Luttikhuizen et al. [101] argued that other decapod such as the shrimp Crangon crangon might also have survived in the Black Sea in brackish waters with salinities of less than 7 ppt.

The marked genetic divergence between both geographically delineated groups (harbouring different Cox1 types) could reflect the impact and likely involvement of historical processes (i.e., Pleistocene climate oscillations; see [4]) in shaping Mediterranean marine diversity. In particular, the significant divergence of type III, being separated by large number of mutational steps from both types I and II, indicates a relatively old separation along the study area and suggests that these mitochondrial clades had been formed and delineated by long-term biogeographic barriers, and their differentiation was probably maintained by restricted historical gene flow. This assumption was confirmed by our dating procedure, based on the use of different models and calibration strategies, placing the divergence between type III and both types I and II at 1.54 Mya to 0.80 Mya. This range of divergence time estimation corresponds approximately to the Early Pleistocene (1.8 to ~ 0.781 Mya according to Riccardi [102]) and coincides with historical episodes of separation between types I and II (1.2 Mya to 0.69 Mya; [12]), providing evidence of simultaneous impact of climate change on genetic structure of the green crab across different parts of the Eastern Mediterranean. It is known that this period was characterized by strong climate shifts and sea-level fluctuations that might have profoundly affected the genetic structure of populations of several Mediterranean marine species (see [4, 103]). Furthermore, the marked gradual transition between type I and type III along the Adriatic, Ionian and Aegean seas may also reflect the impact of historical events in the surveyed region. Genetic clines, such as here observed, could originate from genetic admixture at secondary contact zones, following postglacial recolonization from Pleistocene refugia [104]. Accordingly, we hypothesize that the Ionian Sea could be considered a contact zone between the two defined parapatric divergent groups of C. aestuarii within the Eastern Mediterranean, following episodes of historical isolation between the Aegean Sea and the Eastern Mediterranean, as has been reported in several studies [5, 16, 39,40,41,42,43,44,45]. Further sampling along the eastern Greek coastline is required in order to test this hypothesis.

The three applied neutrality tests (Tajima’s D, Fu’s F S , and Ramos-Onsins and Rozas’s R2) for both divergent genetic groups of C. aestuarii (delineated by SAMOVA) showed significant values, indicating significant deviation from neutrality due to historical factors (such as a population bottleneck or sudden expansion) and contemporary processes (i.e., natural selection). Notably, significant negative D and F S values (interpreted as signatures of historical population expansion, supported by the small and non-significant value of Harpending’s raggedness index rg) were retrieved for both geographic groups. Such patterns clearly indicate that these genetic entities have undergone demographic expansion and refer to a loss of equilibrium among mutation, gene flow and genetic drift [83, 84]. Hence, historical processes, rather than contemporary ones, are supposed to be likely involved in triggering the onset of the retrieved phylogeographic pattern. This assumption could be supported by the outcome of Bayesian Skyline Plots analyses. Indeed, different patterns of demographic history were revealed for both groups of C. aestuarii, highlighted by more recent expansion event in group 1 (35,000 years ago) than that recorded in group 2 (51,000 years ago). Notably, while both demographic events were found to precede the Last Glacial Maximum (between 26,500–20,000 years ago; [105]), the temporal difference between both genetic groups of C. aestuarii points out to the potential intensity and effect of historical processes (i.e., paleoclimate fluctuations) on pattern of genetic diversity evolution. In addition to the observed temporal variation, marked spatial differences in the demographic history of the Mediterranean green crab also have been noticed. While a pattern of sudden population growth was detected for group 1 (composed mainly by types I and II), the species exhibited a pattern of slight and progressive increasing population size in the Aegean Sea and the Sea of Marmara (group 2 harbouring exclusively type III). We may therefore hypothesize that these spatial differences could likely stem from the potential impact of Pleistocene climate fluctuations (alternating glacial-interglacial cycles) on the availability of favourable abiotic conditions (such as suitable temperature and salinity conditions, as well as habitat availability highlighted by suitable ecological niche). In this context, we may attribute the less magnitude of expansion event, recorded in the Aegean and Marmara seas (group 2), to the lack of suitable habitat availability for C. aestuarii which could have probably led to population expansion limitation. Indeed, during glacial cycles of the Pleistocene, particularly during the Middle Pleistocene and the last glacial period, the Aegean Sea experienced massive sea-level regression that probably caused half of its area to become land [106, 107]. Moreover, even after the rising of sea level, at the end of the Pleistocene glaciations, the post-glacial recolonization of the Aegean Sea could not have led to considerable increase in the effective population size of C. aestuarii owing to the limited areas of the Aegean and Marmara seas.

Alternatively, the recorded pattern of parapatric genetic divergence among the two adjacent groups, with the potential existence of contact zone among both groups in the Ionian Sea, suggest the impact of particular past and present oceanographic patterns (such as marine currents and gyres) across the Eastern Mediterranean. Accordingly, we assume that the recorded pattern of population structure could be mainly caused by the potential effects of the hydrographic isolation of the Adriatic-Ionian and Aegean seas [4, 5] on larval dispersal. In addition to the impact of the strong currents impeding the mixing of the different water bodies at each sub-region, selective forces associated with environmental features in each basin could account for a phylogeographic break. Despite the existence of significant genetic differentiation across the Siculo-Tunisian Strait, as already unveiled by previous investigations [12, 23, 35], the partition maximizing genetic variance among groups (more than 75%) was recorded between the Aegean-Marmara seas and the remaining group of populations (as also confirmed by pairwise comparisons of genetic differentiation, MDS plot and SAMOVA). These results are in concordance with those inferred from other studies on Mediterranean marine invertebrates and vertebrates (see [5, 16, 39,40,41,42,43,44,45, 108]), corroborating the isolation of the Aegean Sea from the remaining Mediterranean Sea. In particular, Nikula and Väinölä [39] unravelled genetic subdivision within the bivalve Cerastoderma glaucum in the Eastern Mediterranean, highlighting a major phylogeographic break in mitochondrial Cox1 gene sequences between a group of Ponto-Caspian-Aegean Sea haplotypes and those to the west of the Peloponnese peninsula in the Mediterranean. This remarkable pattern of genetic break, within the Eastern Mediterranean, has also been recorded in crustaceans. Notably, Shemesh et al. [108] and later Pannacciulli et al. [18] revealed marked genetic isolation in the barnacle species Chthamalus montagui between the western-central Mediterranean Sea and Aegean-Black seas. Gene flow restriction between western and eastern Mediterranean populations was shown to be mainly linked to the hydrographic isolation of the Aegean Sea, and found to be more marked than that associated with the Siculo-Tunisian Strait [45].

Genetic divergence within C. aestuarii across the Eastern Mediterranean could also be triggered by historical isolation events resulting in strong bottlenecks. Such isolation events could include Pleistocene hydrographical shifts that allowed repeated isolation and separation of the Aegean Sea from the remaining Eastern Mediterranean [38]. The continuous isolation of both delineated geographic groups (as inferred from SAMOVA) until the present is likely maintained by the Peloponnese anticyclonic front [109,110,111], which would have prevented gene flow between both geographic groups. What is more, the exclusive existence of type III in the Aegean Sea and Sea of Marmara (with the consequent total lack of types I and II) matches perfectly with the specific oceanographic features of the Sea of Marmara-Aegean Sea area (due to the low salinity surface water mass from the Black Sea flowing from the Black Sea to the Aegean on top of the denser saline Mediterranean waters flowing towards the Black Sea), which should prevent larvae flow from the Mediterranean Sea into the Aegean Sea [39, 45, 112]. Limited genetic connectivity was only revealed between the Ionian and Aegean seas, with occasional dispersal events from the latter area to the former. For instance, type III (predominant in the Aegean Sea) was otherwise only found in the very adjacent locations of Lefkada and Navarino (Ionian Sea) and in lower frequencies, indicating a unidirectional dispersal event.

In addition to the previously discussed impact of palaeoecological history as well as past and present oceanographic processes on shaping the genetic variability and population structure of the green crab C. aestuarii, selection on mtDNA haplotypes is another important factor that should be taken into consideration. For instance, the significant genetic patterns, linked to geography within the examined Mediterranean coasts, may strongly suggest that gene flow is not only limited by fluctuating events and neutral processes (as suggested by the observed significant pattern of isolation by distance), but also by environmental factors, i.e., hydrological factors (affecting dispersal) or differential selection (affecting fitness). Reduced dispersal among populations can lead to genetic subdivision of populations and may facilitate local adaptation [113]. The impact of selection on genetic structure has already been suggested in marine invertebrates [114, 115]. Therefore, we may attribute genetic distinctiveness of the eastern group to the specific abiotic features of the Aegean Sea. Water temperatures in the Aegean are known to be influenced by the cold-water masses of low temperature that flow in from the Black Sea through the Dardanelles Straits [116]. The sea surface temperature in the Aegean generally ranges from about 16 to 25 °C. Furthermore, hypersaline Mediterranean water (moving northward along the west coast of Turkey) characterizes the Aegean surface water before being displaced by less dense Black Sea outflow [117]. Hence, we hypothesize that specific environmental features of the Aegean Sea might have exerted selective pressures on the gene pool of C. aestuarii in the Eastern Mediterranean. A significant correlation between sea surface temperature and mitochondrial haplotypes has been already recorded in marine species, i.e., the walleye pollock Theragra chalcogramma in the North Pacific [118]. Nevertheless, with no detected sign of positive natural selection in examined Cox1 sequences of C. aestuarii, this hypothesis remains questionable and would need to be verified.

Overall, the results of the present investigation, along with those already obtained for more western populations of the Mediterranean green crab [12, 35], allow us to postulate the following evolutionary history scenario for C. aestuarii throughout its distribution range: During glaciations periods of the Early Pleistocene, dropping sea levels led to the restriction of biotic exchange across the Gibraltar and Siculo-Tunisian straits. The Eastern Mediterranean Basin was more affected by these environmental shifts and experienced desiccation episodes of greater or lesser importance in different parts [119]. Hence, being isolated from the rest of the Mediterranean, an endemic fauna of the Eastern Mediterranean may have originated and evolved a different genetic composition. In this context, climate oscillations during the Pleistocene may have contributed to the simultaneous onset of different genetic isolates (Cox1 Types II and III) in different parts of the Eastern Mediterranean (eastern Mediterranean endemic isolates). Later, the relative impact of historical and contemporary barriers to gene flow, as well as different patterns of postglacial recolonization of C. aestuarii from Pleistocene refugia in different parts of the distribution area, might have shaped the current genetic diversity and population genetic structure, such as observed by us. Indeed, in an earlier study by Deli et al. [12], a secondary contact between historically isolated types I and II was noticed in the central Mediterranean (i.e., across the African eastern Mediterranean). The outcome of the present study confirms these earlier insights and revealed further separation between both genetic types (I and II) and the highly diverged type III in the Eastern Mediterranean (see Fig. 1).


Our study provides new and important insights into the evolutionary history of a highly dispersive benthic decapod crustacean in the Mediterranean Sea. Notably, results of this investigation allow unravelling a sharp phylogeographic break in the Eastern Mediterranean (matching the well-known and reported sharp genetic break between the Aegean Sea and the remaining Mediterranean) and elucidating historical and contemporary processes susceptible of driving such complex phylogeographic structure. Our finding also stress the importance of investigating peripheral areas in the species’ distribution zone to fully understand the distribution of the genetic diversity and unravel hidden genetic units and local patterns of endemism. Lastly, regardless the mechanism involved in shaping pattern of genetic structure of C. aestuarii, the two discerned geographic groups deserve to be better investigated. In this sense, analysis of nuclear markers (such as microsatellite loci) in areas where both genetic groups occur is required to confirm this particular divergence pattern. It should be noted that the use of such additional markers may provide changes in phylogeographic patterns and consequent interpretations, as reported by Avise [120], yet still providing complete picture on phylogeographic structure and evolutionary history of the species. In addition, further studies including populations from the Black Sea and those located to the south of the Aegean Sea (i.e. the Levantine Basin) are needed to better understand the evolutionary history of the Mediterranean green crab, and fully characterize and delineate other potential genetic groups or isolates. In addition, further sampling across the Ionian Sea would allow confirming and precisely delineating the geographic occurrence of the observed genetic cline in this study.



Degrees celsius


Deoxyribonucleic acid


Mitochondrial DNA


Million years ago


Million years


Polymerase chain reaction


Parts per thousand


Ultra violet


  1. 1.

    Maggs CA, Castilho R, Foltz D, Henzler C, Jolly MT, Kelly J, et al. Evaluating signatures of the glacial refugia for North Atlantic benthic marine taxa. Ecology. 2008;89:108–22.

    Article  Google Scholar 

  2. 2.

    Berg PR, Jentoft S, Star B, Ring KH, Knutsen H, Lien S, et al. Adaptation to low salinity promotes genomic divergence in Atlantic cod (Gadus morhua L.). Genome biol. Evolution. 2015;7(6):1644–63.

    CAS  Google Scholar 

  3. 3.

    Godhe A, Egardt J, Kleinhans D, Sundqvist L, Hordoir R, Jonsson PR. Seascape analysis reveals regional gene flow patterns among populations of a marine planktonic diatom. Proc Royal Soc B. 2013;280:20131599.

    Article  Google Scholar 

  4. 4.

    Patarnello T, Volckaert FAM, Castilho R. Pillars of Hercules: is the Atlantic-Mediterranean transition a phylogeographic break? Mol Ecol. 2007;16:4426–44.

    Article  PubMed  Google Scholar 

  5. 5.

    Pérez-Losada M, Nolte MJ, Crandall KA, Shaw PW. Testing hypotheses of population structuring in the Northeast Atlantic Ocean and Mediterranean Sea using the common cuttlefish Sepia officinalis. Mol Ecol. 2007;16:2667–79.

    Article  PubMed  Google Scholar 

  6. 6.

    Hewitt GM. Genetic consequences of climatic oscillations in the quaternary. Phil Trans R Soc B. 2004;359:183–95.

    CAS  PubMed Central  Article  PubMed  Google Scholar 

  7. 7.

    Galarza JA, Carreras-Carbonell J, Macpherson E, Pascual M, Roques S, Turner GF, et al. The influence of oceanographic fronts and early-life-history traits on connectivity among littoral fish species. Proc Natl Acad Sci U S A. 2009;106:1473–8.

    CAS  PubMed Central  Article  PubMed  Google Scholar 

  8. 8.

    Krijgsman W, Hilgen FJ, Raffi I, Sierro FJ, Wilson DS. Chronology, causes and progression of the Messinian salinity crisis. Nature. 1999;400:652–5.

    CAS  Article  Google Scholar 

  9. 9.

    Lambeck K, Esat TM, Potter EK. Links between climate and sea levels for the past three million years. Nature. 2002;419:199–205.

    CAS  Article  PubMed  Google Scholar 

  10. 10.

    Pannacciulli FG, Bishop JDD, Hawkins SJ. Genetic structure of populations of two species of Chthamalus (Crustacea: Cirripedia) in the Northeast Atlantic and Mediterranean. Mar Biol. 1997;128:73–82.

    Article  Google Scholar 

  11. 11.

    Reuschel S, Cuesta JA, Schubart CD. Marine biogeographic boundaries and human introduction along the European coast revealed by phylogeography of the prawn Palaemon elegans. Mol Phyl Evol. 2010;55:765–75.

    Article  Google Scholar 

  12. 12.

    Deli T, Chatti N, Said K, Schubart CD. Concordant patterns of mtDNA and nuclear phylogeographic structure reveal Pleistocene vicariant event in the green crab Carcinus aestuarii across the Siculo-Tunisian Strait. Mediterr Mar Sci. 2016;17(2):533–51.

    Article  Google Scholar 

  13. 13.

    Vermeij GJ. Biogeography and adaptation. Patterns of marine life. Cambridge: Massachusetts and London: Harvard University Press; 1987.

    Google Scholar 

  14. 14.

    Thiede J. A glacial Mediterranean. Nature. 1978;276:680–3.

    Article  Google Scholar 

  15. 15.

    Arnaud-Haond S, Diaz Almela E, Teixeira S. Vicariance patterns in the Mediterranean Sea: east-west cleavage and low dispersal in the endemic seagrass Posidonia oceanica. J Biogeogr. 2007;14:963–76.

    Article  Google Scholar 

  16. 16.

    Sanna D, Cossu P, Dedola GL, Scarpa F, Maltagliati F, Castelli A, et al. Mitochondrial DNA reveals genetic structuring of Pinna nobilis across the Mediterranean Sea. PLoS One. 2013;8(6):e67372.

    CAS  PubMed Central  Article  PubMed  Google Scholar 

  17. 17.

    Cordero D, Peña JB, Saavedra C. Phylogeographic analysis of introns and mitochondrial DNA in the clam Ruditapes decussatus uncovers the effects of Pleistocene glaciations and endogenous barriers to gene flow. Mol Phylogenet Evol. 2014;71:274–87.

    CAS  Article  PubMed  Google Scholar 

  18. 18.

    Pannacciulli FG, Maltagliati F, de Guttry C, Achituv Y. Phylogeography on the rocks: the contribution of current and historical factors in shaping the genetic structure of Chthamalus montagui (Crustacea, Cirripedia). PLoS One. 2017;12(6):e0178287.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  19. 19.

    Mori M, Mancon R, Fanciulli G. Notes on the reproductive biology of Carcinus aestuarii Nardo (Crustacea, Decapoda) from the lagoon of san Teodoro (island of Sardinia, Italy). Riv Idrobiol. 1990;29:763–74.

    Google Scholar 

  20. 20.

    Behrens Yamada S, Hauck L. Field identification of the European green crab species: Carcinus maenas and Carcinus aestuarii. J Shellfish Res. 2001;20:905–12.

    Google Scholar 

  21. 21.

    Mistri M, Rossi R, Fano EA. Structure and secondary production of a soft bottom macrobenthic community in a brackish lagoon (Sacca di Goro, northeastern Italy). Estuar Coast Shelf Sci. 2001;52:605–16.

    Article  Google Scholar 

  22. 22.

    Darling JA, Bagley JM, Roman J, Tepolt CK, Geller JB. Genetic patterns across multiple introductions of the globally invasive crab genus Carcinus. Mol Ecol. 2008;17:4992–5007.

    CAS  Article  PubMed  Google Scholar 

  23. 23.

    Marino IAM, Pujolar JM, Zane L. Reconciling deep calibration and demographic history: Bayesian inference of post glacial colonization patterns in Carcinus aestuarii (Nardo, 1847) and C. maenas (Linnaeus, 1758). PLoS One. 2011;6(12):e28567.

    CAS  PubMed Central  Article  PubMed  Google Scholar 

  24. 24.

    Almaca C. Sur la distribution géographique du genre Carcinus Leach (Crustacea, Decapoda, Brachyoura). Revista Fac Ci Univ Lisboa. 1962;10:109–13.

    Google Scholar 

  25. 25.

    Furota T, Watanabe S, Watanabe T, Akiyama S, Kinoshita K. Life history of the Mediterranean green crab, Carcinus aestuarii Nardo, in Tokyo Bay, Japan. Crustacean Res. 1999;28:5–15.

    Article  Google Scholar 

  26. 26.

    Geller JB, Walton E, Grosholz E, Ruiz GM. Cryptic invasions of the crab Carcinus detected by molecular phylogeny. Mol Ecol. 1997;6:256–62.

    Article  Google Scholar 

  27. 27.

    Watts RJ, Johnson MS. Estuaries, lagoons and embayments: habitats that enhance population subdivision in fishes. Mar Freshw Res. 2004;55:641–51.

    Article  Google Scholar 

  28. 28.

    Paula J. Larval retention and dynamics of the prawns Palaemon longirostris H. Milne Edwards and Crangon crangon Linnaeus (Decapoda, Caridea) in the Mira estuary, Portugal. Invertebr Reprod Dev. 1998;33:221–8.

    Article  Google Scholar 

  29. 29.

    Lundquist CJ, Oldman JW, Lewis MJ. Predicting suitability of cockle Austrovenus stutchburyi restoration sites using hydrodynamic models of larval dispersal. New Zeal J Mar Fresh. 2009;43:735–48.

    Article  Google Scholar 

  30. 30.

    Largier JL. Estuarine fronts: how important are they? Estuaries. 1993;16:1–11.

    Article  Google Scholar 

  31. 31.

    Cognetti G, Maltagliati F. Biodiversity and adaptive mechanisms in brackish water fauna. Mar Pollut Bull. 2000;40:7–14.

    CAS  Article  Google Scholar 

  32. 32.

    Ward RD, Woodwark M, Skibinski DOF. A comparison of genetic diversity levels in marine, freshwater, and anadromous fishes. J Fish Biol. 1994;44(2):213–32.

    Article  Google Scholar 

  33. 33.

    Pelc RA, Warner RR, Gaines SD. Geographical patterns of genetic structure in marine species with contrasting life histories. J Biogeogr. 2009;36:1881–90.

    Article  Google Scholar 

  34. 34.

    Marino IAM, Barbisan F, Gennari M, Giomi F, Beltramini M, Bisol PM, et al. Genetic heterogeneity in populations of the Mediterranean shore crab, Carcinus aestuarii (Decapoda, Portunidae), from the Venice lagoon. Estuar Coast Shelf Sci. 2010;87:135–44.

    CAS  Article  Google Scholar 

  35. 35.

    Ragionieri L, Schubart CD. Population genetics, gene flow and biogeographic boundaries of Carcinus aestuarii (Crustacea: Brachyura: Carcinidae) along the European Mediterranean coast. Biol J Linn Soc. 2013;109:771–90.

    Article  Google Scholar 

  36. 36.

    Deli T, Said K, Chatti N. Genetic differentiation among populations of the green crab Carcinus aestuarii (Brachyura, Carcinidae) from the eastern and western Mediterranean coasts of Tunisia. Acta Zool Bulg. 2015;67(3):327–35.

    Google Scholar 

  37. 37.

    Deli T, Said K, Chatti N. Morphological differentiation among geographically close populations of the green crab Carcinus aestuarii Nardo, 1847 (Brachyura, Carcinidae) from the Tunisian coast. Crustaceana. 2014;87(3):257–83.

    Article  Google Scholar 

  38. 38.

    Svitoch AA, Selivanov AO, Yanina TA. The Pont-Caspian and Mediterranean basins in the Pleistocene (paleogeography and correlation). Oceanology. 2000;40:868–81.

    Google Scholar 

  39. 39.

    Nikula R, Väinölä R. Phylogeography of Cerastoderma glaucum (Bivalvia: Cardiidae) across Europe: a major break in the eastern Mediterranean. Mar Biol. 2003;143:339350.

    Article  Google Scholar 

  40. 40.

    Costagliola D, Robertson DR, Guidetti P. Evolution of coral reef fish Thalassoma pavo spp. (Labridae). 2. Evolution of the eastern Atlantic species. Mar Biol. 2004;144:377–83.

    Article  Google Scholar 

  41. 41.

    Domingues VS, Bucciarelli G, Almada VC, Bernardi G. Historical colonization and demography of the Mediterranean damselfish, Chromis chromis. Mol Ecol. 2005;14:4051–63.

    Article  PubMed  Google Scholar 

  42. 42.

    Chevolot M, Hoarau G, Rijnsdorp AD, Stam W, Olsen J. Phylogeography and population structure of thornback rays (Raja clavata L., Rajidae). Mol Ecol. 2006;15:3693–705.

    CAS  Article  PubMed  Google Scholar 

  43. 43.

    Peijnenburg KT, Fauvelot C, Breeuwer JA, Menken SB. Spatial and temporal genetic structure of the planktonic Sagitta setosa (Chaetognatha) in European seas as revealed by mitochondrial and nuclear DNA markers. Mol Ecol. 2006;15:3319–38.

    CAS  Article  PubMed  Google Scholar 

  44. 44.

    Zulliger DE, Tanner S, Ruch M, Ribi G. Genetic structure of the high dispersal Atlanto-Mediterranean Sea star Astropecten aranciacus revealed by mitochondrial DNA sequences and microsatellite loci. Mar Biol. 2009;156:597–610.

    CAS  Article  Google Scholar 

  45. 45.

    Borrero-Pérez GH, Gonzalez-Wangüemert M, Marcos C, Pérez-Ruzafa A. Phylogeography of the Atlanto-Mediterranean Sea cucumber Holothuria (Holothuria) mammata: the combined effects of historical processes and current oceanographical pattern. Mol Ecol. 2011;20:1964–75.

    Article  PubMed  Google Scholar 

  46. 46.

    Schubart CD. Mitochondrial DNA and decapod phylogenies: the importance of pseudogenes and primer optimization. In: Martin JW, Crandall KA, Felder DL, editors. Crustacean issues 18: decapod crustacean Phylogenetics. Boca Raton, Florida: Taylor and Francis/CRC Press; 2009. p. 47–65.

    Chapter  Google Scholar 

  47. 47.

    Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol. 1994;3(5):294–9.

    CAS  PubMed  Google Scholar 

  48. 48.

    Technelysium Pty Ltd. Chromas lite ver. 2.1.1. 2012. Available from:

  49. 49.

    Hall TA. BioEdit: a user-friendly biological sequence alignment editor and analysis program for windows 95/98/NT. Nucleic Acids Symp Ser. 1999;41:95–8.

    CAS  Google Scholar 

  50. 50.

    Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetic analysis version 7.0 for bigger datasets. Mol Biol Evol. 2016;33:1870–4.

    CAS  Article  PubMed  Google Scholar 

  51. 51.

    Nei M. Molecular evolutionary genetics. New York: Columbia University Press; 1987.

    Google Scholar 

  52. 52.

    Tajima F. Evolutionary relationships of DNA sequences in finite populations. Genetics. 1983;105:437–60.

    CAS  PubMed Central  PubMed  Google Scholar 

  53. 53.

    Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.

    CAS  Article  PubMed  Google Scholar 

  54. 54.

    Poon AF, Frost SD, Pond SL. Detecting signatures of selection from DNA sequences using Datamonkey. Methods Mol Biol. 2009;537:163–83.

    CAS  Article  PubMed  Google Scholar 

  55. 55.

    Nei M, Gojobori T. Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol Biol Evol. 1986;3:418–26.

    CAS  PubMed  Google Scholar 

  56. 56.

    Clement M, Posada D, Crandall K. TCS: A computer program to estimate gene genealogies. Mol Ecol. 2000;9:1657–60.

    CAS  Article  PubMed  Google Scholar 

  57. 57.

    Hudson RR. Gene genealogies and the coalescent process. In: Futuyama D, Antonovics JD, editors. Oxford surveys in evolutionary biology. Oxford: Oxford University Press; 1990. p. 1–44.

    Google Scholar 

  58. 58.

    Templeton AR, Crandall KA, Sing CF. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics. 1992;132:619–33.

    CAS  PubMed Central  PubMed  Google Scholar 

  59. 59.

    Pons O, Petit RJ. Estimation, variance and optimal sampling of gene diversity. Theor Appl Genet. 1995;90:462–70.

    CAS  Article  PubMed  Google Scholar 

  60. 60.

    Pons O, Petit RJ. Measuring and testing genetic differentiation with ordered versus unordered alleles. Genetics. 1996;144:1237–45.

    CAS  PubMed Central  PubMed  Google Scholar 

  61. 61.

    Petit RJ, Hampe A, Cheddadi R. Climate changes and tree phylogeography in the Mediterranean. Taxon. 2005;54:877–85.

    Article  Google Scholar 

  62. 62.

    Demeusy N. Recherches sur la mue de puberté du décapode brachyoure Carcinus maenas Linné. Arch Zool Exp Gen. 1958;95:253–492.

    Google Scholar 

  63. 63.

    Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969–73.

    CAS  PubMed Central  Article  PubMed  Google Scholar 

  64. 64.

    Posada D, Crandall KA. Modeltest: testing the model of DNA substitution. Bioinformatics. 1998;14:817–8.

    CAS  Article  PubMed  Google Scholar 

  65. 65.

    Yule GU. A mathematical theory of evolution, based on the conclusions of Dr. J. C. Wills, F. R. S. Philos Trans R Soc Lond B. 1924;213:21–87.

    Article  Google Scholar 

  66. 66.

    Nee S, May RM, Harvey PH. The reconstructed evolutionary process. Philos Trans R Soc Lond B. 1994;344:305–11.

    CAS  Article  Google Scholar 

  67. 67.

    Gernhard T. The conditioned reconstructed process. J Theor Biol. 2008;253:769–78.

    Article  PubMed  Google Scholar 

  68. 68.

    Drummond AJ, Ho SYW, Phillips MJ, Rambaut A. Relaxed phylogenetics and dating with confidence. PLoS Biol. 2006;4:e88.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  69. 69.

    Kingman JFC. The coalescent. Stoch Process Their Appl. 1982;13:235–48.

    Article  Google Scholar 

  70. 70.

    Rambaut A, Drummond AJ. Tracer v 1.4.8. Institute of evolutionary biology: University of Edinburg; 2007. Available from:

  71. 71.

    Kass RE, Raftery AE. Bayes factors. J Am Stat Assoc. 1995;90:773–95.

    Article  Google Scholar 

  72. 72.

    Tavaré S. Some probabilistic and statistical problems in the analysis of DNA sequences. Lectures on mathematics in the life sciences. Amer Math Soc. 1986;17:57–86.

    Google Scholar 

  73. 73.

    Burridge CP, Craw D, Fletcher D, Waters JM. Geological dates and molecular rates: fish DNA sheds light on time dependency. Mol Biol Evol. 2008;25:624–33.

    CAS  Article  PubMed  Google Scholar 

  74. 74.

    Rambaut A. FigTree v1.3.1. 2009. Computer program available at: http: //

  75. 75.

    Excoffier L, Smouse PE, Quattro JM. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics. 1992;131:479–91.

    CAS  PubMed Central  PubMed  Google Scholar 

  76. 76.

    Excoffier L, Laval G, Schneider S. Arlequin ver. 3.0: an integrated software package for population genetics data analysis. Evol Bioinformatics Online. 2005;(1):47–50.

  77. 77.

    Tajima F, Nei M. Estimation of evolutionary distance between nucleotide sequences. Mol Biol Evol. 1984;1:269–85.

    CAS  PubMed  Google Scholar 

  78. 78.

    Narum SR. Beyond Bonferroni: less conservative analyses for conservation genetics. Conserv Genet. 2006;7(5):783–7.

    CAS  Article  Google Scholar 

  79. 79.

    Hammer Ø, Harper DAT, Ryan PD. PAST: paleontological statistics software package for education and data analysis. Palaeontol Electron. 2001;4(1):1–9.

    Google Scholar 

  80. 80.

    Mantel N. The detection of disease clustering and a generalized regression approach. Cancer Res. 1967;27:209–20.

    CAS  PubMed  Google Scholar 

  81. 81.

    Miller MP. Alleles in space: computer software for the joint analysis of Interindividual spatial and genetic information. J Hered. 2005;96:722–4.

    CAS  Article  PubMed  Google Scholar 

  82. 82.

    Dupanloup I, Schneider S, Excoffier L. A simulated approach to define the genetic structure of populations. Mol Ecol. 2002;11:2571–81.

    CAS  Article  PubMed  Google Scholar 

  83. 83.

    Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123:585–95.

    CAS  PubMed Central  PubMed  Google Scholar 

  84. 84.

    Fu YX. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147:915–25.

    CAS  PubMed Central  PubMed  Google Scholar 

  85. 85.

    Ramos-Onsins SE, Rozas J. Statistical properties of new neutrality tests against population growth. Mol Biol Evol. 2002;19:2092–100.

    CAS  Article  PubMed  Google Scholar 

  86. 86.

    Harpending HC. Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum Biol. 1994;66(4):591–600.

    CAS  PubMed  Google Scholar 

  87. 87.

    Drummond AJ, Rambaut A, Shapiro B, Pybus OG. Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol. 2005;22:1185–92.

    CAS  Article  PubMed  Google Scholar 

  88. 88.

    Hulva P, Fornuskova A, Chudarkova A, Evin A, Allegrini B, Benda P, et al. Mechanisms of radiation in a bat group from the genus Pipistrellus inferred by phylogeography, demography and population genetics. Mol Ecol. 2010;19:5417–31.

    CAS  Article  PubMed  Google Scholar 

  89. 89.

    Grant WS. Problems and cautions with sequence mismatch analysis and Bayesian skyline plots to infer historical demography. J Heredity. 2015;106:333–46.

    Article  Google Scholar 

  90. 90.

    Yamada SB. Global invader: the European green crab. Corvallis: Oregon Seagrant; 2000.

    Google Scholar 

  91. 91.

    Simon C, Frati F, Beckenbach A, Crespi B, Liu H, Flook P. Evolution, weighting, and phylogenetic utility of mitochondrial gene sequences and a compilation of conserved polymerase chain reaction primers. Ann Entomol Soc Am. 1994;87:651–701.

    CAS  Article  Google Scholar 

  92. 92.

    Spalding MD, Fox HE, Allen GR, Davidson N, Ferdaña ZA, Finlayson M, et al. Marine ecoregions of the world: a bioregionalization of coast and shelf areas. Bioscience. 2007;57:573–83.

    Article  Google Scholar 

  93. 93.

    Almaça C. Evolutionary and zoogeographical remarks on the Mediterranean fauna of the brachyuran crabs. In: Moraitou-Apostolopoulou M, Kirostis V, editors. Mediterranean marine ecosystems. New York: Plenum Press; 1985. p. 347–66.

    Chapter  Google Scholar 

  94. 94.

    Hofricher AA. Das Mittelmeer-Fauna, Flora, Ökologie Band I. Berlin: Spektrum Akademischer Verlag, Heidelberg; 2002.

    Google Scholar 

  95. 95.

    Aksu AE, Hiscott RN, Mudie PJ, Rochon A, Kaminski MA, Abrajano T, et al. Persistent Holocene outflow from the Black Sea to the eastern Mediterranean contradicts Noah's flood hypothesis. GSA Today. 2002;12(5):4–10.

    Article  Google Scholar 

  96. 96.

    Kalkan E, Kurtuluş A, Maraci Ö, Bilgin R. Is the Bosphorus Strait a barrier to gene flow for the Mediterranean mussel, Mytilus galloprovincialis (Lamarck, 1819)? Mar Biol Res. 2011;7(7):690–700.

    Article  Google Scholar 

  97. 97.

    Fratini S, Ragionieri L, Deli T, Harrer A, Marino IAM, Cannicci S, et al. Unravelling population genetic structure with mitochondrial DNA in a notional panmictic coastal crab species: sample size makes the difference. BMC Evol Biol. 2016;16:150.

    PubMed Central  Article  CAS  PubMed  Google Scholar 

  98. 98.

    Fairbanks RG. A 17,000-year glacio-eustatic sea-level record: influence of global melting rates on the younger Dryas event and deep-ocean circulation. Nature. 1989;342:637–42.

    Article  Google Scholar 

  99. 99.

    Svitoch AA, Selivanov AO, Yanina TA. Paleohydrology of the Black Sea Pleistocene basins. Water Resour. 2000;27(6):594–603.

    CAS  Article  Google Scholar 

  100. 100.

    Zaitsev Y, Mamaev VO. Biological diversity in the Black Sea: a study of change and decline, Black Sea environmental series, vol. 3. New York: United Nations Publishing; 1997. p. 208.

    Google Scholar 

  101. 101.

    Luttikhuizen PC, van Bleijswijk JCJ, Peijnenburg KTCA, van der Veer HW. Phylogeography of the common shrimp, Crangon crangon (L.) across its distribution range. Mol Phylogenet Evol. 2008;46:1015–30.

    CAS  Article  PubMed  Google Scholar 

  102. 102.

    Riccardi AC. IUGS ratified ICS recommendation on redefinition of Pleistocene and formal definition of base quaternary. International Union of Geological Sciences; 2009.

    Google Scholar 

  103. 103.

    Hewitt GM. The genetic legacy of the quaternary ice ages. Nature. 2000;405:907–13.

    CAS  Article  PubMed  Google Scholar 

  104. 104.

    Barton NH, Hewitt GM. Analysis of hybrid zones. Annu Rev Ecol Syst. 1985;16:113–48.

    Article  Google Scholar 

  105. 105.

    Clark PU, Dyke AS, Shakun JD, Carlson AE, Clark J, Wohlfarth B, et al. The last glacial maximum. Science. 2009;325:710–4.

    CAS  Article  PubMed  Google Scholar 

  106. 106.

    Lykousis V. Sea-level changes and shelf break prograding sequences during the last 400 ka in the Aegean margins: subsidence rates and palaeogeographic implications. Cont Shelf Res. 2009;29:2037–44.

    Article  Google Scholar 

  107. 107.

    Perissoratis C, Conispoliatis N. The impacts of sea-level changes during latest Pleistocene and Holocene times on the morphology of the Ionian and Aegean seas (SE alpine Europe). Mar Geol. 2003;196:146–56.

    Article  Google Scholar 

  108. 108.

    Shemesh E, Huchon D, Simon-Blecher N, Achituv Y. The distribution and molecular diversity of the eastern Atlantic and Mediterranean chthamalids (Crustacea, Cirripedia). Zool Scr. 2009;38:365–78.

    Article  Google Scholar 

  109. 109.

    Malanotte-Rizzoli P, Bergamasco A. The circulation of the eastern Mediterranean. Part I Oceanol Acta. 1989;12:335–51.

    Google Scholar 

  110. 110.

    Roussenov V, Stanev E, Artale V, Pinardi N. A seasonal model of the Mediterranean Sea general circulation. J Geophys Res. 1995;100:13515–38.

    Article  Google Scholar 

  111. 111.

    Millot C. Circulation in the Mediterranean Sea: evidences, debates and unanswered questions. Sci. 2005;69:5–21.

    Google Scholar 

  112. 112.

    Kalkan E, Karhan SÜ, Bilgin R, Hemond EM. The Turkish straits system as a phylogeographic boundary – a literature review. In: Özsoy E, Çağatay MN, Balkıs NE, Nu B, Öztürk B, editors. The sea of Marmara; marine biodiversity, fisheries, conservation and governance. Istanbul, Turkey: Turkish Marine Research Foundation (TUDAV), Publication No:42; 2016. p. 550–69.

    Google Scholar 

  113. 113.

    Slatkin M. Gene flow and the geographic structure of natural populations. Science. 1987;236(4803):787–92.

    CAS  Article  PubMed  Google Scholar 

  114. 114.

    Rizzo C, Cammarata M, Di Carlo M, Pancucci A, Parrinello N. RAPD profiles distinguish Paracentrotus lividus populations living in a stressing environment (Amvrakikos gulf, Greece). Russ J Genet. 2009;45:499–503.

    CAS  Article  Google Scholar 

  115. 115.

    Penant G, Aurelle D, Feral J-P, Chenuil A. Planktonic larvae do not ensure gene flow in the edible sea urchin Paracentrotus lividus. Mar Ecol Prog Ser. 2013;480:155–70.

    Article  Google Scholar 

  116. 116.

    Skliris N, Sofianos SS, Gkanasos A, Axaopoulos P, Mantziafou A, Vervatis V. Long-term sea surface temperature variability in the Aegean Sea. Adv Oceanogr Limnol. 2011;2(2):125–39.

    Article  Google Scholar 

  117. 117.

    Aksu AE, Yaşar D, Mudie PJ, Gillespie H. Late Glacial-Holocene paleoclimatic and paleoceanographic evolution of the Aegean Sea: micropaleontological and stable isotopic evidence. Mar Micropaleontol. 1995;25(1):1–28.

    Article  Google Scholar 

  118. 118.

    Grant WS, Spies IB, Canino MF. Biogeographic evidence for selection on mitochondrial DNA in North Pacific walleye pollock Theragra chalcogramma. J Hered. 2006;97:571–80.

    CAS  Article  PubMed  Google Scholar 

  119. 119.

    Meijer PT, Krijgsman W. A quantitative analysis of the dessication and re-filling of the Mediterranean during the Messinian salinity crisis. Earth Planet Sci Lett. 2005;240:510–20.

    CAS  Article  Google Scholar 

  120. 120.

    Avise JC. Phylogeography: the history and formation of species. Cambridge: Harvard University Press; 2000.

    Google Scholar 

Download references


We would like to thank Daou Haddoud for his help with crab sampling from the Libyan coast, Valerio Matozzo for sending specimens from the Venice Lagoon, and Plamen Mitov for collecting materials from Alexandroupolis. This study was supported by a grant (No: 09S101) from the Research Fund of Boğaziçi University in Istanbul (Turkey) and a DAAD exchange project (PPP program, D/08/02059) of the University of Regensburg (Germany) with the Institute of Fish Resources in Varna (Bulgaria). Research was furthermore co-funded by the University of Monastir (Tunisia) and the Institute of Zoology (Chair: Prof. Jürgen Heinze) at the University of Regensburg (Germany), and special thanks is due to the latter for extended financial support. We are also very grateful to four anonymous reviewers for their very helpful and interesting comments and suggestions to improve the quality of the manuscript.


This study was supported by a grant (No: 09S101) from the Research Fund of Boğaziçi University in Istanbul (Turkey) and a DAAD exchange project (PPP program, D/08/02059) of the University of Regensburg (Germany) with the Institute of Fish Resources in Varna (Bulgaria). We acknowledge support of the publication fee by the University of Regensburg.

Availability of data and materials

The datasets supporting the conclusions of this article are included within the article and in the Additional file 1: Table S1.

Author information




TD, RB and CDS conceived and funded the study. TD, EK, SÜK, AK and SU participated in laboratory work. TD carried out the statistical analyses. TD and CDS wrote the manuscript. All authors read, edited and approved the final version of the manuscript.

Corresponding author

Correspondence to Christoph D. Schubart.

Ethics declarations

Ethics approval and consent to participate

None of the sampled populations of Carcinus aestuarii is endangered or protected by any international or national legal framework. In addition, no specific permissions were required for sampling activities.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional file

Additional file 1:

Table S1. Geographic distribution of the 106 haplotypes of C. aestuarii recorded at the 22 sampling sites within the Mediterranean Sea. (DOCX 44 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Deli, T., Kalkan, E., Karhan, S.Ü. et al. Parapatric genetic divergence among deep evolutionary lineages in the Mediterranean green crab, Carcinus aestuarii (Brachyura, Portunoidea, Carcinidae), accounts for a sharp phylogeographic break in the Eastern Mediterranean. BMC Evol Biol 18, 53 (2018).

Download citation


  • Crustacea
  • Population genetics
  • Biogeographic boundaries
  • Evolutionary history
  • Mitochondrial DNA
  • Mediterranean Sea