First endemic freshwater Gammarus from Crete and its evolutionary history—an integrative taxonomy approach

The Mediterranean islands are known as natural laboratories of evolution with a high level of endemic biodiversity. However, most biodiversity assessments have focused mainly on terrestrial and marine fauna, leaving the freshwater animals aside. Crete is one of the largest islands in the Mediterranean Basin, with a long history of isolation from the continental mainland. Gammarid amphipods are often dominant in macrozoobenthic communities in European inland waters. They are widely used in biomonitoring and exotoxicological studies. Herein, we describe Gammarus plaitisi sp. nov., endemic to Cretan streams, based on morphological characters and a set of molecular species delimitation methods using mitochondrial cytochrome oxidase subunit I and 16S rRNA genes as well as nuclear 28S rDNA, ITS1 and EF1-alpha genes. The divergence of the new species is strongly connected with the geological history of the island supporting its continental origin.


INTRODUCTION
Due to its complex geological history and unique combination of geological and climatic factors, the Mediterranean Region is recognized as one of the globally most important hotspots of biodiversity and endemism, and is a model system for studies of biogeography and evolution (Woodward, 2009;Poulakakis et al., 2015). The freshwater fauna of the region is still heavily understudied, yet it is estimated that the Mediterranean is inhabited by ca. 35% of Palearctic species, which mean the region contains more than 6% of the world's freshwater species. At least 43% of the freshwater Mediterranean species are considered to be local endemics (Figueroa et al., 2013). Most of these endemics occupy the Mediterranean islands (Myers et al., 2000;Whittaker & Fernández-Palacios, 2007).
Crete is the fifth largest of the Mediterranean islands and the largest of the Aegean islands. At the beginning of the Miocene, Crete was a part of the mainland composed of the Balkan Peninsula and Asia Minor (23-12 million years ago). Around 12 million years ago, the split of the Balkan Peninsula (including Crete) from Asia Minor began. Afterwards, about 11-8 million years ago, the isolation of Crete from Peloponnesus started, due to the rise of sea levels. Later, between 5.96 and 5.33 million years ago, the dessication of the Proto-Mediterranean Sea during the Messinian Salinity Crysis led to the formation of hypersaline deserts around Crete and other islands, and this is the last known land connection between Crete and the mainland (Poulakakis et al., 2015). During the Pliocene, Crete was divided temporarily into at least four islands due to sea level rise associated with the Zanclean flood (Sondaar & Dermitzakis, 1982). At the end of the Pliocene or in the Early Pleistocene, Crete gained its present configuration.
In this paper, we show evidence that the Cretan population of Gammarus pulex pulex is, in fact, a new species and describe it as Gammarus plaitisi sp. nov., based on morphological, ultrastructural and molecular features. We also reconstruct, based on a multimarker dataset, the phylogeny of this species with respect to other lineages of G. pulex to reveal its biogeographic afiliations and possible origin.

