Environmental DNA metabarcoding of Danish soil samples reveals new insight into the hidden diversity of eutardigrades in Denmark

Abstract

Tardigrades are rarely included in large biodiversity surveys, mainly because of the impracticalities that follow larger sampling and identification of these microscopic animals. Consequently, there is a lack of data on their biogeographical distribution. Here, we analyse environmental DNA sequences of eutardigrades obtained with a metabarcoding protocol on Danish soil samples collected during a national biodiversity project (Biowide). Specifically, we aimed to investigate the applicability of the V4 region (~400 bp) of the 18S rRNA marker gene to assign taxonomy to 96 eutardigrade molecular operational taxonomic units (MOTUs), using three different methods (alignment-, tree- and phylogeny-based methods). Tardigrade reference libraries are currently suffering from insufficient taxon coverage, in some cases challenging the interpretation of data based on similarity searches. This can, to some extent, be accounted for by supplementing identification with algorithms that incorporate a backbone phylogeny and infer models of evolution. Together, the present findings suggest that the V4 region of the 18S rRNA gene offers a promising tool to identify unknown MOTUs of eutardigrades to genus or family level and can, in some cases, be used to assign to species level.

Eutardigrades, metabarcoding, environmental DNA, 18S, biodiversity, taxonomic resolution

INTRODUCTION

Tardigrades are ubiquitous in nature and have been found in many biomes on Earth, from mountaintops (Dastych 2004) to the deep sea (Renaud-Mornant 1987, Hansen 2007) and from the cold poles (Møbjerg et al. 2007, Vecchi et al. 2016) to tropical forests (Stec et al. 2018). Nevertheless, little is known about the distribution of tardigrades over different spatial scales or the mechanisms behind their abundances in different environments (McInnes and Pugh 2007, Ramsay et al. 2021). Because of their microscopic size, tardigrades are hard to monitor using traditional biomonitoring methods, and thus they face the same challenges as many other invertebrates in biodiversity assessments (Cardoso et al. 2011): many species remain to be discovered; the global distributions of many described species are unknown; the abundances of tardigrades, including how they change spatially and temporally within natural settings, are understudied; and knowledge on their ecology and their tolerances to natural habitat changes remain unknown for most species. Hence, it is no surprise that the otherwise well-documented phylum is currently suffering from an incoherence of biogeographical data, with most of these data originating from sampling localities of voucher specimens. Extant tardigrades are represented by two major evolutionary lineages, Eutardigrada and Heterotardigrada (Jørgensen et al. 2018, Degma et al. 2022). The present study revolves around semi-terrestrial eutardigrade species presumably inhabiting soil samples.

The bioinformatic field is maturing rapidly, with a wealth of well-understood methods and algorithms becoming increasingly available. With the increased accessibility of simple genetic tools over the last decades, supplementing morphological data with barcoding and experimental data with genomics and transcriptomics has become increasingly popular in studies of tardigrade evolution, phylogeny, systematics, taxonomy, phylogeography (e.g. Sands et al. 2008, Bertolani et al. 2014, Jørgensen et al. 2018, Guil et al. 2019, Gąsiorek and Michalczyk 2020, Møbjerg et al. 2020, Morek and Michalczyk 2020, Stec et al. 2020, 2021a, Morek et al. 2021, Tumanov 2022, Kayastha et al. 2023, Vecchi et al. 2023), adaptation and physiology (e.g. Hashimoto et al. 2016, Kamilari et al. 2019, Arakawa 2022, Møbjerg et al. 2022, Neves et al. 2022, Yoshida et al. 2022). In many ways, DNA barcoding has proved to be a valuable addition to the taxonomic classification and is now considered a standard practice in species descriptions. However, DNA barcoding, using Sanger sequencing, can only add to the taxonomic aspect of establishing a list of species and does not reduce the sampling effort (Coissac et al. 2012). Formal species descriptions rely on phenotypic traits, but tardigrades are small (50–1200 µm) and largely invisible to the naked eye. Traditional methods involve hand-picking microscopic individuals from a sample and mounting these onto microscopy slides for species identification, which requires a high degree of taxonomic expertise.

Reliable identification of organisms is essential to explore their spatial and temporal occurrences and to study their diversification, ecology and dispersal. Organisms typically represented in biodiversity assessments include flowers, insects and larger animals (Thomsen and Willerslev 2015, Brunbjerg et al. 2019), with occurrence data from many organisms publicly available in online databases, such as GBIF (Global Biodiversity Information Facility). Museum collections of tardigrades were recently included in this database (GBIF.org 2022); however, the data are still in their infancy, primarily showing collection sites of voucher specimens. The need for a comprehensive library of tardigrades is recognized, and over the past decade efforts have been made to increase the taxonomic coverage of genomic databases, such as GenBank (Sands et al. 2008, Jørgensen et al. 2010, Guil and Giribet 2012, Bertolani et al. 2014), laying the groundwork for metabarcoding to provide interpretable results.

DNA metabarcoding is a combination of DNA barcoding with high-throughput sequencing (Deiner et al. 2017), which allows simultaneous production of millions of DNA sequences, as opposed to the Sanger sequencing method of manually processing one sequence at a time. Living organisms leave traces of their DNA in the environment in which they live (Willerslev et al. 2003), commonly known as environmental DNA (eDNA). The eDNA can be extracted directly from environmental samples (e.g. sediment, water, air or debris samples), amplified using universal primers and sequenced using next-generation sequencing to generate millions of reads, usually grouped into molecular operational taxonomic units (MOTUs) (Ruppert et al. 2019). Environmental DNA metabarcoding has wide applicability and presents some practical advantages, being cost-effective and non-invasive and allowing for rapid analysis of large samples, with the opportunity to cover a wide range of taxa (Deiner et al. 2017). Accordingly, eDNA metabarcoding has revolutionized the exploration of biodiversity and has already been used to describe the diversity of meiofauna (de Faria et al. 2018), macrofauna (Lynggaard et al. 2022), fungi (Frøslev et al. 2019) and plants (Barnes et al. 2022). It can also be applied in understudied habitats, which are hard to access (Oliverio et al. 2018). Identifications are based on comparisons of obtained environmental sequences Molecular Operational Taxonomic Units (MOTUs) with a taxonomically annotated collection of DNA sequences, such as GenBank’s nucleotide database (Benson et al. 2013). Methods for taxonomic annotation can generally be grouped into three categories: alignment-, tree- and phylogeny-based methods (Holovachov et al. 2017).

