Neogene paleogeography provides context for understanding the origin and spatial distribution of cryptic diversity in a widespread Balkan freshwater amphipod

Sampling and taxonomic identification

A total of 348 sites located in springs, streams, rivers and lakes all over the Balkan Peninsula, including the Peloponnese (Fig. 1) were surveyed in the years 2006–2008, using semi-quantitative effort, based on a kick sampling method with a benthic hand-net (e.g., Jażdżewski, Konopacka & Grabowski, 2004). All the collected material was sorted on-site and gammarids were immediately fixed in 96% ethanol. In the laboratory, animals were identified to species level using morphological characters described in available keys (e.g., Karaman & Pinkster, 1977a; Karaman & Pinkster, 1977b; Karaman & Pinkster, 1987; Pinkster, 1993). Gammarus roeselii was found at 26 sites (Table 1). All the material used in this study has been stored in the permanent collection of the Department of Invertebrate Zoology and Hydrobiology, University of Lodz, Poland.

Molecular procedures (DNA extraction, PCR amplification and sequencing)

Genomic DNA was extracted from pieces of muscle tissue obtained from 1 to 12 individuals from each site, making up a total of 177 individuals (Table 1). All individuals were processed, via PCR, to amplify a ca. 650 bp fragment of the mtDNA cytochrome C oxidase subunit I (COI), using the primer pairs: LCO1490 and HCO2198 (Folmer et al., 1994); UCOIF and UCOIR (Costa et al., 2009); and a newly designed primer pair COIGrF (5′-GCTAGH GCCGTAGGYACATC) and COIGrR2 (5′-RAATARGTGYTGGTACAGAATAGG). Details about which primers were specifically used for each individual are provided in Table S1. After defining molecular operational taxonomic units (MOTU) based on the COI marker (see ‘Results’), a subset of 50 individuals representing all the defined MOTUs was also amplified to provide a ca. 900 bp long fragment of the 28S nrDNA using 28F and 28R primers (Hou, Fu & Li, 2007). An aliquot of each PCR product was used to check the specificity and yield by agarose gel electrophoresis. PCR products were EXO-FastAP (Thermo Scientific) purified prior to sequencing by Macrogen Inc., Korea, using BigDye terminator technology. All the details of the molecular procedures were as described in Mamos et al. (2016).

Sequence data authentication, editing and alignment

All sequences were positively verified as Gammarus DNA via BLASTn searches against GenBank (Altschul et al., 1990). Geneious 6.1.4 (Biomatters Inc, 2013) was used for sequence editing, trimmed to either 530 nucleotides (COI) or 563 nucleotides (28S) and aligned using Geneious alignment algorithm. Haplotypes were identified for each marker with the DnaSP software (Librado & Rozas, 2009). All single marker haplotypic sequences were deposited in GenBank (Accession numbers in Table 1). For COI, Kimura 2-parameter (K2p) distance was calculated (e.g., between all haplotypes, or within and among MOTUs) using MEGA 6.0 (Tamura et al., 2013). Standard error estimates were obtained by a bootstrap procedure (1,000 replicates).

Time-calibrated reconstruction of phylogeny

The time-calibrated phylogeny was reconstructed using the Bayesian inference (BI) in BEAST 1.8.2 (Drummond et al., 2012). Two amphipod species Dikerogammarus villosus (Sovinsky, 1894) (EF570297) and Pontogammarus robustoides (GO Sars, 1894) (JF965990) were used as out-groups for G. roeselii in the following analysis. Given that no molecular clock calibration based on gammarid fossil material is available, the calibration was performed applying three different schemes in parallel: (1) secondary calibration points; (2) standard rate of evolution; and (3) geological calibration points. The first scheme was based on the secondary calibration points from Hou et al. (2011). We constrained the node between Gammarus sp. 3 (JF965959) and Gammarus rambouseki Karaman 1931 (JF965946) to 15.5 Ma. The node between these two species and our sequences of G. roeselii was constrained to 22 Ma. Dating constraints were set to normal distributions with 10% standard deviation. The second calibration scheme applied the “standard” rate of the COI evolution as 0.0115 per site per My (Brower, 1994). The third scheme was based on two geological calibration points: (a) withdrawal of Moesian Sea and subsequent formation of modern Danube river system; (b) origins of ancient Lake Baikal. The calibration points, as well as the prior settings, were described in detail in Mamos et al. (2016). For the first calibrated point (a) we used sequences of G. cf. balcanicus from the Southern Carpathians (KU056394) and the Balkan Mountains (KU056256). The constraint was set to the time between the withdrawal of the Moesian Sea at ca. 2.6 Ma and emergence of the modern Danube system at ca. 1.3 Ma, which connected formerly-isolated rivers of the Southern Carpathians and the Balkan Mountains (Miklos & Neppel, 2010; Mamos et al., 2016). For the second calibration point (b) we used the common ancestor of Acanthogammarids from the ancient Lake Baikal. A permanent lake in the Baikal basin that originated probably ca. 27 Ma, but the conditions in the lake were subsequently changing. Thus, the lacustrine species flocks may be much younger than the lake itself (Mats et al., 2000). According to Mats, Shecherbakov & Efimova (2011), the radiation of Acanthogammarids, represented here by Acanthogammarus brevispinus Dorogostaisky, 1922 (AY926651), Ommatogammarus albinus Dybowsky, 1874 (AY926686) and Eulimnogammarus viridulus Bazikalova, 1945 (AY926665), started at c. 17 Ma. Thus, the constraints were designed to produce results centred on this value but ranging back to the origins of the lake.

