Samples were collected in two ecologically distinct locations along the west coast of Sweden in August 2014.

Hållö island samples: Coarse shell sand was sampled by dredging at 7-8m depth along the north-eastern side of Hållö island near Smögen, Sotenäs municipality, Västra Götalands county (N 58° 20.32-20.38', E 11° 12.73-12.68').

Gullmarn Fjord samples: Soft mud was collected using a Waren dredge at 53 m depth in the Gullmarn Fjord near Lysekil, Lysekil municipality, Västra Götalands county (N 58°15.73', E 11°26.10').

Hållö island. Hållö island samples were extracted in the lab using two different variations of the flotation (decanting and sieving) technique.

Flotation (freshwater): Freshwater was used to induce an osmotic shock in meiofaunal organisms and force them to detach from heavy sediment particles. 200 mL of sediment were placed in a large volume of fresh water and thoroughly mixed to suspend meiofauna and lighter sediment particles. The supernatant was sieved through a 1000 µm sieve to separate the macrofaunal fraction, which was then discarded. The filtered sample was sieved again through a 45 µm sieve to collect meiofauna and discard fine organic particles. This procedure was repeated three times. Meiofauna was then rinsed with seawater from the sieve into large falcon tubes. Twelve sediment samples were processed, ten of them were fixed immediately in 96% ethanol for molecular analysis and stored at -20°C. The other two samples were first screened for live representatives of Xenacoelomorpha, and later preserved in 4% formaldehyde for morphology-based identification of nematodes.

Flotation (MgCl2 solution): A 7.2% solution of MgCl2 was used to anesthetize meiofauna. As above, twelve samples were processed in total, ten of them were decanted through 125 µm sieve and fixed immediately in 96% ethanol for molecular analysis and stored at -20°C, while two samples were decanted through a 125 µm sieve which was subsequently placed in a petri dish with seawater. After 30 minutes, the petri dish as well as the inside of the sieve were searched for Xenacoelomorpha using a stereo microscope. Afterwards they were preserved in 4% formaldehyde for morphology-based identification of nematodes.

Gullmarn Fjord. Meiofauna was extracted from the Gullmarn Fjord samples using two different methods: flotation and siphoning.

Flotation (freshwater): Freshwater was used to induce an osmotic shock in meiofaunal organisms. 2.4 L of sediment were placed in a large volume of freshwater, thoroughly mixed to suspend meiofauna and lighter sediment particles. The supernatant was sieved through a 1000 µm sieve in order to separate macrofauna, which was then discarded. The filtered sample was then sieved three times through a 70µm sieve to collect meiofauna and discard fine organic particles. Meiofauna was then rinsed with seawater from the sieve into a large container and equally divided between 12 falcon tubes. Six samples were fixed in 96% ethanol for molecular analysis and stored at -20°C. Six samples were screened for live representatives of Xenacoelomorpha, and preserved in 4% formaldehyde for morphology-based identification of nematodes.

Siphoning: A total volume of 12 L of sediment was processed as follows: an approximately 5 cm thick layer of mud was placed in a container and covered with 20 cm of seawater.  The sediment was allowed to settle for 20 hours. Half of the sediment area was then siphoned through a 125 µm sieve, the residue in the sieve was immediately fixed in 96% ethanol, large macrofauna was manually removed, and the entire volume was split equally into six samples and placed at -20°C for subsequent molecular analysis. The remaining half of the area was similarly siphoned through a 125 µm sieve, the sieve contents were stored in sea water, large macrofauna manually removed, the entire volume split into six samples, which were screened for live representatives of Xenacoelomorpha, and preserved in 4% formaldehyde for morphology-based identification of nematodes.

Xenacoelomorpha. Four samples from Hållö and 12 samples from Gullmarn Fjord were used for morphology-based assessment of the diversity of Xenacoelomorpha. All samples were stored in seawater and searched for Xenacoelomorpha with a stereo microscope. All specimens found were immediately identified to the lowest taxonomic rank possible using a compound microscope equipped with DIC.

