DNA metabarcoding for high-throughput monitoring of estuarine macrobenthic communities.

Morphology-based profiling of benthic communities has been extensively applied to aquatic ecosystems' health assessment. However, it remains a low-throughput, and sometimes ambiguous, procedure. Despite DNA metabarcoding has been applied to marine benthos, a comprehensive approach providing species-level identifications for estuarine macrobenthos is still lacking. Here we report a combination of experimental and field studies to assess the aptitude of COI metabarcoding to provide robust species-level identifications for high-throughput monitoring of estuarine macrobenthos. To investigate the ability of metabarcoding to detect all species present in bulk DNA extracts, we contrived three phylogenetically diverse communities, and applied four different primer pairs to generate PCR products within the COI barcode region. Between 78-83% of the species in the contrived communities were recovered through HTS. Subsequently, we compared morphology and metabarcoding-based approaches to determine the species composition from four distinct estuarine sites. Our results indicate that species richness would be considerably underestimated if only morphological methods were used: globally 27 species identified through morphology versus 61 detected by metabarcoding. Although further refinement is required to improve efficiency and output of this approach, here we show the great aptitude of COI metabarcoding to provide high quality and auditable species identifications in estuarine macrobenthos monitoring.

Introduction

Macrobenthic invertebrate surveys have been widely used for the assessment of the ecological status of aquatic ecosystems worldwide[1-5]. They are one of the key compulsory components of biological monitoring programs implemented in numerous countries’ environmental directives, such as the European Union Water and Marine Strategy framework directives (WFD 200/60/EC and MSFD 2008/56/EC) or the USA (EPA 841-B-99-002) and Canadian Aquatic Biomonitoring Network[6]. Under the WFD, for example, EU member states are required to establish a regular biological monitoring programme for freshwater systems and transitional waters (e.g. estuaries), which include macrobenthic communities[7]. So far, routine assessments of macrobenthic invertebrates have been carried out using almost exclusively morphology-based approaches for species identifications. This is a time-consuming and skill-dependent approach, which has resulted in low-throughput in processing biomonitoring samples. Very often the specimens cannot be accurately assigned to species, either because morphology-based identifications are inherently difficult, or because of organisms’ bodies that are damaged and missing diagnostic parts. In the case of immature stages, the species level identifications are extremely difficult or nearly impossible[8]. Published data on macrobenthic surveys frequently report specimens assigned only to family or genus level[9,10]. The growing reports of cryptic species among dominant macrobenthic invertebrates further calls into question the precision of morphology-based identifications[11-13]. However, many of the biotic indexes applied to macrobenthic communities require species-level identifications, such as for example the AZTI’s Marine Biotic Index (AMBI), which is based on a list of nearly 8000 species that are assigned to five ecological groups depending on each species’ tolerance to environmental disturbance[14].

Recent DNA-based approaches to species identifications from bulk community samples, such as environmental DNA barcoding or DNA metabarcoding[15,16], have the potential to help circumvent many of the above-described limitations of the morphology-based method. DNA metabarcoding is expected to help improving macrobenthic surveys, by providing a high-throughput approach that generates auditable species-level identifications. Although proof of concept studies have shown the feasibility of the application of metabarcoding approaches for monitoring river macrobenthos[15,17], no equivalent comprehensive studies have been developed specifically for estuarine macrobenthic communities. Most of the published High Throughput Sequencing (HTS) based studies in estuarine ecosystems targeted genomic regions where species level resolution is limited[18] and focused exclusively on meiofaunal communities, which used environmental DNA (eDNA) from the sediment[19-21] rather than bulk communities. So far, only a few studies have applied DNA metabarcoding to marine macrobenthic communities using the standard cytochrome c oxidase I (COI) barcode region[22-25]. Yet, these few studies either did not test the amplification success of different primer pairs in engineered communities of known species composition for methodology validation, and/or or did not target specifically estuarine soft-bottom macrobenthos. Due to their high phylogenetic diversity, marine and estuarine communities may convey additional difficulties in PCR-based approaches due to potential primer mismatch and amplification bias[26-28], therefore specific approaches must be comprehensively tested before conducting full “blind” metabarcoding assessments[22].

Our aim in the present study was a) assess the ability of the metabarcoding approach to detect the diversity of species typically present in estuarine macrobenthic communities, through the use of experimentally assembled communities; b) evaluate whether the metabarcoding approach provides comparable, more or less detailed species inventories compared to the traditional morphology-based approaches; and c) assess the ability of the metabarcoding approach to effectively discriminate among natural macrobenthic communities within an estuary, therefore reflecting environmental conditions at different sites and enabling its use in the assessment of the sediment environmental quality and ecological status of the estuary. By combining experimental and field studies, here we demonstrate for the first time the feasibility of using COI metabarcoding for monitoring estuarine macrobenthic communities, which provides equal or more sensitive data on the species composition compared to morphology.

Results

Metabarcoding of the assembled communities

