Patterns of Cryptic Diversity and Phylogeography in Four Freshwater Copepod Crustaceans in European Lakes

: Comparative phylogeography has become a powerful approach in exploring hidden or cryptic diversity within widespread species and understanding how historical and biogeographical factors shape the modern patterns of their distribution. Most comparative phylogeographic studies so far focus on terrestrial and vertebrate taxa, while aquatic invertebrates (and especially freshwater invertebrates) remain unstudied. In this article, we explore and compare the patterns of molecular diversity and phylogeographic structure of four widespread freshwater copepod crustaceans in European water bodies: the harpacticoids Attheyella crassa, Canthocamptus staphylinus and Nitokra hibernica , and the cyclopoid Eucyclops serrulatus, using sequence data from mtDNA COI and nuclear ITS/18S rRNA genes. The three taxa A. crassa , C. staphylinus and E. serrulatus each consist of deeply diverged clusters and are deemed to represent complexes of species with largely (but not completely) non-overlapping distributions, while in N. hibernica only little differentiation was found, which may however reﬂect the geographically more restricted sampling. However, the geographical patterns of subdivision differ. The divisions in A. crassa and E. serrulatus follow an east–west pattern in Northern Europe whereas that in C. staphylinus has more of a north–south pattern, with a distinct Fennoscandian clade. The deep mitochondrial splits among populations of A. crassa , C. staphylinus and E. serrulatus (model-corrected distances 26–36%) suggest that divergence of the lineages predate the Pleistocene glaciations. This study provides an insight into cryptic diversity and biogeographic distribution of freshwater copepods. visualization, E.K.; A.H.;


Introduction
Freshwater animals generally comprise well-defined populations separated by unsuitable habitats and geographical barriers, and often manifest remarkable intraspecific diversity, facilitated by limited gene flow between populations [1]. Studies on the genetic diversity of widespread freshwater taxa, largely applying DNA barcoding methods, are increasingly revealing complex patterns of cryptic variation. Apart from investigations of freshwater fishes [2] and other macroscopic taxa [3], studies on the phylogeny of smaller freshwater crustaceans, such as cladocerans [4][5][6], isopods [7], amphipods [8][9][10], and copepods [11] have demonstrated the presence of hidden genetic diversity and occurrence of cryptic species across presumed conspecific populations.
The concept of cryptic or hidden species has various definitions in the literature, the most commonly holds that cryptic species are two or more (genetically) distinct species that were classified as a single species due to their morphological similarity [12]. The frequency Here we present a comparative survey of the broad phylogeographic patterns and cryptic diversity in four one of the most widespread freshwater copepods in Northern Europe, Attheyella crassa Sars [53], Canthocamptus staphylinus Jurine [54], Eucyclops serrulatus Fischer [55], and Nitokra hibernica Brady [56]. We use our previous mtDNA sequence data for A. crassa [51], C. staphylinus [52], and E. serrulatus [48] and supplement them with new mtDNA data from several additional populations and with nuclear DNA sequences, and newly analyse mtDNA data for N. hibernica. We analyze the intra-species phylogenies and genetic structuring and compare the patterns of cryptic diversity and distribution to understand how the potential cryptic subdivision is formed, what the contact zones and divergence levels between cryptic species are, and what are the implications on the views of the age of diversity.

