Multiple Recent Colonizations of the Australian Region by the Chydorus sphaericus Group (Crustacea: Cladocera)

Abstract

Biotic introductions are an ongoing disruption for many ecosystems. For passively dispersed freshwater zooplankton, transcontinental introductions have been common but are poorly studied in the southern hemisphere. Here we assess the hypothesis of recent introduction for populations of the Chydorus sphaericus group (Crustacea: Cladocera) in Australia. We analyzed 254 sequences (63 original sequences) from the cytochrome oxidase I region of mitochondrial DNA of Chydorus sp., which included global representation. Three Australian populations were connected with separate clades in the northern hemisphere, suggesting multiple colonization events for Australia. The timescale of the divergences was consistent with recent (Quaternary) dispersal. As Australian populations are exposed to migrating birds from the northern hemisphere, both avian and anthropogenic sources are candidates for dispersal vectors. We concluded that recent cross-hemisphere dispersal in the Chydorus sphaericus group is more common than previously believed.

1. Introduction

The rise of phylogeographic methods brought about a paradigm shift in biogeography that deemphasized population genetics while spawning tree-based approaches such as coalescent modeling [1,2]. These approaches permitted advances in reconstructing the history of dispersal at a regional scale [3,4] or more rarely a global scale [5,6]. Large-scale phylogeography was bolstered by a desire to better understand the “peopling” of Earth [7].

The Cladocera (Crustacea: Branchiopoda) played an important role in recent phylogeographic investigations of freshwater microfauna. Cladocerans are a relatively large group containing more than 850 known species belonging to more than 100 genera, most of them belonging to the order Anomopoda Sars [8]. Anomopod diversity is greatly underestimated [9]. It is an intensively studied group, but the majority of recent authors deal with a single genus Daphnia Müller which is regarded as a very important model for genetics, evolutionary biology, and toxicology [10,11]. Many examples of the “continental endemism” in the Cladocera were detected by morphology-based taxonomists [12,13,14,15,16]. These conclusions were greatly reinforced by studies with molecular markers [17,18,19,20,21]. Several papers were published on trans-Eurasian and trans-Holarctic phylogeography [22,23,24,25,26,27], but truly global phylogeographic reconstructions remain rare [19,28,29,30].

It is now widely accepted that cladocerans are a very old group. At least some “recent” genera (and even the subgenera of Daphnia) were already differentiated by the Jurassic/Cretaceous boundary [31], or, most probably, even earlier [32,33]. However, there may have been near cosmopolitan distributions for some taxa before continental drift sundered the supercontinents. At the same time, some preliminary investigations using mitochondrial genes revealed modest differentiation between multiple species from different continents and some transcontinental dispersal scenarios in cladoceran biogeography [34].

Recent anthropogenic introductions have altered ecosystems and have threatened biodiversity. Thousands of scientific publications have demonstrated their harmful economic effect and destructive consequences for terrestrial, marine, and freshwater ecosystems [35,36,37,38,39]. Many intergovernmental documents, such as the Convention on Biological Diversity, warn of biological invasions and their consequences. It is obvious that biological invasions and the mixing of faunas are international problems, as the invasive species do not recognize international borders [38,40]. Some of these introductions are “cryptic”, defined as “the introduction and spread of non-native lineages within the species’ native range (intra-specific cryptic invasion) or the invasion of non-native species that goes unnoticed due to misidentification as a native or another invasive species (inter-specific cryptic invasion)” [41]. Cryptic invasions are often thought to be less harmful to local ecosystems and are often poorly studied [41,42].

Most recorded invasions of cladocerans are less conspicuous than those of large organisms. However, even among cladocerans, some invasions are more conspicuous than others. For example, cladocerans with conspicuous defensive morphology such as Bythotrephes longimanus, Cercopagis pengoi, Daphnia lumholtzi, and Bosmina coregoni [43,44,45,46,47] are among the best-studied cladoceran invasions. Invasive D. curvirostris is less well studied, in part because it is difficult to distinguish by nonexperts from the native fauna. Other invasions in the Cladocera are even more cryptic because a native taxon was replaced by a closely related congener [48,49,50]. Finally, long-distance introductions below the species level have been revealed in some studies with molecular markers (i.e., in the D. galeata, D. curvirostris, D. pulex, and D. magna groups) [20,26,51,52,53]. So far, few studies of nonindigenous populations have been made for anomopod cladocerans beyond the genus Daphnia [54,55,56]. The underrepresentation of non-Daphnia species applies to genetic and morphological investigations [57]. The dearth of studies does not appear to be based on a lack of potential for introduced anomopod populations [29,58]. The increasing records of mitochondrial DNA from cladocerans in public databases such as BOLD and GenBank have accelerated the discovery of introductions in this group [59].

