An integrated study of the biodiversity within the Pseudechiniscus suillus–facettalis group (Heterotardigrada: Echiniscidae)

Pseudechiniscus is the second most species-rich genus in Heterotardigrada and in the family Echiniscidae. However, previous studies have pointed out polyphyly and heterogeneity in this taxon. The recent erection of the genus Acanthechiniscus was another step in making Pseudechiniscus monophyletic, but species identification is still problematic. The present investigation aims at clarifying biodiversity and taxonomy of Pseudechiniscus taxa, with a special focus on species pertaining to the so-called ‘suillus–facettalis group’, by using an integrated approach of morphological and molecular investigations. The analysis of sequences from specimens sampled in Europe and Asia confirms the monophyly of the genus Pseudechiniscus. Inside the genus, two main evolutionary lineages are recognizable: the P. novaezeelandiae lineage and the P. suillus–facettalis group lineage. Inside the P. suillus– facettalis group, COI molecular data points out a very high variability between sampled localities, but in some cases also among specimens sampled in the same locality (up to 33.3% p-distance). The integrated approach to the study of Pseudechiniscus allows confirmation of its monophyly and highlights the relationships in the taxon, pointing to its global distribution.


INTRODUCTION
The phylum Tardigrada consists of more than 1200 described species (Guidetti & Bertolani, 2005;Degma & Guidetti, 2007;Degma et al., 2018), subdivided into three classes: Eutardigrada, Mesotardigrada and Heterotardigrada (Ramazzotti & Maucci, 1983). The latter class is mainly made up of semi-terrestrial or marine species, and it is subdivided into two orders: Arthrotardigrada and Echiniscoidea (Marcus, 1929). Inside Echiniscoidea, the second most species-rich genus is Pseudechiniscus Thulin, 1911, surpassed only by the genus Echiniscus C. A. S. Schultze, 1840 Gąsiorek et al., 2018c). However, Pseudechiniscus has always been controversial. For example, Maucci (1973-74) and Dastych (1984) suggested a complete revision of the genus. It was finally emended by Kristensen (1987), but two main groups inside the genus were still distinguished: the P. suillus/conifer group, characterized by the absence of trunk cirri; and the P. victor group, identifiable with the presence of trunk cirri. Moreover, the advent of molecular investigations quickly pointed to polyphyly in Pseudechiniscus Guil & Giribet, 2012;Guil et al., 2013). A recent integrative analysis (Vecchi et al., 2016) has led to the erection of the genus Acanthechiniscus, more or less corresponding to the old P. victor group, another step in making Pseudechiniscus monophyletic, with the genus now comprising 40 described species .
T h e h i g h n u m b e r o f d e s c r i b e d s p e c i e s notwithstanding, species identification is still difficult in Pseudechiniscus (Kathman & Dastych, 1990;Fontoura & Morais, 2011). Some of the species are difficult to distinguish from one another, which is due to different issues, such as older, imprecise or incomplete diagnoses (for actual standards) with omissions of important characters, probably linked to the use of poor equipment and, in some cases, the type material being absent or lost (Fontoura & Morais, 2011). Moreover, some characters are subject to different interpretations, e.g. faceted cephalic and terminal plates, fine or coarse dorsal granulation and small cones on lateral positions. In addition, there is still some confusion between P. suillus (Ehrenberg, 1853), the type species, and P. facettalis Petersen, 1951, especially concerning the subdivision of the pseudosegmental plate and the faceting on the terminal plate, a highly variable character (Kathman & Dastych, 1990;Fontoura & Morais, 2011). Interestingly, even in the original description, Petersen (1951) described P. facettalis as a form of P. suillus, stating 'the differences are not so great that it can be apprehended as an independent species'.
The present investigation aims at studying the diversity and taxonomy of Pseudechiniscus taxa, using an integrated approach of morphological (light and scanning electron microscopy) and molecular (18S, 28S and COI genes) investigations in order to elucidate relationships in the genus Pseudechiniscus, with a special focus on species pertaining to the P. suillusfacettalis group.