Sample collection, identification and material deposition
The study material was collected from seven out of 53 sampling sites, including springs, streams, rivers and lakes, visited during two sampling campaigns to Crete in 2011 and 2015 (Fig. 1). Multihabitat sampling was done with rectangular kick sample nets (aperture 25 × 25 cm and 0.5 mm mesh size). The samples were sorted at the site and amphipods were immediately fixed in 96% ethanol. Afterwards, the material was evaluated with a Nikon 800 stereomicroscope. Identification to species was done according to the diagnostic morphological characters described by Karaman & Pinkster (1977a), Karaman & Pinkster (1977b), Karaman & Pinkster (1987) and by Pinkster (1993

Scanning electrone microscope analysis
Individuals used for scanning electron microscope (SEM) analysis were critical point dried and sputter-coated with colloidal gold (10 nm). Pictures were taken with a PHENOM PRO X SEM in the Department of Invertebrate Zoology and Hydrobiology of University of Lodz. The photographs of the composition of the pores on antenna 1 and epimeral plate 2 were taken from three same-sized individuals belonging respectively to G. plaitisi sp. nov. and other populations of G. pulex pulex under four different magnifications.

DNA extraction, PCR amplification, sequencing, haplotype diversity and sequence analysis
About 3 mm 3 of the muscle tissue was taken out from each individual, with a sharp-edged forceps and incubated overnight at 55 • C in a 1.5-ml tube containing 200 µl of Queen's lysis buffer with 5 µl of proteinase K (20 mg ml −1 ) (Seutin, White & Boag, 1991). Total DNA was extracted using the standard phenol/chlorophorm method (Hillis, Moritz & Mable, 1996). Air-dried DNA pellets were resuspended in 100 µl of TE buffer, pH 8.00, stored at 4 • C until amplification and finally longterm stored at −20 • C. At first, 57 individuals from 7 sampling sites were barcoded for cox I gene fragment using LCO1490/HCO2198 (Folmer et al., 1994) and LCO1490 Hou, Fu & Li, 2007). The cleaning of the PCR products was done with exonuclease I (20 U mL-1, Fermentas) and alkaline phosphatase FastAP (1 U mL-1, Fermentas) treatment according to the manufacturer's guidelines. Subsequently, the products have been sequenced using the same primers as at the amplification stage. Sequencing of the PCR products was performed using BigDye terminator technology by Macrogen Inc. All resulting sequences were verified and confirmed as Gammarus DNA via BLASTn searches in GenBank (Altschul et al., 1990) and then assembled and aligned in Geneious software (Kearse et al., 2012). The alignment was performed using MAFFT plugin with G-INS-i algorithm in Geneious software (Katoh et al., 2002).
Additional COI sequences of closely related lineages from Greece and Sweden (geographically nearest to type locality of G. pulex pulex), and outgroup Gammarus species were downloaded from NCBI GenBank and added to analyses to test the monophyly of G. cf pulex group (Table 1). The neighbour-joining tree of all COI sequences, using Tamura-Nei model of evolution with 1,000 bootstrap replicates, was created in MEGA7 software (Kumar, Stecher & Tamura, 2016). Afterwards, at least three individuals per each delimited cluster were amplified for one additional mitochondrial and two nuclear markers for phylogeny reconstruction: (1) mitochondrial 16S rRNA using 16STf and 16SBr markers (Palumbi et al., 1991;MacDonald IIII, Yampolsky & Duffy, 2005) under the following PCR conditions: initial denaturation at 94 • C for 150 s; 36 cycles of denaturation at 94 • C for 40 s, annealing at 54 • C for 40 s, extension at 65 • C for 80 s; and a final extension at 65 • C for 8 min (Weiss et al., 2014); (2) the nuclear 28S rRNA gene amplified with 28F and 28R primers (Hou, Fu & Li, 2007) under following conditions: initial denaturation at 94 • C for 3 min, 35 cycles of denaturation at 94 • C for 20 s, annealing at 55 • C for 45 s, and elongation at 65 • C for 60 s, followed by a final extension for 2 min at 65 • C and 5 min extension at 72 • C; (3) the nuclear ITS1 gene with ITS1F and ITS1R primers (Chu, Li & Ho, 2001)

MOTU delimitation-cryptic diversity
The Molecular Operational Taxonomic Units (MOTUs) were delimited, based on the COI marker, with five methods and two different approaches (as done before by Grabowski, Wysocka & Mamos, 2017): the distance-based approaches, namely Barcode Index Number (BIN) System (Ratnasingham & Hebert, 2013) and barcode gap discovery with the ABGD software (Puillandre et al., 2012) and the tree-based approaches, using two GMYC modelbased methods (Pons et al., 2006) according to Monaghan et al. (2009) and the bPTP procedure described by Zhang et al. (2013).
The BIN method is a distance-based approach, embedded in the Barcode of Life Data systems (BOLD; Ratnasingham & Hebert, 2007). The sequences already deposited in BOLD database are confronted with the newly submitted ones. Afterwards, according to their molecular divergence, the sequences are clustered using algorithms that identify discontinuities between the clusters. A unique and specific Barcode Index Number (BIN) is assigned to each cluster. If the submitted sequences do not group together with already known BINs, a new number is created. Each BIN is registered in the BOLD database.
The Automated Barcode Gap Discovery (ABGD) method uses pairwise distance measures. ABGD clusters the sequences into MOTUs (Molecular Operational Taxonomic Units), in the way that the genetic distance between two sequences belonging to two separate groups will always be greater than an indicated threshold (i.e., barcode gap). In our study, the primary partitions were used as a principal for cluster delimitation, as they tend to remain stable on a wider range of prior values, minimising the oversplitting of the number of groups and are usually the closest to the number of taxa described by taxonomists (Puillandre et al., 2012). The default value of 0.001 was applied as the minimum intraspecific distance. As the maximum intraspecific distance we investigated a set of values up to 0.1, which has been proposed as suggested maximum distance value in amphipods distinguishing two separate species (Costa et al., 2007). The standard Kimura two-parameter (K2P) model correction was used (Hebert, Cywinska & Ball, 2003).
The bPTP approach for species delimitation is a tree based method, utilising nonultrametric phylogenies. The number of substitutions in incorporated into the model of speciation and the bPTP assumes that the probability that a substitution leads to a speciation event follows a Poisson distribution, as the lengths of the branches of the input tree are generated independently according to either to speciation or coalescence, which are two classes of the Poisson processes. In bPTP, the Bayesian support values are added for each delimited cluster (Zhang et al., 2013). As an input tree, the phylogeny was generated using Bayesian inference in Geneious software package using MrBayes plugin (Kearse et al., 2012) with MCMC chain 1 million iterations long, sampled every 2,000 iterations. The TN93 + I + G was chosen as the substitution model, as bestfit based on bModel test (Bouckaert & Drummond, 2017). The consensus tree was constructed after removal of 25% burn-in phase. The analysis itself was done using the bPTP web server (http://species.h-its.org/) with 500,000 iterations of MCMC and 10% burn-in.
The GMYC method identifies the transition from intraspecific branching patterns (coalescent) to typical interspecific branching patterns (Yule processes) on an ultrametric, phylogenetic tree, using the maximum likelihood approach. The estimation of the boundary between coalescent and Yule branching processes can be done using two different GMYC approaches, one using the single threshold and the second one based on multiple threshold model. We have reconstructed an ultrametric tree, which is required for GMYC analyses, in BEAST software, using 20 million iterations long MCMC chain, with TN93+ I + G as the best-fit substitution model. The consensus tree was analysed in the GMYC web server (available at: http://species.h-its.org/gmyc/) using both the single and multiple threshold models.

Time calibration and phylogeny reconstruction
The time-callibrated phylogeny was reconstructed based on data from sequences of COI (586 bp), 16S rRNA (299 bp), 28S rRNA (781 bp), ITS1 (548 bp) and EF1-alpha (602 bp) in BEAST2 software package (Bouckaert et al., 2014) with the use of five MCMC chains of 50,000,000 runs with following models of substitution: TN93 + I + G (for COI), HKY + I + G (for 16S), TN93 + I + G (for 28S), HKY + I + G (for ITS1) and TN93 + I + G (for EF1-alpha). The models for each marker were selected according to bModel test (Bouckaert & Drummond, 2017). The relaxed log-normal clock model was used and based on the selected rate of 0.0115 substitutions (SD 0.0026) per million years for COI according to already established rate (Brower, 1994), which was cross-validated against two other rates (0.0113, 0.0127) established recently for other freshwater members of Gammarus, in the G. roeselii species complex (Grabowski et al., 2017). All other clock rates were set on estimate. For 16S rRNA and EF1-alpha also relaxed log-normal clock was used, whereas for 28S rRNA and ITS1 the strict clock was used. All the models were tested beforehand in MEGA software, using an implemented test for molecular clock model based on Maximum Likelihood phylogeny (Kumar, Stecher & Tamura, 2016). The resulting trees were checked for ESS values in Tracer and two trees with the best ESS values were combined in LogCombiner and annotated in TreeAnnotator. The final output tree was edited in FigTree software (http://tree.bio.ed.ac.uk/software/figtree/). Order: Amphipoda Latreille, 1818Family: Gammaridae Leach, 1814Genus: Gammarus Fabricius, 1775Pinkster, 1970: 179, Karaman & Pinkster, 1977a: 3, Barnard & Barnard, 1983  Island, Greece) Diagnosis: Large species, making a robust impression. Similiar to G. pulex pulex by the characteristic antenna 2 with swollen flagellum, bearing a flag-like dense brush of setae and similar armature of pereiopods. It may be distinguished from G. pulex pulex by the lack of spines on the dorsal surface of the first segment of urosome, the shape of the posterodistal margin of the second and third epimeral plate and by the size and the arrangement of the pores on the cuticle surface. It is also clearly distinguishable from G. pulex pulex on the molecular level, with respect to the COI nucleotide sequence.  (Fig. 2C): First segment unarmed. Second segment with ventral setae: in the proximal part 2-3 setae much shorter than the diameter of the segment, in the distal part 10-13 setae as long as or up to 2.5× longer than the diameter of the segment. Third segment armed with 2 groups of long A-setae, a regular comb of 25-30 D-setae and 5-6 long E-setae. Maxillipeds (Fig. 2D): The maxillipeds with the inner plate armed distally with strong spine-teeth; the outer plate with spine-teeth and long plumose setae; the palp is well developed. Gnathopod I (Fig. 2E): Palm oblique, setose, with one strong medial palmar spine, strong angle spine accompanied by several small spines intermixed with longer setae along the posterior palmar margin with addition of small spines and short setae on the lateral surface. Gnathopod II (Fig. 2F): Propodus trapezoid, widening distally. Palm concave, setose, with one medial palmar spine and three angle spines. Many groups of setae, variable in length, are visible both on the inner and outer as well as the lateral surface of the propodus Pereopod III (Fig. 3A): Anterior and distal margin of coxal plate slightly convex, posterior margin straight. Distal corners rounded. The last three segments of third pereiopod bear groups of long, often curved setae along the posterior margin, usually two to three times longer than the diameter of segments. The anterior margin of merus armed with 1 spine. Dactylus short, robust with one seta at joint of unguis. Pereopod IV (Fig. 3B): Coxal plate dilated distally. Distal corners rounded. The last three segments of fourth pereiopod bear groups of long, often curved setae along the posterior margin, usually two to three times longer than the diameter of segments. The anterior margin of merus armed with 1 spine. Dactylus short, robust with one seta at joint of unguis. Pereopod V (Fig. 3C): Basis with a subrectangular shape, posterior margin slightly concave, posterodistal lobe well developed, posterior margin with 10-12 very short setae, anterior margin with 4-5 spiniform setae. Ischium naked. Merus, carpus and propodus with robust spines on both margins, occasionally intermixed with relatively short setae. Dactylus short, robust usually with one seta at joint of unguis. Pereopod VI (Fig. 3D): Similar to PV, but slightly longer and wider, posterior margin convex, posterodistal lobe less prominent and basis more more elongated with a single, little spine on posterointerior corner. Ischium to propodus armed with robust spines and very few short setae. Dactylus short, robust with one seta at joint of unguis. Pereopod VII (Fig. 3E): Basis wider than in PVI with a single, little spine only at posteroinferior corner and even more elongated. Further articles armed same as in preceding pereopods. Uropod III (Fig. 3F): The inner ramus attains about 2/3 of the length of the outer ramus. Most of setae along the inner and outer margin of endo-and exopodite plumose. Telson (Fig. 3G): Deeply cleft, rather setose. Each lobe with two apical strong spines intermixed with few short and long setae, several short subapical setae present. Epimeral plates (Fig. 3H): First epimeral plate with one spine at the laterodistal margin. Second epimeral plate with one spine at the laterodistal surface, posterodistal margin rounded. Third epimeral plate with three spines at the laterodistal surface, posterodistal margin rounded with the posterodistal corner slightly pointed. Urosome (Fig. 4A): very flat without any elevation. First urosomite lacking any spines on dorsomedial or dorsolateral surface and armed only with a few groups of setae.

Systematics
Second urosomite with dorsomedial and dorsolateral groups of robust spines (2-2-2). Third urosomite only with two groups of dorsolateral spines on each side (3-0-3), and a dorsmedial group of 2-4 setae. Ultrastructure (Figs. 5 and 6) The pores are larger and more distinctly marked in comparison to G. pulex pulex. This pattern holds true for both A1 and E2, however on A1 the difference is more pronounced. On A1 pores form the regular rows for both G. plaitisi sp.nov. and G. pulex pulex, whereas on E2 the rows of pores are much more regular in G. plaitisi sp.nov. compared to those in G. pulex pulex. The distances between rows of pores are always about 1.5 times wider than in G. pulex pulex. Female: Smaller than male. The setation of the peduncle segments of the first and second antennae is longer than in the male. The characteristic brush of second antenna flagellum is absent. The propodi of the gnathopods smaller than in males and the setation of P3 and P4 is less abundant and shorter. Variability: Morphology of G. plaiti is stable with respect to features such as presence of calceoli in males, presence of brush in peduncle of A2, flatness and armature of urosomites.
Larger individuals tend to have higher number of flagellum segments in antenna I and II, as well as more and longer setae on all appendages. The density of the setation and spinulation is also rather variable depending on age of the individual. Such variability is typical for most species of this genus (Karaman & Pinkster, 1977a;Karaman & Pinkster, 1977b;Karaman & Pinkster, 1987).

Haplotype diversity and phylogeny reconstruction
We identified three haplotypes of G. plaitisi sp.nov. in the dataset composed of the forty three COI sequences, with one haplotype being represented only by one specimen. The most common haplotype, H2, was present in the majority of sites, except for locus typicus of the species (Fig. 7). The overall haplotype diversity was quite high (Hd = 0,375 ± 0,076), whereas nucleotide diversity (Pi = 0,00126 ± 0,00075) was low. Generally, the differentiation was very low as the most common haplotype differed from the two remaining ones by a maximum of two mutation steps with intraspecific distance not exceeding the value of 0.005. All MOTU delimitation methods supported distinctness of G. plaitisi, which always formed a single MOTU and was separated from its closest relative by the mean K2P distance of 0.12 (Table S1). It also formed a unique BIN in the BOLD database (BOLD: ADG8205). All the applied MOTU delimitation methods provided constant results with six MOTUs delimited for the G. pulex morphospecies. Only the ABGD method indicated one MOTU less within the Peloponnese group. Both the used GMYC approaches produced the same outcome with the same LR test values. Results of MOTU delimitation methods support high cryptic diversity within Gammarus pulex morphospecies from Greece, as no morphological differences amongst the representatives of respective MOTUs have been found. The topology of the neighbour-joining tree confirms that G. plaitisi sp. nov. is nested within the clade of lineages belonging to the G. pulex morphospecies (Fig. 8). This suggests that G. pulex is, in reality, a paraphyletic group of cryptic and pseudocryptic species.
Multimarker time-calibrated phylogeny indicated that divergence of the whole G. pulex lineages from Peloponnese happened around 15 million years ago, whereas divergence of  G. plaitisi sp.nov. from its continental relatives took place around 9.2 million years ago Moreover, divergence within the continental groups of G. pulex lineages spanned the last 5 million years (Fig. 9). All three rates used for time calibratead reconstruction of Bayesian phylogeny gave congruent results (Table 2).

DISCUSSION
We provided evidence for the existence of new freshwater Gammarus species from Crete, making this the third known freshwater endemic gammarid to Crete. The endemic freshwater species of Gammaridae before this work were, Echinogammarus platvoeti and E. kretensis (Pinkster, 1993), making G. plaitisi sp. nov. the first endemic of the genus Gammarus. The integrative taxonomy approach confirmed the distinctness of the species not only on a morphological basis, but also on a molecular level. This study also stressed the importance of using SEM photography, which may provide additional diagnostic features that are impossible to detect on usually used optical devices (Platvoet et al., 2008).
Despite the presence of G. plaitisi sp. nov. in seven, mostly isolated sites located both in the eastern and western part of Crete, its haplotype diversity is surprisingly low, with only two mutation steps separating the three known haplotypes (Table 3). This pattern suggests a strong founder effect and recent dispersal, probably in the late Pleistocene, as suggested by the time-calibrated phylogeny, possibly due to rearrangement of the local hydrological networks at the end of the last ice age. This is a rather unusual finding considering the fact that Pleistocene glaciations, which strongly affected the river systems, promoted the diversification of various taxa in the Mediterranean (Previšić et al., 2009;Goncalves et al., 2015), including the freshwater gammarids (Grabowski et al., 2017a). However, such a founder effect scenario has also been found in other freshwater members of the genus Gammarus, such as Gammarus minus which inhabits both surface and groundwaters of North America. Gooch & Glazier (1986) confirmed postglacial dispersal of this species from refugia, which resulted in strong decrease in their allele diversity. This scenario is the most plausible one also for G. plaitisi sp. nov., which may have colonised the current distribution area from a single refugium. The distribution of haplotypes (Fig. 7) suggests that the individuals originate from a founding population from the western part of Crete, where all of the known haplotypes are present. Yet another question concerns the way of dispersal between isolated freshwater systems, separated by more than 100 km. One must consider passive dispersal i.e., by birds (Rachalewski et al., 2013), however, groundwater connections cannot be excluded (Harris, Roosa & Norment, 2002). On the other side, there may be still some localities, particularly in the mountains, where the species is present or could have been present in the early Holocene but died out due to climatic changes. We still do not have enough data to reveal the dispersal history of this species.
Our results suggest that G. plaitisi sp. nov. has diverged from the continental lineages of G. pulex around 9 million years ago (Fig. 9). This result has been strongly supported by cross-validation with other substitution rates proposed for freshwater gammarids in earlier studies (Grabowski et al., 2017a). The timescale seems to be convergent with the estimated date of the first isolation of Crete from Peloponnese (Poulakakis et al., 2015). Since that time Crete could be colonized only by overseas dispersal. This finding suggests the continental origin of the newly described species. The molecular data suggest rather the possibility of its dispersal to Crete before first isolation of this island than migration during the temporal land connection during the Messinian Sality Crysis and after its final isolation at around 5 million years ago.

CONCLUSIONS
G. plaitisi sp. nov. is the first endemic insular freshwater Gammarus in the Mediterranean. However, given the scarcity of the sampling in the fresh waters of the Mediterranean islands, there is a high chance there are more representatives of the genus in the Aegean Basin and other Mediterranean islands. The description of this new species using the integrative taxonomy approach not only broadens the knowledge about freshwater diversity of Crete, but also provides a link between the geological history of this island with the evolution of the local freshwater species. The results provide yet another piece of the puzzle in explaining the evolution of the family Gammaridae.