Species detection success through metabarcoding of the assembled macrobenthic communities (AMC) was generally high, ranging from 78% of the species in AMC1 to 83% in AMC3, with 81% of the species detected in AMC2. Two non-target species were detected in the AMC1 although they were included in other AMC - Abra alba (W. Wood, 1802) in AMC2 and AMC3 and Scolelepis (Scolelepis) foliosa (Audouin & Milne Edwards, 1833) included in AMC2 (see Table 1). The cumulative success of recovery, considering all species present in the three AMC, was 83%; it increases to 89% if we consider the non-target species. Primer pairs were not equally effective for species detection in all three AMC. Globally, the most effective primer was D (78%), followed by B (75%), C (61%) and A (44%) (Fig. 1). Primer pairs B and D showed significant differences with the primers A and C (cophenetic correlation coefficient of the Bray-Curtis similarity = 0.73). Notably, the global maximum success rate of species detection was attainable using only two primer pairs, B and D. Six species were not detected with any primer pair, namely three crustaceans (Corophium sp.3, Cyathura carinata (Krøyer, 1847) and Melita palmata (Montagu, 1804)) and one polychaete (Scolelepis sp.). The bivalve (Abra alba) and the polychaete species (Scolelepis (Scolelepis) foliosa) have to be added if the non-target species are considered.

Table 1

Species composition of the three assembled communities and number of specimens per species.

Phylum	Class	Order	Family	Species	AMC1	AMC2	AMC3
Annelida	Polychaeta		Capitellidae	Notomastus profondus (Eisig, 1887)	1	1	3
			Maldanidae	Euclymene santandarensis (Rioja, 1917)	1	1	6
				Euclymene sp1		1	1
				Heteroclymene robusta Arwidsson, 1906	1	1	1
				Leiochone leiopygos (Grube, 1860)	1	1	1
				Leiochone sp1		1	1
				Praxillella praetermissa (Malmgren, 1865)		1	1
				Praxillella sp1		1	1
		Eunicida	Eunicidae	Marphysa sanguinea (Montagu, 1815)		1
			Lumbrineridae	Lumbrineris latreilli Audouin & Milne Edwards, 1834		1
			Onuphidae	Diopatra neapolitana Delle Chiaje, 1841		1
		Phyllodocida	Glyceridae	Glycera alba (O.F. Müller, 1776)		1
				Glycera tridactyla Schmarda, 1861		1
			Nephtyidae	Nephtyidae ni		1
			Nereididae	Hediste diversicolor (O.F. Müller, 1776)	1	1	13
		Spionida	Spionidae	Scolelepis (Scolelepis) foliosa (Audouin & Milne Edwards, 1833)		1
				Scolelepis sp		1
		Terebellida	Cirratulidae	Cirratulidae ni		1
			Terebellidae	Pista cristata (Müller, 1776)		1	1
Arthropoda	Malacostraca	Amphipoda	Ampeliscidae	Ampelisca sp		1	1
			Corophiidae	Chorophium sp1		1	1
				Chorophium sp2		1
				Chorophium sp3		1
				Chorophium sp4		1
			Leucothoidae	Leucothoe incisa (Robertson, 1892)		1
			Melitidae	Melita palmata (Montagu, 1804)	1	1	30
		Decapoda	Alpheidae	Athanas nitescens (Leach, 1813 [in Leach, 1813–1814])		1
			Diogenidae	Diogenes pugilator (Roux, 1829)	1	1	1
			Hippolytidae	Eualus cranchii (Leach, 1817 [in Leach, 1815–1875])		1
			Pilumnidae	Pilumnus hirtellus (Linnaeus, 1761)		1
			Porcellanidae	Pisidia longicornis (Linnaeus, 1767)		1
			Upogebiidae	Upogebia deltaura (Leach, 1815)		1
		Isopoda	Anthuridae	Cyathura carinata (Krøyer, 1847)	1	1	1
Mollusca	Bivalvia	[unassigned] Euheterodonta	Solenidae	Solen marginatus Pulteney, 1799	1	1	2
		Veneroida	Cardiidae	Cerastoderma edule (Linnaeus, 1758)		1	1
			Semelidae	Abra alba (W. Wood, 1802)		1	1

Gray color represents presence of the species in each AMC.

Figure 1

Species detection success for the four primer pairs (A–D). The columns in each primer pair (from left to right) denote assembled macrobenthic communities (AMC): AMC1, AMC2, AMC3 and global result for the three AMC. *Significant differences were detected between primers A/C and B/D (cophenetic correlation coefficient of the Bray-Curtis similarity = 0.72).

Natural communities

The sediments’ types in the 4 sites of the Sado estuary sampled for the natural macrobenthic communities (NMC) varied considerably in their features, ranging from sandy to muddy sediments (Table 2). The sediment at the NMC1 and NMC3 sites had respectively the lowest and highest content of organic matter (TOM) among the four sediments analyzed. Sediments of the NMC2 and NMC3 also had high organic matter content (1.30% and 2.05%, respectively), however, NMC2 had lower fine fraction (FF; particle size < 63 μm) because the site was probably disturbed due to the recent construction of the ferryboat wharf. The results are summarized in Table 2.

Table 2

Sediment features in the four sites of the Sado estuary sampled for the natural macrobenthic communities.

	NMC1	NMC2	NMC3	NMC4
Salinity	34 ± 1	34 ± 1	34 ± 1	34 ± 1
TOM (%)	0.62 ± 0.05	1.30 ± 0.11	2.05 ± 0.20	0.74 ± 0.18
FFa (%)	5.16	6.4	16.93	9

