DNA Metabarcoding Methods for the Study of Marine Benthic Meiofauna: A Review

Meiofaunal animals, roughly between 0.045 and 1 mm in size, are ubiquitous and ecologically important inhabitants of benthic marine ecosystems. Their high species richness and rapid response to environmental change make them promising targets for ecological and biomonitoring studies. However, diversity patterns of benthic marine meiofauna remain poorly known due to challenges in species identification using classical morphological methods. DNA metabarcoding is a powerful tool to overcome this limitation. Here, we review DNA metabarcoding approaches used in studies on marine meiobenthos with the aim of facilitating researchers to make informed decisions for the implementation of DNA metabarcoding in meiofaunal biodiversity monitoring. We found that the applied methods vary greatly between researchers and studies, and concluded that further explicit comparisons of protocols are needed to apply DNA metabarcoding as a standard tool for assessing benthic meiofaunal community composition. Key aspects that require additional consideration include: (1) comparability of sample pre-treatment methods; (2) integration of different primers and molecular markers for both the mitochondrial cytochrome c oxidase subunit I (COI) and the nuclear 18S rRNA genes to maximize taxon recovery; (3) precise and standardized description of sampling methods to allow for comparison and replication; and (4) evaluation and testing of bioinformatic pipelines to enhance comparability between studies. By enhancing comparability between the various approaches currently used for the different aspects of the analyses, DNA metabarcoding will improve the long-term integrative potential for surveying and biomonitoring marine benthic meiofauna.


INTRODUCTION
Marine sediments harbor high, but understudied biodiversity (Snelgrove, 1999;Armonies et al., 2018;Hughes et al., 2021). They are inhabited by a highly diverse and abundant community (Giere, 2009;Balsamo et al., 2010;Appeltans et al., 2012;Sinniger et al., 2016), a substantial fraction of which is part of the benthic meiofauna. Meiofauna usually refers to organisms that pass through a 1 mm sieve but are retained on a sieve of approximately 45 µm mesh size, although no fixed size definition exists (Higgins and Thiel, 1988;Brannock and Halanych, 2015). The appropriateness of mesh sizes for meiofauna collection has been recently investigated, and an alternative minimum mesh size of 20 µm has been proposed (Ptatscheck et al., 2020). Furthermore, some authors refer to organisms with a great anatomical and physiological modification to interstitial life, thus covering a vast taxonomic variety (viz., Zeppilli et al., 2018). Indeed, meiobenthic species are found in at least 24 out of 35 described animal phyla (Giere, 2009;Balsamo et al., 2012) and are characterized by a high abundance of up to 10 6 individuals/m 2 , consisting of (amongst others) nematodes, copepods, gastrotrichs, and platyhelminthes (Higgins and Thiel, 1988;Giere, 2009). Meiobenthic animals affect nutrient cycling, microbial community structure, hydrodynamics, and are a link in food webs between microbial organisms and macrobenthos (Montagna, 1984;Coull, 1999;Schratzberger and Ingels, 2018). Meiofaunal species commonly have fast reproduction and turnover rates, allowing for rapid, but also heterogeneous responses to environmental changes. This makes them potentially valuable as bioindicators to evaluate the quality of marine ecosystems and perform environmental impact assessments, which is important in times of high anthropogenic pressure and climate change (Semprucci et al., 2015;Zeppilli et al., 2015).
Traditionally, meiofaunal communities are studied by morphological identification using light microscopy (Coull, 1999;Giere, 2009). However, light microscopy techniques are time consuming and often impeded by difficult identification of species, due to few (documented) distinctive morphological features and expert taxonomists (Fontaneto et al., 2009;Curini-Galletti et al., 2012). Recent advances in DNA-based methods allow for the evaluation of biodiversity at comparatively low costs, high speed, and on a larger scale compared to classic morphological studies (Bik et al., 2012a;Brannock and Halanych, 2015). In particular DNA metabarcoding (Taberlet et al., 2012a), meaning high-throughput amplicon sequencing targeting specific gene regions like the mitochondrial cytochrome c oxidase I (COI) or the ribosomal nuclear small subunit (nSSU) 18S gene, termed molecular barcodes, is increasingly used to explore meiofaunal community composition in a wide range of ecosystems (Guardiola et al., 2015;Fonseca et al., 2017;Wangensteen et al., 2018;Weigand and Macher, 2018;Müller et al., 2019;Fais et al., 2020a). In addition, meiofaunal community metabarcoding data can be applied in ecological monitoring, as has been done to assess ecosystem quality and recovery (Chariton et al., 2015;Cordier et al., 2017), short-term direct anthropogenic effects (Laroche et al., 2016;DiBattista et al., 2020;Martínez et al., 2020;He et al., 2021) and climate change (Leasi et al., 2021).
Although the potential of metabarcoding of meiofaunal communities for biodiversity and ecological studies has been demonstrated, currently used workflows differ among studies. Moreover, a consensus on an ideal metabarcoding approach for monitoring and assessing benthic meiofauna has not yet emerged. Each step in metabarcoding workflows is a potential source of bias in recovering community diversity and composition, as shown for example for the type of sample used for DNA extraction (i.e., whole sediment cores versus meiofauna isolated from sediment, Brannock and Halanych, 2015;Klunder et al., 2019;Pansu et al., 2021), sampling depth and design (Montagna et al., 2017;Nascimento et al., 2018;Fais et al., 2020b), target genes (Haenel et al., 2017;Cordier et al., 2019;Atienza et al., 2020;Laroche et al., 2020), primer pairs (Cowart et al., 2015;Fais et al., 2020b), and bioinformatic processing (Brannock and Halanych, 2015;Leasi et al., 2018;Antich et al., 2021). This implies that careful consideration and evaluation of the applied methods depending on the research question is crucial, and accurate reporting of adopted methodologies is important for better comparability between studies.
This review focuses on DNA metabarcoding of marine benthic meiofauna, which complements recently published reviews focusing on various types of bulk samples (van der Loos and Nijland, 2021) and freshwater meiofauna (Schenk and Fontaneto, 2020) by evaluating methods used specifically for marine meiofauna in bulk (isolated) as well as the sediment samples. Moreover, since the review on use of metagenetic tools for meiofauna by Carugati et al. (2015), a vast number of experimental studies has been published addressing the performance of specific steps of the meiofauna metabarcoding workflow, leading to new insight. The aim of this review is to assist researchers in selecting a DNA metabarcoding workflow for the study of marine benthic meiofauna by summarizing the currently used methodologies in different steps of the metabarcoding procedure, identifying limitations, and pointing out aspects of the metabarcoding workflow that should be carefully evaluated in future studies.

