Integrative description of a new Tunisian tardigrade species, Macrobiotus azzunae sp. nov. (Eutardigrada, Macrobiotidae, hufelandi group)

In this paper a new tardigrade species, Macrobiotus azzunae sp. nov., from Tunisia, is described. An integrative taxonomic approach was applied by combining morphological, morphometric and molecular data. In particular, light and scanning electron microscopy observations, and four genetic markers, three nuclear (18S rRNA, 28S rRNA and ITS-2) and one mitochondrial (COI) were used. The analysis showed that M. azzunae sp. nov. belongs to the Macrobiotus hufelandi group and is most similar to Macrobiotus sandrae Bertolani & Rebecchi, 1993. It differs from M. sandrae by a more pronounced constriction of the first macroplacoid (hardly visible in M. sandrae) and for the eggshell shape, with thinner wires of the reticulum and meshes around the processes larger than the inter-process meshes in M. azzunae sp. nov., while all meshes are similar in size in M. sandrae. The species is gonochoristic. With this discovery, there are 33 species of tardigrades identified in Tunisia, all nonmarine. This result, compared with nearby Sicily, where more research has been conducted, indicates that there is a considerable potential for identification of new species. Further research will be most informative if multiple habitats are explored and if carried out with an integrated approach as done in this present work.


