Phylogeography of the Chydorus sphaericus Group (Cladocera: Chydoridae) in the Northern Palearctic

The biodiversity and the biogeography are still poorly understood for freshwater invertebrates. The crustacean Chydorus sphaericus-brevilabris complex (Cladocera: Chydoridae) is composed of species that are important components of Holarctic freshwater food webs. Recent morphological and genetic study of the complex has indicated a substantial species diversity in the northern hemisphere. However, we know little of the geographic boundaries of these novel lineages. Moreover, a large section of the Palearctic remains unexamined at the genetic level. Here we attempt to address the biodiversity knowledge gap for the Chydorus sphaericus group in the central Palearctic and assess its diversity and biogeographic boundaries. We sequenced nuclear (ITS-2) and mitochondrial (COI) gene regions of Chydorus specimens across the Palearctic and compared them with already available Holarctic sequences. We detected six main clades in the C. sphaericus group in the Palearctic, of which two of the groups are novel. Three of the more divergent clades are geographically widespread. The central portion of Eurasia (the Yenisey River basin) appears to be a narrow zone of secondary contact for phylogroups that expanded from European and Beringian refugia. As such, the previously unsampled central Palearctic represents an important region for understanding the evolutionary consequences of Pleistocene climatic oscillations on the Chydorus sphaericus group.


Introduction
Cladocerans (Branchiopoda) are well-studied crustaceans in aquatic biology [1]. The best studied cladocerans belong to the planktonic genus Daphnia [2][3][4], but a few recent studies have investigated non-planktonic genera [5][6][7][8][9][10]. Planktonic and littoral cladocerans often differ in their resting egg structure and dispersal strategy. Females of the Order Anomopoda Sars place their resting eggs in ephippia-modified exuvia. There is evidence that the ephippia of various pelagic cladocerans are adapted for easy dispersal from one water body to another by waterfowl or by some abiotic vectors [11][12][13]. However, the ephippia of littoral anomopods are heavy and appear adapted to be kept in the water body of the parents [14][15].
The best known genus of the littoral cladocerans is Chydorus Leach, the type genus of the family Chydoridae Dybowski & Grochowski. This family is the most species-rich among the cladocerans. The first taxon of this genus was described by O. F. Müller [16] as Lynceus sphaericus O.F. Müller, 1776. Leach [17] later established the genus Chydorus Leach, 1816 for a single species C. sphaericus. Since the revision of Leach, many authors noted the presence of the genus in different regions of the world, and many new taxa were described [18][19][20][21][22][23][24][25][26][27]. It soon became apparent that taxa of this genus are among most common cladoceran species in many regions of the world. Moreover, Chydorus was detected in a very wide range of water bodies, from minute puddles to large lakes and rivers, and from oligotrophic alpine lakes to hypertrophic lowland reservoirs [28][29]. Chydorus sphaericus was recorded innumerable times in studies of freshwater from the 19th-21st centuries.
Smirnov [30] provided a major revision of the family Chydoridae. He noted the existence of a great number of taxa being synonyms of earlier described species in some groups, including Chydorus. Twenty three species of Chydorus were found to be valid, while some species were regarded as having several subspecies [30]. But the diversity of the genus remained poorly understood at the global level. This made adequate biogeographic conclusions concerning Chydorus difficult.
Another wave of advances occurred in cladoceran systematic biology with the application of genetic methods. In addition to an external method to evaluate the correctness of the separation and grouping of some taxa, we now have a powerful tool to understand their evolutionary history, reconstructing genealogical relationships within each group based on molecular markers. Such works for chydorids are only beginning to appear [6][7].
The most comprehensive study of chydorid taxonomy, biogeography and evolution based on molecular markers was conducted by Belyaeva & Taylor [7] for the Holarctic Chydorus sphaericus-brevilabris complex. The authors found that there are at least seven taxa in the C. sphaericus group in the Holarctic, including C. biovatus. In agreement with the opinion of Frey [31,37], C. brevilabris group was found to be a sister group to the sphaericus-group based on molecular methods. But Belyaeva & Taylor [7] lacked material from a large geographic area, the Eastern Palearctic (only a few populations from the Eastern coast of Asia were represented). Due to this knowledge gap, it is presently difficult to make detailed conclusions about the biogeography, diversity and phylogeography of Palearctic Chydorus.
Other recent investigations on the diversity of the genus Chydorus resulted from the DNA barcoding of the cladocerans in the New World [53][54][55][56][57]. These data were largely composed of samples beyond the Chydorus sphaericus-group.
Our study aims to address the taxonomic and biogeographic knowledge gap for the Chydorus sphaericus species group in the Eastern Palearctic. We combine original DNA sequences of the mitochondrial cytochrome c oxidase (COI) and nuclear rDNA internal transcribed spacer (ITS-2) genes with those deposited in Genbank to assess phylogeography, biogeographic boundaries and taxonomic representation of this group in the Palearctic.

