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 and Giere 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.

Figure 1.  

Location of each of the sampling sites along the Swedish western and southern coasts. (1) Tjärnö near the Tjärnö Marine Laboratory (formerly Sven Lovén Centre for Marine Sciences Tjärnö), 58°52’41.0 N, 11°06’56.0 E; (2) Fiskebäckskil near the Sven Lovén Centre for Marine Sciences Kristineberg, 58°14’52.1 N, 11°27’05.8 E; (3) Halmstad at Påarp Beach, 56°36’10.0 N, 12°54’48.3 E; (4) Kåseberga, just east of Ales stenar, 55°23’00.6 N, 14°03’46.5 E; (5) Landön western side, 55°58’20.6 N, 14°24’29.6 E.

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.

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, Leasi et al. 2018). 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).

Bioinformatic data processing was performed via QIIME2 (Quantitative Insight Into Microbial Ecology) version 2018.11 (Caporaso et al. 2010) following the procedures presented in Bourlat et al. (2016). Sequence quality was initially assessed using FastQC (Andrews 2010), with primer removal, further quality control, chimera removal and OTU grouping performed via DADA2. Example scripts are presented in Suppl. material 2.

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.

For 18S sequences, preliminary OTU taxonomy was assigned using QIIME’s feature-classifier 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.

Reference alignments and phylogenetic trees were created in order to identify OTUs of the four taxa of interest and were constructed using 18S and COI gene sequences downloaded from GenBank. Datasets include 1) for Acoela, a total of 343 (18S) and 185 (COI) reference sequences, including a combination of sequences downloaded from GenBank, as well as new, previously unpublished sequences; 2) for Gastrotricha, 174 (18S) and 148 (COI) sequences, based on the combined dataset of Kånneby et al. 2013, Kieneke et al. 2012, Kolicka et al. 2018; 3) for Macrostomorpha, 239 (18S) and 95 (COI) sequences, based on Atherton and Jondelius 2019; and 4) for Rhabdocoela, 393 (18S) and 43 (COI) reference sequences representing all Rhabdocoela sequences currently available on GenBank. All reference alignments were carried out via MAFFT v 7 (https://mafft.cbrc.jp/alignment/server; Katoh et al. 2019) and the maximum likelihood (ML) trees were built using RAxML v. 8.2.10 (Stamatakis 2014) with a GTR-GAMMA substitution model and 1000 bootstrap replicates. Specimen information, including GenBank accession numbers, can be found in Suppl. material 3 .

Query sequences were aligned against the appropriate reference database using MOTHUR v 1.39 (Schloss et al. 2009) with default settings. Taxonomic predictions were generated via PPlacer v 1.1 (Matsen et al. 2010). Taxonomic assignments were based either on high likelihood (above 90% threshold) of single placement or on high cumulative likelihood of multiple placements when all are within a single monophyletic clade (see Holovachov et al. 2017). All reference alignments and Placement trees are available as Suppl. materials 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19.

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 (Table 2; Fig. 2a,b). For COI, 708 (31%) OTUs were assigned to Metazoa, 185 (8.1%) to the SAR superphylum, and 13 (0.57%) to Fungi, while a large portion of the OTUs (939, 41%) remained unassigned (Table 3; Fig. 2 c,d).

Figure 2.  

A graphical summary of the higher level composition for all samples combined by (A) the percent of OTUs as determined by 18S sequences; (B) the percent of sequence reads as determined by 18S sequences; (C) the percent of OTUs as determined by COI sequences; (D) the percent of sequence reads as determined by COI sequences.

Tables 4, 5 and Fig. 3 present the OTUs assigned to any Metazoan phylum. For 18S, more than two-thirds of the metazoan diversity based on number of OTUs was split relatively evenly between Arthropoda (624 of 1,639 total OTUs; 38%) and Nematoda (496 OTUs; 30%), with Platyhelminthes accounting for 14% (228 OTUs) and all remaining phyla accounting for only 11% (175 OTUs). A total of 7% (116 OTUs) of the Metazoan OTUs were unable to be identified to phylum level. Such results are consistent with other benthic metabarcoding studies based on the same gene region (Bik et al. 2011, Fonseca et al. 2014, Fonseca et al. 2010, Lallias et al. 2015).

Figure 3.  

A graphical summary of the Metazoa OTUs identified to phylum level for each sample, as well as all samples combined, as shown by (A) the percent of OTUs based on 18S sequences; (B) the percent of sequence reads based on 18S sequences; (C) the percent of OTUs based on COI sequences; (D) the percent of sequence reads based on COI sequences. 1: Tjärnö, 2: Fiskebäckskil, 3: Halmstad, 4: Kåseberga, 5: Landön. Total: all samples combined.

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) based on OTUs and 41% (3,401,317) based on the number of reads. Phyla with high biodiversity in our samples include: Arthropoda with 265 OTUs (37%) across 2,230,455 sequences (27%); Nematoda with 106 OTUs (15%) and 591,410 sequences (7.2%); Platyhelminthes with 82 OTUs (12%) across 1,539,916 sequences (19%); and Mollusca with 57 OTUs (8.1%) and 353,181 sequences (4.3%; Table 5; Fig. 3c,d).

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 Table 2, Kåseberga attained much higher numbers of OTUs from every major category except Metazoa. Metazoan diversity of Kåseberga is relatively low, with only 373 OTUs across ~1.5 M sequences. Pairwise ANOSIM analyses (Table 6) indicated that sample composition differed significantly (p ≤ 0.011, R ≥ 0.31) between all sites, excepting Tjärnö and Fiskebäckskil (p = 0.115, R = 0.27). This is illustrated very clearly in the PCoA emperor plot, which shows four distinct clusters representing Halmstad, Kåseberga, Landön and Tjärnö+Fiskebäckskil (Fig. 4c).