METHODS
A search of scientific literature published between January 2010 and May 2021 was conducted within the Google Scholar database between March and May, 2021. Published articles were searched using terms that designate meiofaunal communities (i.e., "meiofauna, " "meiobenthos") in marine benthic environments (i.e., "beach, " "sediment") with or without the commonly used marker genes (i.e., "COI" and "18S") in different combinations with the terms "marine" and "metabarcoding". Search terms are disclosed in Supplementary Table 1.
Each paper was initially screened to confirm that the study included DNA metabarcoding of benthic meiofauna communities in marine environments (including estuarine and coastal areas). Publications on specific meiofaunal taxa (e.g., Nematoda) were excluded, as the applied methods (e.g., specimen extraction and taxa specific primer pairs) mostly cannot be applied to the study of whole meiofaunal communities. In studies on marine sediment environmental DNA (eDNA) or metabarcoding of complete benthic communities, the parts that focused on non-meiofaunal taxa were not further analyzed.
The information retrieved from each selected publication included: (1) the geographic area of the study; (2) the type of sampling sites; (3) the amount of sediment collected; (4) the number of replicates per sampling point; (5) the sampling instrument; (6) the depth of the layer of sampled sediment; (7) number of sliced layers; (8) the smallest sieve mesh size; (9) sample conservation method prior to DNA extraction; (10) the type of sample for DNA extraction; (11) the DNA extraction protocol, (12) amount of sediment for DNA extraction; (13) implementation of extraction controls; (14) number of extraction replicates; (14) the amplified marker gene(s) and gene region(s); (15)  The geographic area was identified at the national level. The type of sampling site was identified according to four categories: estuary, intertidal, subtidal up to 200 m depth, or at depths of over 200 m. The amount of sediment collected per sampling site was reported in various units in different studies (volume, mass, and surface area). The amount of sediment was converted to mL by squaring the inner radius of the sediment core in cm and multiplying with π * height, if reported. Amounts reported in other units (e.g., mass, surface area) were not converted as this could lead to imprecise estimates. The depth of the layer of sampled sediment was reported in cm starting from the surface sediment (0 cm depth). The type of the sample for DNA extraction was identified according to three categories: from sediment samples, from isolated meiofauna, or from extracellular DNA. Extraction of extracellular DNA was distinguished from extraction of sediment (total) DNA by the use of a saturated phosphate buffer protocol optimized for the recovery of extracellular DNA while avoiding organismal DNA (Taberlet et al., 2012b;Guardiola et al., 2015). The sieve mesh size was reported only for studies in which meiofauna was extracted from the sediment. As MoBio Laboratories was acquired by QIAGEN in 2015, data from identical DNA extraction kits that were first produced by MoBio was combined with QIAGEN data (e.g., MoBio PowerSoil kit and QIAGEN DNeasy PowerSoil kit). The bioinformatic unit was reported as the bins into which sequences are clustered based on their similarity (e.g., Operational Taxonomic Unit, Amplicon Sequence Variant). The bioinformatic processing was categorized according to the similarity threshold or pipeline used for delineation of the respective bioinformatic unit. We report sequencing depth as the total number of raw and quality filtered reads. All other variables were noted in the review as reported in the original paper, unless indicated otherwise. When authors did not mention certain methods, they were noted as "not reported". Graphs were created using the ggplot2 package v3.3.3 (Wickham, 2016) in R v4.0.3 (R Core Team, 2020), Rstudio v1.3.1093 (RStudio .

RESULTS
A total of 67 papers matched the criteria and were included in this review (see Supplementary Table 2 for a list of publications and metadata). The sum of studies (n) for different steps of the metabarcoding workflow sometimes exceeds the total number of analyzed studies, as some studies applied multiple methods. The first study on DNA metabarcoding of benthic meiofauna was published in 2010 (Creer et al., 2010), and the number of meiofaunal community metabarcoding studies increased over the last decade (Supplementary Figure 1).
The reviewed papers covered samples taken from estuarine (n = 12), intertidal (n = 12), subtidal (n = 26), or deep (n = 12) sediments, or a combination of the previous (n = 5, Supplementary Figure 2). Samples were usually collected with a grab (n = 11), corer (n = 23), box corer (n = 7), multicorer (n = 10), or directly into tubes or jars (n = 9). Directly comparing studies regarding the amount of sediment collected was often impossible due to differences in reported units. However, it is apparent that the amount of sediment that was collected per replicate varied greatly amongst studies. Differences of over 1,000-fold were found, ranging from 0.035 to 5 kg and 0.003 to 4.5 L per replicate. Most studies implement 2-5 sampling replicates (n = 49), whereas more than five replicates (n = 5) or one replicate (n = 9) were used less frequently. Approximately a third of studies analyzed collected sediment from the upper sediment layer to a depth of 2 cm or less (n = 20), n = 18 studies collected sediment up to a depth of 3-5 cm, n = 17 to a depth of 10 cm, and sediment cores deeper than 10 cm were used in n = 5 studies. Collected sediment was mostly examined as a whole (n = 64) rather than separately for two (n = 1) or three (n = 3) different depth layers. After sample collection, samples were conserved predominantly frozen (n = 40), in ethanol (n = 14, varying concentrations and temperatures) or in DESS (n = 8).
DNA was most often extracted directly from sediment samples (n = 40), followed by DNA extracted from meiofaunal organisms separated from the sediment (n = 29), or from isolated extracellular DNA (n = 4). n = 27 studies used between 8 and 10.5 g or 10 mL of sediment or meiofauna isolate for DNA extraction, whereas extractions from relatively small amounts of material (0.2-0.35 g) are common as well (n = 11). Amounts were not clearly reported in n = 17 studies. Isolation of meiofauna from sediment was mostly done by sieving (n = 5), or through decantation over a sieve using Ludox solution (n = 10, separation based on density) or the anesthetic MgCl 2 (n = 5, Supplementary Figure 3). Sieve mesh widths ranged from 20 to 500 µm, with 45 µm being used most frequently (n = 15, Supplementary Figure 4).
DNA from sediment samples was primarily extracted using the commercial QIAGEN (previously MoBio) PowerSoil (n = 17) and PowerMax (n = 23) kits (Figure 1). DNA from isolated meiofaunal samples was mostly extracted with these kits as well, with the PowerSoil used more frequently (n = 13) than the PowerMax Soil kit (n = 7). Kits or protocols not specifically designed for sediment samples were used in n = 9 studies (e.g., QIAGEN QIAamp DNA Blood Maxi kit, n = 5). Sediment extracellular DNA was extracted with a saturated phosphate buffer protocol (Taberlet et al., 2012b), in combination with either a commercial DNA extraction kit (n = 3) or chloroform and isoamyl alcohol (IAA) protocol (n = 1). n = 12 studies report on the number of extraction replicates used, varying from one to ten extraction replicates. Negative DNA extraction controls (field controls or extraction blanks) were reported to be included in n = 16 studies, while positive controls (including DNA spiking or mock samples) were reported in n = 8 studies.
Primers targeting the COI gene (n = 21) all amplify parts of the 5 end of the gene, known as the Folmer region (Folmer et al., 1994), which is the standard molecular barcode for many animal taxa. A total of seven different primer combinations were used for the amplification of the COI fragment. In most cases, the mlCOIintF forward primer was combined with the reverse primer jgHCO2198 (n = 7, Leray et al., 2013), which amplifies the 313bp "mini-barcode region" in the 3 of the Folmer fragment. The mlCOIintF primer in combination with the reverse primer LoboR1 (Lobo et al., 2013) was used less frequently (n = 5). Almost two thirds of studies used more than one replicate PCR (n = 40). Of these, PCRs were mostly performed in duplicates (n = 12) or triplicates (n = 25) to as many as eight replicates (n = 2). PCR negative controls (water only and no-template) were included in n = 34 studies. FIGURE 2 | A total of 18 and 7 primer combinations were used that amplify the different variable regions of the 18S rRNA gene or partial segments of the 5 end of the mitochondrial COI, respectively. Bars of primer pairs with a matching color target the identical molecular marker as indicated in the right margin of the figure.
Sequencing was mostly conducted on the Illumina MiSeq sequencing platform (n = 39), followed by Roche 454 pyrosequencing (n = 17; Supplementary Figure 5). Sequencing depth was found to be reported most commonly as the total number of raw reads obtained from sequencing, and the number of reads retained after quality filtering. As the total number of sequenced samples and replicates was not always reported, inferring sequencing depth per sample was often difficult. Therefore, we report the total number of reads per sample to illustrate the wide variation in sequencing depth. We found a range of approximately 125,000-61,000,000 total raw reads per marker, and 50,000-45,000,000 total quality filtered reads per marker. Subsequently, a variety of bioinformatic clustering methods and denoising tools were used to assign sequences into Operational Taxonomic Units (OTUs) or Amplicon Sequence Variants (ASVs; Supplementary Figure 6). OTU clustering (n = 52) was mostly performed using a cut-off similarity threshold (n = 45), with 97% genetic identity being the most commonly used threshold (n = 24), followed by 96% similarity (n = 7). Denoising and ASV formation was done using the DADA2 pipeline (n = 12, Callahan et al., 2016). In other studies, ASVs were subsequently clustered into OTUs using varying methods (n = 6) with a similarity threshold or with SWARM v2 (Mahé et al., 2015). Sequences were subsequently assigned to meiofaunal taxa using, for example, NCBI GenBank or nucleotide database (n = 22, Sayers et al., 2021), SILVA (n = 40, Pruesse et al., 2007), MIDORI (n = 6; Machida et al., 2017), or BOLD (n = 2, Ratnasingham and Hebert, 2007).

DISCUSSION
We found a diverse array of methods in the metabarcoding workflow of meiobenthic community studies analyzed. Variance in metabarcoding procedures amongst studies is shared with other target communities, such as benthic macroinvertebrates . Our results specifically identified discrepancies in pre-treatment of samples and DNA extraction techniques, selected marker genes and regions, preferred primer pairs, sequencing procedures, and bioinformatic processing of sequences. In the following sections, we discuss how bias may be potentially introduced within each of these steps, which conditions would favor the choice for a certain method over another, and indicate which aspects need careful consideration in future metabarcoding studies.

Sample Collection and Pre-treatment
Our results showed that the amount of sampled sediment varied greatly amongst studies and was often not reported clearly. The amount of sampled sediment may be important especially when DNA is extracted from isolated meiofaunal specimens or from a subsample of homogenized sediment and should thus be reported consistently in addition to the amount of sediment used for DNA extraction. In 92% of studies in which the sampling depth of the sediment layer was reported, the upper 10 cm of sediment or less was sampled, corresponding to the most biologically active part of the sediment in which meiofauna generally tends to concentrate (Thiel, 1983;Simpson and Batley, 2016).
Sample preservation is a critical step to avoid DNA degradation and was consistently reported in all but three studies. Our survey highlights that freezing or freeze-drying is the predominant method used, without fixatives such as ethanol or DESS. However, when a fixative was used this mostly comprised isolated meiofauna samples rather than sediment. Consistently, van der Loos and Nijland (2021) report that freezing is often used for sediment samples, and recommend the use of DESS over ethanol for bulk samples as it preserves both DNA and morphology with higher quality.
DNA extraction was most often performed directly from sediment samples compared to specimens isolated from sediment prior to DNA extraction. Compared to isolation of meiofauna, using sediment samples may be easier, it limits processing time, it avoids contamination, and it does not introduce any bias due to different methods of specimen isolation Fais et al., 2020b). Our survey highlights that sediment samples are dominant in studies targeting environments hard to sample, such as the deep sea (e.g., Sinniger et al., 2016;Atienza et al., 2020;Kitahashi et al., 2020), and muddy and fine sediments difficult to elutriate (i.e., separating meiofauna from sediment by specific gravity using a liquid stream, Somerfield et al., 2005;viz., Lanzén et al., 2017;Faria et al., 2018;Salonen et al., 2019;Brandt et al., 2020). In addition, studies including a wider size range of benthic organisms such as micro-and macrofauna assessed through metabarcoding of sediment eDNA (e.g., Xie et al., 2017Xie et al., , 2018Fais et al., 2020a,b;Klunder et al., 2020a,b) might have artificially raised the number of studies in which DNA was directly extracted from sediment samples. This method may lead to partially misleading results on the species composition present, as DNA in sediment can be recent or older extracellular DNA (Thomsen and Willerslev, 2015). This is especially relevant in deep sea environments, where extracellular DNA can constitute for 50-90% of the DNA, and ancient DNA can be very well conserved because of adsorption onto the sediment matrix (Torti et al., 2015). On the contrary, isolation of meiofaunal specimens from the sediment allows researchers to process a larger volume of sediment per sample, thereby possibly reducing sample heterogeneity and increasing recovery of rare species Brannock and Halanych, 2015). Sieving reduces the dominance of large taxa (Cowart et al., 2015;Fonseca et al., 2018;He et al., 2020). Furthermore, elutriating or sieving samples was reported to result in a relatively high abundance of metazoan taxa while reducing the amplification of non-target bacterial or non-benthic taxa (Brannock and Halanych, 2015;He et al., 2020). When isolating meiofauna from the sediment, careful consideration is needed as soft-bodied organisms are easily damaged, whereas relatively heavy hard-bodied taxa such as mollusks may be retained in the sediment during elutriation, possibly leading to an underrepresentation of certain taxonomic groups (van der Loos and Nijland, 2021).
Recently, Castro et al. (2021) compared the recovered community composition from DNA extracted directly from sediment cores with two meiofauna extraction methods commonly used in metabarcoding studies, namely (1) separation of meiofaunal specimens from sediment by flotation and centrifugation with Ludox and (2) separation of meiofaunal specimens from sediment by decantation of specimens with MgCl 2 . The community composition recovered through direct DNA extraction from sediment differed strongly from that obtained with the isolation techniques (Castro et al., 2021). Isolation of specimens with Ludox or MgCl 2 led to more taxa and meiofaunal species compared to DNA extraction directly from sediment (Castro et al., 2021). Differences were also found amongst isolation techniques, as some genera belonging to the phyla of platyhelminthes, annelids, mollusks, xenacoelomorpha, gastrotrichs, and nematodes were found in a lower abundance in samples obtained through MgCl 2 decantation relative to Ludox flotation (Castro et al., 2021).
Hence, as sample type significantly affects both the OTU richness and taxonomic composition of meiofaunal communities (Brannock and Halanych, 2015;Haenel et al., 2017;Nascimento et al., 2018;Castro et al., 2021), a greater consideration should be given to the choice of sample and DNA extraction workflow. Indeed, any extraction method will lead to a bias toward taxonomic groups that are more easily isolated with the chosen technique. Therefore, if the goal is to maximize recovery of taxa, a combination of isolation methods is the best option (Castro et al., 2021).