Collection of samples
Field collection in Russia was carried out by our team or by our colleagues under the fulfillment of the governmental program "Ecology and biodiversity of aquatic ecosystems and invasions of alien species" ( 0109-2014-0008), or during engineering ecology surveys with governmental permission to collect such samples from public property. Sampling in the natural reserves of Russia was conducted with special permission of their Directors (A.L. Strelnikov, Komandorsky State Natural Reserve; T. I. Shpilenok, Kronotsky Biosphere Reserve; Yu.P. Suschitsky, State Khankaisky Biosphere Reserve). All collected samples were listed in special reports to the administration of the reserves. Verbal permissions to collect in private farm ponds were obtained from local owners. Field collections in South Korea were carried out as part of the fulfillment of the program from the National Institute of Biological Resources with governmental permission to collect on public property. Samples from Norway were provided by I. Diemantovica, under permission to collect as a member of a Norwegian government institute. The field studies and collections did not involve endangered or protected species.
Specimens were collected by plankton nets with a diameter of 0.20-0.4 m and a mesh size of 30-50 μm, or rectangular dip nets of the same mesh size but a width of 0.2-0.3 m. Collections were preserved in 90-96% alcohol. Before the start of genetic studies, ten specimens of Chydorus from each locality were preliminarily identified by morphological characters [7,50]. Only members of the C. sphaericus group were incorporated into this study. One to five specimens from each of 22 localities through whole Palaearctic were studied here.

DNA sequencing
Genomic DNA was extracted using a Wizard Genomic DNA Purification Kit (Promega Corporation, Madison, WI, USA) according to the manufacturer's instructions. The polymerase chain reaction (PCR) was used to amplify a 461 bp fragment of the 5' region of the mitochondrial cytochrome c oxidase subunit I gene (COI) using the primers Chy-f and Chy-r [7]. A subsample of 14 individuals was sequenced for the nuclear fragment, which included the complete sequence of internal transcribed spacer 2 (ITS-2) as well as small partial sequences of 5.8S and 28S ribosomal RNA genes, using the primers 5.8SF [58] and D2r [59]. In the text below we refer to the nuclear fragment as ITS-2 for convenience. The 25 μl PCR reaction consisted of 2 μl of genomic DNA, 8.5 μl of double-distilled H 2 O, 1 μl of each primer (10 mM) and 5 μl PCR 5x Taq Screen-Mix-HS (Evrogen, Moscow, Russia). The PCR conditions for the COI amplification were 1 cycle of 5 min at 95˚C, 40 cycles of 40 s at 94˚C (denaturation), 30 s at 50˚C (annealing) and 90 s at 72˚C (extension), followed by 1 cycle of 7 min at 72˚C. The PCR conditions for the amplification of ITS-2 were the same, but the annealing temperature was 60˚C [7]. PCR products were purified by precipitation in ethanol and sequenced bi-directionally on an ABI 3730 DNA Analyzer with ABI PRISM BigDye Terminator v. 3.1 sequencing kit (Applied Biosystems, USA). A single consensus sequence was assembled using the forward and reverse sequences using CodonCode Aligner v. 6.0.2 (CodonCode Corp, USA). DNA sequences were submitted to the NCBI Gen-Bank database (accession Nos. KX431583-KX431669 for COI and KX 448799-KX448812 for ITS-2 sequences).