Alignment-based approaches find regions of local similarity, based solely on the alignment between a MOTU and reference sequences, finding the reference sequence with the highest nucleotide similarity to the specific MOTU (Holovachov 2016a). Alignment-based approaches can be performed directly through the BLASTN function of the NCBI server (Sayers et al. 2020) or other search algorithms, such as VSEARCH (Rognes et al. 2016).

Tree-based approaches use a phylogenetic tree composed of reference sequences (a backbone tree) to construct a phylogeny de novo. The method incorporates the MOTUs into the backbone tree, using phylogeny inference algorithms and models of evolution (Holovachov 2016b). To assign taxonomy, the position of MOTUs within the cladogram is evaluated with respect to the neighbouring taxa and the branch support. If a reliable backbone tree is already known for a subset of reference sequences, it is computationally more efficient to find the position of the MOTUs within this existing tree (Balaban et al. 2020, Smirnov and Warnow 2021).

Phylogeny-based approaches also use a backbone tree; however, instead of incorporating the MOTUs into the tree, the tree topology is fixed, and MOTUs are placed onto existing branches within the tree, using phylogenetic placement algorithms (Czech et al. 2022). The placement of each MOTU is tested across all nodes in the topology, and the algorithm finds the most likely insertion position for every MOTU, based on either distance or maximum likelihood (ML).

Tardigrade diversity in Denmark remains to be described properly, and diversity in Danish soil has not been covered since the study by Hallas and Yeates (1972). In the present study, we analysed data from the metabarcoding eDNA project Biowide, covered in the study by Brunbjerg et al. (2019), which offers a unique opportunity to explore the applicability of metabarcoding (using the V4 region of the 18S rRNA gene) to identify eutardigrade species associated with soil samples collected across Denmark. Given that species and unique taxa constitute a fundamental unit of metabarcoding diversity assessments, the most important step in this process is the taxonomic assignment of MOTUs. Specifically, we evaluated the extent to which MOTUs of eutardigrades associated with Danish soil samples (Brunbjerg et al. 2019) can be assigned to the family, genus or species level using three different assignment approaches (alignment-based, tree-based and phylogeny-based).

MATERIALS AND METHODS

Data collection

The 18S rDNA sequences analysed in the present study originate from the biodiversity project Biowide, which combined traditional species identification with eDNA metabarcoding, covering 130 terrestrial sampling sites in Denmark (Brunbjerg et al. 2019). Specifically, in order to create an eDNA inventory, 81 soil cores were taken from each of these 130 sites at 15 cm depth (see Brunbjerg et al. 2019: appendix B). Large debris and visible plant parts were removed from the cores before they were pooled and homogenized. DNA was subsequently extracted for marker gene amplification and sequencing. For the 18S rRNA marker gene, the V4 region was amplified using universal eukaryote primers [TAReuk454FWD1 (5ʹ-CCAGCASCYGCGGTAATTCC-3ʹ) and TAReukREV3 (5ʹ-ACTTTCGTTCTTGATYRA-3ʹ)] according to Stoeck et al. (2010), and amplicons subsequently underwent 300 bp paired-end sequencing (Illumina MiSeq platform). The 18S metabarcoding sequence data contained a total of 45 468 MOTUs, and the data have previously been used in other publications (Brunbjerg et al. 2019, Frøslev et al. 2022a, 2022b). For details on methods for soil sampling, DNA extraction and sequencing, see Frøslev et al. (2017) and Brunbjerg et al. (2019). For this study, all tardigrade MOTUs retrieved were of eutardigrade origin. Notably, heterotardigrades seem to have large insertions within the 18S rRNA V4 region, making the barcoded fragment too long for MiSeq Illumina sequencing. Hence, the primers from Stoeck et al. (2010) cannot be used for metabarcoding of heterotardigrades.

Data cleaning and filtering

We constructed an 18S reference database composed of 368 sequences representing the phylogenetic diversity of eutardigrades available on GenBank. Retrieved sequences were trimmed to the barcoding region. To extract potential tardigrade sequences from the eDNA 18S eukaryote dataset, we matched all environmental sequences against this reference database using VSEARCH. Specifically, we retained sequences that matched any reference sequence with ≥ 87% sequence identity using dada2 (Callahan et al. 2016). The retrieved sequences (provided in Supporting Information, File S1) were aligned initially with MAFFT v.7.490, using default parameters (Katoh et al. 2009), and inspected in MEGA7 v.7.0.26 (Kumar et al. 2016). Ambiguous sequences were subsequently blasted using BLASTN on NCBI. Forty-two sequences were matched with organisms not belonging to Tardigrada and thus removed from further analyses. The remaining sequences, all of eutardigrade origin, were subsequently filtered, applying a lower threshold of 10 reads, i.e. only MOTUs represented by ≥ 10 reads in a soil sample were included in the analyses. This filtering step was applied because damaged eDNA can acquire the appearance of a new genotype (Deiner et al. 2017). MOTUs with potential errors are likely to have a low number of reads, hence the removal of low-read MOTUs from the dataset before analyses. This filtering removed 83 sequences and left 96 MOTUs for further analyses, with an average length of 382 bp (range: 380–390 bp), representing a total of 20 953 raw reads across the 130 samples. Figure 1, which illustrates the distribution of the 96 MOTUs, was made in R v.4.2 using the packages tidyverse (Wickham et al. 2019) and readxl (Wickham and Bryan 2022).

Figure 1.

Graph showing the number of MOTUs registered at one or more sample sites. Thirty-six of the 96 MOTUs were found at only one of 130 sample sites.

Reference database