Cladocerans are well-known for passive dispersal of the resting egg stage. Gamogenetic females of the Anomopoda produce ephippia–modified molting exuviae containing resting eggs. Ephippia are well-protected against unfavorable environments and are dispersed readily by wind, water currents, and other natural vectors [60]. Waterbirds, however, represent the main vector of natural dispersal of cladocerans and other microscopic invertebrates [61,62]. The ability of waterbirds to carry and transfer the resting stages of different microscopic animals was the main justification for Darwin’s [63] dispersalist paradigm. Recently, however, humans have become important vectors of resting eggs [46,49,61].

The anomopod, Chydorus sphaericus O.F. Müller, has often been regarded as a cosmopolitan taxon [64], but after the rejection of the “cosmopolitanism” concept [12], a series of taxa, earlier regarded as belonging to this group, were separated [13,65]. Even after these revisions, sphaericus-like populations were found to be very common in different regions of the Holarctic and beyond. Chydorus sphaericus is still rare in the tropics and subtropics of the eastern hemisphere, where other species of the genus are dominant [66,67]. In the eastern hemisphere, the C. sphaericus group has a bipolar distribution: it is the dominant taxon in North Eurasia and relatively common in North China, Korea, and Japan [25,68,69], but remains undetected by recent authors in South China [70], SE Asia [71,72,73], and the Malay Archipelago [74,75]. Chydorus sphaericus-like populations are very rare in the Oriental Zone; previous records of this group mostly belong to other taxa [74]. They are also rare in tropical North Australia, while disjunct populations are common in South Australia and Tasmania [66,76].

C. sphaericus was first reported from Australia in the 19th century. King [77] described two species of the genus Chydorus from Australia: C. angustus King and C. leonardi King. Sars [78] concluded that C. angustus was “related perhaps to the European species C. globosus Baird”. However, C. globosus now belongs to another genus, Pseudochydorus Fryer, [60,75]. Sars also stated that C. leonardi was “very near to the well-known European species C. sphaericus (Müller)”. However, the resolution of species identification was very low at that time, and almost all representatives of the genus Chydorus were assigned to C. sphaericus until David Frey [13,79] and Nikolai Smirnov [64,80]. The monograph of Smirnov & Timms [66] verified the presence of C. sphaericus in Australia. Moreover, the species was found to be common and abundant in the southern portion of this continent (including Tasmania) while very rare in tropical Australia. However, molecular methods have shown that even the well-studied Holarctic populations of C. sphaericus can represent a group of several species with local distributions [22,25,81]. Diagnostic morphological characters have been proposed for these regional taxa [22,65]. Sharma & Kotov [82] noticed that the mtCOI sequences obtained from two reservoirs in South Australia were very similar to those obtained from Iceland and Greenland. They concluded that Australian populations were established as a result of an anthropogenic introduction, most probably associated with the mass stocking of fishes from Europe. However, the endemicity status and clade associations of other populations of C. sphaericus in Australia remain unstudied.

In this study, we analyze populations of the Chydorus sphaericus group from Pacific Asia, SE Asia, and Australia with the aim of testing the question of indigenous/nonindigenous status and reconstructing the history of colonization of Australia. We note that a taxonomic revision of the C. sphaericus complex is due but beyond the scope of the present study.

2. Materials and Methods

2.1. Field Collection and Preliminary Morphological Analysis