Phylogenetic analyses
The initial authenticity of the sequences was verified by BLAST (blastn) comparisons. The mitochondrial sequences were translated using the invertebrate mitochondrial genetic code to verify the open reading frame. The original sequences from the present study (S1 Table) and known sequences from GenBank (S2 Table) assigned to the C. sphaericus-brevilabris complex were submitted for alignment. The sequences KC616878, KC617532, KC617533, KC617535 by Prosser et al. [58] were excluded from our analysis as they apparently do not belong to the genus Chydorus. Sequences were edited and assembled in the UGENE v. 1.23 package [60]. The DNA sequences were first automatically aligned using the Clustal Omega algorithm [61] using default options and then manually edited. The following population genetics values were estimated: the number of polymorphic sites (S), number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (Pi) and average number of nucleotide differences (k). All calculations were performed using DnaSP v. 5.1 [62] and MEGA v.7 [63].
For tests of neutrality, we calculated Tajima's D [64] and Fu's Fs [65] using Arlequin v. 3.5 [66] with 10000 permutations. Mismatch distribution was constructed for each geographic population to test a model of exponential population growth [67]. A goodness of fit test was performed to test the validity of the sudden expansion model, using a parametric bootstrap approach based on the sum of square deviations (SSD) between the observed and simulated mismatch distributions, with P-values calculated as the proportion of simulations producing an SSD larger than or equal to the observed SSD. The demographic parameter Tau was estimated using a generalized nonlinear least square approach, and the confidence interval of this parameter was computed using a parametric bootstrap with 10000 replicates in Arlequin 3.5.
The best-fitting models of nucleotide substitution were selected in jModelTest 2.1.7 [68] based on the likelihood scores for 88 different models and the Akaike information criterion [69]. Within-and among-clade distances were calculated in MEGA7 using the Tamura 3-parameter model [70] with the shape parameter estimated by jModelTest for COI and the Kimura 2-parameter model [71] for ITS-2. Phylogenetic analyses were performed separately for each gene.
We used sequences from the chydorid genus Pleuroxus for outgroup rooting of the phylogenetic trees. The maximum likelihood (ML) phylogeny reconstruction was performed with MEGA7 and 1000 nonparametric bootstrap pseudoreplicates. Maximum parsimony (MP) analyses were performed in PAUP Ã 4.0a149 [72]. Heuristic MP searches were done using equal weighting, 10 random sequence addition replicates and tree bisection and reconnection (TBR) branch swapping. Non-parametric bootstrapping was performed to assess the nodal support using 1000 pseudoreplicates for MP. Bayesian analyses (BI) were performed in MrBayes v.3.2.6 [73]. Four independent Markov chain Monte Carlo (MCMC) analyses were run simultaneously for 10 million generations and sampled every 1000 generations. The first 25% of the generations were discarded as the burn-in and a 50% majority rule consensus tree was calculated from the remaining trees. The mean genetic sequence divergence between major phylogroups was calculated in MEGA7 using models estimated by jModelTest with pairwise deletion of gaps.
The program Network 5 was employed to construct a median-joining haplotype network that is based on minimum spanning trees in which median vectors representing extinct central or unsampled haplotypes were added following the maximum-parsimony principle [74]. Central haplotypes were identified by their internal⁄central positions in the network, their high frequencies and the number of low-frequency haplotypes derived from them.
Both data sets were tested for recombination events using GARD software [75]. The model selection tool available on the web-server [76] was used to obtain the input nucleotide substitution model. The other settings were: general discrete model for rate variation with four rate classes.