Nematoda. Two samples from each location/extraction method were used to assess nematode diversity using morphology-based identification. Samples from Hållö (flotation with fresh water and MgCl2) and Gullmarn Fjord (siphoning) were processed whole and samples from Gullmarn Fjord extracted using flotation with fresh water were subsampled by taking 1/10 of the entire sample. Formaldehyde–preserved samples were transferred to glycerin using Seinhorst’s rapid method as modified by De Grisse (1969). Permanent nematode mounts on glass slides were prepared using the paraffin wax ring method. It is common practice to estimate the diversity of marine nematodes by counting a predetermined number (usually 100 or 200) of randomly picked nematodes per sample (Vincx 1996), which may not provide sufficiently detailed results for samples with high diversity. Therefore, all nematode specimens were counted and identified for each analyzed sample. All nematode specimens were identified to genus, and, when possible, to species level.

DNA extraction. 30 samples were processed for total DNA extraction, twelve from the Gullmarn Fjord and eighteen from Hållö island, using 10g of sediment and the PowerMax^® Soil DNA Isolation Kit (MO BIO Laboratories), according to manufacturer’s instructions.

Primer design. Illumina MiSeq reagent v3. produces paired-end reads of 300bp in length, allowing a maximum marker length of 500bp when taking into account a 50 bp overlap. Universal COI primers available for the Metazoa amplify a 658bp region (Folmer et al. 1994), which is too long for most NGS applications.

Accordingly, primers amplifiying a 313 bp fragment of the mitochondrial cytochrome oxidase 1 (COI) gene were used, as described in Bourlat et al. 2016. The primers used for COI are modified from Leray et al.’s ‘mini-barcode’ COI primers (mlCOIintF-dgHCO2198; Leray et al. 2013) by adding the Illumina MiSeq overhang adapter sequences. The Leray et al. ‘mini-barcode’ primers have been shown to amplify up to 91% of metazoan diversity in a sample (Leray et al. 2013). In combination with Leray et al.'s mini barcode forward primer (mlCOIintF), we used Folmer et al.'s COI reverse primer (dgHCO2198; Folmer et al. 1994) as well as a reverse primer developed by Lobo et al., shown to enhance amplification of the COI region in a wide range of invertebrates (Lobo et al. 2013).

For the 18S region, Illumina overhang adapter sequences were appended to the primers from Fonseca et al. (SSU_FO4-SSU_R22; Fonseca et al. 2010), yielding a 364 bp fragment. These primers target a homologous region of the gene and flank a region that is highly divergent, corresponding to the V1-V2 region of the 18S gene (Lindeque et al. 2013, Fonseca et al. 2010).

Sequence overlap in the paired-end reads was calculated in Geneious Kearse et al. 2012. COI shows a sequence overlap of 230 bp and 18S shows an overlap of 190 bp.

All primer sequences used are shown in Table 2.

Illumina MiSeq library preparation using fusion primers. For Illumina MiSeq library preparation, we used a dual PCR amplification method as described in Bourlat et al. (2016). The first PCR, the amplicon PCR, uses amplicon specific primers including the Illumina adapter overhang, as described above. The second PCR, the index PCR, allows the incorporation of Illumina index adapters using a limited number of cycles (Bourlat et al. 2016).

Amplicon PCR. PCR amplifications of the COI and 18S regions were set up as follows. For a 50µl reaction volume, we used 5µl Pfu polymerase buffer (10x), 1µl dNTP mix (final concentration of each dNTP 200µM), 0.5 µl of each primer at 50 pm/µl, 2 µl DNA template (~10 ng), 0.5µl Pfu DNA polymerase (Promega) and 40.5µl of nuclease free water. Each DNA sample was amplified with the 3 primer pairs described above (COI Leray, COI Lobo and 18S). PCR cycling conditions were 2 min at 95°C (1 cycle); 1 min at 95°C, 45 s at 57°C, 2 min at 72°C (35 cycles); 10 min at 72°C (1 cycle). The PCR was checked on a 2% agarose gel. 20µl of each PCR reaction were then purified with Agencourt^® AMPure^® XP paramagnetic beads (Beckman Coulter), allowing size selection of PCR fragments by using different PCR product to bead ratios (Bourlat et al. 2016).