TOM total organic matter, FF fine fraction. ± of salinity indicates standard error of the refractometer. ± of TOM indicates standard deviation of three samples per site. aParticle size <63 μm.

Morphological identification of the specimens was carried out in five corer samples per site, except for NMC2, where no specimens were found after sieving one corer (NMC2.10). Species level identifications were attempted in all specimens except those taxonomically difficult groups (e.g. oligochaetes and nemerteans) and the very damaged or fragmented specimens due to the sieving process. Globally, 55 taxa were identified in the natural communities, 27 of which were identified to species level and the remaining 28 to higher taxonomic ranks. A few specimens of polychaetes (family Cirratulidae and genera Euclymene and Glycera) and amphipods (genus Ampelisca) were classified to higher taxonomic ranks, since these taxa are especially difficult to identify through traditional methods. Considering only the specimens identified to the species level, four phyla were detected in all NMC (Annelida, Arthropoda, Echinodermata and Mollusca), but this number increased to five (plus Nemertea) if we consider specimens identified to a higher taxonomic level. All communities showed a diverse taxonomic composition, comprising between 3 and 5 phyla, except NMC1, which was only composed of polychaetes and mollusks (Supplementary Fig. S1A,B).

Metabarcoding-based identification generated a total of 61 species matches in all 4 natural communities, obtained through searches against both GenBank public database and our own reference library (dx.doi.org/10.5883/DS-3150). The 61 species were distributed among six phyla, the same five reported above from the morphological identification, plus Bryozoa. The variation of the species richness among sites displayed a similar pattern in morphology or metabarcoding-based assessments (NMC2 < NMC1 < NMC3 = NMC4 for morphological identifications and NMC1 < NMC2 < NMC3 = NMC4 for metabarcoding) but the number of species recorded was more than twice using the latter method (Supplementary Fig. S1C). NMC1 was also the less taxonomically diverse together with NMC2, represented only by three phyla. Forty-three of the 61 species were detected by any of the primer pairs used (B and D). Among the remaining 18 species, 10 were recovered exclusively with primers B and 8 exclusively with the primers D. The number of reads assigned to species in each sample of all NMC and primer pair is available as Supplementary Table S1.

Comparison between morphological and metabarcoding species-level identifications in the 4 natural macrobenthic communities resulted in only 23% (range 20–28%) of the species detected simultaneously by the 2 approaches (Fig. 2). Significant differences between both approaches were detected using the Bray-Curtis similarity coefficient (cophenetic correlation = 0.73). On average, as much as 65% of the species were detected exclusively by metabarcoding (range 62–71%), while species detected exclusively by morphology were only 12% in average (range 9–15%). Among the latter, there were 4 species for which there were no reference COI barcodes available when the analysis was performed: Corbula gibba (Olivi, 1792), Ecrobia ventrosa (Montagu, 1803), Spisula solida (Linnaeus, 1758) and Parvicardium pinnulatum (Conrad, 1831). Polychaetes were the dominant taxa in all sites, regardless of whether they were identified by morphology or metabarcoding, except for morphology based identifications in NMC2, that were dominated by molluscs. The second most important groups were arthropods in the case of metabarcoding-based identifications, and molluscs in the case of morphology-based identifications. The detailed list of species identified in each site by the two approaches is available as Supplementary Table S2, while the proportion of taxa in each corer and approach is displayed in Supplementary Fig. S1.

Figure 2

Comparison between morphological and metabarcoding species-level identifications in four natural macrobenthic communities (NMC1-NMC4) of the Sado estuary, Portugal. The upper bar chart shows the distribution of the number of species per phylum obtained either by morphology or metabarcoding, in each macrobenthic community. The circles in the lower part of the figure represent the proportion of species detected exclusively by morphology (white circles), exclusively by metabarcoding (shaded circles), and by both approaches (overlapping circles dashed area). Significant differences were detected between both approaches (cophenetic correlation coefficient of the Bray-Curtis similarity = 0.73).

Figure 3 shows the graphical distribution of the natural communities as a function of their similarities in species composition, obtained by non-parametric MDS, for either the morphology, metabarcoding-based identifications and also the combination of two approaches. Three maps show a similar pattern, NMC1 and NMC2 in the left part of the map and NMC3 and NMC4 in the right side. The results obtained using AMBI also showed a similar pattern between the morphological identification, HTS and combination of both approaches, all calculated using only absence-presence of species. Conversely, the original AMBI index also showed similar results with the AMBI index using absence-presence of species for the morphology-based identifications. The four NMC were classified as slightly disturbed probably because the majority of the species obtained in each natural community through the two approaches was similar. Although NMC1 was the community closer to the EG-III (moderately disturbed) and NMC2 and NMC3 the less disturbed (see Fig. 4).

Figure 3

Non-metric multidimensional scaling (nMDS) for the morphological identification (A), COI metabarcoding (B) and morphological identification plus COI metabarcoding (C) results of the four natural macrobenthic communities (NMC). Similarity index of Bray-Curtis was applied for the absence-presence of the species.

Figure 4