In order to assign taxonomic labels to MOTUs, a local eutardigrade V4 region 18S rDNA reference database was created (Supporting Information, Table S1). Specifically, 18S rRNA gene sequences from various eutardigrades were selected from GenBank on 1 July 2022. We excluded sequences directly submitted (unpublished), not included in any multigene studies, or using primers outside the V4 region from the reference database. The database encompasses 313 sequences, covering 48 genera, 14 families and two orders, with the shortest sequence (Mesobiotus ethiopicus, MF678793.1) covering 65.9% of the V4 region. Specifically, the four evolutionary lineages of Parachela, i.e. Hypsibioidea (four families, 19 genera), Macrobiotioidea (four families, 13 genera), Eohypsibioidea (one family, two genera) and Isohypsibioidea (four families, 13 genera), in addition to the genus Milnesium (order Apochela), are represented in this database. When excluding cryptic species (abbreviated aff., cf. or sp.), the database has representants of 184 of the 884 currently described eutardigrade species (Degma et al. 2022).

Creating a backbone tree

To create a backbone tree, the 313 reference sequences were aligned using MAFFT with the G-INS-i algorithm (% mafft, reorder, maxiterate 2, retree 1, globalpair input) and visually inspected using MEGA. An ML phylogeny was built with IQ-TREE v.1.6.12 (Nguyen et al. 2015a), with gaps and missing characters treated as unknown characters. The tree was rooted with two related lines of Heterotardigrada; Echiniscus testudo (GQ849022.1) and Batillipes mirus (GQ849016.1). GTR+F+I+G4 was chosen by ModelFinder (Kalyaanamoorthy et al. 2017) as the optimal substitution model, according to the Akaike information criterion. Branch support was calculated using 1000 SH-aLRT (Guindon et al. 2010) and 1000 ultrafast (UF) bootstrap replicates (Hoang et al. 2018). Using UF bootstraps in IQ-TREE produces a consensus tree (backbone.fas.contree) with assigned branch supports, in which branch lengths are optimized on the original alignment. This file was subsequently used as a reference topology (see Tree- and Phylogeny-based approach below) to reflect current consensus on eutardigrade phylogeny. The chosen reference topology did not influence any downstream analysis.

Alignment-based approach

The alignment-based approach estimates the similarity between MOTU and reference sequences. Specifically, MOTU sequences were compared with the reference sequences with VSEARCH v.2.21.1 (Rognes et al. 2016), using the --usearch_global command with a 90% identity threshold and ≥ 70% query coverage (--threads 4 --dbmask none --qmask none --rowlen 0 --notrunclabels --userfields query+id+target --maxaccepts 10 --maxrejects 0 --maxhits 10) with the option --iddet 2. A 100% match between a MOTU and a single reference sequence was considered sufficient for species-level determination. Similarity ranges of < 100% to ≥ 98%, < 98% to ≥ 95% and < 95% to ≥ 90% were used for assignments to genus, family and superfamily, respectively, as an adjustment to the threshold suggested by Wu et al. (2015).

Tree-based approach

The tree-based approach evaluates MOTUs by analysing their position relative to reference sequences within the backbone tree. The 96 MOTUs were aligned with the backbone alignment using MAFFT (% mafft -inputorder --keeplength --addfragments fragments --auto input) (Katoh and Frith 2017). An ML phylogeny was built with IQ-TREE, using the consensus backbone tree as topological constraint tree (-g backbone.fas.contree) and GTR+F+I+G4 as the substitution model. Problematic MOTUs that disrupted tree topology were removed before analysis. Branch support was calculated using 10 000 UF bootstrap replicates. The UF bootstrap values ≥ 95% were regarded as significant support (Zawierucha et al. 2020). To achieve an assignment to species, the MOTU had to share terminal nodes with a reference sequence and still show good support (UF bootstrap ≥ 95%). To achieve an assignment to genus, the MOTU should branch within a clade containing only members of this genus with ≥ 95% UF bootstrap support. For assignment to a higher-level taxon (family and subfamily), the MOTU had to branch within a clade containing primarily members of this taxon with ≥ 95% UF bootstrap support. MOTUs were assigned to superfamily level when they were placed within one of the four major monophyletic and highly supported clades.

Phylogeny-based approach

The phylogeny-based approach aims to find the optimal assignment of MOTUs onto a fixed backbone tree based on estimated support for distance-based placements. The phylogenetic placement was performed in Python 3.8.10 using APPLES v.2.0.9 (Hasan et al. 2022). This algorithm bases phylogenetic placement on distance, using a least-squares method to calculate differences between pairwise sequence and patristic distances on the backbone tree (Czech et al. 2022). The optimal placement of a MOTU on the backbone tree is therefore a placement on the branch that minimizes the difference between the calculated MOTU-to-reference distances and the distances present in the tree (Balaban et al. 2020, 2022). Before placement, branch lengths on the consensus backbone tree (backbone.contree) were re-estimated using FastTree v.2.1.11 (Price et al. 2010). The analysis was performed on 22 October 2022; it included fast bootstrapping with 1000 replicates (-F -N 1000) and printed only placements with the minimum least square error (--lse). The minimum least square error (MLSE), pendant lengths (PLs) and support values were evaluated manually from the output file (jplace-file; for details, see Matsen et al. 2012). Good placements were indicated by MLSE and PL of zero paired with support values ≥ 99%. A non-zero MLSE paired with a PL of zero indicated potential misplacements. For better visualization, the branch lengths were ignored on the phylogenetic placement tree produced by APPLES-2 owing to occasional negative branch lengths that follow the re-estimation using FastTree. Taxonomic assignment was performed by labelling each branch of the backbone tree by the most descriptive taxonomic path of its descendants and subsequently assigning each MOTU sequence to these labels based on placement (Czech et al. 2022).

Presentation of results and final taxonomic assignment

Tree illustrations were made with iTOL v.6.5.2 (Letunic and Bork 2021) and Adobe Illustrator CS6 v.16.0.0. Taxonomic assignment of MOTUs (Table 1) was conducted by comparing the results across the three methods, i.e. the alignment-based, tree-based and phylogeny-based approaches (Holovachov et al. 2017, Shchepin et al. 2017), using open nomenclature as suggested by Sigovini et al. (2016) and Horton et al. (2021). Accordingly, ‘inc.’ (incerta) was used to indicate that a MOTU was identifiable to a specific rank, without certainty, whereas ‘indet.’ (indeterminabilis) was used in connection with MOTUs that were unidentifiable below a certain rank. During the final taxonomic assignment, each MOTU was aligned with the best-matching reference sequences according to the three methods, and the alignments were inspected visually using MUSCLE v.3.8 (Robert 2004).

