Rotifera is a group of microscopic animals (< 2 mm) that encompasses approximately 2 000 species (Segers, 2007). They are globally distributed in freshwater and brackish habitats, from large lakes to small ponds, along with a small number of marine habitats (Wallace, Snell, Ricci, & Nogrady, 2006). Rotifers play a pivotal role in aquatic food webs by transferring energy to higher trophic levels (Xiang, Xi, Zhu, & Xu, 2017). The subclass Monogononta, which includes most species in the phylum, comprises planktonic, littoral and benthic organisms (Segers, 2008). Monogonont rotifers reproduce by cyclical parthenogenesis, alternating periods of asexual reproduction with sexual phases (Wallace, 2002). When sexual reproduction is carried out, they produce resting eggs that are resistant to harsh environmental conditions and are adapted to passive dispersal (Mills, Lunt, & Gómez, 2007). These eggs can also accumulate in sediments providing an archive of the species genetic variability (Gómez, Carvalho, & Lunt, 2000).

As with other zooplankton, the taxonomy of rotifers is understudied (Segers, 2007; Fontaneto, Kaya, Herniou, & Barraclough, 2009). The causes of this are well known: 1) there are few taxonomists that can identify a rotifer reliably, 2) rotifers are typically morphologically austere and as such have few recognizable morphological characters and 3) the small size of rotifers makes it difficult to identify some taxa to the species level, this is particularly the case for littoral and benthic rotifers (Mark Welch & Mark Welch, 2005). Among monogonont rotifers, morphology may not be sufficient to differentiate between species due to 1) limited or not morphological variation among cryptic species (Campillo, Serra, Carmona, & Gómez, 2011; Mills et al., 2017), and 2) phenotypically plastic morphological variation within a single species resulting from cyclomorphosis and inducible defenses (Gilbert, 2001; Wallace, 2002). Phenotypic plasticity is an important characteristic of monogonont rotifers, and different morphological variants have in the past erroneously been described as different subspecies, forms or varieties (Segers & De Smet, 2008). Nevertheless, the taxonomic status of these morphotypes remains unclear for most monogonont rotifers (Segers, 2007).

Molecular results of the best studied monogonont rotifers Brachionus plicatilis and B. calyciflorus have shown that each of these two taxa contains at least 15 genetically divergent lineages (Mills et al., 2017; Papakostas et al., 2016). In the case of the B. plicatilis species complex, some cryptic taxa show little interspecific morphological differentiation. Therefore, it seems that extensive genetic and ecological variation in some rotifers can occur along with an apparent absence of morphological differentiation (Serra, Galiana, & Gómez, 1997; Mills et al., 2017). Cryptic species have also been reported in other monogonont rotifer based on genetic data, for example in Epiphanes senta (Müller, 1773) (Schröder & Walsh, 2007), Lecane bulla (Gosse, 1851) (Walsh, Schröder, Wallace, & Rico-Martinez, 2009), Testudinella clypeata (Müller, 1786) (Leasi et al., 2013), Polyarthra dolichoptera Idelson, 1925 (Obertegger et al., 2014), Synchaeta pectinata (Ehrenberg, 1832) (Kimpel, Gockel, Gerlach, & Bininda-Emonds, 2015), Euchlanis dilatata Ehrenberg, 1832 (Kordbacheh, Garbalena, & Walsh, 2017), Keratella cochlearis Gosse, 1851 (Cieplinski et al., 2017), Limnias melicerta Weisse, 1848 and L. ceratophylli Schrank, 1803 (Kordbacheh, Wallace, & Walsh, 2018) and several more (García-Morales & Elías-Gutiérrez, 2013). In addition some studies have used morphological and ecological analysis as well as crossmating experiments to examine cryptic diversity. For example in the B. plicatilis complex, morphology and molecular analysis and mating experiments were carried out (Suatoni et al., 2006). A same approach was used by Schröder & Walsh (2007) to study the Ephiphanes senta complex. Whereas Papakostas et al. (2016) used morphological, ecological and molecular analyses to studied B. calyciflorus species complex. All these studies provide evidence of how molecular approaches are important tools for understanding the nature of genetic differentiation and are also relevant to rotifer taxonomy.

In the most recent revision of Rotifera, the genus Brachionus (Pallas, 1766) contained approximately 62 taxa (Segers, 2007), most of them predominantly inhabiting warm water. Brachionus is diagnosed by a loricated body that is formed by a dorsal and a ventral plate. Dorsal anterior margin of the lorica can contain six, four or two spines, whereas the ventral anterior margin is generally undulated. Some species have spines on the posterior edges of the lorica or at the foot opening. For the genus a malleate trophi is characteristic (Koste & Shiel, 1987). Brachionus are polymorphic in some species groups. One of these polymorphic species is Brachionus quadridentatus Hermann, 1783 which is typically found in the littoral zone of fresh and brackish water (Koste & Shiel, 1987). B. quadridentatus is a cosmopolitan species reported in all continents including Antarctica (Dumont, 1983; Segers, 2007). This species exhibits broad intraspecific morphological variation associated with spine development, which is reflected in different forms, varieties and synonyms of doubtful taxonomic validity (Koste, 1978; Segers, 2007). Two of these forms are B. quadridentatus f. brevispinus (Ehrenberg, 1832) and B. quadridentatus f. cluniorbicularis (Skorikov, 1894), which are common in waterbodies worldwide (Koste & Shiel, 1987). García-Morales and Elías-Gutiérrez (2013) reported a significant difference in mtDNA sequence in some populations of B. quadridentatus collected from southeast Mexico, suggesting the presence of cryptic species.

In the present study we used sequence variation in the mitochondrial cytochrome c oxidase subunit I (COI) and the nuclear ribosomal 18S genes to explore the diversity of Brachionus quadridentatus from inland waterbodies of Mexico. In particular, we investigated whether different morphotypes in this taxon may represent putative cryptic species, with some remarks on the variation in the lineages distribution patterns.

MATERIALS AND METHODS

Sampling and identification: Samples were collected in 12 locations across Mexico (Fig. 1; Table 1 and Table S1). Geographic data from another five water bodies were obtained from a previous study (see Table 1) (García-Morales & Elías-Gutiérrez, 2013). Our study area covered four biogeographic provinces: the Trans-Mexican Volcanic Belt (TMVB) and Baja California Peninsula (temperate zones), as well as the Gulf of Mexico and Yucatan Peninsula (tropical zones) (Fig. 1; Table S1). Samples were collected using a 50 µm mesh size plankton net and preserved in 96 % ethanol. Brachionus quadridentatus specimens were sorted from the samples under a compound microscope, rinsed with distilled water to remove debris, and transferred into PCR tubes for DNA extraction.