Comparison of AMBI for the morphological identification using absence-presence of species (A), morphological identification using abundance of species (B), COI metabarcoding (C) and morphological identification plus COI metabarcoding (D) results of the four natural macrobenthic communities (NMC).

Discussion

The tests performed with assembled communities of known composition, showed that high success rates in species detection are attainable using COI amplicons and employing only two primer pairs. In the field tests, COI metabarcoding generated concordant results with morphology-based assessments, and detected a higher number of species in all stations and samples. The metabarcoding approach was sensitive and able to reflect differences in the species composition among natural communities and therefore it demonstrates the potential for implementing COI metabarcoding in the monitoring of estuarine macrobenthic communities.

In spite of the differences in the proportion of specimens per species, relatively high success rates in species detection were attained in all of the assembled communities (78% to 83%). AMC2, composed of the highest number of species (36), each represented by a single specimen, constituted an extreme test for the robustness of the metabarcoding approach, particularly for the effectiveness of the bulk DNA extraction and amplification procedures. In this community, no sequences were generated only for two species (Corophium sp. 3 and Scolelepis sp.) with any of the four primer pairs tested. Specimens of these species were very small (<5 mm in length) and the possibility that their DNA was not effectively isolated and that not enough template DNA was available for amplification cannot be discarded. Two species of peracaridean crustaceans (Cyathura carinata and Melita palmata) were apparently recalcitrant to amplification generating no reads, although they were present in the three assembled communities. However, the fact that previously we have been able to generate full DNA barcodes for individual specimens by Sanger sequencing using one of the primer pairs (Lobo primers[26]), and that the isopod C. carinata was recovered in the natural communities, excludes the possibility of amplification inhibition in these species. A possible explanation is that these species have a low affinity to the tested primer pair and may be outcompeted by higher affinity DNA templates from other species present in the PCR reaction. This is an important issue when considering primer match for metabarcoding studies and demonstrates the need for primer evaluation using assembled mixtures prior to large-scale analysis of bulk samples. Several reasons could explain the detection of two species (Abra alba and Scolelepis (Scolelepis) foliosa) in AMC1 where they were not included, but not in AMC2 and AMC3, where they were present. Because the organisms were processed in the same collection event, some tissue or body fragment may have been accidentally transported together with other specimens, or they may have even been preyed upon by some of the predator species (e.g. Hediste diversicolor) present in AMC1.

Mismatches between primers and target templates are a key concern in PCR-based metabarcoding, since it can lead to some level of systematic failure in species detection[29]. Because in silico analyses reveal high variability in the actual and potential primer annealing regions within the COI barcode, this marker has been dismissed as appropriate for metabarcoding[30]. Alternative markers, such as the nuclear gene coding for 18S rRNA, with lower variability in priming sites, have been used and proposed for metabarcoding marine macroinvertebrates[22,31], but the species level resolution is substantially lower than using COI[19,22,27]. When compared side by side in a field test of metabarcoding invertebrates of seagrass meadows[22], both markers showed taxonomic bias, with the 18s rRNA recovering a higher number of species (compared to full length COI barcodes (658 bp)) but amplifying preferentially meiofaunal groups such as nematodes. Since species level identification is essential for applying macrobenthic invertebrate indices (e.g. AMBI[14]), and reference libraries for marine invertebrates are available and continuously growing, the standard barcode marker for metazoans is the natural candidate for metabarcoding macroinvertebrate communities. Several studies[32] have shown that shortcomings of PCR bias may be minimized by using enhanced degenerate primers, and multiple amplification primers. The deep sequencing provided by the HTS platforms used in metabarcoding may also improve global primer success compared to what has been found using individual specimen sanger sequencing[32-34]. Dowle et al.[35] also highlighted the potential benefits of using approaches, such as gene enrichment, where the initial PCR step is not necessary.

Despite there being no major differences observed in species detection success rates among the three different assembled communities, there were considerable differences among the 4 primer pairs. The primer pair A, amplifying 658 bp, was the least successful one; hence we conclude that smaller length sizes appear to be more efficient for metabarcoding. Short fragments of COI barcode (mini-barcodes), even of 150 bp, can achieve unambiguous species-level identifications, as it was observed for a diversity of taxonomic assemblages in previous studies[33,34,36,37]. A much better success rate was obtained with primer pairs B and D compared to A and C. The two former primer combinations were tested for the first time, and proved effective in the amplification of more target species from three phyla than the remaining two primers. There was also no indication of a major taxonomic bias in these primers, as they were able to amplify targets from any of the three phyla. This indicates that, despite the large phylogenetic diversity of estuarine communities, a combined approach of degenerate primer design and multiple amplification primers can minimize substantially primer-template mismatch issues.