Index PCR. For dual indexing we used the Nextera XT index kit (96 indices, 384 samples, Illumina) according manufacturers’ instructions. Dual indexing allows an increase in the multiplex level of sequencing per lane, so that more samples can be sequenced on the same flow cell (Fadrosh et al. 2014). It also eliminates cross-contamination between samples and the occurrence of mixed clusters on the flow cell (Kircher et al. 2012). The index PCR was set up as 50µl reactions using 5µl of cleaned up PCR amplicons, 5µl of Nextera XT Index Primer i5, 5µl of Nextera XT Index Primer i7, 25µl of 2x KAPA HiFi HotStart ready mix (Kapa Biosystems) and 10µl of nuclease free water. PCR cycling conditions were: 3 min at 95°C (1 cycle); 30 s at 95°C, 30 s at 55°C, 30 s at 72°C (8 cycles); 5 min at 72°C (1 cycle). A bead purification was carried out after the index PCR with Agencourt^® AMPure^® XP magnetic beads (Beckman Coulter) using a ratio of 0.8, allowing the selection of fragments larger than 200 bp. DNA was quantified before sequencing using a Qubit Fluoremeter (Invitrogen) and average fragment size was verified using Tapestation (Agilent Technologies). Further library normalization and pooling steps are described in Bourlat et al. (2016).

Sequencing. The pooled libraries were sequenced three times independently using Illumina MiSeq Reagent Kit v3, producing in total 24 132 875 paired-end reads of 300 bp in length, of which 15 883 274 COI reads and 8 249 601 18S reads (Table 3).

Most analytical steps were performed using Qiime (Quantitative Insight Into Microbial Ecology) version 1.9.1 (Caporaso et al. 2010) and custom python scripts (Fig. 1).

Figure 1.  

Schematic workflow of bioinformatic analytical steps

The data underpinning the analysis reported in this paper are deposited at the GenBank SRA under project number PRJNA388326 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA388326).

Illumina MiSeq produced at total of 24 132 875 raw reads, of which 15 883 274 COI reads and 8 249 601 18S reads. These were quality filtered (see methods section for details) resulting in 7 954 017 COI sequences and 890 370 18S sequences. These were clustered into 2805 and 1472 representative OTUs respectively, yielding 190 metazoan OTUs for COI and 121 metazoan OTUs for 18S at 97% sequence similarity (see methods, Table 5 & Fig. 2).

Taxonomic assignment of OTUs at a 97% similarity threshold shows community composition of the samples at the phylum level (Fig. 2). Of 2805 COI OTUs, 190 (7%) were assigned to the Metazoa, 22 (1%) to plants and algae, 1 (0%) to Fungi. 2592 OTUs remained unassigned, corresponding to 92% of COI OTUs.

For the 18S dataset, 121 of 1472 OTUs (8%) were assigned to Metazoa, 104 (7%) to plants and algae, 10 (1%) to Fungi, and 8 (1%) to Protozoa. 1229 OTUs remained unassigned, corresponding to 83% of all 18S OTUs.

The large numbers of unassigned OTUs reflect the incompleteness of the databases used for COI and 18S. When unassigned OTUs are disregarded, differences between the taxonomic ocverage of the markers can be observed (Fig. 2, B and D). COI is the ‘standard’ animal barcode and is thus mostly useful for diversity surveys within the Metazoa (Hebert et al. 2003). 18S has on the other hand much larger taxonomic coverage and can be used for biodiversity profiles of whole eukaryotic communities, at higher taxonomic scales.

