One More “Living Fossil”: Mesozoic Lineage Differentiation, Morphological Stasis and Weak Phylogeographic Structure in Eurasia in the Bosminopsis Deitersi Group (Crustacea: Cladocera: Bosminidae)


 Background: 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. Results: Here we use an integrated morphological and multilocus genetic approach to 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. Conclusions: 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).

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 (Supplement Table 2). N -number of sequences; n -total number of sites (excluding sites with gaps or missing data); Snumber of segregating (polymorphic) sites; Hd -haplotype diversity; h -number of haplotypes; Pinucleotide diversity per site; k -average number of nucleotide differences; Fs -Fu's neutrality statistic [49]; D -Tajima D neutrality test [50], the star-sign is indicated statistical signi cance P < 0.05; DemMod -most likely demographic model by DnaSP v.6 coalescent simulation (PBpopulation bottleneck, PD -population decline, n/d -not de ned) based on Theta-W, probability P(Sim < Obs) is indicated as a percentage in brackets.
Three major clades of Bosminopsis deitersi group were revealed from a tree based on the mitochondrial dataset (Supplement Fig. 1A). The rst clade was B. zernowi -widely 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 Bosminopsis sp. (from Thailand) absent from the tree based on nuclear genes. The large clade of B. zernowi was again subdivided into two subclades (1 and 2) with low support. There were some con icts 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 Supplement Fig. 1A and 1B). 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 re ecting historical processes.
The nal tree based on combined mitochondrial and nuclear datasets (Supplement Fig. 1C) was fully congruent with the mitochondrial tree -major clades are well-supported. No con icts were found between ML and BI (with unlinked data) trees for separate genes, mitochondrial and nuclear trees or the nal consensus tree.
The results of all 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. All approaches suggested an independent status of B. zernowi, B. deitersi and Bosminopsis sp. from Thailand. Only the nuclear tree suggested a "Far Eastern" sub-clade based on the information in the hypervariable domains D1 and D2.
There was some evidence that North American and South American populations of B. deitersi form independent species. More sampling of North American populations and loci is warranted to test this hypothesis further. Compared to bGMYC and mPTP, STACEY suggested signi cantly more taxa. However, increased splitting is expected (compared to morphological evidence) with STACEY [51,52]. 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, 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 [53,54]. Again, North and South American populations may be separate species, but more sequences are needed to test this hypothesis.  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 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 xation of the substitutions in mitochondrial and nuclear genomes [56]. Minor heterochrony is not an insurmountable obstacle for phylogenetic reconstructions [57].  Fig. 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 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 signi cantly 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.    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 rst and second segment) setae subequal in size, covered by ne setules.
Limb I large, its corm conically narrowing distally. Epipodite ( Fig. 7D: epp) with two long nger-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 V (Fig. 7I) with a small, ovoid setulated preepipodite and an epipodite supplied with a long ngerlike projection. Exopodite with ve soft setae covered by long setules. The distalmost portion of limv as a densely setulated at lobe, two soft setae near it, two setulated setae of subequal length near gnathobase. Filter plate with two long setae.
Juvenile female. Instar I has a dorsal head pore [46]. Body elongated, head relatively high, elevating over valves, without a cervical incision (Fig. 5E). Carapace with a short posterior spine and a long posterventral 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 modi ed as compared to parthenogenetic female   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 elds 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).   (Fig. 10H, Fig. 13) as in previous species.
Adult male with a dorsal contour of head humped (Fig. 14A); head large ( Fig. 14B-C), with a smooth rostrum and expressed ocular dome, a series of mucro-like spines always present postero-ventrally. Valve as in female (Fig. 14D-E). Basal spine on postabdominal claw shorter than the claw itself (Fig. 14F).
Antenna I free, its distal portion only slightly bent (Fig. 14C, G). A long additional male seta on endopod not curved at tip (Fig. 14H). Limb I with a copulatory hook (Fig. 14I) relatively more massive that in previous species.
Size. 0.25-0.47 mm.   Amur River basin [97]. The opinion of Manujlova [37] that "in the USSR it was not found east to Ob' River" was based on inadequate knowledge of previous literature, i.e. [27]. It is widely distributed in Korea [98] (also see descriptions above) and Japan [24] (and our data), and present in South China and Vietnam (our data). Most probably, it is present on the Paci c coast of Asia from the Amur to the Mekong basins.
Most probably, the tropical countries of Asia are fully populated by other taxa, as in all illustrations the females have a single strong mucro [39,103,104]. Also, a single mucro is illustrated in the gures of Bosminopsis from North America [22], Africa [42] and Australia [38] (Fig. 15). 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.