DNA Extraction and PCR Procedures
We found that when DNA was extracted directly from sediment, researchers used kits specifically designed for soil samples, while these kits were used in only two-thirds of studies that extract DNA from meiofauna isolated from the sediment (Figure 1). Humic substances often present in sediments can be retained during DNA extraction and inhibit enzymes that amplify DNA during PCR, which can lead to false negatives (Matheson et al., 2010). Soil DNA extraction kits specifically remove these substances, thereby alleviating PCR biases caused by enzyme inhibition. Nonetheless, in earlier studies, clean isolated meiofauna samples and good amplification results were obtained using non-soil kits (e.g., Creer et al., 2010;Fonseca et al., 2010). More recently, it was recommended to use a soil-specific kit in samples that may contain remains of sediment to obtain high-purity DNA (van der Loos and Nijland, 2021), which also works well for marine and coastal sediments (Lekang et al., 2015;Taberlet et al., 2018;Fais et al., 2020b). The two most used kits in our analysis (i.e., MoBio/QIAGEN PowerSoil and PowerMax kit) differ in the amount of sediment that can be processed (a maximum of 0.25 and 10 g, respectively). Indeed, we found that in the analyzed studies DNA is mostly extracted from 5 to 10 g or 0.2 to 0.35 g, corresponding to the processing maxima of these kits. As some studies have shown that the number of meiofaunal taxa and OTUs increase with sample size where no elutriation steps are performed (Nascimento et al., 2018;Brandt et al., 2020;Fais et al., 2020b), 10 g of sediment is recommended for eDNA extraction. However, we point out that the QIAGEN PowerMax is costly (about 25€/USD per extraction), thereby limiting the number of samples that can be processed in many projects, as reflected in the number of extraction replicates found in this study. None of the studies that used the PowerMax kit used multiple extraction replicates, while these were reported using the PowerSoil kit (QIAGEN). Ultimately, several replicate DNA extractions and a lysis step are strongly suggested to improve OTU and taxa richness in studies on marine meiobenthos through DNA metabarcoding approaches (Lanzén et al., 2017). Additionally, it is recommended to perform at least three replicate PCR reactions to correct for the effects of PCR stochasticity and to potentially increase the number of OTUs and species detected (Leray and Knowlton, 2015;Bourlat et al., 2016;Alberdi et al., 2018), which was performed in less than half of the studies and not reported in 36% of studies analyzed here. Again, whether PCR replicates improve the results and outweigh the costs depends on the research goal (van der Loos and Nijland, 2021).
Even though negative extraction and PCR controls are essential control measures within metabarcoding, their use was reported in less than half of the studies analyzed here. Whether these controls were truly not implemented in the analyses, or this result only reflects a gap in reporting remains unknown, but either way this indicates an important point of improvement (van der Loos and Nijland, 2021). Field extraction controls (i.e., opening sample storing material at the sampling site without sample collection), no-template DNA extraction controls, notemplate or water-only PCR controls all aim to identify potential contamination introduced in any step of the workflow and should therefore be processed and sequenced along the field samples and controlled for during bioinformatic quality filtering steps (e.g., Brandt et al., 2021a).