Figure 4.  

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.

The majority of the diversity of every sampling location was assigned either to the SAR supergroup or to Metazoa. Halmstad in particular had much more diversity of Metazoans (619 OTUs across ~ 2.17 M reads) as compared to the other sites (358-465 OTUs across ~ 1.47-2.34 M reads), largely due to the abundance and diversity of the Arthropoda, which was higher here than any other site samples (316 OTUs, ~ 1.47 M reads). Table 4 and Fig. 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).

CO1. 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.

Figure 5.  

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.

OTUs unidentified even to the level of larger assemblages accounted for a large amount of the CO1 datasets for every locale, fully dominating at three of the five locations. Otherwise, the majority of the known diversity was assigned primarily to Metazoa, although fewer metazoan phyla, in general, were represented at each site as compared to the 18S dataset. Of the metazoan phyla represented in the CO1 datasets, Arthropoda, Mollusca, Nematoda and Platyhelminthes were present in every sample, with Annelida and Xenacoelomorpha present at every locale (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.

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). Overall, 15 different species of Acoela were identified from the samples, with multiple OTU placements occurring for seven species (2 each for Anaperus tvaerminnensis, Eumecynostomum macrobursalium, Mecynostomum lutheri and Philactinoposthia sp. 3; 3 for Arachaphanostoma sp. 1; 7 for Arachaphanostoma agile and 9 total OTUs identified as Arachaphostoma macrospiriferum). The remaining two OTUs were each placed with accumulated likelihoods above 0.90 in a monophyletic clade representing a single genus (Arachaphanostoma and Mecynostomum).

For COI, 20 of the 25 Acoela OTUs were identified to individual species level (Table 9), with three additional OTUs placed within a single genus clade (Philactinoposthia). Of the remaining two, one was identified to family level only (Mecynostomidae) and the other was unable to be placed within a single monophyletic clade with high cumulative likelihood. Overall, five species (Arachaphanostoma agile, Arachaphanostoma macrospiriferum, Arachaphanostoma sp. 1, Paedomecynostomum bruneum and Philactinoposthia sp. 3) were represented, all with multiple OTU assignments per species.

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 paradoxus (Table 8). Four OTUs were placed with high cumulative likelihood in single-genus 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 single-genus 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.

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, Fonseca et al. 2010, 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 (Elbrecht and 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).