Old Mesozoic group
Our results are consistent with the hypothesis that B. deitersi is, in fact, a species group. We nd no evidence for the nominate species in the Palearctic. However, we did nd 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 bene t from morphological and ecological information. Single gene datasets may disagree with one another and with the species tree for manifold reasons [105,106]. 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 [107]. The antiquity of cladocerans of different ranks, from genera to orders, has been con rmed by the fossil record, i.e., from the Mesozoic [108 -112]. Unfortunately, the Palaeozoic records [113,114] are dubious: the described animals could belong to the Cladocera, but also could be members of other crustacean groups [4]. Kotov and Taylor [111] 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 [115][116][117]. Perhaps the only known calibration point for relaxed molecular clocks is the Daphnia/Ctenodaphnia/Simocephalus split at 145 MYA [111]. Non-calibrated molecular clocks also suggest earlier differentiation of the cladocerans, i.e., differentiation of the subfamilies within Chydoridae in the mid-Palaeozoic [115]. A fast molecular clock estimate gave a divergence time for Daphnia at more than 66 MYA [117]. This value is probably too young given the calibration point of 145 MYA. A more realistic estimate should exceed the minimum fossil calibration [118].
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 [118]. Chydorids are probably of Palaeozoic origin [115]. No Mesozoic bosminids are known to date. Bosmina was one of the rst genera to be studied with paleolimnology. Unlike Bosmina, subfossil remains for Bosminopsis are unknown from the Holocene and Pleistocene bottom sediments [119,120]. 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 [4,121]. Most probably, the Paci c Coast Region ("Far East") was the center of B. zernowi diversi caton, as this region is the richest in mitochondrial haplotypes (Fig. 3), and is often a center of diversity for cladoceran taxa [54].
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 [121] 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 [122] also has phylogeographic support [123,124]. But phylogeographic publications referring to previous epochs with non-Holarctic samples are rare [7,54]. 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 Bosmniposis appears to be weak since the Mesozoic. This divergence involves ne-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 [54] 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 [125]. 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 [126] 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.

Preliminary comments on further taxonomic revision
Further studies are needed to demonstrate that the North American populations form a separate species from South American specimens. If so, then the taxonomic name for North American specimens would be Bosminopsis birgei Burckhardt, 1924. Records of Bosminopsis are infrequent in North America and are mainly from the southeastern USA [30,127]. Recent biotic exchange between North America and South America has occurred for several cladoceran genera [128]. In such cases, there is very little genetic differentiation for mitochondrial markers among continents. The status of populations from East Asia [39] must also be addressed, including that of B. devendrari Rane (a possible proper name for the SE taxon recorded above). To date, we have no information on Bosminopsis cf. deitersi from Africa, i.e., described by Daday [129] from Lake Nyassa as Bosminella Anisitsi var. africana with a single mucro. This taxon is found in different African countries [12,14,42]. The status of Australian populations [38], also having single mucro, remains unclear.
Subtle geographic variation in the morphology of Bosminopsis has been known since the early 1900's. Meissner [71] and Klocke [78], for example, pointed out the numerous stout spines at the postero-ventral corner of the valves from Russian and Japanese Bosminopsis -the most prominent difference between B. deitersi and B. zernowi. Still, subsequent authors failed to recognize this variation as taxonomically valuable. Klocke[78] concluded that there are two species in Japan, B. ishikawai and B. deitersi. He concluded that the former has stronger denticles on antenna I, better-developed reticulation, a posterodorsal projection located more dorsally and longer spines at postero-dorsal angle "making it similar to Ilyocryptus". In reality, all listed differences are characteristic of juvenile females. Therefore, Klocke[78] erroneously regarded the populations with large adults and without large adults those as two separate species. The same mistake was made by Rey

Conclusions
Here we revised only populations from Eurasia. Other taxa discussed above in the genetic section must be reconsidered in the future, and biological differences must be studied in detail. Each of these putative species has a single mucro at a postero-dorsal angle and minimal differences between parthenogenetic females. We expect that comparing males will be the most fruitful for assessing morphological diagnoses, as male morphology tends to differ among species more than female morphology in cladocerans [130,131].

Collecting samples and morpholohical analysis
Samples from different localities were collected via small-sized plankton nets (with mesh size 50 µm) and xed via 4% formaldehyde or 96% ethanol in the elds, using conventional techniques. All samples were initially examined using a stereoscopic microscope LOMO. Individuals of Bosminopsis were initially identi ed via available references using morphological features [46,65,132]. Existing museum samples were used for morphological analysis (see the list of material in Supplement Table 1).
The morphology of populations from the Neotropics and the Palaearctic was examined in detail to assess the presence of taxonomically signi cant characters. Specimens of Bosminopsis from presorted 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, United Kingdom), and examined under a scanning electron microscope (Vega 3 Tescan Scanning Electron Microscope, TESCAN, Czech Republic).

