Fuzzy Species Borders of Glacial Survivalists Revealed Using a Multimarker Approach - A Case Study of Amphipod Crustaceans in the Carpathian Biodiversity Hotspot

doi:10.21203/rs.3.rs-628030/v1

Download PDF

Research

Fuzzy Species Borders of Glacial Survivalists Revealed Using a Multimarker Approach - A Case Study of Amphipod Crustaceans in the Carpathian Biodiversity Hotspot

https://doi.org/10.21203/rs.3.rs-628030/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this older preprint version

Read the latest preprint version →

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 shown 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 G. stasiuki sp. nov. and another we redescribe as G. tatrensis (S. Karaman, 1931).

Animal Science

Behavioral Ecology

Carpathians

key biodiversity hotspots

Alpine orogenesis

Gammarus balcanicus

Pleistocene

morphological examination

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 20th century (e.g., Schäferna 1922; Martynov 1925, 1931; S. Karaman 1931, 1934; Schellenberg 1937, 1943). Subsequently, although 30 years have passed since then, these efforts have been summarized and critically revised in a series of works by G. Karaman and Pinkster (1977a, b, 1987), 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., Jażdżewski 1975a; Jażdżewski & Konopacka 1989; G. Karaman 1989; Iannilli & Ruffo 2002; Alther et al. 2016; Grabowski et al. 2017a; Hupało et al. 2018; Rudolph et al. 2018). 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 (Hou et al. 2011; Hou et al. 2014; Hou & Sket 2016; Sket & Hou 2018; Mamos et al. 2016). Also, recent phylogeographic studies on several geographically widespread Gammarus morphospecies pointed out to, formerly unforeseen, outstanding cryptic and pseudo-cryptic diversity with an ancient divergence between phylogenetic lineages (e.g., Mamos et al. 2014, 2016; Grabowski et al. 2017a,b; Copilaș-Ciocianu et al. 2017, 2018; Hupało et al. 2020; Wattier et al. 2020). 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., Grabowski et al. 2017a; Hupało et al. 2018). However, the latter approach has been disputed for some years. Although it has several shortcomings (for discussion see Meier 2008), it has been appreciated as a tool that helps to overcome the evident global taxonomic impediment (e.g., Coleman 2015; Puillandre et al. 2020). 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 (Puillandre et al. 2020). 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 (Wattier et al. 2020).

The Carpathian Arc is the third-longest and most extensive mountain chain in Europe, stretching for ca. 1,500 km across the central and eastern part of the continent. Being 100–350 km in width, it covers an area of ca. 190,000 km² (Kondracki 1989, Král 2001). This mountain chain extends mainly through Poland, Slovakia, Ukraine and Romania, but also encompasses small areas in Czechia, Austria, Hungary and, disputably, a fragment in Serbia across the Danube valley (Mráz & Ronikier 2016). The Carpathians are known as a major biodiversity hotspot in Europe with a long history and strong tradition of natural history research (Bálint et al. 2011). Nevertheless, modern tools in ecology and evolution, especially those based on molecular markers, were introduced into the study of Carpathian biodiversity only during the last couple of decades, after the disappearance of the local geopolitical constraints (Mráz & Ronikier 2016).

Located north of the Mediterranean Region, already in the Pleistocene permafrost zone and recurrently covered by glaciers in their higher parts, the Carpathians were not perceived, among the major European glacial refugia until recently. However, recently growing evidence that much of the area remained unglaciated during the Last Glacial Maximum (LGM), corroborated by results of palaeoecological and molecular phylogeographic studies, has largely altered such view (for an overview see Bálint et al. 2011; Mráz & Ronikier 2016). Now it is generally accepted, that the Carpathians acted as refugium not only for temperate biota during the LGM but, during warmer episodes, also for the cold-adapted taxa, either at higher elevations (e.g., Schmitt & Varga 2012; Ronikier 2011) or in habitats such as fens and peat bogs at lower elevations (Hájková et al. 2015). Most of the phylogeographic studies in the Carpathians focused on the terrestrial biota; however, the region has also been presumed as an important refugium for aquatic taxa. Already Malicky (1983, 2000), in his “dinodal theory”, predicted that cold-adapted and cold-tolerant freshwater biota would persist in suitable permanent habitats such as springs and headstreams, present in periglacial areas. Subsequent studies have proven this assumption for some aquatic insects, like the dipteran genus Pedicia Latreille, 1809 (Keresztes et al. 2011), riffle beetles (Bozáňová et al. 2020) or several caddisflies (e.g., Bálint et al. 2011).

Most recently, high lineage diversification in the Carpathian springs and streams, indicating long local divergence processes, even in the northernmost parts of the Carpathian Arc (i.e., the Western Carpathians), has been shown for the two gammarid morphospecies widespread in Europe, i.e., Gammarus balcanicus (Mamos et al. 2016; Copilas-Ciocianu & Petrusek 2017) and G. fossarum (Copilas-Ciocianu et al. 2017; Rudolph et al. 2018).

In this study we follow our former large-scale phylogenetic research (Mamos et al. 2016), which has revealed the Miocene phase of the Alpine orogenesis as the main factor influencing lineage divergence of G. balcanicus in the Carpathians. We aim at an in-depth exploration of the spatiotemporal divergence patterns and potential hidden diversity within this morphospecies in headwaters of the northernmost part of the Carpathian Arc, where it has its northern and northwesternmost range limit. In this region, G. balcanicus is the most common and abundant gammarid, particularly at higher altitudes (Grabowski & Mamos 2011).

First, we hypothesise that various lineages survived Pleistocene glaciations in local periglacial refugia in the Carpathians from which, after LGM, some expanded to their present ranges. It should result in a mosaic pattern of lineage distribution with recently colonised areas showing low molecular diversity and the refugial regions with numerous divergent, locally endemic, lineages and high molecular diversity. Second, we hypothesise that some of these divergent lineages may represent cryptic or pseudo-cryptic species that may be defined and described via integrative methods.

1. Material collection, identification and analysis

The amphipods were collected by kick-net sampling from 75 sites (BOLD: DS-GAMNCARP, doi: provided after acceptance) representing the northernmost range of G. balcanicus and initially identified using available keys (Jażdżewski 1975; Karaman & Pinkster 1987). 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 of diagnostic value 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 (2006, 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 5,000x, 10,000x, and 30,000x. 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).

2. DNA processing