Choice of Marker Gene and Region
We found that the nuclear 18S rRNA molecular marker was the most commonly used, especially when only one marker region was amplified (n = 44 for 18S and n = 1 for COI). Most studies did not mention the reason behind marker region choice. Amongst the reasons for the choice of marker, sequencing length limitations of the Illumina platform was sometimes mentioned as a reason for selecting the 18S rRNA V9 region (Brannock et al., 2014(Brannock et al., , 2018Brannock and Halanych, 2015). In other studies, 18S V4 was chosen over COI because this region is shorter and easier to PCR amplify and sequence using Illumina, and amplifies only eukaryotic sequences (Lejzerowicz et al., 2015;He et al., 2021). Furthermore, poor amplification results for marine nematodes with the COI Folmer primers in early barcoding studies (Derycke et al., 2010) have stimulated the use of 18S rRNA gene regions as genetic markers. The more preserved nature of 18S rRNA primer binding sites compared to COI allows for the amplification of a broader range of taxa (Schenk and Fontaneto, 2020). However, the relatively high conservation of 18S priming sites comes at the cost of resolution, often lacking the ability to discriminate at the species level (Tang et al., 2012).
Despite the extensive use of the 18S rRNA marker to study meiofauna, the mitochondrial COI gene is becoming more popular for meiofauna metabarcoding studies (28% of the studies analyzed in this review). Compared to 18S rRNA, the main advantage of COI is its high taxonomic resolution, often allowing for identification of species and even intraspecific variability (Tang et al., 2012;Cowart et al., 2015;Leray and Knowlton, 2016;Turon et al., 2020). Notwithstanding, meiofaunal species are often underrepresented in molecular reference databases (Sinniger et al., 2016;Curry et al., 2018;Wangensteen et al., 2018), leading to high percentages of uncharacterized OTUs. For example, 70.9% of the eukaryotic clusters remained unassignable at the phylum level in a deep-sea study by Atienza et al. (2020). However, it is important to consider that a high number of unassigned reads in COI datasets can be a result of non-target amplification of bacterial and other non-target species (Weigand and Macher, 2018). Therefore, both commonly used markers have their pitfalls and the choice of the most appropriate one depends on methodology, study system and research objective, as is true for all steps in the metabarcoding workflow. Using multiple molecular markers within a single study is advisable to maximize the recovery of meiofaunal taxa (Fais et al., 2020b and references therein), which is now implemented by less than a third of the studies analyzed. Furthermore, efforts should be undertaken to barcode individual species and place molecular references in barcode reference libraries to improve our ability to taxonomically assign metabarcoding data (Wangensteen et al., 2018;Weigand et al., 2019).

Selection of Primer Pairs
Optimal primer pairs for amplification of whole meiofaunal communities ideally amplify a broad range of meiofaunal animals while minimizing off-target amplification. Ideally, extensive reference databases containing sequences for this marker should be available to allow for high taxonomic resolution. Our data showed that a total of 22 different primer combinations were used in the studies analyzed in this review. This high number can be partly explained by studies that test multiple primer pairs, but it also reflects that a consensus on the optimal primers for amplification of marine meiofaunal communities has not been reached. Primer choice is important as primer bias due to mismatches affects the relative abundance of amplicons and species recovery (Elbrecht and Leese, 2015). Thus, conserved sequences facilitate universal primer design. As such, the SSU_F04 forward and SSU_R22 reverse primer set (Blaxter et al., 1998;Fonseca et al., 2010) targeting the 18S rDNA V1-V2 region is highly conserved and amplifies a broad range of meiofaunal organisms (Creer et al., 2010). In meiofaunal taxa, the V1-V2 region amplified by these primers is the most variable nSSU region and its length of 450 bp favors amplification of DNA from living organisms (Lallias et al., 2015), benefiting recovery of present meiofaunal communities. Further modifications, such as shortening and including an ambiguity (R) into the SSU_R22 reverse primer to accommodate the inclusion of more metazoans (called SSU_R22mod, Sinniger et al., 2016), can be used to target specific groups of taxa within meiofauna.
After 18S rRNA V1-V2, the V4 and V9 were commonly used to amplify eukaryotes. The 18S rRNA V4 region is the most variable region in eukaryotes (Nickrent and Sargent, 1991). Indeed, the 18S V4 primer pairs are designed to amplify a broad range of eukaryotic taxa [F566 and R1200 by Hadziavdic et al. (2014); Uni18S and Uni18R by Zhan et al. (2013); E572F and 897R by Comeau et al. (2011) andHugerth et al. (2014)] or eukaryotic microbes [TAReuk454FWD1 and TAReukREV3 by Stoeck et al. (2010)] instead of meiofauna taxa specifically. The same goes for the eukaryotic primer set 1380F and 1510R targeting the V9 region (Amaral-Zettler et al., 2009). The short (100-110 bp) regions that are amplified using the 18S_allshorts primer set targeting 18S rRNA V7 are suitable for amplifying extracellular DNA due to its length, however limiting meiofaunal identification at species level (Guardiola et al., 2015(Guardiola et al., , 2016. The 658-bp fragment amplified by the original primer pair for the mitochondrial COI gene (LCO1490 and HCO2198, Folmer et al., 1994) is too long for commonly used Illumina sequencing platforms. The shorter mlCOIintF and dgHCO2198 (Meyer, 2003;Leray et al., 2013) and mlCOIintF and jgHCO2198 (Geller et al., 2013;Leray et al., 2013) are most commonly used for amplification of marine meiobenthos. However, the degenerate LoboR1 primer (Lobo et al., 2013) is a promising reverse primer in combination with the mlCOIintF forward primer, as it uncovered a higher diversity compared to the degenerate reverse primer dgHCO2198 (Haenel et al., 2017). This combination was used and recommended in several meiofaunal community studies (Haenel et al., 2017;Fais et al., 2020a,b;Castro et al., 2021). Additionally, inosine bases present in the widely used jgHCO2198 reverse primer may impair the performance of highfidelity polymerases commonly used for PCR amplification in metabarcoding studies (Jungbluth et al., 2020). For analyses of meiobenthic communities extracted from the sediment, using highly degenerate primers that amplify the maximum number of taxa, such as the Leray and Lobo primers (Meyer, 2003;Leray et al., 2013;Lobo et al., 2013), may be beneficial. On the other hand, these primers can be suboptimal for amplification and sequencing of DNA extracted directly from sediment, as the primers can amplify a wide range of non-target microbial taxa abundantly present in the sediment preventing detection of target taxa (Alberdi et al., 2018;Weigand and Macher, 2018). Alternative methods such as elutriation or sieving, as suggested before, could reduce non-target amplification. Alternatively, less degenerate, more specific primers can be beneficial, at the cost of potentially losing some target taxa.
Consequently, no optimal primer pair is currently available for marine benthic meiofauna. To increase comparability between studies, researchers should aim to narrow down the amount of different primer combinations that are used in meiofaunal community studies, while considering the best primer pairs for their research objectives.

Sequencing
Primer pair selection eventually affects the choice of sequencing platform, and vice versa. Illumina MiSeq is currently the main sequencing platform used in meiofaunal community DNA metabarcoding studies and can produce 2 × 300 bp long reads, rendering the use of primer pairs that generate longer amplicons unfeasible. The 454 pyrosequencing platform was identified as the main sequencing platform in Carugati et al. (2015). In our review, we found that this platform was mainly used in older papers, as Roche discontinued 454 platform production in 2013 and support by 2016. It was recently shown that Illumina's NovaSeq 6000 outcompetes the MiSeq platform in a seawater eDNA metabarcoding study, where NovaSeq detected a higher diversity at the same sequencing depth. Moreover, where the number of new sequence variants detected by the MiSeq platform leveled off at approximately one million reads, data from NovaSeq did not show such a plateau (Singer et al., 2019). At ten times the costs of MiSeq sequencing, however, the NovaSeq 6000 platform may be inaccessible for most studies (Singer et al., 2019) and was indeed not used in the studies analyzed in this review. However, costs might decrease and the platform, which can produce 2 × 250 bp long reads, might become a cost efficient option for future metabarcoding studies. In contrast, the popular Illumina MiSeq platform has relatively low costs per sequenced base pair. Sequencing techniques such as PacBio and Nanopore allow sequencing of longer amplicons (Callahan et al., 2019) and have been successfully used for metabarcoding of environmental samples (Heeger et al., 2018;Karst et al., 2018;Davidov et al., 2020;Okazaki et al., 2021). Even though trade-offs regarding the higher error rate and read coverage compared to Illumina short-read sequencing exist, longer reads comprising several genes may provide better taxonomic resolution due to the high informational content per read (Weirather et al., 2017;Baloglu et al., 2021). To our knowledge, these sequencing platforms have not been used in marine meiofaunal metabarcoding studies and might thus be valuable options for future studies (Hebert et al., 2018).
Diversity detection increases with sequencing depth in metabarcoding studies, as shown for different target taxa (e.g., Smith and Peay, 2014;Lanzén et al., 2017;Alberdi et al., 2018). Sequencing depth increases the number of unique OTUs per PCR replicate, and thereby decreases similarity between PCR replicates (Alberdi et al., 2018). If a study aims to identify the whole meiofaunal diversity, a high sequencing depth is crucial. However, for general monitoring purposes deep sequencing depth may not be required. Nonetheless, a sufficient sequencing depth should be considered fitting the purpose of the study. Singer et al. (2019) inferred sequencing depth of 10 metabarcoding studies using DNA from various sources and found that an average sequencing depth of 55,000-60,000 reads per sample is common. Indeed, they also found that sequencing depth is often not reported clearly. Similar to the use of highly degenerate primers, detection of more species comes at the cost of an increasing number of artifactual or non-target sequences and should thus be considered during bioinformatic processing (Alberdi et al., 2018;Weigand and Macher, 2018).

Bioinformatic Approach
It has been shown by multiple studies that bioinformatic techniques affect recovered diversity (Brannock and Halanych, 2015;Leasi et al., 2018;Antich et al., 2021). We found a total of 23 different bioinformatic approaches used in studies on marine benthic meiofauna, which can potentially render comparison of results between studies unfeasible. The clustering methods (i.e., combining sequences into representative units such as OTUs or MOTUs) analyzed in this study are mostly based on either Bayesian methods (CROP; Hao et al., 2011), single linkage clustering (SWARM; Mahé et al., 2015), or similarity thresholds. The latter was most commonly applied in the reviewed studies, often with a cut-off value of 97% similarity. Denoising approaches without clustering (i.e., rectifying sequencing errors by merging sequences containing an error with the true sequence) as done in DADA2 (Callahan et al., 2016) or UNOISE (Edgar, 2016) pipelines, yields ASVs or zero-radius OTUs (ZOTUs), which allow studying intraspecific diversity. Similarly, but not widely implemented in meiofauna research, the Deblur pipeline (Amir et al., 2017) obtains single-nucleotide resolution and produces comparable weighted results to DADA2 and UNOISE2 (Nearing et al., 2018) and could be valuable for identification at (sub)species level.
Currently, the debate whether OTUs or ASVs are more appropriate for analyzing metabarcoding data is ongoing. It is argued that ASVs, due to their higher resolution and reusability across studies, should be the standard reporting unit (Callahan et al., 2017), but other researchers argue that clustering is crucial in more variable markers such as COI to conform to the biological species concept (Antich et al., 2021;Brandt et al., 2021b). Furthermore, denoising and clustering programs that were originally developed for rRNA (e.g., CROP, SWARM, DADA2, UNOISE) are now commonly used for the mitochondrial COI gene (Antich et al., 2021). Using these programs on sequences from variable markers such as COI demands for a critical consideration of parameter settings as done for SWARM2 Brandt et al., 2021b) and CROP (Leray et al., 2013;Leray and Knowlton, 2015;Wangensteen et al., 2018). However, not all studies provide the used parameter settings, implying that the default settings were used without further consideration (Antich et al., 2021).
It has been pointed out in many studies that taxonomic assignment of recovered sequences is hampered by significant gaps in reference databases (Sinniger et al., 2016;Curry et al., 2018;Wangensteen et al., 2018;Weigand and Macher, 2018), leading to difficulties especially for the interpretation of metabarcoding data at lower taxonomic levels. The COI gene tends to be more variable than rRNA genes like 18S, meaning that for reliable identification on a low taxonomic level, representatives of the species or at least closely related species have to be present in a reference database (Leray and Knowlton, 2016). Moreover, databases are often limited to few genetic markers (e.g., SILVA, which contains rRNA sequences, BOLD, which is focused on COI), contain many sequences that are not openly available and may contain erroneous or wrongly assigned sequences (e.g., Radulovici et al., 2021). Combining multiple reference databases could therefore improve the number and accuracy of OTUs assigned within a study (Macher et al., 2017). However, this was only implemented by a limited number of studies analyzed here. Almost two-third of studies use the ribosomal RNA sequence database SILVA (Pruesse et al., 2007) for taxonomic assignments, corresponding to 18S being the most used marker in meiobenthic metabarcoding studies. Surprisingly, BOLD was only used in n = 2 studies. Studies targeting COI tend to rely more on GenBank, as some taxonomic groups, including annelids and mollusks, are better represented in GenBank compared to BOLD (Weigand et al., 2019).
Overall, implementation of a mock sample with known species composition could confirm the reliability and accuracy of the metabarcoding workflow, including bioinformatic processing, and can be used not only to verify the detection of mock species but also control for overestimation of OTUs. Of the studies reviewed, only n = 6 studies (e.g., Chariton et al., 2015;Klunder et al., 2019;Brandt et al., 2020) included a mock community, and routinely sequencing such a control sample might be beneficial for future studies.

CONCLUSION
We found that a consensus on optimal methods for metabarcoding of marine meiofaunal sediment communities has not yet emerged, making comparison of study results difficult or impossible. We identified several key aspects that should receive attention in future studies on marine meiobenthos: (1) comparison of sampling methods, i.e., direct extraction of DNA from sediment versus separation of specimens from the sediment, followed by DNA extraction; (2) comparison of different primers and multiple molecular marker regions to maximize taxon recovery; (3) all steps of the sampling and laboratory procedures should be accurately described to allow for replication and better comparability of experiments; and (4) bioinformatic approaches for processing data and assigning species names to sequences should be thoroughly compared and evaluated. Considering the vast number of improvements that have been made over the last decade, we believe metabarcoding will further develop into an accurate and efficient standard tool to monitor marine meiobenthos in the coming decade.

AUTHOR CONTRIBUTIONS
RG designed and performed the research, analyzed data, and drafted the work. J-NM designed the research, interpreted data, and drafted and revised the work. MF, DF, WR, SC, and FC drafted and revised the work. All authors contributed to the article and approved the submitted version.

FUNDING
This work was funded by a BEN (Biodiversity-Ecology-Nature) grant (Number T0206/37197/2021/kg) of the Bauer-Hollmann foundation to J-NM.

ACKNOWLEDGMENTS
We thank the Naturalis Marine Biodiversity Group for support and helpful comments and two reviewers for their suggestions.