DNA extraction, ampli cation and sequencing
Each specimen was accurately identi ed by morphological characters (Supplement Table 1). Genomic DNA was extracted from single adult females using the Wizard Genomic DNA Puri cation Kit (Promega Corp., Madison, WI) and QuickExtract DNA Extraction Solution (Epicentre by part Illumina, Inc., Madison, WI) using manufacturer's protocols. Two mitochondrial and two nuclear markers were investigated here: (1) the 5'-fragment of the rst subunit of mitochondrial cytochrome oxidase (COI) -a protein-coding marker widely used in DNA barcoding [133]; (2) the 5'-fragment of the mitochondrial 12S rRNA gene (12S) with a structure demonstrating an alternation of highly conservative duplexes and variable loops [134]; and (3-4) 5'-fragments of the nuclear ribosomal genes (18S and 28S). Each fragment contains both long conservative portions and two variable domains. Although these nuclear markers are predominantly used for divergent taxa [135], they are informative at the species level for many microcrustaceans [136].
Primers used for ampli cation are listed in Table 4 PCR products were visualized in a 1.5% agarose gel stained with ethidium bromide and puri ed by QIAquick Spin Columns (Qiagen Inc., Valencia, CA). The PCR program included a pre-heating of 3 min at 94°С; 40 cycles (initial denaturation of 30 sec at 94°С, annealing of 40 sec at a speci c temperature, an extension of 80 sec at 72°С); and a nal extension of 5 min at 72°С (Table 4). Each PCR product was sequenced bi-directionally on the ABI 3730 DNA Analyzer (Applied Biosystems) using the ABI PRISM BigDye Terminator v.3.1 kit at the Syntol Co, Moscow. The authenticity of the sequences was veri ed by BLAST comparisons with published cladoceran sequences in mBLAST [137].
The sequences from this study were submitted to the NCBI GenBank database for 16S acc. no. MT757174-MT757231, for COI acc. no. MT757459-MT757473, for 18S acc. no. MT757232-MT757274 and for 28S acc. no. MT757314-MT757388. D2b2 CGTACTATTGAACTCTCTCTT 48-50 [139] Population analysis, alignment and phylogenetic analysis Nucleotide diversity analysis [140] and neutrality tests were carried out using DnaSP v.6.12 [141]. We applied the Fs [49] and D [50] tests to con rm neutrality and describe demographic processes in Bosminopsis population [142,143]. To determine the most probable demographic model for Bosminopsis sequences, we performed a coalescent simulation for each locus (1000 replications) in DnaSP v.6.12 [141] using ve demographic models (Standard Neutral Model, Population Growth, Population Decline, Population Bottleneck, Population Split/Admixture). The best model was selected based on the Theta-W (theta with Watterson) estimator [140,144].
Alignment of sequences from each locus was carried out using the MAFFT v.7 algorithm [145] available on the server of the Computational Biology Research Center, Japan (http://mafft.cbrc.jp). For the proteincoding gene COI, we used the "Translation Align" option 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 [146].
The best-tting models of the nucleotide substitutions for each locus and for linked data were selected  [149]. 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, 3th). In our opinion, independent phylogenetic analysis for each locus is more correct as compared to uniting of different loci based on general models of the nucleotide substitutions. Mechanical lumping leads in some average probabilities on substitutions bases on theoretical models instead of empirical facts. We believe that quality of phylogenetic reconstruction is higher in case of using of the information from all the trees, We used the multi-taxon coalescence model "star" in BEAST2 [150] and an partitioned models method using numerous models [151]. 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 did not exclude sequences with incomplete and missed data as such exclusion can strongly reduce the accurateness of phylogenetic reconstruction [152].
We used the IQ-TREE v. Cybertaxonomic species delimitation based on DNA data Our approach to cybertaxonomic taxon delimitation was described in a previous paper [54]. 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 [161]. Analysis of Multi-rate Poisson Tree Processes (mPTP) [166], which is most useful for datasets with small genetic distances [167] albeit prone to "splitting" [52], 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 Classi cation Estimation, Yarely) [51], was made in BEAST2. Genealogical relationships were estimated by STACEY with four independent generations (50M generations of MCMC, sampling of every 10 k generation) after incorporating the suggestions from an initial run. STACEY log les were examined with  [53,116,169,170]. Apparently this is a very rough estimation [171]. An alternative is a stochastic approach based on the "calibration" of the molecular clocks by paleontological data and subsequent strict proof mathematical analyses.
To estimate the probability of molecular clock-like data, we applied a Maximum Likelihood method implemented in MEGA-X v. 10 with a Yule speciation model as the most proper for datasets with several potential species [173]. Four independent runs of 50M generations were done, with each 100 k tree sampled. Subsequent analysis was performed as above for BI following the recommendations [174].

Phylogeographic reconstructions
To test phylogeographic models, we used the packet BioGeoBEARS [175] with the integrated statistical package of the "R" language in RASP4 v.4.2 [176]. The data set was composed based on the maximum number of geographic localities and representation of all phylogenetic lineages revealed by cybertaxonomic methods. We reconstructed a phylogenetic tree based on sequences of two mitochondrial genes (COI and 16S) for all main locations in the distribution range of Bosminopsis cf. deitersi. A mitochondrial phylogeny has some advantages, such as the absence of recombination.
Objective software limitations allowed us to analyze only 27 sequences from ve 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 [177]. 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 10 k generation Ethics approval and consent to participate

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. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.