Applicability and perspectives for DNA barcoding of soil invertebrates

Introduction

Soils are among the most diverse habitats on earth, harbouring 25% to 50% of the biodiversity on Earth (Decaëns et al., 2006; Decaëns, 2010; Anthony, Bender & van der Heijden, 2023). This biodiversity drives essential processes for life on Earth and provides ecosystem services that impact human wellbeing, such as the decomposition of dead organic material, recycling of nutrients and carbon storage (Wardle et al., 2004; Lavelle et al., 2006; Bardgett & van der Putten, 2014). Characterizing and monitoring soil biodiversity therefore is of general interest to maintain and preserve soil functions (Orgiazzi et al., 2015). However, this is a challenging and time-consuming task due to the enormous taxonomic diversity and cryptic lifestyles of soil-organisms. Molecular methodologies offer great advantages for soil biodiversity assessment in terms of time and cost efficiency, and taxonomic resolution (Antil et al., 2022; Eisenhauer, Bonn & Guerra, 2019).

A large fraction of soil animal biodiversity is represented by microarthropods with body-sizes between 0.1 and 2 mm. Collembola (Hexapoda) and Oribatida (Acari: Sarcoptiformes) are dominant and omnipresent microarthropod taxa, and occur in all soil-related habitats where they reach high abundances of up to 50,000–100,000 individuals per square meter (Bardgett & van der Putten, 2014). Traditionally, Collembola and Oribatida have been described as decomposers, microbivores and fungivores, but studies using stable isotopes showed that they actually cover several trophic levels, demonstrating trophic specialization and functional diversity within these taxa (Schneider et al., 2004; Pollierer et al., 2009; Potapov et al., 2016; Maraun et al., 2023). These microarthropods spend their entire life in the soil matrix or in the litter layer, which makes them interesting candidates as bioindicators of soil quality in monitoring programs (Gulvik, 2007). Collembola are typical r-strategists with fast reproduction cycles, whereas Oribatida are usually considered K-strategists with long-life spans of 1–3 years and low fecundity, but species with shorter life-cycles are also common (Maraun & Scheu, 2000; Pfingstl & Schatz, 2021). The general differences in life-history traits and trophic diversity between Collembola and Oribatida could be informative for monitoring programs. Collembola respond and recover more quickly to disturbances (Ponge et al., 2003; Santorufo et al., 2012) than Oribatida, which have long recovery times and therefore are more sensitive to environmental changes (Zaitsev et al., 2002; Gulvik, 2007; Pfingstl & Schatz, 2021). However, the wide range of functional and life-history traits among different species necessitates species level determination in order to better understand their interactions in the soil system or to use them as bioindicators for changes in soil functions. About 9,000 species of Collembola and 11,000 species of Oribatida are described worldwide, but this likely represents only about 20 % of the expected species (Potapov et al., 2020; Behan-Pelletier & Lindo, 2023). Local species richness of these two taxa can be very high, reaching 60–100 species in forest soils (Rusek, 1998; Schatz & Behan-Pelletier, 2008). High species richness and abundance, and small body sizes of both, Collembola and Oribatida, pose a significant challenge for biodiversity assessments. Molecular applications, such as DNA barcoding and metabarcoding, have great potential to aid specimen identification and biodiversity assessment (Valentini, Pompanon & Taberlet, 2009). These methods utilize a standardized DNA fragment for taxonomic assignment of specimens by matching DNA sequences of undetermined individuals to a reference database (Hebert et al., 2003; Hebert & Gregory, 2005). This enables to automatically assign any taxonomic level and even species names to undetermined individuals. It is applicable to mixed samples of pooled specimens, which significantly reduces workload and costs. Further, molecular identification tools are equally applicable to juveniles that often lack taxonomic characters (Richard et al., 2010; Grzywacz et al., 2021). Automated handling of samples, simultaneous identification of multiple individuals in a single reaction, and the scalability of molecular data to any taxonomic level offers new opportunities for analysing spatial and temporal dynamics of soil-living animals (Arribas et al., 2021; Decaëns, 2021), and thereby provide new perspectives for monitoring of soil biodiversity. The method, however, relies on two preconditions: (1) a representative reference database and (2) a marker (barcoding) gene that reliably separates species. The most common databases are BOLD (“The Barcode of Life Data System”, Ratnasingham & Hebert, 2007) and NCBI (https://ncbi.nlm.nih.gov/). The standard barcoding gene for Metazoa is a 658 bp region of the mitochondrial cytochrome oxidase I gene (COI; Hebert et al., 2003; Hubert et al., 2008). In general, a minimum of 500 bp of COI is required, but shorter fragments can also be used for specimen identification and species discovery (Hajibabaei et al., 2006; Collins & Cruickshank, 2012).

Two types of methods have been developed for molecular species assignment. Distance-based methods, such as the barcoding gap (Hebert et al., 2003), the BIN system of BOLD (Ratnasingham & Hebert, 2013) and Neighbor Joining (Saitou & Nei, 1987; Hebert et al., 2003), which transform sequence alignments into a genetic distance matrices. These genetic distances can be calculated based on the observed differences between sequences, or by including a model of sequence evolution that accounts for mutational processes. Based on a threshold value of similarity, these distances are then used to assign sequences to known species or to identify putatively new species. Alternatively, character-based methods, such as Maximum Likelihood, GMYC (Pons et al., 2006) and PTP (Zhang et al., 2013), rely on a phylogenetic tree. These approaches use DNA sequences directly without a distance matrix and information on character evolution is not getting lost by capturing differences among sequences in a single metric. These methods infer species based on branching frequencies in a time-calibrated tree (GMYC) or by comparing the number of substitutions on branches (PTP) using all characters in an alignment. Reciprocal monophyly of taxa on a Maximum Likelihood or Bayesian inference tree is also appropriate to check if sequences can be assigned to known species. The dependence on a phylogenetic tree, however, makes them computational more demanding compared to distance-based methods. Character-based methods also require a threshold for delineating species boundaries and a model of sequence evolution to infer phylogenetic relationships. In DNA barcoding, the K2P model is the most widely used (Hebert et al., 2003; Nishimaki & Sato, 2019), and can be viewed as a compromise between observed genetic distances that do not account for evolutionary changes and complex models which might overestimate genetic distances in closely related taxa. However, the use of K2P as a standard model has also been criticized, but identification success is hardly affected, even if K2P poorly fits as model for a dataset (Srivathsan & Meier, 2011; Collins et al., 2012).

The success of species delineation based on genetic markers depends on the presence of a threshold value, also known as the barcoding gap, which implies that genetic distances within a species are smaller than genetic variances to congeneric and other species (Meyer & Paulay, 2005). A global threshold of 2%–3% intraspecific sequence divergence, or 10x the mean intraspecific divergence, has been proposed to reliably separate species (Hebert et al., 2004). Such a universal threshold is extremely helpful for automated species assignment of genetic data in bioinformatic pipelines. This threshold seems to be valid for a range of taxa (Hebert, Ratnasingham & de Waard, 2003; Hebert et al., 2004; Barrett & Hebert, 2005), but its universal application has been questioned for other species (e.g., Burns et al., 2007; Chapple & Ritchie, 2013; Elias et al., 2007; Meier et al., 2006; Meyer & Paulay, 2005; Wiemers & Fiedler, 2007). In particular soil-living animals show high intraspecific divergences in the COI gene that commonly exceed the standard barcoding threshold. Examples cover different families of earthworms (King, Tibble & Symondson, 2008; Novo et al., 2010; Martinsson, Rhoden & Erseus, 2016), Collembola (Porco et al., 2012a; Porco et al., 2012b; von Saltzwedel, Scheu & Schaefer, 2016; Zhang et al., 2019) and Oribatida (Rosenberger et al., 2013; von Saltzwedel et al., 2014). These studies question the general effectiveness of COI for specimen identification in these taxa. Moreover, asexual reproduction occurs in 7–10% of all species in several families of Collembola and 10% of all species of Oribatida (Chahartaghi, Scheu & Ruess, 2006; Cianciolo & Norton, 2006; Chernova et al., 2010; Bluhm, Scheu & Maraun, 2016), and asexual species can be dominant in temperate forests (Maraun & Scheu, 2000). According to theory, asexual organisms accumulate mutations over time, until they go extinct due to the accumulation of too many deleterious mutations (Muller’s ratchet: Muller, 1964; Kondrashov’s hatchet: Kondrashov, 1988). This suggests that present day populations of asexual species represent a range of COI haplotypes, while populations of sexual species should represent discrete clusters of similar COI haplotypes, standing for independently evolving lineages that interbreed (Barraclough, Birky jr & Burt, 2003). In consequence, a barcoding gap should not be present in asexual species but rather a continuum of slightly divergent individuals. Further, hybridization events are potential origins of asexual species, and if followed by mitochondrial introgression the detection of a barcoding gap is difficult (Mutanen et al., 2016; Dupont et al., 2016). Altogether, asexual reproduction could blur lines of species identification, and a hybrid species could be wrongly identified as its maternal species.

The aims of this study were to test (i) if the standard barcoding marker gene COI meets the precondition to reliably assign species in Collembola and Oribatida, and (ii) the accuracy of separating species based on a barcoding gap. Many species of these taxa have wide distribution ranges, and European species often occur across Palaearctic or Holarctic regions. Geographic coverage of samples provided in DNA barcode reference libraries can affect species assignment (Hebert et al., 2003). We therefore focused on species with a dense and broad geographic sampling to cover the potential range of intra-specific haplotype variation of COI (Phillips, Gillis & Hanner, 2022). We also included parthenogenetic (asexual) species to test (iii) if the reproductive mode affects the barcoding gap, because parthenogenetic species likely carry a continuum of divergent haplotypes due to the accumulation of mutations and the absence of homogenizing effects of mixis. We analysed the performance of COI for species delimitation using distance- and character-based methods.

Materials & Methods

Results

Datasets included 970 Oribatida and 612 Collembola sequences of COI, and alignments were between 507 bp and 657 bp long (Table 2), and covered the standard barcoding region of COI. Only a few sequences were below 500 bp long, predominantly in the Oribatida genus Steganacarus and the Collembola genus Folsomia.

Barcoding gap threshold detection for different genetic distances

The global datasets (all Oribatida, all Collembola) had relatively high cumulative errors (false positives and false negatives, Fig. 1; Table 3). The optimized local barcoding threshold for the local datasets differed among taxa (Fig. 1; Table 3). One dataset had a narrow threshold without species mismatches (Achipteria, 15%), others had large threshold ranges without mismatches (Nothrus, Platynothrus and Ceratophysella), and for the remaining datasets it was not possible to define a barcoding threshold without any mismatches (Oppiidae, Steganacarus, Folsomia). The optimal barcode thresholds at genus level exceeded the standard barcoding threshold of 2%–3% by 1% (Platynothrus) and up to 6% (Achipteria), except in the two genera Nothrus and Ceratophysella.

Table 3:

Range of barcoding gap thresholds and cumulative errors for all datasets.

The cumulative error is the sum of false positives and false negatives. Except for Nothrus and Ceratophysella the barcoding gap is not present in the investigated species, or exceeds the standard threshold of 2%–3%.

	Optimal barcoding gap threshold
	Cumulative error = 0	Smallest cumulative error
All Oribatida	–	4% (error: 26)
Achipteria	15%	9%–14% (error: 1)
Steganacarus	–	8%–15% (error: 8)
Nothrus	2%–17%	1%, 18% (error: 2)
Platynothrus	>4%	2%–3% (error: 1)
Oppiidae	–	8%–15% (error: 3)
All Collembola	–	6%–11% (error: 9)
Ceratophysella	3%–11%	2% (error: 3)
Folsomia	–	6%–14% (error: 9)

DOI: 10.7717/peerj.17709/table-3

Distance-based barcoding gap with ASAP

Accurate specimen assignment requires a gap between the largest genetic distance within and the smallest distance between species. We compared the distribution of intra- and interspecific distances of nominal species with that of the genetic lineages inferred by ASAP (Fig. 2), selecting the partitions with the least number of subsets. Our analysis consistently demonstrated that genetic distances of COI within and between morphologically assigned species overlap, which makes accurate species assignment impossible. The parthenogenetic oribatid mite genus Nothrus was a single exception, which showed a clear barcoding gap for morphologically assigned species. When genetic lineages (ASAP subsets) were considered, a barcoding gap between intra- and interspecific distances was present. Overall, the assignment of genetic lineages to morphospecies reduced the overlap of intra- and interspecific genetic distances considerably in all datasets, However, the effect was much more pronounced in Collembola than Oribatida and generated a barcoding gap that spanned a range of more than 10% between intra- and interspecific genetic distances. Among Oribatida, the effect of splitting morphotypes into genetic lineages was not strong. However, for two Oribatida datasets, Achipteria and Platynothrus, the choice of using the partition with the lowest number of subsets was too conservative because the resulting barcoding gap was very narrow. The barcoding gap thresholds estimated for the partition with the lowest number of subsets ranged in Oribatida 6.9% (Achipteria), 12.2% (Nothrus), 15.0% (Steganacarus, Platynothrus, all Oribatida) and 15.8% (Oppiidae); in Collembola 8.0% (Folsomia), 11.3% (all Collembola) and 13.8% (Ceratophysella). Notably, the intraspecific distances of the genetic lineages of Platynothrus show three clusters in distribution frequencies (<3%, at 3%–8%, 14%–17%) that likely represent three genetic lineages in P. peltifer (Table 1).

The outliers, i.e., single data points scattered within the range of the barcoding gap, likely belong to sequences that were considerably shorter than the average sequences. Both Collembola datasets were more heterogeneous in sequence lengths than the Oribatida datasets. Only the datasets of Steganacarus and Oppiidae also had very short sequences compared to the median sequence lengths, and both also had outliers after splitting morphotypes into genetic lineages. A few outliers remained for the genus Nothrus, which likely belonged to the species N. palustris and N. pratensis. After splitting both species into two genetic lineages as proposed by ASAP, nearly all outliers disappeared.

Character-based specimen assignment with Maximum Likelihood

Reliability of specimen assignment based on phylogenetic inference and therefore on molecular characters was very poor in Collembola (Fig. 3). The two genera Ceratophysella and Folsomia and the species within each genus were not monophyletic. Species clustered within clades of other species several times. The topology of Oribatida supported monophyly for most genera (Fig. 4). The species in the genus Nothrus and Platynothrus were also monophyletic. Only Nothrus pratensis separated into two highly supported, non-monophyletic clades. Within the genus Steganacarus all species were monophyletic, except S. magnus which formed five clades, two with 100% bootstrap support. The species S. carinatus was monophyletic, but separated into two highly supported clades. One sequence of S. applicatus clustered within S. carinatus while the remaining sequences were monophyletic with very high support. Possibly this single sequence represents a misidentified individual. Most genera of Oppiidae were monophyletic, except for Disorrhina and Oppia. The sequences assigned to Oppia sp. were sister to one clade of O. nitens with very high bootstrap support. It is possible that these sequences belong to the species O. nitens. All remaining species were monopohyletic with 100% bootstrap support. The genus Achipteria was represented by only three species. The two species A. howardi and A. catskillensis were monophyletic, but the sequences of A. coleoptrata were non-monophyletic.

Representativeness of sampling effort

Rarefaction curves (Fig. 5) showed that Oribatida species had more haplotypes than Collembola and that sexual Oribatida species (A. coleoptrata, S. magnus) had more haplotypes than parthenogenetic Oribatida (P. peltifer, N. silvestris). Further, Collembola reached saturation in species diversity at a sampling size of less than 200 individuals (for COI and 28S rDNA), the pattern was similar for the parthenogenetic Oribatida N. silvestris. However, the parthenogenetic Oribatida P. peltifer and both sexual Oribatida species did not reach saturation at a sampling size of more than 600 individuals and the expected diversity exceeded 200 haplotypes.

Discussion

This study tested the validity of a barcoding gap and the applicability of the standard barcoding gene COI for species assignment in two of the most species rich and abundant taxa of soil-living invertebrates, Collembola and Oribatida. The analysed datasets comprised two genera of Collembola with eleven and nine species, respectively. Oribatida datasets comprised four genera with three to six species per genus, and one family-level dataset with ten species in six genera.

Our results showed that correct species assignment was possible within some genera, but not all. However, both distance- and character-based methods were not able to assign species without mismatches when all Collembola or all Oribatida were analysed together. This is likely due to the different ranges of intraspecific genetic distances, demonstrating the absence of a general (global) barcoding gap for COI in these taxa. The genetic divergence separating intra- and interspecific distances differed among taxa and exceeded the standard species threshold of 3% intraspecific genetic distance in all but one species, indicating that taxon-specific thresholds should be applied for correct specimen assignment (Phillips, Gillis & Hanner, 2022). Here, the application of algorithms that dynamically adjust thresholds for sequence clusters, and therefore apply flexible thresholds, could improve species assignment in soil invertebrates (James, Luczak & Girgis, 2018; Chiu & Ong, 2022).

Absence of a global barcoding gap in the COI gene seems to be particularly relevant for soil-living animals and hampers the application of automated specimen assignments in DNA-based biodiversity surveys such as metabarcoding. Absence of a global barcoding gap had also been demonstrated for Annelida, among which earthworm taxa accounted for one third of interspecific comparisons with 0% genetic divergence (Kvist, 2016). In metabarcoding studies, Collembola had a high failure rate and high numbers of false positives for species assignments based on public databases and COI (Recuero, Etzler & Caterino, 2023). Among mites, specimen assignment is in general correct at least to family and order level (Oliverio et al., 2018; Ustinova et al., 2021; Young, deWaard & Hebert, 2021; Young & Hebert, 2022; Recuero, Etzler & Caterino, 2023). General explanations for failures in species assignments include the lack of completeness and misidentified individuals in reference databases, geographic underrepresentation of species and a neglect of assigning genetic lineage identities to sequences in reference databases (Kvist, 2016; Martinsson, Rhoden & Erseus, 2016; Young et al., 2019; Young, deWaard & Hebert, 2021; Phillips, Gillis & Hanner, 2022; Recuero, Etzler & Caterino, 2023). The rarefaction analysis demonstrated that genetic diversity is exceptionally high within morphospecies of soil-living invertebrates, and more genetic diversity is to be expected in additional samples. In particular, rarefaction curves for Oribatida did not reach saturation without sampling hundreds of additional individuals. For Collembola the number of expected COI haplotypes is lower as curves reached saturation at an expected sampling size of about 200 individuals, indicating that required sampling effort can be reached sooner than in Oribatida.

Results of this study provide an additional explanation why molecular species assignment often fails in Collembola and Oribatida. The more detailed analysis of the individual datasets at genus level showed that intra- and interspecific distances of taxa greatly overlapped, demonstrating the absence of a barcoding gap between species for all taxa, except for the parthenogenetic Oribatida genus Nothrus. The automated partitioning of datasets based on genetic distances (ASAP) suggested that each morphospecies (except most species within Nothrus and Platynothrus) consists of several genetic lineages, indicating the presence of putative or cryptic species. After assigning individuals according to genetic lineages, a barcoding gap between intra- and interspecific distances became apparent, but it still exceeded the standard threshold of 3%. Alternative partitions in the species assignment analyses also opted for smaller thresholds, but resulted in even more genetic lineages. From a conservative approach, the two sexual Oribatida species A. coleoptrata and S. magnus comprised 12 and 18 genetic lineages, respectively, with the relatively high barcoding threshold estimates of 6.9% (A. coleoptrata) and 15.0% (S. magnus). The Collembola species C. denticulata consisted of seven genetic lineages (barcoding gap of 13.8%) and F. quadrioculata of 24 genetic lineages (barcoding threshold of 8.0%). Notably, the parthenogenetic Oribatida species O. nova separated only into nine genetic lineages, P. peltifer into three and N. silvestris remained a single species which was consistent with morphological assignments.

In contrast to our hypothesis, the detection of a barcoding gap and thus species delimitation worked well for the parthenogenetic, but not for the sexual taxa. Species boundaries of Nothrus were clear and unequivocal. However, intra- and interspecific distances among Platynothrus overlapped, likely due the presence of three genetic lineages in P. peltifer. This is consistent with previous studies that identified seven genetic lineages in P. peltifer based on a transcontinental sampling, and demonstrated that lineages are consistent with species based on the 4x rule of parthenogenetic speciation (Heethoff et al., 2007; Birky & Barraclough, 2009; Birky et al., 2010).

Detection of deeply divergent genetic lineages in morphological consistent species is a common phenomenon and detection rate of cryptic species accelerated with the application of molecular identification tools (Bickford et al., 2007; Pfenninger & Schwenk, 2007; Skoracka et al., 2015; Struck et al., 2018). However, it remains important to consider these putative species carefully based on barcoding approaches, as delimitation is based only on a single genetic marker. The putative genetic lineages were highly congruent with nuclear haplotype diversity in C. denticulata, but not in F. quadrioculata and O. nova. Interestingly, genetic variance of nuclear and mitochondrial genes was opposite in the two latter species. In F. quadrioculata a single nuclear haplotype comprised many COI lineages, but in O. nova a single COI lineage comprised several nuclear haplotypes. This suggests that different selective forces might act on mitochondrial and nuclear genes in the two species. The higher mutation rate of mitochondrial compared to nuclear genes explains the higher diversity in COI in F. quadrioculata, indicating relatively recent divergence of lineages that had not yet been accompanied by variation in the nuclear gene. By contrast, in O. nova the mitochondrial gene shows relatively little variation, which likely is related to stronger purifying selection in parthenogenetic Oribatida (Brandt et al., 2017; Brandt et al., 2021). The two other parthenogenetic species (N. silvestris and P. peltifer) also show very little genetic variation, unfortunately, no additional genes were available for these taxa.

Our results demonstrate, that soil-living microarthropods comprise deeply divergent genetic lineages. Barcoding or metabarcoding studies based on single genes will therefore likely result in high numbers of unassigned reads or overestimate species numbers and consequently misrepresent species richness in communities. Potential species status should therefore be corroborated with an integrative taxonomic approach using multiple genetic markers and, if possible, re-examination of morphotypes (Schäffer, Kerschbaumer & Koblmüller, 2019; Lienhard & Krisper, 2021). However, morphological differences are often subtle, making traditional determination of soil microarthropods even more challenging. The nuclear 28S rDNA gene has been proposed as secondary barcoding marker for Oribatida (Lehmitz & Decker, 2017), but its applicability in a wider geographic range and different habitats has not been tested. Alternatively, metagenomic studies provide multiple genes per specimen which likely improves accuracy in specimen assignment. However, similar to metabarcoding based on single genes such as COI and/or 28S rDNA, successful application of metagenomics depends on representative reference databases. Notably, single reference genomes, or one/few barcodes per species will not cover intraspecific variation. Species with high intraspecific genetic variance would require “pan-barcodes” i.e., multiple barcodes from individuals that were sequenced across the range of a species to cover the extent of its intraspecific genetic variance.

The limited taxon sampling in this study demonstrates that even for the relatively intensive sequenced COI gene, databases do not provide taxonomic breadth for reliable species delimitation of Collembola and Oribatida. It is possible that species assignment will improve with a better reference database, but it is also important to understand the mechanisms that explain the barcoding gap, i.e., the substantial genetic divergence of COI sequences between closely related Collembola and Oribatida taxa. It is unknown if genetic variance is neutral or adaptive, or if mitonuclear or environmental interactions (Hill, 2020) generate the genetic structure in soil-living microarthropods. Fixation of neutral variance is one likely mechanism in the investigated taxa. The high numbers of haplotypes and nucleotide diversity suggest that COI is already highly saturated in these species. Many Collembola and Oribatida species are very abundant in local communities, suggesting high effective population sizes. This could enable the maintenance of neutral allelic variation and blur a barcoding gap in order to maintain the highly conserved protein sequence of COI. Repeated episodes of extreme population bottlenecks can also generate a barcoding gap between species. However, this is unlikely because high genetic variance in general argues against repeated population bottlenecks. However, the Oribatida species N. silvestris shows exceptionally low genetic variance compared to the other taxa, and consists of a single genetic lineage. It is possible that the low genetic variance resulted from a bottleneck this species experienced during Quaternary glaciations (<2.6 mya). Molecular divergence times among genetic lineages in the other species are several million years old, most date back to the Miocene (23-5 mya) and support the accumulation of neutral variance by genetic drift in Oribatida and founder events in Collembola (Rosenberger et al., 2013; von Saltzwedel, Scheu & Schaefer, 2016). Directional selection on mitochondrial genotypes and disrupted gene flow can lead to rapid divergence among populations. Collembola and in particular Oribatida are poor active dispersers due to their small body size, which reduces gene flow among populations and is a possible explanation for mitochondrial lineages corresponding with nuclear 28S rDNA haplotypes and sampling locations in C. denticulata (Porco et al., 2012b; von Saltzwedel, Scheu & Schaefer, 2016). However, reduced gene flow seems unlikely in F. quadrioculata due to the low genetic variance in the nuclear 28S rDNA gene compared to the highly variable mitochondrial COI gene. Genetic distances among lineages suggest maintenance of relatively ancient divergences, which argues against rapid divergence and disrupted gene flow. Further, this explanation does not apply for parthenogenetic species. Apparently, different mechanisms seem to account for the genetic variance in COI within species of Collembola and Oribatida. This is not surprising considering that the species in this study likely are separated by tens to hundreds of millions of years, each having its own evolutionary trajectory (Schaefer et al., 2010; Schaefer & Caruso, 2019; Leo et al., 2019; Katz, 2020; van Straalen, 2021).

This study showed that metabarcoding using the standard gene COI is problematic when investigating biodiversity of soil invertebrates. Advances in second- and third-generation sequencing technologies can significantly contribute to improve the reliability of barcodes for genetically diverse and potentially cryptic species. Proposed as an alternative to small barcoding fragments, low coverage shotgun sequencing and genome skimming offer increased species discrimination by covering entire organellar genomes and ribosomal sequences (Coissac et al., 2016). PacBio sequencing technology generates reads of approximately 3 kb with very low error rates. This enables sequencing of nearly full-length marker genes and their flanking regions, which improves taxonomic resolution and reduces spurious Operational Taxonomic Units (OTUs) (Tedersoo & Anslan, 2019). Notably, genomes of Collembola and Oribatida typically range between 350 and 500 Mb, enabling to obtain reasonable sequencing read depth at moderate prices. Further, wet-lab protocols for genome sequencing of small, non-model invertebrates have been developed (Collins et al., 2023) and the results underscore the importance of taking intragenomic variance into account in order to integrate genetic and morphological species boundaries. We propose that characterizing pan-genomes is crucial for identifying species in soil invertebrates (Tettelin et al., 2005). This approach will also contribute to develop informative barcoding genes (pan-barcodes) in soil invertebrates that lack a distinct barcoding gap. A pan-genome includes the complete set of genes shared by all individuals within a species and consists of conserved (core) and variable (accessory) gene regions (Golicz et al., 2020). The core genome covers all genes that are present in all individuals and the accessory genome includes the genomic regions that are variable among species. This variance is often due to ecological, geographical or reproductive boundaries (Reno et al., 2009). Accordingly, pan-genomes offer a holistic view of a species’ genome, allowing to identify both conserved and variable regions that are suitable for designing robust barcoding markers, in particular in taxonomically challenging organisms.

Conclusions

This study demonstrated that intra- and interspecific genetic divergences in the standard barcoding gene COI overlap in several species of Collembola and Oribatida. This is violating the assumption of a barcoding gap, which is a precondition for molecular species assignment and questions the applicability of the standard barcoding gene COI for soil-living microarthropods. Further, the presence of deeply divergent genetic lineages within morphologically consistent species emphasizes that (meta-)barcoding results solely based on a single genetic marker should be interpreted carefully. Based on COI, morphologically consistent species comprised numerous cryptic species. Without additional genetic and morphological data, the taxonomic status of these cryptic species is questionable. The assignment of genetic lineages to sequences in reference databases and application of flexible or species-specific thresholds could improve specimen assignment. However, the strong discrepancy between morphological conservativeness and genetic variance of many soil invertebrates calls for a more general approach. We are promoting to develop barcoding approaches with alternative sequencing technologies that generate more genetic data than metabarcoding, such as low-coverage shotgun sequencing of genomes (e.g., genome skimming and metagenomics) or long-read sequencing of marker genes using third generation sequencing technologies. Further, we advocate for the construction and analysis of pan-genomes to understand genetic species boundaries and to develop reliable barcoding markers that cover the whole range of genomic variance of species (pan-barcodes). Regardless of the approach taken, it is essential for reference databases to cover the intraspecific variability of a species throughout its geographic range.

Supplemental Information