Mitonuclear discordance across taxa is increasingly recognized as posing a major challenge to species delimitation based on DNA sequence data. Integrative taxonomy has been proposed as a promising framework to help address this problem. However, we still lack compelling empirical evidence scrutinizing the efficacy of integrative taxonomy in relation to, for instance, complex introgression scenarios involving many species. Here, we report remarkably widespread mitonuclear discordance between about 15 mitochondrial and 4 nuclear Brachionus calyciflorus groups identified using different species delimitation approaches. Using coalescent-, Bayesian admixture-, and allele sharing-based methods with DNA sequence or microsatellite data, we provide strong evidence in support of hybridization as a driver of the observed discordance. We then describe our combined molecular, morphological, and ecological approaches to resolving phylogenetic conflict and inferring species boundaries. Species delimitations based on the ITS1 and 28S nuclear DNA markers proved a more reliable predictor of morphological variation than delimitations using the mitochondrial COI gene. A short-term competition experiment further revealed systematic differences in the competitive ability between two of the nuclear-delimited species under six different growth conditions, independent of COI delimitations; hybrids were also observed. In light of these findings, we discuss the failure of the COI marker to estimate morphological stasis and morphological plasticity in the B. calyciflorus complex. By using B. calyciflorus as a representative case, we demonstrate the potential of integrative taxonomy to guide species delimitation in the presence of mitonuclear phylogenetic conflicts.
Integrative taxonomy aims to delimit species using knowledge acquired from multiple complementary perspectives such as morphology, patterns of mitochondrial and nuclear DNA diversity, and ecology ( Dayrat 2005 ; Padial et al. 2010 ). Resolving conflicts between perspectives in integrative taxonomy allows for the prioritization of criteria to obtain rigorous species-level taxonomies ( Schlick-Steiner et al. 2009 ; Andujar et al. 2014 ). By doing so, we can reciprocally inform the description of the morphology or the ecology of the species and achieve greater accuracy or avoid biases ( Wielstra and Arntzen 2014 ). Furthermore, the opportunity then arises to study evolutionary phenomena associated with these conflicts that would otherwise have remained cryptic ( Schlick-Steiner et al. 2014 ).
Introgressive hybridization (in short, introgression), often revealed by discordant patterns between mitochondrial and nuclear phylogenies (mitonuclear discordance), is now considered more prevalent than was previously thought ( Petit and Excoffier 2009 ; Toews and Brelsford 2012 ). Introgression can affect species integrity, obscure species characters, lead to phylogenetic conflict, and thus mislead species identifications ( Petit and Excoffier 2009 ). In the light of frequent introgression across disparate taxonomic groups, a pressing question emerges: how can we critically recognize evolutionary significant units of diversity?
Zooplankton organisms may offer suitable models to assess the problem of mitonuclear discordance in species delimitation. The frequently complex life cycles, high dispersal capacities, and rapid local adaptations of planktonic species can, for instance, facilitate interspecific gene flow ( Cristescu et al. 2012 ). Cyclical parthenogens such as freshwater monogonont rotifers and cladocerans have long been known to be prone to hybridization ( Hebert 1985 ). Most of the work in this direction has been informed by genetic data ( Taylor and Hebert 1992 ; Xu et al. 2013 ). Molecular phylogenetics has been a valuable tool in understanding cryptic diversity in these microscopic organisms ( Fontaneto 2014 ), the cryptic species complex of the marine rotifer Brachionus plicatilis being a well-known case ( Gómez et al. 2002 ). Remarkably, very little is known about the degree of hybridization between cryptic species of Brachionus rotifers. Experimental crosses between B. plicatilis cryptic species have resulted in F1 hybrids in some cases ( Suatoni et al. 2006 ), although no evidence for the occurrence of such hybrids in nature has been found ( Mills et al. 2016 ). Cases of mitonuclear discordance have been observed between cryptic species of the complex of the freshwater rotifer Brachionus calyciflorus ( Xiang et al. 2011 ). However, whether these cases are the result of hybridization or incomplete lineage sorting (ILS) remains unknown.
In this study, we set out to explore the degree of mitonuclear discordance in the B. calyciflorus complex ( Gilbert and Walsh 2005 ). We based our analyses on rotifer individuals sequenced for the mitochondrial cytochrome oxidase subunit I gene (mtCOI) and the nuclear internal transcribed spacer 1 locus (nuITS1)—currently the two most widely used genetic markers in Brachionus rotifers. We also sequenced parts of the nuclear 18S (nu18S) and 28S (nu28S) ribosomal RNA genes as a way to validate our findings using the nuITS1 marker. In the absence of a single best approach to species delimitation ( Carstens et al. 2013 ; Flot et al. 2015 ; Fontaneto et al. 2015 ), units of diversity of putative species status were delineated using different criteria. To investigate the mechanism driving the observed discordance, we mainly employed methods based on coalescent simulations and Bayesian modeling of genetic admixture. For the latter, we genotyped our rotifers using a set of 12 recently developed microsatellite markers for B. calyciflorus ( Declerck et al. 2015 ). Furthermore, we aimed to address phylogenetic incongruence from both a morphological and an ecological perspective. We measured standard morphological characters in rotifer clonal lines and examined whether mitochondrial or nuclear delimitations are better predictors of the observed morphological variation. We carried out a competition experiment under a variety of growth conditions to test whether putative cryptic species show systematic differences in competitive strength. By so doing, this study assesses the utility of integrative taxonomy in overcoming mitonuclear discordance and guiding species delimitations.
M aterials and M ethods
Resting Egg Collection and Establishment of Clonal Lines
Rotifer sediment samples were collected from 22 sites in the Netherlands (Supplementary Appendix 1: Table S1, available on Dryad at http://dx.doi.org/10.5061/dryad.8rc4r ). Brachionus sp. resting eggs were separated from the sediment using a sugar flotation technique ( Gómez and Carvalho 2000 ) and hatched under light in Petri dishes using double-distilled H 2 O. Upon hatching, females of Brachionus calyciflorus Pallas, 1766 were identified and isolated under a stereoscope then used to set up clonal lines in the laboratory (two replicate clonal cultures per clonal line). Clonal cultures were maintained throughout the study by transferring about half of the culture (∼20 ml) to a clean tube with 20 ml of fresh culture medium each week. We used Chlamydomonas reinhardtii from nutrient-sufficient phytoplankton chemostats as food source, as described in Declerck et al. (2015) .
DNA Extraction, PCR Amplification, DNA Sequencing, and Genotyping
DNA was extracted from single rotifers using the HotSHOT method ( Montero-Pau et al. 2008 ). PCR amplification of a part of the mtCOI region (amplicon size: 642 bp) was performed using a set of specific primers for B. calyciflorus (LCOmodBc: 5′-GTCAACAAATCATAAAGATATTGGAACTC-3′, HCOmodBc: 5′- GGGTGACCAAAAAATCAAAATAARTGTT-3′). These primers were based on the LCO1490 and HCO2198 universal primers ( Folmer et al. 1994 ), redesigned following the same principles as described in Vasileiadou et al. (2009) . The complete nuITS1 region (between 296 bp and 313 bp in length) was amplified using the primers III: 5′-CACACCGCCCGTCGCTACTACCGATTG-3′, and VIII: 5′-GTGCGTTCGAAGTGTCGATGATCAA-3′ ( Palumbi 2006 ). Parts of the nu18S and nu28S ribosomal RNA genes (amplicon sizes: 588 bp and 824 bp, respectively) were amplified using the primers 1F: 5′-TACCTGGTTGATCCTGCCAGTAG-3′, and 4R: 5′-GAATTACCGCGGCTGCTGG-3′ ( Giribet et al. 1996 ) for nu18S, and the primers D1F: 5′-GGGACTACCCCCTGAATTTAAGCAT-3′, and Rd4b: 5′-CCTTGGTCCGTGTTTCAAGAC-3′ ( Crandall et al. 2000 ; Park and O'Foighil 2000 ); for nu28S. PCR conditions are provided in Supplementary Appendix 2: Section A, available on Dryad. Macrogen Europe (Amsterdam, The Netherlands) performed Sanger sequencing in both directions. Genotyping of 12 microsatellite loci developed for B. calyciflorus ( Declerck et al. 2015 ) was conducted with an ABI Prism 3130 DNA Analyzer (Applied Biosystems, CA) and the GeneMapper v.4.0 software (Applied Biosystems, CA). Microsatellite primer sequences and multiplex PCR conditions are also provided in Supplementary Appendix 2: Section A (see also Supplementary Appendix 2: Table S1), available on Dryad.
Alignment and Phylogenetic Inferences
Alignment of the mtCOI dataset was performed using the MUSCLE algorithm ( Edgar 2004 ) as implemented in the MEGA v.6 software ( Tamura et al. 2013 ). The mtCOI dataset consisted of 203 newly sequenced samples and 685 sequences downloaded from GenBank (Supplementary Appendix 1: Tables S2 and S3, available on Dryad). Coamplified mitochondrial pseudogenes located in the nucleus (numts) may be a confounding factor to species delimitations for mitochondrial datasets ( Song et al. 2008 ). Absence of numts in the mtCOI dataset was verified by inspecting the alignment of the translated sequences and conducting blastp searches against the UniProtKB database ( www.uniprot.org ). This step confirmed that the studied sequences were free of frameshift mutations and stop codons, which may be diagnostic of numts ( Bensasson et al. 2001 ).
Alignment of the nuITS1 dataset was performed using the mlocarna function of the LocARNA v.1.8.7 tool with default settings. LocARNA aligns noncoding RNAs by considering both sequence and secondary structure similarities ( Will et al. 2012 ). As the secondary structure of internal transcribed spacer regions has a role in the maturation of ribosomal RNA genes, it has been suggested that accounting for structural information when aligning ITS1 may improve the phylogenetic utility of this marker (e.g., Gottschling et al. 2001 ; Goertzen et al. 2003 ). The nuITS1 dataset consisted of 176 newly sequenced samples and 485 sequences downloaded from GenBank (Supplementary Appendix 1: Tables S2 and S3, available on Dryad). The nuITS1 alignment also included 4 bp and 82 bp from the neighboring 18S and 5.8S ribosomal RNA gene regions, respectively.
Alignments of the nu18S and nu28S datasets were performed using the MUSCLE algorithm. Secondary structure may also be important for phylogenetic analysis using ribosomal RNA genes ( Dixon and Hillis 1993 ), but the low levels of polymorphism found in the sequenced parts of these genes precluded the need for a structure-based alignment procedure. Because these datasets were used to validate the findings based on the nuITS1 marker, sequencing was performed on selected clonal lines that represented distinct units of diversity according to nuITS1. While the nu18S dataset consisted of 25 sequences due to the lack of polymorphism, we sequenced 93 samples in the case of the nu28S gene to account for potential intra-group variation.
To account for heterozygous individuals in the case of the nuclear markers, the forward and reverse chromatograms were also aligned using Sequencher 4 (Gene Codes) then inspected for double peaks indicative of heterozygosity ( Flot et al. 2006 ). The markers were phased using the approach described in Fontaneto et al. (2015) : this was trivial when chromatograms had a single double peak, but in the case of length-variant heterozygotes, which may contain many double peaks ( Flot et al. 2006 ), co-occurring sequences were separated using the program CHAMPURU ( Flot 2007 , available online at: http://seqphase.mpg.de/champuru/ —last accessed January 17, 2016). Because no phase information was available for the sequences downloaded from GenBank, this analysis was applied only to the samples sequenced in this study. All chromatograms are available in the Supplementary Material, available on Dryad.
The absence of substitution saturation in each dataset was assessed using the test of Xia et al. (2003) as implemented in the program DAMBE v.5 ( Xia 2013 ). Upon alignment, identical sequences were collapsed to single sequences. These unique sequence types, hereafter treated as haplotypes, were found for each marker using the program DNAsp v.5 ( Librado and Rozas 2009 ) and designated with a number prefixed with “Hap_” for those sequences also present in the dataset downloaded from GenBank or with “Hap_O” for those sequences found only in the newly sequenced samples. Haplotype alignments are found in the Supplementary Material, available on Dryad.
Bayesian and maximum likelihood phylogenetic inferences were performed using the program MrBayes v.3.2.5 ( Ronquist et al. 2012 ) and a rapid-bootstrap version of the RAxML v.7.7.1 algorithm for web servers ( Stamatakis et al. 2008 ), respectively. In the case of Bayesian analyses, two independent runs were carried out for 20 million generations, with one cold and three heated chains, and with a tree being sampled every 2000 generations. To avoid the problem of long-tree solutions ( Marshall 2010 ; Brown et al. 2010 ), the prior for branch length was set to brlensp = unconstrained:Exp(100). Input files for MrBayes with all the parameters can be found in the Supplementary Material, available on Dryad. Convergence was assessed by (i) examination of the potential scale reduction factor (PSRF, ; maximum 1.012), (ii) examination of the average standard deviation of split frequencies ( 0.014), as well as (iii) using Tracer v.1.6 ( Rambaut et al. 2014 ) by requiring an effective sample size (ESS) values above 200 for all parameters. Trees were summarized using the sumtrees command in DendroPy v.3.12.0 ( Sukumaran and Holder 2010 ) after discarding the first 20% of the trees as burn-in. The RAxML analysis was run with 100 bootstraps. A B. plicatilis sequence was used as an outgroup for the mtCOI phylogeny (GenBank accession: AY785179; Suatoni et al. 2006 ). For nuITS1, fitting an outgroup to the B. calyciflorus alignment [e.g., using the add and keeplength functions of the MAFFT v.7.266 program ( Katoh and Standley 2013 )] proved a difficult task. Due to this reason, and because this study focuses on phylogenetic relationships within the B. calyciflorus complex, we omitted the use of an outgroup in the nuITS1 dataset. In this case, for visualization purposes, in trees used for ultrametric tree conversions, and in the case of the Poisson tree process (PTP) method (see below), we applied the method of midpoint rooting using the program FigTree v.1.4.2 (available from: http://tree.bio.ed.ac.uk/software/figtree/ —last accessed January 17, 2016). In the case of the nu28S dataset, the lack of corresponding 28S sequence for B. plicatilis led us to use a sequence from Plationus patulus (GenBank accession: FR729700; Stelzer et al. 2011 ) as an appropriate outgroup ( Reyna-Fabian et al. 2010 ). The absence of polymorphism in the sequenced part of the 18S gene did not allow for any phylogenetic analysis. The best-fitting models of nucleotide substitution were determined following the Bayesian Information Criterion (BIC; Schwarz 1978 ) using jModelTest 2 ( Darriba et al. 2012 ). Overall, 88 models and 11 substitution schemes were tested, and the likelihoods of the models were calculated using maximum-likelihood topologies resulting from heuristic searches using the subtree pruning and regrafting algorithm as implemented in PhyML v.3.0 ( Guindon et al. 2010 ).
Species delimitation was based on three main types of approaches ( Flot et al. 2015 ). First, we used two different tree-based coalescent methods, the generalized mixed Yule-coalescent model (GMYC; Pons et al. 2006 ), and the PTP method ( Zhang et al. 2013 ). Second, we performed automatic barcode gap discovery (ABGD), a distance-based method ( Puillandre et al. 2012 ). Third, we investigated species boundaries using a haploweb ( Flot et al. 2010 ), which is an allele sharing-based approach. GMYC uses ultrametric trees (i.e., trees whose branch lengths are proportional to time) to calculate the most likely threshold between interspecific (modeled as a Yule process) and population-level branching rates (modeled as a coalescent process), thereby delineating evolutionary significant units akin to species ( Fontaneto et al. 2015 ). PTP does not require ultrametric trees and distinguishes between population- and species-level processes by assuming that intra- and interspecific substitutions follow two distinct Poisson processes ( Fontaneto et al. 2015 ). ABGD uses pairwise genetic distances to determine the gap between intra- and interspecific divergence and delimit primary species hypotheses ( Puillandre et al. 2012 ). Haplowebs delineate reproductively isolated gene pools using information from haplotypes found co-occurring in heterozygous individuals ( Flot et al. 2010 ). Unlike GMYC, PTP, and ABGD that can be performed on haploid as well as diploid markers, haplowebs require diploid nuclear markers in which heterozygous individuals are detected as having double peaks in the sequencing chromatograms ( Flot et al. 2006 ). In all cases, species delimitations were performed on the ingroup. Whenever applicable, to remove the outgroup from the trees we used the drop.tip function of the R package “ape” ( Paradis et al. 2004 ). In the case of tree-based methods, GMYC and PTP, this step ensured that no processes of diversification between higher taxa were involved in the investigated tree topology.
GMYC was employed on ultrametric trees calculated with three different methods for each marker. First, summarized chronograms were generated using the program BEAST v.1.8.2 ( Drummond et al. 2012 ), which is currently the best recognized practice for GMYC ( Tang et al. 2014 ). Second, and in order to control for potential methodological biases, we also applied GMYC on conversions of the summarized Bayesian trees (obtained from MrBayes) to ultrametric trees using the penalized likelihood criterion as implemented in the program r8S v.1.7 ( Sanderson 2003 ). Third, we applied GMYC on conversions of the summarized Bayesian trees to ultrametric trees using the program PATHd8 according to the mean-path length method ( Britton et al. 2007 ). BEAST was run for 60 million generations and a tree was sampled every 6000 generations with the substitution parameters suggested by the best fitting model for each marker under a lognormal relaxed (uncorrelated) clock [following Monaghan et al. 2009 and Wertheim et al. 2009 ], with a constant-size coalescent tree prior. BEAUti files with all the BEAST parameters for each marker can be found in Supplementary Material, available on Dryad. Convergence was assessed with Tracer v.1.6 by inspecting the ESS of all parameters. The summarized trees were calculated using the program TreeAnnotator v.1.8.2 (which is part of the BEAST package) by keeping the node heights of the highest log clade credibility identified after discarding the first 20% of the trees as burn-in. For ultrametric tree conversions with the programs r8s and PATHd8, the age of the most recent common ancestor at the root of the tree was arbitrarily set to 100, while polytomies were resolved randomly with zero-branch length dichotomies using the multi2di function of the R package “ape”. The smoothing parameter for the calculations using penalized likelihood in the r8s program was set to 1.50 for the mtCOI and nu28S markers and to −2.00 for nuITS1 (following crossvalidation of values ranging from −6.00 to 6.00 in increments of 0.50). All calculations with the GMYC model were performed using the single-threshold method, which in a recent simulation study exhibited a bias toward overlumping rather than oversplitting ( Dellicour and Flot 2015 ).
To account for uncertainty in tree space and in the parameters of the GMYC model, we also employed a Bayesian implementation of the model (bGMYC) using trees from the BEAST analysis. The bGMYC method allows analysis of multiple post-burn-in topologies and performs Markov Chain Monte Carlo (MCMC) simulations to examine the posterior distribution of the GMYC model ( Reid and Carstens 2012 ); in the simulation study of Dellicour and Flot (2015) , bGMYC was found to perform better than the single-locus GMYC model and to present a tendency toward oversplitting rather than overlumping. The bgmyc.multiphylo function of the R package “bGMYC” v.1.0.2 was run using 100 trees sampled from the BEAST run as follows: a tree was sampled every 480 thousand generations using LogCombiner v.1.8.2 (part of the BEAST package) and then the first 20% of the trees were discarded as burn-in. This step ensured a uniform sampling of the post-burn-in trees. bGMYC was then run for 110,000 iterations, with burn-in set to 10,000, and sampling every 100th step. The length of the burn-in period was decided after observing convergence in less than 1000 generations using the bgmyc.singlephy function and the highest clade credibility tree—identified with TreeAnnotator—of the BEAST run. To increase the accuracy of the method, we used estimates of the upper limit of number of species, t2 , based on the result of the other species delimitation methods. Convergence was assessed by inspecting the posterior probability of the simulations against the number of generations. In the end, the function bgmyc.point was used to assess conspecificity at the posterior probability level , while the function plot.bgmycprobmat was used to plot the matrix of probability of conspecifity onto the summarized Bayesian phylogeny.
Calculations with the PTP method were performed using the “PTP” Python package (v.2.2; date: 14 February 2014) using the default settings, notably requiring a -value of 0.0001 for the presence of distinct intra- and interspecific branch length classes to be considered. Because the PTP method does not require ultrametric trees, we applied the method directly on both the best scoring tree obtained from the RAxML analysis and the highest log clade credibility tree identified using TreeAnnotator from the MrBayes analysis. The same two trees were also employed for determining the barcode gap with the ABGD method. Pairwise genetic distances were calculated directly from the trees using the cophenetic function of the R package “ape,” which allowed calculations to be made accounting for the models used to build each of the trees.
Haploweb analyses were performed as described in Flot et al. (2010) . All identified haplotypes from the chromatograms were used to construct a median-joining haplotype network using the program Network v.4.613 (available online at http://www.fluxus-engineering.com/sharenet.htm —last accessed 17 January 2016). Haplotype frequencies were estimated at the level of clonal lines, and the pattern of co-occurrence of alleles in heterozygotes was used to determine fields for recombination of putative cryptic species status ( Doyle et al. 1995 ; Flot et al. 2010 ).
Lastly, we evaluated the degree of agreement between all different species delimitation methods and tried to come up with a consensus for species delimitations. To account for the tendency of some phylogenetic delimitation criteria to oversplit ( Fontaneto et al. 2015 ; but see Dellicour and Flot 2015 ), we also considered more conservative criteria for species number in the case of the tree-based delimitation methods. For each of the GMYC-based approaches (using BEAST, r8s, and PATHd8), we considered the delimitations at the lower limit of the 95% confidence interval. In bGMYC we calculated the delimitations at a lower conspecificity probability threshold ( ). Finally, for the PTP method we considered the lower number of species estimates obtained from the Bayesian and maximum-likelihood trees.
Testing for Hybridization
To assess hybridization as a potential mechanism responsible for the observed mitonuclear discordance, we employed two different approaches. First, we used the DNA sequence data and simulated gene trees under the multispecies coalescent model according to the method implemented in the program JML v.1.3.0 ( Joly 2012 ). With JML we tested whether the minimum interspecific genetic distance observed using mtCOI sequences was significantly smaller than that predicted from the posterior distribution of species trees calculated using either the nuITS1 or the nu28S marker. In this way, we tested the null hypothesis that ILS was sufficient to explain the observed discordance ( Joly et al. 2009 ). Second, we used the microsatellite genotypes with two Bayesian-based methods for the detection of genetic admixture in our samples as implemented in the programs STRUCTURE v.2.3.4 ( Pritchard et al. 2000 ) and NewHybrids v.1.1 ( Anderson and Thompson 2002 ).
Tests with JML were performed separately for two different marker combinations, mtCOI-nuITS1 or mtCOI-nu28S. To obtain accurate results we used only rotifer samples that were sequenced for both markers in each case, and considered every haplotype combination only once (Supplementary Appendix 1: Table S2 for mtCOI-nu28S and Appendix 1: Table S4 for mtCOI-nuITS1, available on Dryad). We acknowledge that nuITS1 and nu28S are physically connected and do not represent independent replicates of nuclear markers; yet, given their different evolutionary rates, our analyses implicitly tested for the potential effects of such differences by repeating the analyses for each of the two nuclear markers. Because ILS is more likely for nuclear than for mitochondrial markers (e.g., Zink and Barrowclough 2008 ), coalescent simulations were performed using each of the nuclear markers, while mtCOI sequences were used to estimate the observed minimum interspecific genetic distances. Species were delimited according to the nuITS1 marker as suggested by our integrative taxonomic approach (also supported by the nu28S gene, see below).
Coalescent simulations were obtained by running *BEAST ( Heled and Drummond 2010 ) for 150 million generations with a sample taken every 15,000 generations. Two independent runs were performed for each nuclear marker and the two runs were merged using the LogCombiner program of the BEAST package, after discarding 20% of each run as burn-in. In each case, we also considered either a Yule (pure birth) or a birth–death process as a prior for the species tree, and the results obtained with each process were compared using the program Tracer v.1.6 according to the harmonic mean of the combined likelihood trace of the two runs with 100 bootstrap replicates. The nucleotide substitution model was set as suggested by jModelTest2, with the clock model set to “lognormal relaxed clock,” and the population size model set to “piecewise constant” as suggested by the author of the JML program. BEAUti files with all the settings of the *BEAST runs are available in the Supplementary Material, available on Dryad.
For the JML tests, due to computer memory limitations, the post-burn-in merged simulations were thinned to one-third, that is, one sampled tree every 45,000 *BEAST generations, or a total of 5334 simulations after 20% burn-in. We used the best fitting model of nucleotide substitution in each case, and the relative locus mutation rate as was estimated from the mean locus rate of *BEAST runs that combined each of the nuclear loci with the mtCOI gene in the same run. Population sizes for the simulations were scaled using the appropriate relative heredity scalar of one-fourth, given that the effective population size of nuclear loci is typically four times that of mitochondrial genes ( Birky et al. 1989 ). We acknowledge this may be a debatable assumption in the case of the multi-copy nuclear markers (such as nuITS1 and nu28S) due to gene conversion. For this reason, we validated our use of JML with nuITS1 and mtCOI sequences from 14 species of the B. plicatilis complex for which no mitonuclear discordance has been observed ( Mills et al. 2016 ; Supplementary Appendix 1: Table S5, available on Dryad). We hypothesize that the JML test in this dataset, if not strongly influenced by the effect of gene conversion ( Hartfield et al. 2016 ), would not reject the null hypothesis of ILS for any of the pairwise species comparisons. Analyses for the B. plicatilis complex were performed as described for the B. calyciflorus dataset (detailed in Supplementary Appendix 2: Section B, available on Dryad). To further account for false positives, a potential confounding factor in JML tests ( Heled et al. 2013 ), we used the program QVALUE and calculated the false positive rate at the applied significance threshold ( Storey and Tibshirani 2003 ). -values provide an extension of the false discovery rate describing the proportion of false positives incurred within a set of significant features ( Storey and Tibshirani 2003 ). Furthermore, because the JML method assumes no recombination within loci, we confirmed the absence of recombination within each marker with a difference of sum of squares analysis ( McGuire and Wright 2000 ) using a sliding window of 80 bp and a step size of 10 bp, as implemented in the program TOPALi v.2.5 ( Milne et al. 2009 ).
STRUCTURE was run under the admixture model, to simultaneously estimate the number of populations in the sample and identify individuals with ancestry from more than one population. The model was run assuming numbers of clusters ( ) from to , using a burn-in of 20,000 followed by 100,000 MCMC samples. For each individual, the estimated proportion of ancestry from each cluster (q) and 95% confidence limits of this estimate were returned. Runs were repeated three times for each , and all other variables were left at their default parameters. The most likely number of clusters was inferred according to the rate of change in the log likelihood of data between successive values ( ; Evanno et al. 2005 ).
NewHybrids was informed by the results obtained with STRUCTURE and was run assuming two hybridizing groups and two generations of hybridization, resulting in six genotypic classes: two pure parental classes, F1 and F2 hybrids, and backcrosses in the direction of each parent. The program was run three times using different starting seeds, each time with a burn-in of 20,000 followed by 200,000 sweeps, default parameters for all other variables, and no reference genotypes provided. To investigate whether the employed set of microsatellite loci conferred sufficient power to discriminate between different hybrid classes, we also simulated a second-generation hybrid population using HYBRIDLAB v.1.0 ( Nielsen et al. 2006 ) and then used NewHybrids to assign simulated individuals. Eighty-eight individuals, inferred by STRUCTURE to have ancestry from a single cluster, were used as the F0 generation; 300 parental and F1 offspring were generated from these individuals, and these were in turn used to generate the F2 hybrid generation.
Morphometric Measurements and Variation Partitioning Analysis
Morphometric analysis was performed on formalin-fixed females. A selection of 23 clonal lines was used to provide the most efficient contrast between groups identified by mtCOI and nuITS1. Twenty, whenever possible, randomly picked individuals were measured from each clonal line. For each individual, microphotographs were taken under a LeitzLaborlux S optical microscope. Morphometric measurements were made using ImageJ (available online at: http:/rsbweb.nih.gov/ij/ —last accessed 10 January 2016). A total of 19 lorica traits were measured. These traits included those measured by Fu et al. (1991) , Ciros-Pérez et al. (2001) , and Proios et al. (2014) , with additional ones on the anterodorsal and anterovental sides (Supplementary Appendix 2: Section C, available on Dryad). Two traits of the anterodorsal side, namely “d” and “f,” were not included in further analysis due to high distortion of the placement of the anterior spines during preservation. Schematic representations of each of the measurements can be found in Supplementary Appendix 2: Figures S1–S4, available on Dryad. All the rotifer microphotographs analyzed are available in the Supplementary Material, available on Dryad.
To evaluate which species delimitation explained best the morphometric variation, we applied variation partitioning on redundancy analysis models (RDA). RDA is a linear regression technique designed to test the power, adjusted R 2 , of variables in explaining variation in a multivariate response variable matrix (here the morphometric variable matrix). The technique of variation partitioning ( Peres-Neto et al. 2006 ) allows one to estimate the unique explanatory contribution (conditional effect) of each explanatory variable as well as the amount of explained variation that it shares with the other explanatory variables in the model (collinear effect). We contrasted the performance of two types of delimitations (according to the consensus between all different methods): a nuclear-based delimitation based on the nuITS1 marker (the nuclear sequences of nu28S and nu18S showed lower levels or absence of genetic diversity, respectively), and a mitochondrial-based delimitation based on the mtCOI gene. Species delimitations in each case were coded as dummy variables, that is, each delimited species was represented in the explanatory matrix by a column in which each case (row) that corresponded to the respective species was coded “1” and the rest “0.” The significance levels of marginal and conditional effects were assessed with 999 random permutations. Given the distributional properties of morphometric data, we performed the analyses on untransformed data. To perform RDA and variation partitioning, we used the functions rda and varpart of the package “vegan” ( Oksanen et al. 2015 ) in R, respectively.
We tested for ecological differentiation between rotifers belonging to two nuITS1-delimited groups, namely “B” versus “C”, and used rotifers found in sympatry in the location coded “69” (Supplementary Appendix 1: Table S1, available on Dryad). Rotifer clones were genetically characterized for both the mtCOI and nuITS1 markers before the start of the experiment (Supplementary Appendix 1: Table S6, available on Dryad). The competition experiment was performed in semi-continuous batch cultures ( ) following a two-factorial randomized block design with stoichiometric food quality and dilution rate as experimental factors. To avoid the possibility that the competition outcome would be contingent on specific clonal combinations, we replicated each factorial combination eight times, with each replicate representing one of eight unique combinations of clones available from our cultures (cf. the blocks in the design; Supplementary Appendix 1: Table S6, available on Dryad). Food was derived from chemostat-grown phytoplankton ( Chlamydomonasreinhardtii ) and stoichiometric quality consisted of three levels: C:N:P = 33:4:1, C:N:P = 67:3:1, and C:N:P = 548:56:1, that is, nutrient-sufficient, nitrogen-limited, and phosphorus-limited food, respectively. These elemental ratios are within the natural range of freshwater habitats ( Elser et al. 2000 ). Two dilution treatments were tested, low dilution (4%.day −1 ) and high dilution (15–25%.day −1 ) (additional information in Supplementary Appendix 2: Section D.1, available on Dryad). The experiment lasted for 30 days.
At the end of the experiment, we randomly selected 10 rotifers from each culture and determined their nuITS1 group using the restriction endonuclease Dra I. On an agarose gel, Dra I digests the nuITS1 amplicon and generates fragments that can distinguish nuITS1 B from C rotifers (details in the Supplementary Appendix 2: Section D.2; Appendix 2: Table S2, available on Dryad). If both groups have equal competitive abilities, their relative abundances at the end of the experiment are expected to be equal. We tested for deviations from this one-to-one expectation using chi-square tests on the pooled data of each multifactorial combination separately. To test for an effect of food quality, of the dilution rate and of their interaction on the relative performance of both species, we applied a generalized linear mixed model (GLMM; e.g. Bolker et al. 2009 ) on the counts of both species in the experimental units, using a binomial error distribution with logit link function. In this model, food quality and dilution rates were specified as fixed factors and clone combination as random blocking factor. The analyses were performed in R using the “lme4” package ( Bates et al. 2015 ). To investigate whether hybrids had been produced during the experiment, we carried out microsatellite genotyping and mtCOI sequencing on all the rotifers from three (out of the eight) randomly selected experimental blocks, that is, starting clonal combinations (Supplementary Appendix 1: Table S6, available on Dryad).
We found 404 mtCOI haplotypes, 194 nuITS1 sequence types (with 49 heterozygous clonal lines among the 131 clonal lines inspected; Supplementary Appendix 1: Table S2, available on Dryad), nine nu28S sequence types, and a single nu18S haplotype for B. calyciflorus (Supplementary Appendix 1: Tables S2 and S3, available on Dryad). Because of this lack of polymorphism, the nu18S gene was not considered further. New sequences were deposited in GenBank (Accession numbers: KT729841-KT730043 for mtCOI, KT729547-KT729722 and KU364083-KU364144 for nuITS1, KT729748-KT729840 for nu28S, and KT729723-KT729747 for nu18S). Alignment length was 544 bp for mtCOI, 415 bp for nuITS1, and 461 bp for nu28S following trimming of low-quality nucleotide calls and, in the case of mtCOI and nuITS1, inclusion of the sequences downloaded from GenBank. Alignment files and data input and output files from all the analyses are available in the Supplementary Material, available on Dryad.
The best fitting models were found to be TVM+G for mtCOI ( Posada 2003 ), TPM3uf+G for nuITS1 ( Kimura 1981 ), and JC+I for nu28S ( Jukes and Cantor 1969 ). Because TVM and TPM3uf models are not presently available in MrBayes, the settings of the next best fitting models available were used: GTR+G for mtCOI ( Lanave et al. 1984 ; Tavaré 1986 ) and HKY+G for nuITS1 ( Hasegawa et al. 1985 ). Likewise, BEAST and *BEAST analyses were run using the GTR+G model for mtCOI, the HKY+G model for nuITS1, and the HKY+I model for nu28S.
Bayesian and maximum-likelihood phylogenetic inferences were very similar in their general topology and yielded good branch support (above or equal to 0.85 posterior probability or 75% bootstrap support) that largely corroborated the results of species delimitations ( Fig. 1 ). Different species delimitation methods yielded slightly different results, but there was a clear consensus of about 15 mtCOI haplotype groups, henceforth referred to as “1” to “15,” and of four main nuITS1 groups, labeled “A” to “D” ( Fig. 1 ). Even by the more conservative species estimates (shown as black bars in Fig. 1 ), the consensus estimate of species number remained essentially unchanged. Due to these findings, the upper limit for species number in the bGMYC analyses was set to species for mtCOI and for nuITS1 and nu28S. As expected from the literature ( Dellicour and Flot 2015 ), the bGMYC approach yielded a higher number of putative species than the single-locus GMYC approach performed on the BEAST tree ( Fig. 1 ). Among the methods tested, ABGD performed best on the mtCOI dataset, yielding delimitations in perfect agreement with the consensus delineation, but performed worst on the nuITS1 dataset (where it failed to detect any species). For the nuITS1 dataset the haploweb approach performed best, yielding a delimitation that was in perfect agreement with the consensus of the other approaches ( Fig. 1 ; Supplementary Appendix 2: Section E and Fig. S5, available on Dryad).
In the case of the nu28S gene, the very low levels of polymorphism—for example, just 8 or about 1.7% of the studied sites were parsimony informative—did not allow confident delimitations. For instance, GMYC-based delimitations yielded nonsignificant solutions ( ). We also did not find any heterozygous individuals to employ the haploweb approach. Regardless, distinct groups of nu28S haplotypes were recognized—each with fixed mutations—that corresponded to the delimited nuITS1 groups (Supplementary Material available on Dryad). As a result, the phylogenetic placement of the nu28S haplotypes in both Bayesian and maximum-likelihood trees totally matched the consensus delimitation scenario for the nuITS1 marker with Bayesian posterior probabilities 0.80 and likelihood bootstrap support (Supplementary Appendix 2: Section F; Appendix 2: Fig. S6, available on Dryad).
The microsatellite amplification pattern was also consistent with the nuITS1 delimitations (Supplementary Appendix 1: Table S7, available on Dryad). In practice, all the 12 microsatellite loci were amplified in nuITS1 C rotifers, while 9 of them were amplified in nuITS1 B rotifers (regardless of different mtCOI delimitations; Supplementary Appendix 1: Table S7, available on Dryad). For the rest of nuITS1 groups, “A” and “D,” the studied microsatellites seemed not to work; amplification and genotyping were difficult and inconsistent (Supplementary Appendix 1: Table S7; Appendix 2: Table S3, available on Dryad). For the nine co-amplified loci, levels of genetic diversity also supported a distinction between nuITS1 B and C rotifers (F ST = 0.533, from 1000 permutations; Supplementary Appendix 2: Fig. S7, available on Dryad). Allelic polymorphism was somewhat higher for nuITS1 C rotifers—but without accounting for sample heterogeneity (Supplementary Appendix 2: Table S4, available on Dryad), and group-specific alleles could be observed in each case (Supplementary Appendix 1: Table S7; Appendix 2: Section G, available on Dryad).
Using rotifer individuals that had been sequenced for both the mtCOI and nuITS1 markers, we identified many instances of mitonuclear discordance. In six of the 10 considered mtCOI-delimited groups, we observed coexistence of at least two distinct nuITS1-delimited groups ( Fig. 2 ). In the most extreme cases, rotifers of the mtCOI groups “8” or “15” were found to harbor sequence types from three different nuITS1 groups. Conversely, rotifers of, for example, the nuITS1 C group had mtCOI from seven different groups ( Fig. 2 ; Supplementary Appendix 1: Table S4, available on Dryad). In other words, discordance was widespread, occurring between several delimited groups across the mtCOI and nuITS1 phylogenies ( Fig. 2 ), even using more conservative delimitations ( Fig. 1 ).
Testing for Hybridization
Both DNA sequence- and microsatellite-based analyses supported the hypothesis that hybridization was a driver of the observed mitonuclear discordance. JML rejected the null hypothesis that ILS was solely responsible for the discordant pattern observed with for two out of the three pairwise comparisons for either the nuITS1 or the nu28S marker ( Table 1 ). Even in the remaining case in each marker, significance was quite high ( between “A” and “C” for nu28S) or hybridization was eventually supported otherwise [although between “B” and “C” for nuITS1, “B”/“C” hybrids were observed both in the wild samples and in the competition experiment using the microsatellites ( Fig. 3 b; Supplementary Appendix 2: Fig. S7, available on Dryad) or in the haploweb (Supplementary Appendix 2: Fig. S5, available on Dryad)]. At , the false positive rate was estimated to be low, at about 3% for nuITS1 and < 1% for nu28S. In the case of the B. plicatilis complex, JML did not reject the null hypothesis at for any of the species comparisons [but at and for two out of the 91 pairwise species comparisons (Supplementary Appendix 1: Table S8, available on Dryad)].
Note: Tests were conducted using posterior species tree simulations based on either the nuITS1 (lower diagonal) or the nu28S marker (upper diagonal). Cases with are shown in bold. At the false positive rate has been estimated at for nuITS1, and at for nu28S.
Because microsatellite loci amplified only in the nuITS1 B and C rotifers, STRUCTURE and NewHybrids analyses were limited to these two groups. We used the nine loci that could be amplified in both nuITS1 B and C (Supplementary Appendix 2: Tables S3 and S4, available on Dryad). Bayesian estimates of admixture proportions using STRUCTURE suggested two genotypic clusters, , as the most likely solution. The two clusters corresponded perfectly to each of the nuITS1 B- and C-delimited rotifers but were incongruent with the mtCOI delimitations, confirming once again the mitonuclear discordance (Supplementary Appendix 2: Fig. S7, available on Dryad). Three of the wild-sampled rotifers were confirmed as “B”/“C” hybrids (including a clonal line, 7C, identified as “B”/“C” hybrid in the haploweb; Supplementary Appendix 2: Fig. S5, available on Dryad) as the estimated 95% confidence intervals of ancestry did not overlap with either of the two parental groups (Supplementary Appendix 2: Fig. S7a, available on Dryad). One additional individual, sample 22BQ1, also appeared to have mixed ancestry, but wide confidence intervals precluded its definite identification as a hybrid. Using the same approach, five other “B”/“C” hybrids were detected at the end of the competition experiment (Supplementary Appendix 2: Fig. S7b, available on Dryad).
These findings were in complete agreement with the results from the NewHybrids analyses ( Fig. 3 ). NewHybrids assigned all but five wild-sampled rotifers with >95% probability to one of the two pure parental classes nuITS1 B or C; the remaining specimens were assigned with >30% probability to one or more of the hybrid classes ( Fig. 3 a). Five rotifers in the competition experiment were also assigned to at least one hybrid class, most frequently the F1 ( Fig. 3 b). In agreement with these findings, some of the hybrids of the competition experiment also produced a “hybrid restriction pattern” at nuITS1, that is, the Dra I enzyme produced fragments diagnostic for both “B” and “C” rotifers (Supplementary Appendix 2: Fig. S8, available on Dryad). In silico simulations with the HYBRIDLAB program indicated that, although the nine microsatellite loci could be used to positively identify hybrids among our samples, the precise mating events giving rise to these hybrid individuals could not be confidently inferred. The test with simulated data assigned all “pure” individuals with >80% probability to their correct parental class, with several of them also being partly assigned to the corresponding backcross class. However, the nine microsatellite loci were unable to assign simulated hybrids unequivocally to their correct hybrid class (F1, F2, or backcross).
Hybridization between “B” and “C” rotifers was also observed in the haploweb as a rare co-occurrence event of two otherwise abundant haplotypes characteristic of the “B” and “C” groups (Supplementary Appendix 2: Fig. S5; Appendix 1: Tables S2 and S4, available on Dryad). Flot et al. (2010) predicted that haplowebs could be used to allow hybrid detection, but to the best of our knowledge this is the first time that this is confirmed on an empirical dataset.
Each of the two kinds of delimitations, nuITS1- and mtCOI-based, significantly explained variation in the morphometric dataset ( Fig. 4 a). Analyzed separately, the mtCOI consensus delimitation explained 36% of the morphometric variation, while the nuITS1 consensus delimitation explained 71%. The effects of the mtCOI delimitations became insignificant when the nuITS1 delimitation was accounted for, whereas the conditional effect of the nuITS1 delimitation still amounted to 39%. All variation explained by the mtCOI delimitation was shared with the nuITS1 delimitation (34%). The above is also reflected in the results of a principal component analysis of which the two first axes represent 88% of the total observed morphological variation ( Fig. 4 b). None of the identified hybrids was included in the morphometric analysis (Supplementary Appendix 1: Table S9, available on Dryad).
In 11 of the 48 experimental units, rotifer populations had gone extinct before the end of the experiment. All remaining units proved strongly dominated by nuITS1 C rotifers. Compared to the abundance of nuITS1 B, the relative abundance of nuITS1 C averaged 96% and ranged between 67% and 100% (Supplementary Appendix 1: Table S10, available on Dryad). Therefore, chi-square tests showed very significant deviations from the expected 1:1 ratios for each of the 6 multifactorial combinations ( ). We found no significant effects of dilution rate, stoichiometric food quality, or the interaction of the two on the relative performances of the species (the GLMM was not significant).
As mentioned, Bayesian estimates of admixture proportions using the programs STRUCTURE and NewHybrids revealed five cases of admixed rotifers between nuITS1 C and D in the 30-day duration of the competition experiment ( Fig. 3 b; Supplementary Appendix 2: Fig. S7b; Appendix 1: Table S10, available on Dryad). Interestingly, the mtCOI identity of these hybrid rotifers was either “8” (two cases) or “10” (three cases). Based on the mtCOI identity of the rotifers combined at the start of experiment, “8” derived from nuITS1 B rotifers and “10” from nuITS1 C rotifers ( Fig. 3 b; Supplementary Appendix 2: Fig. S7b, available on Dryad) suggesting bidirectional hybridization.
In this study, we have shown the utility of integrative taxonomy to inform species delimitations in the B. calyciflorus cryptic species complex. This is despite a remarkable degree of mitonuclear discordance across species comparable only to a few other known cases ( Toews and Brelsford 2012 ), and limited expected phenotypic and ecological differentiation between sister Brachionus species ( Ortells et al. 2003 ; Fontaneto et al. 2007 ; Papakostas et al. 2013 ). By first conducting a comprehensive species delimitation analysis using different approaches (as shown in Fig. 1 ), we were able to account for the potential limitations of individual methods and to detect 15 mtCOI-based or 4 nuITS1-based putative species in our analyzed datasets ( Fig. 1 ). By then focusing on rotifer individuals sequenced for both markers, we demonstrated widespread discordance between the mtCOI and nuITS1 delimitations: for example, all three examined (out of the four identified) nuITS1-delimited species shared same mtCOI-delimited groups ( Fig. 2 ). Using morphological measurements and a competition experiment, we were able to demonstrate that nuITS1 delimitations were better predictors of morphological variation and of competitive abilities of rotifers than mtCOI delimitations (e.g., Fig. 4 ).
Finding evidence for hybridization as the driver of mitonuclear discordance in B. calyciflorus introduces a critical new dimension to the study of the evolution of this species complex. Hybridization may have varied impacts on evolutionary and speciation processes [as outlined in Barton (2001) and Abbott et al. (2013) ]. We provided three kinds of evidence in support of hybridization: first, JML tests rejected ILS as the sole driver of the discordance ( Table 1 ). Second, estimates of genetic admixture recognized many instances of hybrids both in the wild-derived samples and in the competition experiment ( Fig. 3 ; Supplementary Appendix 2: Fig. S7, available on Dryad). Importantly, observing hybrid genotypes over the 30-day course of the competition experiment provided convincing empirical evidence for the occurrence of bidirectional hybridization between nuITS1 B and C rotifers; mtCOI “8” haplotypes from nuITS1 B rotifers and mtCOI “10” haplotypes from nuITS1 D rotifers were both found in admixed individuals ( Fig. 3 b; Supplementary Appendix 2: Fig. S7b, available on Dryad). Third, hybridization between “B” and “C” rotifers was also supported by the nuITS1 haploweb (Supplementary Appendix 2: Fig. S5, available on Dryad). Notably, hybridization was not only supported by different approaches but also by different sources of molecular information: nuITS1 and nu28S sequences for the JML tests, microsatellite genotypes for the Bayesian admixture analyses, and nuITS1 sequences for the haploweb. NewHybrids analyses ( Fig. 3 ) suggested that hybridization in the wild populations had progressed further than the F1, although simulations indicated that more genetic information will be required to confirm these backcrosses.
We acknowledge that our use of the JML method can be disputed as we employed multi-copy markers, nuITS1 and nu28S, for the coalescent simulations. Due to gene conversion, multi-copy markers are prone to shorter coalescent times compared with single-copy genes ( Hillis and Dixon 1991 ; Hartfield et al. 2016 ), which may increase the false positive rate of the JML test. Brachionus rotifers are also facultative sexual organisms, and low rates of sex may enhance the effect of gene conversion on coalescent times ( Ceplitis 2003 ; Hartfield et al. 2016 ). For this reason, we validated our JML tests using nuITS1 sequences from a recently compiled dataset from the B. plicatilis complex, for which no mitonuclear discordance and thus no evidence of hybridization was found ( Mills et al. 2016 ). Assuming no strong differences (e.g., in effective population sizes) between B. plicatilis and B. calyciflorus , the lack of support for hybridization, at for any of the B. plicatilis species comparisons (Supplementary Appendix 1: Table S8, available on Dryad) suggests that gene conversion had a minor influence on the outcome of our JML tests. We also did not find any evidence for recombination within our studied markers (tested with the program TOPALi v.2.5—Materials and Methods section), which further suggests a minor effect of intra-locus gene conversion ( Wiuf 2000 ). In general, the baseline of gene conversion rate is low, for instance, at about 10 −5 to 10 −6 per site per generation for single-copy genes in the asexual bdelloid rotifer Adineta vaga ( Flot et al. 2013 ). However, as we cannot estimate the conversion rate for multi-copy markers in Brachionus rotifers, future work should also verify our findings using single-copy nuclear markers.
Assuming widespread hybridization, it is intriguing that we did not observe a significant breakdown of species boundaries in the B. calyciflorus phylogenies ( Fig. 1 ). Ecological specialization and meiosis suppression have been identified as mechanisms involved in generating and maintaining divergence in hybridizing cyclical parthenogenetic Daphnia species ( Cristescu et al. 2012 ; Xu et al. 2013 ). Different hypotheses predict barriers to gene flow between hybridizing species: restricted recombination of particular genomic regions in hybrids ( Noor and Bennett 2010 ) or, as in Daphnia , the occurrence of traits subject to divergent selection that would also contribute to some degree of reproductive isolation between species ( Servedio et al. 2011 ; Smadja and Butlin 2011 ). As such, B. calyciflorus may be a candidate model for future investigations on the controversial topic of reticulate evolution, particularly speciation, in the face of high levels of gene flow.
Perhaps one of the most notable findings of this study is that morphological variation in B. calyciflorus rotifers was best explained by the nuITS1 delimitations ( Fig. 4 ). Morphology has been known to conflict with mitochondrial delimitations in cases of introgressive hybridization ( Sullivan et al. 2004 ). Our results may have profound implications for the interpretation of morphological stasis and morphological plasticity−two components that are considered critical to understand cryptic species diversity ( Schlick-Steiner et al. 2007 ; Flot et al. 2011 ). It is striking that, based solely on mtCOI information and without knowing that nuITS1 delimitations are more accurate predictors of morphological variation, we would have overestimated levels of both morphological stasis and morphological plasticity in our samples. For example, we would have assumed high levels of morphological stasis between the two types of rotifers mtCOI “9” and “10”, were it not for the fact that, in this particular case, these rotifers belong to the same nuITS1 C group ( Fig. 4 b). Also, we would have assumed high levels of morphological plasticity among mtCOI “8” rotifers, were it not for the fact that, in this particular case, these rotifers belong to two different nuITS1 A or B groups ( Fig. 4 b). Admittedly, some of the studied cases were represented only by one clonal line, warranting further investigation, but the overall pattern is very strong and consistent ( Fig. 4 ). One would indeed expect that the larger size and number of genes of the nuclear genome exert a greater effect on the morphology than the mitochondrial genome. Intriguing testable hypotheses for further research include whether F1 hybrids are more similar to one of the parental species than to the other, and how the morphology of the hybrids is influenced by repeated backcrossing with one of the parental species ( Mallet 2005 ). As we were unaware of hybrid identity when we selected rotifer clonal lines for morphometric measurements, our selection of lines did not include any hybrids and these hypotheses remain therefore to be tested.
In conclusion, we have shown that integrative taxonomy is an extremely helpful framework to manage conflicting species delimitations in the challenging case of a rotifer cryptic species complex. By contrasting molecular-based species delimitations with information about the morphology and the ecology of the species, we were able to resolve mitonuclear discordance and to draw objective conclusions regarding the levels of morphological stasis and plasticity between species. The use of multiple nuclear markers—also single-copy—will still be needed, if not a genome-wide approach ( Seehausen et al. 2014 ), to fully understand the role of hybridization in the evolutionary history of the B. calyciflorus complex ( Mallet 2007 ; Abbott et al. 2013 ) and/or to identify what mechanisms maintain species integrity despite interspecific gene flow ( Kulathinal et al. 2009 ; Noor and Bennett 2010 ; Cruickshank and Hahn 2014 ; Krause and Whitaker 2015 ). Integrative taxonomy will likely be an important component of future speciation genomics studies ( Seehausen et al. 2014 ) aiming to better understand the mechanisms of speciation and of species diversity.
S upplementary M aterial
Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.8rc4r .
This work was supported by the Academy of Finland (grant number 258048). J.R. acknowledges funding by the Finnish Kone Foundation.
We thank Jean-François Flot for his critical help with the identification of heterozygous individuals and with the implementation of the haplowebs method, as well as for his suggestions that greatly improved this study. Frank Anderson and Robb Brumfield as well as an anonymous reviewer also provided constructive comments that further improved the manuscript. Matthieu Bruneaux helped with some of the bioinformatics work included in this article, Dennis Waasdorp provided technical assistance, and Pavlos Pavlidis helped understand better gene conversion and coalescence.