The late Pleistocene origin of the Italian and Maltese populations of Potamon fluviatile (Malacostraca: Decapoda): insights from an expanded sampling of molecular data

Abstract

Evidence available for most inland water and terrestrial organisms highlights the significant role played by southern Italy, Sicily and the Maltese islands as refuges during Pleistocene climatic fluctuations. However, to date, the hypothesis that these areas may have acted as Pleistocene refugia for the freshwater crab Potamon fluviatile has not been explicitly tested, and a recent origin of local P. fluviatile populations was proposed on the basis of a small set of analysed molecular data. We have thus expanded the currently available data set on the population genetic structure of P. fluviatile through dedicated samplings in Sicily (Italy, 18 specimens), the Maltese Islands (Malta, 15 specimens) and the island of Corfu (Greece, seven specimens), with the explicit aim of testing the role of southern Italy, Sicily and the Maltese archipelago as refugia where early Pleistocene or pre-Quaternary populations of the freshwater crab could have persisted to date. Based on the analysis of both novel and published cytochrome c oxidase subunit I (COI) sequence data, we gathered evidence supporting a recent origin of the southern Italian, Sicilian and Maltese populations of Potamon fluviatile from the Balkans, thus rejecting the hypothesised Pleistocene persistence in the area of pre-existing local populations of the same species.

Keywords:

• Freshwater crab
• mtDNA Cytochrome c oxidase subunit I (COI)
• range expansion
• refuge areas

Introduction