Sequence variation and alignment
The overall sequencing success rate for the mitochondrial gene was about 85%, negative results appeared due to improper sample preservation. Of 144 specimens from 22 localities (Fig 1) sequenced for COI, 64 unique haplotypes were detected. A summary of the diversity of this marker is presented in Tables 1 and 2, but only unique sequences were included in the phylogenetic analyses. The COI alignment for all ingroups and Pleuroxus was 461 bp long, unambiguous and contained no indels or reading frame disruptions. There were 87 variable and 75 parsimony informative sites, with 93 total substitutions. The average base composition was as follows: A-T = 61%. The observed A-T content around 60% is common for the COI region of chydorids [6]. The translated amino-acid alignment had 153 characters, of which 8 were variable and 6 parsimony informative. The overall transition/transversion bias was R = 6.554.
In ITS-2 sequences, there were both very variable regions, and very conserved regions which displayed no variation. The variable regions were difficult to align and we routinely observed all individual sequences in the alignment. We controlled visually and edited manually all dubious substitutions, returning to the original electropherograms in each dubious case. Only 14 ITS-2 sequences were homozygous in ITS-2. These gave good quality sequences with no or few ambiguities. The rest of the sequences were unreadable due to within-individual variation in the sequence and its length, resulting in ca. 60% of double electropherogram peaks. Such ITS-2 heterozygotes were not associated with a particular clade, but rather they

COI phylogenetic tree
The best-fit model selected by jModeltest for the COI data set was Tamura 3-parameter (T92 +G+I) with a gamma distribution and a parameter for the proportion of invariable sites. The T92 distances among COI ingroup taxa varied between 0% and 12%.
Original sequences together with the GenBank sequences could be associated with 11 large clades of Chydorus in ML search (Fig 2). Trees constructed by different methods were congruent, but the clade support differed among analyses. In general, the clades having strong support in ML also had strong support in BI. Deep branches had moderate to strong support in ML and BI, which is not usual for trees based on COI alone. Among the apparent 11 phylogroups, eight belong to the C. sphaericus-brevilabris complex, while three clades (C. pubescens and two unassigned taxa from Genbank) are beyond this species complex. We named the main clades from Eurasia according to Belyaeva & Taylor [7]: A1-A3. The North American clades are also occurring predominantly in Eurasia (although the A3 clade was also found in northwestern North America, east to Manitoba). Two latter clades form C. sphaericus group, and our work is concentrated on the third group (A1 + A2 + A3) only. The clade A1 (Chydorus sphaericus s.str.) is distributed in Europe and Western Siberia, with two sub-clades: widely distributed A1_1 and locally distributed A1_2 which has been recorded to date in Iceland and Greenland only (although it could be more widespread than is presently accepted)-A1_2 was introduced recently to Australia [77]. The clade A2 is present in Northern Europe and Western Siberia, and widely distributed in Eastern Siberia and Far East, including Kamchatka Peninsula and Bering Island. There are three sub-clades (Fig 3): widely distributed A2_1, A2_2 from Kamchatka and Bering Island only, and A2_3 represented to date by only a single population in Yakutia (Eastern Siberia).
The clade A3 is distributed in Far East of Russia, Eastern Siberia, Korea and Japan-genetic distances between specimens from very distant territories within this clade are very small.

ITS-2 gene phylogenetic tree
The best-fit model selected by jModeltest for the ITS-2 fragment data set was Kimura 2-parameter (K80+G) with a gamma distribution with a lightest AIC weight. GARD detected no significant recombination breakpoints by either criterion used. The ITS-2 sequence divergences (K80) varied between 0% and 23% among ingroup taxa.
Only the group A1+A2+A3 is discussed here, because it contains our original sequences, see the discussion of other clades in Belyaeva & Taylor [7]. There are two main clades in this group: A3, and A1+A2 grouped in this tree with American clade A4 (Fig 4). But the support of such grouping is low, the nuclear gene does not support differentiation of the latter two mitochondrial clades. In this large group, the mitochondrial clade A2_3 is not represented (it was among bad sequences, see above), while the clade A2_2 forms a unique branch. Overall, the nuclear phylogeny provided even less resolution than the mitochondrial phylogeny.