Samples were collected by plankton nets with a diameter of 0.20–0.4 m and a mesh size of 30–250 µm and immediately preserved in 70–96% alcohol in the field (90–96% alcohol after laboratory sorting). Before the start of genetic studies, a few specimens of Chydorus from each locality were preliminarily identified by morphological characters. Only members of the C. sphaericus group were incorporated into this study. They were identified based on the following characters of the parthenogenetic females: (1) subglobular shape of body; (2) elongated labral keel with an acute tip; (3) absence of denticles at posteroventral portion of valve; (4) aesthetascs located only terminally on tip of antenna I; (5) sensory seta of antenna I located in its middle; (6) antenna II with setal formula: 0–0–3/0–1–3; (7) inner distal lobe of limb I bearing the largest seta additionally chitinised (Figure 1).

2.2. Genetic Analysis

DNA extraction was performed from the specimens accurately identified by morphological characters. Genomic DNA was extracted from single adult parthenogenetic females using the Wizard Genomic DNA Purification Kit (Promega Corp., Madison, WI, USA). We used as a genetic marker the 5′-fragment of the first subunit of mitochondrial cytochrome oxidase (COI)—a protein-coding marker widely used in routine DNA barcoding [83] and monitoring of biological invasions [84]. Amplification was performed using general primers Chy_COI-F1 5′-AAT GTA ATT GTC ACT GCA CAY GC-3′ and Chy_COI-R1 5′-AAA AGT ATG GTA ATA GCY CCT GC-3′, with a PCR-product size of about 0.5 kbp. For partially degraded (i.e., improperly fixed and kept in collections for extended periods) samples, we used the internal primers for a conservative portion of the COI fragment: Chy_COI-Fin 5′-GCK GGY ACA GGT TGA ACT GTT TA-3′ and Chy_COI-Rin 5′-RAT CAA CTG AAG CTC CTG AAT GA-3′ designed specifically for the C. sphaericus group [25]. A combination of Fin+R1 and F1+Rin primers gave two PCR products of about 0.3 kbp length each, with an overlap zone of about 70 bp. Polymerase chain reactions (PCR) were carried out in a total volume of 20 μL, consisting of 4 μL of genomic DNA solution, 1 μL of each primer (10 mM), 4 μL of double-distilled H₂O, and 10 μL 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 diaGene Spin Columns (Dia-M Co., Moscow, Russia) according to the manufacturer’s protocol. The PCR program included a preheating of 3 min at 94 °С; 40 cycles (initial denaturation of 30 s at 94 °С, annealing of 40 s at 52 °С, an extension of 60 s at 72 °С) and a final extension of 5 min at 72 °С. Each PCR product was sequenced bidirectionally on the ABI 3730 DNA Analyzer (Applied Biosystems, Waltham, MA, USA) using the ABI PRISM BigDye Terminator v.3.1 kit at the Syntol Company, Moscow, Russia. Initial analysis of the chromatograms, formation of contigs and their subsequent editing was performed with the Sanger Reads Editor in Unipro uGENE v.37 [85]. The authenticity of the sequences was verified by BLASTn comparisons with published cladoceran sequences in NCBI BLAST nr/nt database [86].

The alignment was carried out using the MAFFT v.7 algorithm [87] on the server of the Computational Biology Research Center, Tokyo, Japan [88], and the FFT-NS-i strategy. An open reading frame was verified by translation. Substitution models (partitioned by codon position) were selected using ModelFinder v.1.6 [89] at the Center for Integrative Bioinformatics Vienna web-portal, Vienna, Austria [90,91]. Substitution models were selected based on minimal values of Bayesian information criterion, BIC [92]. All selected models demonstrated convergence according to the BIC and AICc information criteria, providing additional evidence for model fit [92].

To reduce the influence of population structure on our results, we a priori subdivided all specimens into large population groups (clades) following Belyaeva & Taylor [22], Kotov et al. [25], and Wang et al. [81]. Nucleotide diversity analysis [93], as well as neutrality tests Fs [94] and Tajima’s D [95], were carried out using DnaSP v.6.12 [96]. Genetic distances were calculated in MEGA-X v.10.2 [97], using the recommended p-distances for barcoding purposes [98]. Standard error estimates were obtained by a bootstrap procedure (100 replicates). Analyses were conducted using the maximum composite likelihood model [99]. This analysis involved 254 nucleotide sequences. All ambiguous positions were removed for each sequence pair (pairwise deletion option). There was a total of 453 positions in the final dataset. Testing of the “isolation-by-distance” model using a Mantel test with 1000 permutations was performed in PASSaGE v.2.0 [100], as this test is found to be useful for detecting nonindigenous specimens in a general pool [58]. A matrix of p-distances for this test was formed in MEGA-X, while a matrix of geographic distances was made in PASSaGE.