DNA was isolated from a total of 286 individuals. The isolation and amplification followed the procedure from Mamos et al. (2016). 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. Subsequently, to check the quality and length of the obtained PCR products, 2 µl aliquot of each reaction mixture was visualised in MidoriGreen-stained (Nippon Genetics) 1.0 % agarose gels. Sequencing of the PCR products was performed using BigDye terminator technology by Macrogen Inc., Korea. To explore the level of heterozygosity, the nuclear markers were sequenced in both directions.

3. Analysis of DNA data

Obtained sequences were tested for contamination via BLASTN (Altschul et al. 1990) and were verified as Gammaridae. Sequences of each respective gene fragment were assembled, aligned and trimmed to the same length in GENEIOUS 6.1.4 (Biomatters 2013). Their alignment was performed using MAFFT with automatic settings to select the best algorithm for each data set (Katoh et al. 2002). All newly generated sequences were deposited in GenBank (accession numbers provided after acceptance) and BOLD dataset DS-GAMNCARP (doi provided after acceptance). The data set was supplemented by the already published sequences from Mamos et al. (2016) and Copilas-Ciocianu & Petrusek (2017). Altogether, our data set consisted of 327 sequences for COI, 78 for 16S, 81 for 28S, 80 for H3 and 70 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 (Xia 2013) using the index proposed by Xia et al. (2003). In order to visualise distances between sequences, the phylogeny was reconstructed for each single-marker data sets in MEGA 7 (Kumar et al. 2016) using the Neighbor-Joining (NJ) method (Saitou & Nei 1987) and Kimura 2-parameter (K2p) distance (Kimura 1980) with 1,000 replicates (Felsenstein 1985). 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.

4. 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 seven different methods. Two were distance-based: (i) Barcode Index Number (BIN) System (Ratnasingham & Hebert 2013) and (ii) barcode-gap approach using the Automatic Barcode Gap Discovery (ABGD) software (Puillandre et al. 2012). The remaining five methods were based on phylogeny reconstruction: general mixed Yule coalescent (GMYC) model-based approach (Pons et al. 2006), using both, (iii) single (sGMYC) and (iv) multiple threshold models (mGMYC); (v) approach using the Bayesian implementation of the Poison Tree Processor (bPTP) (Zhang et al. 2013) and (vi) multi-rate PTP (mPTP) (Kapli et al. 2017). All the mentioned methods are based on a single marker only, COI in this case, therefore we also used the multimarker delimitation method (vii) STACEY (species tree and classification estimation, yarely) (Jones 2017).

4.1 Distance based methods

(i) The BIN method is implemented as part of The Barcode of Life Data Systems (BOLD; Ratnasingham & Hebert 2007). Newly submitted sequences are compared with each other 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 (Barcode Index Number or BIN), already available or newly created if submitted sequences do not cluster with already known BINs. Each BIN is registered in BOLD. (ii) 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 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 (Puillandre et al. 2012). The standard K2p distance correction was applied. The default values of 0.001 to 0.1 were explored as intraspecific distances and gap values from 1 to 1.5 (already tested approach, e.g., Grabowski et al. 2017).

4.2 Phylogeny based methods

As a proxy for phylogeny-based methods, Bayesian tree was reconstructed in BEAST 2.5.2 (Bouckaert et al. 2019). The site model was set up with bModelTest (Bouckaert & Drummond 2017). The tree prior was set to Birth-Death following Bayes factors (> 2 to next model) selection through Path Sampling in BEAST 2.5.2. Four runs of Markov chain Monte Carlo (MCMC) were performed each 20 M generations-long, sampled every 2,000 generations. Runs were examined for convergence in Tracer 1.7 (Rambaut et al. 2018). All runs reached the effective sample size (ESS) above 200 and were combined using LogCombiner 2.5.2. The final tree was summarised with TreeAnnotator 2.5.2, all being part of BEAST 2.5.2 package. The outgroup (G. balcanicus from locus typicus) was removed in all species delimitation analysis.