Specimens were identified morphologically using the features of the lorica, in particular its overall shape, the shape of its posterior margin, length and shape of the median anterior spines and the posterior caudal spines. We identified the specimens following the taxonomic keys of Koste (1978) and Koste & Shiel (1987). Specimens from different locations and morphogroups were illustrated using a camera lucida attached to an optical microscope.

DNA extraction, amplification, sequencing, and alignment: DNA extraction and amplification for the COI gene were conducted following the protocol described by García-Morales and Elías-Gutiérrez (2013). Details of the 18S rRNA amplification profile are provided in Table S2. The primers for the COI gene were LCO1490: 5´-GGT CAA CAA ATC ATA AAG ATA TTG G-3´ and HCO2198: 5´-TAA ACT TCA GGG TGA CCA AAA AAT CA-3´ (Folmer, Black, Hoeh, Lutz, & Vrijenhoek, 1994). The primers for the 18S gene were 18SF new: 5´-AGA TTA AGC CAT GCA TGT CT-3´ and 9R: 5´-GAT CCT TCC GCA GGT TCA CCT AC-3´ (Tang et al., 2012). PCR products were sequenced in the High Throughput Genomics Center in Seattle, USA. Obtained sequences were deposited in GenBank under accession numbers XXXX. To increase the geographic range of the samples we included 22 mitochondrial sequences of B. quadridentatus from La Escondida, Vernet, Charco 1, Parque Ecologico and El Rosario locations (Fig. 1; Table 1 and S1), reported by García-Morales and Elías-Gutiérrez (2013), and downloaded from GenBank (accession numbers JX216559-JX216569, JX216571-JX216573, JX216575-JX216580, JX216582-JX216583). COI gene sequences were aligned in MEGA v. 6.0 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013) using ClustalW with default settings, whereas 18S gene sequences were aligned with MAFFT v. 7 using the Q-INS-I algorithm regarded as the optimal strategy for ribosomal markers (Katoh, Asimenos, & Toh, 2009). This last aligment was carried out on the MAFFT webserver http://mafft.cbrc.jp/alignment/server/index.html (accessed June 9, 2017). Final alignments were checked visually.

Species delimitation: We applied five methods of species delimitation for mtDNA, rRNA, and the concatenated datasets:

1. The algorithm of Statistical Parsimony (SP) implemented in TCS 1.2.2 (Clement, Posada, & Crandall, 2000), that was designed to construct intraspecific relationships in haplotype networks, but has also been used to detect candidate species (Jörger, Norenburg, Wilson, & Schrödl, 2012). We used this algorithm applying the 95 % connection probability limit, with gaps coded as missing.