Table 1.

Identification of MOTUs based on the three methods, i.e. alignment-, tree- and phylogeny-based analyses (for details, see Figs 2, 3 and Supporting Information, Figs S1.1-1.4, S2.1-2.4; Tables S2 and S3)

Taxon	MOTU: Soil region (number of reads at distinct sampling sites)
Milnesiidae
Milnesium indet.	#074: Sand EJ (32)
Eohypsibioidea: Eohypsibiidae
Eohypsibius indet.	#044: Organic EJ (48), ZL (72)
Bertolanius indet.	#051: Sand EJ (20); Organic FLM (58); Clay WJ (10); #073: Organic WJ (32)
Hypsibioidea: Hypsibiidae
Hypsibiinae indet.	#001: Sand NJ (60, 25, 77), WJ (16, 94), ZL (124), FLM (19); Organic NJ (197), EJ (20, 39), FLM (10, 121, 134, 58); Clay WJ (407), EJ (38, 201), FLM (82); Sand clay ZL (45, 14, 17, 27)
	#009: Sand NJ (20), ZL (10), FLM (67); Organic NJ (15), WJ (35), ZL (61), FLM (13); Clay FLM (433); Sandclay ZL (27); #020: Sand ZL (80); Organic ZL (108), FLM (104); #041: Sand NJ (128); #045: Sand WJ (29, 43); Sandclay FLM (28) #052: Organic EJ (82); #063: Sand EJ (34); Clay FLM (22); #105: Sand ZL (14)
Hypsibius indet.	#002: Sand NJ (65, 13, 103, 18, 27), WJ (14, 67), EJ (14), ZL (15); Organic NJ (59), WJ (177, 24), ZL (28, 170, 37, 103), FLM (243, 205, 169, 86); Clay EJ (31, 38, 12); Sandclay WJ (36); #047: Sand NJ (26); Organic EJ (32); Clay WJ (17), EJ (13)
Diphascon indet.	#004: Sand NJ (28), WJ (40, 48); Organic NJ (130), WJ (43, 23), EJ (174), ZL (36, 49), FLM (173, 186, 83, 80); Clay EJ (16); Sandclay NJ (12), FLM (87); #008: Sand NJ (70, 10, 67), EJ (44); Organic EJ (31), ZL (31), FLM (91, 53, 27); Clay EJ (73); Sandclay WJ (52), EJ (19), ZL (230), FLM (11); #085: Sand WJ (11, 13)
Diphascon pingue inc.	#023: Sand NJ (28, 34, 48, 11), WJ (18); Organic NJ (23), FLM (66); Sandclay EJ (22), ZL (80)
Itaquasconinae indet.	#027: Sand NJ (42), WJ (20), EJ (31); Organic NJ (60), ZL (70, 18); #053: Sand WJ (33); Organic ZL (42); #056: Clay EJ (69); #072: Sandclay NJ (33); #084: Organic ZL (24)
Adropion indet.	#021: Sand NJ (18, 35, 67), WJ (74, 36); Organic EJ (58), ZL (86); #059: Sand NJ (62); #100: Sand NJ (16)
Adropion scoticum inc.	#017: Sand NJ (95, 36, 39, 37), WJ (51); Organic ZL (25), FLM (18); Sandclay WJ (29, 27, 62), EJ (40)
Astatumen indet.	#014: Sand NJ (38, 53), WJ (13, 31, 139), EJ (82); Sandclay WJ (55), EJ (30, 38), FLM (38); #066: Clay EJ (20), FLM (31)
Pilatobius indet.	#018: Sand NJ (20, 28, 32), ZL (21, 42); Organic FLM (24, 174); Sandclay ZL (13), FLM (33, 11); #019: Clay WJ (278), EJ (22); Organic ZL (91); Sandclay EJ (20); #037: Organic ZL (49); Sandclay NJ (10), EJ (18), FLM (30, 40); #087: Sand WJ (22); #106: Organic FLM (14) #112: Clay EJ (12)
Notahypsibius indet.	#036: Sand NJ (52), FLM (53); Organic EJ (31)
Mixibius indet.	#068: Sand NJ (29), ZL (19); #102: Sand ZL (15)
Mixibius saracenus inc.	#003: Sand NJ (23, 336, 109, 44); Organic WJ (276, 46, 120, 64), EJ (57), FLM (132)
Hypsibioidea: Microhypsibiidae
Microhypsibiidae indet.	#089: Organic ZL (21)
Microhypsibius indet.	#013: Sand NJ (25); Organic NJ (69), WJ (49), EJ (30), ZL (186), FLM (115); Sand clay NJ (52); #055: Organic FLM (18, 52); #093: Organic ZL (19); #114: Organic NJ (11)
Microhypsibius bertolanii inc.	#006: Sand NJ (26, 71, 72), WJ (48, 22), EJ (15), ZL (16), FLM (343); Organic WJ (44), EJ (38), ZL (15, 33), FLM (52, 53); Sand clay EJ (12), FLM (37)
Hypsibioidea: Ramazzottiidae
Ramazzottius indet.	#070: Sand EJ (19); Organic WJ (25)
Isohypsibioidea indet.	#033: Sand NJ (24), WJ (28); Organic NJ (82), ZL (16, 14)
Isohypsibioidea: Doryphoribiidae
Doryphoribiidae indet.	#026: Organic ZL (102, 105), FLM (41)
Doryphoribius indet.	#057: Clay ZL (66); #104: Clay ZL (12)
Doryphoribius macrodon inc.	#040: Clay ZL (24), FLM (65); Sand clay FLM (39)
Grevenius indet.	#010: Organic WJ (12, 48), EJ (28), ZL (16, 161); Sand clay NJ (25), FLM (295); #031: Sand NJ (32), WJ (29); #086: Organic FLM (23); #095: Organic FLM (18)
Thulinius indet.	#029: Sand ZL (25); Organic NJ (116), FLM (88); #054: Organic WJ (36), FLM (44); #061: Sand WJ (35); Clay WJ (26); #096: Organic EJ (17)
Thulinius augusti inc.	#042: Sand WJ (22); Organic EJ (67, 31)
Isohypsibioidea: Hexapodibiidae
Hexapodibius indet.	#028: Sand NJ (134), WJ (103)
Isohypsibioidea: Isohypsibiidae
Isohypsibiidae indet.	#015: Organic ZL (116, 404); #050: Organic FLM (89); #092: Organic FLM (19); #098: Sand NJ (16); #107: Organic FLM (13); #111: Organic FLM (13)
Dianea indet.	#024: Sand NJ (13), WJ (22), ZL (16); Organic EJ (20), ZL (19), FLM (29, 97, 30); Sand clay ZL (58), FLM (31); #067: Clay WJ (49); #069: Organic FLM (23, 22); #076: Sand EJ (12), ZL (13); #091: Organic ZL (14); #103: Organic FLM (14)
Dianea sattleri inc.	#035: Sand NJ (10), WJ (16); Organic NJ (11), FLM (51); Clay FLM (71)
Eremobiotus indet.	#079: Sand NJ (29)
Isohypsibius indet.	#071: Organic ZL (33); #077: Clay EJ (26)
Isohypsibius dastychi inc.	#038: Sand NJ (86, 23), FLM (49)
Isohypsibius cfr. dastychi inc.	#046: Sand NJ (106)
Ursulinius indet.	#011: Sand NJ (18); Organic NJ (30), EJ (47, 70), FLM (327); Clay ZL (34), FLM (54)
Ursulinius lunulatus inc.	#080: Sand ZL (13, 15)
Macrobiotoidea indet.	#065: Organic FLM (30, 24)
Macrobiotoidea: Adorybiotidae
Crenubiotus indet.	#058: Sand WJ (63)
Macrobiotoidea: Macrobiotidae
Macrobiotus indet.	#062: Clay EJ (58)
Macrobiotus caelestis inc.	#060: Organic EJ (62)
Macrobiotus hufelandi complex indet.	#030: Sand NJ (178), EJ (31)
Macrobiotus persimilis complex indet.	#082: Organic FLM (20)
Mesobiotus indet.	#081: Sand ZL (20); #083: Organic EJ (26)
Mesobiotus harmsworthi group inc.	#032: Sand NJ (61), WJ (26); Organic FLM (25); Clay EJ (13); Sand clay NJ (10), WJ (13, 11)
Minibiotus indet.	#005: Organic FLM (441, 513)
Paramacrobiotus indet.	#039: Sand WJ (52); Organic ZL (55); Sand clay FLM (12)
Macrobiotoidea: Murrayidae
Murrayidae indet.	#025: Sand NJ (31), WJ (12); Organic WJ (64), EJ (70), ZL (38), FLM (81)