Phylogenetic reconstructions based on the maximum likelihood (ML) and Bayesian (BI) methods were made based on “unlinked” data on each nucleotide position in the codon. ML analysis was conducted using IQ-TREE v.2.1 [101] with support calculated using 10k replicas of the bootstrap test in UFboot2 [102]. A Bayesian phylogeny was reconstructed in BEAST2 v.2.6 [103], with all of the parameters of the substitution model determined using BEAUti [104]. In each analysis, we conducted four independent runs of Markov chain monte carlo (MCMC, 50M generations, with selection of each 10k generation), with effectiveness control in Tracer v.1.7 [105]. A consensus tree based on the maximum clade credibility (MCC) was obtained in TreeAnnotator (part of BEAST2) with a burn-in of 25%. Posterior probabilities were used as support values [106]. Because the main clades for BI and ML were congruent, we presented the BI trees, with branch support for key nodes.

For COI clade delimitation, we used the divergence threshold optimizing and clustering approach—locMin [107]. All calculations were performed using the algorithm of [108] in “Microsoft R-Open and MKL” 64-bit v.3.5 [109], as it gives good results for single locus studies. The use of parameter-rich models and many loci has been applied to more complicated species complexes [20].

The general mixed Yule-coalescent (GMYC) model was made to assign individuals to species according to ultrametric time trees derived from single-locus data [110]. We used the Bayesian GMYC model in the ‘bGMYC’ package for “R” [111] to discover species in the MCMC trees. Given that bGMYC is prone to over-split trees containing identical alleles (i.e., zero-length branches) into species [111], we dropped zero-length tips from the MCMC tree before the analysis, then ran ‘bGMYC’ using the multiple-threshold models in the ‘bGMYC’ package with sequences of Chydorus pubescens Sars, 1901, as outgroups. Ultrametric trees for each locus were evaluated in BEAST2 using substitution models, (see the phylogenetic analysis above), a strict clock model, and a Yule process. We sampled 20 million MCMC generations with retention of one tree every 10k generations and a 20% burn-in parameter. For the bGMYC analysis, we randomly selected 100 ultrametric trees after burn-in from BEAST2 trees-file. For each of the 100 trees selected, the analyses consisted of 200k generations with a burn-in of 20%. We set the threshold parameter priors, t1 and t2, to 2 and 100, respectively, and the starting parameter value was set at 100. The results were accepted as statistically significant at a modified p > 0.99 level. This p-value should significantly reduce the likelihood of excessive ‘fragmentation’ in the designation of taxonomic structure.

A single COI gene tree was analyzed applying the multirate Poisson Tree Processes algorithm (mPTP). This software has been shown to be the most accurate method for genetic species delimitations (better than previous PTP and popular distance-based methods) because the mPTP allows each species to have different evolutionary rates. We evaluated species limits with mPTP using the IQ-TREE COI gene tree as input. This analysis was performed on a web server at the Heidelberg Institute for Theoretical Studies, Germany [112]. For the data set, we performed 10 different runs. The convergence of the independent runs was assessed through the average standard deviation of delimitation support values and the overall support for the ML estimate calculated computing the mean of the average support values over the 10 runs, following [113].

For estimation of possible divergence age of different clades, we used the molecular clocks with paleontological calibration [114]. To test the fit of our data to molecular clock models, we used a maximum likelihood test in MEGA-X [97,115]. Nucleotide substitution parameters (using a maximum likelihood substitution model statistical method) were also made in MEGA-X based on lowest BIC (Bayesian information criterion) scores. The null hypothesis of equal evolutionary rate throughout the tree was not rejected at a 5% significance level.

For determination of the relative rate of substitutions, we used both paleontological information [116] and points based on molecular phylogenetic data [32,117]. 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; plus an additional calibration point within the outgroups: Lepidurus/Triops 122 MYA [32,116,118]. The age of lineage differentiation according to a strict molecular clock model was estimated in BEAST2 with the optimized relaxed clock model in ‘ORC’ package for BEAST2 [119] with a Yule process as most properly describing our dataset [120]. Four independent runs of 50M generations were completed, with each 100k trees sampled. Further analysis was conducted following the recommendations of Barido-Sottani et al. [121].