MATERIAL AND METHODS
Mosses were sampled in five different localities in Eurasia (Table 1). After soaking a fragment of the moss for half an hour, tardigrades were extracted from samples by repeatedly washing and squeezing the fragment through consecutive 500-μm and 38-μm sieves. Tardigrades pertaining to the genus Pseudechiniscus were individually isolated using a needle and a glass pipette under a stereomicroscope. The remaining fractions of the samples are stored at -80 °C in the Biobank of the Evozoo Lab at the Department of Life Sciences (University of Modena and Reggio Emilia, Italy -UNIMORE).
Tardigrades (53 specimens) were mounted in Faure-Berlese mounting medium for observations under light microscopy (LM), while the other 28 animals were prepared for scanning electron microscopy (SEM), following the protocol described in Bertolani et al. (2011b). Observations with LM were carried out with both phase contrast (PhC) and differential interference contrast (DIC) up to the maximum magnification (objective 100× with oil immersion) with a Leica DM RB microscope, equipped with a digital camera Nikon DS-Fi 1, at UNIMORE. Observations with SEM were carried out with a Nova Nano SEM 450 (FEI Company -Oxford Instruments), available at the 'Centro Interdipartimentale Grandi Strumenti' at UNIMORE. All slides are deposited in Bertolani's collection at the Department of Life Sciences, UNIMORE. Unfortunately, for all studied populations, morphometric data according to Fontoura & Morais (2011) were impossible to collect, because most specimens on slides were in an unsuitable position for the measurement of the scapular plate. The latter was almost never in only one focal plane, preventing its measurement. Moreover, in the few slides in which the scapular plate was completely visible in only one focal plane, the specimens were so flattened that the edge of all the other dorsal plates were not identifiable and, therefore, impossible to measure.
A total of 27 specimens were utilized for molecular analyses (Table 1; Supporting Information, Table S1). Before total DNA isolation, single Pseudechiniscus specimens were temporally mounted in water on a slide, and observed and identified at LM up to 100× magnification, following the protocol described in Cesari et al. (2011). Voucher images of taxonomically relevant structures (i.e. cuticle, papillae, spines, buccalpharyngeal apparatus and claws) were recorded for each specimen and are available upon request. Genomic DNA was extracted from single tardigrades using the MasterPure Complete DNA and RNA Purification kit (Epicentre -Illumina, Madison, WI, USA), following the manufacturer's protocol. A region of the nuclear ribosomal small subunit gene (18S) was amplified with the primers and protocols described in Vicente et al. (2013), whereas a region of the nuclear ribosomal large subunit gene (28S) was amplified with the primers and protocols as described in Guidetti et al. (2014). A portion of the mtDNA cytochrome c oxidase I gene (COI) was carried out using primers LCO (5'-GGT CAA CAA ATC ATA AAG ATA TTG G-3′) and HCOout (5´-CCA GGT AAA ATT AAA ATA TAA ACT TC-3'), with the following protocol: initial denaturation at 94 °C for 5 min, followed by five cycles of 94 °C for 1 min,  Table S1). For 18S and 28S gene analysis, nucleotide sequences were aligned with the MAFFT algorithm (Katoh et al., 2002), as implemented in the MAFFT online service (Katoh et al., 2017) and checked by visual inspection. Sequences pertaining to a Florarctus sp. (Heterotardigrada, Arthrotardigrada; GenBank acc. nos.: GQ849017 for the 18S gene, GQ849034 for the 28S gene; Jørgensen et al., 2010) specimen were used as outgroup. Other Echiniscoidea sequences from GenBank were also included in the analysis for appropriate comparisons (Supporting Information, Table S2). A Bayesian inference dendrogram was computed with the program MrBayes v.3.2.6 (Huelsenbeck & Ronquist, 2001;Ronquist & Huelsenbeck, 2003) as implemented in CIPRES v.2.0 (Miller et al., 2010). Best-fitting model evaluations were performed taking into account the Corrected Akaike Information Criterion (jModelTest v.2.1.10; Guindon & Gascuel, 2003;Darriba et al., 2012), which identified the TPM1+I+G model for the 18S gene and the GTR+G model for the 28S gene as the most suitable. The TPM+I+G model was replaced by the GTR+I+G model in the MrBayes' analysis. Two independent runs, each of four Metropolis-coupled Markov chains Monte Carlo method, were launched for 30 × 10 6 generations and trees were sampled every 1000 generations. Convergence of runs was assessed by tracking average standard deviation of split frequencies between runs and by plotting the log likelihood of sampled trees in TRACER v.1.5 (Rambaut & Drummond, 2007) and the first 3 × 10 6 sampled generations were discarded as burn-in. A maximum likelihood (ML) analysis was performed with the program RAxML v.7.2.4 (Stamatakis, 2006), as implemented in CIPRES, using the models described above. Bootstrap resampling with 1000 replicates was undertaken via the rapid bootstrap procedure of Stamatakis et al. (2008) to assign support to branches in the ML tree.
For COI gene analysis, chromatograms obtained were checked for the presence of ambiguous bases. Sequences were translated into amino acids by using the invertebrate mitochondrial code implemented in MEGA7 (Kumar et al., 2016) in order to check for the presence of stop codons and, therefore, of pseudogenes. Nucleotide sequences were aligned with the MAFFT algorithm, as described above, and checked by visual inspection. For appropriate molecular comparisons, COI sequences from GenBank pertaining to other Echiniscidae specimens (Supporting Information, Table S2) were considered in the analysis. Pairwise nucleotide sequence divergences between scored haplotypes were computed using p-distance with MEGA7. The relationships among COI haplotypes were estimated using a parsimony network by applying the method described by Templeton et al. (1992), as implemented in TCS v.1.21 (Clement et al., 2000) and visualized using tcsBU (Santos et al., 2015). A 95% connection limit was used as a useful general tool in species assignments and discovery (Hart & Sunday, 2007). Putative species were also inferred by using the Poisson Tree Process (PTP; Zhang et al., 2013) and the Automatic Barcode Gap Discovery method (ABGD; Puillandre et al., 2012). PTP is a coalescentbased species delimitation method that use nonultrameric gene trees as input, and utilizes heuristic algorithms to identify speciation events relative to numbers of substitutions. The PTP method produces robust diversity estimates, in some cases more robust than those estimated under the generalized mixed Yule coalescent model (Tang et al., 2014). The starting gene tree was a ML tree computed using RAxML under the GTR+G model with the same procedure as described above for the 18S gene. In the distancebased ABGD method, the sequences are sorted into hypothetical species based on the barcode gap (i.e. whenever the divergence among organisms belonging to the same species is smaller than divergence among organisms from different species). The method first detects the barcode gap as the first significant gap beyond a model-based one-sided confidence limit for intraspecific divergence, and then uses it to partition the data. ABGD settings for the COI dataset were: prior minimum divergence of intraspecific diversity (P min ) = 0.001; prior maximum divergence of intraspecific diversity (P max ) = 0.1; Steps = 10 and gap width = 1.5. The analysis was performed on the ABGD website (http://wwwabi.snv.jussieu.fr/public/abgd/ abgdweb.html).

RESULTS
The morphological analyses with LM and SEM of all investigated Pseudechiniscus specimens evidence the presence of pedal (leg) plates, the absence of the spine on the first leg pair and a dome-like cephalic papilla (Figs 1A, B, 2A, B, 3A, B, 4A, B). The latter character points to their belonging to the 'suillus-facettalis' species group (Petersen, 1951). Although the cephalic and caudal plates' faceting should allow a clear attribution of the specimens to either the P. suillus (no plate faceting) or to the P. facettalis (plate faceting present; Ramazzotti & Maucci, 1983;Gąsiorek & Degma, 2018), this character is not clear and does not help in unambiguously identifying specimens (Figs 1C, D, 2C, D, 3C, D, 4C, D). All sampled localities feature male and female specimens. The female gonopore is made up of six cuticular folds that appear as a roselike structure (Fig. 1E, F), while the male gonopore is round (Fig. 3E, F). These data point to the bisexual condition of analysed specimens.
The specimens pertaining to the Slovak sample (C3019; Fig. 1) have a peculiar ventral sculpture with respect to all other analysed specimens. The sculpture is made up of the heads of the intracuticular pillars that appear as dotted areas forming a reticulate pattern (Fig. 1E, F). This peculiar reticulation is very similar to that of Pseudechiniscus xiai Wang et al., 2018. The only differences found between these species are the division of median plate 2 in the Slovak sample (undivided in P. xiai) and the not uniform terminal plate in the Slovakian specimens, which presents notches, as depicted in fig. 10, p. 45 of Petersen (1951). In P. xiai, the terminal plate is without notches, nor faceted.
All other analysed specimens (Figs 2-4) present no evident differences in the cuticle ventral sculpture, which is made up of the heads of the intracuticular pillars forming a fine uniform granulation that becomes coarser between the legs and in correspondence of the gonopore (Figs 2E, F, 3E, F, 4E, F).
The phylogenetic inference based on 1652 bp of the 18S + 28S genes (Fig. 5) shows Echiniscoides branching of the base with respect to a group with all other Echiniscoidea taxa. Inside the latter cluster, the Parechiniscus genus is sister to five main lineages: the Hypechiniscus + Testechiniscus + Diploechiniscus + Echiniscus cluster, three different lineages representing the genera Mopsechiniscus, Oreella and Bryodelphax, respectively, and a large group containing four other genera. Inside this latter cluster, two main lineages can be observed: one groups Proechiniscus, Cornechiniscus and Acanthechiniscus, while the other clusters all Pseudechiniscus specimens. The Pseudechiniscus lineage is subdivided into two highly supported clusters: one groups a sequence from GenBank misidentified as Echiniscus sp. with Pseudechiniscus titianae Vecchi et al., 2016 (characterized by an elongated cephalic papilla, typical of the Pseudechiniscus novaezaelandiae group), while all sequences grouped inside the second cluster belong to specimens with a dome-like cephalic papilla, which characterizes the species pertaining to the P. suillus-facettalis group. This second cluster is divided into two main groups: (1) one sequence from a Mongolian specimen belonging to the first cluster, together with P. facettalis specimens from Greenland and (2) all other sequences obtained from newly analysed specimens belonging to a highly supported group comprising also Norwegian and Spanish sequences from GenBank, with the sequences from Slovakian specimens sister to all. No clear-cut division between P. suillus and P. facettalis is found; both species are comprised in the latter cluster.
The analyses of 675 bp of the COI gene of the Pseudechiniscus specimens show a high variability. Ten new haplotypes are found among 15 specimens. The genetic p-distances between scored haplotypes (Supporting Information, Table S3) are in the range of 0.2-33.3%. Inside each sampled location, p-distances are generally low (Portugal: 0.0-1.7%; Slovakia: 0.0-0.2%; Italy: 0%), with the exception of both Mongolian samples, where a very high variability (up to 33.3% p-distance) is found among specimens. The PTP analysis (Fig. 6) shows nine putative species clusters: one each for Portuguese (together with Spanish P. facettalis specimens from GenBank), Italian, Slovakian and Greenlandic specimens. On the other hand, the Mongolian specimens are subdivided into five different putative species. This subdivision is confirmed by both the ABGD and the haplotype network analysis (Fig. 6).

DISCUSSION
The integrated approach to the study of Pseudechiniscus allows us to confirm and clarify its taxonomic position. Our new data confirm the subdivision between Acanthechiniscus and Pseudechiniscus identified by Vecchi et al. (2016) and also points out that all specimens sampled in European, Asian, Arctic and Antarctic regions belong to Pseudechiniscus. This highlights the global distribution of Pseudechiniscus, differing from what was recently discovered for another echiniscid genus (i.e. Testechiniscus; Gąsiorek et al., 2018a). While all Pseudechiniscus specimens clustered together in a well-supported group, their relationships with other members of the recently erected subfamily Pseudechiniscinae (Guil et al., 2019) appear less evident. In fact, our analysis finds good support for the Pseudechiniscus + Proechiniscus + Cornechiniscus + Acanthechiniscus lineage, but it does exclude all Mopsechiniscus specimens from this evolutionary lineage, also negating the tentative including of this genus inside the Pseudechiniscini tribe (Guil et al., 2019). Therefore, more in-depth analyses will be needed in order to confirm or disprove this new Echiniscidae systematics.
Two species-groups were identified inside Pseudechiniscus, recognizable by both morphological and molecular approaches, i.e. the P. novaezealandiae group (characterized by an elongated cephalic papilla) and P. suillus-facettalis group (characterized a domelike cephalic papilla; Fig. 5). This datum is also confirmed by a morphological phylogenetic analysis by Gąsiorek et al. (2018c), in which P. novaezealandiae and P. suillus are two clearly separated lineages. These species groups could be described as new subgenera, but we think that more integrative data are necessary in order to erect new taxa, especially for the P. novaezealandiae group.
However, no clear-cut subdivision could be found between P. suillus and P. facettalis specimens. Even though some previous authors (Ramazzotti & Maucci, 1983;Gąsiorek & Degma, 2018) used the presence of faceting in the terminal plate to discriminate between these two species, this character was not detectable and, therefore, it was not possible to unambiguously identify specimens. Considering that also Dastych (1988) synonymized P. suillus and P. facettalis, the taxonomic status of these two species is not well defined and needs revision.
The integrative analysis allows us a clear-cut identification of the specimens from Slovakia. They all have a peculiar ventral sculpture, characterized by a reticulated pattern (Fig. 1E, F), similar to that exhibited by P. xiai, sampled in the Liupan Mountains (Central China;Wang et al., 2018). The molecular results of Slovakian specimens clearly point to a welldefined taxon (Fig. 6) but, unfortunately, no molecular data are presently available for the type population of P. xiai. Given the large geographical distance and the small morphological differences between Slovakian and Chinese specimens (divided median plate 2, plate IV slightly faceted in some individuals), presently there are no definitive elements to assign the Slovakian specimens to P. xiai or to a new species.
The present study reveals the presence of semicryptic species in the genus Pseudechiniscus (Fig. 6). At least seven different putative species are found: while for some of them there is at least a geographical subdivision (e.g. Italy and Iberian Peninsula), in other cases four and two different putative species are found in the same gathering (Mongolia C2592 and Mongolia C2595, respectively), with one putative species shared between the two sampling sites. The presence of cryptic species in tardigrades has been known since Faurby et al. (2008) found high genetic distances in a single morphotype of Ramazzottius oberhaeuseri (Doyère, 1840). Even though some criticisms have recently arisen about the cryptic status of these molecularly identified species , the fact that taxonomically important morphological characters are scarce in tardigrades must not be underrated. Actually, the application of an integrated approach has allowed the definition of many tardigrade species (Cesari et al., 2009(Cesari et al., , 20112016b;Bertolani et al., 2011a, b;Guidetti et al., 2013Guidetti et al., , 2014Vicente et al., 2013;Stec et al., 2015Stec et al., , 2017aStec et al., , b, c, 2018bGąsiorek et al., 2016Gąsiorek et al., , 2017aGąsiorek et al., , b, c, 2018bMøbjerg et al., 2016;Morek et al., 2016;Zawierucha et al., 2016Roszkowska et al. 2017Roszkowska et al. , 2018Buda et al., 2018;Kaczmarek et al., 2018;Nowak & Stec, 2018;Perry et al., 2018), but it also pointed out that many difficult taxonomic situations still exist in both Heterotardigrada and Eutardigrada, and it evidenced the surprisingly large presence of putative cryptic species in the phylum (Guidetti et al., , 2016 Bertolani et al., 2011b;Faurby et al., 2011Faurby et al., , 2012Faurby & Barber, 2015;Cesari et al., 2016a;Stec et al., 2018a;Morek et al., 2019). Presently, the most sensible solution would be to place species difficult to distinguish into the category of cryptic or semicryptic species (complexes that may display minor, but still detectable, morphological differences; Korshunova et al., 2017) until new diagnosable characters are found (Bickford et al., 2007). Evolutionary lineages identified only by molecular methods should not be described as new species, but temporarily be considered unidentified candidate species (UCS; Padial et al., 2010). This will allow the flagging of potentially delicate situations and the call for the use of different approaches (e.g. karyological, ecological, ethological, etc.) and possible new genetic markers. The present taxonomic situation in the P. suillus-facettalis group appears intricate due to several factors: (1) the original description of P. suillus is dated and incomplete, (2) the discrimination between P. suillus and P. facettalis is problematic and (3) recent morphological (e.g. ventral cuticular pattern) and morphometric characters are still not applied to a large number of taxa in order to be reliable for a clear intra-and interspecific discrimination. Therefore, considering the low morphological variability and the lack of morphometric data in our studied populations, we designate eight Pseudechiniscus (semi)cryptic species as UCS. Following Padial et al. (2010), these UCS should be defined with the combination of the binomial name of the most similar species, followed (in square brackets) by the abbreviation 'Ca' (for candidate) with an attached numerical code referring to the particular candidate species and terminating with the author name and year of publication of the article (i.e. this paper) in which the lineage was first discovered. In order to be more precise, we also have decided to designate these lineages as part of a species group by adding the affinis designation (aff.). One UCS (Portuguese specimens) is designed as P. aff. facettalis because the obtained sequences match with four GenBank sequences assigned to P. facettalis (FJ435811-2; JX683830-1). Another UCS (Slovak specimens) is attributed to P. aff. xiai because it is the most similar described species, even though some small morphological differences are blue-coloured branches to red-coloured branches. Newly scored haplotypes are in bold. The scale bar shows the number of substitutions per nucleotide position. Centre: orange rectangles represent specimens grouped by ABGD analysis. Right: haplotype network of COI gene in Pseudechiniscus. Circles represent haplotypes, while circle surface denotes haplotype frequency. Small black squares indicate missing/ideal haplotypes. Networks falling below the value of the 95% connection limit are disconnected. Grey rectangles delimitate putative species. T h e i n t e g r a t e d a p p r o a ch t o t h e s t u d y o f Pseudechiniscus allows confirmation of its monophyly and highlights the relationships in the taxon. The present study also, for the first time, reports the presence of semicryptic species in the P. suillus-facettalis group, and underlines the taxonomic problems related to this group. New studies will be performed in the future to further highlight the taxonomic relationships and the biodiversity of this peculiar tardigrade taxon.

SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article at the publisher's web-site.