Molecular operational taxonomic units (MOTUs) were assigned to various taxa using the open nomenclature provided by Horton et al. (2021). The 130 sampling sites (see Brunbjerg et al. 2019) are clustered into five regions: Northern Jutland (NJ); Eastern Jutland (EJ); Western Jutland (WJ); Zealand (ZL); and Funen/Lolland/Møn (FLM). Four different soil types were covered in each region (Sand, Organic, Clay and Sand clay). The number of reads at each site within each region is provided in parentheses. All taxa and species (except Macrobiotus caelestis) have previously been observed confidently in neighbouring (European) countries.

RESULTS

The majority of MOTUs are represented by relatively few reads (< 1000; Table 1) at a few sampling sites (Fig. 1), probably reflecting that tardigrades are very small and thus produce little eDNA. The 130 Biowide sampling sites cover a variety of different soil environments. Eutardigrade MOTUs were registered at all but 21 sites. The sites that apparently lacked eutardigrades were agricultural lands, dry environments and nutrient-poor wet environments. The diversity (number of MOTUs) is generally highest in wet organic soils, especially on the island of Lolland, with 47.9% of the total reads being from organic soils.

VSEARCH provided 10 reference sequences with the highest sequence similarity (based on percentage sequence similarity) to each of the 96 MOTUs. The best match, with the highest similarity and ≥ 70% query coverage, are reported in the Supporting Information (Table S2). Twenty MOTUs have a sequence similarity of 100%, with eight of these being to multiple reference sequences; 60 MOTUs have sequence similarities between < 100% and ≥ 98%; 14 MOTUs between < 98% and ≥ 95%; and two MOTUs have similarities ≤ 95%. Among the eutardigrade genera present in the reference database, some are represented by species with highly similar or identical sequences in the barcoding region (e.g. Milnesium, Dactylobiotus, Crenubiotus, Paramacrobiotus, Macrobiotus and Tenuibiotus).

The backbone tree generated in the present study is largely in agreement with previous 18S rDNA phylogenies (e.g. Guil and Giribet 2012, Guidetti et al. 2019a), in addition to those based on multiple markers (e.g. Bertolani et al. 2014, Dabert et al. 2014, Jørgensen et al. 2018, Gąsiorek et al. 2019, Gąsiorek and Michalczyk 2020, Guidetti et al. 2019b, Guil et al. 2019, Morek and Michalczyk 2020, Stec et al. 2020, 2022, Zawierucha et al. 2020, Tumanov 2022). As noted above, problematic genera included Milnesium, Dactylobiotus, Crenubiotus, Paramacrobiotus and Tenuibiotus, for which it was not possible to distinguish between single species within each genus. Moreover, the phylogeny of Mesobiotus was not in agreement with any recently published literature (Stec 2021, Stec et al. 2021b, 2022).

Below, we describe the results of the tree- and phylogeny-based approaches. As an example of the results obtained with these approaches, we show trees pertaining to Isohypsibioidea (Figs 2, 3). For the sake of simplicity, we do not show trees pertaining to the remaining taxa, but these are available for download in the Supporting Information (Figs S1.1–1.4, S2.1–2.4).

Figure 2.