In 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, two different GMYC approaches, one using the (iii) single threshold and the other one (iv) multiple-threshold model are used to estimate the boundary between intra- and interspecific branching patterns. Both analyses were performed in R software package 'SPLITS' (Species Limits by Threshold Statistics) (Ezard et al. 2009) in R v3.1.0 (R Core Team 2020). The presence of significant differences between the two models was tested using the likelihood ratio test (LRT) in R package 'spiderDEV'. (v) 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). The analysis was performed on the bPTP webserver (available at https://species.h-its.org) with 500,000 iterations of MCMC and 10 % burn-in. (vi) The mPTP method incorporates different levels of intraspecific genetic diversity deriving from differences in either the evolutionary history or sampling of each species. The method 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. Multimarker species delimitation was performed using the (vii) Species Tree And Classification Estimation, Yarely - STACEY v.1.2.1 (Jones 2015), in BEAST 2.5.2 (Bouckaert et al. 2019). All markers were used in the analysis; each COI haplotype not represented by a nuclear marker was removed from the analysis. The minimal clusters were determined by BINs. The nucleotide substitution models were set up with the bModelTest (Bouckaert & Drummond 2017). The Birth-Death model was used to estimate the species tree, other priors were set according to authors guidelines and initially tested [priors: Collapse Height = 0.001, Collapse Weight = 0.5 using a beta prior (1.1), bdc growth rate log-normal (4.6, 2) population prior scale in inverse gamma (2.2), the relative death rate (1.1)]. Ploidy was equal to 2 for nuclear genes and 0.5 for mtDNA genes. Two runs were performed with 100 M generations of MCMC sampled every 10,000 generations. The delimitation results were analysed with speciesDA (Jones 2015) with different settings for the CollapseHeight parameter. Moreover, we tested other species delimitation hypotheses, based on BINs, ABGD, and morphological delimitation (Morphological form A and B) employing multispecies coalescent model: *BEAST package in BEAST 2.5.2. Species trees were constructed utilising each of the mentioned species hypotheses and using the same settings and data set as for the STACEY analysis. The marginal likelihood of each hypothesis was calculated using path sampling in BEAST 2.5.2 with 25 steps and 10 M MCMC generations. Bayes factors were calculated for each hypothesis.

5. Time calibrated reconstruction of phylogeny

For estimating and visualising the temporal framework of G. balcanicus and its sister lineages, we reconstructed the time-calibrated species tree using COI haplotypes and all available markers combined with COI haplotypes in the *BEAST package of BEAST 2.5.2 (Bouckaert et al. 2019). The provisional “species” were defined as a priori according to the BINs provided by BOLD. The molecular clock was calibrated using the COI rate of 0.0166 substitutions/site Ma^− 1 (SD: 0.0022). The value is an average of rates from two independent studies focused on G. balcanicus morphospecies: Mamos et al. (2016) (0.0167 Ma^− 1, SD: 0.0026) and Copilas-Ciocianu & Petrusek (2017) (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 (Bouckaert & Drummond 2017). Four runs of the MCMC, each 20 M generations long and sampled every 2,000 generations, were performed and examined for convergence in Tracer 1.7 (Rambaut et al. 2018). All runs reached the effective sample size (ESS) above 200 and were combined using LogCombiner 2.5.2. The final tree was summarised with TreeAnnotator 2.5.2.

6. Speciation rate through time

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 (Rabosky 2014). As input, we used the maximum clade credibility species chronogram generated in *BEAST. First, the priors were preselected using the R package 'BAMMtools' (Rabosky et al. 2014). 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'.

7. 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 (Tajima 1989), Fu’s Fs (Fu 1997) and Ramos-Onsins and Rozas R2 (Ramos-Onsins & Rozas 2002) indices were calculated using DNAsp6 software (Rozas et al. 2017). Their statistical significance was evaluated using coalescent simulations with 1,000 replications. Second, the extended Bayesian skyline plot (eBSP) (Heled et al. 2008) 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 (http://www.r-project.org) after a 10 % burn-in phase. All demographic analyzes were 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 (Bandelt, Forster, & Röhl 1999) in POPART 1.7 (Leigh & Bryant 2015). 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 (Flot et al. 2010). 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 (Stephens, Smith, & Donnelly 2001) 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 (Spöri & Flot 2020).

Results of morphological and molecular analyses provided support for a redescription of Gammarus tatrensis and description of new species - G. stasiuki sp. nov. Both species belong to the G. balcanicus species complex.

MOTUs delimitation

Seven 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 method divided the data into 18 MOTUs (G1-G18) with P stable for 0.0269–0.0437 (distribution of pairwise distances given in Fig. S2). 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² = 4.654517, df = 6, p = 0.5888235). Both methods supported the hypothesis of multiple species (sGMYC: likelihood ratio = 44.497, p = 2 x 10^− 10, mGMYC: likelihood ratio = 49.152, p = 2x 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 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 Miocene. Soon after, the lineage of G. tatrensis, endemic to Ukrainian Lowlands, diverged from the others. The remaining four lineages started to diverge in Late Miocene/Early Pliocene, while the intra-lineage divergence is dated to Late Pliocene/Pleistocene. G. tatrensis is now the most widely distributed species in the northern Carpathians (Fig. 2b). Its distribution reaches Western/Southern Carpathians, Apuşeni Mts. and Transylvanian Plateau. The divergence of G. stasiuki sp. nov. into three lineages is dated to Late Miocene/Early Pliocene. Its distribution is limited to several locations in the Eastern Carpathians. Except for a few locations in Poland the species was found in northern Romania.

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.

Network reconstructions

The haploweb of EF1-alpha revealed the existence of a complex pattern among the studied lineages (Fig. 4). All lineages belonging to G. tatrensis, with exception of those from the Ukrainian Lowlands (G11, G12), share haplotypes. Only partial distinctiveness can be observed for the lineage G1. The haplotypes of G. stasiuki sp. nov. are, in general, intermingled with G. tatrensis, with partial distinctiveness of haplotypes belonging to BIN ADK0679. A similar pattern is observed in haploweb reconstructed for the H3 nuclear marker, in which almost all MOTUs share haplotypes, including those from the Ukrainian Lowlands. The network of 28S rDNA shows a separation of G. tatrensis from G. stasiuki, the distinctiveness of most of the G. tatrensis MOTUs, and clear divergence of the lineage inhabiting the Ukrainian Lowland from the others within G. tatrensis.

Species description

Gammarus stasiuki sp. nov Jażdżewski K., Mamos T., Grabowski M.

(Taxon name and the article after being published will be registered with Zoobank, type and paratype as well as DNA isolates will be deposited in Museum and Institute of Zoology Polish Academy of Sciences)

Gammarus balcanicus form B; Jażdżewski, 1975, pp. 51, 90–91; Fig. 6L, M, Fig. 8D, E, F, I, K, Fig. 23. Locality – Bieszczady Mts. Gammarus balcanicus form B; Jażdżewski & Konopacka, 1988, p. 76; Fig. 1D, E, F, Fig. 2. Locality – Bieszczady Mts. G. balcanicus form B; Zieliński, 1995, pp. 61, 64, 69; Fig. 1. Locality – Dwernik = Prowcza stream, Bieszczady Mts., 1020 m a.s.l. G. balcanicus form B; Zieliński, 1998, p. 131, Fig. 1, p. 137, Fig. 6. Locality – four streams in south-easternmost Poland. Gammarus balcanicus form B; Konopacka, Jażdżewski & Jędryczkowski, 2001, pp. 36–39; Fig. 2. Eleven localities in Bieszczady National Park.

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.

Material examined

The morphologically examined material consists of 4 samples collected from Eastern Carpathians (Poland and Romania). Type sample was collected in the locality Mała Rawka Mt., 1140 m a.s.l., on 23.04.2016, by M. Grabowski and T. Mamos.

Description

The morphological analysis supported discrimnation of G. stasiuki sp. nov. and G. tatrensis based on setation of the second antenna (Fig. 1, Fig. 2e). The examination of ultrastructure through SEM did not show any features delimiting the taxa or MOTUs (Fig. S4).

Male (Fig. 6, 7, 8): max. length observed 13 mm. A I of the length of head and 4 pereon segments, A II somewhat shorter – head + 3 pereon segments. Head lateral lobe rounded, eyes medium size, oval or reniform; eye length +/- equal to the A I basal width (Fig. 5a, b). A I main flagellum with 20–25 articles and accessory flagellum with 3–4 articles (Fig. 5b). A II flagellum with 10–12 articles; in adult males (over 10 mm) with 4–5 calceoli on 4 to 8 flagellum articles, often poorly visible (Fig. 5c). Third article of mandibular palp with a brush of over 20 D-setae forming even row (Fig. 5 ).

Lower margins of basal A II articles 4 and 5 richly setose, with 5–6 groups of 2–4 setae. The length of longest setae +/- equal to the width of 4th article and somewhat longer than the width of 5th article (Fig. 5b). Hind margins of gnathopod’s carpus and propus richly setose (Fig. 5D, E). P III merus hind margin with 4–5 groups of numerous setae as long or a bit longer than merus width (Fig. 5a). Basis of P VI and P VII with hind 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 VI and P VII +/- rectangular. Hind margin of P VI basis somewhat concave, hind margin of P VII basis slightly convex (Fig. 5e,f).

Second epimeral plate 2 (E2) rather characteristic, with lower margin distinctly convex ending with a small tooth; hind 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 III biramous, endopodite length +/- 2/3 of exopodite length (Fig. 6c). Outer margin of the first article of U III 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 III apically with 1–2 spines and 3–4 setae, some over 2 times longer than spines. Inner margin of U III exopodite with several groups of 1–3 setae, ca. half of these setae are feathered. Outer margin of U III endopodite with 3–5 groups of spines and setae, of which 3–4 are feathered. Inner margin of U III 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 (Fig. 8, 9): max. length observed 11 mm. Clear sexual dimorphism in the setation of appendages: setae on the lower margin of A II peduncular articles 4 and 5 distinctly longer than in males; in 4th article setae are 1,5x longer than this article width and in 5th article two times longer than this article width (Fig. 8c). Similarly, in females, P III merus is more setose, with setae twice as long as the merus width (Fig. 8g). In the U III exopodite outer margin longest setae are 2 times longer than spines. Also, apical setae of U III 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.

Gammarus tatrensis (S. Karaman, 1931) - redescription

(Taxon name and the article after being published will be registered with Zoobank, neotype and paratype as well as DNA isolates will be deposited in Museum and Institute of Zoology Polish Academy of Sciences)

Rivulogammarus tatrensis; S. Karaman 1931, pp. 97–98, Figs. 4a, b; localities: 1) Ďumbier Mt. (Low Tatra Mts, Slovakia, that time Czechoslovakia); 2) Kuzy, Carpathian Mts, western Ukraine (probably not G. tatrensis)

Gammarus (Rivulogammarus) balcanicus tatrensis (S. Karaman) 1931 (sic!); Straskraba 1953, pp. 218–222, Figs. 4A-D, locality – Tatra Mts, several streams in eastern Slovakia (that time Czechoslovakia)

Gammarus (Rivulogammarus) balcanicus subsp. tatrensis; Straskraba 1957, p.256, several streams in western Ukraine, Carpathian Mts ( probably not G. tatrensis)

Gammarus (Rivulogammarus) balcanicus tatrensis S. Karaman 1931 (sic!); Micherdziński 1959, pp. 562–567, Tab. LXXVI Fig. 3, Tab.LXXVII Figs. 3–4.

Gammarus balcanicus tatrensis (S. Karaman, 1931); Straskraba 1962, Fig. 2; numerous localities in northern Carpathian Mts in Slovakia (that time Czechoslovakia), i.a. Vyšná Boca near Ďumbier Mt, 7.07.1954.

G. b. tatrensis (S. Kar.) and G. balcanicus form A; Jażdżewski 1975b, Fig. 8a,b; numerous localities in Polish Carpathian Mts and sub-Carpathian region.

Gammarus balcanicus tatrensis and G. balcanicus form A; Konopacka and Jażdżewski 1985, pp. 375–376 (probably not G. tatrensis)

Material examined

The morphologically examined material consists of 54 samples collected all along the Carpathian Arch and Ukrainian Lowlands. Type sample was collected in the locality Vysna Boca at the Ďumbier Mt, 1030 m a.s.l., on 15.05.2015, by all the authors.

Redescription

Male (Fig. 10, 11, 12): max. length observed 14 mm. Head lateral lobe rounded, eyes oval or reniform (Fig. 10a); eye length equal to the A I basal width. A I of the length of head and over 4 pereon segments, A II of the length of head and over 3 pereon segments. A I flagellum with articles 20–30, accessory flagellum with 3–4 articles. Flagellum of A II with 10–12 articles. Calceoli set on basal articles are present only in males larger than 9 mm. Mandibular palp with a brush of over 20 D-setae in even row. Lower margin of A II basal articles (4 and 5) poorly setose: 4th article with 2–3 groups of short setae (ca. 0.5 of the article width) and 5th article with 3–4 groups of setae, their length equal to the article width (Fig. 10c ). Hind margin of merus of pereopod III with 3–4 groups of setae, the longest is equal to merus width (Fig. 11a). Hind margin of basis of peropod VI and VII crenulated with small setules; this margin in P VI slightly concave, in P VII slightly convex. Distal part of these basis articles 1,5 times wider than ischium width; their distaloposterior lobe rectangular (Fig. 12b,c). Lower margin of epimeral plate 2 slightly convex, hind margin straight or slightly concave; posterodistal tooth of E 2 medium size. Epimeral plate 3 with strongly concave hind 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 III 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 III biramous, endopodite 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 a 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 III 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).

