Biodiversity between sand grains: Meiofauna composition across southern and western Sweden assessed by metabarcoding

Abstract The meiofauna is an important part of the marine ecosystem, but its composition and distribution patterns are relatively unexplored. Here we assessed the biodiversity and community structure of meiofauna from five locations on the Swedish western and southern coasts using a high-throughput DNA sequencing (metabarcoding) approach. The mitochondrial cytochrome oxidase 1 (COI) mini-barcode and nuclear 18S small ribosomal subunit (18S) V1-V2 region were amplified and sequenced using Illumina MiSeq technology. Our analyses revealed a higher number of species than previously found in other areas: thirteen samples comprising 6.5 dm3 sediment revealed 708 COI and 1,639 18S metazoan OTUs. Across all sites, the majority of the metazoan biodiversity was assigned to Arthropoda, Nematoda and Platyhelminthes. Alpha and beta diversity measurements showed that community composition differed significantly amongst sites. OTUs initially assigned to Acoela, Gastrotricha and the two Platyhelminthes sub-groups Macrostomorpha and Rhabdocoela were further investigated and assigned to species using a phylogeny-based taxonomy approach. Our results demonstrate that there is great potential for discovery of new meiofauna species even in some of the most extensively studied locations.


Introduction
Meiofauna (i.e. animals that pass through a 1 mm mesh, but are retained on a 45 µm mesh; Higgins and Thiel 1988) have been relatively unexplored (Eisenhauer et al. 2019), despite the fact that they are an important component of benthic habitats and play key roles in their environments (Aller and Aller 1992, Schratzberger and Ingels 2018, Snelgrove and Smith 2002. Assessments of meiobenthic diversity have often been performed manually using traditional morphological approaches, whereby researchers pick individual animals from a sample and attempt to identify the animal to species level (Giere 2009, Higgins andThiel 1988). Such explorations are time-consuming and require high levels of taxonomic expertise and are therefore, by-necessity, focused on a limited number of taxonomic groups. Moreover, there is potential for great amounts of missed diversity, as animals that are particularly small, delicate, or simply unusual may be missed or ignored, while the presence of asexual or juvenile animals could make identification impossible.
Nevertheless, there is a long history of studies on meiofaunal taxa in the North Sea and Baltic waters surrounding Sweden. For instance, Westblad (Westblad 1940, Westblad 1942, Westblad 1945, Westblad 1946, Westblad 1948) published a series of papers describing Acoela and Platyhelminthes from the Gullmarn Fjord and surrounding Swedish and Norwegian waters, while Swedish marine Rhabdocoela (Platyhelminthes) were assessed on several occasions by, for example, Karling (Karling 1956a, Karling 1956b, Karling 1963, Karling 1974), Luther (Luther 1962, Luther 1963) and Westblad (Westblad 1954). For a long time, the Kristineberg Marine laboratory was especially the focal point of many studies of meiofaunal diversity (e.g. Boaden 1960, Elofson 1944, Jägersten 1952, Odhner 1937 and many others) and broad meiofaunal surveys have more recently occurred at the the Sven Lovén Centre for Marine Sciences on the island of Tjärnö (Willems et al. 2009, Curini-Galletti et al. 2012 . As a consequence, the Swedish meiofauna is comparatively well known.
Further in 2002, the Swedish Species Information Centre (ArtDatabanken) established a 20-year project to improve the taxonomical knowledge of all Swedish multicellular organisms and new opportunities were provided for modern day surveys of many different meiofaunal groups, including such taxa-of-interest as Acoela (e.g. Kånneby et al. 2015), Gastrotricha (e.g. , Kånneby et al. 2012, Kånneby 2011 and Platyhelminthes (e.g. Atherton and Jondelius 2018a, Atherton and Jondelius 2018b, Atherton and Jondelius 2019, Larsson et al. 2008, Larsson and Willems 2010. The intense sampling activities during these surveys resulted in discovery of species new to science and new records for Sweden. In addition, molecular barcoding sequences for both new and known species were attained and deposited in public databases. The availability of these and other DNA sequences allowed for the creation of reference databases, which are a critical component in the recently developed metabarcoding approach to biodiversity assessment. High-throughput DNA sequencing of a subset of gene(s) extracted from environmental or bulk samples (i.e. metabarcoding; TABERLET et al. 2012) is a promising tool that allows for rapid and comprehensive examination of community compositions (Brannock andHalanych 2015, Leray andKnowlton 2016). Metabarcoding techniques have been previously used to assess metazoan colonies along autonomous reef monitoring structures (Leray and Knowlton 2015), planktonic assemblages (de Vargas et al. 2015, Lindeque et al. 2013) and marine benthic meiofaunal communities (Brannock and Halanych 2015, Cowart et al. 2015, Fonseca et al. 2014. Two prior studies used metabarcoding approaches to examine meiobenthic communities in Sweden (Haenel et al. 2017, Holovachov et al. 2017) with a primary focus on improving and testing methods of metabarcoding, including examining the effects of sediment extraction methods, the usefulness of different primer pairs and gene loci and the optimal strategies to assign and identify operational taxonomic units (OTUs).
In this study, five locations across the southern and western coast of Sweden were sampled and a high-throughput metabarcoding approach was employed to provide a snapshot of their meiobenthic communities using fragments of the mitochondrial COI and the nuclear 18S rDNA genes. How efficient was the metabarcoding approach in capturing previously-recorded species diversity? Will metabarcoding reveal hitherto undetected species, even in localities that have been extensively sampled with traditional methods? To address these questions, four meiofaunal groups previously surveyed as part of the Swedish Taxonomy Initiative and hence with available sequence reference databases were selected for further assessment and OTUs assigned to these groups were identified to species level.

Sampling
Marine sandy sediments were collected from five locations within Sweden during May-June 2016 (Fig. 1). Sites 1 and 2 are situated near the Marine labs at Tjärnö and Kristineberg where, especially at site 2, there is a long history of taxonomic meiofauna studies (see Coull andGiere 1988, Swedmark 1964 on the history of meiofauna research). The remaining three sites span the Swedish southern coastline where little previous sampling of the meiofauna has occurred.
At each site, two or three replicates of 500 ml sediment were collected at 1.5 m depth, placed in jars and transported to the laboratory for processing. Sediments were placed in a 7.2% solution of MgCl initially for two minutes and then again for an additional five minutes. After both time periods, the samples were thoroughly mixed to suspend meiofauna and lighter sediment particles and the supernatant decanted through a 125 µm sieve. Samples were then fixed immediately in 95% ethanol for molecular analyses and stored at -20°C.

Library preparation and sequencing
DNA was extracted using Qiagen's DNeasy PowerSoil Kit following the manufacturer's instructions.
Primers were selected from Haenel et al. (2017) to target a 370 bp fragment of the 18S gene corresponding to the V1-V2 region, which have been shown to be effective for the amplification and identification of metazoan diversity in previous metabarcoding studies (Haenel et al. 2017, Holovachov et al. 2017). In addition, the COI reverse 'Folmer' primer (dgHCO2198, Folmer et al. 1994, as well as a newly-designed forward primer, were used to target a 390 bp fragment of the mitochondrial cytochrome oxidase 1 (CO1) gene encompassing the 'mini-barcode' region (Leray et al. 2013). Following the dual PCR amplification method of Bourlat et al. (2016), Illumina overhang adapter sequences were appended to the primers for the first PCR (the amplicon specific PCR), while a second PCR (the index PCR) was then performed to incorporate Illumina index adapters. Suppl. material 1 lists all primer information.
The amplicon PCR reactions were performed using 0.2 ml PuReTaq Ready-To-Go PCR Beads (GE Healthcare) with 5 pmol each forward and reverse primers and 3 µl DNA and cycling conditions of 5 min at 95°C, 35 cycles of (30s at 95°C, 90s at 50°C, 60s at 72°C) and 10 min at 70°C. PCR reactions were checked on a 2% agarose gel and purified with Agencourt AMPure XP magnetic beads (Beckman Coulter).
The index PCR reactions were performed using 0.2 ml PuReTaq Ready-To-Go PCR Beads (GE Healthcare) with 5 pmol each of Nextera XT Index Primer i5 and Nextera XT Index Primer i7 and 13 µl DNA. Cycling conditions consisted of 5 min at 95°C, 10 cycles of (30s at 95°C, 30s at 62°C, 30s at 72°C) and 10 min at 70°C. PCR reactions were again checked on a 2% agarose gel and purified with Agencourt AMPure XP magnetic beads (Beckman Coulter) using a 0.8 ratio to select for > 200bp fragments.
In order to minimise random sampling error, three libraries were created and run independently, each multiplexed to include amplicons from every sample. Libraries were pooled to equimolar concentration and sent to SciLifeLab (Stockholm, Sweden) for sequencing via the Illumina MiSeq platform with v3 chemistry. A total of 41,410,434 and 20,484,405 paired end reads of appropriate length were produced for 18S and CO1, respectively (Table 1) Table 1.
Total number of reads per marker at each sampling location before and throughout the DADA2 quality control, chimera removal and OTU grouping process.

Alpha and Beta diversity analyses
Overall species richness was calculated for each location in QIIME2 using the nonparametric Chao1 index (Chao 1984) with rarefied datasets to correct for bias due to unequal sampling size. One sample of the COI dataset with a very low yield (433 total sequences and 7 total OTUs) was excluded prior to all analyses. Rarefaction was performed without replacement to a depth of 141,711 for 18S sequences and 269,894 for COI sequences and was equal to the number of sequences in the smallest sample for each dataset. Alpha diversities amongst sites were compared by Analysis of Variance (ANOVA).
Principal coordinate analysis (PCoA) and Analysis of Similarities (ANOSIM) tests were performed to assess Beta diversity and differences in sample community composition. Distance matrices of the rarefied 18S and CO1 sample datasets were calculated in QIIME2 based on the Jaccard index (Jaccard 1901) and Analysis of Similarities (ANOSIM) tests were performed both pairwise and across all groups with 999 permutations. Example scripts for all alpha and beta diversity analyses are presented in Suppl. material 2.

Preliminary OTU assignment
For 18S sequences, preliminary OTU taxonomy was assigned using QIIME's featureclassifier classify-sklearn with SILVA release 128 at 99% (Quast et al. 2013) and default settings. This identifies query sequences to phyla, based on similarity levels of 80% and to species at 97%. A total of 281 OTUs were identified as sequences belonging to four, previously relatively well-studied taxa of interest (Acoela, Gastrotricha, Macrostomorpha and Rhabdocoela) and further processed for phylogeny-based taxonomic assignment.
COI sequences were blasted against the full NCBI nucleotide database (http:// blast.nlm.nih.gov/Blast) and MEGAN v 6.17 (MEtaGenomics Analyzer;Huson et al. 2011) was used to parse the results and assign preliminary OTU taxonomy (min. support 1; min. score 100; top percent 10; min. complexity 0). As with the 18S OTUs, a total of 58 COI OTUs were identified as sequences of interest and processed for phylogeny-based taxonomic assignment.

Overall community composition
Illumina MiSeq produced at total of 41,410,434 raw reads of 18S and 20,484,405 raw reads of CO1, which were reduced following the quality filter step to 9,691,423 and 9,043,200 reads, respectively (Table 1). Reads were clustered into a total of 3,615 18S and 2,276 COI representative OTUs and an 80% similarity threshold was used to determine phylum level taxonomy. All MiSeq data underpinning the analyses reported in this paper are deposited at the GenBank SRA under project number PRJNA627723 (https:// www.ncbi.nlm.nih.gov/bioproject/PRJNA627723).
For the 18S dataset, the majority of OTUs were assigned to either Metazoa (1,639 or 45%) or the SAR superphylum (1,313 or 36%), while 247 (6.8%) were unidentified (  A summary of the higher-level identification results from each sample based on the 18S gene sequences. The numbers of OTUs and reads of each taxon overall and at each sampling location are listed. Table 3. A summary of the higher-level identification results from each sample, based on the COI gene sequences. The numbers of OTUs and reads of each taxon overall and at each sampling location are listed.   When assessing total number of reads for the 18S dataset, Arthropoda attains more than half the numbers of reads for all of Metazoa (4,621,137 of 9,086,841 total number of reads; 51%), with Nematoda (2391247; 26%) and Platyhelminthes (1,487,606; 16%) following in abundance (Table 4; Fig. 3a,b). All other phyla represent only 6.4% (476,080) of the total number of metazoan reads, with only 0.1% (9,886) remaining unidentified.
For COI, a large amount of the biodiversity could not be assigned to a phylum: approximately 19% (134)

Community composition at different locations across Sweden
18S. Fig. 4 shows the alpha diversity rarefaction and box plots for the 18S sequences of each sampling site. Though alpha diversity was consistent across all other sites, Kåseberga had a significantly higher number of OTUs and Chao1 diversity (p < 0.05) than the remaining locations. As can be seen in  3 show the Metazoan composition of each location according to amounts of 18S OTUs and reads. Within the metazoan phyla, Arthropoda, Nematoda, Platyhelminthes and Xenacoelomorpha were present in every sample, with Annelida, Cnidaria, Gastrotricha and Rotifera further present at every locale. Such diversity of interstitial taxa is expected in shallow, sandy marine waters.
A number of 18S OTUs were common along the sampled part of the Swedish coast (Table  7). Overall, a total of 37 OTUs were present in every locality, with two of these (one Arthropoda, one Nematoda) present in every single sample. The majority of these shared OTUs (28/37) were identified as Metazoa with the remaining either members of the SAR super assemblage (4), Archaeplastida (3) or Fungi (2). Diversity measurements, based on 18S gene sequences, including (A) rarefaction curves; (B) Chao1 diversity box plots; (C) principal coordinate analysis (PCoA) emperor plots for each sample site.

Gene
Group 1 Fig. 5 gives the alpha diversity rarefaction and box plots based on the CO1 sequences for each sampling site. Although the Chao1 diversities of Tjärnö and Fiskebäckskil were overall slightly higher and that of Halmstad slightly lower, the alpha diversities based on CO1 sequences did not significantly differ between any of the sampling sites (p ≥ 0.05) and results from the pairwise ANOSIM analyses (Table 6) found the composition of each location differed significantly only between Landön and Tjärnö/ Fiskebäckskil and Kåseberga and Tjärnö/Fiskebäckskil. As with the 18S results, the PCoA emperor plot found four distinct clusters (Fig. 5c) representing Halmstad, Landön, Kåseberga and Tjärnö+Fiskebäckskil.  Table 7.

CO1.
OTUs of 18S and COI sequences that were present in at least one sample from all five localities. **OTUs that were further present in every sample.  (Table 5; Fig. 3c,d). Ten CO1 OTUs were present in every location sampled, of which eight were identified as Metazoa (5 Arthropoda, 2 Platyhelminthes, 1 Mollusca) with the remaining two unable to be identified (Table 7). One Platyhelminthes OTU (further identified as a Proseriate) was recorded from every sample.

Species identification within taxa of Interest
Acoela. All but two of the 37 18S-based OTUs initially assigned to Acoela could be positively identified to species level following the phylogeny-based taxonomy assessment (Table 8) Table 8.
Alignment-based taxonomy assignments of taxa of interest based on 18S sequences. Species of Acoela, Gastrotricha, Macrostomorpha and Rhabdocoela that were identified are listed along with the total number of OTUs and number of reads assigned to each species from each sampling locality and overall. Species in bold are identified species with at least one OTU present at more than one sampling location. **Species that were previously found within Sweden.   Gastrotricha. For 18S sequences, 10 of the 16 OTUs assigned to Gastrotricha were identified to the species level, with eight OTUs being placed with high likelihood as Tubanella cornuta and the remaining two OTUs as Macrosdasys sp. 2 and Halichaetonotus Table 9.
Alignment-based taxonomy assignments of taxa of interest based on COI sequences. Species of Acoela, Gastrotricha, Macrostomorpha and Rhabdocoela that were identified are listed along with the total number of OTUs and number of reads assigned to each species from each sampling locality and overall. Species in bold are identified species with at least one OTU present at more than one sampling location. **Species that were previously found within Sweden.
paradoxus (Table 8). Four OTUs were placed with high cumulative likelihood in singlegenus clades (one each in Chaetonotus and Paraturbanella; two in Halichaetonotus) and one OTU was placed within a clade representing a monophyletic family (Turbanellidae). The last OTU could only be assigned to the order Chaetonotida.
Only three COI sequences were assigned to the phylum Gastrotricha and none could be identified to species level (Table 9).

Macrostomorpha.
Fourteen of the 228 18S-based OTUs initially assigned to Platyhelminthes were identified as belonging to Macrostomorpha and, from these, four could be identified to species level (one OTU each Dolichomacrostomum uniporum and Psammomacrostomum sp. 1; two OTUs placed as Microstomum crildensis), with an additional eight placed in single-genus clades (Macrostomum, Microstomum). The remaining two OTUs could not be placed within a single, lower-level monophyletic clade (Table 8).
Only two of the 82 OTUs based on COI sequences were identified as Macrostomorpha.
One was identified as Bradynectes sterreri, while the other was placed in a clade of Microstomum (Table 9).
Rhabdocoela. The majority (141 of 228) of the Platyhelminthes 18S OTUs were identified as Rhabdoceola. Of these, 78 were identified to species and 21 were placed within singlegenus clades, with the remaining 42 only able to be identified to family or higher taxonomic levels (Table 8). Nineteen total species were represented, with the highest number of OTUs (19) placed with high likelihood as Schizorhynchoides caniculatus and 10 and 9 OTUs placed as Placorhynchus octaculeatus and Cicerina tetradactyla, respectively.
For COI, 29 of the 82 OTUs were identified as Rhabdocoela, and all of these were placed with high likelihood as a single species (Table 9). All but one of the 29 OTUs were identified as Austrorhynchus pacificus. The last was placed as Carcharodorhynchus sp. 24.

General Notes
This study supplements previous metabarcoding and traditional surveys of the Swedish coastal fauna. Results of both loci found high abundances of Arthropoda, Nematode and Platyhelminthes, consistent with metabarcoding studies utilising the same gene regions (Bik et al. 2011, Fonseca et al. 2014, Haenel et al. 2017, Lallias et al. 2015. Additionally, numerous other interstitial taxa such as Annelida, Gastrotricha and Rotifera were present and reflect the diversity patterns expected of shallow, sandy sediments. While results between COI and 18S genes were proportionally consistent, there were nevertheless some discrepancies. Primer choice is well known to influence results in a substantial manner. Primer bias driven by mismatches with their target has been shown to skew the relative abundance of amplified DNA from mock communities and, in certain cases, prevent rare species detection Leese 2015, Piñol et al. 2015). The COI gene is limited in usefulness because the lack of conserved primer binding sites makes it challenging to design primers that amplify consistently across taxa. The primer pair in this study was designed to amplify a broad range of invertebrate taxa (Haenel et al. 2017) with some emphasis on Platyhelminthes. It is maybe for this reason that smaller phyla, such as Kinorhynchia and Tardigrada, were not evident in the COI dataset.
However, the largest discrepancy in results from the 18S and COI datasets was between the amounts of OTUs that remained unassigned even to higher-level taxa. While the unidentified OTUs for the 18S sequences remained below 10%, nearly half (41%) were unplaced for COI. COI remains the most common DNA marker for animals, but it is less well-used with non-metazoan organisms (Hebert et al. 2003, Leray et al. 2013). There are gaps even within Metazoa, especially in less well-studied taxa, and the Genbank database still lacks COI references for numerous families of marine invertebrates (Curry et al. 2018). Thus, unassigned OTUs may be primarily attributed to the incompleteness of the COI database, although COI barcode misidentifications in GenBank or methodological artifacts (e.g. PCR and sequencing errors or amplification of pseudogenes) also potentially contribute (Bucklin et al. 2016). As the OTUs were assigned to phylum based on an 80% similarity cutoff and COI is most appropriate for higher resolution identifications (Erpenbeck et al. 2006), it is not surprising that a large number of the COI OTUs were unable to be assigned to even a phylum, and our results are consistent with or better than previous metabarcoding studies (e.g. Cowart et al. 2015, Haenel et al. 2017, Leray and Knowlton 2015. Finally, there were large discrepancies between the number of reads of a particular taxa and the number of OTUs assigned to that taxon. This can clearly be seen when looking at the 18S results (Tables 4, 5; Fig. 3), where Arthropoda and Nematoda attained fairly similar numbers of OTUs (624 or 38%, and 496 or 30%, respectively) but where Arthropoda almost doubled Nematoda in number of reads (4,621,137 or 51% compared with 2,391,247 or 26%). The influence of specimen biomass on sequence read abundance has been previously examined (Elbrecht and Leese 2015, Elbrecht et al. 2017, Thomsen et al. 2016) and a distinct positive correlation between the number of MiSeq reads and the size of an individual animal was found, such that much larger specimens were over-represented in metabarcoding studies when extracted in bulk together with smaller organisms. It may be, then, that taxa composed of smaller individuals are actually under-represented in the data (but see also Lim et al. 2016).

Community Composition at Five Locales
This study provides a snapshot of the meiofauna at five different loacations on the Swedish coast. We detected significant community composition differences amongst Locales at higher taxonomic levels, suggesting that metabarcoding even to this small extent can be useful for describing coarse level biodiversity trends. Rarefaction curves for both 18S (Fig.  4a) and COI (Fig. 5a) indicated that sequence variation within the samples was exhaustively captured. For 18S, the asymptote was reached after 20000 reads with each sample yielding between 358 and 619 OTUs and a total of 1639 metazoan OTUs for our five sample sites. Haenel et al. (2017) investigated the effects of sampling methods and primer choice on estimates of Swedish marine meiofauna. While they found lower sublittoral diversity, 266 OTUs in sublittoral sand samples from a single location, they also did not reach the asymptote after 20000 reads of 18S sequences. The alpha diversity of meiofauna along the Swedish west coast is high from a global perspective. The five sites sampled in this study yielded a higher number of OTUs than that detected by Fonseca et al. (2014) when using 18S markers to metabarcode samples from 23 littoral locations in the UK and mainland Europe or by  in 19 sites at sandy tropical beaches.
Results from the beta analyses identified four statistically-supported clusters corresponding to four distinct and diverse communities, with samples from Tjärnö and Fiskebäckskil grouping together to form a single cluster. The most strikingly different locale was Kåseberga, which unlike all other locations, was dominated by OTUs identified as belonging to the SAR super assemblage (549 18S OTUs) instead of to Metazoa (373 18S OTUs). Indeed, compared to the other Locales, Kåseberga had a much higher number of 18S OTUs from every major category except Metazoa, for which the diversity was comparatively low (Table 2). Further, though the overall biodiversity profiles within Metazoa were driven by commonly occurring taxa, several phyla (e.g. Bryozoa, Cnidaria, Echinodermata, Porifera, Nemertea) occurred predominantly or uniquely within one or two Locales (Tables 4, 5).
Alternatively, a number of 18S OTUs and COI OTUs displayed broad geographic distributions and were present at every locality. Amongst the widely-distributed Metazoa taxa, the two COI OTUs identified as Platyhelminthes (one Proseriate, one otherwise unidentified flatworm) were perhaps the most unexpected, since free-living Platyhelminthes species are thought to have limited dispersal abilities (Palmer 1988). However, without further identification, more should not be inferred at this time.

Taxa of Interest
A phylogeny-based taxonomy approach was used to identify OTUs that were preliminarily assigned to four meiofaunal taxa: Acoela, Gastrotricha, Macrostomorpha and Rhabdocoela. These four taxa of interest were selected because their biodiversity within the Swedish littoral sands is arguably well-known compared to most places worldwide. Each has been the focus of previous taxonomic research (e.g. Boaden 1960, Karling 1963, Karling 1974, Luther 1962, Luther 1963, Van Steenkiste et al. 2013, Westblad 1948, Westblad 1954, as well as modern surveys funded by ArtDatabanken and the STI (e.g. Atherton and Jondelius 2019, Kånneby et al. 2012, Larsson and Willems 2010 and DNA sequences of numerous Swedish species are available for each on GenBank. Yet, despite previous efforts, our results found OTUs from all four taxa of interest that could not be identified to species level, suggesting the possibility of new, unknown species within Sweden still remains.
Acoela. The reference alignment and subsequent tree used to identify OTUs of Acoela included a total of 343 18S and 185 COI sequences, with numerous new, unpublished sequences and sequences specifically attained from Swedish specimens collected over 25 years. Out of all four taxa of interest, the reference database for Acoela was the most complete and most targeted towards the present study. As such, it is unsurprising that so many of the 18S and COI OTUs could be identified to species level and it demonstrates the potential of the metabarcoding technique to monitor species, map their distributions and otherwise be of use in biodiversity surveys when a comprehensive reference database is available.
It is of particular interest to note, then, that despite the high amounts of reference data, there were OTUs that could not be identified to species level. The unidentified species of Arachaphanostoma and Mecynostomum based on the 18S gene could potentially be attributed to the fact that 18S does not necessarily differentiate between closely-related species (Tang et al. 2012). Archaphanostoma occulta and A. sublitorallis, for instance, have identical 18S gene sequences, though they are clearly distinguishable through 28S and COI gene sequences, as well as morphologically (Kånneby et al. 2015 There were several occasions where multiple OTUs were identified as the same species of Acoela. The DADA2 pipeline of QIIME2 clusters OTUs from Exact Sequence Variants (ESV), such that each sequence haplotype is a different OTU (i.e. even a single nucleotide difference makes a new OTU) and assigning multiple OTUs to a single species may simply reflect the intraspecific genetic variation of that species. Every one of the five species of Acoela identified in this study included multiple OTUs based on the COI gene, demonstrating both the expected variability of this gene, as well as the ability of the metabarcoding method and phylogenetic placement technique to capture and assess that variability.
However as previously stated, the 18S locus, though useful for discerning deeper phylogenetic relationships (e.g. , Holterman et al. 2006, Mallatt et al. 2004, may have limited ability to discriminate closely-related or cryptic species (Tang et al. 2012). Little or no intraspecific variation within the 18S genes of acoel species has been measured (Kånneby et al. 2015) and was expected, yet multiple 18S OTUs identified as the same species occurred seven times for Acoela, with as many as seven and nine different OTUs identified as Arachaphanostoma agile and Arachaphostoma macrospiriferum, respectively. Even if taking a conservative approach and attributing differences in OTU sequences to unusually abundant intraspecific variation or potential sequencing error, such differences should still be limited to only a few bp (Pfeiffer et al. 2018), but larger than expected variability nevertheless occurred in 18S sequences of, for instance, OTUs assigned to Mecynostomum lutheri, which differed by 15 bp (4.29%) and those assigned to Eumecynostomum macrobursalium, which differed by 19 bp (5.97%, Suppl. material 20). Further, though the majority of the nine OTUs assigned to Archaphanostomum macrospiriferum differed from each other by 1-3 bp, one OTU differed from the others by up to 35 bp (9.44%), on par with the number of differences between it and OTUs assigned to A. agile and higher than the number of differences, for instance, between OTUs assigned to A. ylvae and A. fontaneti (13-14 bp). The presence of COI OTUs that could not be assigned to any species as well as a set of 18S OTUs assigned to the same nominal species despite high sequence divergences indicates the presence of hitherto undetected diversity in the Swedish acoel fauna.
A total of 48 species of Acoela (41 described and 7 undescribed species) collected from Sweden were represented in the 18S reference alignment (currently 66 species of Acoela are known from Sweden, with the majority-63/66-recorded from the west coast; Kånneby et al. 2015). Of the 48 species represented, 25 (20 described, 5 undescribed) are known to occur at a comparable depth and habitat (interstitial sand, sublittoral ≤ 20 m depth) as the collected samples and thus may be reasonably expected to be present in the metabarcoding results. In total, 12 of these species (including 31 total OTUs) were detected in our samples (Table 8). Four other species of Acoela previously unrecorded from Sweden were also identified, including Actinoposthia sp. 8 (reference sequence collected from Helgoland), Paramecynostomum sp UJ0853 (New Caledonia), Paraproporus sp. 3 (Italy) and Paratomella unichaeta (unknown location).
Sites 1-3 on the west coast showed a much higher acoel 18S biodiversity, each having between 15-18 OTUs (9 or 10 species), as compared to Sites 4 and 5 along the southern coast (7 and 5 OTUs/species, respectively; Table 8). However, despite the high diversity, Sites 1 and 2 together were composed almost entirely of OTUs identified as species already known to Sweden, with only a single potential new record (OTU from Site 1 identified as Paramecynostomum sp. UJ08-53), a somewhat unsurprising result given the high degree of previous sampling efforts that have occurred at and near the marine labs at Kristineberg and Tjärnö. Three of the four OTUs identified as species not previously known from Sweden as well as both OTUs that could not definitively be identified were found only in the three less-well-studied southern Locales. The present study is limited in scope and the patterns are tentative, but the results from the 18S sequences indicate that the coast near Site 3/Halmstad, may represent a particularly good area for future sampling efforts, attaining both high amounts of overall acoel biodiversity (18 OTUs/9 species), as well as real potential for new species.
Gastrotricha. Gastrotricha of Sweden was relatively-recently surveyed (Curini-Galletti et al. 2012, Kånneby 2011, Kånneby et al. 2012, Todaro et al. 2011, Willems et al. 2009) and records, including some DNA sequences, were presented and published from each location. Currently, 46 species of Gastrotricha are known from Sweden, of which 35 inhabit marine or brackish environments. The 18S reference alignment included representatives from 19 of these 35 marine species (+15 of the 21 species known only from freshwaters), but only three species total, two described (Halichaetonotus paradoxus, Turbanella cornuta) and one undescribed (Macrodasys sp. 2), known from Sweden, were identified from the OTUs assigned to Gastrotricha. Six additional OTUs were found but could not be assigned to any of the Gastrotricha species in the reference alignment. As with Acoela, none the six unidentified OTUs of Gastrotricha was found in Sites 1 and 2, reflecting again how well-studied these sites are compared the southern coast.  tested the genetic variation of the 18S and COI gene in four species of Turbanella collected from the Baltic and North Seas, including from two locations on the southern Swedish coast. They found six distinct genetic clusters of T. cornuta with much larger between-group genetic distances (p-distance 12.1-18.8%) than within-group distances (maximum 1.5%) for the COI locus. Though only from a single gene, the results nonetheless reflect some potential for cryptic species. Our results further corroborate such potential, with eight OTUs based on 18S sequences all being identified as Turbanella cornuta. Five of the eight OTUs differed by no more than 5 bp from each other, as well as the three T. cornuta reference sequences, but the other two OTUs differed by 15 to 18 bp. At the very least, as with Acoela, our results suggest that there remains potentially undiscovered species of Gastrotricha within the littoral sands of Sweden.

Macrostomorpha.
There have been relative few taxonomic studies of macrostomorphs in Sweden apart from a report by Westblad (1953) and a more recent series on marine and freshwater Swedish Microstomum (Atherton and Jondelius 2018a, Atherton and Jondelius 2018b, Atherton and Jondelius 2019). Nevertheless, the majority of the OTUs identified as belonging to this group of flatworms could be identified at least to genus level and records are overall consistent with that which has been previously reported. The 18S gene has been demonstrated to be poor at distinguishing between closely-related macrostomorph species (Atherton and Jondelius 2018b), which may be the most likely explanation for the inconclusive identification of so many single-genus OTUs; however, the unidentified Microstomidae OTU based on the COI locus suggests some-albeit limited-missed biodiversity.

Rhabdocoela.
A little more than half (78 of 141) of the Rhabdocoela OTUs based on 18S data were identified to species level with a further 21 OTUs identified to a single genus. Most of the species (15 of 19) were identified from multiple OTUs and many of these (12 of 15) were collected from more than two distinct locations. The almost ubiquitous widespread geographic ranges and high genetic 18S diversity of the identified species suggests lots of untapped potential for new diversity within this group.
Juxtaposing the relatively well-curated reference database of Acoela, the Rhabdoceola COI reference database was clearly not sufficient for accurate OTU identification. Freeliving Rhabdocoela have been recorded from the littoral and sublittoral waters of Sweden (Curini-Galletti et al. 2012, Karling 1974, Van Steenkiste et al. 2013, Willems et al. 2007, Willems et al. 2009), but very few Rhabdocoela COI sequences have been submitted to Genbank (only 43 at time of writing with no sequences of Swedish specimens). For COI, 28 of the 29 OTUs were identified as Austrorhynchus pacificus, while the last placed as Carcharodorhynchus sp. 2, but further inspection showed that none of these OTUs matched the Austrorhynchus pacificus COI reference sequence by over 90%. Undoubtedly, such results are more likely due to the paucity of the COI reference data than a true identification. Metabarcoding-based species identification requires a reference database of DNA sequences that is taxonomically complete and geographically comprehensive for each species and gene region (Bucklin et al. 2016) and such is currently severely lacking for this group of animals.

Conclusions
This study contributes to the growing body of research surrounding the methodology and applications of metabarcoding. Here, Illumina MiSeq technology was utilised to examine the metazoan community composition at five locations along the Swedish coastline and to assess the biodiversity of four meiofaunal taxa therein.
Our results provide a snapshot of the meiofauna communities in the sampled localities. It is evident that metabarcoding is an effective and efficient method for assessing biodiversity, but it is contingent on the availability of a comprehensive reference database and different representations of genes within public databases can thus affect the quality of the results. Of the four taxa examined in detail, Acoela had the most complete database. Consequently, we were able to identify nearly all OTUs initially designated as acoels to the species level. We can also conclude that unknown acoel species still exist even within the most well-studied parts of Sweden. In juxtaposition, the limited depth of the Rhabdocoela COI reference database meant that even though the rhabdocoel fauna of the sampled area has been extensively studied in the past using traditional methods, metabarcoding of the COI gene failed to provide any further information beyond a basic and ambiguous OTU count.
The differing data and results between 18S and COI also underline the importance of using multiple markers in biodiversity assessments. Within each of Acoela, Gastrotricha, Macrostomorpha and Rhabdocoela, our results showed numerous instances of COI OTUs that could not be identified as a single species, as well as multiple 18S OTUs assigned to a single species. Both instances are interpreted as potential evidence of new species and, taken together, suggests that knowledge of meiofaunal biodiversity is yet incomplete, even in those areas where taxa can be considered best known. the staff at the DNA Lab of the Naturhistoriska riksmuseet, very much especially Dr. Rodrigo Esparza-Salas and Dr. Bodil Cronholm, for their efforts and help throughout the entirety of this project. This research was supported by a grant from the Swedish Taxonomy Initiative (no. DHA 2014-151 4.3) to UJ.

Funding program
This research was supported by a grant from the Swedish Taxonomy Initiative (no. DHA 2014-151 4.3) to UJ.

Grant title
Diversity and Taxonomy of Swedish Macrostomorpha