Bosminopsis deitersi (Crustacea: Cladocera) as an ancient species group: a revision

Water fleas (Crustacea: Cladocera) of the Family Bosminidae have been studied since the founding of paleolimnology and freshwater ecology. However, one species, Bosminopsis deitersi, stands out for its exceptional multicontinental range and broad ecological requirements. Here we use an integrated morphological and multilocus genetic approach to address the species problem in B. deitersi. We analyzed 32 populations of B. deitersi s. lat. Two nuclear and two mitochondrial loci were used to carry out the bGMYC, mPTP and STACEY algorithms for species delimitation. Detailed morphological study was also carried out across continents. The evidence indicated a widely distributed cryptic species in the Old World (Bosminopsis zernowi) that is genetically divergent from B. deitersi s.str. We revised the taxonomy and redescribed the species in this complex. Our sampling indicated that B. zernowi had weak genetic differentiation across its range. A molecular clock and biogeographic analysis with fossil calibrations suggested a Mesozoic origin for the Bosminopsis deitersi group. Our evidence rejects the single species hypothesis for B. deitersi and is consistent with an ancient species group (potentially Mesozoic) that shows marked morphological conservation. The family Bosminidae, then, has examples of both rapid morphological evolution (Holocene Bosmina), and morphological stasis (Bosminopsis).