No relationship was found in this study between the number of specimens and the number of reads. Indeed, for phylogenetically diverse assemblages such as macrobenthic communities, comprising organisms varying widely in size, biomass and anatomically (thus varying also in the amount of DNA template available in a bulk extraction), the possibility of quantitative inferences from the number of reads data was not anticipated. For example, the polychaete species Hediste diversicolor, represented by 1 specimen in the AMC1 and 14 specimens in the AMC3, produced 8 and 23194 reads respectively, using the primer pair B. However, the polychaete species Notomastus profondus, represented by 1 specimen in the AMC1 and 3 specimens in the AMC3, produced 4601 and 3161 reads respectively, using the primer pair D. Also, in the AMC2 where all species were represented by 1 specimen, 6131 and 21576 reads of the similar-sized decapod specimens of Pilumnus hirtellus and Upogebia deltaura were respectively obtained, among various other examples of deep mismatches between the number of reads and organisms abundance and size patent in our results. Empirical relationships between specimen numbers, body size or biomass and the number of reads have been found occasionally, usually in studies targeting a closely related and more or less homogeneous group (e.g. chironomids[38,39]). In comprehensive tests performed by Elbrecht and Leese[40] such relationships were still elusive, probably because the primer efficiency is highly species-specific, preventing straightforward inference of species abundance in the assembled communities. This could be possbily circumvented by using recent techniques, such as gene enrichment, where no initial PCR step is necessary[35].

Marine macrobenthic communities are complex, highly diverse communities, where morphology-based species identifications can be rather challenging. In our study many specimens could not be identified to the species level due to uncertain species identity, mostly when they were immature stages, belonged to difficult taxa or were missing diagnostic body parts as a result of sieving, handling and ethanol. Such difficulties are common, even when a group of experts is available, as reported in numerous studies[41]. We have found many fragments of organisms, especially annelids, as a result of sieving and handling process and therefore many species could not have been identified using morphology, although they were present in the samples. A comprehensive search of over 138 published reports and inventories of benthic communities has found that approximately one third of the specimens were not identified to species level when using morphological methods, although this proportion of missed species identifications varies greatly between different taxonomic groups[42]. The morphology-based macrobenthic community profile that we obtained for the four sites in the Sado estuary, provided similar results to previous morphology-based surveys in nearby and similar ecosystems, both regarding the species richness and species-specific composition (e.g. Tagus estuary[43]).

In our study, metabarcoding approaches for identifying species composition in communities indicated that the species richness would be underestimated if we used only morphological methods. Similar findings have been reported in a study made on seagrass associated benthic communities[22], where HTS-based species inventories were considerably richer compared to morphology-based ones. The advantage of using DNA barcodes for metabarcoding approaches is that the reference libraries are being established and improved for all major groups of eukaryotic organisms. Thereby, it is possible to verify the species attribution of the samples. Contrary to the morphological approach, HTS enabled the recovery of sequences from damaged specimens, immature stages, difficult taxonomic groups, fragments of organisms and even from endoparasites, namely the copepod Mytilicola orientalis Mori, 1935 and the decapod Pinnotheres pisum (Linnaeus, 1767). M. orientalis, native to East Asia, occurs in the intestinal tracts of bivalve species and has been recorded as an alien species in European waters[44,45]. Metabarcoding could be used as a tool for early detection of invasive species[46]. P. pisum, living in the mantle cavity of bivalves, is also a parasite[47]. In addition, six species of algae were recovered: Pyropia haitanensis (T.J.Chang & B.F.Zheng) N.Kikuchi & M.Miyata, 2011, Ceramium secundatum Lyngbye, 1819, Durvillaea sp., Leathesia marina (Lyngbye) Decaisne, 1842, Petalonia fascia (O.F.Müller) Kuntze, 1898 and Scytosiphon lomentaria (Lyngbye) Link, 1833 (see Supplementary Table S2). Although the algae were not a targeted taxonomic group, this illustrates that studies with different scopes are possible, even when using the primer pairs applied in this work.

As patent in Fig. 4, the four sampling sites displayed fairly distinct sediment types and macrozoobenthic species compositions. NMC1 and NMC2 appeared to be more similar to each other (on the left side of the graphic) and the same for NMC3 and NMC4 (right side of the graphic), agreeing with their geographic vicinity, on either the north or south margin of the estuary. The species richness was consistently higher in NMC3 and NMC4 for both morphological and HTS approaches. According to the original AMBI and p/a AMBI indexes, the four NMC were classified as slightly disturbed. Sado estuary is globally considered a slightly disturbed ecosystem due to its high hydrodynamics and multiple anthropogenic activities, although the south margin is less disturbed than the north one[48]. Contrary to these global patterns, our AMBI and species richness results indicate NMC1, located in the south margin, as the most disturbed community. The NMC1 sampling site had the lowest sediment’s organic matter and fine fraction content. which would also suggest less disturbed communities. Therefore, in this case, it appears that the sediment type was not the main factor determining the macrozoobenthic assemblages composition. Regular dredging operations and construction works, together with a strong hydrodynamics, can affect and promote sudden changes in macrobenthic communities in the Sado estuary[49,50] and may help to explain these results. However, the key finding is that either morphological or metabarcoding approaches produced similar global outcomes (AMBI classifications and species richness ranks), and metabarcoding consistently outperformed morphology in the ability to detect a higher number of species and to provide species level identifications, despite the still incipient state of completion of the reference libraries for macrobenthic invertebrates.

