The authenticity of the sequences was verified by BLAST comparisons with published D. magna sequences [49]. COI Sequence JF821194 from Turkey was excluded from the analysis because a preliminary tree suggested that this sequence belongs to a separate clade, which, however, could be an artifact and needs further confirmation. Sequences were edited and assembled in uGene v.1.26 [50]. The original sequences of the present study (S1 Table) and sequences from previous publications [51,52,53,54,55,56,57,58,59,60] deposited to the GenBank (S2 Table) were used for multi-sequence alignments. The sequences were first automatically aligned using the T-Coffee algorithm [61] with default options of the uGene package. Due to the strong variability of the COI flanking regions in Crustacea, no universal primers for barcoding exist, in contrast, for example, to fishes [62]. As a consequence, sequences from GenBank had been amplified with varying primers, explaining the variation in fragment lengths. Hence, all COI sequences were cropped to the minimal overlapping length of 563 bp.

To analyze the genetic variation among samples, we estimated the following parameters for each gene fragment and for each species separately: number of haplotypes (Nh), number of variable (polymorphic) sites (Nv), number of parsimony informative sites (Np), haplotype diversity (Hd), nucleotide diversity (Pi), and average number of nucleotide differences (k) [63,64]. We also carried out neutrality tests and assessed mismatch distributions. All analyses were performed in DnaSP v.5.1 [65] and MEGA v.7 [66].

The best-fitting model of nucleotide substitution was selected using the ModelFinder web application [67], based on likelihood scores for 154 different models and the Bayesian information criterion [68]. Within- and among-clade distances were calculated in MEGA using p-distance [69]. Phylogenetic trees were constructed using D. similis and D. sinensis as outgroup species. A maximum likelihood (ML) phylogenetic reconstruction was performed using the IQ-TREE web server [70] and ultrafast bootstrap [71] resampled 10000 times. Maximum parsimony (MP) analyses were performed in PAUP*4.0a152 [72]. A heuristic MP searches was done using equal weighting, 10 random sequence addition replicates and TBR branch swapping. Non-parametric bootstrapping was performed to assess the nodal support, using 1000 pseudoreplicates for MP. Bayesian analyses (BI) were performed in BEAST v2.4.6 [73]. Six independent Markov chain Monte Carlo (MCMC) analyses were run simultaneously for 50 million generations and sampled every 1000 generations. The first 50% of the generations were discarded as burn-in. A 50% majority rule consensus tree was generated from the remaining trees, and the posterior probability of each node was estimated as the percentage of trees recovering any particular node. MP and BI analyses were performed on the computer cluster CIPRES Science Gateway v.3.3 [74].

To establish whether there was any evidence for cryptic species, we performed a Bayesian generalized mixed Yule coalescent (bGMYC) analysis [75] with a threshold of 0.5, using the last 100 trees of the BEAST MCMC file and the GMYC-web server [76]. Furthermore, we used our haplotypes to construct a network in the program PopART [77] to show relationships among the individuals sampled from different locations. The TCS algorithm was selected for this network, based on the implemented statistical parsimony [78]. Second, the Network package (version 5, http://www.fluxus-technology.com) was used to construct a median-joining haplotype network with Steiner maximum parsimony post-processing [79], putting equal weight on each variable nucleotide site. Before constructing this network, data were processed using the star contraction algorithm [80] to reduce of number of similar haplotypes and hence to reduce the complexity of the network. As an additional test for presence of distinct clades, we used the program ABGD [81] to assess p-distances [69] with 100 replicate steps in this reduced-complexity network.

The reduced-complexity network was also used for a nested clade analysis, which identifies significant phylogroups, and to test for historical demographic processes within the major phylogroups. The nested clade phylogeographic analysis (NCPA) was performed in the program ANeCA v.1.2 [82] with integrated module GeoDis v.2.6 for calculation based on a new approach to the interpretation of biological processes according to [83]. Subsequently, to test for neutrality within individual clades, Tajima’s D tests [84] and Fu’s FS tests [85] were performed in Arlequin v.3.5 [86] with 1000 permutations. We also performed R2-tests [87] in DnaSP, again with 1000 permutations. Mismatch distributions were investigated for some groups to evaluate a model of exponential population growth [88]. For the latter, a goodness of fit test was performed, using a parametric bootstrap approach based on the sum of squared deviations (SSD) between the observed and simulated mismatch distributions [89]. The demographic parameter Tau was estimated using a generalized nonlinear least square approach, and the confidence interval of this parameter was computed using parametric bootstrap with 1000 replicates in Arlequin.

Divergence times were estimated using a relaxed molecular clock approach with uncorrelated lognormal distributions of branch rates in the program BEST v.1.6.1, following methods and calibrations of [90]. A second variant of this analysis was done using the additional calibration point of [91] for the Dapnia/Ctenodaphnia split (145 MYA), which suggests that the speed of nucleotide substitutions is significantly lower than suggested by [90].