Novel integrative data for two Milnesium Doyère, 1840 (Tardigrada: Apochela) species from Central Asia


 Tardigrada
 are a phylum of microscopic animals inhabiting a variety of ecosystems, both aquatic and terrestrial, being recognised for their remarkable abilities to withstand tough environmental conditions. The order Apochela groups exclusively carnivorous species, with the vast majority representing the genus Milnesium Doyère, 1840. Representatives of this genus are characterised by simplified morphology, therefore possessing an extremely limited set of taxonomically meaningful morphological traits. Nevertheless, the taxonomy of Milnesium is mostly based on classical data: observations and measurements in light microscopy with the majority of descriptions lacking integrative data, most importantly DNA barcodes, but also scanning electron microscopy photographs and developmental variability analysis. Hence, re-descriptions that include novel integrative data are urgently needed. In this contribution, we provide new taxonomic data for two species described from Central Asia, Milnesium almatyense (a single population) and Milnesium reductum Tumanov, 2006 (five populations): morphometrics, DNA barcodes, SEM observations and description of developmental variability. As a result, we amend the description of both species and reveal phylogenetic relationships of those species and other sequenced congeners. The integrative data confirm the validity of the two species and include them in the growing set of Milnesium species associated with DNA sequences.


Introduction
Tardigrades, commonly named as water bears, are a phylum of cosmopolitan invertebrates, inhabiting almost all environments across the world (Nelson et al. 2018). Currently around 1300 tardigrade species are recognised (Degma et al. 2009(Degma et al. -2019 with 44 belonging to the order Apochela Schuster et al., 1980 which is characterised by a set of unique traits, i.e. a separation of the primary and secondary claw branches or the presence of peribuccal lamellae. Additionally, in contrast to parachelans, apochelans lack important morphological characters commonly utilised in their taxonomy, for example, placoids and septulum in the muscle pharynx or ornamented egg shells. The great majority of extant species (41 of 44) within the family Milnesiidae Ramazzotti, 1962, the only one in order Apochela, are classified within genus Milnesium Doyère, 1840(37 listed in Degma et al. 2009, excluding two nomina dubia: M. dujiangensis and M. tardigradum trispinosa, but with a further four taxa described later on by Kaczmarek et al. 2019, Moreno-Talamantes et al. 2019, and Surmacz et al. 2019. Due to the limited number of taxonomically meaningful traits the classification of this genus is challenging and, so far, has been based mostly on classical morphological data that were shown to be at least partially insufficient in resolving the taxonomy of these tardigrades, as pseudocryptic species have recently started to be recognised (e.g. Morek et al. 2016a;Morek et al. 2019a;Surmacz et al. 2020). Moreover, it has been demonstrated that some Milnesium species exhibit developmental variability, in which the immature life stages, i.e. hatchlings and juveniles (1 st and 2 nd instar, respectively; Morek et al. 2019b) may exhibit different morphologies than the adults (3 rd instar onwards), for example, in the number of points on secondary branches of claws (so called claw configuration, CC; Michalczyk et al. 2012a, b), dorsal cuticle sculpturing or buccal tube morphology (Morek et al. 2016a;Morek et al. 2019b;Morek et al. in press).
So far, the tardigrade fauna of Kazakhstan and Kyrgyz Republic have not been thoroughly investigated, with only a few contributions published in the current century, focusing mostly on descriptions of new species (e.g.: Tumanov 2003Tumanov , 2005Tumanov , 2006Tumanov , 2007Kaczmarek et al. 2018;Zawierucha et al. 2018or Coughlan et al. 2019. In this region, three Milnesium species were described and two of them, M. almatyense Tumanov, 2006 andM. reductum Tumanov, 2006, are characterised by rare morphological traits: the [2-3]-[2-2] CC (M. almatyense) and the lack of accessory points (M. reductum). Although the descriptions are accurate and detailed, they lack the integrative data, such as scanning electron microscope imaging, information on ontogenetic variability and, most importantly, DNA sequences. Additional phenotypic data and genetic data for these species are needed, not only to aid their identification and to test against potential synonymies, but also to pinpoint the phylogenetic positions of both taxa. The phyletic relationships of M. almatyense and M. reductum are particularly interesting because of the rare phenotypic traits mentioned above. Specifically, in the first extensive phylogeny of Milnesium in Morek and Michalczyk (2020) only two species with a [2-3]-[2-2] CC and no species without accessory points were included.
Thus, the aim of this study was to collect integrative data for two species, M. almatyense and M. reductum, based on the material collected in the proximity of their loci typici in Kazakhstan and Kyrgyz Republic, respectively. The integrative analysis of six populations of the two species provided novel morphological and developmental traits and DNA barcodes, amending the descriptions of these taxa and allowed us to identify their kin.

Sampling and specimens
We analysed six samples containing M. almatyense (one sample) and M. reductum (five samples) originating from Kazakhstan and Kyrgyz Republic (detailed collection data are listed in Table 1). The samples were examined according to the protocol described by Stec et al. (2015). The extracted specimens were afterwards split into four analyses: (i) imaging and morphometry in phase-contrast light microscopy (PCM), (ii) DNA extraction and sequencing, (iii) imaging in scanning electron microscopy (SEM) and (iv) developmental variability analysis (Morek et al. 2016a). The exact numbers of specimens per population utilised for each analysis are provided in Table 1.
In addition to the new populations, paratypes of both species were examined: a single specimen of M. al-matyense (slide No. 199) and five specimens of M. reductum (slide No. 192). The slides were loaned from the Zoological Institute of the Russian Academy of Sciences, St. Petersburg, Russia.

Microscopy, imaging and morphometry
The specimens were mounted on permanent microscope slides according to the method by Morek et al. (2016b). The measurements of buccal tube follow Michalczyk et al. (2012a, b), whereas the remaining traits were measured following Tumanov (2006). The branch height ratio is a ratio of the secondary to the primary branch height and it is expressed as a percentage. The pt is a ratio of a given structure to the length of the buccal tube (Pilato 1981), expressed as a percentage and, in the text and tables, is given in italics. The number of measured specimens follows the recommendations by Stec et al. (2016). The morphometric data were handled using the Apochela spreadsheet ver. 1.3., available from Tardigrada Register, www.tardigrada.net (Michalczyk and Kaczmarek 2013). The measurements and photographs were taken with an Olympus BX53 PCM associated with an Olympus DP74 digital camera. For deep focus structures that could not be focused in a single photo, a series of up to 24 photos were taken and then merged into one focused image using Corel Photo-Paint X8. Specimens were processed for SEM imaging according to the protocol by Stec et al. (2015) and examined under high vacuum with Versa 3D DualBeam SEM at the ATOMIN facility of the Jagiellonian University.

Developmental variability analysis
Cultures were established from live specimens and viable eggs deposited in exuviae, which were incubated at standard rearing conditions described by  with rotifers Lecane inermis Bryce, 1892, as a food source. In order to test for the potential presence of developmental variability, hatching was applied because developmental tracking (Morek et al. 2016a) failed (hatchlings died before reaching the juvenile stage). As ontogenetic variability in CC in both analysed species was detected, we further applied the analytical method described by Surmacz et al. (2020) to assess whether our morphometric dataset encompassed juveniles in addition to hatchlings and adults and to determine whether the species exhibit early or late CC change.
Moreover, the analytical method by Surmacz et al. (2020) was also applied to data available for M. berladnicorum Ciobanu et al., 2014 from the Tardigrada Register and an additional, topotypic hatchling (exhibiting the [2-2]-[2-2] CC), to test whether the species exhibits the early or late change. The data unequivocally indicated the early positive shift to [2-3]-[2-2] CC (see also Suppl. material 1).

Genotyping
Genomic DNA was extracted from four individuals from each of the six analysed populations (except for a single population of M. reductum where only one specimen was sequenced). The extraction method follows Chelex ® 100 resin (Bio-Rad) protocol by Casquet et al. (2012), with modifications by Stec et al. (2015). Prior to extraction, the specimens were mounted on temporary water slides to determine their CC. Whenever possible, voucher specimens were obtained from the extraction vial (in total 13 vouchers out of 21 extractions). Four genetic markers were sequenced, three nuclear: the small ribosomal subunit (18S rRNA), large ribosomal subunit (28S rRNA), Internal Transcribed Spacer 2 (ITS-2); and one mitochondrial, Cytochrome Oxidase C subunit I (COI). The PCR protocols follow Stec et al. (2015), primers and PCR programmes with relevant references are listed in Table  2. The chromatograms were manually checked using BioEdit ver. 7.2.5 (Hall 1999). As the COI is a protein coding gene, the obtained fragments were translated into amino acids using MEGA 7 (Kumar et al. 2016) to test against pseudogenes. All sequences are deposited in GenBank (the accession numbers are listed in Table 1).

Phylogenetic analysis
To uncover the phylogenetic relationships of the two analysed species, the recent dataset by Morek and Michalczyk (2020) Stec et al. 2018a). In next step, the obtained alignments were visually checked in BioEdit and trimmed to 1047 (18S rRNA), 865 (28S rRNA), 648 bp (ITS-2) and 580 bp (COI). Afterwards, the four obtained alignments were concatenated in SequenceMatrix (Vaidya et al. 2011). We used the Bayesian Information Criterion in PartitionFinder version 2.1.1 (Lanfear et al. 2016) in order to find the most suitable substitution model for posterior phylogenetic analysis. As the COI is a protein coding fragment, before partitioning, the alignment was always divided into three data blocks constituting three separated codon positions. The analysis was run to test for all possible models implemented in the programme (for further, Bayesian Inference, BI). The best fit-models for six partitions in BI were: GTR+I+G for 18S rRNA, 28S rRNA, the second and the third codon positions in COI, whereas the GTR+G was indicated for ITS-2 and the first codon position in COI. BI marginal posterior probabilities were calculated using MrBayes v3.2 (Ronquist and Huelsenbeck 2003). Ran- dom starting trees were utilised and the analysis was run for ten million generations, sampling the Markov chain every 1000 generations. An average standard deviation of split frequencies of < 0.01 was used as an indicator that the two independent analyses had converged. We used Tracer v1.3 (Rambaut et al. 2014) to ensure Markov chains had reached stationarity and to determine the correct "burn-in" for the analysis (the first 10% of generations). A consensus tree was obtained after summarising the resulting topologies and discarding the "burn-in". In the BI consensus tree, clades recovered with posterior probability (PP) between 0.95 and 1.00 were considered well supported, those with PP between 0.90 and 0.94 were considered moderately supported and those with lower PP were considered unsupported. The consensus tree was visualised in FigTree v.1.4.3, available from http://tree.bio.ed.ac.uk/software/figtree. As the dataset analysed herein is only slightly larger than the dataset by Morek and Michalczyk (2020), we did not perform the Maximum Likelihood analysis, as previously it resulted in weakly-supported nodes.

Genetic diversity of analysed populations
Milnesium almatyense was represented by single haplotype in each gene, whereas M. reductum exhibited three 18S rRNA haplotypes, four 28S rRNA haplotypes, three ITS-2 haplotypes and four haplotypes of COI. The maximum genetic distance between the analysed populations of M. reductum was as follows: 0.2% in 18S rRNA, 0.5% in 28S rRNA, 1.3% in ITS-2 and 2.9% in COI, thus there was very little structuring within the M. reductum clade. The summary of haplotypes and the matrices of genetic distances between the populations are available in Suppl. materials 2.

Novel morphological traits
Thanks to the observations in high quality PCM and SEM, as well as morphometric measurements of multiple specimens, including sexually-immature instars (the analysis with the algorithm developed by Surmacz et al. 2020, confirmed the presence of the first three life stag-es), we were able to identify novel traits and update differential diagnoses for both species.
Specifically, for M. almatyense ( Fig. 2A): 1. Morphometric data missing in the original description, such as the dimensions of peribuccal and lateral papillae, anterior and posterior widths of buccal tube, II internal, I-III external and anterior claws lengths, are provided in Table 3. Morphometric ranges have expanded in the majority of traits due to the increased sample size (six in the original description and 13 herein), but mostly because of the inclusion of immature life stages, which were not present in the original description. The only discrepancy in the measurements was the standard buccal tube width both in absolute and relative dimensions (14.1-16.3 µm and 36.5-44.0 in the original description vs. 6.5-13.1 µm and 26.2-33.1 herein). The gap in absolute dimensions can be explained by the corresponding gap in body lengths between the type series and the new population analysed herein: 672-798 µm vs. 311-665 µm, respectively (Tumanov 2006 measured only larger specimens, hence the skewed body size in the original description). The gap in relative values is harder to explain, but (i) the pt index does not remove the allometric effects entirely (Bartels et al. 2011) and (ii) buccal tube width is susceptible to deformation under cover-glass pressure and even minor changes in pressure may considerably affect pt values (Morek et al. 2016b). Importantly, the susceptibility to deformation increases with the size of the buccal tube (Morek et al. 2016b), thus, given that larger specimens were present in the type series than in the new population, any differences in cover-slip pressure would magnify the effect, resulting in a discrepancy in buccal tube width. Thus, we consider the gaps in buccal tube width between the type series and the new population analysed herein as an artefact caused by differences in body size, possibly magnified by differences in cover-slip pressure.

The analysis of ontogenetic variability in CC
showed that hatchlings (Fig. 1A) exhibit a [2-2]-[2-2] CC (Fig. 1B), whereas juveniles and adults have . 4. Observation of the mouth opening in SEM confirmed that the species has six peribuccal lamellae, but it also indicated that the lamellae are of unequal size, i.e. the pair of dorsal and ventral lamellae are larger than the two lateral lamellae, i.e. the configuration is 4+2 (Fig. 2C).

Observations under a high quality PCM and SEM
indicated that M. almatyense has a sculptured cuti-cle, which is visible both, in the population analysed herein (Fig. 3A, C, E, F), as well as in the analysed paratype (Fig. 3B, D). The sculpturing has a form of delicate, irregular wrinkles, similarly to M. berladnicorum and M. variefidum (Morek et al. 2016a). Pseudopores and dorsal pseudoplates are also clearly visible in the analysed paratype (Fig. 3B) and the new population (the scheme of the pseudoplate arrangement is depicted in Fig. 2B, based on the new population). Row I is not visible, whereas for rows II and III, the determination of exact shape was impossible, due to the obscured outline of the pseudoplates (numbering according to Moreno-Talamantes et al. 2019). The cuticle in hatchlings is similar, but the dorsal sculpture is more evident and regular, forming a delicate reticulation (Fig. 1C).
Similarly, for M. reductum (Fig. 4A): 1. Morphometric data missing in the original description, such as the dimensions of peribuccal and lateral papillae, anterior and posterior widths of buccal tube, II internal, I-III external and anterior claws  Table 4. Morphometric ranges have slightly expanded in the majority of traits due to the increased sample size (six in the original description and 34 herein) and mostly because of the inclusion of immature life stages, which were not recognised in the original description.

SEM observations confirmed that the species has
six peribuccal lamellae, but it also indicated that the lamellae are of unequal size, i.e. the pair of dorsal and ventral lamellae are larger than the two lateral lamellae, i.e. the configuration is 4+2 (Fig. 4B).

Observations under a high quality PCM of all new
populations, as well as of the analysed paratype, revealed that this species possesses pseudopores visible in the dorsal cuticle (Fig. 4D, E). A weak outline of a single pseudoplate is also visible (Fig.  4C, E). The SEM imaging further confirmed the species has s smooth dorsal cuticle, with a barely visible pseudoplate in row VIII (Fig. 4C). As only a single pseudoplate is observable both in PCM and SEM, a scheme of pseudoplate arrangement is not provided.  On the other hand, M. reductum is represented by five closely-related populations forming a well-supported subclade that is in a sister relationship with all remaining taxa of "clade A" in Morek and Michalczyk (2020).
In the currently available dataset, the species with the closest affinity to M. almatyense is M. berladnicorum, thus we compared the genetic distances between the two species, which are as follows: 0.2% in 18S rRNA, 0.9% in 28S rRNA, 4.0% in ITS-2 and 6.9% in COI. M. reductum is indicated as the sister species for the entire clade (genetic distances between M. reductum and the remaining species of clade A are as follows: 0.5-2.4% in 18S rRNA, 3.4-5.7% in 28S rRNA, 10.2-16.3% in ITS-2 and 13.8-16.3% in COI).
Considering the similarity and problematic taxonomy of M. tardigradum and M. pseudotardigradum, as well as of M. almatyense, M. berladnicorum and M. variefidum, we introduce names for these two species clusters: the tardigradum complex and the almatyense complex, respectively (Fig. 6). The names are derived from the first species described within respective complexes.   al. 2019b; Surmacz et al. 2019;Morek et al. in press;present contribution), which is a worryingly low percentage, considering the challenging taxonomy of this genus reflected by the low number of useful phenotypic traits for species differentiation and the widespread developmental variability. Moreover, within the 14 species for which integrative data are available, only in 10 species has the ontogenetic variability been analysed. This shows how much there is still to be done to stabilise Milnesium taxonomy, both in terms of re-descriptions and encouraging taxonomists to publish descriptions supported by ontogenetic variability and genetic analyses. Especially, that recent studies revealed a considerable undescribed species diversity within the genus  and demonstrated evidence for pseudocryptic Milnesium species (Morek et al. 2019b;Surmacz et al. 2019), the identification of which would not have been possible with the sole use of classical tools.

Discussion
Although the number of species that were tested against developmental variability in CC is not too high (20 species, including the undescribed species that were analysed by Morek and Michalczyk 2020), some patterns seem to emerge. Specifically, amongst the 10 tested species with adult [3-3]-[3-3] CC, ontogenetic variability was observed only in a single species (an undescribed Milnesium sp. nov. 10 PH.014 in Morek and Michalczyk 2020).
In contrast, all 10 tested species with adult CC different from [3-3]-[3-3], exhibited developmental variability in CC. Moreover, in the majority of cases (7 of 10 species; 70%) an increase in the number of secondary claw branches is observed (i.e. positive CC change). Importantly, in the vast majority of cases (9 of 10 species; 90%), the shift takes place between hatchlings and juveniles (i.e. early CC change). Given that hatchlings are most likely to be missing in the type series, taxonomists should be especially cautious when describing and analysing taxa in which adults exhibit the CC other than [3-3]-[3-3].  Figure 6. Fragment of the Bayesian phylogenetic tree ("clade A"), based on the analysis of concatenated 18S rRNA + 28S rRNA + ITS-2 + COI nucleotide sequences, obtained by Morek and Michalczyk (2020) with an addition of six new populations, showing the positions of the two Milnesium species analysed in this contribution (shaded in grey): M. almatyense Tumanov, 2006 andM. reductum Tumanov, 2006. The numbers at nodes represent Posterior Probability (PP) supports, with the values in grey font indicating poor support (conventionally recognised as polytomy). The scale bar shows the number of substitutions per site. For the remaining parts of the tree, please see Morek and Michalczyk (2020).
In contrast to the original description of M. almatyense, our analysis of the type and new material showed that the dorsal cuticle is not smooth, but covered with fine sculpturing and pseudoplates. These characteristics were missed due to the insufficient quality of the light microscope used by Tumanov (2006) and because such traits were not recognised at the time. This amendment, along with the exact same pattern of ontogenetic CC change (i.e. early positive) makes M. almatyense extremely similar to M. berladnicorum. Specifically, in the original differential diagnosis of M. berladnicorum (Ciobanu et al. 2014), the species was discriminated from M. almatyense by the sculptured cuticle and ontogenetic variability in CC was not known to occur in either species. Moreover, measurements of all morphometric traits presented herein (Table 3) and in the original description (Table 3 in Tumanov  However, the absence of eyes was observed by Tumanov (2006) in specimens fixed on microscope slides, thus it is not known whether the eyes were absent or they were present in living animals, but they dissolved in acetic acid and Faure medium (D. Tumanov, pers. inf.). Importantly, the eyes were present in all living individuals in the M. almatyense population analysed herein and dissolved in Hoyer's medium in hatchlings. Thus, it could be assumed that M. almatyense exhibits eyes as M. berladnicorum does, making this trait invalid in the differentiation of these two species. We have, however, found a single morphological trait that could potentially differentiate M. almatyense and M. berladnicorum: pseudoplate arrangement, specifically in M. almatyense pseudoplate IV is singular (Fig. 2B), but it is double in M. berladnicorum ( fig. 6A in Morek et al. 2016a). However, the taxonomic value of dorsal pseudoplates in Milnesium has been recently questioned by Moreno-Talamantes et al. (2019). Thus, further studies on this trait are needed to validate its value in species differentiation.
Moreover, the new genetic data, presented herein, show that M. almatyense and M. berladnicorum are a pair of a very closely-related species (Fig. 6), with genetic distances in the analysed variable markers between the two species being moderate in both ITS-2 (4.0%) and COI (6.9%). Considering that there are Milnesium species with high intraspecific genetic distances (e.g. M. tardigradum with p-distances up to 4.0% in ITS-2 and 11.4% in COI, Morek et al. 2019bor M. eurystomum Maucci, 1991 with respective maximum distances of 1.7% and 8.9%, Morek et al. in press), it is not clear whether M. almatyense and M. berladnicorum are a single or distinct species. Thus, a more detailed analysis, for example, with the use of new genetic markers, is needed to verify whether M. berladnicorum is a good species or a junior synonym of M. almatyense. If the taxonomic status of M. berladnicorum is supported, then this pair of species should be considered pseudocryptic, as the phenotypic differences are difficult to spot and are not straightforward. In other words, a confident identification of these species may require DNA barcoding.
The recent phylogenetic analysis of the genus Milnesium suggested that species generally cluster by geography . Thus, as both M. almatyense and M. reductum originate from the Palearctic realm (sensu Holt et al., 2013), we expected them to be placed either in clade A or clade B shown by Morek and Michalczyk (2020) and the current analysis confirmed this prediction by demonstrating that both species represent clade A (Fig. 6). Therefore, the correlation between molecular phylogeny and geographic origin of populations is maintained and strengthened. M. reductum is so far the only species characterised by the lack of accessory points, for which its phylogenetic position is known. Therefore, it is difficult to speculate, whether this trait was lost only once in Milnesium or this is yet another example of convergent evolution (as has been shown for many other traits, such as cuticular sculpturing, claw or lamellae configuration; Morek and Michalczyk 2020). However, given that generally Milnesium species seem to exhibit limited geographic ranges and since the other four known species that lack the accessory points originate from different zoogeographic realms, i.e. India (M. longiungue Tumanov, 2006) and USA (M. alabamae Wallendorf & Miller, 2009;M. swansoni Young, Chappell, Miller &Lowman, 2016 andM. zsalakoae Meyer &Hinton, 2010), it would not be surprising if this trait was lost independently in different lineages.

Supplementary material 1
The relationships be tween body length and buccal tube length in Milnesium berladnicorum Authors