Taxonomical comment

S. Karaman (1931) has described Rivulogammarus tatrensis based on the samples collected at the Ďumbier Mt. (Low Tatra Mts, Slovakia, Czechoslovakia at that time) and in the locality Kuziy Massif (western-most Ukraine - Eastern Carpathians). In this second location, according to our recent examination of numerous samples collected along the whole northern arch of Carpathian Mts, the true Gammarus tatrensis s.s. would not occur. More probably, this was another member of the Gammarus balcanicus complex. Similarly, the records of Gammarus balcanicus tatrensis (Straškraba 1953, 1957, 1962, Micherdziński 1959) from numerous localities of the northern Carpathian Mts in Slovakia and Poland (except the Tatra Mts) had to do with two different species of G. balcanicus complex.

Survival through the Ice Age

The Carpathian Arch is recognised as one of the most important extra-Mediterranean European glacial refugia (Provan & Bennett 2008). However, so far the majority of studies suggested that it concerned mostly the Southern and Eastern Carpathians (for a summary see Bálint et al. 2011). Moreover, some studies rejected the existence of northern “cryptic” glacial refugia (Tzedakis et al. 2013). 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 (Copilas-Ciocianu et al. 2017), as well as for G. leopoliensis, one of the lineages from within the G. balcanicus complex (Mamos et al. 2016; Copilas-Ciocianu et al. 2018). The existence of glacial refugia in the northern Carpathians is evidenced by Pleistocene fossils of land forest-dwelling molluscs (e.g., Orcula dollium, Harl et al. 2014), which supports the hypothesis of broadleaf forest micro-refugia in the region (Juřičková et al. 2014). Organic debris, such as leaves, is known as the primary diet component for amphipods of the genus Gammarus (Väinölä et al. 2008) and could provide a food source even through the Last Glacial Maximum. Moreover, compared to other European species, the G. balcanicus complex consists of amphipods that can be usually found in higher altitudes, especially in the northern regions of the Carpathians (Grabowski & Mamos 2011), therefore being adapted to well-oxygenated and low-temperature habitats. Thanks to such adaptation, survival in cold, northern refugia could be promoted for the species complex.