In summary, our study demonstrates the aptitude of COI metabarcoding using HTS approach for implementation in biodiversity assessments of estuarine macrobenthic communities. High-throughput metabarcoding may enable more frequent and spatially detailed biomonitoring with higher information content[6,17], concomitantly reducing time and cost constraints in the monitoring of benthic communities. By virtue of the generation of readily comparable DNA sequence data, the metabarcoding approach can provide species-level information of high quality, with reduced ambiguity and susceptible to scrutiny in the future. The ability to provide data on parasite occurrence, for example, and to enable early detection of alien species, or to discriminate cryptic species, constitute highly relevant additional benefits of this approach. Nevertheless, further refinement is still required, to improve its overall efficiency and output, namely the improvement of the recovery rates through the refinement of primers and testing of alternative combinations, especially for the recalcitrant species. Given that the direct measurement of species abundance is still not attainable, further studies are required to generate large datasets, which will allow extensive comparison of the performance of morphology and metabarcoding-based approaches. Lastly, since the reference libraries of DNA barcodes are still incipient for marine invertebrates, the continuing completion will be decisive to fully materialize the potential of metabarcoding.

Methods

Experimental design

The sampled macrobenthic fauna does not include any protected or endangered species. This study was designed in two main sequential phases. The first phase focused on analysis of macrobenthic communities with known composition, while the second phase comprised natural field-collected macrobenthic communities. In the first phase, the ability of four combinations of primer pairs to successfully amplify fragments of the COI barcode region between 250 to 658 base pairs (bp), was tested in three assembled macrobenthic communities. The assembled communities included a different number of species and individuals per species. The two most efficient primer pairs were then used in the second phase. A schematic overview is presented in Fig. 5.

Figure 5

Schematic overview of the experimental design.

In the second phase, morphology-based taxonomic identification of the species composition in the natural macrobenthic communities was directly compared with the species inventory obtained from HTS of COI amplicons generated from bulk community DNA extractions, using two primer pairs selected among the four previously tested. This comparison was applied to NMCs collected in four separate sites in the Sado estuary, Portugal, encompassing distinct sediment features and levels of anthropogenic impact (Fig. 6). In each site, half of the replicate samples were used for conventional morphology-based identification while the remaining half was used for metabarcoding community analyses. Because data generated through metabarcoding does not provide a direct measure of specimen abundance, we used a biotic index based solely on the presence and absence of species to compare morphology and metabarcoding approaches.

Figure 6

Map of the study area showing the collection sites. (A) for the creation of artificial communities and (B) for natural communities. NMC = natural macrobenthic communities. The map was drawn using the a free and open source image editor GIMP v. 2.8, at: https://www.gimp.org/.

To enable species-level DNA based identifications from the NMC, a reference library of DNA barcodes was compiled for dominant groups of Atlantic European macrobenthic invertebrates. The reference library (dx.doi.org/10.5883/DS-3150) comprises GenBank-published[13,26,51-53] and unpublished DNA barcodes of marine invertebrates of southern European Atlantic coast, plus the DNA barcodes generated for the specimens used in the AMC study.

Sediment and specimen collection

Assembled macrobenthic communities (AMC)