Demographic history: mitochondrial median haplotype network
The median-joining network for 103 COI haplotypes (Fig 5) corroborated the ML and BI phylogenetic trees, but the A1 cluster is not monophyletic: clades A1_1 and A1_2 form unique groups. Although all A2 sub-clades form a monophyletic group, this clade is also clearly subdivided into three separate sub-clades, A2_1, A2_2 and A2_3, what also is an evidence of a strong isolation of the sub-clades of A2 phylogroup. The networks for the A1_1, A2_1, A2_2 and A3 clades showed star-like topology, suggesting recent population expansion events.
The center of the haplotype network for A1_1 clade is occupied by the haplotype H1, present in Germany, several regions of southern half to center of European Russia (Moscow Area,  Vladimir Area, Saratov Area, Krasnodar Territory, Kalmykia Autonomous Republic) and Khakassia Autonomous Republic in SE portion of Western Siberia. A Haplotype H17 (from Germany and Norway) is in the center of a smaller star near H1.
The center of the haplotype network for A2_1 clade is occupied by the haplotype H76, present in the subarctic portion of European Russia (Nenets Autonomous Area and the Arkhangelsk Area). The center of the haplotype network for A2_2 clade is occupied by the haplotype H42, present in Kamchatka Peninsula only, but absent in Bering Island.
The center of the haplotype network for A3 clade is occupied by the haplotype H99, present in the Asian Beringia (Chukotka Peninsula, Kamchatka Peninsula), South of Russian Far East (Primorie Territory), Eastern Siberia (Krasnoyarsk Territory, Yakutia Autonomous Republic), American Beringia (Alaska) and Arctic Canada (Manitoba).   Tajima's D and mismatch distribution analysis supported the population expansion hypothesis for all these clades.