In all three calibration schemes the strict clock (vs relaxed lognormal and relaxed exponential), general time-reversible (GTR) model of evolution with gamma-distributed rate heterogeneity (G) and a proportion of invariable sites (I) model of evolution (vs all other models available in BEAST) as well as the Birth-Death speciation process (vs Yule) were set as priors, following the Akaike information criterion (AIC) through the Markov chain Monte Carlo (MCMC) method of moments estimator (Baele et al., 2013) in Tracer 1.6 (Rambaut et al., 2014). Four MCMC, each with 20 M iterations, were run and sampled every 2,000 iterations. The Effective Sampling Size (ESS) of each parameter was verified to be above 200 in Tracer 1.6 (Rambaut et al., 2014). Runs were combined with LogCombiner 1.8.2 (Drummond et al., 2012) with 25% burn-in, and the maximum clade credibility chronogram was annotated using TreeAnnotator 1.8.2 (Drummond et al., 2012). Since all the time-calibrated phylogeny reconstructions provided congruent results (Table S2) we have used the tree obtained from the first calibration scheme as a proxy for following analysis of MOTU delimitation and Bayesian reconstruction of ancestral states. Our interpretation of the obtained results in the context of local paleogeography follows the established and commonly accepted intercalibrations of direct and stratigraphic dating provided by Gradstein, Ogg & Smith (2004) and by Piller, Harzhauser & Mandic (2007), which were summarised for the area by Harzhauser & Mandic (2008).

In order to provide additional support for the BI topology we have also reconstructed a phylogeny using the Maximum Likelihood approach (ML) in MEGA6 software (Tamura et al., 2013). The GTR + G + I model of evolution was used with 1,000 bootstrap replicates (Felsenstein, 1985).

Analysis of diversification rates

The history of diversification was explored using the lineage through time (LTT) plot generated from 1000 BI trees in Tracer 1.6. In order to avoid artificial influence of coalescence on terminal part of the plot we have used BI trees generated with the *BEAST algorithm (Heled & Drummond, 2010) in BEAST 1.8.2. Species delimited for the analysis were coded according to the consensus MOTUs. The BI analysis was conducted using the same priors and MCMC settings as in the first calibration scheme of the time calibrated reconstruction of phylogeny. To identify the potential shifts in diversification we used Laser 2.3 package (Rabosky, 2006) written in R (R Core Team, 2012), a set of maximum likelihood-based methods for analyzing lineage diversification rates. We estimated the fit of six speciation models (function fitdAICrc), assuming constant (pure birth, birth death) and variable (Yule 2-rate, Yule 3-rate, exponential and logistic density dependent) diversification rates, to the reconstructed BI phylogeny.

Cryptic diversity—MOTU delimitation

To explore the number of MOTUs that could represent putative cryptic species within the morphospecies G. roeselii, we applied five different methods: (i) a distance based Barcode Index Number (BIN) System (Ratnasingham & Hebert, 2013); (ii) a distance-based barcode-gap approach using the Automatic Barcode Gap Discovery (ABGD) software (Puillandre et al., 2012); (iii) the tree-based, phylogenetic approach using the Bayesian implementation of the Poison Tree Processor (bPTP) (Zhang et al., 2013) and (iv and v) the tree-based general mixed Yule coalescent (GMYC) model-based method (Pons et al., 2006), using either single and multiple threshold models.

The BIN method is implemented as part of The Barcode of Life Data systems (BOLD; Ratnasingham & Hebert, 2007). Newly submitted sequences are compared together as well as with sequences already available in BOLD. Sequences are clustered according to their molecular divergence using algorithms aiming at finding discontinuities between clusters. Each cluster is ascribed a globally unique and specific identifier (aka Barcode Index Number or BIN), already available or newly created if newly submitted sequences do not cluster with already known BINs. Each BIN is registered in BOLD.

The ABGD method is based upon pairwise distance measures. With this method the sequences are partitioned into groups (MOTUs), such that the distance between two sequences from two different groups will always be larger than a given threshold distance (i.e., barcode gap). We used primary partitions as a principal for group definition, as they are typically stable on a wider range of prior values, minimise the number of false positive (over-split species) and are usually close to the number of taxa described by taxonomists (Puillandre et al., 2012). The default value of 0.001 was used as the minimum intraspecific distance. As there is currently no consensus about which maximum intraspecific distance is reflecting delimitation of species, neither based on morphology (Costa et al., 2007; Weiss et al., 2014, or Katouzian et al., 2016) nor on reproductive barrier (Lagrue et al., 2014) we explored a set of values up to 0.1. The standard Kimura two-parameter (K2p) model correction was applied (Hebert et al., 2003).

