Diversification rates in Tardigrada indicate a decreasing tempo of lineage splitting regardless of reproductive mode

Understanding the dynamics of speciation and extinction events is one of the most interesting subjects in evolutionary biology that relates to all life forms, even the smallest ones. Tardigrades are microscopic invertebrates that attracted public and scientific attention mostly due to their ability to enter into the diapause stage called cryptobiosis and in such stage resist extremely harsh environmental conditions. However, although recent research solved a considerable number of phylogenetic uncertainties and further uncovered physiological mechanisms of cryptobiosis, not much attention is given to the evolutionary forces shaping tardigrade diversity. Here, we investigated the effect of reproductive mode on diversification rates in tardigrades using three groups: macrobiotids, echiniscids and milnesids, which represent low, moderate and high levels of parthenogenesis, respectively. Our results showed a decreasing tempo of diversification events for each of the studied groups without any differences that could be ascribed to reproductive mode. We discussed the observed lack of effect in tardigrades acknowledging deficiencies in available data sets and encouraging further studies to understand whether our results can be considered reliable.


Introduction
Tardigrades are a phylum of ubiquitous microscopic invertebrates with a body size smaller than one millimetre. Although most of the 1380 currently recognized species have been found to live in mosses and lichens, there are also groups inhabiting exclusively marine or freshwater environments (Degma & Guidetti, 2007;Degma et al., 2021;Guidetti & Bertolani, 2005;Nelson et al., 2019). Tardigrades are mostly known for their ability to withstand harsh and unfavourable environmental conditions by entering cryptobiosis, a state in which measurable metabolic processes stop, preventing life functions to occur (Møbjerg & Neves, 2021). Along with sexual reproduction, a wide occurrence of asexual (parthenogenetic) reproduction is well known and documented in this phylum, whereas hermaphroditism was reported for only a few taxa (Rebecchi et al., 2000). Although there was an obvious increase in interest in studying tardigrades during the last three decades (mostly thanks to their cryptobiotic abilities), the most intensive studies towards describing their diversity burst about ten years ago. It was the time of the first reliable DNA taxonomy research in this group that showed a promise for enhancing diversity studies by utilizing genetic data (e.g. Cesari et al., 2009;Bertolani et al., 2011). From that point, genetic data in the form of DNA sequences started to accumulate in public data bases mostly thanks to the taxonomic studies conducted within the integrative taxonomy framework . Nevertheless, although each year brings another piece of the puzzle to clarify some systematic and phylogenetic uncertainties, not much attention has been given so far to the evolutionary forces shaping the observed taxa diversity.
Given that the two opposite types of reproduction are almost equally common within Tardigrada (sexual vs. asexual), and that the reproduction type may influence the diversification rates (Barraclough et al., 2003;Fontaneto et al., 2012;Tang et al., 2014), our working hypothesis was to check for the effect of reproductive mode on diversification rates in tardigrades. We focused on three large and different limno-terrestrial tardigrade clades that, besides differences in the proportion of sexual and asexual taxa within each of them, constitute also groups that have the largest coverage in terms of DNA sequences that are publicly available. The selected groups represent two classes distinguished within Tardigrada, the armoured Heterotardigrada Marcus, 1927 andthe soft-bodied Eutardigrada Richters, 1926. The first class is represented by the speciose family Echiniscidae Thulin, 1928 that currently groups about 300 herbivorous species within 21 distinct genera and is characterised by a moderate proportion (about 50%) of parthenogenetic taxa (Gąsiorek et al., 2019(Gąsiorek et al., , 2020(Gąsiorek et al., , 2021Piotr Gąsiorek pers. Com.). Eutardigrades are represented by two distinct evolutionary lineages belonging to different orders, namely the family Macrobiotidae Thulin, 1928 and the genus Milnesium Doyère, 1840. The family Macrobiotidae currently comprises about 300 omnivorous species grouped within 14 genera and is characterised by a high proportion of taxa (about 80%) with sexual reproduction (Bryndová et al., 2020;Stec et al., 2021a). The last group, the genus Milnesium, is the only representative of the order Apochela Schuster et al., 1980 for which genetic data are available. This order has been recently proposed as a separate class, Apotardigrada , not supported by following studies . Notably, the current systematics for some parts of the tardigrade tree is based on a descriptive interpretation more than a confidently confirmed scenario (Fleming & Arakawa, 2021). The genus Milnesium comprises 47 nominal species that are exclusively carnivorous and is characterised by an extreme level of morphological stasis with the presence of many putatively cryptic taxa (Morek et al., 2021;Roszkowska et al., 2016). This group practically constitutes an equivalent of a proper tardigrade family with a high proportion of taxa (about 80%) that reproduce parthenogenetically Morek et al., 2021;Witold Morek pers. com.).
We chose these three groups and reconstructed their phylogenesis in order to investigate their evolutionary rates as they (i) represent three very distinct evolutionary lineages within Tardigrada that (ii) differ in the proportion of sexually and asexually reproducing taxa (namely low, moderate and high level of parthenogenesis) and (iii) are characterised by good genetic data coverage that is publicly available. As not all sequences used in the study were properly assigned to valid taxa, we used modern methods of DNA taxonomy to delineate entities akin to species and at the same time checked their performance with large tardigrade data sets. Based on previous studies on rotifers comparing sexual and asexual taxa, we predicted that also tardigrades that reproduce predominantly via parthenogenesis would exhibit a higher accumulation of diversification events towards the root of the tree compared with their relatives reproducing mostly sexually.

Data set
The whole data set analysed in this study comprised sequences of a fragment of the cytochrome oxidase c subunit I gene (COI), belonging to three tardigrade clades: Echiniscidae, Macrobiotidae, Milnesium (Table 1). All these sequences were downloaded from GenBank and include at least part of the standard 658-bp-long Folmer's region of the COI gene. Genetic diversity within each of three groups is comparable: two groups represent families and the genus Milnesium has a vast amount of genetic diversity with apparent phenotypic stasis (Morek et al., 2021;Roszkowska et al., 2016), making also this group comparably diverse to the other two.
The degree of occurrence of parthenogenetic taxa within each group was catagorised based on published literature and anecdotal evidence from colleagues, namely with Echiniscidae assumed to have at least 50% of the species that are able to reproduce sexually (Gąsiorek et al., 2019(Gąsiorek et al., , 2020(Gąsiorek et al., , 2021Piotr Gąsiorek pers. com.), Macrobiotidae a high proportion of taxa, about 80%, with sexual reproduction (Stec et al., 2021a), and Milnesium being the group with asexuality being known in most of the species, about 80% (Morek &

Alignments preparations
The downloaded sequences were aligned using the AUTO method of MAFFT ver. 7 (Katoh & Toh, 2008;Katoh et al., 2002) and visually checked against potential problems such as gaps, ambiguously aligned positions in BioEdit ver. 7.2.5 (Hall, 1999) and stop codons in MEGA7 ver. 7.0 (Kumar et al., 2016). Only a few extremely short sequences and old sequences that contained gaps were eliminated from the data set. Then the alignment was used to calculate neighbour joining trees in MEGA7 to check for strangely behaving sequences. Such sequences were seen in Macrobiotidae (JX865312, MF503451-54) as well as Milnesium (FJ435810) data sets, and after BlastX confirmation of their affinity to different tardigrade groups, they were discarded from further analyses. The final alignments for Echiniscidae, Macrobiotidae and Milnesium comprised 548, 344 and 114 sequences, which were trimmed to 570, 624 and 569 bp, respectively. Some sequences contained IUPAC ambiguity codes (e.g. Y, R, S) that were converted in DnaSP (Librado & Rozas, 2009) to replace ambiguity codes by Ns. Then, we collapsed the alignments to haplotypes with DnaSP and used only these reduced datasets to reconstruct the phylogenies for species delimitation analysis. We used datasets with unique haplotypes as a more conservative solution to avoid biased and false coalescence processes caused by replicated haplotypes during the analyses. All DNA sequences converted in DnaSP and analysed in this study are given in supplementary materials (Online Resource 1).

Species delimitation and phylogenetic analyses
Since many complexes of cryptic or pseudo-cryptic taxa are present in tardigrades and their taxonomy is far from being satisfactory untangled (e.g. Grobys et al., 2020;Guidetti et al., 2019;Morek et al., 2021;Stec et al. 2018aStec et al. , 2020aStec et al. , b, 2021b we utilised four tools/approaches of DNA taxonomy to access a number of independent evolutionary entities, akin to species in our data sets. Specifically, we implemented two tree-based approaches, the Generalized Mixed Yule Coalescent (GMYC, Fujisawa & Barraclough, 2013) and the Poisson Tree Processes with a Bayesian upgrade (bPTP; Zhang et al., 2013) as well as two distance-based approaches, the Automated Barcoding Gap Discovery (ABGD; Puillandre et al., 2012) and the Assemble Species by Automatic Partitioning (ASAP; Puillandre et al., 2021). To delimit distinct evolutionary entities, tree-based approaches look for the transition point between speciation/extinction events and intraspecific coalescence processes whereas distance-based methods look for a barcoding gap within a given data set. For all phylogenetic analyses conducted in this study, an appropriate model of sequence evolution, as well as the best partitioning scheme, were chosen by using PartitionFinder ver. 2.1.1 (Lanfear et al., 2016) under the Akaike Information Criterion (AIC). The results of these searches are given in supplementary materials (Online Resource 2). As the input data for GMYC, we used ultrametric trees calculated in BEAST v2.5 (Bouckaert et al., 2019) through the CIPRES online portal (Miller et al., 2010) from the data sets collapsed to unique haplotypes. Parameters in BEAST were set as default, except for appropriate models from Partion-Finder2, relaxed clock log normal and 100,000,000 generations with a burn-in of 10%. In order to obtain the number of GMYC entities, the resulting ultrametric trees were loaded in R (R Core Team, 2020) and analyzed using the package 'splits' ver. 1.0.19 (Ezard et al., 2021). For bPTP analyses, we used Maximum-likelihood (ML) trees constructed using RAxML v8.0.19 (Stamatakis, 2014). The strength of support for internal nodes of ML construction was measured using 1000 rapid bootstrap replicates. Then the trees were used as input data for the bPTP performed on the web server (https:// speci es.h-its. org/) with default settings except for the number of MCMC generation that was set for maximum (500,000). As the result of the bPTP, we considered two numbers of recovered distinct entities that received the highest support in Maximum likelihood and Bayesian solution. For ABGD and ASAP analyses we used alignments with all sequences (not collapsed to haplotypes) to allow more reliable calculations for distance-based methods with sequences of unequal lengths (Fregin et al., 2012). The analyses were run on the respective servers (https:// bioin fo. mnhn. fr/ abi/ public/ abgd/; https:// bioin fo. mnhn. fr/ abi/ public/ asap/ asapw eb. html) with default settings, except for Milnesium ABGD analysis where X (relative gap width value) had to be lowered to 0.5 as the default settings returned only one partition and the lowering was recommended by the software. As the final output of all these delimitation approaches in each data set, each sequence was assigned to a specific independent evolutionary entity given by the specific approach. This information is provided in supplementary materials (Online Resource 1). To choose the final number of the distinct entities present in each data set and investigate the delimitation results, we visualised them as heatmaps with the R package 'heatmap3' v1.1.9 (Zhao et al., 2021). For each data set, all squares akin to species were positioned diagonally and their number was equal to the most conservative results from all tested delimitation approaches that was further considered by us as the species delimitation approach for further analyses. The heatmaps are provided in supplementary materials (Online Resource 3).
To avoid biased results by artificial inflation of the number of branching events towards the tips of the tree, for diversification rates analyses, we included only one representative sequence per each delimited species/entity from each of the three analysed data sets (Online Resource 1). From these reduced data sets, the final ultrametric trees for further diversification analyses were calculated in BEAST ver. 2.5 through the CIPRES online portal, as described above with prior search of the best model and partition scheme using PartionFinder2 (Online Resources 2 and 4).

Diversification rates analyses
To analyse shifts in diversification rates, we used the ultrametric trees constructed based on reduced data sets containing singular unique sequence per species/entity. Significant departures from the constant diversification rate model were tested following the rationale of Pybus and Harvey (2000). Their γ statistic compares the relative positions of nodes in phylogeny to those expected under a constant diversification rate model, under which the statistic follows a standard normal distribution. The situation in which nodes are closer to the tips than expected under the constant rate model (i.e. there has been an apparent increase in diversification rate toward the present) is reflected by positive values of γ. The opposite situation, when γ values are negative, signify an apparent deceleration that is visible as nodes (branching events) that are positioned more towards the tree root than the tips. Therefore, we obtained γ statistics for each of the three clades to test whether net diversification rates changed over time and how. Analysing diversification events with their shifts by using the γ statistics may be biased towards negative values if the phylogeny does not comprise all the species of a given group. In order to solve this problem, Pybus and Harvey (2000) developed the Monte Carlo constant rates test (MCCR). Their test accounts for taxonomic sample size and undersampling in the phylogeny (Fordyce, 2010) to recover estimates of whether the observed negative values are indeed significantly negative. As mentioned above, the high proportion of unresolved taxonomic issues in many tardigrade groups makes it impossible to assess the amount of undersampling that could be then implemented in the MCCR test. Nevertheless, despite the taxonomic obstacle, by using the distribution of the number of sequences that represent each delineated taxon in our data sets, we could predict the potential total number of species that should be present for each analysed group. To do that, we used the R package 'vegan' ver. 2.5-6 (Oksanen et al., 2020) to estimate the number of these potentially unseen species. The output is given as four estimations (Chao, first-order jackknife, secondorder jackknife and bootstrap) for each of the three analysed data sets (Online Resource 1). The γ statistics and the MCCR tests were computed by using the R package 'phytools' ver. 0.6.99 (Revell, 2012) whereas lineages-through-time plots to visualise diversification rates were obtained by using the R package 'ape' ver. 5.3 (Paradis and Schliep 2019). To increase the confidence of our MCCR test, the final tree for each of the analysed tardigrade group was tested twice-with the highest and the lowest estimated number of species that would maximise and minimise the bias caused by undersampling, respectively.

Statistical analyses
Since we analysed only three tardigrade groups, a proper statistical comparison of the eventual confounding factors would be unreliable. Therefore, we included data on sexual and asexual clades of rotifers from Fontaneto et al. (2012) in order to compare different diversification metrics and potential effects of confounding factors between dioecious and parthenogenetic groups. As our tardigrade groups cannot be easily classified as sexual vs. asexual, we made assumptions and considered Macrobiotidae and Echiniscidae as sexual groups whereas the genus Milnesium as asexual. For Macrobiotidae and Echiniscidae, we used the same rationale that is used for monogonont rotifers: males are known for at least half of the species of the group and sex with gene flow can thus be considered an actual drive of evolutionary diversification in these groups (Fontaneto et al., 2012). We used linear models (LM) in R to check the effect of differences between phylum (Tardigrada, Rotifera), the type of reproduction (sexual, asexual) and the number of entities/ species included in the analysis on the values of γ statistics. We tested also for interactions between predictors and the interaction terms were retained in the final model only if significant and improving model fit. Model fit was visually assessed with the R package 'performance' ver. 0.5.0 (Lüdecke et al., 2021).
Additionally, to assess the drivers of the difference in the number of delimited tardigrade species between tree-based and distance-based methods, we performed a paired t-test comparing mean values of the delimited entities for each of the approaches and each of the studied tardigrade group (see Online Resource 1 for raw input data).

Results
Out of the 1006 tardigrade COI sequences, 480 (ca. 48%) of them constituted unique and distinct haplotypes ( Table 1). The total number of tardigrade species/entities of the three analysed groups was 192, representing the most conservative solution of the DNA taxonomy delimitation methods (the lowest number of delimited entities per each group; Table 1). Importantly, paired t-test comparison indicated a significant difference between tree-based and distance-based methods, with tree-based methods demonstrating over-splitting tendency (t = −6.642, df = 2, p = 0.0219; Fig. 1).
The lineages-through-time plots obtained from the reconstructed phylogenies of the delimited species did not show any obvious difference between Echiniscidae, Macrobiotidae and Milnesium (Fig. 2). In all of the three analysed groups, diversification events seem to accumulate towards the root of the trees, as revealed by the significantly negative values of the γ statistics, supported also by the MCCR test to account for potential undersampling (Table 1 and Fig. 2).
The analyses performed including Tardigrada and Rotifera revealed that there is no significant influence of the phylum and the type of reproduction on the γ statistics (LM, phylum: t = −1.419, p = 0.199; type of reproduction: t = 0.948, p = 0.375). However, the number of species/entities had significant and negative effect on the γ statistics (t = −2.856, p = 0.0245; Fig. 3).

Discussion
The main result of our analyses is that we could not confirm our initial hypothesis that tardigrades that reproduce predominantly via parthenogenesis would exhibit a higher accumulation of diversification events towards the root of the tree than their relatives reproducing mostly sexually. All three analysed tardigrade groups with a low, moderate and high level of parthenogenesis (Macrobiotidae, Echiniscidae and Milnesium, respectively) revealed significantly negative values of the γ statistics suggesting initial high speciation rate followed by a decrease in the tempo of lineage splitting, without any differences that could be ascribed to reproductive mode.
Such results are in contrast to what is known in rotifers, where the putatively asexual bdelloids have different diversification rates than the cyclically sexual monogononts (Fontaneto et al., 2012;Tang et al., 2014). Rotifers and tardigrades live in the same habitats, from sediments of aquatic environments to limno-terrestrial interstices in mosses, lichens and soils; they also share microscopic body size, a variety of reproductive modes, desiccation tolerance, survival capabilities in the outer space and other unusual features (e.g. Jönsson et al., 2019;Nelson et al., 2019;Schill & Hengherr, 2019;Song et al., 2021;Wallace & Snell, 2001;Wallace et al., 2006). Notwithstanding the many similarities, the effects of reproductive mode on evolutionary trajectories are different in the two phyla, with apparently a strong effect in rotifers (Fontaneto et al., 2012) and a lack of effect in tardigrades (this study). Unfortunately, it is still premature to understand whether such lack of effect in tardigrades could be considered an actual biological pattern or if it simply reflects biases in the analysed data set. For example, whereas in the tested monogonont rotifers the occurrence of sexual reproduction has been found and described for all the taxa included in the analyses and for bdelloid rotifers the lack of sexual reproduction is assumed for all the analysed taxa (Wallace et al., 2006;Flot et al., 2012), in tardigrades none of the analysed clades can be considered as completely asexual or completely sexual. Only a gradient of occurrence of sexual reproduction could be inferred Fig. 1 Number of species/ entities delimited by tree-based and distance-based single-locus delimitation methods used in this study for three tardigrade groups: family Echiniscidae, family Macrobiotidae, genus Milnesium from the scattered literature (e.g. Altiero et al., 2019) and the lack of occurrence of males in the analysed populations. One might also argue that the simple distinction between sexual vs. asexual can be complicated in tardigrades considering hermaphrodites and heterogony (alternation of a dioecious with a parthenogenetic generation). Yet, both are extremely rare in Fig. 2 Bayesian phylogenetic reconstructions for three distinct evolutionary lineages of Tardigrada with respective lineage-through-time plots. The trees were calculated based on the reduced datasets with singular sequence representing a given species/entity delimited in this study with multiple DNA taxonomy approaches (see the "Mate-rial and methods" section for details). The tardigrade photos from the left to the right represent: Testechiniscus spitsbergensis tropicalis Gąsiorek et al., 2018, Macrobiotus shonaicus Stec et al., 2018band Milnesium variefidum Morek et al., 2016. Scale bar = 50 µm Fig. 3 The negative relation between values of γ statistics and the number of entities based on which they were calculated. Circles indicate Tardigrada; triangles indicate Rotifera; black indicates asexual reproduction; white indicates sexual reproduction tardigrades. Hermaphroditism has been documented only for six tardigrade species so far and out of which only two come from the family Macrobiotidae. These are two closely related pseudocryptic Macrobiotus taxa namely Macrobiotus hannae Nowak &Macrobiotus joannae Pilato &Binda, 1983. The remaining four hermaphroditic taxa (Bertolanius weglarskae (Dastych, 1972), Borealibius zetlandicus (Murray, 1907), Grevenius granulifer (Thulin, 1928) and Ursulinius lunulatus (Iharos, 1966)) belong to families that were not investigated in our study. Heterogony is not well documented and is even less frequent than hermaphroditism as it was reported only once for an unidentified Milnesium species (Suzuki, 2008). Although the rarity of these two unusual reproduction types indicates no influence on our study, the lack of effect in tardigrades could be due to our inability to clearly define monophyletic clades with unique reproduction modes, masking a potential effect of reproductive mode on evolutionary rates.
The confounding factor of overall diversity within each analysed clade seems to be the only significant driver of the estimated values of γ statistics (Fig. 3). Such a strong effect of sampling bias, even if accounted for in the calculations of the probabilities (Fordyce, 2010), may mislead the inference on the actual effect of differences in reproductive mode. Unfortunately, with the data at hand, we cannot assess the magnitude of such effect and additional data on the already analysed taxa and other taxa within rotifers and tardigrades is strongly needed to find an answer to the issue. We hope that our study and caveats of the currently available data will stimulate the community and other researchers to further test our or related hypotheses. There are still many interesting things awaiting their exciting discoveries in tardigrades and meiofauna in general, but for this better sampling and high resolution of integrated data will be crucial to produce confident answers.

DNA taxonomy in tardigrades
Nowadays, we are facing an extremely rapid accumulation of genetic data thanks to facilitated access to cost-and timeeffective sequencing technologies. On the other hand, we are dealing at the same time with a global biodiversity crisis triggered by the fast decrease of species diversity caused by environmental devastation and lack of human power that would describe this diversity before it disappears. Therefore, molecular species delimitation approaches may constitute a handy tool to taxonomists, ecologists and other biologists for accelerating species discovery and establishing the number of species in the sample through DNA taxonomy (Tautz et al., 2003). In our study, we conducted delimitation analyses using four methods commonly used in DNA taxonomy. Two of them are treebased approaches (bPTP, GMYC) whereas the other two are distance-based approaches (ABGD, ASAP). For the very first time, these methods have been tested with large tardigrade data sets enabling us to check their performance, namely we asked if distance-based and tree-based approaches might differ in the number of delimited entities. We found a significant difference between both approaches with three-based methods delimiting more species compared with distance-based methods. The results are in line with previous studies that systematically evaluated the performances of these methods (e.g. Dellicour & Flot, 2018;Magoga et al., 2021). It has been shown that generally, distancebased methods outcompete tree-based methods giving more conservative solutions with a lower number of delimited taxa that better match morphologically diagnosed species whereas tree-based methods have over-splitting tendencies (Magoga et al., 2021). Although single-locus delimitation approaches are quite straightforward, easily accessible and time effective, there is no single ideal approach and the general recommendation for integrative approaches advocated for other taxa (Dayrat, 2005;Dellicour & Flot, 2018;Padial et al., 2010) applies also to tardigrades. The use of several methods, looking for the best congruency between them covering data coming from DNA sequences, morphology, ecology, geography, etc., should become a common approach in tardigrade taxonomy.

Acknowledgements
We are especially grateful to our colleagues Piotr Gąsiorek and Witold Morek for sending us tardigrade microphotographs and consultation on tardigrade groups in which they are specialists as well as to Matteo Vecchi for his valuable advice in handling DNA data with R. Finally, we are greatly indebted to Krzysztof Miler and an anonymous reviewer for their precious comments and improvements to the earlier draft of our manuscript. During this study, DS was supported by the Foundation for Polish Science (FNP). Data availability All data generated or analysed during this study are included in this published article and its supplementary information files.

Declarations
Ethics approval Not applicable.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.