Overview of the Studied Taxa
The harpacticoid copepod Attheyella crassa is widespread in Europe in different types of freshwater water bodies, and the overall distribution is Palearctic [43,57,58]. The short life cycle (4-6 weeks) has facilitated its use as an experimental organism in ecotoxicological studies, which have shown that exposure to toxic substances can significantly affect the body length, individual fertility, and genetic diversity of populations [59,60]. Preliminary studies of structuring among several European populations revealed significant differentiation in morphometric characters and marked genetic divergence (c 20% between two main clades [51]. Canthocamptus staphylinus is one of the most widely distributed freshwater harpacticoid copepod taxa in the Palearctic [43,61]. It is considered a stenothermal and psychrophilic species inhabiting water bodies with temperatures from 0 • C to 19 • C [62,63]. As part of their life cycle, adult C. staphylinus rest encysted in the bottom mud. They occur in a wide range of habitats, from small spring pools to large lakes and rivers. Populations are both geographically and ecologically variable in several morphological characteristics, including the structure of the fifth pair of thoracal legs, numbers of spinules on the anal operculum, form of the spermatophore, and the development of the aesthetask on the fourth segment of the antennule [43,52,61]. Two distinct mitochondrial clades (22-23% divergence in COI) were recently reported [52], which raises the question whether the different forms described as C. staphylinus in the literature actually represent a group of closely related species.
Eucyclops serrulatus is a widespread cyclopoid copepod, which was considered a cosmopolitan taxon until recent splits and a comprehensive revision of the E. serrulatus group [64]. In the Palearctic, the range of this broader E. serrulatus group extends from North Africa to continental Europe and central Asia to Eastern Siberia. All records from Australia, North America, and other zoogeographical regions are likely the result of recent invasions [64]. Recent studies from East European mountains and European water bodies showed significant differences in morphology of setae on caudal rami and swimming legs, as well as in mtDNA. Eight clades with high mitochondrial divergences (around 30%) have been recognized [25,48].
Nitokra hibernica is a benthic harpacticoid with broad salinity tolerance, occurring in fresh, brackish, coastal, and estuarine waters. It is usually associated with macrophyte beds occurring in nearshore zones of large rivers and lakes. It has a broad distribution in Eurasia, e.g., in the Black and Caspian Seas, the coasts of the Atlantic and Baltic Seas, inland waters in the South, East, and Central Europe as well as in Central Asia and the Caucasus [65,66]. N. hibernica was probably introduced in ballast water to the Great Lakes of North America, where it has become one of the most common nearshore harpacticoids [67,68]. No previous information on mitochondrial diversity is available.

Sampling
Our sampling comprised a total of 27 European freshwater bodies with 11 populations of A. crassa, 13 of C. staphylinus, 11 of E. serrulatus, and 5 of N. hibernica (Table S1). The locations include individual lakes, rivers or reservoirs in Norway, Sweden, Finland, Estonia, Switzerland, Germany, Austria, Belarus, France, Italy, and various regions of Russia (see later in Results and Discussion). The material includes samples of three species from our previous studies [48,51,52] augmented by new samples from Germany (lakes Plau am See and Stechlin), Austria (lakes Fernsteinsee, Zell am See and Hallstatt), Switzerland (Lake Zurich) and Russia (the Severnaya Dvina and Lena rivers), as well as new data on 5 populations of the previously unassessed harpacticoid N. hibernica. Samples were collected by horizontal sweeps at 1 m depth in near-shore areas using a 100 µm mesh size hand net. Living samples were sorted under a stereomicroscope; copepods were then individually preserved in 96% ethanol. mtDNA data on five additional specimens from four populations were obtained from GenBank (see Table S1).

Molecular Analysis
Genomic DNA was extracted using different protocols for the different taxa. For E. serrulatus these included the salt extraction method of Aljanabi & Martinez, 1997 [69], and use of the Diatom™ DNA Prep 200 Kit (Galart, Russia); for C. staphylinus and N. hibernica the Chelex protocol outlined in Walsh et al., 1991 [70] and described in Kochanova et al., 2018 [52]; and for A. crassa the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany). The primers used in this study are listed in Table S2 in the supplementary material. The PCR mixes and amplification protocols for the mitochondrial COI gene for A. crassa and N. hibernica are described in Kochanova & Gaviria, 2018 [51], for C. staphylinus in Kochanova et al., 2018 [52] and for E. serrulatus in Sukhikh & Alekseev, 2015 [48]. In addition, sequences of nuclear markers, either ITS or 18S rRNA, were obtained for the three species C. staphylinus, A. crassa, and E. serrulatus. The amplification of ITS1-5.8S-ITS2 region of the rRNA operon (for short, ITS) involved initial denaturation at 95 • C (30 s), followed by 38 cycles of denaturing at 95 • C (30 s), annealing at 50 • C (30 s), and extension at 72 • C (70 s), and then a final extension at 72 • C for 7 min. The PCR program for the 18S gene was: 15 min HotStarTaq activation step at 95 • C; 35 cycles of 1 min at 94 • C, 1 min at 55 • C (fragment 1), 59 • C (fragment 2) or 57 • C (fragment 3), 2 min at 72 • C; and 10 min extension hold at 72 • C. PCR products were visualized by electrophoresis in 2% agarose gel and purified with the ExoSap-IT PCR Product Clean-Up kit (Applied Biosystems, Waltham, MA, USA). Sequencing was carried out in both directions, using the BigDye Terminator v3.1 (Life Technology) reagent kit on the ABI PRISM 310 Genetic Analyzer (Applied Biosystems, Waltham, MA, USA).