Specimens were collected in the Sado (Geographical coordinates 38.49/−8.84) and Lima (Geographical coordinates 41.68/−8.82) estuaries, west coast of Portugal (Fig. 6A) during April, May, September and October 2012. Sediment samples were collected using a corer sampler (110 mm diameter, 495 mm height) and sieved through a 0.5 mm screen in order to separate the macrobenthic invertebrates. Sieved samples were transported refrigerated to the laboratory where they were individually separated from the debris and stored in absolute ethanol at −20 °C until processing. Morphological identifications to the lowest possible taxonomic level were carried out employing a stereomicroscope, using identification keys[54-57]. Species’ names were checked in the online databases World Register of Marine Species (http://www.marinespecies.org) and Integrated Taxonomic Information System (www.itis.gov). A total of 112 specimens belonging to 36 morphospecies (25 of which identified to species level) were assembled, comprising 3 mollusks, 13 crustaceans and 20 annelids species, therefore representing the 3 most dominant taxa in typical estuarine macrobenthic communities. These specimens were distributed in 3 groups in order to originate the following assembled macrobenthic communities: AMC1 was composed by 9 morphospecies of 9 specimens (one of each) (5 annelids, 3 crustaceans and 1 mollusk), AMC2 by 36 morphospecies of 36 specimens (one of each) (19 annelids, 14 mollusks and 3 crustaceans) and AMC3 by 67 specimens of 18 morphospecies (10 annelids, 5 crustaceans and 3 mollusks) (Table 1).

Natural macrobenthic communities (NMC)

Natural communities were sampled in four sites (NMC1, NMC2, NMC3 and NMC4) of the Sado estuary, west coast of Portugal (Fig. 6B) in May 2014. Geographical coordinates of each location were: 38.48/−8.88 for NMC1, 38.47/−8.86 for NMC2, 38.50/−8.84 for NMC3 and 38.49/−8.82 for NMC4. The macrobenthic communities of the Sado estuary have been extensively studied and the diversity of the soft bottom habitats and environmental impacts provides an appropriate test case for this study. NMC1 and NMC2 are situated in the Tróia Peninsula, near the protected area of the “Sado Estuary Nature Reserve”, and are generally less exposed to direct contamination sources of anthropogenic origin, except for ferryboat wharf located near NMC2. These Tróia Peninsula sites are more influenced by tidal hydrodynamism and have lower water residence time[58]. NMC3 and NMC4 are located on the north margin of the estuary, near the industrial zone close to the city of Setúbal which harbours a number of potential sources of pollution such as factories for the production of pesticides, fertilizers and pulp mill, a thermoelectric power plant, shipyards, etc.[58]. As opposed to NMC1 and NMC2, these sites have a lower hydrodynamism, therefore facilitating the retention of contaminants and sediment’s fine particles from the upper estuary. Eleven sediment samples were collected from each site (44 samples in total) using a corer sampler (110 mm diameter, 495 mm height). One sample was used for sediment’s physico-chemical characterization and the remaining 10 samples were used for macrobenthic community assessment. The latter were sieved in situ through a 0.5 mm screen in order to separate the macrobenthic invertebrates, transported refrigerated to the laboratory and stored in absolute ethanol at −20 °C until processing. Five samples were then randomly chosen for morphology-based identifications and the other 5 used for metabarcoding-based identifications. Morphology-based identifications were carried out in individually separated specimens as described in the previous section. Specimens for the metabarcoding approach were processed collectively as a bulk natural community, as described below.

Sediment characterization

At each site two sediment’s features were determined: (a) organic matter content (extrapolated from total combustible carbon, TOM): sediments were dried at 60–80 °C and combusted at 500 ± 25 °C for 4 h; and (b) fine fraction (particle size < 63 μm): determined by sieving after treating the samples with hydrogen peroxide and disaggregation with pyrophosphate.

DNA barcoding and HTS analyses

Assembled communities

Standard COI barcodes were obtained for every specimen used in the AMC study, and included in the compiled reference library. A small piece of tissue (1–2 mm) from each specimen was used for DNA extraction employing Nucleospin® Tissue kit (Macherey-Nagel Inc., Bethlehem, PA, USA) according to manufacturer’s protocols. COI was amplified using the primers LoboF1 and LoboR1 (see Table 3)[26]. PCR reactions were assembled in a 25 μL volume [2 μL DNA template, 17.5 μL molecular biology grade water, 2.5 μL 10 × Invitrogen buffer, 1 μL 50 × MgCl2 (50 mM), 0.5 μL dNTPs mix (10 mM), 1.5 μL forward primer (10 μM), 1.5 μL reverse primer (10 μM) and 0.5 μL Invitrogen Platinum Taq polymerase (5 U/μL)]. The amplification cycle was: 95 °C for 5 min; 5 cycles of 94 °C for 30 s, 45 °C for 1 min 30 s, 72 °C for 1 min; 45 cycles of 94 °C for 30 s, 54 °C for 1 min 30 s, 72 °C for 1 min; final extension at 72 °C for 5 min. PCR products were sequenced bidirectionally using an ABI 3730XL DNA sequencer.

Table 3

Primer pairs used to amplify cytochrome c oxidase I barcode fragments from bulk samples.

Pair	Primer	Direction (5′–3′)	Reference	Fragment length (bp)	PCR thermal cycling conditions
A	LoboF1	(F) KBTCHACAAAYCAYAARGAYATHGG	Lobo et al.[26]	658	(1) 94 °C (5 min); (2) 5 cycles: 94 °C (30 s), 45 °C (1 min 30 s), 72 °C (1 min); (3) 45 cycles: 94 °C (30 s), 54 °C (1 min 30 s), 72 °C (1 min); 4) 72 °C (5 min).
	LoboR1	(R) TAAACYTCWGGRTGWCCRAARAAYCA	Lobo et al.[26]
B	LoboF1	(F) KBTCHACAAAYCAYAARGAYATHGG	Lobo et al.[26]	250	(1) 94 °C (5 min); (2) 35 cycles: 94 °C (30 s), 48 °C (1 min 40 s), 72 °C (1 min); (3) 72 °C (5 min).
	250 R	(R) CTTATRTTRTTTATICGIGGRAAIGC	Shokralla et al.[33]
C	ArF2	(F) CCIGAYATRGCITTYCCICG	Gibson et al.[32]	310	(1) 94 °C (5 min); (2) 35 cycles: 94 °C (30 s), 48 °C (1 min 40 s), 72 °C (1 min); (3) 72 °C (5 min).
	ArR5	(R) GTRATIGCICCIGCIARIACIGG	Gibson et al.[32]
D	ArF2	(F) CCIGAYATRGCITTYCCICG	Gibson et al.[32]	418	(1) 94 °C (5 min); (2) 35 cycles: 94 °C (30 s), 48 °C (1 min 40 s), 72 °C (1 min); (3) 72 °C (5 min).
	LoboR1	(R) TAAACYTCWGGRTGWCCRAARAAYCA	Lobo et al.[26]

A reference Sanger based DNA barcode library was built using the COI sequences obtained for all 112 specimens. An in silico analysis was carried out based on Hajibabaei et al.[15], in order to evaluate the species level discrimination ability of the various fragments sizes. Sequences from all 112 specimens were aligned using the program MEGA v.6.0[59]. Phenograms were constructed for the complete fragment (658 bp) and two fragments of 200 bp (1–200 bp and 458–658 bp) with the Neighbor-Joining (NJ) method[60] using the Kimura 2-parameter (K2P) substitution model[61] and 1000 bootstrap replicates. Results demonstrated that unambiguous species level identifications (intraspecific divergences below 3%) are possible even for short fragments of 200 bp (data not shown).