A MP haplotype network was constructed in popART v.1.7 [122] using new ultrafast TNT v.1.5 [123] for reconstructing initial most parsimonious nine trees.

For a preliminary testing of the phylogeographic models for the clade A3, we analyzed aligned sequences in BioGeoBEARS [124] with the integrated package “R” in the RASP v.4.2 [125]. Program limitations allowed us to take into the analysis 40 sequences from 9 large geographic regions. For this analysis, we separately reconstructed a phylogeny of this clade only in BEAST2 using sequence JN233874 (clade B5) as the outgroup and running four independent MCMC chains of 10M generations with the retention of each thousandth tree. For RASP files, 100 trees were randomly selected from the BEAST2 set through the script in “R” package [126], with the outgroup removed. A condensed tree was obtained with 50% burn-in in BEAST2. Six biogeographic models (standard dispersion-vicariant models and those with a correction to the speciation events (+J)) were evaluated according to the Akaike information criterion, AICc_wt [127]. A DIVALIKE+J model, which permitted new lineage formation as a result of direct colonization by a founder without a transitional widely distributed ancestor [128] had the best fit to our dataset. We reconstructed the clade history based on the DIVALIKE+J model in RASP4 with three possible distribution ranges and no a-priori hypotheses of species structure. We used the RASP v.4b module of the Bayesian inference for discrete areas—BayArea [129] as an alternative and accurate stochastic approach for the description of the distribution-range evolution as a continuous process. All input files and localities for this analysis were the same as for DIVALIKE+J. We performed 10M MCMC generations, selecting each 10k tree, with a burn-in of 20% of the first generations. Four independent runs of BayArea were made to increase the effective sample size; the results of runs were combined and visualized through the “Combine result” tool in RASP4. The results of DIVALIKE+J and BayArea analyses were identical.

3. Results