The tree based bPTP method for species delimitation uses non-ultrametric phylogenies. The bPTP method incorporates the number of substitutions in the model of speciation and assumes that the probability that a substitution gives rise to a speciation event follows a Poisson distribution. The branch lengths of the input tree are supposed to be generated by two independent classes of the Poisson events, one corresponding to speciation and the other to coalescence. Additionally the bPTP adds Bayesian support values (BS) for the delimited species (Zhang et al., 2013). For the input tree we generated phylogeny using Bayesian inference in MrBayes 3.2 (Ronquist et al., 2012). Prior to tree construction, the out-group was removed from the alignment. The GTR + G + I model was selected using AICM method of moments estimator (Baele et al., 2013) in Tracer 1.6 (Rambaut et al., 2014). Four MCMC, each with 20 M iterations, were run and sampled every 2,000 iterations. This was enough to obtain potential scale reduction factor values equal 1 for all parameters. Stable convergence was determined by inspecting the log likelihood values of the chains and the split frequencies that reach, after 10 million generations, value above 0.01. The consensus tree was constructed after removal of 25% burn-in phase. Analysis was performed on the bPTP web server (available at: http://www.species.h-its.org/ptp/) with 500,000 iterations of MCMC and 10% burn-in.

With the GMYC method, species boundaries are assessed based on sequences in a maximum likelihood framework by identifying the switch from intraspecific branching patterns (coalescent) to typical species branching patterns (Yule process) on a phylogenetic tree. First, a log-likelihood ratio test is performed to assess if the GMYC mixed models fit the observed data significantly better than the null model of a single coalescent species. Providing there is evidence for overlooked species inside the phylogenetic tree tested, two different GMYC approaches, one using the single threshold and the other one on multiple threshold model are used to estimate the boundary between intra- and interspecific branching patterns. Since the GMYC-method requires an ultrametric tree, we have used the already reconstructed Bayesian time calibrated phylogeny. The consensus tree was loaded into the R software package ‘SPLITS’ (Species Limits by Threshold Statistics) (Ezard, Fujisawa & Barraclough, 2009) in R v3.1.0 (R Core Team, 2012) and analysed using the single and multiple threshold models. Presence of significant differences between the two models was tested using likelihood ratio test (LRT).

Bayesian reconstruction of the ancestral states

The program BayesTraits (Pagel, Meade & Barker, 2004) was used to implement the MCMC Bayesian method, aiming at reconstructing ancestral states of a given trait along a phylogenetic tree. Two traits, habitat (river vs lake) and geographic coordinates (longitude and latitude) were explored. The BI maximum clade credibility chronogram was used as the input tree and the actual habitat type and geographic coordinates of its tips were used as input data for the actual state of the trait. For each ancestral node along the tree the results were reported as (i) a chance (in percentage) for the ancestor to occur in a given habitat or (ii) possible coordinates of the ancestor’s locality. In both cases the analysis was performed in two steps. Firstly, using the mentioned input data, a distribution of models was estimated and saved into a file. Secondly, the input data with file containing models distribution was used in the final analysis in which four independent, 10 M iterations long, MCMC runs were performed and sampled every 1,000 iterations, with 5% burn-in. Continuous and multistate reconstruction options with equal priors models, were used for geographic coordinates and habitat, respectively. All runs in each analysis produced convergent results that were combined and summarised (see Table S3 in Supporting Information). Mean and standard deviations for the ancestral geographic coordinates were calculated in MS Excel.

Phylogenetic tree of 28S haplotype and relationship with COI MOTUs

To display the relationships between the nuclear 28S rDNA haplotypes a phylogenetic Neighbor-Joining (NJ) tree based p distance and node robustness being evaluated by 1,000 bootstraps was constructed using MEGA 6.0 (Tamura et al., 2013). To display relationship (i.e., discrepancy between phylogenies) between nuclear and mitochondrial phylogenies COI MOTUs were plotted on 28S NJ tree.

Within MOTU diversity, divergence, differentiation and historical demography: a Pleistocene tale

Molecular genetic diversity, divergence, differentiation and historical demography based on mtDNA COI sequences were, wherever possible, estimated for each MOTU. The relative frequency, geographic distribution and molecular divergence between COI haplotypes, were graphically presented as a median-joining (MJ) network (Bandelt, Forster & Rohl, 1999) implemented in the program NETWORK 4.6.1.2 (http://www.fluxusengineering.com). Given the differences in substitution rates, we applied a weight of 1 for transition and 2 for transversions. The topology was obtained at the homoplasy level parameter default value (e = 0).

Diversity was assessed as the number of haplotypes (k), haplotypic diversity (h) and nucleotide diversity (π) (Nei, 1987) with the DnaSP software (Librado & Rozas, 2009). Molecular divergence was estimated as average Kimura 2—parameters (K2p) distance between haplotypes using MEGA 6.0 (Tamura et al., 2013). In case of MOTUs present in more than one site (MOTUs A, C, E, G and K) and given a minimum sampling size per site of 3 individuals, the differentiation in haplotype frequencies was estimated as overall F_ST through an Analysis of Molecular Variance (AMOVA, Excoffier, Smouse & Quattro, 1992) with 10,000 permutations. In addition, pairwise differentiation between locations was estimated through F_ST (Weir & Cockerham, 1984) with 1,023 permutations, and through Exact Tests of Sample Differentiation (ETSD, Raymond & Rousset, 1995) with 30,000 Markov steps. Historical demographic expansion within MOTUs at both the scale of MOTUs and location was examined based on Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1996) neutrality tests with 1,000 replicates. All the above analyses were performed in Arlequin 3.5.1.3 (Excoffier & Lischer, 2010).

Species geographic distribution and associated habitat

We collected Gammarus roeselii from 26 sites on the Balkan Peninsula, including 18 riverine and 8 lacustrine sites belonging to 12 river basins draining to the Black, Adriatic, Ionian and Aegean seas as well as the isolated Prespa lake (Table 1, Fig. 1). The altitudinal range of the sites spans from the sea level to nearly 870 m.a.s.l. The longitudinal range of sites in which G. roeselii was found extended from the Adriatic/Ionian coast of the Balkan Peninsula (15.7°E) on the west to the Struma river system (24°E), and the latitudinal range extended from the Drava River (46°N) down to the Corinthian Bay (38.5°N). Despite extensive search, G. roeselii was not found neither in the eastern part of the Balkan Peninsula nor on the Peloponnese Peninsula. Other areas where we did not encounter presence of G. roeselii included central Serbia, Bosnia & Herzegovina and coastal regions of Croatia.

COI haplotypes: diversity, distribution and divergence

Out of 177 individuals sequenced, a total of 74 haplotypes were identified (Accession numbers KP780673 –KP780746) (Table 1). Only seven haplotypes were shared between a pair of locations (Table 1). Each of the remaining 67 haplotypes was specific to only one location (Table 1). Minimum, average and maximum K2p distance between haplotypes were, respectively, 0.0019 (SE = 0.0017), 0.1463 (SE =0.0108) and 0.2405 (SE =0.0108). Such high average K2p distance suggests presence of cryptic diversity within the morphospecies (see also Table S4).

Time calibrated reconstruction of phylogeny and analysis of diversification rate

To infer the chronology of evolutionary events of diversification in G. roeselii, we performed a Bayesian time-calibrated reconstruction of phylogeny based on the COI marker. Both the Bayesian inferred (BI) chronogram and the Maximum Likelihood (ML) tree constructed in parallel had merely the same topology and high node supports (Fig. 2). The mean substitution rate associated with BI chronogram was 0.0113 in substitutions per site, per My, (SD 0.0016, 95% HPD 0.0083–0.0144) for first calibration scheme and 0.0129 substitutions per site, per My (SD 0.0030, 95% HPD 0.0077–0.0193) for third calibration scheme (Table S2).

The spatio-temporal pattern of lineage diversification within G. roeselii in the Balkan Peninsula appeared to be very complex (Fig. 2A, Table S1). The diversification started ca. 18 Ma and produced a high number (>10) of extant lineages. Spatial distribution of these lineages is heterogeneous, most of them are local endemics, two are relatively widespread and few occur in sympatry (Fig. 2B). The initial divergence (ca. 18 Ma) produced two groups of lineages—one spread over the northern part of the Balkan Peninsula, while the other one colonised southern part of this region. Each of the two groups is composed of a set of lineages that diverged before the beginning of Pleistocene glaciations, namely in Miocene (ca. 20.0–5.3 Ma) and Pliocene (ca. 5.3–2.6 Ma). The lineages through time plot (LTT, Fig. 3) shows that the accumulation of lineages was rather constant until Pleistocene. The analysis of diversification rate change through time showed that constant rate pure birth model and density-dependent speciation rate model with logistic variants fit the data best, with only a little difference between them (dAIC =0.21).

Cryptic diversity/MOTU definition

The thirteen MOTUs (A–M, Fig. 2) representing a set of potential cryptic species, we retained, are based on the combination of the information brought by (i) five species delimitation methods (BINs, ABGD, bPTP, GMYC single threshold and GMYC multiple threshold model) and (ii) the divergence age on the BI chronogram (Fig. 2).

All the 177 sequences submitted to BOLD clustered into 27 BINs (details in Table S1).

The ABGD distance-based approach (barcode-gap analysis), resulted in eight primary partitions of the analysed COI sequences. Four of them were stable over a wide range of P values (0.001–0.10) and defined 12 MOTUs, being the same in each partition (Fig. 2).

The estimated number of MOTUs in bPTP ranged between 14 and 29, out of which 20 have Bayesian support values above 0.46. Eight of them are represented only by a single haplotype (Fig. 2).

In the GMYC, the null model of a single species was rejected in both, single and multiple threshold models (GMYC single LRT < 0.002, GMYC multiple LRT < 0.0001). Results of the chi-square test showed that there was no statistically significant difference between the two models (LRT =6.305, p = 0.39). Results of the single threshold analysis suggested the presence of 19 MOTUs (confidence interval: 12–30), out of which five consisted of single haplotypes. The results of the multiple threshold analysis suggested 24 MOTUs (confidence interval: 17–28), out of which seven consisted of single haplotypes (Fig. 2).

The 13 MOTUs retained, corresponded to the MOTUs identified by the ABGD method with the only exception that we decided to split one ABGD MOTU in two MOTUs (D + E). They diverged at ca. 4 Ma, already in the Pliocene, and the K2p distance between these MOTUs is ca. 6%. Such a distance might have been associated with reproductive barrier observed in other Gammarus species (Lagrue et al., 2014) (Table S4, see Discussion).

Chronology of divergence and geographical distribution of MOTUs

The identified MOTUs formed two well defined primary clades. The first clade grouped MOTUs (A–H) occupying the northern and the central part of the Balkan Peninsula. The second one included MOTUs (I–M) present predominantly in the southern part of the peninsula. The clades appear to have diverged in the early Miocene (ca. 18 Ma) (Fig. 2A).

Divergence within both clades started ca. 13 Ma, leading to separation of MOTU H from A to G in the northern clade and I–J from K to M in the southern one. Next, within the northern clade, at ca. 8 Ma, the group including MOTUs A–G separated into two subgroups A–E and F–G. The former, at ca. 7 Ma, gave rise to MOTU sets A–C and D–E. The MOTU C separated from the others at ca. 4.5 Ma, while A and B diverged at ca. 4 Ma. MOTUs D and E separated at the same time. Of the other subgroup, MOTUs F and G diverged at ca. 7 Ma. Within the southern clade, MOTUs I and J diverged already at ca. 10 Ma. Then MOTU K separated from L to M at ca. 7 Ma and, finally, MOTU L diverged from M at ca. 3.3 Ma (Fig. 2A).

Currently, for the northern clade, MOTU A had the widest distribution. It was present in seven sites, including rivers and lakes in the central Balkan Peninsula from the Aliakmonas drainage on the southwest, to the Vardar drainage on the north, and the Struma drainage on the east (Fig. 2B). MOTU A coexisted with MOTU D in the Vardar River drainage, the only place where MOTU D was found, and with MOTU E in the Struma River. The lacustrine populations of MOTU A inhabited Kastoria, Vegoritis and Petron lakes; all belonging to the Aliakmonas drainage. MOTU B was found only at one site in the southernmost part of the South Morava drainage, a major Danube tributary (Fig. 2B). MOTU C showed the northernmost distribution, being present only in rivers of the Danube drainage system (Black Sea basin), on the Pannonian Basin. MOTU E was another widely distributed lineage, found in the central part of the Balkan Peninsula. Interestingly, it was always found in sympatry with other MOTUs—with A in the Struma River, with G in the Shkumbin River, and with J in the upper Devoll River (a tributary of the Seman River). MOTU G, present in six sites, was predominantly associated with the Drin River drainage, including also the Ohrid and Prespa lakes, but also the nearby Shkumbin River, draining directly to the Adriatic Sea. MOTU H was found at only one isolated site in a small river on the Adriatic side of the Hellenic Mountains (Fig. 2B).

The MOTUs belonging to the southern group are spatially scattered and have very localised distributions (Fig. 2B). Each of them is restricted usually to only one waterbody or drainage area. For example, the highly divergent MOTU I occupies the southernmost locality. It was found only in the Trichonis Lake, on the Ionian side of the Hellenic Mountains. Its sister MOTU J was present only in a remote locality on the upper Devoll River, in the north-eastern outskirts of the Hellenic Mountains. MOTU K was found in the Sperchios drainage on the Thessalian Plain and in the neighbouring Sofaditikos drainage in the Pelion Mountains. MOTU L was found in the Pinios drainage, also in Thessaly. Interestingly, its sister MOTU M was found only in the remote Thyamis drainage, on the Adriatic side of the Hellenic Mountains.

Reconstruction of the ancestral longitudes/latitudes and ancestral habitats of MRCAs

Reconstruction of the geographical origin of the most recent common ancestors (MRCAs) (Fig. 4A, Table S3) allowed to project the lineage diversification within G. roeselii on the present geography of the Balkan Peninsula and to interpret it in the paleogeographic context. The analysis showed that the early Miocene diversification of G. roeselii lineages might have started in the western part of the central Balkan Peninsula. From there, the northern clade could have spread to the neighbouring areas of the present Danube, Drin, Vardar, Shkumbin, Aliakmonas and Struma drainages; where it diversified throughout Miocene and Pliocene. The southern clade had colonised the southwestern parts of the Balkan Peninsula, diversifying in the Miocene and Pliocene, most probably on the eastern side of the Hellenic Mountains, before they effectively crossed the mountain chain and colonised now isolated watersheds. Finally, in the Pleistocene, a number of local diversification events occurred in MOTUs of both clades all over the species distribution on the Balkan Peninsula.

According to the reconstruction of the ancestral habitats of MRCAs (Fig. 2A), most of the Miocene ancestors could, with equal probability, have live both in lacustrine and riverine habitats. In several cases (MOTUs A, F and G), the extant riverine populations derived most probably directly from lacustrine ancestors. In case of MOTUs A and G the lacustrine ancestors of riverine populations were of Pleistocene age, while in the case of MOTU G, the presumed ancestor lived in the late Miocene.

Projection of the reconstructed longitudes and latitudes on paleomaps (Figs. 4B–4D) suggests that the early MRCA of extant G. roeselii MOTUs occurred in the area covered by the Neogene lake system occupying large parts of the western Balkan Peninsula. The reconstructed positions of the Miocene/Pliocene regional diversifications coincide with the major orogenic uplifts in the area, causing, among others, a decline of the lake systems and complete rearrangement of the local hydrological area.

28S: diversity, phylogeny and evidence for introgression?

For each of the 13 MOTUs identified by the COI data, one to 18 individuals (50 in total) were sequenced for the nuclear 28S marker, producing 23 haplotypes (N1–N23) (Table 1). The topology of the Neighbor-Joining tree (Fig. 5) showing relationships between the 28S haplotypes was partially concordant with the results derived from the mtDNA COI marker. For example, MOTUs with high COI divergence (ca > 15% K2p distance) were also divergent in the case of 28S. Only MOTU E showed a peculiar pattern as its individuals, associated with three 28S haplotypes (N2, N4 and N5), and belonged to two highly divergent clades of the 28S N–J tree (Fig. 5). Interestingly, MOTU E was always sympatric (sites 11, 13, 20) with another distantly related MOTU G, J and A, respectively. In two cases the individuals with COI haplotypes of MOTU E had a nuclear haplotype either close (site 11, N2 vs N3 and N9) or even identical (site 13, N4 and N5) to the nuclear haplotypes associated with the other MOTU present at the respective site.

COI: diversity and demography within MOTUs

Analysis of genetic diversity within all 13 defined MOTUs (Fig. 6), revealed a generally rather high level of haplotypic and nucleotide diversity as well as K2p distance from 0.002 to 0.029 (Table 2). This was also reflected by the fact that in numerous cases various species delimitation methods have subdivided the MOTUs we defined into smaller units. In all cases this intra-MOTU divergence is younger than 2.5 Ma. Such sub-MOTUs have very localized distributions, e.g., they are restricted to one lake (e.g., Prespa, Vegoritis/Petron) or one lake system (Kastoria, Ohrid). This suggests that the Pleistocene glaciations might have promoted diversification at the local scale. Only MOTU M showed a clear sign of demographic expansion (Table 2, Table S6). On the other side, results of AMOVA (Table 2), pairwise F_ST and Exact Test of Sample Differentiation (ETSD) revealed significant differentiation between most populations of the same MOTUs inhabiting various locations (Table S7).

High level and pre-Pleistocene onset of cryptic diversity

Our findings point out the presence of 12–29 different MOTUs, depending on the delimitation methods, that may represent plausible species. The higher estimates by GMYC and bPTP. which are methods based on phylogeny and coalescence, may be associated with their general tendency to overestimate the number of MOTUs in the presence of numerous highly divergent haplotypes, or when MOTUs appear to be unevenly sampled; with some containing low, and others high, within-MOTU genetic variability (Talavera, Dincă & Vila, 2013; Zhang et al., 2013). On the other hand, aggregating the haplotypes into BINs is based on a genetic distance threshold of 2.2% (Ratnasingham & Hebert, 2013) that may be too low for Gammarus (for details see Lagrue et al., 2014). Thus, it seems that the most conservative results of the ABGD analysis provide more realistic estimates. A detailed study upon the taxonomy of Gammarus roeselii is beyond the scope of our paper. In the only taxonomic review of European Gammarus made so far, and based exclusively on morphological characters, Karaman & Pinkster (1977a), Karaman & Pinkster (1977b) and Karaman & Pinkster (1987) synonymised nine species and subspecies and recognised them as populations of G. roeselii. Our results are in clear opposition to their conclusion and suggest that, in reality, G. roeselii is a conglomerate of several locally endemic species. For example, our study revealed that MOTUs A and D were found in the Treska River of the Vardar system, while MOTUs A and E were found in the Struma River. Then, it is quite likely that at least some previously synonymised species, such as Gammarus vardarensis (Karaman, 1929), originally described from the Vardar River system, and Gammarus strumicae (Karaman, 1974), described from the Struma River system, correspond to some of these MOTUs. This points that another revision of G. roeselii should be made in the future, implementing the integrative taxonomy approach, which is the only reasonable option in the case of diverse animal groups with conservative morphology or poorly defined morphological variability (Pante, Schoelinck & Puillandre, 2015). High cryptic diversity with numerous divergent, locally endemic MOTUs representing possibly cryptic and/or yet undescribed species within other European freshwater amphipods was already evidenced for the morphospecies G. balcanicus, which is widely spread from the eastern Carpathians through the Balkan Peninsula and to the eastern Alps (Mamos, 2015; Copilaş-Ciocianu & Petrusek, 2017). Accordingly, a high level of cryptic diversity was found in case of G. fossarum inhabiting Carpathians, Balkans, Alps and other upland areas in Europe (Westram et al., 2011; Weiss et al., 2014; Copilaş-Ciocianu & Petrusek, 2017).

Mountainous areas with steep elevational gradients, such as the Balkan Peninsula, have often been reported as local endemism and diversity hotspots of biota (Hou, Li & Li, 2014; Luebert & Muller, 2015; Zhou et al., 2015; Katouzian et al., 2016). In southern Europe it was often a result of speciation by isolation in consequence of severe alterations and fragmentations of habitats during Pleistocene climatic oscillations (Hewitt, 2004; Salicini, Ibanez & Juste, 2013; Goncalves et al., 2015). However, the presence of pre-Pleistocene lineages within conventionally recognised morphospecies was reported for several terrestrial and fewer freshwater taxa, including: reptiles; amphibians; and fish, but also insects and molluscs (e.g., Albrecht et al., 2007; Marková et al., 2010; Falniowski et al., 2012; Garcia Munoz et al., 2014; Previšić et al., 2014; Pabijan et al., 2015).

Our study revealed that the level of genetic divergence between most of the Gammarus roeselii lineages in the Balkan Peninsula is too high to be attributed to climatic fluctuations in the Pleistocene. Similarly, Mamos et al. (2016) as well as Copilaş-Ciocianu & Petrusek (2017) found that the divergence within G. balcanicus started some 15–20 Ma, in early/middle Miocene. Copilaş-Ciocianu & Petrusek (2015) found that G. fossarum diversified in the south-western Carpathians 10–17 Ma. A high cryptic diversity of pre-Pleistocene origin was also recently discovered within the freshwater isopod, Asellus aquaticus (Linnaeus, 1758), widely distributed in Europe including the Balkan Peninsula. The MOTUs of A. aquaticus diverged over Miocene and a number of them are localised in various parts of the Balkan Peninsula (Verovnik, Sket & Trontelj, 2005; Sworobowicz et al., 2015), which is similar to the spatial distribution pattern we have obtained for G. roeselii.

The impact of Neogene and Quaternary paleogeography on the diversification of Gammarus roeselii

Historical framework of cryptic diversity within the Balkan freshwater taxa has been poorly studied, with a clear exception of taxa inhabiting ancient lakes, such as: the Ohrid; Prespa; or Trichonis (Albrecht et al., 2008; Wysocka et al., 2013; Wysocka et al., 2014). This is probably a result of the dynamic paleogeography and geology of the area, which may restrain a historical interpretation of phylogenetic relationships between MOTUs and their observed spatial distribution patterns, making it an uneasy endeavour.

Our BI chronograms based on two different calibration schemes provided estimations of mean COI substitution rates, in substitutions per site, per million years, which did not differ substantially (0.0113; SD 0.0016 vs 0.0129; SD 0.0030). They were also quite close to the COI mutation rate (0.0165; SD 0.0018) obtained for the cryptic lineages of G. balcanicus by Copilaş-Ciocianu & Petrusek (2017) and to the rate (0.0115) routinely used for arthropods (Brower, 1994). Thus, we can assume that our estimation of divergence times were reasonable and reliable in comparison to other studies. According to our results, the divergence of G. roeselii MOTUs took place within similar time frame as obtained in the other studies (reported above) upon Gammarus. The two main groups of MOTUs split at ca. 18 Ma (early middle Miocene) in the Neogene, and this continued their diversification to individual MOTUs until less than 1 Ma (Pleistocene) in the Quaternary. During this time frame, the region was already a part of continental Europe; undergoing rapid changes of terrain due to the uplift of the Dinaric and Hellenic mountain chains. This included high tectonic activity of the area and recurrent reorganisations of local hydrological networks (Popov et al., 2004; Dumurdzanov, Serafimovski & Burchfiel, 2005). In the case of other gammarids (Copilaş-Ciocianu & Petrusek, 2015; Copilaş-Ciocianu & Petrusek, 2017; Mamos et al., 2016), we can assume that such landscape changes driven by Alpine orogenesis had a profound impact on the diversification of G. roeselii in the Balkans. A major feature of the Balkan landscape was large freshwater lacustrine systems, well supported by sediment and fossil data, existing throughout the Neogene and undergoing dynamic changes due to local orogenic movements (Harzhauser & Mandic, 2008; Krstić, Savić & Gordana, 2012). Unfortunately, we could not unambiguously define the ancestral habitats of Miocene common ancestors for the extant MOTUs; they could have lived in either rivers or lakes. However, reconstruction of the geographical origin of the Miocene and Pliocene ancestors suggests that they lived most probably in the areas occupied contemporarily by large freshwater lakes (Popov et al., 2004; Harzhauser & Mandic, 2008; Krstić, Savić & Gordana, 2012). This allows us to speculate that they could have been lacustrine organisms that diversified gradually along with the regression of the lakes and formation of new hydrological networks.

At least some of the extant ancient, or presumably ancient, Balkan lakes are probably relicts of these paleolake systems (Albrecht & Wilke, 2008). Populations of G. roeselii were found in several of them, such as: Ohrid; Prespa; Trichonis etc. (see Table 1, Fig. 1). The lacustrine populations showed various levels of divergence from their sister riverine counterparts, dating from as early as ca. 10 Ma (Miocene) to as late as ca. 1 Ma (Pleistocene). In a couple of cases, populations of both types shared the same or closely related haplotypes. Our results for the ancestral habitat reconstruction suggest that for MOTUs A, F and G, the extant riverine populations derived from lacustrine ancestors. For example, the closest relatives of haplotypes found in rivers of Northern Greece (MOTU A) are those from the local lakes of Kastoria, Petron and Vegoritis, which are connected to the Aliakmonas River system. Similarly, in MOTU G, haplotypes from the Shkumbin and Drin rivers in Albania are most closely related to those from the ancient lakes Ohrid and Prespa. Some individuals from the Ohrid Lake share haplotypes with those from the River Drin. The lake drains to Drin and, possibly, was historically connected to the nearby Shkumbin River. Interestingly, for haplotypes from the currently isolated Lake Prespa, they form a sister clade to all the other haplotypes of MOTU G, from which they diverged 1.5 Ma. Similar sister relationships are also known e.g., within the gastropod genus Radix inhabiting both lakes (Albrecht et al., 2008). It is also known that Lake Ohrid and Prespa were at once part of the large Dessarete lake system (Albrecht & Wilke, 2008). The most recent studies of the lake sediment cores indicate that Lake Ohrid is at least 2 My old (Jovanovska et al., 2015). The age of Lake Prespa remains unknown but it may be at least as old (Karaman, 1971; Radoman, 1985). A split of that clade may reflect a separation of the Dessarete lakes and a relict character of the population inhabiting rivers derived from the Dessaretes. An interesting case is whether MOTU I is endemic to the ancient Lake Trichonis. This lake harbours several endemic species of molluscs and fish (Doadrio & Carmona, 1998; Albrecht et al., 2009; Almada et al., 2009). Its exact age is not known, yet it is hypothesized to originate in the Pliocene, at ca. 3.4–1.8 Ma, or in the early Pleistocene (Economidis & Miller, 1990; Khondkarian, Paramonova & Shcherba, 2004). The Pleistocene origin of MOTU I diversity in the lake is supported by our results (Fig. S1, Table S3). Yet it is worth noting that it diverged from its sister, riverine MOTU J, at ca. 10 Ma. Given the presumably Pliocene/Pleistocene age of Lake Trichonis, and results of our ancestral state analysis, we may assume that MOTU I could derive either from a riverine or lacustrine ancestor. In fact, almost all the Miocene ancestors in our analysis had a rather equal chance to be either lacustrine or riverine. This may suggest that multiple habitat switches are evident in the evolutionary past of G. roeselii, and this points out the potential importance of lacustrine systems for their diversification history. Climate aridisation that had started in the Miocene and intensified greatly during the Pleistocene glaciations led to a periodical depletion of the riverine network in southern Europe, as seen in the present day genetic diversity of freshwater biota (Tzedakis, 2007; Macklin, Lewin & Woodward, 2012; Gonzalez, Pedraza-Lara & Doadrio, 2014). In such climatic conditions, the lakes could serve as local microrefugia, working as sources for the recolonisation of river systems, restored during more humid periods. That is particularly plausible for the Pleistocene, when recurrent glaciations resulting in droughts intermingled with humid interglacials. Indeed, our results indicate that the direct lacustrine ancestors of extant riverine populations for MOTUs A and G are of Pleistocene age. Similar conclusions may be derived from the results of the F_ST and ETDS analyses showing significant differentiation between localities. Also, the aggregation of haplotypes into BINs may provide support for the Pleistocene speciation by an isolation process, as observed in other taxa (Salicini, Ibanez & Juste, 2013; Goncalves et al., 2015). It seems also that at least some MOTUs are versatile in changing habitats. They are present in both lakes and rivers, and we may also assume that changing between such habitats played an important role in the evolutionary history of G. roeselii.

Evidence of potential introgression between currently divergent MOTUs

In several examples, we observed incongruence between phylogenies based on the nuclear (28S) and mitochondrial (COI) markers. Similar conflict between nuclear and mitochondrial genetic diversity patterns was evident for the Balkan cyprinid fish, suggesting an exchange of genes after initial isolation (Marková et al., 2010). Nuclear haplotypes shared randomly among numerous mitochondrial lineages were observed for A. aquaticus and interpreted as a result of incomplete lineage sorting (Verovnik, Sket & Trontelj, 2005; Sworobowicz et al., 2015). However, in the case of G. roeselii such interpretation is implausible. Two nuclear haplotypes were shared only between MOTUs E and J, which were not direct relatives and diverged from each other ca. 18 Ma. A third nuclear haplotype of MOTU E was a close relative of a haplotype from MOTU G; while, based on the COI sequence, the two MOTUs belonged to different clades that had separated ca. 8 Ma ago. Moreover such a situation was observed only in the localities where the MOTUs occurred in sympatry; suggesting that the female harboring nuclear haplotypes of MOTU E were fertilised by males from the other locally present MOTUs. According to Morando et al. (2004) such phylogenetic patterns indicate secondary contact, resulting in hybridisation and/or introgression. It is hard to speculate what enabled such secondary contacts in the case of G. roeselii. MOTU E emerged at ca. 3.5 Ma, already in the Pliocene, and is the most diversified and internally divergent (five COI BINs) of all the defined MOTUs, and is also present in three distant localities. Thus, two main scenarios would be possible. One would assume a much wider presence of the MOTU before the Quaternary as well as an extinction event for much of the former distribution area due to Pleistocene climate aridisation resulting in a depletion of the river systems. According to this scenario, secondary contact with other MOTUs would be possible via recurrently changing watersheds during contemporary climate oscillations. A second scenario would assume a possible overland transport by birds. Such possibility was proved for freshwater amphipods living in lakes and shallow lowland rivers occupied by waterfowl (e.g., Segerstråle, 1954; Rachalewski et al., 2013). However it cannot be excluded, we believe it is less probable given that MOTU E was found exclusively in rivers in rather mountainous terrain.