Analysis of Sequences and Reconstruction of Phylogeny, Haplowebs and Population Structure
The obtained nucleotide sequences were aligned using the ClustalW algorithm [71]. In the nuclear dataset, alleles were not distinguished and double peaks in chromatograms were coded using ambiguity codes. The nucleotide and haplotype diversities, number of polymorphic sites, and mutations were assessed in DnaSP 5.1 [72]. The best-fitting models of nucleotide substitution for both mitochondrial and nuclear datasets were selected in jModelTest v. 2.1.7 based on the likelihood scores for 88 different models and under the Akaike Information Criterion (AICc) [73]. The model selected for COI was GTR+Г (general time reversible model with inter-site rate variation modeled with a gamma distribution) in all species but N. hibernica, for which HKY+Г was chosen (see Results). The model for all nuclear datasets was HKY.
To display the relationships of individual haplotypes, separate median-joining (MJ) haplotype networks for the COI marker and ITS/18S markers were constructed using POPART 1.7 [78]. The homoplasy level parameter (ε) was set at the default value (ε = 0).
An analysis of molecular variance (AMOVA) among and within populations of A. crassa, C. staphylinus, E. serrulatus and N. hibernica was conducted using Arlequin v.3.5.2.2 [79] and presented in Supplementary Material (Table S3). Tajima's D [80] and Fu's F S [81] were calculated for populations with more than four sequences and tested for significant departure from demographic stability with 10,000 coalescent simulations in DnaSP 5.10 [72].