In total, 254 COI sequences (63 original sequences and 191 from GenBank) of Chydorus sp. were analyzed (Supplementary Materials, Table S1). Among them, 251 sequences belonged to the C. sphaericus complex (clades A+B of [22,25], and three sequences belonged to other species (С. pubescens and C. sp.). The C. sphaericus complex was represented by 228 sequences belonging to the C. sphaericus-group (clade A of [22]) with 23 sequences being assigned to the C. brevilabris-linguilabris group (clade B of [22]) (Figure 2).

Here we describe in detail only the major clade A, namely the C. sphaericus group (Figure 3). It consisted of four large subclades:

(1)

A1 which is found in Europe, Greenland, western Siberia, with a single population in Australia;

(2)

A2 which is present in northern Europe, Greenland, Siberia, northern Far East, and Tibet;

(3)

A3 which has a wide distribution area, from eastern Siberia to northern and southern Far East and northern North America, having some populations in Australia; and

(4)

A4 which is present in northern North America (and is not discussed below as it was represented in our analysis by few populations and is not known from Australia).

Figure 3. Phylogenetic ultrametric tree for Chydorus sphaericus group. Major clades are marked as A and B (black bars), internal subclades separated by locMin, bGMYC, and mPTP algorithms are marked by grey bars. Sequences from Australia are marked by red arrows. All shapes marking subclades are as in Figure 2. Phylogenetic analysis of the mitochondrial data set with relaxed molecular clock estimates based on fossil calibration points. The bars depict the 95% highest probability density (HPD) interval of the estimated divergence times. Statistic support of branches is coded by the color gradient from pink (low) to violet (high). Node supports for key points are: UFboot2 (ML) and posterior probabilities (BI) for first and second digits.

Figure 3. Phylogenetic ultrametric tree for Chydorus sphaericus group. Major clades are marked as A and B (black bars), internal subclades separated by locMin, bGMYC, and mPTP algorithms are marked by grey bars. Sequences from Australia are marked by red arrows. All shapes marking subclades are as in Figure 2. Phylogenetic analysis of the mitochondrial data set with relaxed molecular clock estimates based on fossil calibration points. The bars depict the 95% highest probability density (HPD) interval of the estimated divergence times. Statistic support of branches is coded by the color gradient from pink (low) to violet (high). Node supports for key points are: UFboot2 (ML) and posterior probabilities (BI) for first and second digits.

Each clade A1–A3 was represented by several subclades, including widely (A1_1, A2_1, A3_1) and locally distributed ones. A1_2 is present only in Greenland, Iceland, and Australia. A1_3 is represented by a single population in Korea and a single population in European Russia. A2_2 is present in Yakutia only. A2_3 is found in Kamchatka and Bering Island. A3_2 is present in the Eastern Plain in China (a single sequence is analyzed there, but the subclade was found in two localities by Wang et al. [81]) and a single locality in Australia. A3_3 is represented by a single population in Vietnam.

Therefore, sequences from Australia and neighboring islands belong to the subclades A1_2 and A3 (within the latter, we detected isolated haplotypes belonging to the subclades A3_1 and A3_2). Polymorphism metrics for major clades of the C. sphaericus group are represented in

Table 1. Polymorphism of COI gene in Chydorus major clades and subclades of the clade A.

Group	N	S	p	h	Hd	Pi	k	theta	D	Fs
A1	56	52	35	43	0.977	0.014	6.34	0.024	−1.498 *	−38.5
A2	59	53	27	53	0.996	0.017	8.05	0.025	−1.086 *	−56.84
A3	107	27	22	21	0.711	0.009	2.51	0.021	−1.692 *	−8.55
A4	6	45	2	6	1	0.035	15.93	0.046	−1.446 *	−0.29
A total	228	68	54	92	0.931	0.063	16.31	0.043	0.548	−34.65
B total	26	115	91	17	0.951	0.111	38.95	0.107	−0.162	3.09

N—number of sequences; S—number of segregating (polymorphic) sites; p—number of parsimony informative sites; Hd—haplotype diversity; h—number of haplotypes; Pi—nucleotide diversity per site; k—average number of nucleotide differences; theta—parameter of nucleotide diversity by (ps/a1), see [93]; Fs—Fu’s neutrality statistic [94]; D—Tajima’s D neutrality test [95], the one star-sign is indicated statistical significance p < 0.10.

. Evolutionary divergence estimates are represented in Supplementary Materials, Table S2. A Mantel test rejected the “isolation-by-distance” model (p < 0.05) for both whole major clade A, and subclades A1 and A3.

The MP network for COI haplotype groups (Figure 4) agreed with the close associations of Australian populations to northern populations found by the ML and BI trees. Widely distributed subclades of the major clades A1, A2, A3, namely A1_1, A2_1, and A3_1 are represented by star-shape patterns, which may suggest recent population connections (see discussion by Kotov et al. [25]). All smaller clades (A1_2, A1_3, A2_2, A2_3, A3_2, and A3_3) are well-differentiated in the network.

Our molecular clocks suggested a Late Mesozoic differentiation of the genus Chydorus, Late Mesozoic-Earlier Caenozoic differentiation of the C. sphaericus complex (major clades A + B) and clade A and a Mid-Caenozoic differentiation of the subclades A1, A2, A3 and A4. Note that some Australian sequences are not identical to those from the northern hemisphere. Their appearance in Australia most probably had a Quaternary origin. At the same time, additional estimates based on fixed mutation rates (Supplementary Materials, Table S3) suggest significantly younger clade differentiation, i.e., c.a. 16MYA for A/B separation, c.a. 10 MYA for A1/A2/A3 separation, 0.2–2 MYA for the subclades, and very recent separation of the Australian populations. A similar time of phylogeographic split between Europe and western Asia vs. eastern Asia is known for other cladoceran species [24,26].

Phylogeographic reconstruction for clade A3 of the С. sphaericus group based on an unrooted tree is represented in Supplementary Materials, Figure S1. The result of RASP4 analysis suggested a complicated biogeographic history with a combination of dispersal and vicariant events. Both Bay Area and DIVALIKE analyses are consistent with the origin of Australian populations as dispersal events.

4. Discussion

Our results indicate that cases of intercontinental dispersion of C. sphaericus are more numerous than previously believed, and these cases are recent or relatively recent. The consequences of such invasion for ecosystems are unknown and need to be studied in detail: there is a chance that C. sphaericus replaced some native taxa (as it is very common in the temperate portion of Australia now). Both human and avian vectors may play a role in these connections. Common migratory flyways run approximately north–south, and “bird-mediated transport to tropical islands from tropical latitudes is expected to be less common than dispersal from boreal or temperate regions” [132]. The distribution area of the A3 clade is strongly congruent with zone of “the East Asian-Australasian Flyway (EAAF)” of bird’s migrations [133,134] (see, also Supplementary Materials, Figure S1d, lower panel). Many of the migratory birds make long-distance trips from Pacific Eurasia to Australia and represent potential vectors of the invertebrate resting stage introduction [135,136,137]. In particular, the EAAF migratory wader (or shorebird) species such as the Bar-tailed Godwit (Limosa lapponica), Eastern Curlew (Numenius madagascariensis), and Japanese or Latham’s Snipe (Gallinago hardwickii) [138,139,140] have been recorded on Lord Howe Island [141,142]. Although avian migration is ancient, most recent migratory routes are of “geologically recent origin”, with Pleistocene glaciation promoting their formation [143].

Calibrated molecular clocks suggest a Late Mesozoic-Earlier Caenozoic differentiation of the C. sphaericus complex (major clades A+B) and clade A (C. sphaericus group). Moreover, we can hypothesize a Laurasian origin of the major clades based on the DIVALIKE+J and BayArea analyses. Such a conclusion agrees with the antiquity hypothesis suggested by many recent investigators of the Cladocera [12,116,144] and large branchiopods [145].

Several unusual distribution ranges are revealed among the subclades in this study: Greenland + Iceland + South Australia (A1_2) (like [82]), Pacific Region + South Australia + Tasmania (A3_1), and China Eastern Plain + Lord Howe Island (A3_2). Interpretation of such distribution areas as a result of vicariant events within each subclade could be regarded as technically impossible as the molecular clock calculations suggest a recent (Quaternary) age of such divergences. However, during the Late Caenozoic, Australia was separated from Eurasia and North America for tens of millions of years or more. Therefore, the aforementioned pattern is a result of a series of at least three separate dispersal events in relatively recent times. As such, recent dispersal is the favored hypothesis for the origin of Australian populations by DIVALIKE+J and BayArea analyses.

The “standard” rate for crustaceans and other arthropods usually suggests ~1–2% mutations per million years [146]. Thus, most divergences would have occurred much more recently (~2–4 million years for the major clades A1–A3, and younger age of the subclades, see Supplementary Materials, Table S3), somewhere in the Holarctic. However, this rate differs substantially among animals [26,84,147], and calibrated molecular clocks are perhaps more informative for estimating clade ages. In any case, the molecular clocks yield approximate estimates (see the very wide bars in Figure 3). Note that in each case for the Australian population, the error bars include 0 years, and very recent (human-mediated) colonization of Australia could not be excluded. Indeed, both types of molecular clock agree on this aspect.

Studied populations of the C. sphaericus group exhibited a high level of genetic polymorphism, i.e., a high nucleotide diversity together with a relatively high haplotype diversity (Table 1). At the same time, strong differences between Pi and Theta values suggest the absence of a panmixia together with an imbalance between mutation level and genetic drift [93]. Strong differences between Tajima’s D and Fu’s Fs values could be explained by a strong population heterogeneity. Negative values of the neutrality test point to a high probability of expansion processes, although this may be a feature of this group as well as for other animals [148]. Based on the concept of expansion processes, however, it is easy to explain the negative result of the Mantel test (p < 0.05) for the whole major clade A and subclades A1 and A3 containing Australian sequences. Such a pattern could not have originated as a result of “isolation-by-distance”. We hypothesize that some long-term dispersal events took place in the past in the history of A1 and A3 clades.

Similar conclusions could be made based on a comparison of differentiation levels within the major clade A, including sequences from Australia (Supplementary Materials, Table S3). Apparently, representatives of clades A1 and A3 have colonized Australia independently. The first group (represented by single sequence KC020624 [82]) belongs to the northeastern subclade A1_2. The Australian sequence is relatively well-differentiated (more than 2%) from other sequences from this subclade. In any case, the time of the subclade A1_2 differentiation excludes a vicariant pattern in the appearance of an Australian population. It apparently appears as a result of a dispersal event. All other Australian sequences belong to A3 clade. Most of them belong to a central group of haplotypes of the subclade A3_1 widely distributed in Beringian region. A single sequence belongs to the subclade A3_3 being endemic in China, with minimal differences from a Chinese population (less than 1% difference). Representatives of these two subclades apparently have appeared in Australia relatively recently and independently. However, is their penetration connected with human activity? A final conclusion on this topic requires further study.

We agree with Sharma & Kotov [82] that such an anthropogenic scenario seems to be preferable for the Australian populations from clade A1_2. Indeed, the chance of a natural colonization of Australia by the North Atlantic invaders is minimal [34,82]. The Australian sequence is definitively nested within the cluster of European sequences in the phylogenetic tree. Australian populations of this group demonstrate a relatively low divergence level from the Eurasian populations (Supplementary Materials, Table S2). Moreover, some haplotypes both in Australia and in the North Atlantic could be missed in our study. At the same time, the Pacific Region is relatively well-studied, and we are confident that representatives of the A1_2 clade are rare or absent there. It is difficult to propose a vector of such invasion rather than an anthropogenic one. At the same time, we cannot exclude dispersal of this clade via birds (even earlier than in A3) and subsequent accumulation of unique mutations.

The A1_3 subclade (Moscow + South Korea) is a candidate for recent (anthropogenic) introduction within Eurasia. The exact direction of this invasion is unknown; this is a rare clade represented by a single population in each of two different, well-studied regions. We hypothesize that Korea was a donor of this invasion because the population in Moscow is found in the pond in the Botanical Garden of M.V. Lomonosov Moscow State University. Botanical gardens are frequently populated by invasive species. For example, Daphnia ambigua was described initially from the lake at Kew Gardens in London [149] and then found to be an invader from North America [150].

All southern hemisphere sequences from the A3_1 subclade belong to the Pacific Eurasian-North American haplotype group., Moreover, most Australia-Tasmanian sequences belong to a central haplotype of this subclade (Figure 4). Genetic distances between Australia-Tasmanian and Palaearctic populations of the A3 clade are also small. A haplotype that has been introduced by humans is often a common widely distributed haplotype, as in case of A3_1. Ultimately, the revealed pattern could have both an anthropogenic or an older than anthropogenic = natural (bird-mediated Pleistocene?) origin, although such a scenario was ignored by Sharma & Kotov [82]. The South Australian reservoirs, Hindmarsh Island, and Tasmania are sites where introduced fish and macrophyte taxa are recorded [82,151,152], but no invasive fishes have been found on Lord Howe Island to date [153]; direct anthropogenic invasion in its freshwater seems unlikely.

In the case of the C. sphaericus group, we favor a series of long-distance dispersal events over the “isolation-by-distance” model. However, using the COI dataset, we are unable to determine an exact age of the A3 pattern revealed here. Very recent anthropogenic and natural Pleistocene dispersal events are similar in some aspects. Moreover, the consequences of recent anthropogenic disturbance and the collapse of freshwater communities are similar to those due to the climatic changes in the Pleistocene (cooling, aridification). Previously, several cases of the appearance in the Australian Region of nonindigenous microcrustaceans were interpreted as anthropogenic invasions [154,155,156,157]. As some of these taxa are common in the Pacific Region or the northern hemisphere (on Eurasian or/and North American Pacific Coast), we also cannot be fully certain that such colonization is strictly anthropogenic: dispersal by avian vectors remains a viable hypothesis [158].

5. Conclusions

Our results indicate that Quaternary and even pre-Quaternary movement between Holarctic and Australian populations of C. sphaericus has been more prevalent than previously believed. Introductions can be facilitated by both natural dispersal and anthropogenic processes [159]. For example, the formation of eutrophic ponds, reservoirs, and embayments by humans has almost certainly aided in the dispersal of “weedy” cladocerans [160]. Indeed, C. sphaericus has often been associated with eutrophic habitats [161,162]. As freshwater ecosystems are increasingly threatened, it remains important to understand the sources of recent biotic change and what diversity is lost by introductions. As some cladoceran groups leave a sediment record (and in some cases an ancient resting egg bank), we expect that the challenge of distinguishing anthropogenic from recent passive dispersal will be aided by paleolimnological investigation and high-resolution genetic studies.