Fuzzy species borders of glacial survivalists in the Carpathian biodiversity hotspot revealed using a multimarker approach

The Carpathians are one of the key biodiversity hotspots in Europe. The mountain chain uplifted during Alpine orogenesis and is characterised by a complex geological history. Its current biodiversity was highly influenced by Pleistocene glaciations. The goal of the current study was to examine the phylogenetic and demographic history of Gammarus balcanicus species complex in the Carpathians using multiple markers as well as to delimit, using an integrative approach, and describe new species hidden so far under the name G. balcanicus. Results showed that divergence of the studied lineages reaches back to the Miocene, which supports the hypothesis of their survival in multiple micro refugia. Moreover, the increase of their diversification rate in the Pleistocene suggests that glaciation was the driving force of their speciation. The climatic changes during and after the Pleistocene also played a major role in the demography of the local Carpathian lineages. Comparison of diversity patterns and phylogenetic relationships of both, the mitochondrial and nuclear markers, provide evidence of putative hybridisation and retention of ancient polymorphism (i.e., incomplete lineage sorting). The morphological examination supported the existence of two morphological types; one we describe as a G. stasiuki sp. nov. and another we redescribe as a G. tatrensis (S. Karaman, 1931).

The Carpathians are one of the key biodiversity hotspots in Europe. The mountain chain uplifted during Alpine orogenesis and is characterised by a complex geological history. Its current biodiversity was highly influenced by Pleistocene glaciations. The goal of the current study was to examine the phylogenetic and demographic history of Gammarus balcanicus species complex in the Carpathians using multiple markers as well as to delimit, using an integrative approach, and describe new species hidden so far under the name G. balcanicus. Results showed that divergence of the studied lineages reaches back to the Miocene, which supports the hypothesis of their survival in multiple micro refugia. Moreover, the increase of their diversification rate in the Pleistocene suggests that glaciation was the driving force of their speciation. The climatic changes during and after the Pleistocene also played a major role in the demography of the local Carpathian lineages. Comparison of diversity patterns and phylogenetic relationships of both, the mitochondrial and nuclear markers, provide evidence of putative hybridisation and retention of ancient polymorphism (i.e., incomplete lineage sorting). The morphological examination supported the existence of two morphological types; one we describe as a G. stasiuki sp. nov. and another we redescribe as a G. tatrensis (S. Karaman, 1931).
European freshwaters are inhabited by numerous amphipod species classified in the genus Gammarus Fabricius, 1775. Extensive studies upon the diversity and taxonomy of this speciose genus in Europe were undertaken in the early twentieth century (e.g., [1][2][3][4] ). Subsequently, although 30 years passed since then, these efforts were summarized and critically revised in a series of works by G. Karaman and Pinkster [5][6][7][8] , with the result of the revision synonymizing most of the previously described species. Since that time, however, more than a dozen of new species were described from inland waters of Europe and adjacent regions (e.g., [9][10][11][12][13][14][15][16]. Furthermore, during the recent decade, phylogenetic studies based on extensive molecular data sets have entirely altered traditional views on the taxonomy of the family Gammaridae and also on the definition and validity of genera belonging to this family [17][18][19][20] . Also, recent phylogeographic studies on several geographically widespread Gammarus morphospecies pointed out to outstanding cryptic and pseudo-cryptic diversity with an ancient divergence between phylogenetic lineages (e.g., 14,15,[20][21][22][23][24][25] ). Unrevealing such a high level of hidden diversity stimulated further taxonomic studies, resulting in the description of new taxa, based either on a modern integrative approach using morphological, ultrastructural and molecular characters or, in the case of species for which the morphological distinction was impossible, based on molecular descriptive characters only (e.g., 14,15 ). However, the latter approach has been disputed for some years. Although it has several shortcomings (for discussion see 26 ), it has been appreciated as a tool that helps to overcome the evident global taxonomic impediment (e.g., 27,28 ). Although the hypotheses of the species described in this way are viewed as interim, particularly if they are based only on a single molecular marker, they are as valid as those founded on morphological features. Such hypotheses can subsequently be verified using a multiproxy integrative approach 28 . Importantly, nowadays they may be used in broadly defined nature conservation, e.g., for fast, effective and precise quantification of biodiversity and detection of endemic lineages in the presumed diversity and endemism hotspots 25 . MOTUs delimitation. Eight molecular operational taxonomical units (MOTU) delimitation methods provided partially similar results (Fig. 1). The BOLD clustering method revealed 21 BINs within G. tatrensis s.l. and 3 BINs within the G. stasiuki sp. nov. Current study provided also 12 additional BINs belonging to the sister lineage of the studied Gammarus species. The ABGD and ASAP provided similar results, among 10 resulting partitions P was stable for 0.0269-0.0437 (distribution of pairwise distances given in Fig. S2), with ASAP species threshold for preferred partition being: 0.0447 p-distance, the data was divided into 18 MOTUs (G1-G18). Six MOTUs fell within the G. tatrensis (G1-G4, G11, G12) and one within the G. stasiuki sp. nov. Sister lineages were classified in 12 MOTUs. There was no statistically significant difference between single and multiple GMYC methods (Chi 2 = 4.654517, df = 6, p = 0.5888235). Both methods supported the hypothesis of multiple species (sGMYC: likelihood ratio = 44.497, p = 2 × 10 -10 , mGMYC: likelihood ratio = 49.152, p = 2 × 10 -11 ). The sGMYC and mGMYC revealed 25 and 23 MOTUs respectively, for G. tatrensis, 3 for G. stasiuki sp. nov. and 14 for sister lineages. The PTP methods revealed different values for G. tatrensis: 21 (bPTP) vs. 11 (mPTP) as well as 14 bPTP vs. 9 mPTP in sister lineages. In the case of G stasiuki sp. nov., the number of MOTUs (3) was consistent for these two methods. The lowest posterior probability for MOTU grouping within bPTP was 0.43 while for mPTP it was 0.73. Multilocus species delimitation (STACEY) suggested the presence of 36 MOTUs (fraction 0.81) identical to BINs used as minimal clusters, from which 21 fell into G. tatrensis and 3 to G. stasiuki sp. nov. The same results were reported in additional analysis employing path sampling on multilocus data. The BIN delimitation had the highest marginal likelihood in comparison to delimitation using ABGD/ASAP or morphology (Table S2).
Time calibrated reconstruction of phylogeny and distribution. The saturation test revealed no significant loss of phylogenetic signal for any of the tested markers (ISS < ISS.cSym, p < 0.00). The Neighbor-Joining distance tree reconstructed for each marker separately resulted in weakly resolved topologies for nuclear markers with few well-supported clades and strong phylogenetic signal for mitochondrial markers, especially COI. The analysed markers did not convey conflicting phylogenetic signals (Fig. S1). The multi-marker reconstruction of phylogeny (Fig. 2a) showed that both G. tatrensis and G. stasiuki sp. nov. form well supported monophyletic clades. Both species separated from other lineages of G. balcanicus species complex already in the Middle Speciation rate changes through time. Bayes factors showed that the model with one rate shift is preferred (BF > 6) over the null model with 0 shifts. The shift of diversification occurred at the divergence of G. tatrensis Carpathian lineages in the Late Miocene (Fig. 2c, Fig. S3). The BAMM plot illustrating changes in diversification rates shows an increase in the Late Miocene and the following fluctuation ended at the beginning of the Pliocene with strong increment (Fig. 2c). These increments are also visible on the LTT plot (Fig. S3), however, the pick of diversification there is not so prominent.
Demographic analysis. All MOTUs of G. tatrensis revealed by ABGD show some signs of postglacial demographic expansion (Fig. 3). This is particularly evident for MOTU G1 where population growth after putative bottleneck event is supported statistically. MOTU G4 shows only mild population size increment but the significance of neutrality tests suggests also an impact of bottleneck effect on its molecular diversity. Population size growth after decline shown on eBSP of MOTU G2 is not supported by neutrality tests. MOTU G2 is showing a possible bottleneck during the most recent glacial maximum. Lack of statistical support for demographic growth is also observed for MOTU G3. Gammarus stasiuki went through a bottleneck during the Last Glacial Maximum, and the population size has not fully recovered, but these changes are not significant.   Etymology. We name this species in honour of Andrzej Stasiuk, a very successful and internationally acclaimed Polish writer, journalist and literary critic. By this we pay tribute to his travel literature and essays that describe the natural and cultural environment of Eastern Europe, including the Carpathians, where he has chosen to settle.   www.nature.com/scientificreports/ Lower margins of basal A 2 articles 4 and 5 richly setose, with 5-6 groups of 2-4 setae. Length of longest setae sub-equal to the width of 4th article and somewhat longer than the width of 5th article (Fig. 5b). Posterior margins of gnathopod's carpus and propus richly setose (Fig. 5d,e). P 3 merus posterior margin with 4-5 groups of numerous setae as long or a bit longer than merus width (Fig. 5a). Basis of P 6 and P 7 with posterior margin crenulated, in crenules a small setule. Distal part of basis of these pereopods 1.5 times wider than ischium width. Distoposterior lobe of basis P 6 and P 7 sub-rectangular. Posterior margin of P 6 basis somewhat concave, posterior margin of P 7 basis slightly convex (Fig. 5e,f). Second epimeral plate 2 (E2) rather characteristic, with lower margin distinctly convex ending with a small tooth; posterior margin of E2 slightly to distinctly convex (Fig. 6a). Urosomites 1-3 with two medial and two lateral groups of 1-3 spines and/or 1-2 setules (Fig. 6b). Uropod 3 biramous, endopodite length ca. 2/3 of exopodite length (Fig. 6c). Outer margin of the first article of U 3 exopodite with 3 groups of 1-2 spines and 1-4 setae, some a bit longer than spines. Apically this first article of exopodite with 3-5 spines and several setae, some are longer than spines. Second exopodite article apically with 3 setae, the longest as long as this article. Endopodite of U 3 apically with 1-2 spines and 3-4 setae, some over 2 times longer than spines. Inner margin of U 3 exopodite with several groups of 1-3 setae, ca. half of these setae are feathered. Outer margin of U 3 endopodite with 3-5 groups of spines and setae, of which 3-4 are feathered. Inner margin of U 3 endopodite with 2-3 groups of setae, several setae are feathered (Fig. 6c). Feathered setae are sometimes broken. Telson lobes with apical 2-3 spines and 3-4 setae, some a bit longer than spines. On the surface of telson lobe 1-4 subapical and/or subbasal setae (Fig. 6d), very rarely 1 spine. Female (Figs. 8,9): max. length observed 11 mm. Clear sexual dimorphism in the setation of appendages: setae on the lower margin of A 2 peduncular articles 4 and 5 distinctly longer than in males; in 4th article setae are 1. 5 × longer than this article width and in 5th article two times longer than this article width (Fig. 8c). Similarly, in females, P 3 merus is more setose, with setae twice as long as the merus width (Fig. 8g). In the U 3 exopodite outer margin longest setae are 2 times longer than spines. Also, apical setae of U 3 endopodite can be 3 times longer than accompanying spines (Fig. 9c). The body surface, plate margins and appendages in juvenile specimens are less dressed with spines and setae than in adults.  (Fig. 10c). Posterior margin of merus of pereopod 3 with 3-4 groups of setae, the longest is equal to merus width (Fig. 11a). Posterior margin of basis of peropod 6 and 7 crenulated with small setules; this margin in P 6 slightly concave, in P 7 A B C www.nature.com/scientificreports/ slightly convex. Distal part of these basis articles 1. 5 times wider than ischium width; their distoposterior lobe rectangular (Fig. 12b,c). Lower margin of epimeral plate 2 slightly convex, posterior margin straight or slightly concave; posterodistal tooth of E 2 medium size. Epimeral plate 3 with strongly concave posterior margin and its posterodistal part is distinctly produced (Fig. 11c). Each urosomite segment dorsally with 4 groups of spines and/or setules; two central groups in urosomite 3 are near each other. Lateral groups usually with 2 spines and 1-3 setules, central groups usually with 1 spine and/or 2-3 setules (Fig. 11e). Uropod 3 biramous, endopodite www.nature.com/scientificreports/ length ca 2/3 of exopodite length. Outer margin of the first exopodite article usually with 3 groups of spines accompanied by 1-2 short setae, sometimes longer than spines; rarely on this exopodite margin sub-apically there is fourth group of only setae. Exopodite first article apically with 2-3 spines and several short setae; exopodite second article apically with 2-4 short setae. Inner margin of U 3 exopodite with several groups of setae, several of these setae are feathered. Endopodite apically with 1-3 spines and 3-4 setae; both inner and outer margin of endopodite with several groups of setae, sometimes accompanied with a spine; several setae are feathered (Fig. 11f). Telson lobes apically with 1-3 spines and 2-3 setules. On telson lobes one subbasal and sometimes one subapical spine and/or 1-2 setules (Fig. 11d). Female (Figs. 13, 14, 15): max. length observed 13 mm. Setae on the lower margin of A 2 peduncle articles 4 and 5 longer than in males; the longest setae of the 4th article are bit longer than this article width; longest setae of 5th article are up to 1. 5 × times longer than this article width (Fig. 13b). Posterior margin of pereopod 3 merus richly setose-in 3-4 groups altogether 10-15 setae (Fig. 14a); the longest 1. 5 × longer than merus width. DISTRIBUTION and ECOLOGY: The species has wide distribution in Western Carpathians but also can be found in some locations on Eastern, Southern Carpathians and Ukrainian Lowlands. Species can be found on different altitudes but predominantly in mountain streams with rock bottom and abundant leaf litter substrate. www.nature.com/scientificreports/ Eastern Carpathians). In redescription we are using material from the first location noted as locus typicus in the original description. However we did not have material from the second location (Ukraine) mentioned by S. Karaman, according to our examination of numerous samples collected along the Carpathian Mts, it is highly possible that another MOTU of Gammarus tatrensis may occur there.

Discussion
The Carpathian Arch is recognised as one of the most important extra-Mediterranean European glacial refugia 50 . However, so far the majority of studies suggested that it concerned mostly the Southern and Eastern Carpathians (for a summary see 31 ). Moreover, some studies rejected the existence of northern "cryptic" glacial refugia 51 . In the current study, we provide clear evidence that two species, belonging to the G. balcanicus species complex, survived through the Ice Age in northern refugia within the Carpathians. Moreover, within the most widely distributed species, G. tatrensis, several distinct MOTUs could be distinguished, which survived the Ice Age in separate refugia. The MOTUs, whose evolutionary history reaches the Miocene, are mostly endemic to the northern area of the Carpathians. Only MOTU G2 is widely distributed around the Danube drainage and is characterised by high genetic diversity, the presence of locally endemic haplogroups and private haplotypes. The range of this MOTU extends through the periglacial region. The most striking case of survivalist is MOTU G1, whose present distribution lies within the area covered, in the highest parts, by a glacier during the Last Glacial Maximum. Such a survival of old lineages in the northern Carpathian refugia was also reported for G. fossarum morphospecies 39 , as well as for G. leopoliensis, one of the lineages from within the G. balcanicus complex 20,23 . The existence of glacial refugia in the northern Carpathians is evidenced by Pleistocene fossils of land forest-dwelling molluscs (e.g., Orcula dollium 52 ), which supports the hypothesis of broadleaf forest micro refugia in the region 53 . Organic debris, such as leaves, is known as  (Figs. 13, 14, 15): max. length observed 13 mm. Setae on the lower margin of A 2 peduncle articles 4 and 5 longer than in males; the longest setae of the 4th article are a bit longer than this article width; longest setae of 5th article are up to 1. 5 × times longer than this article width (Fig. 13b). Posterior margin of pereopod 3 merus richly setose-in 3-4 groups altogether 10-15 setae (Fig. 14a); the longest 1. 5 × longer than merus width. In general, we deduce a relatively stable population size through the Pleistocene with a post-Pleistocene increase of the population size (particularly MOTUs G1, G2 and G3) that could be related to postglacial dispersal and colonisation processes. The relatively constant population size, revealed for the lineage G4, is probably related to the fact that it is present in low altitudes and beyond the glacial reach in the Pleistocene. The absence of a decrease in population size and postglacial increment was already observed in the case of other gammarids: G. fossarum 39 and G. jazdzewskii 16 . The proliferation of populations in periglacial regions was also reported in the case of Asellus aquaticus all over Europe 56 . In contrast to G. tatrensis, G. stasiuki has not gone through a noticeable demographic expansion after the Pleistocene and has remained limited geographically. We cannot discard that this pattern is a result of a competitive displacement by G. tatrensis in their common ecological niche. Definition of species borders has been a topic of discussion for decades and still is controversial. In our study, we decide not to tackle the different concepts but to base species hypothesis on the integrative morphological and molecular approach. Results of the study provides evidence for a species new to science and for resurrecting another one (see above for taxonomic discussion) within the Gammarus balcanicus complex in the northern Carpathians. However, our results show that there are more MOTUs that could represent putative cryptic or pseudo-cryptic species. Hopefully, in the future, new methods based on anatomical features, ecology, ethology or secreted pheromones will help in delimitation, identification and definition of such putative species. So far, our study reveals that molecular analysis, primarily based on COI barcoding, is a fast and handy, even if only www.nature.com/scientificreports/ provisional, method of initial species delimitation within such "taxonomically difficult" animal groups. It should be used as a first step of the species delimitation process to propose primary species hypotheses (see also 28 ). Our study is improving an online reference library of DNA barcodes accessible through BOLD 57 . Provision of wellcurated data sets (available within reference libraries), such as the one from the current study, fills gaps in the knowledge of cryptic biota. Such activities are of utter importance in founding the base for future biodiversity assessments (for an overview see 58 ). In recent years, substantial cryptic diversity in numerous taxa, especially among gammarid crustaceans, was discovered using molecular markers (e.g., 22,25,59 ). The urgency for the classification of newly recognised cryptic species remains largely unchallenged while being essential for planning rational and effective conservation strategies 60,61 . On the other hand, we lack a universal and widely acceptable definition of cryptic species, followed by proper delimitation and diagnostic methodology.
The results of our study show that integrative approach, can reveal important patterns in diversity. We demonstrate that the two species, G. tatrensis and G. stasiuki sp. nov., are extremely similar to each other in terms of morphology and there is only one reliable feature that allows to distinguish them. Unfortunately, using SEM to detect possible differences in cuticle ultrastructure does not provide any additional information. Instead, the molecular species delimitation methods show a plethora of MOTUs that may, hypothetically, represent separate www.nature.com/scientificreports/ species. However, while the mtDNA shows categorical and geographically well-structured patterns of diversification, nuclear markers show a less clear picture. Despite that, in all nuclear markers the divergence between haplotypes is noticeable, in the case of EF1-alpha and H3 these haplotypes are shared between mitochondrial MOTUs. In contrast, the third nuclear marker, 28S rDNA, corroborates the delimitation of species and MOTUs. Such discrepancies between gene histories and putative species borders are commonly attributed to gene flow between species, i.e., hybridization and to incomplete lineage sorting-retention of ancient polymorphism 62,63 . Retention of ancestral polymorphisms between species is especially well documented in the case of an adaptive radiation process (review in 64 ). In the light of our findings, it seems that the mechanism may be more frequently found in various taxa 63 . However, proper identification of the drivers of speciation in the studied group will require the use of wide genome sequencing of multiple lineages.
On one hand, sharing of haplotypes is common feature of slowly evolving markers, when related taxa are considered. On the other hand, COI cannot be said to be a universal solution either. In many taxa (e.g., Coleoptera) "COI-delimitation" is widely accepted. This study suggests that strict division of gammarids into species according to COI can lead to overestimation of species number (BIN analysis indicated the presence of up to 21 species only within G. tatrensis). In the case of more intensive use of molecular features in the delimitation and description of new taxa, a broader discussion will certainly be needed to provide generally acceptable rules.

Conclusions
Molecular diversity of the studied G. balcanicus complex revealed the existence of several lineages that survived Pleistocene glaciations in local periglacial refugia in the northern Carpathians. These lineages show a complex pattern of survival and elevated diversification through the Pleistocene, suggesting that the Ice age was a driving force of speciation. We observed two schemes of molecular diversity spatial patterns: firstly, lineages showing low molecular diversity over its geographic range, including the recently colonised areas, and, second, a pattern of numerous locally endemic lineages. All the described lineages can be aggregated in distinguishable MOTUs that represent putative separate species or cryptic and pseudo-cryptic species, whose existence is to be validated via integrative methods. In most cases, these putative cryptic species show postglacial population growth. We revealed a contrasting pattern of mitochondrial vs. nuclear diversification that probably is a result of the preservation of ancient polymorphism in extant lineages of the G. balcanicus complex. Additionally, we redescribed species G. tatrensis, described a new species G. stasiuki sp. nov. and improved barcode reference database by providing new data that can be used in future recognition and biodiversity assessment through molecular methods.

Material and methods
Material collection, identification and analysis. The amphipods were collected by kick-net sampling from 75 sites (BOLD: https:// doi. org/ 10. 5883/ DS-GAMNC ARP) representing the northernmost range of G. balcanicus and initially identified using available keys 8,41 . Individuals of G. balcanicus were selected for the following morphological and molecular analysis. Selected sexually mature individuals of both sexes were dissected, and all appendages except of pleopods were stained with lignin pink (Azophloxin, C18H13N3Na2O8S2) and mounted with Euparal (Carl Roth GmBH, 7356.1) on microscope slides. Afterwards, they were photographed and drawn according to the protocol described by Coleman (2006Coleman ( , 2009. Twelve specimens were used for Scanning Electron Microscopy (SEM). SEM pictures were taken using dried specimens with 10 nm gold coating under Phenom ProX microscope in the Department of Invertebrate Zoology and Hydrobiology of the University of Lodz. Three magnifications were used 5000×, 10,000×, and 30,000×. Pictures were taken from two same-sized individuals from three molecularly distinct populations of G. tatrensis and one of G. stasiuki sp. nov. (see results).
DNA processing and initial analysis. DNA was isolated from a total of 286 individuals. The isolation, amplification and sequencing followed the procedure from 20 . Altogether, five markers were amplified: two mitochondrial including the barcoding fragment of the cytochrome oxidase subunit I (COI, ca. 650 bp) and 16S ribosomal RNA (ca. 350 bp) as well as three nuclear markers: fragments of 28S ribosomal RNA (ca. 900 bp), histone H3 (ca. 300 bp) and elongation factor EF1-alpha (ca 500 bp). The set of primers used in the study is provided in Table S1. Obtained sequences were tested for contamination via BLASTN 65 and were verified as Gammaridae. Sequences of each respective gene fragment were assembled, aligned and trimmed to the same length. Their alignment was performed using MAFFT with automatic settings to select the best algorithm for each data set 66 . All newly generated sequences were deposited in GenBank and BOLD dataset DS-GAMNCARP. The data set was supplemented by the already published sequences from 20,24 . Altogether, our data set consisted of 325 sequences for COI, 80 for 16S, 83 for 28S, 82 for H3 and 72 for EF1-alpha. Heterozygous sites, observed for the EF1-alpha and H3, were coded as ambiguous nucleotides according to IUPAC code. To assess a potential loss in the phylogenetic signal, each single-marker data set was tested for saturation of substitutions with DAMBE 5.3 67 using the index proposed by 68 . In order to visualise distances between sequences, the phylogeny was reconstructed for each single-marker data sets in MEGA X 69 using the Neighbor-Joining (NJ) method 70 and Kimura 2-parameter (K2p) distance 71 with 1,000 replicates 72 . Ambiguous sites and gaps were treated as complete deletions (Fig. S1). For COI, only sequences above 500 bp long were used. Additionally, NJ tree was reconstructed following the same procedure using COI haplotypes.

Species delimitation methods.
To identify the number of molecular operational taxonomic units (MOTUs) that could represent putative cryptic/pseudo-cryptic species within G. balcanicus, we applied eight different methods. Three were distance-based: (1) Barcode Index Number (BIN) System 73 , (2) barcode-gap approach using the Automatic Barcode Gap Discovery (ABGD) software 74 and ASAP: assemble species by automatic partitioning 28 . In ABGD we used primary partitions as a principal for group definition, as they are typically stable on a broader 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 74 . Both in ABGD and ASAP t he standard K2p distance correction was applied. The default values of 0.001 to 0.1 were explored as intraspecific distances and in ABGD gap values from 1 to 1.5 were applied (already tested approach, e.g., 22 ). The remaining five methods were based on phylogeny reconstruction. As a proxy for phylogeny-based methods, Bayesian tree was reconstructed in BEAST 2.5.2 75 81 was performed on the bPTP webserver (available at https:// speci es.h-its. org) with 500,000 iterations of MCMC and 10% burn-in. (6) The multi-rate PTP (mPTP) 82 implements MCMC sampling that provides a fast and comprehensive evaluation of the inferred delimitation. Five runs of 100 M MCMC generations-long chain with burn-in of 10% were performed on the local server. All the mentioned methods are based on a single marker only, COI in this case, therefore we also used the multimarker delimitation method (7) 24 (0.0165 Ma −1 , SD: 0.0018). The material and strategies to calibrate the molecular clock with known geological events were different but resulted in almost identical rates. Additionally, these calibration schemes were cross-validated in the mentioned studies. The ploidy model was set according to the marker type (mtDNA vs nDNA) and the substitution model was selected via bModelTest. Four runs of the MCMC, each 20 M generations long and sampled every 2,000 generations. Results were processed the same way as in case of phylogeny reconstruction for species delimitation. To explore and visualise putative changes of speciation rates through time, and to interpret them in a spatiotemporal context, we have performed two analyses. First, the history of diversification was visualised as a lineage through time (LTT) plot generated in Tracer 1.7.1 from a subset of 1,500 Bayesian chronograms for the species trees generated by *BEAST. The subset of trees was generated using LogCombiner 2.5. Second, we modelled the macroevolutionary dynamics of diversification across the phylogeny with the program Bayesian Analysis of Macroevolutionary Mixtures-BAMM 85 . As input, we used the maximum clade credibility species chronogram generated in *BEAST. First, the priors were preselected using the R package 'BAMMtools' 86 . Four chains of MCMC were used, each 10 M generations long and sampled and chain swap proposed every 1,000 generations. The ESS was checked using R package ' coda' (Plummer et al. 2006) and proved to be > 200. Additionally, visual inspection of MCMC confirmed convergence. Post-run analysis and visualisation were performed using the R package 'BAMMtools' .
Demographic analysis and haplotype network visualisation. The historical demographic patterns were explored using the COI data employing two approaches. First, to test for a recent demographic expansion, Tajima's D, Fu's Fs and Ramos-Onsins and Rozas R2 (Ramos-Onsins & Rozas 2002) indices were calculated using DNAsp6 software 87 . Their statistical significance was evaluated using coalescent simulations with 1,000 replications. Second, the extended Bayesian skyline plot (eBSP) 88 in BEAST 2.5.2 was used to visualise demographic changes through time. The clock model, rate and priors on substitution models for each group were determined in the same way as for the time-calibrated phylogeny, and the population model was set to 0.5. Two MCMC chains were run to ensure convergence for 40 M iterations, sampled every 20,000 iterations, or both values were doubled to provide good ESS values (> 200). One run for each data set was used to plot the eBSP in an R script 80 after a 10% burn-in phase. All demographic analyses were done also for ABGD delimited MOTU's with a sufficient number of samples available. To visualise the haplotype and MOTUs relationships for nuclear markers, haplotype networks were reconstructed. Relationships between haplotypes of 28S (excluding outgroups) were reconstructed using median-joining (MJ) network in POPART 1.7 89 . The homoplasy level parameter (ε) was set at the default value (ε = 0). Relationships for EF1-alpha and H3 nuclear markers, that show the presence of heterozygosity, were reconstructed using the haploweb approach 90 . This method allows showing additional connections between haplotypes found co-occurring in heterozygous individuals. All sequences containing ambiguous sites coded with IUPAC code were phased using the software PHASE 91 according to author guidelines. The full EF1-alpha and H3 data sets, excluding outgroups, were used to generate a network through the Median Joining algorithm using the HaplowebMaker tool 92 .
Ethics approval and consent to participate. All applicable international, national and institutional guidelines for the care and use of animals were followed. All procedures performed in studies involving animals were in accordance with the ethical standards of the institution at which the studies were conducted. www.nature.com/scientificreports/