Of all OTUs classified as Metazoa, a detailed breakdown per phylum is presented in Table 5 and Fig. 3. Annelida (30% of CO1 metazoan OTUs and 23.97% of 18S metazoan OTUs) and Arthropoda (27.37% of CO1 metazoan OTUs and 11.57% of 18S metazoan OTUs), were the most OTU rich phyla identified in all samples combined, a similar pattern as observed in a recent study on coastal seagrass meadows in Brittany, France (Cowart et al. 2015).

As well as Annelida and Arthropoda, other phyla represented by a high number of OTUs in our samples include Mollusca (13.68% of COI metazoan OTUs and 4.96% of 18S metazoan OTUs), Platyhelminthes (10,74% of 18S metazoan OTUs and 0% of CO1 metazoan OTUs) and Nematoda (8.26% of 18S metazoan OTUs and 0% of CO1 metazoan OTUs) (Table 5 & Fig. 3). Other benthic metabarcoding studies based on the 18S V1-V2 region, found Nematoda and Platyhelminthes as the most OTU rich phyla represented (Fonseca et al. 2014, Fonseca et al. 2010), or Nematoda and Annelida (Bik et al. 2012b), alternatively Nematoda and Arthropoda (Bik et al. 2012a, Lallias et al. 2015).

Taxonomic community composition at both locations surveyed is illustrated in Fig. 4. The bar plots in Fig. 4 take into account the read counts for each OTU, whereas Table 5 and Fig. 3 do not take these into account.

In Fig. 4, clear differentiation in biodiversity between the two habitat types (soft mud versus coarse shell sand) can be observed, as expected. Echinodermata (such as Ophiurida, Echinoidea and Asteroidea), Mollusca (Bivalvia, Gastropoda), Annelida and Arthropoda are represented by higher numbers of reads in samples from the muddy sediments in the Gullmarn fjord samples (grain size 100 μm approx.).

In coarse shell sand in shallow areas, such as in the Hållö island samples, Annelida and Arthropoda are represented by higher numbers of reads, followed by Chordata (cephalohordata such as Branchiostoma sp., ascidians and various fish species such as Gobius sp., Ctenolabrus rupestris, Solea solea) with in addition a larger diversity of small taxa such as Bryozoa, Gnathosthomulida, Gastrotricha, Tardigrada, Rotifera, Sipuncula and Phoronida, reflecting the high diversity of insterstitial taxa found in sandy sediments.

A greater number of phyla were uncovered in the Hållö Island samples than in the Gullmarn Fjord samples (Fig. 4A and 4B) and this observation was corroborated by the alpha diversity rarefaction plots showing that Hållö Island samples (in red) present a higher diversity than the Gullmarn Fjord samples (in blue) (p-value = 0.001) regardless of the marker used (Fig. 5A and 5B). Within the same location, choice of extraction method does not have a significant impact on sample diversity (p-value ~ 1) (Fig. 5C and 5D, Table 6). However, for the 18S dataset, the flotation method seems to be more effective for extraction of nematodes than the siphoning method in the Gullmarn Fjord samples (Fig. 4A and 4B). Moreover, the beta diversity PCoA results highlight the fact that sample composition is influenced by the choice of extraction method for both COI and 18S datasets (p-value = 0.001) leading to four different clusters (Fig. 6and 6B, Table 6). For the COI dataset, in addition to extraction method as a factor of divergence, choice of primer (COI Leray or COI Lobo) also influences the grouping of the samples (p-value = 0.003 excluding unassigned OTUs and 0.001 including unassigned OTUs), in particular for the Hållö Island samples (Fig. 6C). Moreover, the COI Lobo primer seems to uncover a higher diversity of taxa than the COI Leray primer Fig. 5E) even if the results are considered to be non significant (p-value = 0.585 excluding unassigned OTUs and 0.111 including unassigned OTUs) (Table 6, Table 7).