The increase of diversification rate of G. balcanicus complex in the northern Carpathians through Pleistocene suggests that Ice Age could be a driving force of local lineages speciation. Such proliferation of their diversity could be possible due to climate oscillation and survival of different lineages in multiple cryptic refuges in the periglacial zone during adverse conditions.. In fact, large areas of the northern range of the Carpathians was a periglacial area, not covered permanently by ice sheet. It could thus provide suitable habitats for Gammarus sp. survival (Zasadni & Kłapyta 2014 and references therein). It also seems that glaciations were the main factor behind demographic changes of the local lineages. These changes can be visible in both studied species and their MOTUs. For example, in the Western Carpathians, MOTU G1 of G. tatrensis and G. stasiuki sp. nov. experienced a putative population bottlenecks at the end of the Pleistocene (G. stasiuki sp. nov. only according to eBSP). These lineages are present most northerly and, therefore, probably more affected by the Pleistocene glaciations. According to neutrality tests, also MOTU G4 could putatively go through the bottleneck effect. However, this is not supported by the eBSP plot. On the contrary, distributed only in the Eastern Carpathians, MOTU G3, does not have any signs of the negative impact of glaciation on its population size. In general, we can observe 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 also observed in the case of G. fossarum (Copilas-Ciocianu et al. 2017) and G. jazdzewskii (Rudolph et al. 2018). The proliferation of populations in periglacial regions was also observed in the case of Asellus aquaticus all over Europe (Sworobowicz et al. 2020). 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 exclude that this is a result of losing competition with G. tatrensis for the same ecological niche.

The reference library

Our study, based on the integrative approach, 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 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 that is to propose primary species hypotheses (see also Puillandre et al. 2020). Our study is improving an online reference library of DNA barcodes accessible through BOLD (Ratnasingham & Hebert 2007). Provision of well-curated data sets (available within reference libraries), such as the one from the current study, fills gaps in knowledge of cryptic biota. Such activities are of utter importance in founding bases for future biodiversity assessments (for an overview see Weigand et al. 2019).

Species borders

In recent years, substantial cryptic diversity in numerous taxa, especially among gammarid crustaceans, was discovered using molecular markers (e.g., Katouzian et al. 2016; Grabowski et al. 2017; Wattier et al. 2020). The urgency for the classification of newly recognised cryptic species remains largely unchallenged while being essential for planning rational and effective conservation strategies (Bickford et al. 2007; Delić et al. 2017). On the other hand, we still miss a universal and widely acceptable definition of cryptic species, followed by proper delimitation and diagnostic methodology.

The results of our study show that using an integrative approach, surprising patterns in diversity can be revealed. 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 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 can be attributed to gene flow between species, i.e., hybridization and to incomplete lineage sorting - retention of ancient polymorphism (Maddison 1997; Nosil 2008). Retention of ancestral polymorphisms between species is especially well documented in the case of an adaptive radiation process (review in Berner & Salzburger 2015). In the light of our findings, it seems that the mechanism is more frequently found in various taxa (Nosil 2008). However, proper identification of the drivers of speciation in the studied group will require the use of wide genome sequencing of multiple lineages.

The sharing of haplotypes in slowly evolving markers is logical. 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 in order to find a generally acceptable rules.

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: one, when a lineage shows low molecular diversity over its geographic range, including the recently colonised areas, and, second, when numerous locally endemic lineages can be found. 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 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.

Acknowledgments:

We would like to thank all the people that helped during field sampling: Aleksandra Jabłońska, Misel Jelic, Jana Bozanova, Andrzej Zawal, Aleksandra Bańkowska, Przemysław Śmietana, Agnieszka Szlauer-Łukaszewska.

Founding:

The study was funded by the Polish National Science Centre grant Miniatura 2017/01/X/NZ8/01607 as well as by the statutory funds of the University of Łódź. Additional funding for exchange and field work was provided by the Slovak R&D Agency (APVV) bilateral project Slovakia – Poland 2015 SK-PL-2015-0042 and Polish Ministry of Science and Higher Education. Tomasz Mamos was supported by the Scholarship of the Polish National Agency for Academic Exchange (NAWA) Bekker Programme (project nb. PN/BEK/2018/1/00225).

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.

The material was collected in accordance with the permits issued by The District Office, Department of Environmental Care Trenčín (OU-TN-OSZP1-2015/001937-12/Du) and the Ministry of the Environment of the Slovak Republic (5198/2015-2.3, 3735/2015-2.3).

Consent for publication:

Not applicable.

Competing interests:

The authors declare that they have no competing interests.

Availability of data and materials

GenBank accession numbers provided after acceptance. All metadata and sequences are stored in BOLD Dataset DS-GAMNCARP (doi will be provided after acceptance of the MS). Taxon name and the article after being published will be registered with Zoobank, type and paratype as well as DNA isolates will be deposited in Museum and Institute of Zoology Polish Academy of Sciences.

