An upgraded comprehensive multilocus phylogeny of the Tardigrada tree of life

Providing accurate animals’ phylogenies rely on increasing knowledge of neglected phyla. Tardigrada diversity evaluated in broad phylogenies (among phyla) is biased towards eutardigrades. A comprehensive phylogeny is demanded to establish the representative diversity and propose a more natural classification of the phylum. So, we have performed multilocus (18S rRNA and 28S rRNA) phylogenies with Bayesian inference and maximum likelihood. We propose the creation of a new class within Tardigrada, erecting the order Apochela (Eutardigrada) as a new Tardigrada class, named Apotardigrada comb. n. Two groups of evidence support its creation: (a) morphological, presence of cephalic appendages, unique morphology for claws (separated branches) and wide‐elongated buccopharyngeal apparatus without placoids, and (b) phylogenetic support based on molecular data. Consequently, order Parachela is suppressed and its superfamilies erected as orders within Eutardigrada, maintaining their current names. We propose a new classification within the family Echiniscidae (Echiniscoidea, Heterotardigrada) with morphological and phylogenetic support: (a) subfamily Echiniscinae subfam. n., with two tribes Echiniscini tribe n. and Bryodelphaxini tribe n.; (b) subfamily Pseudechiniscinae subfam. n., with three tribes Cornechiniscini tribe n., Pseudechiniscini tribe n. and Anthechiniscini tribe n.; and (c) subfamily Parechiniscinae subfam. n., with two tribes Parechiniscini tribe n. and Novechiniscini tribe n. Reliable biodiversity selection for tardigrades in broad phylogenies is proposed due to biased analyses performed up to now. We use our comprehensive molecular phylogeny to evaluate the evolution of claws in the clawless genus Apodibius and claw reduction across the Tardigrada tree of life. Evolutionary consequences are discussed.

The main objective of this study is better understanding internal relationship within the Tardigrada phylogeny through a more comprehensive analysis. Secondary objectives will be: (a) evaluate monophyletic status from orders to genera considering classification changes, if needed; (b) provide tardigrade taxa selection for future metazoans' phylogenies; and (c) infer evolutionary traces of claws in the clawless genus Apodibius and claw reduction by means of the upgraded Tardigrada phylogeny.
Cycle sequencing with AmpliTaq DNA polymerase was as described by Guil and Giribet (2012). Cycle-sequenced products were cleaned using a standard protocol with ethanol, sodium acetate and formamide. The BigDye ® -labelled products were directly sequenced using an automated ABI PRISM 310 Genetic Analyzer. Chromatograms obtained from the sequencer were read, and contigs assembled using the sequence editing software SEQUENCHER version 4.1.4 (Gene Codes Corporation, Ann Arbor, MI). Assembled sequences were edited with BioEdit version 2007 (Hall, 1999), to identify fragments based on internal primers and conserved regions, as in a previous work (Guil & Giribet, 2012). All new sequences have been deposited in GenBank under accession numbers MH079453 to MH079475 for 18S rRNA, and MH079494 to MH079516 for 28S rRNA (Tables 1 and Supporting information Table S1).