INTRODUCTION
demonstrated morphological stasis for the water fleas (Cladocera) based on paleolimnological records from the Quaternary. He later based the paradigm of "noncosmopolitanism" (Frey, 1982(Frey, , 1987b on this apparent long-term stability in morphology. According to non-cosmopolitanism, geographic differentiation occurred mainly due to vicariant events related to the disruption of Pangaea and the dispersal barriers imparted by subsequent continental drift. The process (often termed "continental endemism") now has strong support among "traditional" taxonomists (Van Damme & Kotov, 2016;

Ethics statement
Field collection in public property in Russia does not require permission. Samples from South Korea were collected in the frame of cooperation between A.A. Kotov and the National Institute of Biological Resources of Korea and does not require special permission. The sample from Arkansas, USA was obtained from collections resulting from a previous NSF grant. The samples from Japan, China, and Thailand were provided by our colleagues having permissions to collect them due to their scientific activity in the governmental institutes of the corresponding countries. Formalin samples from Brazil were kept in the Collection of Zoological Museum of Moscow State University for a long time, they were collected before Brazil introduced very strict regulations for sampling. The species were not assessed as endangered at the time of collection and are currently not subject to specific regulations, however all efforts were taken to ensure that the collection and preservation of animals was performed with due consideration of their welfare. The number of individuals taken did not represent a significant proportion of the population present at each site.

Sample collection and morphological analyses
Samples from different localities were collected via small-sized plankton nets (with mesh size 50 µm) and fixed via 4% formaldehyde or 96% ethanol in the fields, using conventional techniques. All samples were initially examined using a stereoscopic microscope LOMO. Individuals of Bosminopsis were initially identified via available references using morphological features (Kotov, 1997a(Kotov, , 1997bRogers et al., 2019). Existing museum samples were used for morphological analysis (see the list of material in Table S1).
The morphology of populations from the Neotropics and the Palaearctic was examined in detail to assess the presence of taxonomically significant characters. Specimens of Bosminopsis from presorted formalin and alcohol samples were selected under a binocular stereoscopic microscope LOMO. They were then studied in toto under optical microscopes (Olympus BX41 or Olympus CХ 41) in a drop of glycerol formaldehyde or a glycerol-ethanol mixture. Then, at least two parthenogenetic females, two ephippial females, and two males (if available) from each locality were dissected under a stereoscopic microscope for appendages and postabdomens. Drawings were prepared using a camera lucida attached to the optical microscopes.
Some individuals from Neotropical and Palaearctic localities were dehydrated in an ethanol series (30%, 50%, 70%, 95%) and 100% acetone and then dried from hexametyldisilazane. Dried specimens were mounted on aluminum stubs, coated by gold in S150A Sputter Coater (Edwards, West Sussex, United Kingdom), and examined under a scanning electron microscope (Vega 3 Tescan Scanning Electron Microscope; TESCAN, Czech Republic).

DNA extraction, amplification and sequencing
Only alcohol samples were used for the genetic analysis. Each specimen was identified by morphological characters (Table S1). Genomic DNA was extracted from single adult females using the Wizard Genomic DNA Purification Kit (Promega Corp., Madison, WI, USA) and QuickExtract DNA Extraction Solution (Epicentre by part Illumina, Inc., Madison, WI, USA) using manufacturer's protocols. Two mitochondrial and two nuclear markers were investigated here: (1) the 5′-fragment of the first subunit of mitochondrial cytochrome oxidase (COI)-a protein-coding marker widely used in DNA barcoding ; (2) the 5′-fragment of the mitochondrial 16S rRNA gene (16S) with a mosaic of highly conservative duplexes and variable loops (Yang et al., 2014); and (3)(4) 5′fragments of the nuclear ribosomal genes (18S rRNA and 28S rRNA). Each fragment contains both long conservative portions and two variable domains. Although these nuclear markers are predominantly used for divergent taxa (Hovmoller, Pape & Kallersjo, 2002), they are informative at the species level for many microcrustaceans (Karabanov et al., 2018).
Primers used for amplification are listed in Table 4. Polymerase chain reactions (PCR) were carried out in a total volume of 20 ml, consisting of 2 ml of genomic DNA solution, 1 ml of each primer (10 mM), 6 ml of double-distilled H 2 O and 10 ml of ready-to-use PCR Master Mix 2X solution (Promega Corp., Madison, WI, USA). PCR products were visualized in a 1.5% agarose gel stained with ethidium bromide and purified by QIAquick Spin Columns (Qiagen Inc., Valencia, CA, USA). The PCR program included a pre-heating of 3 min at 94 С; 40 cycles (initial denaturation of 30 s at 94 С, annealing of 40 s at a specific temperature, an extension of 80 s at 72 С); and a final extension of 5 min at 72 С (Table 4). Each PCR product was sequenced bi-directionally on an ABI 3730 DNA Analyzer (Applied Biosystems, Foster City, CA, USA) using the ABI PRISM BigDye Terminator v.3.1 kit at the Syntol Co, Moscow. The authenticity of the sequences was verified by BLAST comparisons with published cladoceran sequences in mBLAST (Boratyn et al., 2013).

Population analysis, alignment and phylogenetic analysis
Alignment of sequences from each locus was carried out using the MAFFT v.7 algorithm (Katoh, Rozewicki & Yamada, 2019) available on the server of the Computational Biology Research Center, Japan (http://mafft.cbrc.jp). For the protein-coding gene COI, we used translation alignment with the FFT-NS-i strategy. For alignment of the ribosomal-coding loci, we used the Q-INS-i strategy (secondary structure is considered by this algorithm). Linking sequences and their partitioning for subsequent analyses were made in SequenceMatrix v.1.8 (Vaidya, Lohman & Meier, 2011).
The best-fitting models of the nucleotide substitutions for each locus and for linked data were selected using ModelFinder v. 1.6 (Kalyaanamoorthy et al., 2017) at the Center for Integrative Bioinformatics Vienna web-portal, Austria (http://www.iqtree.org) (Trifinopoulos et al., 2016) based on minimal values of the Bayesian information criteria (BIC) (Schwarz, 1978). The BIC model parameters were almost identical to those obtained using the corrected Akaike's information criterion, AICc.
For the COI locus, the substitution model was partitioned by the nucleotide position of codons (1st, 2nd, 3rd). We used the multi-taxon coalescence model "star" in BEAST2 (Heled & Drummond, 2010) with partitioned models (Chernomor, Haeseler & Minh, 2016). Phylogenetic reconstructions based on the maximum likelihood (ML) and Bayesian (BI) methods were made for each gene separately, the full set of mitochondrial genes, the full set of nuclear genes, and for all "unlinked" genetic data. We included sequences with incomplete or missing data as exclusion can reduce the accuracy of phylogenetic reconstruction (Molloy & Warnow, 2018).
We used the IQ-TREE v.1.6.9 algorithm (Nguyen et al., 2015) via a web-portal CIBIV, Austria for ML tree estimation. Each set of sequences was analysed based on the best model found automatically by the W-IQ-TREE (Trifinopoulos et al., 2016). To estimate the branch support values, we used UFboot2 (Hoang et al., 2018). The Topology of the ML trees was evaluated based on PhyML SH-like tests (Shimodaira, 2002), performed in the block Building Phylogenetic Tree in uGENE v.34 (Okonechnikov, Golosova & Fursov, 2012). BI analysis was performed in BEAST2 v.2.6.2 (Bouckaert et al., 2019), with all of the parameters of the substitution model using the BIC criterion from BEAUti v.2.5.2 (Drummond et al., 2012). In each analysis, we conducted four independent runs of MCMC (40M generations and a sampling interval of 10k generations), with effectiveness control in Tracer v.1.7 . A consensus tree based on the maximum clade credibility (MCC) was obtained in TreeAnnotator v.2.5.2 (Drummond et al., 2012) with a burn-in of at least 20%. Because the main clades for BI and ML were congruent, we presented the BI trees, with ML branch support/ BI posterior probabilities for key nodes.
A haplotype network was constructed for Bosminopsis zernowi (the most sampled taxon in this study) in popART v.1.7 with the Integer Neighbor-Joining Network algorithm (Leigh, Bryant & Nakagawa, 2015) and minimal reticulation tolerance.

Cybertaxonomic species delimitation based on DNA data
Our approach to cybertaxonomic taxon delimitation was described in a previous paper (Kotov et al., 2021). An integrated approach based on genetic species delimitation combined with "traditional" morphological taxonomy was used to estimate the species richness. We used the bGMYC, mPTP and STACEY algorithms for species delimitation (Carstens et al., 2013).
The general mixed Yule-coalescent model (GMYC) was made to assign analyzed individuals to the species according to ultrametric time trees derived from single-locus data (Pons et al., 2006). But the "classical" GMYC has limitations (Lohse, 2009). We used the Bayesian GMYC model in the 'bGMYC' package (Reid & Carstens, 2012) for the statistical language "Microsoft R-Open and MKL" 64-bit v. 3.5.3 (http://mran.microsoft. com). Ultrametric trees for each mitochondrial and nuclear datasets were evaluated in BEAST2. For MCMC, we used 50M generations with a sampling interval of 50k trees. We used Tracer v.1.7 to evaluate the convergence of parameters (based on ESS>200). Sequences of Triops and Bosmina were used as outgroups. Sorting, re-rooting of the trees and outgroup deletion was carried out in "R" according to the script of Sweet et al. (2018). For the bGMYC analysis, we randomly selected 100 ultrametric trees from the 1,000 trees after burn-in from BEAST2. The results were accepted as statistically significant at a modified P > 0.99 level.
Analysis of Multi-rate Poisson Tree Processes (mPTP) (Kapli et al., 2017) was performed on the web-server of Heidelberg Institute for Theoretical Studies (http://mptp. h-its.org/). As the input trees, we used the phylogenetic BI trees from BEAST2 and the ML-tree obtained using W-IQ-TREE for each locus. Delimitation results were congruent for separate loci and were composed of mitochondrial and nuclear datasets.
The combined species tree estimation and species delimitation analysis for STACEY (Species Tree And Classification Estimation, Yarely) (Jones, 2017), was made in BEAST2. Genealogical relationships were estimated by STACEY with four independent generations (50M generations of MCMC, sampling of every 10k generation) after incorporating the suggestions from an initial run. STACEY log files were examined with Tracer v.1.7 to evaluate the convergence of parameters based on ESS > 400. Supports for the tree topologies estimated by STACEY were examined by constructing a maximum clade credibility tree using TreeAnnotator v2.6 (part of the BEAST2) after discarding half of all estimated trees. Species delimitations based on the trees estimated by STACEY were assessed using the Jones' java-application speciesDA: http://www.indriid.com/software. html.

Phylogenetic reconstruction and molecular clock
Two approaches were used for molecular clock estimation. A strict molecular clock (Drummond & Bouckaert, 2015). was based on the assumption of a relatively regular mutation rate in mitochondrial genes. The speed of mutation accumulation differs among organisms. For the crustaceans the rate is ca. 0. 11-2.4% per MYA (Knowlton & Weigt, 1998;Schubart, Diesel & Hedges, 1998;Schwentner et al., 2013;Bekker et al., 2018). Apparently this is a very rough estimation (Schwartz & Maresca, 2006). An alternative approach uses paleontological data to calibrate "molecular clocks".
To estimate the probability of molecular clock-like data, we applied a Maximum Likelihood method implemented in MEGA-X v. 10.1.8 (Kumar et al., 2018). A Maximum Likelihood substitution model was estimated for each locus (separately for each nucleotide position for translated genes, and jointly for non-translated fragments). We used the best model as selected by the lowest BIC (Bayesian Information Criterion) scores and an ML tree with Bosmina as the outgroup.
For determination of the relative rate of substitutions, we used both paleontological information (Kotov & Taylor, 2011) and points based on molecular phylogenetic data (Schwentner et al., 2013). As calibration points (with 15% standard deviations), we used the following estimates: Triops/ all groups 340 MYA, Daphnia/ Simocephalus 145 MYA, Cyclestheria groups 120-70 MYA, and Bosmina / Bosminopsis without an exact date. The age of lineage differentiation according to a strict molecular clock model was estimated in BEAST v. 1.10.4 (Suchard et al., 2018) with a Yule speciation model as the most proper for datasets with several potential species (Gernhard, 2008). Four independent runs of 50M generations were done, with each 100k tree sampled. Subsequent analysis was performed as above for BI following the recommendations (Barido-Sottani et al., 2018).

Phylogeographic reconstructions
To test phylogeographic models, we used the packet BioGeoBEARS (Matzke, 2013) with the integrated statistical package of the "R" language in RASP4 v.4.2 . The data set was composed based on the maximum number of geographic localities and representation of all phylogenetic lineages revealed by cybertaxonomic methods. We estimated a mitochondrial phylogenetic tree based on sequences of two genes (COI and 16S) for Bosminopsis cf. deitersi. Objective software limitations allowed us to analyze only 27 sequences from five phylogenetic lineages and seven main geographic regions.
We tested six biogeographic models in BioGeoBEARS (standard dispersion-vicariant, and those with a correction for speciation events +J), estimated according to the AICc_wt criterion (Matzke, 2014). A phylogeny for RASP4 was reconstructed in BEAST2. In each analysis, we conducted four independent runs of MCMC (40M generations, with a sampling interval of 10k generation). The best model according to the maximum AICc_wt value was DIVALIKE+J, which takes into consideration new lineage origin upon colonization by a founder without the existence of a widely distributed ancestor (Clark et al., 2008). For estimates of the age of historical processes, we used an outgroup "calibration". Palaeoreconstruction was performed in GPlates v.2.2  with PALEOMAP PaleoAtlas v.3 by Christopher R. Scotese (https://www.earthbyte.org/ paleomap-paleoatlas-for-gplates).

ABBREVIATIONS
Abbreviations for collections MGU MLcollection of Zoological Museum of Moscow State University.

Phylogenetics and phylogeography
We analyzed 118 specimens from 32 populations belonging to the B. deitersi group. The specimens originated mainly from Eurasia ( Fig. 1), but a single population from North America and a single population from South America were also analyzed (Table S1).
We obtained 58 original sequences of 16S, 15 sequences of COI, 43 sequences of 18S, and 75 sequences of 28S. Populations had a relatively high genetic polymorphism (Table 1). In contrast, the number of haplotypes was small. Each locus had a differing optimal substitution model (Table S2).
Three major clades of the Bosminopsis deitersi group were revealed from a tree based on the mitochondrial dataset (Fig. S1A). The first clade was B. zernowiwidely distributed in Eurasia and represented by two sub-clades (1 and 2). The second clade was B. deitersi distributed in the Americas, it is represented by two sub-clades: 1 in South America (B. deitersi s.str.) and 2 in North America. Both geographic subclades had modest support. Further study is necessary to examine the independent status of North American populations. A third clade (Bosminopsis sp.) was detected from a single population in Thailand.
The tree based on the nuclear dataset (18S + 28S) had a similar topology to the mtDNA tree, but note that nuclear gene sequences were unavailable for Bosminopsis sp. from Thailand. The large clade of B. zernowi was again subdivided into two subclades (1 and 2) with low support. There were some conflicts between mitochondrial and nuclear sequences. Some specimens from the mitochondrial sub-clade 1 belonged to the nuclear sub-clade 2 (they are marked by asterisks in Figs. S1A-S1B). As support for both mitochondrial and nuclear subclades was low, we do not discuss this below. The 18S locus was almost identical in all populations, suggesting the locus is most informative at the genus level. The 28S locus demonstrated substantial variability in the D1 and D2 variable domains and appeared to contain information for taxa within the genus. Based on the neutrality tests and coalescent simulations in DnaSP v.6, we concluded that the most probable demographic model was a "bottleneck" model reflecting historical processes.
The final tree based on combined mitochondrial and nuclear datasets ( Fig. S1C) was fully congruent with the mitochondrial tree-major clades were well-supported. No conflicts were found for ML and BI (with unlinked data) trees among genes or with the consensus tree.
The results of the phylogenetic reconstructions suggested a deep demographic subdivision of the B. deitersi group. The tests of neutrality were consistent with such a division (Fu's Fs<=0 with Tajima D>>0). The most probable demographic process in this group evolution was an expansion with a strong founder effect resulting in strong differentiation between populations. We further explored the genetic diversity within each group and addressed the taxonomic uncertainty within these lineages.
For cybertaxonomic taxon delimitations (Fig. 2), both bGMYC and mPTP (for both mitochondrial and nuclear genes) suggested a deep divergence within the B. deitersi group.  To estimate divergences among selected OTUs, we calculated "simple" uncorrected p-distances for the best sampled locus, 16S (Table 2). Bosmina was the outgroup. Distances among outgroups are ca. two times greater than the maximum distances within the B. deitersi groups. Groups "Eurasia", "Thailand", and "Americas" are well-differentiated, Figure 2 BI multi-locus tree based on the COI + 16S + 18S + 28S sequences, with a summary of results of the cybertaxonomic species delimitation by different methods. Analyses referring are based on mitochondrial (mit.), nuclear (nuc.) and multi-locus datasets (STACEY). Node supports are: UFboot2 (ML) and posterior probabilities (BI), in percent for mitochondrial genes in the numerator and nuclear genes in the denominator. Dashes indicate branches that were not supported by a method.
Full-size  DOI: 10.7717/peerj.11310/fig-2 while differences between two sub-clades of Eurasia are less than 0.5%. The two subclades may result from moderately separated mitochondrial lineages, which are common in cladocerans (Bekker et al., 2018;Kotov et al., 2021). Again, North and South American populations may be separate species, but more sequences are needed to test this hypothesis. A network of 16S mitochondrial haplotypes (Fig. 3) revealed that all populations from Northern Eurasia belonged to only four haplotypes (Fig. 2): haplotype H1 included 73% of studied specimens and distributed from the Volga basin in European Russia to Pacific coast including Sakhalin Island (but not in Korea and Japan). H2 haplotype was endemic to the Yenisey Basin in Eastern Siberia and seemed to be a derivate of H1. Another well-represented haplotype (H3) differed by two substitutions from H1and was associated with a rare haplotype H4. H3 and H4 were detected only in Japan and Korea. Overall, haplotypic differentiation within groups was weak.
The Maximum Likelihood tests (Table 3) suggested that the hypothesis of molecular clocks is not rejected for each locus tested. In general, the topology of the tree constructed for the molecular clock calculations (Fig. 4A) is congruent with the multilocus tree Note: Nnumber of sequences; ntotal number of sites (excluding sites with gaps or missing data); Snumber of segregating (polymorphic) sites; Hdhaplotype diversity; hnumber of haplotypes; Pinucleotide diversity per site; kaverage number of nucleotide differences; Fs -Fu's neutrality statistic (Fu, 1997) ; D -Tajima's D neutrality test (Tajima, 1989), the star-sign is indicated statistical significance P < 0.05; DemModmost likely demographic model by DnaSP v.6 coalescent simulation (PBpopulation bottleneck, PDpopulation decline, n/dnot defined) based on Theta-W, probability P(Sim<Obs) is indicated as a percentage in brackets. described above. Differences in the divergence pattern of the American (B. deitersi) and Thailand (B. sp.) clades may be explained by heterochrony in the appearance and fixation of the substitutions in mitochondrial and nuclear genomes (Vawter & Brown, 1986).
Minor heterochrony is not an insurmountable obstacle for phylogenetic reconstructions (Allio et al., 2017). "Simple" p-distances between B. zernowi and B. deitersi (based on the COI locus), gave a divergence time estimate of 17-30 MYA (p = 0.241). The estimate used an average divergence time for crustaceans of ca. 0.8-1.4 % per 1 Myr, which is similar to the age of divergence based on the coalescent model (Fig. 4A). Based on the time of divergence of the outgroup (Bosmina), we estimated the divergence of the B. deitersi group at around 200 MYA. Using RASP4, the taxa under consideration had an origin consistent with Laurasia (ancestral distribution range (C(ABCGEF)G), see Figs. 4B, 4C). However, note that Gondwanan populations from Africa, Australia, and India were not studied here. An alternative explanation is that the divergence of Bosminopsis sp. could be explained by its Gondwanan proto-range and subsequent colonization of Eurasia (i.e., due to India's continental drift). Bosminopsis sp. (ancestral range (G)) was separated in the Early Table 3 Estimates of evolutionary divergence over sequence pairs between groups. The number of base differences per site from averaging over all sequence pairs between groups are shown. Standard error estimates are shown above the diagonal and were obtained by a bootstrap procedure (100 replicates). This analysis involved 61 nucleotide 16S sequences. All ambiguous positions were removed for each sequence pair (pairwise deletion option). There was a total of 507 positions in the final dataset. Evolutionary analyses were conducted in MEGA-X (Kumar et al., 2018 (Hasegawa, Kishino & Yano, 1985) ; GTRgeneral time reversible, nst = 6 (Rodriguez et al., 1990). Non-uniformity of evolutionary rates among sites may be modeled by using a discrete Gamma distribution (+G) with 5 rate categories and by assuming that a certain fraction of sites are evolutionarily invariable (+I). Cretaceous from a Euro-Asian-American population of the group (CE(AB)GF)). Subsequent history was probably related to the disruption of Laurasia into Eurasian (CA (BCG)) and American (EF) groups of populations also in Cretaceous. Separation of North (F) and South American (E) populations had no ready explanation-it may be associated with the Gondwana-Laurasia split in the Cretaceous or a significantly more recent dispersal event (Neogene). A strong founder effect could then explain the genetic differences between Neogene populations of the two continents. In any event, the divergence of the entire Bosminopsis group is likely very ancient (at least the Early Cretaceous) and potentially affected by the split of proto-continents.

Morphological analysis
Order Checklist of the formal taxa in the genus Bosminopsis  parts of antennae I, mucro, rostrum, ventral valve edge, base of caudal spine covered with small spinules. Ephippial female with egg chamber sculpture represented by large polygons. A strong medial keel on dorsum, strong paired lateral keels well distinguishable from the dorsal view. Adult male with dorsal contour of head humped, head large, with a smooth rostrum and expressed ocular dome, a short mucro always present posteroventrally. Postabdominal claw bears a basal spine comparable in size with the latter. Antenna I free, remarkably curved distally. A relatively long (somewhat longer than exopod itself), curved at tip additional male seta on endopod apical segment in position of  (Fig. 6I) with a coxal portion bearing a long seta and a short seta on a conical elevation, elongated basal segment and short three-segmented exopod and endopod, antennal formula: setae 0-0-3/1-1-3; spines 0-0-1/0-0-1, but apical spines greatly reduced in size. All apical and lateral (on endopod first and second segment) setae subequal in size, covered by fine setules.
Limb I large, its corm conically narrowing distally. Epipodite (Fig. 7D: epp) with two long finger-like projections. Outer distal lobe (Fig. 7D: odl) with two setae of different size, feathered by sparse, long, robust setules. Inner subdistal lobe (in terms of Kotov, 1997a) ( Fig. 7D: isl) with a single seta, densely fringed by delicate setules. On inner limb edge, three soft setae. A bunch of long setules is located near these setae. Two robust ejector hooks (Fig. 7D: ejh) strongly different in size, armed with short denticles. The maxillar process (Fig. 7D: mxp), a derivative of gnathobase I, with a single long, densely setulated seta, at base of the limb.
Limb II relatively small, with epipodite supplied by a finger-like projection. Inner limb portion with an anterior row of 6 setae (homologs of "scrapers" of the chydorids, see Fryer, 1968) ( Fig. 7E: 1-6) and disjuncted posterior row two setae (Fig. 7E: pos): a seta near gnathobase and another one near the proximal end of the limb. Gnathobase II with distal armature (Fig. 7E: dag) of three setae of different armature. Filter plate consists of five long setulated setae (Fig. 7E: fpl).
Juvenile female. Instar I has a dorsal head pore (Kotov, 1997b). Body elongated, head relatively high, elevating over valves, without a cervical incision (Fig. 5E). Carapace with a short posterior spine and a long postero-ventral mucro, supplied with minute denticles. Antenna I relatively longer than in female. Free and fused parts of antennules, mucro, rostrum, ventral valve edge, base of caudal spine covered with relatively small spinules.
Ephippial female. Only dorsal portion of valves modified as compared to parthenogenetic female (Figs. 5B-5D, 6J-6K). Ephippium yellowish, ovoid, not clearly demarked from ventral and lateral portions of valves. Egg chamber with a single egg, elongated, its sculpture represented by large polygons well visible under light microscope with very clearly, minute wrinkles and tubercles in each polygon. A strong medial keel on dorsum, strong paired lateral keels well distinguishable from the dorsal view. From the dorsal view, keels projected laterally out of body contour.  (Figs. 8I-8J). Antennae I fused to rostrum, but their bases are not fused together (Fig. 9D). Limb I with a short, thick copulatory hook (Figs. 8K-8M).
Adult male. Shape significantly different from that in female, body short (body height/ length ratio about 0.65), dorsal contour of head humped, dorsal contour of carapace straight, valve anterior portion with few setae anteriorly, ventral margin convex, with setules as in female (Fig. 9A). Head large, with a smooth rostrum and expressed ocular dome, compound eye large (Figs. 9B-9C). Valve armature as in female (Fig. 9D), but a short mucro always present postero-ventrally (Fig. 9E). Postabdomen similar with that in female, its ventral margin slightly comvex, preanal margin slightly to moderately concave. Anal margin almost straight, postanal angle absent. Postabdominal claw bears a basal spine comparable in size with the latter (Fig. 9F), both claw and basal spine slightly and regularly bent.
Antenna I free, remarkably curved distally (Fig. 9G). Frontal sensory seta long, located at middle of antennular body, a short male seta somewhat anteriorly to that, several fields of short spinules located at antenna I anterior face. Long aesthetascs located subterminally, two of them are located on the tip of antenna I, the others located on its lateral surface in two rows. Antenna II with apical and lateral setae as in female. A relatively long (somewhat longer than exopod itself), curved at tip additional male seta on endopod apical segment in position of a rudimentary spine in female (Fig. 9H). Limb I with outer distal lobe bearing two setae strongly unequal in size, copulatory hook relatively large and regularly thick, its tip blunt, not expanded bearing small denticles (Fig. 9I).  (5) additional seta on apical segment of male antenna II curved at tip. Morphological differences from other taxa revealed above genetically are not studied yet.
Distribution and ecology. Widely distributed in the Neotropical zone. Records from Mexico (Elias-Gutierrez et al., 2008) and Central America (Collado, Fernando & Sephton, 1984) need to be checked as they could belong to B. deitersi s.str. or the poorly described B. birgei.
Type material. Lost.
The opinion of Manujlova (Manujlova) that "in the USSR it was not found east to Ob' River" was based on inadequate knowledge of previous literature, i.e. (Pirozhnikov, 1937).
It is widely distributed in Korea (Cho & Mizuno, 1977) (also see descriptions above) and Japan (Ueno, 1937b) (and our data), and present in South China and Vietnam (our data). Most probably, it is present on the Pacific coast of Asia from the Amur to the Mekong basins.

Old Mesozoic group
Our results are consistent with the hypothesis that B. deitersi is, in fact, a species group. We find no evidence for the nominate species in the Palearctic. However, we did find evidence for a genetically divergent and morphologically differentiated Old World lineage. Notably, the strong genetic divergences that we observed and our ancient age estimates were unaccompanied by strong morphological divergences. With our integrated approach, we hoped to mitigate some of the limitations of molecular datasets. As coalescent analyses can oversplit taxa, multigene data benefit from morphological and ecological information. Single gene datasets may disagree with one another and with the species tree for manifold reasons (Fisher- Reid & Wiens, 2011;Hailer et al., 2012). In the present analysis, the topological disagreements (e.g. subclades) were weakly supported, indicating that random error may play a role.
A mature biogeography is only possible with an understanding of timescale (Rosen, 1978). The antiquity of cladocerans of different ranks, from genera to orders, has been confirmed by the fossil record, i.e., from the Mesozoic (Smirnov, 1992;Kotov & Korovchinsky, 2006;Kotov & Taylor, 2011;Liao et al., 2020). Unfortunately, the Palaeozoic records (Anderson, Crighton & Hass, 2004;Womack et al., 2012) are dubious: the described animals could belong to the Cladocera, but also could be members of other crustacean groups (Van Damme & Kotov, 2016). Kotov & Taylor (2011) demonstrated that extant genera of the Daphniidae and even the subgenera of the genus Daphnia existed at the Jurassic/Cretaceous boundary, ca. 145 MYA. More fossil calibrations are possible for the group.
Efforts to use fossil calibrations with molecular data have been limited for Cladocera (Sacherova & Hebert, 2003;Schwentner et al., 2013;Cornetti et al., 2019). Perhaps the only known calibration point for relaxed molecular clocks is the Daphnia/Ctenodaphnia/ Simocephalus split at 145 MYA (Kotov & Taylor, 2011). Non-calibrated molecular clocks also suggest earlier differentiation of the cladocerans, i.e., differentiation of the subfamilies within Chydoridae in the mid-Palaeozoic (Sacherova & Hebert, 2003). A fast molecular clock estimate gave a divergence time for Daphnia at more than 66 MYA (Cornetti et al., Figure 15 Schematic representation of distribution of two major morphotypes of the juvenile parthenogenetic female: with several mucro-like spines (red) and a single mucro (blue). Visualisation was made in free software DIVA-GIS7.5.0 https://www.diva-gis.org using free spatial GIS data from http://www.naturalearthdata.com as the layers. Symbols are inserted manually.
Full-size  DOI: 10.7717/peerj.11310/ fig-15 2019). This value is probably too young given the calibration point of 145 MYA. A more realistic estimate should exceed the minimum fossil calibration (Kotov, 2013). The Family Bosminidae contains only the genus Bosmina and the genus Bosminopsis. Our very rough estimation (see Fig. 4A) suggests that the Bosminidae could be even older than the Daphniidae. Such a conclusion agrees with the hypothesis that bosminids are a sister group to Chydoridae (Kotov, 2013). Chydorids are probably of Palaeozoic origin (Sacherova & Hebert, 2003). No Mesozoic bosminids are known to date. Bosmina was one of the first genera to be studied with paleolimnology. Unlike Bosmina, subfossil remains for Bosminopsis are unknown from the Holocene and Pleistocene bottom sediments (Austin, 1942;Hofmann, 1984). It seems unlikely that a detailed fossil record will be found for Bosminopsis. So molecular clocks are the only method to estimate the time of its differentiation. We estimate that the differentiation of the main Bosminopsis lineages took place in the Cretaceous, and coincided with the disruption of Pangaea, or later disruption of Laurasia. Mesozoic lineages survived in SE Asia and elsewhere in Eurasia (the exact location is unclear) after a mass extinction during the mid-Caenozoic (Korovchinsky, 2006;Van Damme & Kotov, 2016). Most probably, the Pacific Coast Region ("Far East") was the center of B. zernowi diversificaton, as this region is the richest in mitochondrial haplotypes (Fig. 3), and is often a center of diversity for cladoceran taxa (Kotov et al., 2021).
While there is strong genetic divergence between New World and Old World lineages, a more detailed assessment of the divergence time awaits further geographic and genomic sampling for the New World. Within B. zernowi, our results suggest a mitochondrial differentiation in the mid-late Caenozoic (or even Quaternary), but the divergence of its mitochondrial haplotypes was weak.
Korovchinsky (2006) postulated that extant cladocerans are relicts of a mass extinction in the mid-late Caenozoic. For Cladocera, Pleistocene mass extinction in the Holarctic due to glaciation and aridization (Hewitt, 2000) also has phylogeographic support (Taylor, Finston & Hebert, 1998;Cox & Hebert, 2001). But phylogeographic publications referring to previous epochs with non-Holarctic samples are rare (Xu et al., 2009;Kotov et al., 2021). For Bosminopsis, our results suggest a Mesozoic differentiation of the lineages and then survival of only two main lineages in the mid-Caenozoic. We failed to detect divergences consistent with the Quaternary. Our results are consistent with contintental endemism and longterm morphological stasis.
Morphological divergence in Bosminopsis appears to be weak since the Mesozoic. This divergence involves fine-scaled characters such as the mucro-like spine number, male basal spine, and antenna I appearance and armature. Such subtle differences among species are known in other cladoceran groups (Kotov et al., 2021) but can only rarely be associated with a timescale.
There are no known fossil records for the globally distributed B. deitersi group. However, Bosminopsis may be a "living fossil" sensu Darwin (Darwin, 1859). The B. deitersi group has survived with very little morphological change since the Mesozoic despite profound abiotic and biotic changes to the continental water bodies over this timescale. Our results indicate that the occupation of differing climates has also left a weak morphological signature. While the concept of "living fossil" is somewhat ambiguous (Casane & Laurenti, 2013) there are several groups that appear to have undergone morphological stasis since the Mesozoic. Our evidence is consistent with Frey (1962) who expected stasis to account for continental endemism in Cladocerans.

Field Study Permissions
The following information was supplied relating to field study approvals (i.e., approving body and any reference numbers): Field collection in public property in Russia does not require permissions. Samples in South Korea were collected in the frame of cooperation between A.A. Kotov and the National Institute of Biological Resources of Korea and does not require special permission. The sample from Arkansas, USA was obtained from collections resulting from a previous NSF grant. The samples from Japan, China, and Thailand were provided by our colleagues having permissions to collect them due to their scientific activity in the governmental institutes in the corresponding countries. Formol samples from Brazil were kept in the Collection of Zoological Museum of Moscow State University for a long time, they were collected before the time when Brazil introduced very strict regulations for sampling. The species were not assessed as endangered at the time of collection and are currently not subject to specific regulations, however all efforts were taken to ensure that the collection and preservation of animals was performed with due consideration of their welfare. The number of individuals taken did not represent a significant proportion of the population present at each site.

DNA Deposition
The following information was supplied regarding the deposition of DNA sequences: All data generated or analysed during this study are included in Open Science Framework project (https://osf.io/4fjnm/) and this published article. All sequences are deposited at the NCBI GenBank accs. no. MT757174-MT757274, MT757314-MT757388, MT757459-MT757473. Specimen series from which DNA was extracted are deposited at Zoological Museum of Moscow State University.

Data Availability
The following information was supplied regarding data availability: Samples

Supplemental Information
Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.11310#supplemental-information.