The freshwater crab family Potamidae is classified into two subfamilies: Potaminae (distributed across Europe, North Africa, Yemen, the Middle East and India) and Potamiscinae (Eastern and South-eastearn Asia) (Yeo & Ng 2003). A revision of the genus Potamon Savigny, 1816 in the West-Palearctic region lists 16 species, divided into four mostly allopatric or parapatric subgenera that were distinguished through morphological characters, mostly related to the gonopods (Brandis et al. 2000). In particular, the subgenus Euthelphusa Pretzmann, 1962, occurring in the western Mediterranean area and on the Balkan peninsula, includes the species Potamon fluviatile (Herbst, 1785), P. pelops Jesse et al., 2010 and P. algeriense Bott, 1967 (Brandis et al. 2000). P. fluviatile is distributed across Italy, Malta, Croatia, southern Dalmatia, Albania and Greece (Brandis et al. 2000), usually segregated from other autochthonous macro-decapods (Barbaresi et al. 2007; Mazza et al. 2017). Adults may reach a carapace length of about 50 mm over a 10-year-long lifespan (Gherardi et al. 1987; Micheli et al. 1990; Scalici et al. 2013); they are omnivorous, feeding primarily on invertebrates, amphibians and macrophytes (Cumberlidge 2008). Since 2008, P. fluviatile has been assessed as “Near threatened” on the International Union for Conservation of Nature (IUCN) Red List of Threatened Species ( http://www.iucnredlist.org/details/134293/0 ), because the populations of the species historically declined as a result of pollution, over-harvesting and mismanagement of habitats (Matthew & Reynolds 1995; Gherardi & Holdich 1999).

In the last decade, some studies were conducted on P. fluviatile in order to investigate its phylogenetic relationships with the other species within the subgenus Euthelphusa, and to elucidate its phylogeography (see Jesse et al. 2010 and references therein). Based on the available molecular data, Jesse et al. (2009, 2010, 2011) proposed a non-human-mediated late Pleistocene origin from the Balkans of the Italian and Maltese populations of the species. However, it is well known that, during Pleistocene glaciations, southern Italy, Sicily and the Maltese archipelago acted as refugia for different species which persisted in situ (Hewitt 2000, 2004), with a resulting genetic differentiation of these “southern” strains, for both aquatic (e.g. Canestrelli et al. 2007; Canestrelli & Nascetti 2008; Kindler et al. 2013; Vamberger et al. 2015) and terrestrial (e.g. Michaux et al. 2003; Bezerra et al. 2015) taxa.

The purpose of our study was thus to explicitly test the possible role played by southern Italy and Malta as refugia for P. fluviatile during the Pleistocene, through the implementation of an expanded sampling protocol of the genetic diversity of different populations of the species in Sicily and in the Maltese archipelago. If this geographic area did indeed act as a Pleistocene refugium for the species, we expect to find the presence of endemic haplotypes, possibly mixed with widespread ones shared with those of the Balkan populations of the species, as in the case of incomplete lineage sorting or as in the case of a secondary contact between the populations inhabiting the “southern Italian” and “Balkan” refugia during the late Pleistocene.

Material and methods

In total, 40 novel specimens (1–7 per site; see Table I) were analysed, taking one leg from each individual, and releasing them back in the same sampling site where they were collected. The loss of one leg does not constitute harm to the sampled individuals, since they are able to routinely regenerate any lost appendages after the next molt (cf. Das 2015 and references therein). Sampled legs were fixed in situ in 90% ethanol. In addition, 21 haplotype sequences, corresponding to 104 individuals collected throughout the known distribution range of the species, were downloaded from GenBank (accession numbers: EU908226–EU908246) and analysed along with our sequences. Further information on the sampling sites and number of sequenced individuals is given in Table I. The map of the sampling sites (see Figure 1) was plotted using the QGIS software v. 2.18.2 ( http://www.qgis.org ).

Table I. Location of the sampling sites.

Location	North	East	n.	Haplotype(s) found
Corfu (Greece)
Valanio	39.7325°	19.7443°	7	NH2
Gozo (Malta)
Xlendi	36.0332°	14.2273°	4	NH2; NH11
Malta (Malta)
Bahrija	35.8947°	14.3385°	6	NH1; NH2
Ghajn Zejtuna	35.9657°	14.3672°	3	NH1; NH2
San Martin	35.9329°	14.3762°	2	NH2
Sicily (Italy)
Fiume Cassibile	36.9882°	15.0268°	1	NH2
Fiume Tellaro	36.8837°	14.9504°	9	NH2
Torrente Frattina	37.8613°	13.3030°	4	NH2
Imera Settentrionale	37.8583°	13.8970°	1	NH10
Fiume Irminio	36.7886°	14.6020°	1	NH2
Torrente Calcinara	37.1402°	15.0291°	2	NH2

Coordinates are expressed as decimal degrees (Map Datum: WGS84); n: number of sampled individuals in each site.

Figure 1. Map of the sampling sites of Potamon fluviatile.

Genomic DNA was extracted from the muscle tissue of each leg using the Genomic DNA Extraction Kit (RBC BioScience) following the manufacturer’s protocol. The selective amplification of a 502-bp-long fragment from the gene encoding for the cytochrome oxidase subunit I (COI) was carried out by the polymerase chain reaction (PCR). The PCR mix consisted of 18.05 µL double-distilled water, 2.5 µL Buffer 10X including 15 mM MgCl₂ solution, 0.25 µL dNTPs (10 mM of each), 0.9 µL of each primer (10 µM), 0.4 µL Taq Polymerase 5 u/µL and 2 µL of DNA template, for a total volume of 25 µL. The thermal cycle (35 cycles; denaturing 95°C for 45 s/annealing 46°C for 45 s/extension 72°C for 50 s) was conducted with the primers COL6b (5ʹ-ACA AAT CAT AAA GAT ATY GG-3ʹ) and COH6 (5ʹ-TAD ACT TCD GGR TGD CCA AAR AAY CA-3ʹ), described by Schubart and Huber (2006). After PCR, 5 µL of each PCR product were separated by electrophoresis on a 2% agarose gel at 90 V for 20 minutes and visualised with a UV Transilluminator. When PCR products showed a clear and single band of the correct expected length, they were purified using the Exo-SAP-IT® kit (Affymetrix USB) and sequenced by Macrogen Inc. (Seoul, South Korea) with ABI 3130xL (Applied Biosystems) sequencer. The forward primers were used for direct sequencing of the PCR product. The quality of the obtained chromatograms was checked through the measurement of their Phred scores (Richterich 1998). Only sequences with continuous reads of high-quality bases Quality Value scores ((QV) > 20) were used. Sequences were analysed and manually proofread with the DNA sequencing software Chromas v. 2.6.2 (Technelysium, Pty. Ltd. 1998, Queensland, Australia) and aligned with ClustalX v. 2.1 (Larkin et al. 2007). MEGA v. 7.0 (Kumar et al. 2016) was used to translate the COI sequences into amino acids in order to check for the possible presence of frameshifts or stop codons, which would indicate the presence of sequencing errors or pseudogenes.

Upon aligning our novel sequences with those downloaded from GenBank, we had to trim out 35 bp from the sequences downloaded from GenBank, which were originally 537 bp long. This trimming led to the merging of some haplotypes (see Table S1 for correspondence of our haplotype codes with those described by Jesse et al. 2009).

A haplotype network for the identified sequences was constructed through the software Popart (v. 1.7, http://popart.otago.ac.nz ), using the median-joining method (Bandelt et al. 1999) (Figure 2).

Figure 2. Median-joining haplotype network based on a 502-bp-long fragment of the mitochondrial DNA gene Cytochrome oxidase subunit 1 (COI). Dashes indicate substitution steps. Each circle represents a haplotype and its size is proportional to its frequency.

We partitioned a priori the sampling sites into three geographical groups: (i) Balkans; (ii) Central Italy; and (iii) southern peninsular Italy, Sicily, and the Maltese Islands (hereafter referred to as “ISM”) in order to detect population differentiation through the analysis of their molecular variance (AMOVA; Excoffier et al. 1992). Tajima’s D test and Fu’s FS test (Tajima 1989; Fu 1997) were performed in order to evaluate the population demography of this species. A mismatch distribution analysis was performed to estimate the parameters of demographic expansion. AMOVA, neutrality tests and mismatch analysis were performed with the software Arlequin (v. 3.5.2.2, Excoffier & Lischer 2010). In addition, haplotype diversity and nucleotide diversity were computed with the software DnaSP v. 5 (Librado & Rozas 2009). In order to estimate the significance of the pair-wise differences in haplotype diversity between the three geographical groups, we used a series of functions, developed by Thomas et al. (2002) and available at http://www.ucl.ac.uk/mace-lab/resources/software (“TEST_h_DIFF”), through the statistical software “R” (v. 3.2.1, https://www.R-project.org ).

Past population demography was inferred using a Bayesian skyline plot model (Drummond et al. 2005), as implemented in BEAST v. 1.8.0 (Drummond et al. 2012) and visualised in TRACER v. 1.6 (Rambaut et al. 2014). The coalescent approach uses standard Monte Carlo Markov chain (MCMC) sampling procedures to estimate a posterior distribution of gene genealogies and population parameters under a TN93 + I substitution model (Tamura & Nei 1993). The choice of the correct evolutionary model was conducted using jModelTest (v. 2.0, Guindon & Gascuel 2003; Darriba et al. 2012) according to the Akaike information criterion (AIC; Akaike 1974). A log-normal model that relaxes the molecular clock hypothesis was used (Drummond et al. 2006). The hyper-parameter “m” (number of grouped intervals) was set to 10. Three independent MCMC analyses of 10⁷ steps were performed, sampling every 10,000th generation with the burn-in set at 10⁶ generations. To calibrate the molecular clock, we used the same split as described by Jesse et al. (2009), i.e. the one between Potamon fluviatile and P. algeriense, which most likely took place at the end of the Messinian salinity crisis (Brandis et al. 2000).

Results

The obtained sequences aligned unambiguously and translated without stop codons, thus suggesting that no pseudogenes were erroneously amplified. A total of 11 haplotypes, including both the haplotypes observed by Jesse et al. (2009) and two new haplotypes recorded within the framework of the current study (“NH11” and “NH10”; GenBank accession numbers KY449461 and KY449462; Table I), were obtained. The three geographical groups showed the following haplotype diversity values (Hd): 0.681 (Balkans), 0.040 (Central Italy) and 0.527 (ISM). The following nucleotide diversity (Pi) values were obtained: 0.0042 for the group “Balkans”, 0.0008 for the group “Central Italy”, and 0.0298 for “ISM”. In terms of pair-wise haplotype diversity differences between the three geographical groups, we obtained significant p values (< 0.0001), whereas the comparison between the Balkan and the ISM cohorts was insignificant (p value = 0.057). The haplotype network obtained shows that there are two common and widespread haplotypes shared by all three groups, plus some rarer haplotypes shared by a single or by a few individuals which are restricted to single areas (Figure 2). The AMOVA test (Table II) shows that the percentage of variation between the three geographical groups is 41.90% while the variation within the groups is 58.10%, with an Fixation index (F_ST) value of 0.419 and a p value < 0.0001.

Table II. AMOVA and F_ST results based on 502-bp-long fragments of the mtDNA COI gene.

Source of variation	df	Sum of squares	Variance components	Percentage of variation
Among groups	2	51.335	0.537 Va	41.90
Within groups	141	105.013	0.744 Vb	58.10
Total	143	156.368	1.281

Fixation Index F_ST: 0.419; p value: 0.000.

“df” are the degrees of freedom, “Va” is the variance among groups, “Vb” is the variance among populations within groups.

The neutrality tests are not significant for the “Balkans” and “ISM” groups, whilst the Fu’s FS statistics p value is only significant for the group “Central Italy”. Specifically, Tajima’s D has a value of −1.185 for “Balkans”, −1.100 for “Central Italy” and 0.367 for “ISM”. Conversely, Fu’s FS values of −2.736 for the “Balkans” and of −1.652 for the “Central Italy” group were recorded, with a p value and a corresponding Fu’s FS value of 0.037 and 1.246, respectively, for the “ISM” group (Table III).

Table III. Neutrality test results based on 502-bp-long fragments of the mtDNA COI gene.

	Statistics	Balkans	Central Italy	Southern Italy – Sicily – Malta	Mean	SD
Tajima’s D Test
	Sample size	32	50	62	48.000	15.099
	S	23	1	6	10.000	11.532
	Pi	3.772	0.040	1.493	1.768	1.881
	Tajima’s D	−1.185	−1.102	0.413	−0.624	0.900
	Tajima’s D p value	0.103	0.121	0.716	0.317	0.345
Fu’s FS test
	Real no. of alleles	13	3	5	7.000	5.291
	Orig. no. of alleles	32	50	62	48.00	15.099
	Theta_pi	3.772	0.040	1.493	1.768	1.881
	Exp. no. of alleles	8.954	1.176	6.143	5.424	3.938
	FS	−2.736	−1.636	1.308	−1.021	2.091
	FS p value	0.121	0.037	0.771	0.309	0.401

SD: Standard deviation; S: Segregating sites; Pi: Nucleotide diversity; Theta_pi: Average number of pairwise differences.

The mismatch analysis carried out on the mtDNA COI sequences showed a skewed unimodal curve with the observed number of pair-wise differences, which differs non-significantly from the expected values (Sum of square deviations (SSD) p value = 0.020; raggedness index (RI) = 0.317; raggedness p value = 0.010; Figure 3).

Figure 3. Result of the mismatch distribution analysis based on a 502-bp-long fragment of the mitochondrial DNA gene Cytochrome oxidase subunit 1 (COI). The dotted line represents the observed frequency distribution of pairwise differences; the curve represents the expected one.

The estimated change in effective population size over time with 95% confidence intervals as computed with the Bayesian skyline plot analysis (BSP) is shown in Figure 4. As indicated by the BSP analysis, the demographic expansion began approximately 12,000 years before present (BP), with a 95% confidence interval ranging from 6000 to 24,000 years BP.

Figure 4. Bayesian skyline plots: evolution of population sizes through time based on Bayesian coalescence. Relevant historic events are indicated. LGM: last glacial maximum.

Discussion

No apparent associations between observed haplotypes and their geographic locations were found, i.e. individual haplotypes could not be ascribed to exclusive geographical locations. However, the haplotype diversity values observed for the three geographical groups indicate the occurrence of a highly significant differentiation for the pairs “Balkans vs. Central Italy” and “ISM vs. Central Italy”; conversely, no significant difference was scored for the pair “Balkans vs. ISM”.

The performed analyses show that the studied Potamon fluviatile populations are not geographically structured, and that the highest haplotype diversity is present in the Balkan area, despite the relatively low number of samples studied from this region. Our results are thus in agreement with the hypothesis formulated by Jesse et al. (2009) that, during the Pleistocene, southern Italy and the Maltese Islands did not act as refugia where early Pleistocene or pre-Quaternary P. fluviatile local populations survived the glaciation periods, but rather they proved to be a “land of conquest” for the Balkan P. fluviatile populations, which could colonise these new lands made accessible by virtue of the marine regressions linked with Pleistocene glaciations. Such a hypothesised dispersal pattern finds a counterpart in other freshwater crustaceans, such as in the anostracan Chirocephalus kerkyrensis Pesta, 1921, whose demographic expansions consistently coincided with the interglacial warm phases of the Pleistocene, while its spatial expansion (and colonisation of Italy from a Balkan origin) was restricted to the Pleistocene cold periods (Ketmaier et al. 2012). Later on, at the end of the Pleistocene, when temperatures rose gradually, Potamon fluviatile might have actively colonised more central and northern regions of Italy from its southern enclaves situated along the same peninsula. This secondary northward migration of the species was denoted by relatively low haplotype diversity values in analysed central Italy populations, despite the relatively high number of sampled individuals from this region, which in turn might be symptomatic of a founder effect.

The existence of a recent genetic exchange between the different populations within the studied areas (or of recent colonisation events via extensive dispersal with no founder effect) is suggested by the AMOVA test results. In the light of the F_ST values obtained and the results from the mismatch analysis, the species is inferred to be in a process of ongoing demographic expansion.

In the Bayesian skyline plot analysis, the upper 95% confidence interval suggests that this inferred demographic expansion could have started even before the last glacial maximum (LGM), when the freshwater crab is assumed to have reached southern peninsular Italy, and subsequently Sicily and the Maltese archipelago, through a trans-Adriatic route (Jesse et al. 2009). Moreover, considering that there is evidence of maritime activities in the Western Mediterranean dating back to 6000 years BP (Broodbank 2006 in Jesse et al. 2009), from the lower 95% confidence interval we cannot exclude the possibility of an anthropogenic introduction of Potamon fluviatile in peninsular Italy, which would thus be an example of “cryptic invasion” (e.g. Marrone et al. 2011). However, areas recently colonised through a human-mediated introduction are usually characterised by substantially lower haplotype richness values when compared, for example, to the haplotype richness values observed in the current study for the group “ISM”. Consequently, and in accordance with Jesse et al. (2009), we consider the anthropogenic introduction of the species in Italy and Malta as rather unlikely.

In the light of the likely autochthonous status of the species throughout its current distribution range, a renewed conservation effort aimed at maintaining the populations of this freshwater decapod should be made, especially since they are currently facing multiple pressures (Gherardi & Holdich 1999; Füreder et al. 2002; Svoboda et al. 2014).

Acknowledgements

This study was funded by the Research Committee of the University of Malta, to which the authors are greatly indebted. The Environment and Resources Authority (ERA) of Malta is acknowledged for having provided the sampling permit to AD (File ref.: NP 27/15) authorising the sampling of Potamon fluviatile individuals within four Maltese sites. Francesco Sacco (Palermo, Italy) is acknowledged for the support he provided.

Supplemental data

Supplemental data for this article can be accessed here

Additional information

Funding

This work was supported by the Research Committee of the University of Malta.