Species Delimitation
For a mechanistic "OTU delimitation" from sequence data of individual loci, we applied three commonly used approaches (the ABGD distance method, coalescence models in the 'splits' and 'bGMYC' packages, and through the Poisson analysis of mPTP) and one (STACEY) for the combined mtDNA and rRNA gene dataset.
ABGD, or Automatic Barcode Gap Discovery (ABGD) [19], was calculated on the webserver Atelier de BioInformatique, France (https://bioinfo.mnhn.fr/abi/public/abgd/ abgdweb.html access on 13 September 2021). Single values for both mitochondrial loci were selected by us: Pmin = 0.001, Pmax = 0.1, Steps = 100, X =10, Nb = 25. GMYC instead applies the coalescence approach based on the general mixed Yule coalescent model [82] and was performed with the 'splits' package [17]. As the input tree, we used an ultrametric tree obtained for each locus in BEAST 2.5.2 programme [83]. The PTP method of species delimitation in turn is based on the Poisson tree processes, for which we used multi-rate Poisson Tree Processes, mPTP [18] on the web server of Heidelberg Institute for Theoretical Studies, Germany (http://mptp.h-its.org/ access on 13 September 2021). As the input tree, we used the phylogenetic ML-tree obtained for each locus.
STACEY v.1.2.4 [84] was used as a method for multilocus phylogenetic taxonomy. This method is a version of the multi-species coalescence model used in the BEAST 2.5.2. program [83]. Bayesian species (partitioned) trees were first constructed in BEAST 2.5.2 using the STACEY template. Final phylogenetic relationships were estimated in four independent runs for the whole data set. Each run consisted of 50 million iterations sampling every 10,000th and discarding the first 50.00 iterations (burn-in). STACEY log files were examined in Tracer v.1.7.1 [85] to assess whether the runs have reached the stationary phase and converged on model parameters (ESS > 400). Support of topologies was evaluated in STACEY by constructing a tree of maximum reliability in TreeAnnotator after the rejection of half of all estimated trees. Species delineation (based on the trees evaluated in STACEY) was carried out using a Java-application 'speciesDA' written by Graham Jones ( http://www.indriid.com/software.html access on 13 September 2021).

Attheyella crassa
In total, an alignment of 37 sequences (individuals) of mtDNA COI gene (545 bp), and one of 17 sequences of rRNA ITS gene (461 bp) were used in the analysis (Table 1). Two distinct mitochondrial clades termed the E (Eastern) and W (Western) clades, were identified from 13 populations in the ML tree ( Figure 1a). The E clade included the single water bodies studied from Finland, Estonia, Sweden, and Belarus, and from Komi, Nenetsia (Pechora and Vychegda Rivers), the lake Glubokoe in Moscow region, and the Rybinsk reservoir in the Yaroslavl Region of Russia. The W clade was made up of individuals occurring in the Swiss Lake Geneva, two lakes in Norway and the Rybinsk reservoir in Russia. Additional sequences of closely related species A. dentata and A. (N.) nordenskioldii formed the separate clade on the COI gene tree. The same division into E and W clades was observed in the ITS gene tree from 17 sequences, including one from the Norwegian Store Le not analyzed for COI ( Figure 1b). The presence of individuals from both E and W clades in the Rybinsk reservoir was identified congruently with the two mitochondrial and nuclear markers. All three single-locus species delimitation methods (ABGD, GMYC and PTP) separated E and W clades in both markers. Further, GMYC and PTP approaches on ITS indicated the Rybinsk and Vänern populations in E clade as additional species, while they showed negligible divergence. The multi-locus approach (STACEY), which was based on the partitioned Bayesian species tree, only recognized the E and W clades as species.
The same division into E and W groups of haplotypes (with 84 mutations between them in COI gene and 23 mutations in ITS gene) is illustrated in the median-joining networks (Figure 1c,d). The average model-corrected distances between the E and W clades were 26.1% for COI and 5.5% for the nuclear ITS marker. For comparison, the divergence between A. crassa and A. dentata and A. (N.) nordenskioldii were about 35%. The intraclade diversity in E and W clades were 6.5% and 1.3% correspondingly for the COI gene and less than 1% for the ITS gene ( Table 2).
The haplotype and nucleotide diversities were higher in the E clade (with the highest values in Pääjärvi lake and lower in Rybinsk reservoir in both clades), values of Tajima's D and Fu's Fs indexes were calculated only for two populations, they were insignificant but negative (Table 3). As an outgroup, sequences of Elaphoidella bidens (accession numbers MW851202-MW851206) for COI gene tree and Tigriopus californicus (KF851276) for ITS gene tree were used. The bars between the trees indicate the "species" identified by different species delimitation methods for each gene (ABGD, GMYC, PTP) and for both genes together (STACEY). Median joining networks of (c) COI and (d) ITS haplotypes from A. crassa. Black dots indicate median vectors (ancestral, unsampled, or extinct haplotypes). Black hatch marks denote the number of mutational steps between haplotypes. Circles are proportional to the haplotype frequency.

Canthocamptus staphylinus
Overall, 59 sequences of COI (606 bp aligned) and 10 of 18S (1576 bp) were obtained ( Table 1). The COI gene tree of C. staphylinus, representing 13 populations, also showed two major clades, nominally C (Continental) and F (Fennoscandian) (Figure 2a). The C clade comprised of ten populations, from Switzerland, Germany, Austria, Estonia, Belarus, and Russia (St. Petersburg). The F clade comprised of three populations, from Finland, Sweden, and Northern Germany. In the broader comparison, the two clades did not cluster together but the sequences of the Asian species C. kitaurensis and C. (B.) longifurcatus were positioned between them. The European C. microstaphylinus made the basal branch of the tree. The 18S gene, analyzed from five populations (Finland, Sweden, Estonia, Switzerland, and Russia) showed almost no variation (one common haplotype and three rare ones separated from them) and no consistent division between the COI clades ( Figure 2b). As an outgroup, sequences of Attheyella crassa (accession numbers MH477663-MH477668) for COI gene tree and A. crassa (EU380307) for 18S were used. A dashed line indicates that the sequence downloaded from GenBank contains ambiguities which are unlikely to affect its position in the tree. The bars between the trees indicate the "species" identified by different species delimitation methods for each gene (ABGD, GMYC, PTP) and for both genes together (STACEY). Median joining networks of (c) COI and (d) 18S haplotypes from C. staphylinus. Black dots indicate median vectors (ancestral, unsampled, or extinct haplotypes). Black hatch marks denote the number of mutational steps between haplotypes. Circles are proportional to the haplotype frequency.
In the COI dataset, the ABGD and GMYC methods identified the C and F clades as different species (as well as C. kitaurensis, C. (B.) longifurcatus and C. microstaphylinus), while the PTP method would further separate the Finnish and Swedish populations in F clade, and the isolated Russian St. Petersburg from the other populations in C clade. No division was implied from the practically invariable 18S gene data. The STACEY approach, considering the two genes together, recognized the C and F clades.
The division into C and F groups of haplotypes (with 102 mutations between them) was similarly illustrated in the COI gene median-joining network, but again not in 18S (Figure 2c,d).
The average genetic distances between clades were 33.2% and 0.07% for the COI and 18S genes, respectively. The divergence among all Canthocamptus species ranged from 15.3% to 46.4% of the COI gene. The average intraclade distances in the F and C clades were 2.8% and 3.6% (COI) and 0.007% (18S) ( Table 4). Haplotype and nucleotide diversities were higher in the C clade than in the F clade, with the highest values in the Orlov pond and lowest in Pääjärvi lake. Tajima's D and Fu's Fs were calculated for four populations: three of them were significant (p < 0.05) a and fourth had negative value of Tajima's D (Table 3).

Eucyclops serrulatus
Here, 38 sequences of COI (654 bp) and 12 sequences of ITS (628 bp) were analyzed ( Table 1). The COI gene tree of E. serrulatus, from ten populations, again displayed two main clades: W (Western) and E (Eastern) (Figure 3a). The W clade comprised of various ponds in Russia (two in St Petersburg), Ukraine, Norway, France, and Italy. The E clade was mainly Russian and recorded in the ponds of St. Petersburg along with samples from three rivers, i.e., Karkalaika (Udmurtia), Severnaya Dvina (Arkhangelsk) and Lena (Yakutia, Siberia), and from the Dniester Liman in Ukraine. The closely related species E. taiwanensis, E. prionophorus, and E. cuatrocienegas formed another clade on the tree. The ITS gene tree only comprised three populations from Russia (two ponds in St. Petersburg and a river in Udmurtia) but showed the same division into two clades (Figure 3b). Both mtDNA and rRNA markers showed that individuals from the two clades co-existed in both ponds from St. Petersburg.
The ABGD and GMYC approaches identified the E and W clades as different species from the mtDNA data, while the PTP method also separated a group of four populations from Russia. From the ITS dataset, all three single-locus methods indicated the three species from Russia: two of them in ponds of St. Petersburg and the third from the river in Udmurtia. STACEY separated two groups, i.e., the E and W clades.
There were 70 mutations between the E and W groups of haplotypes in the COI median-joining network (Figure 3c). In the ITS gene network the specimen from the Tauride population belonging to E and W clades were separated in different haplotypes with 6 mutations between them (Figure 3d). The average COI genetic distances between the two clades were 36.6%, with 2.4% in the ITS gene. The divergence between other Eucylops species ranged from 38.5% to 48.1%. The E and W intraclade distances were 2.1% and 1.3% in COI gene and less than 0.5% in ITS ( Table 5).
The nucleotide and haplotype diversities were higher in the E clade, with the highest values in the Orlov pond. Values of Tajima's D and Fu's Fs for two populations and were insignificant with one instance of negative Tajima's D in Karkalai creek (Table 3). Outgroup sequences were obtained from GenBank, including Cyclops insignis (accession numbers KC627291-KC627294) and C. abysorum (KC627287-KC627290) for COI gene tree and Mesocyclops ogunnus (GQ848502) and M. pehqeiensis (GQ848500) for ITS gene tree. The bars between the trees indicate the "species" identified by different species delimitation methods for each gene (ABGD, GMYC, PTP) and for both genes together (STACEY). Median joining networks of (c) COI and (d) ITS haplotypes from E. serrulatus. Black dots indicate median vectors (ancestral, unsampled, or extinct haplotypes). Black hatch marks denote the number of mutational steps between haplotypes. Circles are proportional to the haplotype frequency.

Nitokra hibernica
This analysis involved 10 COI gene sequences (554 bp) from five populations ( Table 1). The diversity of N. hibernica was lower than in the other taxa, with no deep clades and COI distances around 3% ( Table 6). The basal split was between an Austrian population and another clade comprising the samples from Belarus, Estonia, Finland, and Germany (Figure 4a). N. lacustris formed a sister clade on the tree and Nitocrella achaigae formed another (Figure 4a). The divergence between Nitokra species was 47.3% and between N. hibernica and N. achaigae was 53.4% (Table 6). In the median-joining network the Austrian population was distinct by 15 substitutions. Other populations were closer but showed clear inter-population differences contrasted by within-population homogeneity (maximum 1 substitution within population).

Cryptic Diversification
Three of the four studied taxa, Attheyella crassa, Canthocamptus staphylinus, and Eucyclops serrulatus, show a division into two European clades with strong mitochondrial divergence (26.1-36.6%). The ranges of these sister clades are largely parapatric, although sympatric occurrences were also documented in putative contact zones, in E. serrulatus in St Petersburg, and in A. crassa in Rybinsk Reservoir, Russia. Four sequence-based species delimitation procedures unanimously suggested these sister clades as distinct species. One of the methods (PTP) would even have delineated some populations inside the main clades as species. The fourth taxon studied, Nitokra hibernica, does not show such a deep split, and the divergence between its main clades (3% in COI) rather corresponds to the intra-clade diversity in A. crassa, C. staphylinus, and E. serrulatus (1.3-5.5%). The geographical sampling of N. hibernica was however not as extensive as of the others and not comprehensive in relation to the overall range, and important subdivisions similar to those of other taxa could therefore have been missed.
The presence of clearly diverged sister clades and cryptic species, and in general high levels of intra-and interspecific divergence, have been a common observation in studies of marine copepods [44][45][46][47][86][87][88][89]. Phylogeographic studies of freshwater copepods have demonstrated similar patterns, and their interpretation in the framework of geographical factors and geological events may have been more easily connected with the role of isolation in the speciation process. For instance, in the calanoid Eudiaptomus hadzici four isolated cryptic species were identified from relatively closely situated lakes or groups of lakes in the Western Balkans [21]. Similarly, four geographically separate phylogroups of Neutrodiaptomus tumidus were revealed in Taiwan [76]. In a study of Cyclops, Holynska and Wyngaard [49] pointed out that the distinct distributions of the morphological and genetic sister species C. lacustris, C. bohater, and C. divergens suggest a speciation history connected with geographical isolation. Further, molecular data indicated the subdivision of the freshwater calanoid Hemidiaptomus ingens in the Mediterranean area into three distinct lineages and divided H. roubaui into two geographically isolated species: H. roubaui and H. maroccanus, highlighting ample unknown diversity in morphologically similar widespread species [26]. Similarly, among North American populations of Skistodiaptomus pallidus, three genetically divergent clades were identified [77].
In terms of quantitative divergence estimates, the ranges of interspecific and interclade (cryptic) COI variation found in this study are similar to those in other freshwater copepods: the model-corrected distances of the above-mentioned sister lineages of E. hadzici, N. tumidus, H. ingens, and S. pallidus ranged from 15.2% in N. tumidus to 56.3% in H. ingens (Table 7). For the comparability of studies these estimates have been recalculated using the GTR+Г substitution model similarly to our own data, and aiming to make them proportional to divergence time [27]. The observed p-distances are often almost half as small (also in Table 7), yet also illustrate the parallels among taxa. It should be noted that copepod populations also overall show extreme mitochondrial DNA diversity [90]. Interestingly, this is not only for geographically isolated freshwater populations, but also for marine species, where populations from relatively close locations can demonstrate high variability [11,20,[44][45][46][47]50,[86][87][88][89][90][91][92]. These "intra-species" levels are extraordinary in comparison with vertebrate and invertebrate taxa, which commonly show <3% intraspecific and >2% interspecific variation [93][94][95]. Rawson and Burton [96] explained that the unexpected genetic diversity in copepod mtDNA can be caused by clonal inheritance which reduces the effective population size and increases the likelihood of deleterious mutation fixation. The influences of short life span, fast generation time and peculiar life history of microcrustacean species have also been implicated. As yet there is no theoretical or experimental evidence of the relevance of these factors as accounting for the observations in freshwater taxa. Table 7. Model-corrected COI genetic distances and observed p-distances between cryptic lineages of different freshwater copepod taxa. The rate variation parameter α was estimated for each dataset separately. The "split age" is a transformation of the model distance under a range of COI divergence rates between 1.4%-4% Ma −1 , applied in some other crustacean studies [97][98][99][100][101]. The lower divergence in the nuclear markers is in line with data from previous copepod studies [24,44] and with the generally comparatively slow rates in ribosomal RNA genes [102,103]. The 18S rRNA in particular is too conservative for phylogeographic investigations [24], and showed no informative differentiation in our study. However, the ITS gene was more variable and distinguished the major groups congruently with mtDNA, although with lower sequence divergence (3-7%). It appears a suitable nuclear marker for intraspecific population structure and species delimitation [11]. Most importantly, the congruent information of mtDNA and nuclear markers of the presence of two lineages sympatrically in one lake (in each of A. crassa and E. serrulatus) supports a taxonomic distinction, although more data is desirable for proving reproductive isolation.
The terming of the major clades of A. crassa, C. staphylinus, and E. serrulatus as cryptic species also accords with the absence of phenotypic differences among them. In our previous studies, morphological variability of certain characters among populations of these species was observed [48,51,52]. In A. crassa there was statistical variation in total body length, relative length and width of caudal rami, and the ratio of length and widths of exopod and baseoendopod of the fifth legs [51]. Among populations of C. staphylinus, variable characters included ratios of spine length of the fifth pair of thoracic legs, numbers of spinules on the anal operculum, and the development of the aesthetask of the antennule [52]. Populations of E. serrulatus showed differences in morphology of setae on caudal rami and swimming legs [48]. However, these morphometric characters did not split the groups congruently with the genetic clades and remained variable across all populations in each of three species. No morphological differences were found in populations of N. hibernica (Kochanova, unpubl.) The slow rate of morphological evolution relative to molecular evolution is emerging as a general pattern in copepods [41,91,[104][105][106], and indicates that modern morphological resemblance does not necessarily mirror cladogenetic events [11,77].
The observations on the phylogenetic structure and deep splits among populations of A. crassa, C. staphylinus, and E. serrulatus, together with previous studies on cryptic speciation among freshwater copepods, illustrates the broad underestimation of the true taxonomic diversity. As a rule, and as a clear antithesis to the previous "cosmopolitanism paradigm", molecular studies on putatively widespread freshwater copepods (and other planktonic and meiobenthic organisms [3][4][5][6][7][8][9][10]) show that they indeed represent a high number of different species or complexes of cryptic species, while often remaining indistinguishable by morphological characters. Rather, "regionalism" and "provincialism" of freshwater taxa are replacing the concept of global homogeneity of widespread species [26].

Comparative Phylogeography
Comparative phylogeography is an important approach to understanding the biotic response to past challenges and can provide a compass for contemporary and future changes [28]. For instance, Soltis et al. [107] examined the phylogeographic structures of 157 species of fish, amphibians, birds, reptiles, mammals, molluscs, crustaceans, insects, algae, fungi, moss, trees, and plants from eastern North America and recognized six main recurrent patterns which may generally be attributable to isolation and differentiation during Pleistocene glaciations, although may in some cases be older (Pliocene). Interestingly, many species exhibited distinct patterns that reflect the unique, rather than the shared, aspects of species' phylogeographical histories. Similarly, the comparative phylogeography of oceanic species of cetaceans, fishes, turtles, and plankton [108] showed that the main common phylogeographic patterns were influenced by glacial cycles and changes in sea levels; however, they also demonstrated alternative pathways, such as ecological speciation and parapatric (semi-isolated) divergences. Bowen et al. [108] showed that the general concordance between biogeography and phylogeography indicates that the population-level genetic divergences observed between provinces are a starting point for macroevolutionary divergences between species. This can be crucial for recognizing the hidden cryptic diversity in conservation and natural resource protection and management [109].
At the same time, comparative phylogeographic approaches among closely related and ecologically similar species should be more sensitive to reflect shared processes ('phylogeographic parallelism') [110] and to point out the factors underlying discordant phylogeographic patterns, since the analysis is less biased by major biological differences of the taxa [111]. We know of no previous studies on comparative phylogeography of copepods. But a study of two closely related species of metallic show skink in Tasmania [111] demonstrated remarkably similar phylogeographic patterns suggesting that these species responded similarly to Plio-Pleistocene climate cycling, and that glacial cold and aridity forced them into similar lowland refugial regions throughout the area. Likewise, comparison of phylogeographic patterns of 12 species from desert-dwelling fauna in the Mojave and Sonoran deserts identified six hotspots of high genetic divergence and diversity and suggested that most of the species diversified into distinct lineages prior to Pleistocene climatic changes [31].
In our study, the main lines of subdivision in three European freshwater copepod taxa do not exactly concur, but may be divided into two patterns: Eastern-Western (E-W) in A. crassa and E. serrulatus and Continental-Fennoscandian (C-F) in C. staphylinus ( Figure 5).
A. crassa and E. serrulatus are from different copepod orders but both are thermophilic and prefer summer seasons for active stages of life cycles [43,64]. Similar contacts of eastern and western lineages are known from several terrestrial and aquatic taxa, and they constitute a general pattern in North European biogeography. Generally, the North European fauna has been regarded as a mixture of Siberian and European elements, which have immigrated from non-glaciated refugial areas in the east and the south [4,112,113]. The division of A. crassa runs between Swedish and Norwegian populations, which relates to the hypothesis of Fennoscandian colonization from both north-eastern and south-western directions [114]. The unexpected occurrence of two haplotypes from distant clades in the Rybinsk reservoir might be linked with migration assisted by water birds [115] as the Black Sea Flyway connects the Basin of Volga River with the area of the Western clade [116]. The contact zone between clades of E. serrulatus passes through waterbodies in St. Petersburg, Russia: in two distant ponds (Orlov pond, the type locality of the species, and the pond in Tauride garden), there are representatives of two separate groups.
The structure in C. staphylinus appears as a north vs. south subdivision. This species, in contrast to others, is psychrophilic and can be found in an active form even under ice and at a temperature of 0 • C [62]. It also differs from others by the capability to form resting cysts and encapsulate during unfavorable environmental conditions, for example when the temperature goes below 4 • C or above 15 • C [62]. The limit between the two clades is essentially between Fennoscandia and the Continental part of Europe, in the Baltic Sea. The nucleotide and haplotype diversities in Fennoscandian populations are smaller than on the continent ( Table 3). The obviously consistent north-south pattern may reflect the existence of two climatically adapted taxa whose distributions have been adjusted in parallel through the recurring climatic cycles. It will be interesting to see whether there are real contact areas where these taxa occur together and how their putative ecological and genetic distinctness would be maintained. In Nitokra hibernica, no significant division or divergence was observed which could reflect the more limited sampling of populations in this species. On the other hand, N. hibernica is the only species in the analysis that can inhabit both freshwater and brackish environments as it has a broad salinity tolerance [43]. Such plasticity helps it to occupy a wide variety of water bodies, similar to what has been shown for the calanoid copepod Eurytemora affinis [117]. However, at this stage, we are not able to conclude the phylogeographic structure for this species and further studies are needed.

Demographic History
Apart from the broad-scale divisions and divergence, the intra-clade structure in haplotypic distribution (e.g., 'demographic tests') are informative of the population history, particularly of a post-glacial time scale. In the haplotype networks of A. crassa (Figure 1c,d) E. serrulatus (Figure 3c,d) and particularly of C. staphylinus (Figure 2c) individual populations often comprised of single haplotypes or of a core haplotype with adjacent singletons one step away, and there was little or no overlap between populations (unique coalescence of lake-wise genealogy). Such homotypic or star-like patterns are indicative of population expansion since the (post-glacial) colonization of the lake (a population bottleneck), and are reflected in the negative Tajima's D and Fu's Fs values [80,118,119]. Overall, this kind of structuring is congruent with other copepod species [120]. Indeed, given the environment and reproductive dynamics, population expansions from an initially small stock are not unexpected among copepods or any small aquatic invertebrates [121]. More intriguing is the question of the nature of the initial bottlenecks (rapid coalescence), whether they were connected only with the colonization processes or also with subsequent or preceding population dynamics (cycles of expansion and bottlenecks), or eventually with selective sweeps in the genome.

Age of Diversity
Under general concepts of mitochondrial substitution rates, the deep divergences among the 'conspecific' copepod clades (26-36%) would in the first place suggest distinctly pre-Pleistocene origins for each lineage. While no properly calibrated rate estimates are available for copepods in particular, for other crustaceans estimates of COI rates in the range 1.4% to 4% divergence Ma -1 have been obtained and used as a basis of discussion [97][98][99][100][101]. With these rates, all the three pairs of cryptic species in our data would have diverged already in the Miocene, from 6 to 26 Ma (Table 7). While this appears unexpected in the context of European biogeography, in the context of copepod systematics these estimates are not exceptional: similar or even deeper cryptic lineages are frequently observed [20,21,26,44,45,[89][90][91]104,105]. Several studies on the phylogeography of freshwater copepods (including European populations) have pointed to the Early Miocene as a possible time of initial divergence [21,26,76]. Similarly, Miocene isolation events were implicated for the initial divergences of the North American copepod Skistodiaptomus [77]. Pre-Pleistocene origin is also shown in the freshwater isopod Asellus aquaticus, whose species diversification was interpreted as dating back to the Middle Miocene and was triggered by the Parathethys regression, continuing with several lineages surviving and diversifying through the Pliocene and Pleistocene glaciations and populating the areas free of ice during the interglacials [7,23,101]. We need to stress caution in making generalizations regarding the potentially exceptional dynamics of copepod mtDNA evolution and lack of proper calibration. Nevertheless, even with considerably higher rates (e.g., 10% Ma −1 ) the divergences would appear pre-Pleistocene (Table 7). An implication is that the origin of the newly revealed cryptic diversity was not connected with Pleistocene events. Regarding putative non-monophyly of (e.g.,) C. staphylinus with respect to Asian taxa (Figure 2), the diversification and invasion events should also be viewed in their broader geographic perspective. But at all events the distributions of the component taxa in Europe must have been (repeatedly) adjusted through the Pleistocene dynamics of extirpation, isolation, expansion, and colonization. The data so far may involve unexpectedly little traces of events, e.g., further radiations, at that time scale.

Conclusions
Our study describes contrasting phylogenetic and phylogeographic structures within the four widespread freshwater copepod species Attheyella crassa, Canthocamptus staphylinus, Eucyclops serrulatus, and Nitokra. hibernica in Europe. The first three taxa each represent a complex of cryptic species with an overall high level of genetic diversity but low phenotypic differentiation. N. hibernica in turn shows only low genetic divergence between populations. The two main patterns of distribution of the recognized sister lineages were Eastern-Western (for A. crassa and E. serrulatus) and Continental-Fennoscandian (for C. staphylinus). The deep mitochondrial splits suggest that the lineages could have diverged already in the Micocene or at least before the Pleistocene, prior to major glacial cycles. Their currently largely allopatric distributions point to differential refugial and colonization histories in later Pleistocene times. The pattern of lake-wise haplotype and haplotype cluster distributions in turn reflect more recent colonization dynamics, probably of post-glacial colonization bottlenecks and population expansions.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/d13090448/s1, Table S1: Sampling localities, number of obtained sequences and GenBank accession numbers, Table S2: Primer sequences for PCR amplification and sequencing, Table S3: Analysis of molecular variance (AMOVA) of the studied species among different populations.  Table 1.