| Phylogenetic analyses
We used available tardigrade sequences in GenBank, coincident with fragments analysed in the present study (Supporting information Table S1), to perform a more comprehensive analysis. We used four outgroups as in Guil and Giribet (2012) (Table 2). Disparity of genetic markers used for phylogenetic analyses of the Tardigrada phylum and taxa with those markers made us to perform three parallel analyses with: (a) 18S rRNA (fragment delimited by primers 18S a2.0 and 18S 9R), (b) 28S rRNA (fragment delimited by primers 28Sa and 28S 5b) and (c) a combined analysis with specimens where both genes, 18S and 28S, were successfully sequenced (Table 1).
Parallel analyses of maximum likelihood (ML) and Bayesian analyses (BI) were performed. Prior to likelihood analysis, jModeltest 2.1.1 (Darriba, Taboada, Doallo, & Posada, 2012) was executed to choose the best-fit model of nucleotide substitution for each gene (18S and 28S) and combined matrices, under the Akaike information criterion (AIC). For the 18S data set, the model 012343+I+G+F was obtained (with corrections for gamma distributions, proportion of invariable unchanging sites and the equilibrium base frequencies in the sequences are estimated by observing the occurrence in the data). For 28S, the model TIM2+I+G (transition model) T A B L E 1 (Continued) was resulted (with corrections for gamma distributions and proportion of invariable unchanging sites). Combined analyses were performed with partition data and their respective model described for each one. ML analyses were conducted using the program IQ-Tree (Nguyen, Schmidt, Haeseler, & Minh, 2015) in the web server version (https://iqtree.cibiv.univie.ac.at/), adapting model obtained with jModeltest. Nodal support was evaluated with 100 bootstrap replicates. BI was performed with MrBayes version 3.1.2 (Huelsenbeck & Ronquist, 2001;Ronquist & Huelsenbeck, 2003). Substitution model was specified in each case with parameters specifications as obtained with MrModeltest (Nylander, 2004) and separated models configured in combined analyses. Burn-in times were assessed by first running shorter analyses and graphing the Bayesian log likelihoods (LnL); these burnin times were subsequently confirmed by comparison with the complete log likelihood graphs of all analyses after 15,000,000 generations. Using Tracer version 1.5, burn-in times in a log likelihood graphs of all analyses were visualized, discarding 50,000 trees in each analysis. Support for nodes is expressed as posterior probabilities, calculated on a maximum clade credibility tree of the post-burn-in sample.

| R E S U LTS
We have sequenced 45 specimens from 26 taxa, obtained from moss and lichen samples collected in 18 localities widely distributed (Table 1 and Supporting information Table S1). This study included a large tardigrade diversity, as it covered over 80% of tardigrade families and subfamilies and 53% of genera (Table 3), making relevant conclusions achieved. Sequenced eutardigrade species represent all eutardigrade superfamilies, 92% of families and 59% genera (Table 3 and  Supporting information Table S1). Seven species and one genus (Austeruseus; Trygvadóttir & Kristensen, 2001) were newly sequenced for these molecular analyses (Table 1 and  Supporting information Table S1).
ML and BI analyses have been congruent between them irrespective genes used, being BI support stronger than ML bootstraps (Figures 1-3 for 18S; Figures 4,5 for 28S). Analyses combining 18S and 28S complete data sets agreed with analyses including one gene (18S or 28S) ( Figure 6). Information from the 18S rRNA solved nodes at different levels within the phylogeny (from classes to genera), while 28S rRNA solved deep (classes) and terminal nodes (genera and groups of genera) but not middle nodes. The two classes (Heterotardigrada and Eutardigrada) were supported with 18S, 28S and combined phylogenies, as well as eutardigrade orders Apochela and Parachela (Figures 1-6 Echiniscidae was supported by the three analyses (18S, 28S and combined), and order Echiniscoidea was only recovered with combined analysis (Figures 1-6 (Figures 1, 4 and 6). Mopsechiniscus remained in a doubtful position within the family Echiniscidae.
The family Milnesiidae (Apochela, Eutardigrada) showed two phyletic lines ( Figure 6): (a) Milnesium eurystomum (Spain) with Milnesium tardigradum from Denmark, Greenland and Spain, and (b) Milnesium tardigradum from Spain. Within parachelans, four phylogenetic lineages corresponding to superfamilies were supported (by 18S rRNA and combined analyses; Figures 2, 3 and 6): Hypsibioidea, Eohypsibioidea, Macrobiotoidea and Isohypsibioidea. At the level of parachelan superfamilies and families, 28S rRNA information showed no resolution ( Figure 5). The family Eohypsibiidae confirmed its monophyly incorporating a new genus, Austeruseus (Figures 3, 5 and 6). Within F I G U R E 1 Bayesian phylogram obtained with the nuclear 18S a2.0-9R data set (Supporting information Table S1). First number above branches is posterior probabilities obtained in the BI. Second number is bootstrap support values from ML. Taxa are named following Supporting information Table S1. Parachelan superfamilies are represented in detail in Figures 2-3. Classes, orders, families, superfamilies, genus and group of genera are indicated. Squares in different grey scales and dot limited squares highlight supported clades at different node levels. Scale bar = number of substitutions/site F I G U R E 2 Bayesian phylogram obtained with the nuclear 18S a2.0-9R for the superfamilies Hypsibioidea and Isohypsibioidea (Supporting information Table S1). First number above branches is posterior probabilities obtained in the BI. Second number is bootstrap support values from ML. Taxa are named following Supporting information Table S1. Orders, families, subfamilies genus and group of genera are indicated when monophyletic. Squares in different grey scales and dot limited squares highlight supported clades at different node levels. Scale bar = number of substitutions/site | 127 GUIL et aL.

F I G U R E 3
Bayesian phylogram obtained with the nuclear 18S a2.0-9R for the superfamilies Macrobiotoidea and Eohypsibioidea (Supporting information Table S1). First number above branches is posterior probabilities obtained in the BI. Second number is bootstrap support values from ML. Taxa are named following Supporting information Table S1. Orders, families, genus and group of genera are indicated when monophyletic. Squares in different grey scales and dot limited squares highlight supported clades at different node levels. Scale bar = number of substitutions/site Macrobiotoidea, three phylogenetic lineages can be detected corresponding to families with the combined analysis ( Figure  6): (a) Murrayidae, (b) Adorybiotus (maybe representing the family Richtersiidae, also supported with 28S; Figure  5) and (c) Macrobiotidae. Information with 18S rRNA only showed support to family Murrayidae but not to Richtersiidae (Richtersius and Adoryiotus) and Macrobiotidae (due to the inclusion of Eohypsibiidae; Figure 3). The family Macrobiotidae can be subdivided into four phyletic lines (   information Table S1).
First number above branches is posterior probabilities obtained in the BI. Second number is bootstrap support values from ML. Taxa are named following Supporting information Table S1. Classes, orders, families, genus and group of genera are indicated. Squares in different grey scales and dot limited squares highlight supported clades at different node levels. Scale bar = number of substitutions/site | 129 GUIL et aL. (Figures 1-6)

| Towards a natural classification of Tardigrada
The main purpose of Tardigrada phylogenies has been supporting, modifying or rejecting current tardigrade classification on the phylogenetic basis. We present a more comprehensive Tardigrada phylogeny, which reliability relies on the inclusion of 63 tardigrade genera out of the 119 described (Tables 3 and Supporting information Table S1).
The three classes within Tardigrada (i.e., Heterotardigrada, Mesotardigrada and Eutardigrada) were created at the beginning of the XX century, being Mesotardigrada questioned in several occasions (Grothman et al., 2017;Ramazzotti & Maucci, 1983). Eutardigrada monophyly has also been examined resulting dependent on the selection of outgroups for analyses (Guil & Giribet, 2012). In that study, the order Apochela was independent of class Eutardigrada. Table S1). First number above branches is posterior probabilities obtained in the BI. Second number is bootstrap support values from ML. Taxa are named following Supporting information Table S1. Orders, superfamilies, families, genus and group of genera are indicated. Squares in different grey scales and dot limited squares highlight supported clades at different node levels. Scale bar = number of substitutions/site Morphological differences among Tardigrada classes included: presence of appendages over the body, and morphology of claws and buccopharyngeal apparatuses (Bertolani, et al., 2014;Kristensen, 1987;Ramazzotti & Maucci, 1983). Heterotardigrada includes heterotardigrade (marine and terrestrial) claws and buccopharyngeal apparatus (Figure 7a,b) with a great variety of appendages in head and body, while Mesotardigrada shows heterotardigrade (Echiniscoidea) claws (Figure 7a), eutardigrade buccopharyngeal apparatus and cirrus A on head (Kristensen, 1987;Pilato & Binda, 2010;Ramazzotti & Maucci, 1983). Contrary, within the class Eutardigrada can be found claws and buccopharyngeal apparatuses of apochelan and parachelan types ( Figure  7c-f), while head appendages are present only in apochelans (peribuccal and cephalic papillae; Figure 8 and Schuster, Nelson, Grigarick, & Christenberry, 1980) (parachelans showed in some cases sense organs but not appendages). So, differences between orders Apochela and Parachela include head appendages and claw morphology used to differentiate classes within Tardigrada. In addition, phylogenetic evidences show strong support to class Heterotardigrada, and current orders Apochela and Parachela ( Figure 6). If considering class level as indicated in Figure 6, a new configuration with three classes (and doubtful Mesotardigrada) is evidenced as in other studies (Bertolani et al., 2014;Guidetti et al., 2009;Guil & Giribet, 2012). So, two groups of evidences support the creation of a new class for the current order Apochela: (a) a unique morphology for claws and buccopharyngeal F I G U R E 6 Bayesian phylogram obtained combining nuclear genes 18S rRNA and 28S rRNA data set (Supporting information Table S1).

F I G U R E 5 Bayesian phylogram obtained with the nuclear 28S a-5b for the class Eutardigrada (Supporting information
First number above branches is posterior probabilities obtained in the BI. Second number is bootstrap support values from ML. Taxa are named following Supporting information Table S1. Classes, orders, families, superfamilies, genus and group of genera are indicated. Squares in different grey scales and dot limited squares highlight supported clades at different node levels. New node level for classes proposed is indicated with a vertical line. Scale bar = number of substitutions/site apparatus (Figure 7e,f) together with the presence of cephalic appendages (peribuccal and cephalic papillae; Figure  8) and (b) molecular support from Bayesian and likelihood analyses with 18S rRNA and 28S rRNA information ( Figure  6). Consequently, we propose a new tardigrade class named Apotardigrada, following the former order name (Apochela) that indicates separate primary and secondary branches on claws. Within this new class Apotardigrada, the order Apochela is included, containing the family Milnesiidae, and genera and species composing this family as specified in Degma et al. (2018). Consequently, the class Eutardigrada diagnosis is amended excluding the cephalic appendages and claws with main and secondary branches separated. Since only parachelans remain within Eutardigrada, we propose to erect current superfamilies (Eohypsibioidea, Macrobiotoidea, Hypsibioidea, Isohypsibioidea) as orders within the class and suppression of order Parachela. Detailed taxonomic information is available in the Systematics section. Composition and diagnosis for former superfamilies (now orders) and families are as in Bertolani et al. (2014), Cesari et al. (2016), Guidetti et al. (2016) and Vecchi et al. (2016).
Here, we propose a new internal classification for the family Echiniscidae, with subfamilies and tribes (named after type genera) based on plates' presence and composition and shape of buccal sensory organs. We propose to create three subfamilies (Echiniscinae subfam. n., Pseudechiniscinae subfam. n. and Parechiniscinae subfam. n.) supported by molecular ( Figure 6) and morphological information based on the presence of pseudosegmental and neck plates (see Systematic section for details). Subfamily Echiniscinae subfam. n. is divided into two tribes on the basis of the shape of cirri A, external and internal buccal cirri and phylogenetic information with molecular data ( Figure  6): Echiniscini tribe n. and Bryodelphaxini tribe n. Three tribes organize internally the subfamily Pseudechiniscinae subfam. n. based on specific presence of pseudosegmental plates and phylogenetic support with molecular information ( Figure 6): Cornechiniscini tribe n., Pseudechiniscini tribe n. and Anthechiniscini tribe n. And two tribes are described within the subfamily Parechiniscinae subfam. n. on the basis of the presence of third median and/or head plate and phylogenetic support with molecular data (Figure 6): Parechiniscini tribe n. and Novechiniscini tribe n. Detailed taxonomic information, composition and diagnosis are available in the Systematics section.
In this sense, Apodibius inclusion, a clawless genus, within Isohypsibioidea ( Figure 6; Dabert et al., 2014) allows hypothesizing its claw evolution from an original isohypsibius claws from an isohypsibioidean ancestor until claw lost in current Apodibius. Claws' modification in the soil-dwelling Apodibius could be related to its association with soil and related environments, with tiny spaces between soil grains, where a worm-like shape would favour their movement. Hohberg and Lang (2016) (Figure 6), confirming a homoplastic evolution of the buccopharyngeal apparatus and its structures (Guil, Machordom, et al., 2013). Maybe, diversification to different feeding habits within distinct phylogenetic lineages, and so homoplastic evolution of the buccopharyngeal apparatus, can be related to guarantee of food roles execution within ecosystems (Guil & Sanchez-Moreno, 2013;Guil, Jørgensen, et al., 2013). These hypotheses, relating claw and buccopharyngeal apparatus evolution with ecology, open a new research line within tardigrades that need of further genetic, developmental, taxonomical and ecological information to be clarified. Tardigrada Doyère, 1840 Class Mesotardigrada Rahm, 1937 nomen dubium (diagnosis as in Ramazzotti &Maucci, 1983 andGrothman et al., 2017) Diagnosis: Cirri A present. With heterotardigrada-like spines. Heterotardigrade-like claws with no differentiation in main and secondary branches. Pharyngeal bulb with Eutardigrada-like macroplacoids. Class Apotardigrada (Schuster et al., 1980) Kristensen, 1987) Diagnosis: Tardigrada with cephalic, trunk and leg appendages. Gonopore separated from anus. Malpighian tubules lacking. Placoids consisting of three CaC0 3 elements or three delicate, bar-shaped cuticular structures. Composition: Taxonomic accounts and classification as in Kristensen, 1987, and Fontoura, Bartels, Jørgensen, Kristensen, & Hansen, 2017. Order Arthrotardigrada Marcus, 1927(classification as in Degma et al., 2018 Order Echiniscoidea Richters, 1926 (description as in Kristensen, 1987) Diagnosis: Heterotardigrada without toes on the legs. Median cirrus absent. Family Echiniscidae Thulin, 1928 (description as in Kristensen, 1987