2. The Generalized Mixed Yule-Coalescent model (GMYC) (Pons et al., 2006), that uses a maximum likelihood approach to optimize the threshold identifying the shift in the branching patterns of the gene tree from interspecific branches (Yule model) to intraspecific branches (coalescent). This method has been previously used on rotifers (Fontaneto et al., 2009). For the GMYC model we generated ultrametric trees from the three datasets using BEAST v. 2.1.3 (Bouckaert et al., 2014). The parameters comprised a GTR+G+I substitution model, a relaxed lognormal clock, and a birth-death prior (Tang, Obertegger, Fontaneto, & Barraclough, 2014). The analyses were run with 50 000 000 MCMC for the COI and 18S genes, and 70 000 000 MCMC for the concatenated dataset, sampling every 1 000 generations. Because of the absence of a molecular clock specific to Rotifera, we used calibration clocks for COI of 1.76 % sequence divergence per Myr tested in aquatic invertebrates (Wilke, Schultheiß, & Albrecht, 2009) and 0.02 % per Myr for 18S, both recently applied to rotifers (Tang et al., 2014). The MCMC runs were checked for convergence in Tracer v.1.6 (Rambaut, Suchard, & Drummond, 2013). Trees were combined in TreeAnnotator v.2.1.2 using a maximum credibility tree, with the first 10 % discarded as burn-in. The GMYC model was run through the GMYC webserver, using the multiple threshold option (http://species.h-its.org/gmyc/; accessed Juanary 8, 2017).

3. Bayesian implementation of the Poisson Tree Processes model (PTP) (Zhang, Kapli, Pavlidis, & Stamatakis, 2013), that uses a rooted phylogenetic tree as input and model speciation or branching events in terms of number of substitutions, adding Bayesian support values to delimited species on that tree. We used the ML phylogenetic trees (see below), and the analyses were run with 500 000 MCMC generations on the PTP webserver (http://species.h-its.org/). 4) The Automatic Barcode Gap Discovery (ABGD) model clusters sequences based on genetic distances by detecting differences between intra and interspecific distances. This difference between intra and interspecific distances produce a barcoding gap that can be observed whenever the genetic distance among individuals belonging to the same species is smaller than distance among individuals from different species (Puillandre, Lambert, Brouillet, & Achaz, 2012). Therefore, groups of individuals united by smaller genetic distances were considered to be equivalent to species. ABGD was applied only to the COI gene as this method is based on genetic distances calculated in a single marker, therefore was not used in the concatenated data. The 18S rRNA ABGD did not yield any result due to the poor variation of the gene. ABGD analysis was carried out through its online webserver (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html) using default settings. 5) a multispecies coalescent model was performed using *BEAST v.2.0 (Heled & Drummond, 2010) as implemented in the software package BEAST 2 (Bouckeart et al., 2014). All the specimens for which both COI and 18S rRNA sequences were available were used in the analyses, consequently some specimens show in the COI phylogenetic tree that did not have available 18S sequences were not used. In order to test all well differentiated monophyletic groups found in the phylogenetic results two independent hypothesis were tested: 1) three species, based on the 18S phylogeny and 2) six species based of the monophyletic groups found in the concatenated phylogenetic analyses (since the Charco 1 genetic group is not available for 18S rRNA). Species tree analysis for the two assignments were performed as follows; the best fit models of nucleotide evolution for the model of molecular evolution that best fit the mtDNA (TVM+G) and rRNA (TrN) were used with four gamma categories. We used calibration clocks for COI of 1.76 % sequence divergence per Myr (Wilke et al., 2009) and 0.02 % per Myr for 18S (Tang et al., 2014). For species tree priors and population size models, a Birth-Death process and a constant species tree population size were chosen. Three independent MCMC runs were performed for 500 000 000 generations, sampling every 5 000 generations. Chains convergence was assessed by visualizing the sampled parameter values in Tracer 1.6 (Rambaut et al., 2013). Ten percent of the generations were discarded as burn-in and pooled using the module Log-Combiner in the BEAST 2 package. The maximum clade credibility species tree was obtained by the module Tree-Annotator in the BEAST 2 package.

Phylogenetic relationships: Bayesian inference (BI) and maximum-likelihood (ML) analyses were used to infer phylogenetic relationships among populations of B. quadridentatus, with mitochondrial DNA (mtDNA COI) and nuclear ribosomal DNA (rRNA 18S) analyzed separately and concatenated. The BI and ML analyses were conducted using MrBayes v. 3.1.2 (Ronquist & Huelsenbeck, 2003) and RaxML v. 1.3.1 (Stamatakis, Hoover, & Rougemont, 2008), respectively. For the BI analysis, the jModelTest v. 2.1.1 (Darriba, Taboada, Doallo, & Posada, 2012) was used to identify the model of molecular evolution that best fit the mtDNA (TVM+G), rRNA (TrN), and concatenated (GTR+I+G) data under the Akaike information criterion (AIC). The settings for the BI analyses for each molecular dataset were four simultaneous Markov Chain Monte Carlo (MCMC) runs for 4 000 000 generations for mitochondrial and nuclear data, and 6 000 000 for the concatenated dataset, with trees sampled every 100 generations, and a burn-in of 25 %. Finally, a majority-rule consensus tree was obtained. For ML analyses we used GTR+G (mtDNA data), GTR (rRNA data), and GTR+I+G (concatenated) models as implemented in RaxML v. 1.31. Maximum likelihood analyses were run with 10 000 bootstrap replications.

Genetic diversity and differentiation: We measured genetic diversity and genetic differentiation in order to provide some biogeographic remarks of the B. quadridentatus lineages found using mainly COI gene data. The number of haplotypes, haplotype diversity (Hd), and nucleotide diversity (π) were calculated using DnaSP v. 5.10 (Librado & Rozas, 2009) for mtDNA data for each location, as well as for each lineage discriminated by the phylogenetic analyses of the COI gene, and for two of the three lineages found in the phylogenetic trees for 18S gene. Genetic distances (uncorrected p-values) for the COI gene within and among the lineages discriminated by the phylogenetic analyses for the mitochondrial and the nuclear genes were calculated in MEGA v. 6.0. To assess genetic differentiation among the lineages discriminated by the phylogenetic analyses for COI gene (see results), a hierarchical Analysis of Molecular Variance (AMOVA) was performed in Arlequin v. 3.5 (Excoffier, Smouse, & Quattro, 1992), and its significance was tested using 20 000 permutations. To test the hypothesis that genetic differentiation among these mitochondrial lineages increases with geographic distance, the correlation between the FST distances and geographic distance (km) was obtained using a Mantel test (Mantel, 1967) in Arlequin v. 3.5. Significance of the correlation was determined via 10 000 permutations.

RESULTS

Sequence variation: The final alignment for the mitochondrial COI gene had 635 bp, with 242 polymorphic sites and 233 parsimony informative sites, defining 38 unique haplotypes from 176 sequences (154 de novo and 22 from GenBank). For the nuclear 18S gene 830 bp were obtained, with eight polymorphic sites and seven parsimony informative sites defining seven haplotypes from 125 sequences (all de novo). Most mtDNA haplotypes were unique for a single population, but two haplotypes were shared between two populations that belong to the same watershed (COLVEG1 and VERCH1-1) (Fig. 1; Table S1). Cuitzeo Lake had the highest number of mitochondrial haplotypes (8) followed by La Colorada Lake (6) (Table 1). Cuitzeo was also the site with the highest number of specimens. The highest number of nuclear haplotypes (2) was found in the Cuitzeo, La Zarcita and Ixtlahuaca sites (Table 1).

Species delimitation: The number of putative cryptic species based on COI gene sequences ranged from seven to 13 (Table 2, Fig. S1). The ABGD method discriminated seven, the GMYC method revealed 11 putative species (confidence interval: 5-13, likelihood ratio: 11.15828, LR test: 0.00377), whereas SP and PTP methods delimited 13 putative species (Fig. S1). For 18S rRNA, estimates of putative species ranged from one to four (Table 2, Fig. S1). Statistical Parsimony revealed a single species, whereas the PTP model delimited four putative species. The GMYC model found three putative species (Table 2, Fig. S1) (confidence interval: 1-5, likelihood ratio: 0.23773, LR test: 0.88792). The concatenated dataset detected similar number of putative species as the COI analyses, with exception of those samples in which 18S sequences were not available.

The species tree analyses reveal that the most supported hypothesis is that of 18S rRNA that supports the presence of three putative species with high branch support > 95 % (see Fig. 2). We called these three putative species BqI, BqII and BqIII (see Fig. 2 and Fig. 4). For the six putative species hypothesis (see Fig. 3 and Fig. S2), the branch support was < 75 % for some of the branches (Fig. S2).

Phylogenetic relationships and specimen identification: The mtDNA trees produced by BI and ML were identical for all the analyses. The mitochondrial COI gene retrieved seven major lineages (A-G) which are genetically divergent (> 10 %) (Fig. 3). Lineage A was formed by only one individual from Charco 1 (Fig. 3, see also Fig. S3). Lineage B is represented by 20 samples from four locations (La Colorada, La Vega, Vernet and Parque Ecológico) (Fig. 3). Lineage C includes eight samples from La Zarcita spring, being the sister group with high support of the lineage D. Lineage D includes 43 samples from six locations (Escondida, La Zarcita, Uranden, Magdalena, La Colorada and Yuriria) (Fig. 3). Lineage E form a group composed by 11 specimens in one haplotype from La Joya Maar Lake. Finally, two genetically divergent lineages, F and G are sister clades. The lineage F that includes 58 samples from Cuitzeo and Lerma and the lineage G that is represented by 35 specimens collected in five locations (Atotonilco, Ixtlahuaca, Vernet, El Rosario and Teremendo) (Fig. 3). Uncorrected p distances for mtDNA gene within the seven lineages ranged from 0-3.1 %, whereas, pairwise distances among the seven lineages ranged from 12-18 % (Table S3).

All specimens belonging to the lineages F and G were identified as B. quadridentatus f. cluniorbicularis, and all individuals of lineage B were identified as B. quadridentatus quadridentatus. Lastly, lineages A, C, D, and E were identified as B. quadridentatus f. brevispinus (Fig. 3). In addition, seven clearly morphologically differentiated groups were found. These morphogroups were congruent with the seven lineages found based on COI gene sequences (Fig. 3). Within all the specimens identified as the brevispinus form we found four different morphotypes, which belong to the COI lineages A, C, D and E. Moreover, with the specimens identified as the cluniorbicularis form we discriminated two morphotypes, which are coincident with the COI lineages F and G.

The topologies reconstructed by BI and ML methods for the 18S gene were identical and revealed a polytomy comprising three lineages that we called BqI, BqII and BqIII (see above) (Fig. 4). Percentages of uncorrected distances among the three ribosomal lineages ranged from 0.3-0.5 %. Lineage BqI is formed for all specimens identified as B. quadridentatus f. cluniorbicularis and include the mitochondrial lineages F and G. Lineage BqII was identified as B. quadridentatus f. brevispinus and is concordant with the mitochondrial lineage E. Finally, lineage BqIII encompasses all specimens belonging to lineages B, C, and D in mtDNA analyses. This latter lineage is a mix of specimens identified as B. quadridentatus quadridentatus and B. quadridentatus f. brevispinus (Fig. 4).

In the concatenated analyses, the topologies reconstructed by BI and ML were also identical, revealing the same lineages as the mitochondrial gene analyses (see Fig. 3 and Fig. S4), but with some differences in relationships of the lineages B and E compared to the mtDNA phylogenetic analyses. The lineage A is absent due the lack of sequences for the 18S gene (Fig. S4).

Morphological differences between the seven lineages based on COI gene sequences include the length, width, and shape of lorica, length and shape of the caudal spines, length of median anterior spines, shape of foot opening, and size. Specimens from lineage B are the larger in size than specimens from other lineages. Also lineage B presents median anterior spines and caudal spines longer than the other lineages (Fig. 3, Table S4). Specimens from lineage C are smaller than lineages D and E. Length of the median anterior spines in lineage C are similar to those in lineages D and E, but the caudal spines in lineage C are shorter than the other lineages (Fig. 3, Table S4). In lineage C the body shape is barrel-like with the caudal spines directed outward. Specimens from lineages D and E are similar in size and in the length of the median anterior spines, but differ in body width, with lineage E being wider than lineage D (Fig. 3, Table S4). Caudal spines in lineage D are slightly longer than in lineage E, and their base is narrower than the base of the spines in lineage E (Fig. 3). Lineage D and E also differ in the shape of the foot opening and posterior margin of the lorica. In lineage D the foot opening has two terminal spines of the same size and width, whereas in lineage E the two spines are of different sizes and one is wider. Also in lineage D the posterior margin of the lorica is more concave than in lineage E. Lineage F and G have no caudal spines. Specimens in lineage G are larger and wider than specimens in lineage F. Length of the median anterior spines are similar in lineage F and G. Posterior margin of the lorica in lineage G is rounded in its distal part whereas in lineage F the lorica is rounded laterally and almost straight in its postero distal part (Fig. 3).

Genetic diversity and differentiation: Haplotype diversity (Hd) and nucleotide diversity (π) for most locations were low, with the exception of La Colorada Lake (Hd = 0.810 and π = 0.08757) (Table 1). The Hd and π of the mtDNA lineages B, D and G were high, whereas for the lineages C and F were low (Table S5). Lineage E has only one haplotype in 11 sequences. The haplotype and nucleotide diversity in the 18S rDNA lineages BqI and BqIII were low: Hd = 0.253, π = 0.00050 and Hd = 0.062, π = 0.00016, respectively. AMOVA indicated a high level of genetic differentiation for the seven lineages found in the COI phylogenetic analyses for B. quadridentatus complex (FST = 0.9819, FCT = 0.8602, P < 0.00) (Table S6). The Mantel test did not show correlation between genetic distances and geographic distances among mtDNA lineages (r = 0.0882, p = 0.5846), thus providing no evidence for an isolation by distance pattern in the B. quadridentatus complex.

DISCUSSION

Despite the economic, evolutionary and ecological importance of the Brachionus genus its systematics and taxonomy remains unclear (Segers, 2007). This is true for B. quadridentatus, where its high level of phenotypic variation confuses the taxonomy of the taxon (Segers, 2007). Our study demonstrates that morphology-based taxonomy is not effective at providing an accurate assessment of the diversity of B. quadridentatus. Based upon the results of analysis of the 18S gene at least three putative species exist within B. quadridentatus, however as many as seven species may be present based upon the results of analysis of the COI gene and morphological analysis. Analyses of a larger number of molecular markers as well as breeding and ecological studies will be required to further resolve the taxonomy of B. quadridentatus in this study region. In the case of the B. plicatilis species complex, accompanying the morphology-based taxonomy with the use of two or more molecular markers has been demonstrated to be a more reliable approach to delimit species than use of morphology alone (Fontaneto, Flot, & Tang, 2015; Mills et al., 2017).

Using morphology-based taxonomy, we clustered all the sequenced specimens into three recognized morphological variants: B. quadridentatus quadridentatus (long spines), brevispinus form (short spines) and as cluniorbicularis form (without spines). For the first putative species (BqI) identified by GMYC (Fig. S1) and species tree analyses (BqI lineage in 18S Fig. 2), all the 51 specimens belonging to the group lacking spines. Therefore, these specimens were identified as B. quadridentatus f. cluniorbicularis. This morphotype without spines was exclusive to this lineage (Fig. 4), showing congruence between the genetic and morphological approaches despite the high geographic distance among sampling location. In the case of BqII, that was identified as a putative species in the GMYC and species tree analyses (Fig. 2, Fig. 4 and Fig. S1), all 11 specimens were identified as B. quadridentatus f. brevispinus. However, this morphotype with short spines was also found in some specimens belonging to the BqIII lineage. It is noteworthy that specimens from BqIII had differences in the body shape, length of the anterior and caudal spines and shape of the foot opening. The lineage BqIII in 18S tree (Fig. 4), also delimited as a putative species by the GMYC (Fig. S1) and species tree analyses (Fig. 2), is represented by 52 specimens with short spines, thus classified as B. quadridentatus f. brevispinis and 11 with long spines identified as B. quadridentatus quadridentatus. As such, this putative species is a mixture of morphotypes.

Although the most conservative result found in the 18S gene was three putative species within each of those putative species, there was genetic variation in the COI gene and morphological differentiation. The 18S gene is therefore not suitable to delimit species within the B. quadridentatus species complex, as this gene underestimates diversity (Tang et al., 2012). We found seven morphotypes that display clear diagnostic morphological differences (e. g. length of the median anterior and posterior spines, body shape and shape of the foot opening) (Fig. 3). However, using morphological features alone to differentiate groups has been demonstrated to be problematic. This is especially true of animals that exhibit few morphological attributes as is the case of the B. plicatilis species complex (Mills et al., 2017) and meiofaunal Pontohedyle slugs (Jörger et al., 2012), or that are morphologically plastic (due cyclomorphosis and predators), for instance corals (Ow & Todd, 2010), and the rotifers Euchlanis dilatata, Brachionus havanaensis and Keratella cochlearis (Koste & Shiel, 1989; García-Morales & Elías-Gutiérrez, 2013; Cieplinski et al., 2017). Thus, accompanying morphology-based taxonomy with the use of molecular markers has been demonstrated to be a more reliable approach to identify such diversity (Fontaneto et al., 2015; Mills et al., 2017). Accordingly, in the phylogenetic analysis, ABGD and AMOVA based on COI gene sequences, we also found seven monophyletic lineages that are divergent.

The uncorrected genetic distances calculated among the seven lineages and morphologic variants of the B. quadridentatus species complex are high (12-18 %) in comparison with distances known in other animals (Hebert, Cywinska, Ball, & De Waard, 2003), exceeding the values usually found between congeneric species (Gómez et al., 2002). These genetic distances are similar to the values found between distinct species within the related B. plicatilis (12-23 %, Gómez et al., 2002; 11.9-23.2 %, Mills et al., 2017) and the B. calyciflorus species complexes (9-13 %, Gilbert & Walsh, 2005; 11-22 %, Xiang, Xi, Wen, Zhang, & Ma, 2010; 10.4-26.8 %, Xiang et al., 2011). Moreover, the seven observed morphotypes and monophyletic lineages are congruent. Therefore, we propose the existence of seven putative species within B. quadridentatus from Mexico based on the COI gene and morphological data.

In addition, we compared our COI sequences from Mexico with 23 COI sequences of B. quadridentatus from other parts of the world (Iberian Peninsula, United Kingdom and China) deposited in GenBank (see accessions in Table S7) to determine if they are similar or not. We found that sequences from Iberian Peninsula, United Kingdom and China formed five independent lineages. Genetic distances among these five COI lineages and our seven COI lineages ranged from 13-24 %. Thus increasing the number of putative species in B. quadridentatus. Only five sequences from the Iberian Peninsula reported by Moreno et al. (2017) fell within mitochondrial lineage G of this study (see Fig. S5).

We found evidence for phylogenetic mitonuclear discordance between the mitochondrial COI and the nuclear 18S genes, as the nuclear marker revealed fewer putative species than mitochondrial marker. Other studies of rotifers have also reported varying levels of mitonuclear discordance. COI and ITS1 markers in the B. calyciflorus species complex, shown a high level of discordance (Papakostas et al., 2016), and the Keratella cochlearis, Polyarthra dolichoptera and Synchaeta pectinata species complexes presented intermediate levels of mitonuclear discordance among the same genetic markers (Obertegger et al., 2018). By contrast, in the B. plicatilis species complex a complete mitonuclear concordance between the COI and ITS1 markers was found (Mills et al., 2017). The mitonuclear discrepancies reported in these previous studies have been attributed to processes such as hybridization or incomplete lineage sorting (Papakostas et al., 2016; Obertegger et al., 2018). However, it has also been reported that mitonuclear discordance may be attributed to other factors such as horizontal gene transfer, androgenesis and unresolved phylogenetic polytomies (Obertegger et al., 2018). In our study, phylogenetic analyses conducted on the 18S gene recovered a polytomy. However, we suggest that incomplete lineage sorting is the most likely cause of the mitonuclear discordance observed in our work. This is because differentiation has not been completed in the nuclear gene (18S) and therefore nuclear lineages share haplotypes, retaining ancestral polymorphism. There are lineages living in sympatry (COI lineages: A and G, B and D, B and G) but there is not evidence that hybridization has occurred. Future research will be needed to determine what process is responsible for the mitonuclear discordance observed in the B. quadridentatus species complex by the sequencing of a greater number of molecular markers.

Currently two explanations for the geographic patterns observed in zooplankton species have been outlined in the literature. Firstly, the high dispersal capacity of resting eggs allows high levels of gene flow across large areas resulting in widely separated areas with similar haplotypes and morphotypes. In this case, zooplankton displays little or no geographic structure (Xiang et al., 2010). Long-distance dispersal and colonization of distant habitats, even transcontinental, have been recorded in several zooplankton taxa (De Gelas & De Meester, 2005; Mills et al., 2007; Campillo et al., 2011). Secondly, a rapid adaptation of resident populations to local environmental conditions after a historical colonization event reduces effective gene flow among populations, and cause high levels of genetic differentiation and population subdivision (De Meester, Gómez, Okamura, & Schwenk, 2002).

High dispersal and colonization capacity are evident in B. quadridentatus populations investigated in this study. This is demonstrated by the presence of highly divergent haplotypes and morphotypes within La Colorada, La Zarcita, Vernet and Charco 1 locations (Fig. 1, Fig. 3), which contain specimens belonging to different lineage separated by > 18 mutation steps and that are also morphologically differentiated. This pattern can be explained as a secondary contact between distinct lineages. Other examples that may have stemmed from high dispersal capacity in the group are the presence of specimens showing little genetic and morphological differences within the lineages B, D and G (see Fig. 3, Table S1), but separated by large distances and biogeographic barriers. For example, the most widespread lineage was G, encompassing north, central and southeast Mexico and unexpectedly also found in the Iberian Peninsula (see Fig. S5). Closely genetic and morphologically related haplotypes in lineage G were collected in El Rosario, Baja California and Charco 1, Tabasco, locations that are separated by distances > 1900 km. This suggests a high colonization and adaptation ability of the lineages. A similar scenario was also found in the B. calyciflorus species complex (Xiang et al., 2010), in which wide geographic distribution from one of the cryptic species was attributed to its high capacity for adaptation and colonization of several habitats in East China (Xiang et al., 2010).

In the present study, mitochondrial analysis recovered seven (A-G) deeply divergent lineages, which also exhibited exclusive haplotypes and morphotypes. This may suggests a low level of gene flow among the seven lineages, even between those found in the same water body. Therefore, the high level of genetic structure and the high number of private haplotypes suggests a history of geographic isolation between these lineages. In addition, the high level of genetic differentiation among lineages was not related to the geographic distance with some highly divergent lineages found in geographic proximity. In this case other factors may influence the genetic differentiation among lineages of B. quadridentatus. Similarly, several molecular surveys conducted on rotifers, including Brachionus, Synchaeta and Euchlanis have found high levels of genetic differentiation despite high dispersal ability of the studied species (Gómez et al., 2000; 2002; Xiang et al., 2010; Obertegger, Fontaneto, & Flaim, 2012; Kimpel et al., 2015; Kordbacheh et al., 2017).

For instance, the lineages E and F, which contain specimens from La Joya and Cuitzeo Lakes, respectively, have a high genetic and morphological differentiation with respect to other lineages. These lineages are separated by < 50 km from other highly divergent lineages, thus the observed differentiation is unlikely to be due to dispersal limitations. La Joya and Cuitzeo Lakes display salt levels up to 5 g/L (Sarma, Nandini, Morales-Ventura, Delgado-Martínez, & González-Valverde, 2006), therefore having very distinct environmental conditions from other studied water bodies. Cuitzeo is the second largest lake in Mexico (± 400 km2), and despite the high number of samples collected all around the lake, it is surprising that only one genetic and morphotype group with low genetic diversity was found. This could be explained by founder effects and local adaptation to the environmental conditions of the lakes as has been shown in other studies (e. g. Ishida & Taylor, 2007; Xiang et al., 2010). In addition, the haplotype and nucleotide diversity (Table 1), as well as the star like shape of the haplotype network of the Cuitzeo Lake (Fig. S6), suggest that population experienced a bottleneck with a subsequent expansion. Past bottleneck events and subsequent range expansions have also been documented for other zooplankton taxa (Ishida & Taylor, 2007; Crease, Omilian, Costanzo, & Taylor, 2012).

In summary, this is the first study conducted on the littoral rotifer B. quadridentatus integrating molecular data from COI and 18S genes to examine cryptic diversity in this taxon. The diversity delimited in our study was higher in the mitochondrial COI gene than the nuclear 18S gene, thus providing evidence of mitonuclear discordance. This incongruence might be due to differences in mutation rate between markers, resulting in incomplete lineage sorting. Results from analysis of the COI gene along with the morphotypes observed here support the existence of seven putative species. Moreover, five additional COI lineages were found from sequences obtained from Europe and China (Table S7). Therefore the results presented herein support B. quadridentatus as a species complex with a high diversity comparable to that of the B. plicatilis and B. calyciflorus species complexes. Integrative taxonomic frameworks, including additional nuclear markers, as well as morphological, reproductive, and ecological analyses including samples from all over the world are needed in order to assess the species richnes, species boundaries and mechanisms driving the speciation in B. quadridentatus.

Ethical statement: authors declare that they all agree with this publication and made significant contributions; that there is no conflict of interest of any kind; and that we followed all pertinent ethical and legal procedures and requirements. All financial sources are fully and clearly stated in the acknowledgements section. A signed document has been filed in the journal archives.

ACKNOWLEDGMENTS

This study is part of the postdoctoral tenure of Alma E. García-Morales at the Universidad Michoacana de San Nicolás de Hidalgo and was supported by the Consejo Nacional de Ciencia y Tecnología (CONACYT) grant. Financial support was provided by grant numbers: CIC-UMSNH-2014 and CIC-UMSNH-2015. We grateful Nadia Jogee for the helpful comments to a revision of the English version. We thank Tom Jameson who provided constructive criticism and a revision of the English of this manuscript. GenBank numbers: COI: MN627001-MN627154; 18S: MN636060-MN636184.

RESUMEN

Diversidad molecular críptica en el rotífero morfológicamente variable Brachionus quadridentatus (Rotifera: Monogononta). Brachionus quadridentatus es una especie morfológicamente variable, distribuida por todo el mundo. Su taxonomía es confusa debido a las numerosas variantes infrasubespecíficas descritas en este taxón. Con la taxonomía basada en la morfología, B. quadridentatus tiene tres variantes reconocidas: B. quadridentatus quadridentatus, B. quadridentatus f. brevispinus and B. quadridentatus f. cluniorbicularis. En este estudio, exploramos la diversidad genética entre algunas poblaciones de B. quadridentatus, usando secuencias de los genes COI ADNmt y 18S ARNr. El análisis de delimitación de especies coalescente usando el gen 18S apoya la presencia de al menos tres especies putativas dentro del complejo B. quadridentatus. Estos resultados estuvieron en concordancia con los análisis filogenético y GMYC usando el gen 18S. Sin embargo, se encontró variación en morfología y secuencias del gen COI dentro de cada una de las tres especies putativas. Se encontraron siete linajes delimitados por las secuencias del gen COI usando el método de delimitación ABGD, que además están morfologicamente diferenciadas. Se encontró discordancia mitonuclear entre la filogenia del gen COI y la del gen 18S. La incongruencia entre el marcador mitocondrial y el nuclear puede ser explicada por sorteo incompleto de linaje.

Palabras clave: biodiversidad; citocromo c oxidasa 1; 18S; monogononto; árbol de especies.

REFERENCES

Bouckaert, R., Heled, J., Kuehnert, D., Vaughan, T., Wu, C. H., Xie D., … Drummond, A. J. (2014). BEAST 2: A software platform for Bayesian Evolutionary Analysis. PLOS Computational Biology, 10: e1003537.

Campillo, S., García-Roger, E. M., Martínez-Torres, D., & Serra M. (2005). Morphological stasis of two species belonging to the L-morphotype in the Brachionus plicatilis species complex. Hydrobiologia, 546, 181-187.

Campillo, S., Serra, M., Carmona, M. J., & Gómez, A. (2011). Widespread secondary contact and new glacial refugia in the halophilic rotifer Brachionus plicatilis in the Iberian Peninsula. PLoS ONE, 6, e20986.

Cieplinski, A., Weisse, T., & Obertegger, U. (2017). High diversity in Keratella cochlearis (Rotifera, Monogononta): morphological and genetic evidence. Hydrobiologia, 796, 145-159.

Clement, M., Posada, D., & Crandall, K. A. (2000). TCS: a computer program to estimate gene genealogies. Molecular Ecology, 9, 1657-1659.

Crease, T. J., Omilian, A. R., Costanzo, K. S., & Taylor, D. J. (2012). Transcontinental phylogeography of the Daphnia pulex species complex. PLoS ONE, 7, e46620.

Darriba, D., Taboada, G. L., Doallo, R., & Posada, D. (2012). jModelTest 2: more models, new heuristics and parallel computing. Nature Methods, 9, 772.

De Gelas, K., & De Meester, L. (2005). Phylogeography of Daphnia magna in Europe. Molecular Ecology, 14, 753-764.

De Meester, L., Gómez, A., Okamura, B., & Schwenk, K. (2002). The monopolization hypothesis and the dispersal-gene flow paradox in aquatic organisms. Acta Oecologica, 23, 121-135.

Dumont, H. (1983). Biogeography of rotifers. Hydrobiologia, 104, 19-30.

Excoffier, L., Smouse, P., & Quattro, J. (1992). Analysis of molecular variance inferred from metric distances among DNA haplotypes: applications to human mitochondrial DNA restriction data. Genetics, 131, 479-491.

Folmer, O., Black, M., Hoeh, W., Lutz, R., & Vrijenhoek, R. (1994). DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Molecular and Marine Biology and Biotechnology, 3, 294-297.

Fontaneto, D., Kaya, M., Herniou, E. A., & Barraclough, T. G. (2009). Extreme levels of hidden diversity in microscopic animals (Rotifera) revealed by DNA taxonomy. Molecular Phylogenetics and Evolution, 53, 182-189.

Fontaneto, D., Flot, J. F., & Tang, C. Q. (2015). Guidelines for DNA taxonomy with focus on the meiofauna. Marine Biodiversity, 45, 433-451.

García-Morales, A. E., & Elías-Gutiérrez, M. (2013). DNA barcoding of freshwater Rotifera in Mexico: evidence of cryptic speciation in common rotifers. Molecular Ecology Resources, 13, 1097-1107.

Gilbert, J. J. (2001). Spine development in Brachionus quadridentatus from an Australian billabong: genetic variation and induction by Asplanchna. Hydrobiologia, 446/447, 19-28.

Gilbert, J. J., & Walsh, E. J. (2005). Brachionus calyciflorus is a species complex: mating behavior and genetic differentiation among four geographically isolated strains. Hydrobiologia, 546, 257-265.

Gómez, A., Carvalho, G. R., & Lunt, D. H. (2000). Phylogeography and regional endemism of a passively dispersing zooplankter: mitochondrial DNA variation in rotifer resting egg banks. Proceedings of the Royal Society of London, Series B: Biological Sciences, 267, 2189-2197.

Gómez, A., Serra, M., Carvalho, G. R., & Lunt, D. H. (2002). Speciation in ancient cryptic species complex: evidence from the molecular phylogeny of Brachionus plicatilis (Rotifera). Evolution, 56, 1431-1444.

Hebert, P. D. N., Cywinska, A., Ball, S. L., & De Waard, J. R. (2003). Biological identifications through DNA barcodes. Proceedings of the Royal Society of London, Series B: Biological Sciences, 270, 313-321.

Heled, J., & Drummond, A. J. (2010). Bayesian inference of species trees from multilocus data. Molecular Biology and Evolution, 27, 570-580.

Ishida, S., & Taylor, D. J. (2007). Mature habitats associated with genetic divergence despite strong dispersal ability in an arthropod. BMC Evolutionary Biology, 7, 52.

Jörger, K. M., Norenburg, J. L., Wilson, N. G., & Schrödl, M. (2012). Barcoding against a paradox? Combined molecular species delineations reveal multiple cryptic lineages in elusive meiofaunal sea slugs. BMC Evolutionary Biology, 12, 245.

Katoh, K., Asimenos, G., & Toh, H. (2009). Multiple alignment of DNA sequences with MAFFT. Methods in Molecular Biology, 537, 39-64.

Kimpel, D., Gockel, J., Gerlach, G., & Bininda-Emonds, O. R. P. (2015). Population structuring in the monogonont rotifer Synchaeta pectinata: high genetic divergence on a small geographical scale. Freshwater Biology, 60, 1364-1378.

Kordbacheh, A., Garbalena, G., & Walsh, E. J. (2017). Population structure and cryptic species in the cosmopolitan rotifer Euchlanis dilatata. Zoological Journal of the Linnean Society, 181, 757-777.

Kordbacheh, A., Wallace, R. L., & Walsh, E. J. (2018). Evidence supporting cryptic species within two sessile microinvertebrates, Limnias melicerta and L. ceratophylli (Rotifera, Gnesiotrocha). PLos ONE, 13(10): e0205203. Retrieved from https://doi.org/10.1371/journal.pone.0205203

Koste, W. (1978). Rotatoria. Die Rädertiere Mitteleuropas. Ein Bestimmungswerk, begründet von Max. Vol 2 Tafelband. Germany: Voigt Überordnung Monogononta.

Koste, W., & Shiel, R. J. (1987). Rotifera from Australian inland waters II Ephiphanidae and Brachionidae (Rotifera: Monogononta). Invertebrate Taxonomy, 1, 949-1021.

Koste, W., & Shiel, R. J. (1989). Rotifera from Australian inland waters III. Euchlanidae, Mytilinidae and Trichotriidae (Rotifera: Monogononta). Transactions of the Royal Society of South Australia, 113, 85-114.

Leasi, F., Tang, C. Q., De Smet, W., & Fontaneto, D. (2013). Cryptic diversity with wide salinity tolerance in the putative euryhaline Testudinella clypeata (Rotifera, Monogononta). Zoological Journal of the Linnean Society, 168, 17-28.

Librado, P., & Rozas, J. (2009). DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics, 25, 1451-1452.

Mantel, N. (1967). The detection of disease clustering and a generalized regression approach. Cancer Research, 27, 209-220.

Mark Welch, D. B., & Mark Welch, J. L. (2005). The potential of genomic approaches to rotifer ecology. Hydrobiologia, 546, 101-108.

Mills, S., Lunt, D. H., & Gómez, A. (2007). Global isolation by distance despite strong regional phylogeography in a small metazoan. BMC Evolutionary Biology, 7, 225.

Mills, S., Alcántara-Rodríguez, J. A., Ciros-Pérez, J., Gómez, A., Hagiwara, A., Hinson Galindo, K., … Walsh, E. J. (2017). Fifteen species in one: deciphering the Brachionus plicatilis species complex (Rotifera, Monogononta) through DNA taxonomy. Hydrobiologia, 796, 39-58.

Moreno, E., Conde-Porcuna, J. M., & Gómez, A. (2017). Barcoding rotifer biodiversity in Mediterranean ponds using diapausing egg banks. Ecology and Evolution, 7, 4855-4867.

Obertegger, U., Fontaneto, D., & Flaim, G. (2012). Using DNA taxonomy to investigate the ecological determinants of plankton diversity: explaining the occurrence of Synchaeta spp. (Rotifera, Monogononta) in mountain lakes. Freshwater Biology, 57, 1545-1553.

Obertegger, U., Flaim, G., & Fontaneto, D. (2014). Cryptic diversity within the rotifer Polyarthra dolichoptera along an altitudinal gradient. Freshwater Biology, 59, 2413-2427.

Obertegger, U., Cieplinski, A., Fontaneto, D., & Papakostas, S. (2018). Mitonuclear discordance as a confounding factor in the DNA taxonomy of monogonont rotifers. Zoologica Scripta, 47, 122-132.

Ow, Y. X., & Todd, P. A. (2010). Light-induced morphological plasticity in the scleractinian coral Goniastrea pectinata and its functional significance. Coral Reefs, 29, 797-808.

Papakostas, S., Michaloudi, E., Proios, K., Brehm, M., Verhage, L., Rota, J., … Declerck, S. A. J. (2016). Integrative taxonomy recognizes evolutionary units despite widespread mitonuclear discordance: Evidence from a rotifer cryptic species complex. Systematic Biology, 65, 508-524.

Pons, J., Barraclough, T. G., Gómez-Zurita, J., Cardoso, A., Duran, D. P., Hazell, S., Kamoun, S., Sumlin, W. D., & Vogler, A. (2006). Sequence based species delimitation for the DNA taxonomy of undescribed insects. Systematic Biology, 55, 595-609.

Puillandre, N., Lambert, A., Brouillet, S., & Achaz, G. (2012). ABGD, Automatic Barcode Gap Discovery for primary species delimitation. Molecular Ecology, 21, 1864-1877.

Rambaut, A., Suchard, M. A., & Drummond, A. J. (2013). Tracer v1.6. Retrieved from http://beast.bio.ed.ac.uk/Tracer

Ronquist, F., & Huelsenbeck, J. P. (2003). Mrbayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics, 19, 1572-1574.

Sarma, S. S. S., Nandini, S., Morales-Ventura, J., Delgado-Martínez, I., & González-Valverde, L. (2006). Effects of NaCl salinity on the population dynamics of freshwater zooplankton (rotifers and cladocerans). Aquatic Ecology, 40, 349-360.

Schröder, T., & Walsh, E. J. (2007). Cryptic speciation in the cosmopolitan Epiphanes senta complex (Monogononta, Rotifera) with the description of new species. Hydrobiologia, 593, 129-140.

Segers, H. (2007). Annotated checklist of the rotifers (Phylum Rotifera), with notes on nomenclature, taxonomy and distribution. Zootaxa, 1564, 1-104.

Segers, H. (2008). Global diversity of rotifers (Rotifera) in freshwater. Hydrobiologia, 595, 49-59.

Segers, H., & De Smet, W. H. (2008). Diversity and endemism in Rotifera: a review, and Keratella Bory de St Vincent. Biodiversity and Conservation, 17, 303-316.

Serra, M., Galiana, A., & Gómez, A. (1997). Speciation in monogonont rotifers. Hydrobiologia, 358, 63-70.

Stamatakis, A., Hoover, P., & Rougemont, J. (2008). A rapid bootstrap algorithm for the RAxML web servers. Systematic Biology, 57, 758-771.

Suatoni, E., Vicario, S., Rice, S., Snell, T., & Caccone, A. (2006). An analysis of species boundaries and biogeographic patterns in a cryptic species complex: the rotifer Brachionus plicatilis. Molecular Phylogenetics and Evolution, 41, 86-98.

Tamura, K., Stecher, G., Peterson, D., Filipski, A., & Kumar, S. (2013). MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Molecular Biology and Evolution, 30, 2725-2729.

Tang, C. Q., Leasi, F., Obertegger, U., Kieneke, A., Barraclough, T. G., & Fontaneto, D. (2012). The widely used small subunit 18S rDNA molecule greatly underestimates true diversity in biodiversity surveys of the meiofauna. PNAS, 109, 16208-16212.

Tang, C. Q., Obertegger, U., Fontaneto, D., & Barraclough, T. G. (2014). Sexual species are separated by larger genetic gaps than asexual species in rotifers. Evolution, 68, 2901-2916.

Wallace, R. L. (2002). Rotifers: exquisite metazoans. Integrative and Comparative Biology, 42, 660-667.

Wallace, R. L., Snell, T. W., Ricci, C., & Nogrady, T. (2006). Rotifera volume 1: Biology, ecology and systematics. In H. J. Dumont, R. L. Wallace, T. W. Snell, C. Ricci, & T. Nogrady (Eds), Guides to the identification of the microinvertebrates of the continental waters of the world (2nd Edition, pp. 1-20). The Netherlands: Kenobi Productions and Backhuys Publishers.

Walsh, E. J., Schröder, T., Wallace, R. L., & Rico-Martinez, R. (2009). Cryptic speciation in Lecane bulla (Monogononta: Rotifera) in Chihuahuan Desert waters. Verhandlungen des Internationalen Verein Limnologie, 30, 1046-1050.

Wilke, T., Schultheiß, R., & Albrecht, C. (2009). As time goes by: a simple fool´s guide to molecular clock approaches in invertebrates. American Malacological Bulletin, 27, 25-45.

Xiang, X. L., Xi, Y. L., Wen, X. L., Zhang, J. Y., & Ma, Q. (2010). Spatial patterns of genetic differentiation in Brachionus calyciflorus species complex collected from East China in summer. Hydrobiologia, 638, 67-83.

Xiang, X. L., Xi, Y. L., Wen, X. L., Zhang, G., Wang, J. X., & Hu, K. (2011). Genetic differentiation and phylogeographical structure of the Brachionus calyciflorus complex in eastern China. Molecular Ecology, 20, 3027-3044.

Xiang, X. L., Xi, Y. L., Zhu, L. Y., & Xu, Q. L. (2017). Comparative studies of the population genetic structure of the Brachionus calyciflorus species complex from four inland lakes in Wuhu, China. Biochemical Systematics and Ecology, 71, 69-77.

Zhang, J., Kapli, P., Pavlidis, P., & Stamatakis, A. (2013). A general species delimitation method with applications to phylogenetic placements. Bioinformatics, 29, 2869-2876.

TABLE 2

Results of the different species delimitation methods conducted on Brachionus quadridentatus

Method	COI	18S	Concatenated
SP	13	1	8
PTP	13	4	10
GMYC	11 (5-13)	3 (1-5)	9 (7-13)
ABGD	7	NA	NA

NA means that the test cannot be performed on the dataset.