Eutardigrada phylogenetic tree of the superfamily Isohypsibioidea reconstructed de novo with the obtained eutardigrade MOTUs (written in bold) and using the backbone tree as a constraint tree topology. The phylogeny was computed using IQ-TREE under the GTR+F+I+G4 model with 10 000 UltraFast bootstrap replicates. Only UltraFast bootstrap values ≥ 95% are shown. The scale bar and branch length refer to the maximum likelihood analysis.

Figure 3.

Phylogenetic placement within the superfamily Isohypsibioidea on a re-estimated backbone phylogeny using FastTree. All branch lengths were ignored for better visualization. Phylogenetic placement of MOTUs was performed using APPLES-2, which places the MOTUs based on calculated distances, followed by a phylogenetic correction using the Jukes–Cantor model. The MOTU identities are provided above each placement.

A ML phylogeny was constructed de novo using the 96 MOTUs aligned with the backbone alignment and using the backbone tree as a constraint topology (Fig. 2; Supporting Information, Fig. S1.1–1.4). Some MOTUs (#012, #015, #033, #050, #053, #068, #102 and #112) were problematic and caused severe topological changes when included in the analysis. To increase the support of the placement of surrounding MOTUs, these MOTUs were removed from the analysis and assigned to the family or superfamily level. When assigning taxonomy, MOTUs were considered to belong to a given taxon if the placement showed a high support in the tree (UF bootstraps ≥ 95%). If a MOTU was placed within a clade without sufficient branch support, the assigned taxonomy was raised to the higher-level taxon of the next well-supported branch in the tree.

The phylogenetic placement algorithm, APPLES-2, places MOTUs onto the backbone tree based on the calculated distances between the reference sequences and each MOTU, followed by a phylogenetic correction using the Jukes–Cantor model. Only placements with the highest bootstrapping support values are printed on the tree, and branch lengths are ignored (Fig. 3; Supporting Information, Fig. S2.1–2.4). The MLSE, bootstrap support values and PLs for each MOTU are provided in the Supporting Information (Table S3). The MSLE illustrates good placements at zero (Balaban et al. 2022), and 26 MOTUs show an MLSE and PL of zero. Of these, eleven have support values between 100% and 99%; nine have support values of ≥ 75%; and the remaining six have support values ≤ 75%. Seven MOTUs show non-zero MLSE and PL and ≥ 75% support value. Forty-four MOTUs are placed on inner branches, which suggests that these MOTUs were hard for the algorithm to place. Fifteen MOTUs showed PLs of zero together with a non-zero MLSE, which can be a sign of placement error (Balaban et al. 2022).

Final taxonomic assignment of MOTUs was found by comparing the results across the three approaches, with subsequent manual validation of sequence identity (Table 1). The results uncover a broad diversity within the Danish soil samples, with well-supported assignments to Adropion scoticum Murray, 1905, Dianea sattleri Richters, 1902, Doryphoribius macrodon Binda, Pilato & Dastych 1980, Diphascon pingue Marcus, 1936, Isohypsibius dastychi Pilato, Bertolani & Binda, 1982, Macrobiotus caelestis Coughlan, Michalczyk & Stec, 2019, Mixibius saracenus Pilato, 1973, Microhypsibius bertolaniiKristensen, 1982, Ursulinius lunulatus Iharos, 1966, Thulinius augusti Murray, 1907 and the genera Astatumen, Bertolanius, Crenubiotus, Eohypsibius, Eremobiotus, Grevenius, Hexapodibius, Notahypsibius, Milnesium, Mesobiotus, Minibiotus, Paramacrobiotus, Pilatobius and Ramazzottius.

When comparing the results with traditional surveys in neighbouring countries, all the proposed species (except Macrobiotus caelestis), genera and families found in the Danish eDNA have been observed confidently, which would lend support to their presence also in Denmark (Dastych 1988, Guidetti and Bertolani 2001, Guil and Giribet 2012, Bertolani et al. 2014, Guidetti et al. 2015, DeMilio et al. 2016, Bingemer and Hohberg 2017, Meier 2017, Vuori et al. 2020, Topstad et al. 2021). As noted above, of the 96 MOTUs, 10 MOTUs had well-supported assignments to species level in all three analyses (see Table 1). Specifically, these species, with the exception of Macrobiotus caelestis (seeCoughlan et al. 2019), have all been observed confidently in one or more of the following neighbouring countries: Sweden (Guidetti et al. 2015), Norway (Meier 2017, Topstad et al. 2021), Finland (Vuori et al. 2020), Poland (Dastych 1988), Ireland (DeMilio et al. 2016), Italy (Guidetti and Bertolani 2001, Bertolani et al. 2014) and Spain (Guil and Giribet 2012). Additionally, most of these species have also been recorded previously in soil or sediment samples (Kristensen 1982, Bingemer and Hohberg 2017). Generally, the results indicate that some eurytopic eutardigrades (Guil et al. 2009) inhabit the Danish soil, including Dianea sattleri, Adropion scoticum and a member of the Macrobiotus hufelandi complex. Freshwater eutardigrades, such as Mixibius saracenus and Isohypsibius dastychi, have been collected before in sediment samples (Bertolani et al. 2014); therefore, it is not unlikely for these freshwater tardigrades to occur in moist soil. Generally, the results from this study represent an unprecedented addition to the otherwise poorly confirmed list of tardigrades in Denmark.

DISCUSSION

The most important step in the process of metabarcoding is the taxonomic assignment of MOTUs. Within the metabarcoding field it is not uncommon that MOTUs cannot be identified, and thus could potentially constitute undescribed species, for which no references are available. This is especially evident in meiofaunal surveys, for which large parts of the communities remain to be described properly, and the available reference databases do not exhibit sufficient taxon coverage (Tang et al. 2012). The 18S rDNA database generated in this study aimed to cover as many eutardigrade taxa as possible, whilst keeping a good quality of the reference sequences. When excluding cryptic sequences included in the database (i.e. sequences annotated as ‘sp.’, ‘aff.’ or ‘cf.’), the database covers ~21% of the currently described eutardigrades (Degma et al. 2022). This poses an obvious problem for taxonomic assignment of eutardigrade MOTUs based on 18S rDNA. It follows that an effective application of eDNA metabarcoding on tardigrades requires that more reference sequences are retrieved (18S rDNA and other markers) (Topstad et al. 2021). Hence, efforts should be made to sequence a much greater number of the described tardigrade species and include these in the database. Importantly, heterotardigrades have large insertions within the V4 region (Jørgensen et al. 2010), probably explaining the lack of heterotardigrades within the present dataset. This implies that heterotardigrades are best represented by heterotardigrade-specific primers or metabarcoding of short sequences of COI mtDNA.

Taxonomic resolution depends on the assignment approach, and all three methods applied in the present study suffered loss of taxonomic knowledge when applying the threshold parameters to the results. The three approaches generally agreed on the higher taxonomic levels (family and subfamily) but did not always agree at the genus and species level. Sometimes, a genus-level identity could not be obtained by one approach but was suggested by another approach, making the results harder to interpret. In most cases, assignments to genus could not be obtained because: (1) the support fell outside the threshold parameters; (2) some taxa have low interspecific diversity in the barcoding region, resulting in MOTUs with multiple equally good identifications; and (3) MOTUs placed within para- and polyphyletic clades are hard to evaluate.

The alignment-based approach is fast; however, the results are limited by the number and variation of references in the database. Specifically, many eutardigrade reference sequences are similar or nearly identical in the barcoding region, resulting in MOTUs with multiple and equally good matches (Supporting Information, Table S2). Specifically, selected taxa are represented by species with completely identical barcoding regions (e.g. Grevenius, Paramacrobiotus, Milnesium and Murrayidae), whereas others are distinguishable, suggesting that the level of variability of the barcoding region differs among various eutardigrade taxa. Applied similarity cut-offs were based on a previous study on the 18S rRNA V4 region (Wu et al. 2015), and MOTUs with similarities < 100% were accordingly assigned to genus or a higher taxon (Table 1). To improve this method, similarity thresholds should be based on knowledge of interspecific diversity within each of the various eutardigrade clades (Bik et al. 2012). This highlights the importance of a critical evaluation of the output when applying alignment-based approaches to place MOTUs to supraspecific taxa (genus, subfamily or family). For MOTUs with lower sequence similarity (< 99%), it is not unusual for the most similar reference sequence to not always represent the most closely related species (Koski and Golding 2001). This is evident from our results, where the similarity search was not always in agreement with the methods inferring phylogeny (e.g. MOTU #026 with the highest similarity to Apodibius confusus KC582830.1; for details, see Supporting Information, Table S2). Furthermore, similarity searches can fall short, because sequence similarities do not reveal whether MOTUs with identical sequence similarity to a given reference sequence are highly related to each other, or merely exhibit and equal relationship to the same reference sequence. Reliance on similarity searches for assigning identity to MOTUs can also be problematic in cases where MOTUs show the same sequence similarity to multiple non-identical references (e.g. MOTU #072). In these cases, using methods that infer phylogeny can be an advantage to determine the evolutionarily most likely match.

Phylogenetic reconstructions have greater computational requirements than similarity searches. Specifically, phylogenetic reconstruction present results in a tree-based format, allowing a comparison of the relationship between all sequences, including subsets of MOTUs grouped into distinct clades. Interestingly, our results reveal a well-supported clade consisting of only MOTUs (MOTU #041, #045, #052 and #063, respectively, for details, see Supporting Information, Figure S1.2). A serious downside of using phylogenetics on metabarcoding data is the low phylogenetic signal in MOTUs owing to their short length. The barcodes are chosen purposely because they are short in length and cover a variable region flanked by conserved regions (Holovachov 2016b), and these can be difficult to align using progressive alignment algorithms (Holovachov 2016b). Maximum likelihood phylogenies are very sensitive to changes in the multiple sequence alignment (Nguyen et al. 2015b), and small changes are likely to happen during the alignment of the references and MOTUs (Holovachov et al. 2017). The amount and composition of MOTUs can have a big effect on tree topology and, consequently, branch support, causing neighbouring MOTUs to be misidentified. Some MOTUs had to be excluded from the present ML analysis because they caused severe disruptions in the tree topology and branch support. Another likely result of the sensitive alignment is that some MOTUs, with high similarity to a specific reference sequence, were not placed accordingly on the tree (e.g. MOTU #027). This supports the recommendation by Holovachov (2016a) to use more than one combination of alignment and phylogeny inference algorithms to identify MOTUs reliably. In the present study, the backbone alignment created using MAFFT showed a high frequency of small gaps, which might decrease the accuracy in subsequent ML computed trees (Nguyen et al. 2015b). In future analysis, it is worth considering re-alignment of the reference database using some of the more computationally heavy alignment algorithms, such as UPP (Nguyen et al. 2015b), and subsequently performing de novo reconstruction using RAxML (Stamatakis 2014). Alternatively, IQ-TREE also offers several options of branch support, based on computationally heavier analysis. An obvious extra approach could be to validate the ML method with Bayesian inference.

Nevertheless, for MOTU sample identification, heavy de novo reconstructions might not necessarily be needed, because simpler placements on a backbone tree can be sufficient when supraspecific taxa are well resolved and supported in the backbone tree. The phylogeny-based approach allows for fast estimation of the most likely position of each MOTU within a constrained backbone tree (Holovachov et al. 2017, Czech et al. 2019), incorporating models of sequence evolution, thus providing additional information in comparison to similarity searches (Balaban et al. 2020). The tree-based approach adds new branches to the backbone tree and resolves the phylogenetic relationships between the MOTUs, whereas the phylogeny-based approach keeps the backbone tree fixed, mapping MOTUs onto existing branches in the backbone tree (Balaban et al. 2020, Czech et al. 2022). The quality of each MOTU does not affect the accuracy of the results, i.e. potential erroneous sequences will not change placement of other MOTUs (Holovachov et al. 2017, Hasan et al. 2022). Phylogenetic placement algorithms, such as APPLES-2, thus treat all MOTUs independently, allowing fast analysis and taking only seconds to complete. The downside is that the phylogenetic relationships between individual MOTUs will not be resolved. Hence, as for alignment-based approaches, the phylogeny-based approach only compares MOTUs with references included in the database and does not consider the relationships among the MOTUs. MOTUs belonging to taxa not included in the backbone tree are forced onto branches based on MLSE. This underlines the fact that the accuracy of the results can only be as good as the accuracy of the reference database. Phylogenetic placement algorithms cannot infer evolutionary relationships below the taxonomic level of the backbone tree (Czech et al. 2022), and the taxon composition and sequence quality determine which taxa can be identified and which taxa cannot (Holovachov et al. 2017). Although the phylogeny-based approach has its disadvantages, it has the potential to supplement the commonly used alignment-based approach and to replace heavy de novo reconstructions once the backbone phylogeny of tardigrades includes broader taxonomic coverage, making comparisons between studies more approachable.

In summary, we find that the commonly adopted alignment-based approach is limited by the low interspecific variation found within the V4 barcoding region of the 18S rRNA gene among selected eutardigrade taxa. Both alignment- and phylogeny-based approaches currently suffer from incomplete taxon coverage. Applied thresholds directly affect the results obtained from all approaches, but taxonomical levels are artificial units, and it is therefore impossible to find a perfect genetic similarity threshold adequate for all situations (Shchepin et al. 2017). The tree-based approach was least affected by the insufficient taxon coverage in the database and was able to resolve relationships between closely related MOTUs. However, this approach also represents the most computationally heavy analysis, and it can easily be affected by poor alignment and quality and quantity of the input MOTU data, causing disruptions in the backbone tree and, eventually, uninterpretable results. The optimal assignment of MOTUs was achieved through a combination of the approaches used in the present study. Finally, we note that the MOTUs analysed in the present study derive from a larger dataset aimed at exploring the diversity of various eukaryotes across Denmark. The primers used are universal for eukaryotes, and the original dataset is dominated by reads that belong to taxa other than tardigrades. Moreover, semi-terrestrial tardigrades often inhabit moss and lichens, which are not represented in the present dataset. This is likely to explain the low read count and that the different tardigrade MOTUs often were found at a few sampling sites. Thus, although the present study provides the first insight into the hidden diversity of eutardigrades in Denmark, future studies would benefit from exploration of a wider range of habitats and an implementation of primers and sampling protocols developed specifically for tardigrades.

CONCLUSIONS

Owing to their microscopic size, tardigrades are hard to monitor using traditional biomonitoring methods. eDNA has the potential to fill in some of these gaps and to broaden our knowledge of species boundaries and geographical ranges (Topstad et al. 2021, Saccò et al. 2022). Future studies could supplement the present V4 data by barcoding soil samples using other 18S regions (e.g. V9) and well-known markers, such as ITS2 and COI. Ideally, field studies should test whether the taxa provided in Table 1 can be found at the sampling sites (Brunbjerg et al. 2019) using traditional methods.

Investigation of community correlations could lay the basis for solid information on which tardigrades are widely distributed and abundant and which tardigrades are more likely to be endemic, with narrow distributions. Identifications to higher taxonomic levels (genus and family) could, potentially, provide enough information to group eutardigrades into trophic and functional groups for further ecological assessments, which would be likely to reveal the otherwise unknown diversity and distribution of different tardigrade taxa. More specifically, future studies could analyse the community composition at the different sampling sites to see whether this is correlated with the ecological gradients covered in Biowide (soil moisture, soil fertility and succession) and how it relates to plant communities and other environmental factors. Tardigrades frequently inhabit mosses, lichens and dead organic matter, hence expanding the dataset to include organic matter from such sample sites would be likely to yield a much higher number of reads of each MOTU, and thus provide a higher reliability to some of the dubious MOTUs.

The optimal assignment of MOTUs was best achieved through a combination of taxonomic assignment approaches. Although most MOTUs could not be identified to the species level, these results still provide new information that adds to our current understanding of tardigrade biodiversity in Denmark. With additional studies, this combined method of MOTU assignment has a great potential to serve as a model for the investigation of tardigrade distribution in similar and different environments in other parts of the world.

Supplementary data

Supplementary data is available at Zoological Journal of the Linnean Society online.

File S1. Excel file with detailed information on MOTUs and collection sites (Biowide).

Table S1. Appendix of Genbank references used in the reference database.

Table S2. Results from VSEARCH, showing matches with the highest similarity to MOTUs.

Table S3.Results from APPLES-2, showing Minimum Least Square Error, Pendant length, and Bootstrap support values.

Figures S1.1-1.4. Maximum likelihood phylogenetic tree of Parachela, Eohypsibioidea, Hypsibioidea and Macrobiotoidea reconstructed de novo with the obtained eutardigrade MOTUs (written in bold) and using the backbone tree as a constraint tree topology. The phylogeny was computed using IQ-tree under the GTR+F+I+G4 model with 10.000 UltraFast bootstrap replicates. Only UF bootstrap values ≥ 95% are shown. Scale bar and branch length refer to the Maximum Likelihood analysis.

Figures S2.1-2.4. Phylogenetic placement within Parachela, Eohypsibioidea, Hypsibioidea and Macrobiotoidea on a re-estimated backbone phylogeny using FastTree. All branch lengths were ignored for better visualization. Phylogenetic placement of MOTUs was performed using APPLES-2, that place the MOTUs based on calculated distances followed by a phylogenetic correction using the Jukes-Cantor model. MOTU-IDs are provided on each placement.

ACKNOWLEDGEMENTS

We thank the Biowide team for providing the data. We thank Metin Balaban, Torbjørn Rognes and Henrik Riskær for helping with the algorithms used in this study. We also thank the reviewers for their suggestions and comments. Special thanks are owed to Łukasz Michalczyk, chair of the 15th International Symposium on Tardigrada.

CONFLICT OF INTEREST

The authors confirm there are no conflicts of interest.

DATA AVAILABILITY

Environmental data related to each of the 130 Biowide sites are available at: https://ecos.au.dk/forskningraadgivning/temasider/data. The 18S metabarcoding dataset is available at Github: https://github.com/tobiasgf/sample_storage/tree/main/asv_tables_referencedata (file: bw_tar_seqtab.nochim_Both.rds). The script used for extracting potential tardigrade sequences from the full 18S dataset (filtering_of_sequences.r) and raw data files are available at: https://github.com/fridalp/eutardigrade_eDNA