Alther R, Fišer C, Altermatt F. Description of a widely distributed but overlooked amphipod species in the European Alps. Zoological Journal of the Linnean Society. 2016. https://doi.org/10.1111/zoj.12477.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic Local Alignment Search Tool. J Mol Biol. 1990;215:403–10.
Bálint M, Ujvarosi L, Theissinger K, Lehrian S, Meszaros N, Pauls SU. The Carpathians as a Major Diversity Hotspot in Europe. Biodiversity Hotspots. Distribution and Protection of Conservation Priority Areas. Heidelberg: Springer; 2011. pp. 189–205.
Bickford D, Lohman DJ, Sodhi NS, Ng PK, Meier R, Winker K, Ingram KK, Das I. Cryptic species as a window on diversity and conservation. Trends Ecol Evol. 2007;22:148–55.
Bouckaert R, Vaughan TG, Barido-Sottani J, Duchêne S, Fourment M, Gavryushkina A, Heled J, Jones G, Kühnert D, De Maio N, Matschiner M, Mendes FK, Müller NF, Ogilvie HA, du Plessis L, Popinga A, Rambaut A, Rasmussen D, Siveroni I, Suchard MA, Wu C-H, Xie D, Zhang C, Stadler T, Drummond AJ. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. Plos Comput Biol. 2019;15:e1006650.
Bouckaert RR, Drummond AJ. bModelTest: Bayesian phylogenetic site model averaging and model comparison. BMC evolutionary biology. 2017;17:42.
Coleman CO. Taxonomy in Times of the Taxonomic Impediment – Examples from the Community of Experts on Amphipod Crustaceans. J Crustacean Biol. 2015;35:729–40.
Copilaş-Ciocianu D, Petrusek A. Phylogeography of a freshwater crustacean species complex reflects a long-gone archipelago. J Biogeogr. 2017;44:421–32.
Copilas-Ciocianu D, Rutova T, Paril P, Petrusek A. Epigean gammarids survived millions of years of severe climatic fluctuations in high latitude refugia throughout the Western Carpathians. Mol Phylogenet Evol. 2017;112:218–29.
Copilaş-Ciocianu D, Zimţa A-A, Grabowski M, Petrusek A. Survival in northern microrefugia in an endemic Carpathian gammarid (Crustacea: Amphipoda). Zoolog Scr. 2018;47:357–72.
Ezard T, Fujisawa T, Barraclough TG. 2009. SPLITS: SPecies' LImits by Threshold Statistics. R package version 1.0–18/r45 Available from: URL http://R-Forge.R-project.org/projects/splits/.
Felsenstein J. Confidence limits on phylogenies: an approach using the bootstrap. Evolution. 1985;39:783–91.
Flot JF, Couloux A, Tillier S. 2010. Haplowebs as a graphical tool for delimiting species: a revival of Doyle's "field for recombination" approach and its application to the coral genus Pocillopora in Clipperton. BMC evolutionary biology 10.
Fu YX. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147:915–25.
Grabowski M, Mamos T. Contact Zones, Range Boundaries, and Vertical Distribution of Three Epigean Gammarids (Amphipoda) in the Sudeten and Carpathian Mountains (Poland). Crustaceana. 2011;84:153–68.
Grabowski M, Mamos T, Bącela-Spychalska K, Rewicz T, Wattier RA. Neogene paleogeography provides context for understanding the origin and spatial distribution of cryptic diversity in a widespread Balkan freshwater amphipod. PeerJ. 2017;5:e3016.
Grabowski M, Wysocka A, Mamos T. Molecular species delimitation methods provide new insight into taxonomy of the endemic gammarid species flock from the ancient Lake Ohrid. Zoological Journal of the Linnean Society. 2017;20:1–14.
Hájková P, Horsák M, Hájek M, Jankovská V, Jamrichová E, Moutelíková J. Using multi-proxy palaeoecology to test a relict status of refugial populations of calcareous-fen species in the Western Carpathians. The Holocene. 2015;25:702–15.
Harl, J., Duda, M., Kruckenhauser, L., Sattmann, H., Haring, E. 2014. In Search of Glacial Refuges of the Land Snail < italic > Orcula dolium</italic> (Pulmonata, Orculidae) - An Integrative Approach Using DNA Sequence and Fossil Data. PloS one 9, e96012.
Heled J, Drummond A. Bayesian inference of population size history from multiple loci. BMC evolutionary biology. 2008;8:289.
Hou Z, Sket B. A review of Gammaridae (Crustacea: Amphipoda): the family extent, its evolutionary history, and taxonomic redefinition of genera. Zoological Journal of the Linnean Society. 2016;176:323–48.
Hou Z, Sket B, Fiser C, Li S. Eocene habitat shift from saline to freshwater promoted Tethyan amphipod diversification. Proc Natl Acad Sci USA. 2011;108:14533–8.
Hou Z, Sket B, Li S. 2014. Phylogenetic analyses of Gammaridae crustacean reveal different diversification patterns among sister lineages in the Tethyan region. Cladistics 30.
Hupało K, Karaouzas I, Mamos T, Grabowski M. Molecular data suggest multiple origins and diversification times of freshwater gammarids on the Aegean archipelago. Sci Rep-Uk. 2020;10:19813.
Hupalo K, Mamos T, Wrzesinska W, Grabowski M. First endemic freshwater Gammarus from Crete and its evolutionary history-an integrative taxonomy approach. PeerJ. 2018;6:e4457.
Iannilli V, Ruffo S. Apennine and Sardinian species of Gammarus, with the description of Gammarus elvirae n. sp. (Crustacea Amphipoda, Gammaridae). Boll Acc Gioenia Sci Nat. 2002;35:519–32.
Jażdżewski K. 1975 a. Remarks on Gammarus lacustris G.O. Sars,1863, with description of Gammarus varsoviensis n. sp. Bijdragen tot de Dierkunde 45, 71–86.
Jażdżewski K. 1975 b. Morfologia, taksonomia i występowanie w Polsce kiełży z rodzajów Gammarus Fabr. i Chaetogammarus Mart. (Crustacea, Amphipoda). Acta Universitatis Lodziensis: 185 pp.
Jażdżewski K, Konopacka A. 1988. Notes on the Gammaridean Amphipoda of the Dniester River Basin and Eastern Carpathians. Crustaceana. Supplement, 72–89.
Jażdżewski K, Konopacka A. 1989. Gammarus leopoliensis nov. sp. (Crustacea, Amphipoda) from Eastern Carpathians. Bulletin Zoölogisch Museum, University of Amsterdam 11, 185–196.
Jones G. Algorithmic improvements to species delimitation and phylogeny estimation under the multispecies coalescent. J Math Biol. 2017;74:447–67.
Juřičková L, Horáčková J, Ložek V. Direct evidence of central European forest refugia during the last glacial period based on mollusc fossils. Quatern Res. 2014;82:222–8.
Kapli P, Lutteropp S, Zhang J, Kobert K, Pavlidis P, Stamatakis A, Flouri T. Multi-rate Poisson tree processes for single-locus species delimitation under maximum likelihood and Markov chain Monte Carlo. Bioinformatics. 2017;33:1630–8.
Karaman G, Pinkster S. Freshwater Gammarus species from Europe, North Africa and adjacent regions of Asia (Crustacea Amphipoda). Part II. Gammarus roeseli-group and related species. Bijdragen tot de Dierkunde. 1977;47:165–96.
Karaman G, Pinkster S. Freshwater Gammarus species from Europe, North Africa and adjacent regions of Asia (Crustacea Amphipoda). Part I. Gammarus pulex-group and related species. Bijdragen tot de Dierkunde. 1977;47:1–97.
Karaman G, Pinkster S. Freshwater Gammarus species from Europe, North Africa and adjacent regions of Asia (Crustacea-Amphipoda). Part III. Gammarus balcanicus-group and related species. Bijdragen tot de Dierkunde. 1987;57:207–60.
Karaman GS. 1989. New species of the family Gammaridae from Ohrid Lake basin, Gammarus sketi, n. sp., with emphasis on the subterranean members of genus Gammarus Fabr. (Contribution to the knowledge of the Amphipoda 191.). Glasnik Odjeljenja prirodnih nauka, Crnogorska akademija nauka i umjetnosti 7, 53–71.
Karaman SL. 1931. Beitrag zur Kenntnus der Susswasseramphiopden. Bull. Soc. Scien. Skoplje IX, 93–107.
Katoh K, Standley DM. 2013. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Molecular biology and evolution 30.
Katouzian A-R, Sari A, Macher JN, Weiss M, Saboori A, Leese F, Weigand AM. Drastic underestimation of amphipod biodiversity in the endangered Irano-Anatolian and Caucasus biodiversity hotspots. Sci Rep-Uk. 2016;6:22507.
Keresztes L, Kolcsár L-P, Török E, Dénes A-L. The spring dwelling dipteran genus Pedicia Latreille in the Carpathian area: diversity, divergence and distribution – case studies. In: Keresztes L, Markó B, editors. The Carpathians as speciation centres and barriers: from case studies to general patterns. Cluj University Press; 2011. p. 168.
Kimura M. 1980. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. Journal of molecular evolution 16.
Kondracki J. 1989. Karpaty. WSiP, Warszawa.
Konopacka A, Jażdżewski K. Stream ecosystems in mountain grassland (West Carpathians). 13. Gammarid species. Acta Hydrobiologica. 1985;27:369–78.
Konopacka A, Jażdżewski K, Jędryczkowski W. 3. Skorupiaki (Crustacea). 3.1. Pancerzowce (Malacostraca) Bieszczadów. In: Razowski J, editor. Monografie Bieszczadzkie t. VII.; 2001. pp. 35–47.
Král V. Fyzicka geografie Evropy. Praha: Academia; 2001.
Kumar S, Stecher G, Tamura K. 2016. MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. Molecular biology evolution 33, 1870–4.
Maddison WP. Gene Trees in Species Trees. Syst Biol. 1997;46:523–36.
Malicky H. Chorological patterns and biome types of European Trichoptera and other freshwater insects. Archiv für Hydrobiol. 1983;96:223–44.
Malicky H. Arealdynamik und Biomgrundtypen am Beispiel der Köcherfliegen (Trichoptera). Entomologica Basiliensia. 2000;22:235–59.
Mamos T, Wattier R, Burzyński A, Grabowski M. The legacy of a vanished sea: a high level of diversification within a European freshwater amphipod species complex driven by 15 My of Paratethys regression. Molecular ecology. 2016;25:795–810.
Mamos T, Wattier R, Majda A, Sket B, Grabowski M. Morphological vs. molecular delineation of taxa across montane regions in Europe: the case study of Gammarus balcanicus Schäferna, 1922 (Crustacea: Amphipoda). Journal of Zoological Systematics Evolutionary Research. 2014;52:237–48.
Meier R, Wheeler QD. DNA sequences in taxonomy - Opportunities and challenges. In: Wheeler QD, editor. The New Taxonomy. Boca Raton: CRC Press; 2008. p. 256.
Micherdziński W. Kiełże rodzaju Gammarus Fabricius (Amphipoda) w wodach Polski. Acta Zoologica Cracoviensia. 1959;4:527–637.
Mráz P, Ronikier M. Biogeography of the Carpathians: evolutionary and spatial facets of biodiversity. Biol J Lin Soc. 2016;119:528–59.
Nosil P. Speciation with gene flow could be common. Molecular ecology. 2008;17:2103–6.
Plummer M, Best N, Cowles K, Vines K. CODA: Convergence Diagnosis and Output Analysis for MCMC. R News. 2006;6:7–11.
Pons J, Barraclough TG, Gomez-Zurita J, Cardoso A, Duran DP, Hazell S, Kamoun S, Sumlin WD, Vogler AP. Sequence-Based Species Delimitation for the DNA Taxonomy of Undescribed Insects. Syst Biol. 2006;55:595–609.
Provan J, Bennett KD. Phylogeographic insights into cryptic glacial refugia. Trends Ecol Evol. 2008;23:564–71.
Puillandre N, Brouillet S, Achaz G. ASAP: assemble species by automatic partitioning. Molecular ecology resources. 2021;21:609–20.
Puillandre N, Lambert A, Brouillet S, Achaz G. ABGD, Automatic Barcode Gap Discovery for primary species delimitation. Molecular ecology. 2012;21:1864–77.
Rabosky DL. Automatic Detection of Key Innovations, Rate Shifts, and Diversity-Dependence on Phylogenetic Trees. PloS one. 2014;9:e89543.
Rabosky DL, Grundler M, Anderson C, Title P, Shi JJ, Brown JW, Huang H, Larson JG. BAMMtools: an R package for the analysis of evolutionary dynamics on phylogenetic trees. Methods Ecol Evol. 2014;5:701–7.
Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior Summarization in Bayesian Phylogenetics Using Tracer 1.7. Syst Biol. 2018;67:901–4.
Ramos-Onsins SE, Rozas J. 2002. Statistical Properties of New Neutrality Tests Against Population Growth. Molecular biology and evolution 19, 2092–2100.
Ratnasingham S, Hebert PD. A DNA-based registry for all animal species: the barcode index number (BIN) system. PloS one. 2013;8:e66213.
Ronikier M. Biogeography of high-mountain plants in the Carpathians: An emerging phylogeographical perspective. Taxon. 2011;60:373–89.
Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, Sánchez-Gracia A. 2017. DnaSP 6: DNA Sequence Polymorphism Analysis of Large Data Sets. Molecular biology and evolution 34, 3299–3302.
Rudolph K, Coleman CO, Mamos T, Grabowski M. Description and post-glacial demography of Gammarus jazdzewskii sp. nov. (Crustacea: Amphipoda) from Central Europe. Syst Biodivers. 2018;16:587–603.
Saitou N, Nei M. The neighbor-joining method: a new method for reconstructing phylogenetic trees. Molecular biology evolution. 1987;4:406–25.
Schmitt T, Varga Z. 2012. Extra-Mediterranean refugia: The rule and not the exception? Frontiers in zoology 9.
Schäferna K. Amphipoda balcanica, spolu s pozn_amkami o jin_ ych sladkovodn_ıch Amphipodech. Mem Soc R Sci Boheme Prague. 1922;12:1–111.
Sket B, Hou Z. 2018. Family Gammaridae (Crustacea: Amphipoda), mainly its Echinogammarus clade in SW Europe. Further elucidation of its phylogeny and taxonomy, 61.
Spöri Y, Flot J-F. HaplowebMaker and CoMa: Two web tools to delimit species using haplowebs and conspecificity matrices. Methods Ecol Evol. 2020;11:1434–8.
Stephens M, Smith NJ, Donnelly P. A New Statistical Method for Haplotype Reconstruction from Population Data. The American Journal of Human Genetics. 2001;68:978–89.
Straškraba M. Předběžná zpráva o rozšíření rodu Gammarus v ČSR. Věstník Československé Společnosti Zoologické. 1953;17:212–27.
Straškraba M. 1957. Beitrag zur Kenntnis der Amphipodenfauna Karpatenrusslands (USSR). Věstník Československé Společnosti Zoologické, 21, 256–272.
Straškraba M. 1962. Amphipoden der Tschechoslovakei nach den Sammlungen von. Prof. Hrabě. I. Věstník Československé Společnosti Zoologické, 26, 117–145.
Sworobowicz L, Mamos T, Grabowski M, Wysocka A. Lasting through the ice age: The role of the proglacial refugia in the maintenance of genetic diversity, population growth, and high dispersal rate in a widespread freshwater crustacean. Freshwater Biology n/a; 2020.
Tajima F. Statistical-Method for Testing the Neutral Mutation Hypothesis by DNA Polymorphism. Genetics. 1989;123:585–95.
Theissinger K, Bálint M, Feldheim KA, Haase P, Johannesen J, Laube I, Pauls SU. Glacial survival and post-glacial recolonization of an arctic–alpine freshwater insect (Arcynopteryx dichroa, Plecoptera, Perlodidae) in Europe. J Biogeogr. 2013;40:236–48.
Tzedakis PC, Emerson BC, Hewitt GM. Cryptic or mystic? Glacial tree refugia in northern Europe. Trends Ecol Evol. 2013;28:696–704.
Wattier R, Mamos T, Copilaş-Ciocianu D, Jelić M, Ollivier A, Chaumot A, Danger M, Felten V, Piscart C, Žganec K, Rewicz T, Wysocka A, Rigaud T, Grabowski M. Continental-scale patterns of hyper-cryptic diversity within the freshwater model taxon Gammarus fossarum (Crustacea, Amphipoda). Sci Rep-Uk. 2020;10:16536.
Weigand H, Beermann AJ, Čiampor F, Costa FO, Csabai Z, Duarte S, Geiger MF, Grabowski M, Rimet F, Rulik B, Strand M, Szucsich N, Weigand AM, Willassen E, Wyler SA, Bouchez A, Borja A, Čiamporová-Zaťovičová Z, Ferreira S, Dijkstra KB, Eisendle U, Freyhof J, Gadawski P, Graf W, Haegerbaeumer A, van der Hoorn BB, Japoshvili B, Keresztes L, Keskin E, Leese F, Macher JN, Mamos T, Paz G, Pesic V, Pfannkuchen DM, Pfannkuchen MA, Price BW, Rinkevich B, Teixeira MAL, Varbiro G, Ekrem T. 2019. DNA barcode reference libraries for the monitoring of aquatic biota in Europe: Gap-analysis and recommendations for future work. The Science of the total environment 678, 499–524.
Xia X, Xie Z. DAMBE: Software package for data analysis in molecular biology and evolution. J Hered. 2001;92:371–3.
Xia X, Xie Z, Salemi M, Chen L, Wang Y. An index of substitution saturation and its application. Mol Phylogenet Evol. 2003;26:1–7.
Zasadni J, Kłapyta P. The Tatra Mountains during the Last Glacial Maximum. Journal of Maps. 2014;10:440–56.
Zhang J, Kapli P, Pavlidis P, Stamatakis A. A general species delimitation method with applications to phylogenetic placements. Bioinformatics. 2013;29:2869–76.
Zieliński D. Life History of Gammarus balcanicus Schäferna, 1922 from the Bieszczady Mountains (Eastern Carpathians, Poland). Crustaceana. 1995;68(1):61–72.
Zieliński D. 1998. Life Cycle and Altitude Range of Gammarus leopoliensis Jażdżewski & Konopacka, 1989 (Amphipoda) in South-Eastern Poland. Crustaceana 71.