Using a sequence similarity search at 97% similarity allowed us to identify 213 COI OTUs and 243 18S OTUs to species level (Table 8 and Suppl. material 1). For the COI dataset, 81 species (of which 70 metazoans) were found in both locations, 36 (of which 35  metazoans) were found in the Gullmarn fjord only and 96 (of which 85 metazoans) were found in Hållö island only. For the 18S dataset, 108 species (of which 48 metazoans) were found in both locations, 44 (of which 21 metazoans) were found in the Gullmarn fjord only and 91 (of which 52 metazoans) were found in Hållö Island only (Suppl. material 1). These species observations from metabarcoding represent 'molecular occurrence records' that could be used in monitoring and other types of biodiversity surveys, in the same way as physical observations, such as for mapping species distributions (Bohmann et al. 2014, Lawson Handley 2015).

Five alien species were detected in in the sample, of which two are considered invasive (in bold; Table 9), and the other three are on alert lists. The two invasive species (Acartia tonsa, a copepod, and Alexandrium ostenfefeldii, a dinoflagellate) could easily be overlooked in routine monitoring programs. Species within the genus Acartia are difficult to distinguish (Jensen 2010) and the invasive species can be confused with other native species. Also A. ostenfeldii is easily misidentified as other Alexandrium species; detailed thecal plate observation is often necessary for proper identification (Balech 1995).  This shows the potential of molecular techniques for monitoring  invasive species, and points to problems using traditional identification techniques. Many invasive species arrive in an area as spores, larvae or juveniles - all life stages that may be easily overlooked and problematic to identify to species level. Target barcoding of environmental DNA (eDNA) shows a great promise for detecting species without the need of costly sampling schemes. This would also allow for more random sampling in an area, increasing the probability of actually finding a species even when they occur in low numbers.

Comparison of morphology-based assessment of Xenacoelomorpha diversity with metabarcoding using taxonomic assignments to the phylum level (with 80% similarity threshold; Suppl. materials 2, 3), shows that extraction procedures have strong impact on the effectiveness of morphology-based identification (Tables 10, 11). Using freshwater for extraction of Xenacoelomorpha rendered most of them unrecognizable and unidentifiable, but left their DNA intact and suitable for metabarcoding. No identifiable Xenacoelomorpha were found in the Hållö samples extracted using flotation with fresh water, while all specimens found in Gullmarn Fjord were treated together as one taxon "Acoela sp." for the lack of better alternative. Metabarcoding, on the other hand, recovered between 6 and 15 taxa (OTUs) from the Hållö samples  extracted using flotation with fresh water (Table 11), and up to 13 taxa (OTUs) from the same type of samples from the Gullmarn Fjord site (Table 11), depending on the barcoding region used. Just like for nematodes (see below), 18S barcodes always gave higher overall estimates of diversity (number of OTUs) compared to COI (Table 11). 18S also gave higher diversity estimates, compared to morphology-based identification for the Hållö samples extracted using flotation with MgCl2 (11 versus 7), but lower for the Gullmarn Fjord site samples extracted using siphoning (9 versus 15). COI Leray primers were less effective compared to the COI Lobo primers that recovered 2-6 OTUs more in all samples (Table 11). The most numerous of the morphologically identified species, Mecynostomum tenuissimum, was present with 120 specimens in the manually sorted samples, but was not detected at all in the 18S samples. Note that the 18S and COI sequences for all of the species identified in the visually sorted samples are present in the reference database. This raises the question of the efficiacy of using the SSU_FO4-SSU_R22 18 S fragment for metabarcoding of acoelomorphs. A recent study found a number of unknown xenacoelomorph taxa while data mining metabarcoding sequences from surveys of pelagial and deep benthic habitats (Arroyo et al. 2016). Unknown xenacoelomorph species may exist also at the moderate sampling depths we sampled in the Gullmarn Fjord. Our siphoning technique relies on migration of specimens to the sediment surface in response to hypoxia. It is possible that there are xenacoelomorphs with high tolerance for hypoxia that are not captured by the siphoning method, and thus would not be found in the manually sorted samples, but could be detected by metabarcoding of unprocessed samples. It should be noted that the extraction method used on the Hållö samples does not rely on migration of specimens to the surface.