Diversity and phylogeny of the Chydorus sphaericus group in the Northern Palearctic
Our study markedly expands and clarifies the biogeographic boundaries of the three large phylogroups [7] of the Chydorus sphaericus complex in Northern Eurasia. Continents are imperfect predictors of boundaries. Instead, several geographic boundaries are apparent within continents. Clades A1 and A3 (which also occurs in the western Nearctic), for example are broadly distributed in the Palearctic with a narrow geographic overlap at the center of the continent. While no new large clades were detected, regional subclades were found within the A1 and A2 clades. It is clear now that even a more narrow regionalism has persisted in this cladoceran species complex subsequent to the last glaciation as it was suggested before [7].
The existence of unique clades in the Eastern Palearctic has been detected in several other cladoceran groups [9,[78][79][80]. So the apparent increased diversity of this region seems relevant in both littoral and planktonic cladocerans. Interestingly, the transitional zone between the western A1 and eastern A3 clades occurs at the Yenisey River basin (grey shading in Fig 1). A similar genetic transition zone has been detected for the cladoceran genus Moina [10]. This region of the Palearctic is a candidate for a suture zone [81][82] for the cladoceran taxa, and even for all freshwater fauna. Surprisingly, specimens from three species (A1, A2_1 and A3) have been detected here in a single oxbow pond of the River Abakan, part of the Yenisey River system. It appears that the Yenisey River, separating Western and Eastern Siberia, is an approximate border line between eastern and western cladoceran faunistic zones-morphological evidence had previously supported this boundary [83]; see also the "line Baikal-Taimyr" [84].
It remains to be seen how taxonomically widespread the Palearctic boundary found here is in other freshwater zooplankton. Earlier it was concluded that the cladoceran fauna in Western Siberia is very similar to that of Europe [83][84], and some genetic data confirmed this opinion [85]. At the same time, some specific taxa were detected by morphological methods in the Eastern Siberia and preliminarily identified as "American" ones [50,86]. The Eastern Siberian-Nearctic connection has been detected in previous genetic studies in other cladocerans [7,[78][79][80], but Eastern Siberia is still poorly represented. Thus far, only Daphnia of this region is just beginning to be studied with the use of molecular markers [87][88].
Our current findings mean that more work must be carried out to reconcile morphological and genetic information for the Chydorus sphaericus group. For example, Klimovsky & Kotov [50] conducted a morphological analysis of Central Yakutian populations and failed to detect C. cf. sphaericus while detecting two separate taxa: C. belyaevae  and the C. cf. biovatus, previously found only in the Nearctic zone. Those findings need to be reconciled with the three clades that we detected here by a joint morphological-genetic study. Also we know little of the reproductive boundaries of the subclades that we discovered. These subclades can co-occur (e.g., A2_1 and A2_3 both occur in Yakutia) and so are not geographic subspecies. More nuclear markers and morphological study are required to test for reproductive boundaries in the overlap zone that we have identified. We have but a few sequences of the nuclear ITS gene, and we cannot exclude the possibility of hybridization between different phylogroups, as revealed in other groups of the Cladocera [89][90][91][92][93]. More efforts are necessary to reveal hybrids between different phylogroups of the C. sphaericus group, or to confirm their absence. No conflict of the mitochondrial and nuclear gene was detected by us, but (1) we have a limited number of sequences of the nuclear gene; (2) and in general ITS-2 lacks differentiation of the A2 and A1 clades.

Phylogeography of the Chydorus sphaericus group in Northern Palearctic
Understanding the timescale of biogeographic events in the Cladocera has been severely limited by a dearth of fossils for calibrating molecular clocks (but see Kotov [15,94]). Because the characters that separate the clades (male morphology and number of resting eggs) are unlikely to be detected in the fossil or in the Holocene subfossil record [95][96], calibration options for the Chydorus sphaericus-brevilabris complex remain poor. Calibrations from mutation accumulation experiments are rare and those for the Cladocera appear to result in severe underestimates for deeper divergences [97]. Belyaeva & Taylor [7] calibrated a molecular clock for the clades of Chydorus studied here based on a general arthropod calibration for 16S rRNA molecular clocks. They estimated the common ancestor of the clade A1+A2 as 2-10 MYA, from the late Pliocene-early Pleistocene, and the common ancestor of clades A1+A2 and A3 to be 5-15 MYA, from the late Miocene. Dates estimated with calibrations from other invertebrate lineages can also have pronounced error, i.e. a ten-fold difference between sister clades can occur [98]. Also, such rates decrease as a function of time [99][100][101].
Even with the error of the existing molecular clock estimations we can be reasonably certain that clades A2_1, A2_2 and A2_3, are younger than the geological events that disrupted the supercontinents [11,102] and caused the mass extinction of the cladoceran fauna in the middle Caenozoic [103]. But, it is clear from our study that the groups were old enough to have been affected by Pleistocene climatic oscillations.
The star-like patterns found in the mitochondrial haplotype network may be due to recent postglacial recolonization events [92,98,104]. Such star-like patterns in C. sphaericus group are obviously associated with very distant regions of Eurasia. The network suggests that the complex survived recent glaciation in separate refugia. H1, the central haplotype of the A1_1, clade is distributed predominantly in the territories of Europe, affected by glaciation during the Quaternary maximum, but not glaciated during the Last Glacial Maximum [105]. Presumably, the A1 phylogroup survived in a refugium in the southern portion of European Russia (and, probably, in the southern portion of "non-Russian" Europe, which is poorly sampled).
Interpretation of the fate of A2 phylogroup is more complicated. The central haplotype of A2_1 is present now in the territory (Nenets Autonomous Area and Arkhangelsk Area) with a complicated geological history. This region was covered by an ice shield in the Middle Pleistocene [105][106]. According to the most modern reconstructions [106], the region was covered by an ice sheet in the Early Weichselian, but free of ice in the Late Weichselian. Earlier authors proposed contrasting reconstructions [105][106][107][108]. But, in any case some "cryptic northern refugia" [109][110] may have existed for A2_1 in the late Pleistocene. Remarkably, distinct mitochondrial lineages of the Daphnia longispina complex were found in the same area, which also suggests the existence of a local/regional refugium in this area [111].
The Weichselian ice sheet was associated with huge proglacial lakes (ice-dammed lakes) whose drainage changed from south-north to a north-south direction [112][113]. Such water movement may have enabled dispersal of cladocerans in a southern direction. Our data demonstrate that this region of proglacial lakes may have acted a source for recolonization. A similar scenario has been proposed for the cladoceran genus Polyphemus by Xu et al. [79].
The interpretation of biogeography of clades A1_2 (Greenland and Iceland) and A2_3 (Eastern Siberia), as well as some other small clades (i.e. East Siberian-with upper position in our network-branch of the subclade of A2_2), is not so obvious. They could be relict clades, largely erased by most recent Pleistocene glacial cycle(s), but testing of such hypothesis requires further studies, i.e. intensive sampling in Greenland, Iceland and Eastern Siberia.
The A2_2 and A3 phylogroups appear to have survived in the Beringian refugium. Therefore, a northern refugium was very important for survival of diversity in the complex. The Beringian refugium has been proposed for many other species of cladocerans [79,114] and other freshwater and terrestrial animals [82,[115][116].
It is remarkable that the inferred central haplotype of the A2_2 clade remains undetected on Bering Island, which has a very scarce cladoceran fauna [117]. According to our network, this island was colonized independently by several (at least three) lineages. Although Bering Island is located relatively close to Kamchatka Peninsula, it is separated by a very deep strait. Bering Island has existed at least since the Eocene [118], and was never incorporated into the Beringian land bridge. Therefore, the dispersal hypothesis appears to be supported by the geological history of this region. Earlier it was shown that this island was also colonized by the cladocerans of the genus Eurycercus from the Nearctic, possibly by migratory waterfowl [119]. Chydorus cf. sphaericus of the clade A3 may have colonized Bering Island by the same mode.

Conclusions
We found six different phylogroups (A1_1, A1_2, A2_1, A2_2, A2_3 and A3) of the C. sphaericus group with different phylogeographic histories in northern Eurasia. Among them four (A1_1, A2_1, A2_2 and A3) passed through a relatively recent population expansion, probably associated with the recolonization of territories previously affected by glaciation. Europe was an apparent source for recolonization of two clades (A1 and A2_1), and Beringia for two other clades (A2_2 and A3). The central portion of Eurasia (the Yenisey River basin) appears to be a zone of secondary contact for European and Beringian phylogroups, where we detected their co-occurrence in the same water body. This zone may prove useful in testing reproductive boundaries between clades that had been previously isolated in different refugia. Our results reveal that the northern Palearctic is an important region for understanding the evolutionary outcomes and recolonization routes following Pleistocene glaciation for Holarctic freshwater invertebrates.
Supporting Information S1 Table. Complete list of sequences obtained in this study with information on locality and the GenBank accession for COI and ITS-2 sequences number provided for each specimen. Clade designations correspond to those in other tables. AR in list of states = Autonomous Republic. (DOC) S2 Table. List of sequences from the GenBank that were used in our study. (DOC) Biosphere Reserve; Yu.P. Suschitsky, State Khankaisky Biosphere Reserve), we also thank the Administration of these Natural Reservoirs for the assistance during the sampling; field collections in South Korea were carried out in frames of fulfillment of the program from the National Institute of Biological Resources of Korea.