Fig.S2histogramBarcodeGap.pdf
Fig. S2 The gap analysis. Distribution of pairwise k2distances obtained from ABGD.
Fig.S3LTT.jpg
Fig. S3 Lineages through time plot obtained from Bayesian analysis (see Material and Methods)
Fig.S1.pdf
Fig. S1 Neighbour-Joining distance trees based on sequence data from (a) COI, (b) 16S, (c) 28S, (d) H3, (e) EF1-alpha. Bootstrap values are provided at the respective nodes. Scale bars indicate the number of substitutions per site according to the Kimura 2-parameter distance model. Whenever possible, specimens names are followed by BOLD BINs (see DS-GAMNCARP).
Fig.S4.pdf
Fig. S4 Images of ultrastructure through SEM showing epimeral plate 2 of selected MOTUs. Magnification in each column from the left: 3000x, 10000x, 30000x.
Supportingtables.pdf

Download PDF

Version 1

posted

You are reading this older preprint version

Read the latest preprint version →

Fuzzy Species Borders of Glacial Survivalists Revealed Using a Multimarker Approach - A Case Study of Amphipod Crustaceans in the Carpathian Biodiversity Hotspot

Archived Versions:

Version 1

Abstract

Figures

Introduction

Material And Methods

1. Material collection, identification and analysis

2. DNA processing

3. Analysis of DNA data

4. Species delimitation methods

4.1 Distance based methods

4.2 Phylogeny based methods

5. Time calibrated reconstruction of phylogeny

6. Speciation rate through time

7. Demographic analysis and haplotype network visualisation

Results

MOTUs delimitation

Time calibrated reconstruction of phylogeny and distribution

Speciation rate changes through time

Demographic analysis

Network reconstructions

Discussion

Conclusions

Declarations

References

Supplementary Files

Archived Versions:

Version 1