After tissue subsampling from individuals for building Sanger based barcoding library, the rest of the specimens were then grouped in three AMC as described in Table 1, and each bulk sample (containing whole organisms) was homogenized in 95% ethanol using a conventional blender. The homogenates were incubated at 56 °C for approximately two hours to evaporate residual ethanol. Total genomic DNA of each AMC’s homogenate was extracted using Nucleospin tissue kit (Macherey-Nagel Inc.) according to manufacturer’s instructions. Four primer pairs (A, B, C and D) were used for independent amplification of either multiple fragments of CO1 barcoding region, ranging from 310 bp to 658 bp (see Table 3). PCR thermal cycling conditions for each primer pair are also presented in Table 3.

The generated amplicons from each assembled community were purified using Qiagen MiniElute PCR purification columns and eluted in 30 μL molecular biology grade water. The purified amplicons from the first PCR were used as templates in the second PCR with the same amplification condition used in the first PCR with the exception of using Illumina-tailed primers in a 30-cycle amplification regime. PCR products were visualized on a 1.5% agarose gel to check the amplification success. All generated amplicons were dual indexed and pooled into a single tube and sequenced on a Miseq flowcell using a V2 Miseq sequencing kit (250 × 2) (FC-131-1002 and MS-102-2003). All PCRs were undertaken using Eppendorf Mastercycler ep gradient S thermalcyclers and negative control reactions (no DNA template) were included in all experiments. All sequencing data generated have been deposited in the Dryad data repository (https://doi.org/10.5061/dryad.2686c).

The Illumina generated reads from all COI fragments were merged with SEQPREP software (https://github.com/jstjohn/SeqPrep) requiring a minimum overlap of 25 bp and no mismatches for all primer pairs, except for primer pair A (658 bp fragment) resulting in paired-end reads. For primer pair A, the forward and reverse sequences were quality filtered and then concatenated in a single file before taxonomic assignment. The paired-end reads were filtered for quality using PRINSEQ software[62] with a minimum Phred score of 20, window of 10, step of 5, and a minimum length of 100 bp. USEARCH v6.0.307[63] with the UCLUST algorithm was used to dereplicate and cluster the remaining sequences using a 99% sequence similarity cutoff. This was done to denoise any potential sequencing errors prior to further processing. Chimera filtering was performed using USEARCH with the ‘de novo UCHIME’ algorithm[64]. At each step, cluster sizes were retained, singletons and doubletons were not included for further analysis. Usable reads were compared against the reference Sanger based DNA barcode library (112 specimens) and assign to a species when displaying ≥98% similarity.

DNA extraction, amplification, and HTS of each natural community was carried out as described above for assembled communities. For each of the 20 bulk community samples (4 sites × 5 samples per site), two independent amplifications were performed using the primer pairs B and D. These two primers pairs were selected among the 4 previously tested in the AMC step, because the results together achieved were sufficient to obtain the maximum species recovery rates observed (see below). Amplicons obtained for each of the five samples per site were tagged separately and submitted to HTS in an Illumina MiSeq platform as described in the AMC section. After quality and size filtering, usable reads were first compared against our local barcode reference library and assign to a species when displaying ≥97% similarity. Reads without matching sequences in the reference library were then compared against GenBank using the same minimum threshold for taxonomic assignment. Only reads with a species match, either against the reference library or GenBank, were used in the remaining data analyses.

Community analyses

Non-metric multidimensional scaling (nMDS) was conducted, using PAST version 3.07[65], to show the spatial distribution of the four NMC. Bray-Curtis’s similarity index for absence-presence of species was used in order to compare morphological identification and metabarcoding data, avoiding affecting the number of null values between samples. PRIMER v7.0.13 (Primer-E Ltd, Plymouth, UK) was also used for the Bray-Curtis similarity measurement, performed to identify significant differences between the four primers used in the first step of this study and between the morphological and metabarcoding results obtained in the second step.

AZTI’s Marine Biotic Index (AMBI)[14] is a widely used biotic index to assess the quality of benthic macroinvertebrate communities considering five ecological groups (EG) to which the benthic species are allocated. EG-I: species very sensitive to organic enrichment and present under unpolluted conditions; EG-II: species indifferent to enrichment; EG-III: species tolerant to excess organic matter enrichment; EG-IV: second-order opportunistic species; and EG-V to first-order opportunistic species (V). Because the calculation of the original AMBI index requires species abundance data, an alternative AMBI based only on presence (p) and absence (a) data (p/a AMBI) must be applied when using metabarcoding-derived species inventories, as described in Aylagas et al.[27]. The classifications obtained are somewhat similar using either p/a AMBI or the original AMBI, meaning that species relative abundance does not appear to greatly affect the outcome of the benthic assessments using this biotic index[27]. Since in our study species abundances were not obtained from the metabarcoding results, we applied AMBI to the data from the morphology-based identification, metabarcoding and the combination of both methodologies, using the presence and absence of species to enable results’ comparison. In addition, the original AMBI index based on the abundance of specimens was also applied to the morphology-based identifications in order to validate the results.