Impact of mechanical Arenicola dredging on the benthic fauna communities: assessed by a morphological and molecular approach

: Lugworm Arenicola spp. dredging affects the intertidal benthic community in the Dutch Wadden Sea. Previous studies have found contradicting results regarding the recovery rates of targeted and non-targeted macrozoobenthic species, and meiobenthic communities have been neglected in these studies. The current study explores the short-term effects of dredging on both the macrofaunal and the meiofaunal communities using both a morphological and molecular approach. Benthic samples were collected right before and regularly after dredging for up to 1.5 yr at both control and fished transects. Significant differences between the control and fished transects were found by morphological and molecular approaches. Ordination analysis suggests distinct community compositions between control and fished transects in the first months after dredging. Opportunistic species with short life cycles, typically more than 1 generation yr −1 , thrived more in the fished transects during the spring/summer season compared to these species in the control transects in the same season, whereas recovery for long-lived species was slow. Both ap proaches showed similar results; however, compared to the morphological approach, the molecular approach was more sensitive to the effects due to a larger set of benthic taxa.


INTRODUCTION
Intertidal soft-sediment systems provide wideranging ecosystem services (Levin et al. 2001). Stretching across Dutch, German and Danish coastlines, the Wadden Sea represents both the largest and one of the last remaining relatively undisturbed intertidal areas in the world and appears on UNESCO's list of world heritage sites. The Wadden Sea is a shallow coastal region with large tidal flats consisting of productive soft sediment. The ecosystem is mainly driven by benthic primary and secondary production (Wolff & Binsbergen 1985, Compton et al. 2013, Christianen et al. 2017). On these mudflats, the lugworms Arenicola marina and, in smaller numbers, A. defodiens, are widespread and dominant burrowing polychaetes that engineer the sediments (Luttik huizen & Dekker 2010). Adult specimens grow up to 10−20 cm in length and can live up to 6 yr (Bijkerk & Dekker 1991). The lugworms live in a 20−40 cm deep U-shaped burrow with a funnelshaped shaft through which surface sediment slides down. Sediment particles are ingested and are then defecated at the surface. This ecosystem engineer is a known key species for the benthic mudflat community, as it is a major source of bioturbation and bio -irrigation on the intertidal mudflats of the Wadden Sea (Volkenborn et al. 2007).
Arenicola spp. are harvested by hand or mechanically in some areas as bait for sport fishing. Mechanical dredging involves digging gullies (0.4 m deep, 1 m wide, 200−500 m long) at high tide. Dredged sediment is mixed with water and sieved to harvest the larger organisms; the smaller organisms are discarded in and around the gullies (van den Heiligenberg 1987, Beukema 1995, Leopold & Bos 2009). This mechanical lugworm dredging impacts benthic communities, including lugworms as well as non-target species. Lugworm biomass and abundance are lower at dredging sites shortly after mechanical harvesting (van den Heiligenberg 1987). Different post-harvest recovery rates for the lugworm Arenicola spp. as well as for other macrozoobenthic species have been observed (McLusky et al. 1983, van den Heiligenberg 1987, Beukema 1995. Effects of dredging on the benthic communities could either result from the relocation and disturbance of sediment and smaller species or from the complete removal of lugworms in the dredged area. As the lugworm affects its surroundings by bioturbation of the sediment, the effects of completely removing this species extend beyond its absence and can be complicated. Volkenborn & Reise (2006) demonstrated a positive effect on the biomass of several benthic species shortly after their removal. However, in the long-term, lugworm presence was proven to be essential for a sustainable community (Volkenborn & Reise 2006, Volkenborn et al. 2007.
Research into the effects of lugworm dredging in the Wadden Sea have so far focused solely on macrozoobenthic species and have yielded contradicting re sults. Drent (2013) reported no differences in macro zoobenthos between dredged areas and reference areas at the Vlakte van Kerken, but at the Balgzand area in the same tidal basin, Beukema (1995) found that the biomass in dredged areas only recovered after several years. However, macrozoobenthos might not be the most appropriate indicator for the effects of dredging on the benthic community. Meiofauna species are key species in the marine food web and are known to be good indicators of ecosystem health (Balsamo et al. 2012, Zeppilli et al. 2015. With new metabarcoding techniques, it is now feasible to expand the scope of such research on the impact of lugworm dredging to also more easily incorporate meiofauna. Including this part of the community might provide better insight into the recovery rates of the local benthic community as a whole and the processes that occur during recovery. This study explored the short-term effects of mechanical lugworm dredging on benthic meiofaunal and macrofaunal communities. Benthic and sediment samples for grain size analysis were taken be fore and after digging in fished and control (nonimpacted) sites. Species composition of the prefishing and post-fishing samples was analysed by both morphological taxonomy and by using a molecular (metabarcoding) approach. The aims of the study were to (1) examine whether mechanical Areni cola spp. dredging impacts the composition of benthic communities; (2) explore recovery rates of lugworm abundances and non-targeted fauna following dredging and (3) compare morphological and molecular approaches in their abilities to detect biological change and recovery following an anthropogenic disturbance.

Field sampling
An experimental plot was created on tidal flats in the western part of the Dutch Wadden Sea, Vlakte van Kerken, north-east of the isle of Texel. The plot comprised 18 sampling stations, divided over 3 control transects and 3 experimental transects (Fig. 1). All stations were sampled once before the experimental transects were fished by the company Arenicola B.V. (T0) and another 16 times (T1−T16) after the disturbance across a temporal range of 1.5 yr (Table 1). During each sampling event, all transects were photographed, and 3 cores were taken at each station: one larger core (177 m 2 surface area, 25− 30 cm depth) for traditional morphological identification of macrofauna and 2 smaller cores (5.60 cm 2 area, 10 cm depth) for the molecular identification and grain size analysis, respectively. All cores were taken simultaneously. We chose the current method of 3 separate cores rather than splitting them, as it applies fewer treatments to the sample which might introduce contamination (Aylagas et al. 2016, El brecht et al. 2017. Grain size analysis was performed following Compton et al. (2013) and Klunder et al. (2019a).

Morphological identification
The macrofaunal communities of the intertidal Wadden Sea have been studied for decades within our research group at NIOZ. Sampling strategies and morphological identification within this research builds on the knowledge obtained throughout these years (Beukema 1974, Beukema & Cadée 1997, Beukema & Dekker 2012, Compton et al. 2013. Samples for morphological identification were washed in the field over a 1 mm mesh sieve. Sieve residues were stored in a jar in a cool and dark place and were processed within 36 h after sampling. Spe-cies retained in sieve were sorted by hand and identified while alive by experienced taxonomists based on morphological characteristics according to Hartmann-Schröder (1996) and Hayward & Ryland (1995). Molluscs, crustaceans and polychaetes were identified to the genus level, whereas Oligochaeta were identified to the phylum level.

Molecular identification
Prior to DNA extraction, extracellular DNA was removed from the samples by rinsing with a saturated phosphate buffer (Na 2 HPO 4 ; 0.12 M, pH ≈ 8). Subsequently, the samples were cryodesiccated and ground in liquid nitrogen. DNA was extracted from 10 g of the homogenized sediments using the Pow ermax Soil™ DNA isolation kit (Qiagen). The universal F04 and R22 primer pair for the metazoan community (Sinniger et al. 2016) was used to amplify a 450 bp part of the 18S gene, and both the forward and reverse primers were extended with 12 nt unique barcodes for sample identification in the bioinformatics process. PCR amplifications were performed in triplicate in a 50 μl reaction, containing 0.6 μM of each primer, 0.2 μM dNTP, 800 ng μl −1 BSA, 1 U Phusion ® High-Fidelity DNA Poly merase (Thermo Scientific), 1× Phusion ® HF Buffer (Thermo Scientific) and 5 μl of DNA extract. The thermal cycle started with an initial cycle of 30 s at 98°C, followed by 27 cycles, each comprising 10 s at 98°C, 20 s at 60°C and 30 s at 72°C, followed by a single elongation step at 72°C for 5 min. The PCR products were purified by both excising 450 bp bands from a 1% agarose gel and the Qiaquick Gel Extraction Kit (Qiagen). The final samples were quantified with a Qubit™ 3.0 fluorometer (Qiagen) and pooled in equimolar quantities together with blank PCR controls. Amplicons were submitted for sequencing at Useq (Utrecht) on an Illumina MiSeq using the 2× 300 bp V3 kit.

Bioinformatics
Raw sequences were quality filtered (≤ 30 score over 75% of the positions), de-replicated and clustered at a 98% similarity cut-off. Singleton clusters were discarded and the remaining operational taxo-19 100 meters  nomic unit (OTU) clusters were taxonomically assigned using the RDP Classifier against the SILVA 18S rRNA database (release 128; Pruesse et al. 2007) and our local reference database (Genbank accession numbers MZ709983-MZ710042). Taxonomic assignment was performed at a 0.8 and 0.5 minimum confidence threshold. OTUs identified as macrofauna taxa were extracted from the 0.8 data set; the rest were extracted from the 0.5 data set. Taxonomicassigned OTUs were clustered at the genus level if possible. If no genus-level assignment was possible, higher level assignments were used. Raw Illumina sequences were deposited in the European Nucleotide Archive (ENA accession number: PRJEB46793).

Data analysis
The mud content was relatively stable between 0.8 and 1.2%; therefore, percentage silt was used as a further measure for grain size analysis. Silt content was compared between fished and control sites using ANOVA, including the fished and control transect pairs as a random effect.
For the morphological approach, the abundance (counts) of each genus per sample was assessed. Large variations were observed in the count data of both Oligochaeta and Urothoe sp. Variations were related to observer and sampling bias, respectively and were removed from the community analysis. For the molecular approach, the relative read abundance (RRA) data per OTU was used as a measure of abundance. This transformation was chosen over rarefaction as it preserves valid read abundance data (McMurdie & Holmes 2014, Lanzén et al. 2016. Also, the RRA approach was chosen over a frequency of occurrence approach (FOO), as the sampled environment was very homogeneous, resulting in over 90% presence for many abundant species. This homogeneity would hinder the FOO ap proach in detecting differences between fished and control sites (Deagle et al. 2019). Both the count data for the morphological approach and the RRA were Hellinger-transformed (Legendre & Gallagher 2001), and dissimilarity distances were calculated using the Bray-Curtis equation. Benthic community composition was explored using multi-dimensional scaling (MDS) ordinations based on the dissimilarity distances. Subsequently, permutational multivariate analysis of variance (PERMANOVA) was employed to determine the effect of fishing on the benthic community, and a SIMPER analysis was implemented to determine which taxa were responsible for this effect.
All calculations, statistical analyses and data visualisations were performed in R version 3.5.2 (R Core Team 2018).

Field observations
A selection of the photos taken of each transect during the sampling events is shown in Fig. 2. A dredged gully was visible in the experimental transects immediately after mechanical dredging (T01). The gullies filled quickly and were covered by a layer of diatoms in the first 2 mo (spring season) (T04 and T06). The gullies remained visible as wet transects within the mudflats throughout the summer season (T08). Although only vaguely visible in the photograph, the layer of diatoms appeared again the year after, in the spring season (T13). Unfortunately, porosity of the sediment was not measured, but our observations described below are consistent with a difference in porosity between fished and control transects. All control transects consisted of solid sed-20 Fig. 2. Photographs of the most southern fished transect taken at T01, T04, T06, T08 and T13, respectively iment, whereas anecdotally, the sediment in the fished transects felt looser during sampling from T1−T13. No differences were observed after T13. It was not possible to stand on the fished transects for the entire first year of the study. The average percentage of silt for all samples within either the controlled or fished transects for each sampling event are shown in Fig. 3. Silt content was statistically correlated with treatment when time was considered in the formula (ANOVA, F 1,16 = 2.02, p = 0.01).

Community analysis
Benthic community composition was explored using MDS ordinations based on the count data per taxon from the morphological approach and the RRA for OTUs from the molecular approach (Figs. 4, S1 & S2 in the Supplement at www. int-res. com/ articles/ suppl/ m673 p017 _ supp .pdf). The ordinations from both approaches at T00 indicate no difference in composition between the control and fished transects  T00  T01  T02  T03  T04  T05  T06  T07  T08  T09  T10  T11  T12  T13  T14  T15  when assessed by PERMANOVA (Table 2). Significant differences in composition were found at T01, T05, T06 and T07 using the morphological approach and from T01−T06 using the molecular approach. SIMPER analysis for the samples in the morphological approach suggested that the primary drivers for the difference at T01 were a higher abundance of Scoloplos sp. in the fished transects and a lower abundance of Tharyx sp. Later, at T05, T06 and T07, a higher abundance of Capitella sp. and Pygospio sp. in the fished transects was observed. For the molecular approach, Simper analysis suggested that differences in community species composition were mainly due to meiofaunal taxa. At T01, higher RRA of the meiofauna arthropod taxa Podocopida and Harpacticidae was observed. For T02−T06, lower RRA for several nematode taxa within the fished transects was the main driver of the differences in community composition.

Taxonomic composition
Annelida taxa were the predominant taxa recovered by the morphological approach, while arthropods (mostly Urothoe sp. and Crangon sp.) and molluscs (mostly Limecola sp.) were detected less frequently (Table S1). The abundance of adult specimens of Arenicola, the targeted species within the experimental fished transects, was significantly lower in the fished transect (Wilcoxon, W = 10, p < 0.005) (Figs. 5 & S3). Although the abundance of Arenicola sp. was already lower pre-fishing at T00, this difference in abundance doubled between T03 and T07. On the contrary, the abundance of juvenile Arenicola was higher in the dredged transects than in the control transects at all post-fishing sampling times in the first year (Wilcoxon, W = 7.5, p < 0.005), whereas their abundance was higher in the control transects compared to the fished transects predredging. This effect was still present in the second year (T14, T16) (Fig. S3). Other smaller taxa, such as Capitella sp. and Pygospio sp., were observed in increased abundances within the fished transects (Wilcoxon, respectively, W = 16.5, p < 0.05 and Z = 15, p < 0.05; Fig. S4). Both taxa varied most markedly in their abundance between May and July in both the first and second years (respectively T6−T9 and T13−T16).
With the molecular approach, OTUs assigned to the phyla Nematoda, Arthropoda and Annelida contributed the most total number of reads (Figs. 6 & S5−S7). Overall RRA within Phylum Nematoda was lower for fished transects compared to control transects (Wilcoxon, W = 25, p < 0.05). However, none of the taxonomic orders separately showed a significant decline. In the first 2 mo after fishing, Chromadorida and Monhysterida had a lower RRA (T2−T5). RRAs of Araeolaimida were lower later in the summer in dredged transects compared to control transects (T5−T7). By contrast, RRA of all arthropod taxa, especially Podocopida, was higher in dredged transects (Wilcoxon, respectively, W = 21, p < 0.05 and W = 0, p < 0.0005). Likewise, RRA of Annelida was highest in the fished transects (Wilcoxon, W = 23, p < 0.05). Again, none of the taxa separately showed a significant difference; however, the RRA of Orbiniidae (Scoloplos sp.) and Polynoidae were higher in the first 2 mo in dredged transects compared to control transects.

Morphological versus molecular approach
The molecular approach recovered a larger number of OTUs than the morphological approach, respectively 98 versus 40. A large proportion of the extra OTUs detected were meiofaunal taxa (i.e. Nematoda, Platyhelminthes and small Arthropoda). However, from these 98 OTUs, 51 taxa were taxonomically assigned as macrofauna by the molecular approach, 11 more than detected with the morpho- For instance, the higher abundance of Annelida in dredged transects compared to control transects in the first spring season (T6−T8) was comparable between the 2 methods. However, the outcomes of the morphological and molecular approaches were not comparable with respect to time and seasonality. The RRA for Annelida within the spring season was lower than at the sampling times before, whereas the total abundance in the morphological approach was higher than at the earlier sampling times.

DISCUSSION
The present study applied both a morphological and a metabarcoding approach to analyse the shortterm effects of mechanical dredging for Arenicola spp., a bioturbating lugworm, on the locally present benthic community. Both methods indicated a slow recovery of long-lived species after dredging combined with an increased abundance of opportunistic species in the first spring/summer seasons. In line with previous research (van den Heiligenberg 1987, Leopold & Bos 2009), the observed impact was hypothesized as due to complete removal of lugworms from the sediments as well as the physical relocation and disturbance of the sediment. The combined use of the 2 methodological approaches per-23 T00 T04 T06 T08 T13

Effects of mechanical dredging
Mechanical dredging for Arenicola spp. caused prominent digging tracks within the intertidal mudflats. These tracks remained, although slightly visible, for at least 1.5 yr. Unfortunately, we do not have records beyond that period. Grain size patterns were measured as percentage silt for both the tracks as well as the control transects. Silt content was relatively low, with percentages between 2.5 and 6.9% (Klunder et al. 2019a). Significant differences in grain size patterns were measured between the control and fished samples. The percentage of silt fluctuated more in the fished samples compared to the control samples, especially in the second year; an effect which has also been described after cockle dredging (Piersma et al. 2001). Benthic species are both influenced by and influence (bioturbation) sediment characteristics. The removal of Arenicola spp., a taxon that contributes heavily to sediment bioturbation, might be the cause of the differing grain size patterns measured, including the high porosity long after dredging (Volkenborn et al. 2007).
The most prominent differences in species composition between the fished and control transects were found from April through July, the spring/summer season. MDS ordinations for both the morphological and molecular approaches showed significant differences in community composition in April and May of the first year following dredging, whereas trends in abundances for both approaches also showed differences from April through July in the second year. Life cycles of intertidal marine benthic species in temperate areas, such as the Wadden Sea, are known to follow distinct seasonal patterns, which for most species include hibernation in winter and increased activity in spring/summer due to higher availability of food resources (Beukema 1974, Coma et al. 2000. The largest differences in community composition and species abundances for both years were found within the spring/summer season -the season with increased activity -suggesting recovery rates within the disturbed environment are influenced by the spring/summer activity of the species present (Arntz & Rumohr 1982).
The taxa studied within this experiment showed different recovery rates; taxa adapted to rapid colonization -opportunistic species -can recover more quickly than long-lived and slow-growing species (Newell et al. 1998). Within this study, an increased abundance of these opportunistic taxa, such as Pygospio sp. and Capitella sp., was observed within the first spring/summer season with the morphological approach as well as an increased RRA of Spionidae and Capitellidae. These increased abundances shortly after dredging are consistent with predicted capacity for rapid recolonization (Savidge & Taghon 1988, Shull 1997, Newell et al. 1998). However, the second spring/ summer season also showed a higher abundance of these species in the fished tracks compared to the control tracks. This oscillation in species abundances during community succession has been shown before within the first 2 yr after recolonization by Arntz & Rumohr (1982), who reported that while oscillations were in phase with the seasonal oscillations, this only became normalized in the third year.
Recovery rates for long-lived and slow-growing species, such as Arenicola, were slower compared with other annelid taxa. Both count numbers of adult specimens and the RRA for Arenicolidae were lower in the fished transects until 1.5 yr after dredging at nearly all sampling events. This trend of prolonged suppression of lugworm numbers after dredging was also observed in earlier studies (van den Heiligenberg 1987, Beukema 1995. Moreover, higher recruitment of juvenile specimens was observed in the fished transects (van den Heiligenberg 1987). With the molecular approach, it is not (yet) possible to distinguish between juvenile and adult specimens; the RRA is a rough measure of biomass (Lamb et al. 2019). Adult specimens contribute more to the total biomass than juvenile specimens. Therefore, RRA will be biased based on the relative body size of adults compared to juveniles.
Explanations for differences in recovery rates and colonization mechanisms have been discussed widely, with both passive and active transport mech anisms and recruitment proposed (e.g. Arntz & Rumohr 1982, Savidge & Taghon 1988, Shull 1997, Newell et al. 1998. Taxa with larval recruitment (e.g. Capitella sp.) and mobile taxa (e.g. Pygospio sp. and Podocopida) are able to rapidly colonize a disturbed area (Savidge & Taghon 1988, Shull 1997. Once settled, these taxa grow and reproduce quickly due to the absence of competition with larger bodied, slow-growing taxa (Newell et al. 1998). These opportunistic taxa have versatile re productive strategies which allows for rapid recolonization. For instance, Pygospio sp. is able to reproduce both sexually and asexually in different seasons (Gudmundsson 1985).
Field observations and photographs in the present study suggest a high abundance of diatoms within the fished transects during spring. Although diatom layers are observed naturally on the mudflats (Scholz & Liebezeit 2012), the diatom layers on the dredged gullies were clearly visible, whereas these layers were not visible in the surrounding area and control transects. This high observed abundance is possibly related to higher nutrient availability in the fished transects compared to the control transects, resulting from decaying organisms and/or increased sediment oxidation. This in turn might contribute to increased spring− summer benthic invertebrate abundance in these fished transects (Newell et al. 1998). Especially meiofauna taxa such as Desmodorida and Podocopida demonstrated a rapid increase in RRA directly after dredging; possibly advancing from higher food availability (Boyd et al. 2000).

Comparison between methods
Both approaches demonstrated a significant effect of Arenicola spp. fishing. However, there was a difference regarding the time-range in which this effect was found. The molecular approach detected an effect slightly earlier in the sampling scheme than the morphological approach. The most obvious reason for this difference would be the difference in communities sampled (Elbrecht et al. 2017, Klunder et al. 2019b). The meiofaunal community was included with the molecular approach but was not studied using the morphological approach. The earlier response observed with the molecular approach might be due to an earlier detection of larvae and/or juvenile specimens from macrofauna taxa before the morphological approach was able to detect them, as these larvae or small juveniles would have been washed through the sieve. Another possible explanation could be the more rapid response by true meiofaunal taxa such as Nematoda and small Arthropoda because of their shorter life cycles and reproductive strategies (Boyd et al. 2000, Balsamo et al. 2012, Fonseca et al. 2014. A third possible explanation is the power of the PERMANOVA analysis. A larger number of OTUs was obtained with the molecular approach, reflecting a larger sampling ef fort, hence more power for the multivariate analysis (Ander son & Santana-Garcon 2015). The p-values derived from the PERMANOVA analysis for the morphological approach at T03 and T04 (Table 2) were non-significant, but still relatively low, whereas the molecular approach showed a significant effect. Therefore, the difference in these outcomes could be due to a difference in sampling effort (Anderson & Santana-Garcon 2015).
All taxa detected with the morphological approach were also detected at the family level with the molecular approach. However, relative quantities between the taxa differed. For Annelida, Capitella sp. was the most abundant taxon detected in the first spring− summer season in the morphological data; however, RRA of Capitellidae was low compared to other Anne lida families. Possible explanations for this are that the RRA is more directly related to biomass than species counts (Lamb et al. 2019), and the amount of DNA found in the sediment is highly influenced by ecological factors, such as seasonality or tidal movements, driving the release and dispersal of DNA (Barnes & Turner 2016, Stewart 2019. Nevertheless, patterns observed for most taxa are comparable between the 2 methods. For example, a lower abundance of adult Arenicola spp., Capitella sp. or Pygospio sp. in either fished or control transects in the morphological approach is reflected in a lower RRA of, respectively, Arenicolidae, Capitellidae or Spionidae in the same transects in the molecular approach.

Concluding remarks
The impacts of mechanical dredging for Arenicola spp. were found to last at least 1.5 yr after dredging at the fished transects. Small opportunistic taxa were observed to thrive during their known reproductive season in spring−summer within the fished transects, whereas large-bodied, long-lived taxa showed lower abundance after fishing. Both the morphological approach and the molecular approach detected these changes but the power of the latter appeared to be larger.