Introduction
Tardigrades are hygrophilous micrometazoans whose outstanding resistance enables most of them to inhabit a large variety of habitats from the greatest depths of oceans to the highest mountain peaks, as well as extreme environments such as cryoconite holes. In terrestrial environments they colonize mosses, lichens, leaf litter and soil which receive stochastic hydration (Nelson et al. 2015;Schill 2018). Much research on the diversity and distribution of tardigrade fauna has been carried out in recent years in various parts of the world, especially considering the terrestrial environment. This has led to a significant increase in the number of known species: 531 in 1983 (Ramazzotti & Maucci 1983), about 960 in 2005 ) and more than 1300 species in 2020 (Degma et al. 2020). Nevertheless, knowledge of the diversity and distribution of North African tardigrades is still very limited and in particular in Tunisia, for which only three studies have been conducted. The first report on Tunisian tardigrades comes from a survey by Iharos (1978) in Northern Tunisian regions in which eleven species of tardigrades were detected. Then, Binda & Pilato (1987) studied the tardigrade fauna of Salambo (Tunis), Tabarka and Ain Drahem (Jendouba), and recently Gąsiorek et al. (2017) studied the tardigrade fauna of Bni Mtir, Jendouba. Currently, the known Tunisian tardigrade fauna includes 32 species, of which one, Bryodelphax maculatus Gąsiorek, Stec, Morek, Marnissi, Michalczyk, 2017, was originally discovered from the same region considered here, namely the Kroumirie Mountains. The Kroumirie Mountains are located in Northern Tunisia near the Algerian boundary and are characterized by a climate which varies from sub-humid to humid in winter. The Kroumirie forests are covered especially by cork oaks (Quercus suber L.) and Algerian oaks (Quercus canariensis Willd.) (Stambouli-Essassi et al. 2007). The analysis of the tardigrade fauna of a moss collected at the Kroumirie Mountains led us to find a new eutardigrade species (Macrobiotus azzunae sp. nov.) belonging to the Macrobiotus hufelandi group, considered as one of the most common groups of limnoterrestrial tardigrades on the planet (McInnes 1994;Kaczmarek & Michalczyk 2017a;McInnes et al. 2017). In this study we combined classical morphological and morphometric methods with modern molecular techniques in an integrative approach as suggested by Cesari et al. (2009Cesari et al. ( , 2011Cesari et al. ( , 2020, Stec et al. (2018b) and Kayastha et al (2020). Using phase contrast (PhC) light microscopy (LM), differential interference contrast (DIC) and scanning electron microscopy (SEM), we were able to describe the phenotypic characteristics of the new species whereas the amplification of DNA markers (three nuclear, 18S rRNA, 28S rRNA and ITS-2, and one mitochondrial, COI) provided barcodes for the genetic identification of this new species.
Measurements of animal length, claws and buccal-pharyngeal apparatus details were taken according to Pilato (1981). The external buccal tube diameter was measured at the level of the stylet support insertion point. We calculated the pt index which is the percentage ratio between the length of a structure and the length of the buccal tube measured from the medio-dorsal ridge of the buccal armature to the base of the pharyngeal apophyses (Pilato 1981). Single measurements are available in Supp. file 1.
Fifteen additional specimens were prepared for SEM analyses according to Bertolani et al. (2014), and observed with a Philips SEM XL 40, available at the Centro Interdipartimentale Grandi Strumenti (CIGS) of UNIMORE.

Molecular analysis
Before molecular analysis, each specimen was identified and photographed in vivo with LM using the method described by Cesari et al. (2011) in order to obtain voucher specimens. Supp. file 2 shows the buccal-pharyngeal apparatus of a hologenophore voucher specimen, showing how it is possible to obtain good morphological information on alive animals even when the entire animal has to be used for molecular analysis. Table 1 lists the specimens utilized. Genomic DNA was extracted from six Tunisian single animals using a rapid salt and ethanol precipitation following the protocol described by Cesari et al. (2009). We amplified four DNA fragments: the small ribosome subunit (18S rRNA), the large ribosome subunit (28S rRNA), the internal transcribed spacer (ITS-2), and the cytochrome oxidase subunit I (COI), using the primers and protocols described by Bertolani et al. (2014), Cesari et al. (2016), Stec et al. (2018d) and Cesari et al. (2009), respectively. Additionally, a portion of the 18S gene was sequenced for two additional specimens of M. sandrae sampled in the locus typicus of the species (Black Forest, Germany). GenBank accession numbers of obtained sequences are listed in Table 1. The amplified products were gel purified using the Wizard Gel and PCR cleaning (Promega) kit, and fragments were sequenced according to the protocols described in Bertolani et al. (2014). All COI sequences were translated into protein sequences in MEGAX (Kumar et al. 2018) to check for the presence of stop codons, and therefore for the presence of pseudogenes. Sequences of the four genes pertaining to species of the hufelandi group were used for molecular comparisons (Table 2) and aligned with the newly produced sequences using the MAFFT algorithm (Katoh et al. 2002) as implemented in the MAFFT online service (Katoh et al. 2017) and checked by visual inspection (Supp. file 3,Supp. file 4,Supp. file 5, Supp. file 6). Uncorrected pairwise distances were computed using MEGAX (Supp. file 7,Supp. file 8,Supp. file 9,Supp. file 10). Furthermore, relationships between COI haplotypes pertaining to the clade A nested in the M. hufelandi complex sensu Stec et al. (2021) were estimated using a parsimony network, by applying the method described in Templeton et al. (1992), as implemented in TCS ver. 1.21 (Clement et al. 2000) and visualized using tcsBU (Múrias dos Santos et al. 2016). A 95% connection limit was employed, as it has been suggested 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 coalescent-based species delimitation method that uses non-ultrameric 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 maximum likelihood (ML) tree computed using RAxML ver. 7.2.4 (Stamatakis 2006), as implemented in CIPRES (Miller et al. 2010), under the GTR+I+G model, as inferred by using the Akaike Information Criterion on jModelTest2 (Guindon & Gascuel 2003;Darriba et al. 2012). A sequence of Mesobiotus hilariae Vecchi, Cesari, Bertolani, Jönsson, Rebecchi & Guidetti, 2016 (GenBank accession number: KT226108) was used as outgroup. Bootstrap Table 1. List of specimens of M. azzunae sp. nov. and M. sandrae Bertolani & Rebecchi, 1993 utilized for molecular analyses. Abbreviation: N/A = not available. * These specimens were analyzed for the COI gene in Bertolani et al. (2011b).
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. In order to consider different evolutionary models for the three COI codons, a Bayesian tree was computed using the following models, as inferred by MrModeltest ver. 2 (Nylander 2004): SYM+I+G for the first position of the codon, GTR+I+G for the second position of the codon and GTR+G for the third position of the codon. The Bayesian dendrogram was computed with the program MrBayes ver. 3.2.7a (Huelsenbeck & Ronquist 2001;Ronquist & Huelsenbeck 2003), as implemented in CIPRES. Two independent runs, each of four Metropoliscoupled Markov chains Monte Carlo method, were launched for 2 × 10 7 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 ver. 1.7 (Rambaut et al. 2018) and the first 2 × 10 6 sampled generations were discarded as burn-in. In the distance-based 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 (Pmin) = 0.001; prior maximum divergence of intraspecific diversity (Pmax) = 0.1; Steps = 10 and gap width = 1.5. The analysis was performed on the ABGD website (https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html).

Adult specimens
Body white, transparent after mounting in Faure-Berlese, from 162.2 to 410.3 µm in length ( Fig. 2A, Table 3; structures measured only in the animals more than 200 µm in length). Eye spots present, even after mounting. Cuticle smooth but with small round or oval pores, 1-1.5 µm in diameter (Fig. 2B), better visible after fixation in Carnoy and orcein staining (Fig. 3C), scattered randomly on the entire cuticle, including the dorsal surface of all legs. With SEM, pores look oval or in the shape of a seed (Fig. 3A, D) with the largest diameter of 0.7-0.8 µm. Weak cuticular granulation also present on the lateral surface of all legs and specially on legs IV (Fig. 2B, arrow). Only with SEM is it possible to define the shape of the granulation on the legs, which looks as a regular disposition of star-shaped protuberances (about 0.3 µm; Fig. 3F). Six buccal sensory lobes around the mouth, well recognizable with SEM. Mouth antero-ventral; buccal-pharyngeal apparatus of the Macrobiotus type (sensu Pilato & Binda 2010), with ventral lamina and ten small peribuccal lamellae (in the holotype, after mounting, separated from the mouth). Buccal armature, corresponding to oral cavity armature, OCA, according to Michalczyk & Kaczmarek (2003), without an anterior band of teeth visible, corresponding to the first band of teeth according to Michalczyk & Kaczmarek (2003), and to the anterior band of the buccal ring according to Guidetti et al. (2012); posterior band of teeth poorly visible, corresponding to second band of teeth, according to Michalczyk & Kaczmarek (2003), followed by three dorsal and three ventral crests, corresponding to third band of teeth according to Michalczyk & Kaczmarek (2003); the dorsal crests (Fig. 2D) are distinct transverse ridges, whereas the ventral crests (Fig. 2E) appear as two separate lateral transverse ridges and a roundish median tooth. The posterior band of teeth and the transverse ridges are part of the buccal tube, according to Guidetti et al. (2012). Buccal tube narrow; pharyngeal bulb spherical with triangular apophyses, two rod-shaped macroplacoids, relatively short, the first longer than the second and evidently but not deeply narrowed at its middle (Fig. 2C), the second with a not particularly evident subterminal constriction. Microplacoid present. Slender claws of the hufelandi type (sensu Pilato & Binda 2010); the external claw longer than the internal one and the posterior longer than  The population is dioecious (gonochoristic). Males were recognized using orcein staining, which revealed that the testis is filled with spermatozoa with a coiled head (Fig. 3G) and spermatids. No morphological secondary sexual dimorphism, such as gibbosities on legs IV in males, was identified.

Eggs
Eggs are laid freely, and are white, spherical or slightly oval. One egg containing a fully developed embryo showed the shape of the buccal-pharyngeal apparatus (Fig. 4A). Processes of the eggshell are in the shape of inverted goblets (Fig. 4B) with conical trunks and well-defined distal discs as large as the process bases (for measurements see Table 4). Distal discs concave, with a median small protuberance and, using PhC, with border often smooth, or sometimes slightly jagged, or slightly ragged (Fig. 4C), Table 3. Measurements (in μm) and pt of selected morphological structures of individuals of Macrobiotus azzunae sp. nov. mounted in Faure-Berlese. Abbreviations: Nr = number of specimens/structures measured; Range = smallest-largest value; SD = standard deviation. * Three animals below 200 μm in length were not considered in the analysis. Two were not measurable, while data of the only measurable animal below 200 μm is available in the supplementary material (Supp. file 1).   but never clearly jagged, serrated or dentate. Surface among processes of the hufelandi type (sensu Kaczmarek & Michalczyk 2017a), i.e., covered by a very thin grid (Fig. 4D). Meshes around the process bases slightly larger and with slightly thicker wires compared with interbasal meshes. Mesh diameter around 0.5 µm.

Comparisons
Macrobiotus azzunae sp. nov. has eggs with processes as inverted goblets and a reticulate eggshell between the processes. Consequently, a comparison must be done with the Macrobiotus species listed by Kaczmarek & Michalczyk (2017a) with hufelandi type eggshells, excluding the species with processes that are not like inverted goblets, and adding the species with hufelandi type chorion eggs described after that publication. For the shape of the egg Macrobiotus azzunae sp. nov. differs from M. rawsoni Horning, Schuster & Grigarick, 1978 because this species has only one strip of meshes around each egg process (see Kaczmarek & Michalczyk 2017b); it differs from M. serratus Bertolani, Guidi & Rebecchi, 1996 because in this species the egg surface is porous more than reticulated, with pores small and spaced from each other, and its egg processes have a large, often square, distal disc; it differs from M. seychellensis Biserov, 1994 because the distal disc of the egg processes of this species, even though not dentate, has long and very developed lobes.
The remaining nine species of the hufelandi group should be compared singularly.

Macrobiotus almadai Fontoura, Pilato & Lisi, 2008
Macrobiotus azzunae sp. nov. differs from M. almadai in having a posterior band of teeth in the buccal cavity visible with LM (not visible in M. almadai), and distal disc with a jagged margin instead of very small teeth as in M. almadai.

Macrobiotus canaricus Stec, Krzywański & Michalczyk, 2018
With LM the margin of the distal disc of M. azzunae sp. nov., never dentate in this species, looks similar to that of M. canaricus, but the SEM images of the eggs of the latter species evidence the presence of an almost dentate disc. Moreover, the peribasal meshes of the eggshell are larger than interbasal ones in the new species while they do not differ from the interbasal ones in M. canaricus; regarding the animals there are differences in the buccal armature: in M. azzunae sp. nov. the posterior band of teeth is visible with LM (even though poorly) and the three dorsal crests are distinct transverse ridges, while in M. canaricus the posterior band of teeth is visible only with SEM and with LM the dorsal teeth form a transversal ridge weakly divided into three teeth.

Macrobiotus madegassus Maucci, 1993
The new species differs from M. madegassus by the presence of the eye spots (absent in M. madegassus), pores on the cuticle (absent in M. madegassus), presence in the buccal armature of posterior band of teeth, even though weak (fully absent in M. madegassus), buccal tube much larger (pt of the holotypes 15.9 vs 7), insertion of the stylet supports on the buccal tube much more posterior (pt of the holotypes 76.1 vs 68), first and second macroplacoid longer (pt of the holotypes 25.5 and 18.1 vs 21.3 and 12.0), lunules on the hind legs without kerning (crenate in M. madegassus), eggshell processes with distal disc as large as the base (similar range 3.2-5.2 for both measurements) with respect to that of M. madegassus (disc vs base: 4.3-5.4 vs 2.3-2.6).

Macrobiotus martini Bartels, Pilato, Lisi & Nelson, 2009
The cuticular pores in M. azzunae sp. nov. are much smaller than those of M. martini; the distal disc of the egg processes in M. azzunae sp. nov. has a diameter similar to that of the process base, while in M. martini the distal disc is much narrower than the base.

Macrobiotus nebrodensis Pilato, Sabella, D'Urso & Lisi, 2017
Macrobiotus azzunae sp. nov. differs from M. nebrodensis by the absence of the cuticular bar near the lunules on the first three pairs of legs (a faint bar is present in M. nebrodensis). The egg processes of M. azzuane sp. nov. are in higher number on the circumference (29-33) with respect to those of M. nebrodensis (18). In the latter species there are some egg processes very high (up to 20.6 µm), while in the new species process height and shape are more uniform. The difference in the eggshell between meshes around the process base and the others is much less evident in M. azzunae sp. nov. than in M. nebrodensis. Biserov, 1990 The new species differs from M. personatus by the posterior band of the buccal armature less evident, the presence of a clear constriction in the first macroplacoid (Fig. 5A), in the paratype of M. personatus examined by us barely identifiable (Fig. 6A) and, according to Biserov (1990) usually absent in the type material of that species. The pores on the cuticle of M. azzunae sp. nov. are small, approximately 1 µm in diameter, while in M. personatus they are strongly elliptic and about 3 µm in length (Fig. 6B). Lunules on leg IV are always smooth in M. azzunae sp. nov., sometimes indented in M personatus. With respect to the eggs of M. personatus (Fig. 6C-D), the egg processes of M. azzunae sp. nov. (Figs 4C-D, 5C) are clearly shorter, 5.4 ± 0.6 vs 9.5 ± 0.5 (range 4.2-6.4 vs 9-10.5) and with a narrower base and distal disc (both 3.2-5.2 vs 7-10.5 and 7-9 respectively). Males are present in the new species, while in M. personatus only females were found (Biserov 1990), suggesting parthenogenesis in that species.

Macrobiotus sandrae Bertolani & Rebecchi, 1993
The new species differs from M. sandrae for the eggshell shape, with thinner wires of the reticulum and meshes around the processes larger than the inter-process meshes in M. azzunae sp. nov. (Fig. 5C), all meshes similar in size in M. sandrae (Fig. 5D). Figure 5C-D also show a difference in the process base diameter, narrower in M. azzunae sp. nov.With regard to the animals, M. azzunae sp. nov. differs from M. sandrae by a constriction of the first macroplacoid much more pronounced ( Fig. 5A; it is hardly visible in M. sandrae; Fig. 5B). Moreover, in animals of similar size the posterior band of the buccal armature is just less evident in the new species, and lunules on the hind legs are without hint of teeth (but teeth, present in the holotype of M. sandrae, are often difficult to identify in other specimens of that species).

Macrobiotus terminalis Bertolani & Rebecchi, 1993
Macrobiotus azzunae sp. nov. differs from M. terminalis for the absence of granulation on the cuticle (noted only in the redescription of M. terminalis; see Cesari et al. 2011), for the absence of teeth on the lunules, especially evident on the hind legs of M. terminalis, and for the presence of males, absent in M. terminalis (see redescription by Cesari et al. 2011).

Macrobiotus vladimiri Bertolani, Biserov, Rebecchi & Cesari, 2011
With respect to M. vladimiri, animals of M. azzunae sp. nov. reach a shorter length (up to 410.3 µm vs 515.1 µm), in M. azzunae sp. nov. the posterior band of teeth of the buccal armature is less evident and the lunules on the hind legs are not indented. In M. azzunae sp. nov. the egg diameter without processes (64.7-80.6 µm) is less than that of the eggs of M. vladimiri (89.9-92.0 µm); the processes are shorter (4.2-6.4 µm in the new species vs 6.5-8 µm in M. vladimiri). In the new species the base process diameter is narrower (3.2-5.2 µm) than in M. vladimiri (5.1-7.3 µm), the distal disc is weakly or not jagged (clearly jagged in M. terminalis). In M. azzunae sp. nov. males are present, while they are absent in M. vladimiri.

Genetic distances
The ranges of uncorrected genetic p-distances between M. azzunae sp. nov. and the other species of the M. hufelandi group (Supp. file 7, Supp. file 8, Supp. file 9, Supp. file 10), are as follows:

18S
0.1-5.6%, with the most similar being M. sandrae from Germany (present paper) 28S 0.1%, with the only available M. vladimiri from Spain (FJ435751-5) ITS-2 7.7-32.2%, with the most similar being Macrobiotus vladimiri (MN888347) from Finland COI 6.3-25.6%, with the most similar being Macrobiotus sandrae (HQ876574, HQ876577, HQ876578, HQ876579, HQ876581) from Germany The COI dataset is the most complete and informative for species delimitation investigation. Both phylogenetic reconstructions on the COI dataset resulted in the same topology, and thus the ML tree was utilized for the PTP analysis (Fig. 7, (Fig. 7, centre and right). Present molecular data therefore confirms the validity of the erection of M. azzunae sp. nov.

Discussion
With the discovery of M. azzunae sp. nov. the species identified in Tunisia now comes to 33, all nonmarine. This number is much lower than the number of non-marine species, 120, found in nearby Sicily (Pilato et al. 2017(Pilato et al. , 2019, an island extensively studied from a tardigradological point of view. This means that there is considerable potential for further discoveries that could come from the study of the various habitats colonized by tardigrades in Tunisia. Furthermore, the results will undoubtedly be greater and more informative if the research is carried out with an integrated approach, as done in this work. Results of the Poisson tree process analysis are provided using differently coloured branches: putative species are indicated using transitions from 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: haplotype network of COI gene in M. hufelandi complex. Circles represent haplotypes, while circle surface denotes haplotype frequency. Networks falling below the value of the 95% connection limit are disconnected. Right: rectangles denote specimens grouped by ABGD analysis.