Impact of a coastal protection measure on sandy-beach meiofauna at Ahrenshoop (Baltic Sea, Germany): results from metabarcoding and morphological approaches are similar

Sand nourishment is a common and often applied coastal protection measure along the beaches of the German Baltic Sea coast. At the tourist beach of Ahrenshoop, large-scale sand nourishments are completed every five years along a 4-km stretch. The nourishment in question was executed between the end of October 2021 and mid-February 2022. Meiofaunal communities at the beach-water interface were monitored at six stations before (September 2021) and for one year after the impact (March and September 2022, March 2023). To assess the effects of sand nourishment on meiofauna communities, two complementary approaches were employed: whole community metabarcoding and morphological determination and counting of the animals. Both methods produced consistent results regarding community response to disturbance and subsequent recovery. During the period from autumn 2021 to spring 2023, drastic changes in the composition of communities were observed. Directly after the impact, the taxa Acari and Annelida nearly disappeared, abundances of Copepoda decreased substantially, and Platyhelminthes increased notably. The fact that the studied beach held a low diversity may have facilitated that metabarcoding and morphological determination followed the same trend. After initial studies comprising both methods, efficient assessments using only metabarcoding may be considered for comparable sites and situations.


Introduction
Climate change is causing intensified storms and sea-level rise (Carley and Cox 2017; IPCC Intergovernmental Panel on Climate Change 2019).Overwash, flooding and wave-induced currents displace sediments and erode sandy beaches.However, the primary drivers of beach development and erosion are coastal morphology, sand composition, and long-term hydrodynamics (Jarmalavičius et al. 2016;Mattheus et al. 2018).In the southern Baltic Sea, coastal currents transport beach sand in an east-northward direction, forming sandy spits that extend into peninsulas with lagoons at their rear.The spits are affected by erosion, and countermeasures are taken to ensure the integrity of Metabarcoding and Metagenomics 8: 211-234 (2024), DOI: 10.3897/mbmg.8.127688 Iryna Kapshyna et al.: Sandy-beach meiofauna: metabarcoding and morphology the beaches and seaward dune systems (Boniecka and Kubacka 2024).Apart from the threat to human settlement, land use and tourism (Defeo et al. 2009), the absence of protective measures in the Baltic Sea could result in a breach in the coastal defence and the subsequent destruction of the ecologically important coastal lagoon systems.Sand nourishment, also known as sand or beach replenishment, is cost-effective and widely used for coastal protection (Staudt et al. 2021;Kindeberg et al. 2023).However, as emphasized by Defeo et al. (2009) sand nourishment is not only a protection measure but also a stressor for sandy beach ecosystems.
Animals of the meiofaunal size class (32-1000 μm) are the most abundant metazoans in the marine benthos.Meiofauna play an important role in benthic food webs.They have no pelagic larval stages and their whole life cycle is sediment bound.Due to their diminutive size, meiofauna organisms are readily transported by currents when resuspended by wave action.Only few meiofauna taxa actively swim and enter the water column.Changes in meiofauna composition can be used as a bioindicator to identify the effects of various forms of ecosystem disturbances, including those caused by anthropogenic impacts (Zeppilli et al. 2015;Schratzberger and Ingels 2018).
The diversity (and the high abundances of some) of these taxa hamper traditional morphological assessments.Their identification to species level is extremely demanding mostly due to a lack of identification keys or expertise.Therefore, ecological studies of meiofauna often are limited to the determination and counting of higher taxa (e.g., Menn 2002;Kotwicki et al. 2005;Kuhnert et al. 2010;Veit-Köhler et al. 2018).Genetic tools, such as whole community metabarcoding, make it possible to rapidly determine entire communities, sometimes to the species level.However, the outcome derived from metabarcoding is considered qualitative, due to discrepancies and a lack of correlation between the number of reads of similar sequences and the actual counts of individuals.Despite various efforts to enhance the quantifiability of metabarcoding (e.g., Elbrecht and Leese 2015; Lamb et al. 2019;Rossel et al. 2019), this challenge remains unresolved.
We monitored the effects of a large-scale sand nourishment on meiofauna at the beach of the seaside resort Ahrenshoop (Fischland-Darß-Zingst Peninsula, Baltic Sea, Fig. 2).To our knowledge, so far no such study has been carried out on the non-tidal coasts of the Baltic Sea.Following the current coastal protection strategy, approximately every five years sand is pumped onto the beach which is widened up to 40 m.The last replenishment took place from the end of October 2021 to mid-February 2022.
The objectives of this study were as follows: (1) To monitor and compare the meiofauna communities in the swash zone (beach-water interface) along the beach of Ahrenshoop before and after the sand nourishment.(2) To compare the results of the morphological and genetic approaches and provide advice for future rapid assessment of the effects of disturbance.

Study area and sampling stations
The mainland of the seaside resort Ahrenshoop (Fischland-Darß-Zingst Peninsula, Baltic Sea) is protected from adverse influences of the Baltic Sea by a row of dunes and the beach area in front.Approximately in a 5-year interval, large-scale sand nourishments are completed at the beach.Between the end of October 2021 and mid-February 2022 650548 m 3 of sand were transferred from a 10-m deep off-shore location to a 3925 m long section of the Ahrenshoop beach (information from Staatliches Amt für Landwirtschaft und Umwelt Mittleres Mecklenburg).The beach was widened by 30 to 40 m and the dunes reinforced.At the water line, the layer of sand added was approximately 1 m thick.Our sampling stations (AH01-AH05) were located along the Ahrenshoop sandy beach (Fig. 2) from the boundary of the nature reserve in the north east (54.39°N,12.44°E)to a site just north of the breakwater (54.37°N, 12.41°E).These five stations lay within the zone impacted by the sand nourishment.An unaffected reference station was located south of Ahrenshoop (close to Niehagen, 54.35°N, 12.39°E) at the end of the road Pappelallee (PAP).The position of the sampling locations at each station varied due to the artificial extension of the beach (between T0 and T1) and the subsequent erosion of the added sand (T1-T3).In a schematic cross section, we visualise the temporal changes in the Ahrenshoop beach area (Fig. 3).

Sampling and sample processing
Samples were collected at four dates.The first sampling was carried out before the sand nourishment took place (T0: 14 and 16 September 2021).Three samplings were realised after the impact: T1 (23 March 2022), T2 (27 September 2022), and T3 (28 March 2023).Three replicate samples per station were collected at T0 and four replicates per station were collected at T1-T3.Samples were taken at the water line (beach-water interface) in the middle of the beach section between two groynes.Plexiglass cores (inner diameter 5.4 cm, surface area 22.9 cm 2 ) were placed randomly and inserted vertically into the sediment down to 15 cm depth.Cores were closed with a rubber stopper on top, dug out and closed below with another stopper.Each core was sliced in 5 cm-layers (0-5, 5-10 and 10-15 cm).Sediment horizons were stored in 500-ml Kautex vials and preserved in 96-99% denatured ethanol.After 24-48 hours samples were carefully decanted over a 32-µm sieve, re-fixed with fresh ethanol and stored at -20 °C until further analyses.
To extract the meiofauna, the samples were transferred to an Erlenmeyer flask (1 l) which was filled with filtered tap water and closed with a lid.The flask was vigorously shaken and the liquid content decanted over a 32-μm sieve.This extraction process was repeated four times.After processing each sample, all equipment was washed with a 5% bleach solution and subsequently rinsed with tap water to prevent cross-contamination.Furthermore, negative control samples were conducted with tap water instead of sediment using the previously described method.The extracted meiofauna were fixed with 96-99% denatured ethanol and stored at -20 °C.

Morphological determination of meiofauna
Samples of four stations were used for morphology-based analyses of the meiofauna communities.Stations AH01, AH03, and AH05 were selected because of their location at both ends and in the centre of the impacted zone.Station PAP was included as the control site.One core per station and sampling date (T0-T3) was analysed.
The fixed samples were poured over a 32-µm sieve and transferred to tap water for sorting.Adult copepods were picked and mostly transferred to glycerine.Several specimens were stored in 96% ethanol (-20 °C) for further studies.Then the samples were stained with Rose Bengal.A Leica MZ125 stereo microscope was used for determination of meiofauna higher taxa (classification using the keys of Higgins and Thiel 1988;Schmidt-Rhaesa 2020).The group "Copepoda" consists of adult harpacticoid copepods and copepodid developmental stages.Copepod nauplii were counted separately because of the different ecological niche they occupy due to their small body size.Platyhelminthes and Annelida of the genus Diurodrilus Remane, 1925 were grouped as "Platyhelminthes+Diurodrilus sp.".Due to their morphological similarity, they could not be distinguished during sorting.Unidentified organisms were categorised as "others" and were excluded from the analyses.The total number of individuals per higher taxon was counted and is presented as individuals per 10 cm 2 .Figure 3.A schematic cross-section of the Ahrenshoop beach area shows the beach extent before and its changes after the sand nourishment.Sediment samples for the meiofauna monitoring study were collected at the water line.Immediately before the sand nourishment the beach had its smallest extent (time and sampling location T0, 09/2021, beige dune and sand).After the sand nourishment the beach reached its maximum extent (T1, 03/2022, burgundy coloured sand).We followed the changes of the location of the water line for two more sampling dates: T2 (09/2022) and T3 (03/2023) where beach erosion and loss of sand were observed (dark and light pink sand).Graph created with Inkscape 1.0.2-2.Meiofauna specimens were imaged using a confocal laser scanning microscope.Selected individuals were stained with a watery solution of acid fuchsin and Congo red with glycerine (modified after Michels and Büntzow 2010) and kept in the dark for at least 12 hours.The specimens were mounted in glycerine on slides and scanned using a Leica CTR 5000 (DZMB, Senckenberg am Meer) equipped with a Leica DM 5000 B microscope (200-400× magnification) using three visible light lasers (10% Argon, DPSS 561 and HeNe 633).The Leica LAS software (Leica Microsystems) was used for image processing.Adobe Photoshop Elements 14 was applied to correct contrast and brightness and to create Fig. 1.

Metabarcoding of meiofauna communities
The extracted 198 meiofauna samples and 23 negative control samples (see above) were filtered over sterile glass filters with a pore size of 2.7 μm.An additional 8 blank samples were introduced during the filtering process and distributed among the meiofauna samples.After each sample, all equipment was washed with a 5% bleach solution and rinsed with tap water to prevent cross-contamination.Each filter containing meiofauna and each control filter was transferred to a separate DNA-free Eppendorf tube and further dried using a speed-vacuum system at 45 °C for 1 hour.Genomic DNA was extracted from the filters using the DNeasy PowerSoil pro kit (Qiagen), following their protocols.Realtime-PCR was performed to amplify V1&V2, two hypervariable regions of 18S rDNA gene.The amplification was carried out using the universal primers F04 and R22 (Blaxter et al. 1998).Subsequently, an index qPCR was performed to attach Illumina Unique Dual Nextera Indexes and compatible adapters to the amplicons of each sample.The loci qPCR was conducted on meiofauna and control/blank samples in a total volume of 20 μl, including 10 μl SsoAdvanced Universal SYBR Green Supermix from Bio-Rad, 0.5 μl of each primer (10 pM/µl), 2-9 μl genomic DNA (depends on the concentration), and 0-7 µl molecular grade water.The qPCR cycling settings consisted of an initial denaturation step at 98 °C for 2 minutes, followed by 25 cycles of denaturation at 98 °C for 15 seconds, annealing at 57 °C, and elongation at 72 °C for 1 minute.A two-step qPCR was then conducted to bind the Nextera XT Indexes and Illumina adapter overhang to the amplicons.The cycler settings included denaturation at 98 °C, annealing and elongation at 72 °C for 45 seconds for 10 cycles.Both qPCRs contained a melting temperature analysis (from 65-99 °C) in order to check the amplification success in addition to the Cq values and Relative Florescence Units (RFU).No-template controls have been included in both qPCR steps.After the second qPCR, the indexed amplicons were normalized according to the RFU of each sample, pooled and purified using AMPure XP Beads for PCR Purification (Beckmann Coulter) at a volume of 60%.The pooled library was denatured, and 15% genomic PhiX Reference DNA was spiked before performing the sequencing run using the MiSeq Reagent Nanokit v2 (250 cycles paired end) on an Illumina MiSeq platform at the DZMB Metabarcoding lab in Wilhelmshaven, Germany.
As for bioinformatics, the de-multiplexed NGS reads were trimmed from sequencing primers using BBmap (https://sourceforge.net/projects/bbmap/).Subsequently, the DADA2 pipeline (Callahan et al. 2016) was employed to denoise the amplicons.The sequences were further truncated to generate contigs, filtered based on length and quality scores, and checked for chimeras.The resulting high-resolution Amplicon Sequence Variants (ASVs) were then de-replicated.In order to assign taxonomical information to each ASV, a custom script SGN Metabarcoding pipeline (https://github.com/pmartinezarbizu)was used to blast the ASVs against the NCBI database incorporating the BLASTn pipeline.The taxonomic assignments of the ASVs were cross-checked against WoRMS -World Register of Marine Species (WoRMS Editorial Board 2024) and further literature sources to obtain the most up-to-date status.The target meiofauna ASVs were further classified into Operational Taxonomic Units (OTUs) using a maximum likelihood approach with a 3% cut-off threshold using the DECIPHER library package from the statistical software R (R Core Team 2023).

Metabarcoding
High-resolution meiofauna libraries obtained from the V1&V2 gene fragments were utilised for community analyses and graphical comparisons, employing various R packages such as DECIPHER, dada2pp, ggplot2 and vegan.
DECIPHER was used for aligning ASVs.The relative number of reads per higher taxon and core were visualised as bar plot using the R package dada2pp.Using the R package vegan (specaccum), rarefaction curves were calculated and plotted for the (1) number of OTUs per number of reads per core; (2) number of OTUs per number of samples per sampling date.
To compare the numbers of OTUs per core across the four different sampling dates a one-way Analysis of Variance (ANOVA, factor: Time) was performed using the statistical package SigmaPlot 14.This analysis was carried out for cores of the stations AH01-AH05.PAP was left out of this analysis because this site was not influenced by the sand nourishment.The numbers of OTUs per core per sampling date were visualised as a box-whisker plot (created in R using the ggplot2 library).
To determine the similarity between the meiofauna communities of single cores the number of reads per OTU per core (0-15 cm) were square-root transformed.Samples were compared using the Bray-Curtis Similarity.The results were visualised as non-metric Multidimensional Scaling (nMDS, created in R using vegan library).An Analysis of Similarities (ANOSIM; 99999 permutations) for the factor "Time" and subsequent pairwise tests were carried out using the PRIMER v7 software (Clarke and Gorley 2006).Another aim was to compare the results of the metabarcoding with our findings from the morphological observations (see below).Therefore, further subsets of the metabarcoding data were created through step-wise reduction of the number of considered samples.This involved omitting stations AH02 and AH04 (samples of these stations were not determined morphologically), and the use of only one "sample" per remaining station [median of the numbers of reads per OTU of 2 (T0) or 3 (T1-T3) samples per station, respectively].ANOSIMs were carried out and nMDSs were constructed to visualise the results (PRIMER v7).

Morphology
For meiofauna data based on morphological determination multivariate statistical analyses were carried out using PRIMER v7.An nMDS of communities of single cores was based on the Bray-Curtis Similarity of square-root transformed abundance data (0-15 cm).An ANOSIM (99999 permutations) was used to test for significant differences among the communities of the four sampling dates (factor: Time).Subsequent pairwise tests were carried out to test for differences among sampling dates.Another ANOSIM was applied to test for differences between the communities that were sampled in autumn (T0, T2) and spring (T1, T3; factor: Season).

Use of AI
DeepL.Write and the AI assistant from GoatChat were used to correct the English grammar and style.

Metabarcoding
The Illumina MiSeq sequencing run produced 1101158 reads, with 736978 identified as non-chimeric which were further classified into 3427 ASVs following DADA2 analysis encompassing 58 major groups according to the NCBI taxonomic assignments.Among these, 2322 ASVs were identified as meiofauna belonging to 11 groups containing animals, bivalves, crustaceans, flatworms, gastropods, gastrotrichs, nematodes, mites & ticks, rotifers, segmented worms and tardigrades.The meiofauna ASVs were clustered into 843 OTUs using 3% cut-off genetic distances.The OTUs were reviewed and further taxonomic filtering steps were performed to keep only true marine meiofauna OTUs from the target samples.The sample with the identifier A311 (T0, station AH03, 0-5-cm layer) was prepared but failed in the sequencing process.Therefore, the whole core was left out of the following analyses.
Planktonic Copepoda (order Cyclopoida) were excluded from the analyses.The here used group "Acari" includes the true marine mites of the family Halacaridae and species that belong to the family Hyadesiidae but resulted to be living in marine environments (Hyadesia aff.curassaviensis Viets, 1936) (Fain 1981).Marine Rotifera (Encentrum tectipes Remane, 1949) were kept in the data set.Gastropoda, Bivalvia, and Annelida were included in the analyses, although they usually are of macrofaunal size but their developmental stages may temporarily belong to the meiofauna.Additionally, the meiofauna-sized annelid genus Diurodrilus (Diurodrilidae Kristensen & Niilonen, 1982, Annelida incertae sedis) was detected in the morphological samples (see above) and included together with the other Annelida in the metabarcoding analyses.
Furthermore, OTUs were excluded from the analyses when they matched more than 97% with species that are described as fresh-water and/or terrestrial: Gastrotricha Chaetonotus (Hystricochaetonotus) borealis Kolicka, Kotwicki & Dabert, 2018(Kolicka et al. 2018); Platyhelminthes Microstomum lineare (Müller OF, 1773) (Müller 1774).Information on the habitat of the following excluded species was taken from the World Register of Marine Species (WoRMS Editorial Board 2024): Platyhelminthes Protomonotresis centrophora Reisinger, 1924; Nematoda Rhabdolaimus aquaticus de Man, 1880, Anaplectus porosus Allen & Noffsinger, 1968; Rotifera Lepadella ovalis (Müller, 1786).All taxonomical assignments of the OTUs were provided following the NCBI database and were checked against WoRMS and original literature.After performing quality and taxonomical data filtering, no sequences remained in the negative control and blank samples.Consequently, they were excluded from the downstream diversity analyses.
Based on rarefaction curves (Suppl.material 1: fig.S1), the sequencing reads were evenly distributed among samples (except for one sample at T1), reaching a sequencing depth plateau at approximately 7000 to 10000 reads per core.To avoid losing any diversity, no further rarefaction based on sequencing reads was performed prior to diversity analyses.
The curated OTU table and number of reads per OTU per sample are available from the ZENODO database (Kapshyna et al. 2024b).The unprocessed sequence data of this study can be accessed from the European Nucleotide Archive (ENA) using the project accession number PRJEB75380 (https://www.ebi.ac.uk/ena/ browser/view/PRJEB75380).The scripts for assigning taxonomy to sequencing reads and further post processing analyses in R were modified from original scripts available in the github repository https://github.com/pmartinezarbizu.

Metabarcoding: meiofauna community structure
Metabarcoding analyses showed clear differences in the composition of the meiofauna communities of the four sampling dates (Fig. 4).Before the sand nourishment (T0, September), more homogeneous communities at higher taxon level were detected in the area.Relative numbers of reads were more evenly distributed among higher taxa and communities were more diverse (in terms of presence of different higher taxa) compared to those after the impact (T1, March).Communities of the control station PAP mostly differed in composition from those of the impacted stations at the same date and varied over time.Apart from some reads, Acari were not found at PAP.At T0, Acari showed a total of 13075 reads and Annelida were represented by 16056 reads in samples from the impacted zone.At T1 after the impact both groups were almost absent from the impacted zone and Copepoda were highly reduced (T0 = 20567 reads in 9 samples, T1 = 9647 reads in 15 samples).On the other hand, Platyhelminthes were dominant (124320 reads) and relative numbers of reads of Gastrotricha had increased compared to the other meiofauna groups at T1 after the impact.Although the amount of DNA in each sample was normalised prior to sequencing during the NGS library preparation process, the absolute numbers of reads presented here still express the drastic community change induced by sand nourishment (for data see Kapshyna et al. 2024b).
Copepoda were re-established by T2 (September) but Acari were no longer detected during the remainder of the study (except for some reads during T2).Annelida were again detected at T2, but their relative numbers of reads did not reach the values measured before the sand nourishment.For Nematoda the relative numbers of reads were mostly higher at T3 (March) than at the other sampling dates (Fig. 4).

Number of operational taxonomic units
In order to visualise and compare the meiofauna richness per sampling core before and after the nourishment, a box-whisker plot was generated (Fig. 5).For each sampling date, the results from cores from stations AH01-AH05 were grouped into four categories: T0-T3.The samples collected at the control station PAP on each sampling date are represented as individual symbols (circles).The median number of OTUs per core in samples from the impacted zone was highest at T0 (n = 60) and ranged between 44 and 49 for sampling dates T1-T3.The number of OTUs per core in samples from the control station PAP (T1, T2) fell within the range of those from stations AH01-AH05 at the same time.Notably, at T0, samples from the control station PAP held less OTUs than those from AH01-AH05.By T3 the three PAP-data points were situated at the upper end of the box or above the 75% quartile of the data range of the other stations.
A one-way ANOVA (Factor "Time") of the numbers of OTUs per core of the impacted stations showed that there were no significant differences between the sampling dates (F = 2.778, p = 0.051; Suppl.material 1: table S1, Fig. 5).The differences in the mean values among the treatment groups were not great enough to exclude the possibility that they were caused by random sampling variability.Nevertheless, the standard deviations were high and the power of the conducted test (0.435) was below the desired power of 0.800 (indicating a reduced likelihood of detecting a difference when one actually exists).We therefore interpret this negative result with caution.
Rarefaction analysis of meiofauna OTUs across sampling dates (Fig. 6) revealed greater OTU richness as a function of sampling effort before the impact at T0 compared to T1-T3.Additionally, the number of detected meiofauna OTUs did not reach saturation at any of the sampling dates, suggesting that more samples are necessary to fully capture the meiofauna diversity in the studied area.

Morphological determination
The meiofauna (32-1000 µm; Fig. 1) was mostly represented by Copepoda, Nematoda, Platyhelminthes, Gastrotricha, and some Annelida.Furthermore, some rare taxa such as Rotifera, Halacaridae and Tardigrada were present with low abundances (Table 1; Kapshyna et al. 2024a).Only a few small individuals of Bivalvia and Gastropoda were found at T0 and T3.Annelida (apart from Diurodrilus sp. which is of meiofaunal size) were the only representatives of the macrofauna (>1000 µm) in our samples.However, members of Diurodrilus  could not be visually distinguished from Platyhelminthes during sorting.They were therefore counted as such.This "error" cannot have altered the morphological counts greatly because the metabarcoding results showed that the total number of reads for Diurodrilus over the whole study was only 4564, while Platyhelminthes had a total of 363455 reads (Platyhelminthes represented more than half of the total reads of all samples that were taken into account for this study).
We counted 27445 meiofauna individuals in total encompassing 10 higher taxa; 17 specimens could not be determined.The most abundant taxa were Copepoda (in total 30.5%) and Nematoda (in total 20.8%).Copepod nauplii and the combined group "Platyhelminthes+Diurodrilus sp." accounted for 20.8% and 17.0%, respectively.The rest together made up 11.9% of the total community.The composition of the taxa changed over time.At T0 (before the impact) and T2 (one year after T0) Copepoda and copepod nauplii dominated (≥ 33.8% and ≥ 27.1%, respectively).At T1 (directly after the impact) and T3 Nematoda (≥ 28,9%) and "Platyhelminthes+Diurodrilus sp." (≥ 32.9%) were the most abundant taxa (Table 1).While Platyhelminthes dominated in terms of abundance and number of reads at T1, Diurodrilus annelids were absent from the metabarcoding data at this date (except for only 9 reads in samples from PAP).Data of meiofauna counts per 5-cm sediment layer (0-5, 5-10 and 10-15 cm) are available from the ZENODO data base (Kapshyna et al. 2024a).

Community analyses
For the metabarcoding analyses (assessing numbers of reads per OTU per core), a non-metric Multi-Dimensional Scaling (nMDS) was created based on the Bray-Curtis Similarity of the square-root transformed data (Fig. 7A).The samples from the different sampling dates were clearly separated in the nMDS plot.Notably, due to the impact there was a substantial spatial difference observed between T0 and T1.At T2 the distance to the original state (T0) was less pronounced.Furthermore, from T0 to T2 samples from the control station PAP clustered centrally in the plot.A one-way Analysis of Similarities (ANOSIM) with the factor "Time" showed significant differences between the communities (global R = 0.712, p = 0.00001, Suppl.material 1: table S2A).Pairwise tests revealed significant differences among all sampling dates, but the greatest differences were detected between T0 and T1 (R = 0.838, p = 0.00001) and T1 and T3 (R = 0.926, p = 0.00001).Additionally, there was a significant difference between the tested seasons "Spring" (T1, T3) and "Autumn" (T0, T2): R = 0.229, p = 0.00002.
An nMDS was also created for the morphological data based on the Bray-Curtis Similarity of the square-root transformed abundance data (Fig. 7B).The samples from different sampling dates exhibited distinct clustering, resembling the pattern observed in the metabarcoding analysis (Fig. 7A).One-way ANOSIM with the factor "Time" showed significant differences between the communities (global R = 0.516, p = 0.0001; Suppl.material 1: table S2C).Samples collected at T0 were distinctly separated from the other sampling dates.Noteworthy is that the two T2 samples (AH01 and AH03) that are seemingly within the T0 data cluster (Fig. 7B) were actually situated on the same level as the other two T2 samples.This ordination is better visible in the 3D projection (Suppl.material 1: fig.S2).Pairwise tests indicated the greatest differences between T0 and T1 and T0 and T3 (for both pairs R = 0.896, p = 0.029).The difference between T1 and T2 was less pronounced (R = 0.510, p = 0.029).The pairs T0 vs T2 and T1 vs T3 were just at the border of not being significantly different (p = 0.057) which may be attributed to the limited number of samples analysed in each category.The samples taken at T2 were not significantly different from those of T3 (Suppl.material 1: table S2C).A second ANOSIM detected significant differences between the tested seasons "Spring" (T1, T3) and "Autumn" (T0, T2): R = 0.413, p = 0.003.
Bubble plots superimposed on the nMDS of the morphological data depict the abundances of the key taxa including Nematoda, Copepoda, the combined group "Platyhelminthes+Diurodrilus sp.", copepod nauplii, Gastrotricha, and Acari (Fig. 8).Copepoda were the most abundant taxon at T0 (max 804.7 ind.10 cm −2 at AH03), but they were almost absent from the impacted zone at T1 (max 22.8 ind.10 cm −2 at AH03).The trend of copepod nauplii followed that of Copepoda.In contrast, "Platyhelminthes+Diurodrilus sp." exhibited an inverse pattern with the lowest abundances in the impacted zone at T0 (max 37.1 ind.10 cm −2 at AH05) and increased counts at T1 (max 351.1 ind.10 cm −2 at AH05).The group peaked at T3 (max 417.9 ind.10 cm −2 at AH03).Abundances of Gastrotricha had decreased by T1 but recovered by T2 and T3.Although Acari were scarce, they were most abundant at T0 and almost absent thereafter, except for stations AH05 and PAP at T2.Nematoda did not show a consistent trend during the course of the study (Fig. 8, Table 1).

Comparison of outcomes
For the comparisons of morphological and metabarcoding results, the metabarcoding dataset was stepwise reduced to match the number of samples that were visually inspected and morphologically determined.The outcomes of all ANOSIM analyses conducted in this study are presented in Suppl.material 1: table S2.Non-metric multidimensional scaling (nMDS; 2D stress = 0.09) of communities of single cores collected at the sampling dates T0, T1, T2, and T3 based on the Bray-Curtis similarity of square-root transformed abundance data (core depth 0-15 cm).Superimposed bubble plots show the number of individuals per 10 cm 2 for the meiofauna taxa Nematoda, Copepoda, Gastrotricha and Acari.Results for the combined group "Platyhelmin-thes+Diurodrilus sp." and copepod nauplii are also shown.
In order to align with the morphologically investigated stations, data from stations AH02 and AH04 were omitted from the metabarcoding dataset.The adjustment resulted in a decreased global R value of 0.611 (p = 0.00001, Suppl.material 1: table S2B.1), yet the nMDS plot still showed distinct separations between sampling dates (Suppl.material 1: fig.S3).Ranking the pairwise tests based on their R-value revealed a slightly different order compared to the complete data set, with the highest R observed in the T1 vs T3 pair, followed by T0 vs T1 and then T1 vs T2.
Subsequently, the number of metabarcoding cores was reduced to match that of the visually inspected cores.The medians of the numbers of reads per OTUs were calculated for each remaining station and sampling date.Remarkably, despite these adjustments, distinct group separations persisted (Suppl.material 1: fig.S4) with a highly significant global R value of 0.563 (p = 0.00003; Suppl.material 1: table S2B.2), comparable to the morphological determination outcomes (Fig. 7A, Suppl.material 1: table S2C).However, the ranking of pairwise tests' results underwent another change in order and exhibited lower significance values similar to the morphological results (Suppl.material 1: table S2C) due to reduced sample numbers.

Discussion
In the course of this study, we applied two different methods to assess the impact of sand nourishment on meiofauna communities at the beach-water interface of a touristic beach on the Baltic Sea coast in Ahrenshoop: the contemporary high-throughput approach of whole community metabarcoding and the conventional method of morphological determination and counting of the animals.Here, we compare the findings and offer recommendations for efficiently assessing the impact of disturbance in sandy-beach ecosystems in the future.

Anthropogenic impact monitoring using meiofauna as reference organisms
Meiofauna organisms are well-suited for studying environmental impacts.Because of their small body size, ubiquity, and high abundances their communities can be analysed using small sample sizes and high numbers of replicates.Due to their short generation times, the impact on and recovery of communities can be observed within manageable time frames (Zeppilli et al. 2015;Schratzberger and Ingels 2018).
Benthic communities, including meiofauna, are highly impacted by sand nourishment activities.Nourishments had sizeable impacts on the components of several beach ecosystems, including benthic organisms (review by Speybroeck et al. 2006).Many studies showed that abundances decreased after sand nourishments (e.g., Peterson et al. 2000;Menn et al. 2003;Schlacher et al. 2012).These findings are in line with our results that were obtained directly after the sand nourishment.We observed a sharp decline in total abundance and changes in community composition including a drastic decline in Copepoda and the almost complete disappearance of some groups such as Acari and Annelida.However, meiofauna abundance and richness can even increase after an impact (Fegley et al. 2020).
Van Dalfsen and Essink (2001;North Sea) found that it may take four years for benthic communities to fully recover after a sand extraction.They noted that recovery of the total abundance took place within two years (except for crustaceans), but it was not sufficient for the biomass and diversity of the communities to re-establish completely.This is consistent with the findings by Menn et al. (2003) who recommend at least a 3-year interval between nourishment activities (North Sea).A full recovery of abundances of some specific species may take place faster.It took several Polychaeta species only one year to recover after a beach nourishment (Leewis et al. 2012;North Sea).And there are a number of studies that reported faster, but incomplete recovery regarding only abundance or only diversity.Schlacher et al. (2012;Gold Coast on the east coast of Australia) noticed that the species richness, but not the abundances of macrofauna, recovered after five months.The process of recovery depends on different factors, such as specific life-cycle characteristics of species (Schlacher et al. 2012), sediment mismatches (Peterson et al. 2000), or the volume of the nourished sediment (Menn et al. 2003;Herman et al. 2021).
Our findings revealed a significant impact of sand nourishment on meiofauna communities and a substantial trend towards recovery within six months post-impact.The initial notable difference was reduced (Figs 7, 8) and the composition of meiofauna communities developed towards their pre-impact state (Fig. 4, Table 1).However, the communities of the different sampling dates remained distinct until the last sampling one year after the nourishment (Suppl.material 1: table S2) and never fully reverted to the composition observed at T0 (Table 1).This observation may be attributed to the relatively short monitoring period, the substantial quantity of nourished sediment, or a potential sediment mismatch caused by grain size variations due to the nourishment (as indicated for the Ahrenshoop event by Glueck 2023).Increases in the quantity and heterogeneity of coarser sediment are known to positively affect post-impact meiofauna abundances and richness (Fegley et al. 2020).However, we noted a more uniform sediment composition after the nourishment.Results on grainsize composition will be detailed as part of an upcoming publication.
Organisms of the meiofauna size class were the only noteworthy metazoans down to 15 cm depth at the beach-water interface of the Ahrenshoop beach.Recovery of meiofauna can be influenced by passive resuspension and active migration processes (Menn et al. 2003).Communities may re-establish through resuspension caused by wave action occurring elsewhere (especially after storm events) followed by a current-driven transport.Thus, new organisms (adults and developmental stages) can reach the Ahrenshoop beach by coast-parallel transport from the south-west to the north-east along the impacted zone.Organisms residing in the swash zone demonstrate high adaptability to regular resuspension.They can actively return to the sediment, may float or engage in active swimming in the water column for a certain time before re-entering the sediment.Once settled, reproduction is another decisive factor for a successful re-colonisation.
Seasonal effects cannot be excluded for certain taxa.For Platyhelminthes the highest counts and numbers of reads were observed in part of the spring samples from T1 and T3.This is in line with Dittmann and Reise (1985) who found mass developments of certain platyhelminth species on a North Sea tidal flat as early as February.Platyhelminth abundances in the autumn samples from Ahrenshoop were generally low, regardless of the time since impact.The replenishment started at the end of October 2021 at the southern end of the beach section and was finished in the middle of February 2022 at the northern end.Thus, at station AH05 (southern end), four months had passed between the nourishment and T1 sampling, whereas at station AH01 (northern end), the post-impact period was one month.In the following, we therefore focus on the comparison of the two methods, irrespective of the observed recovery and possible seasonal effects.

Impact assessment approaches: Metabarcoding versus morphology
Meiofauna are frequently neglected in the context of monitoring marine ecosystems although they may be used as key indicator organisms for ecosystem health (e.g., Baguley et al. 2015;Amorri et al. 2022).Their high diversity is challenging and many taxa lack easy-to-detect distinctive morphological features.Therefore, specialised expertise is mandatory for their determination.In the past decade, the emergence of rapid monitoring techniques such as metabarcoding has facilitated the study of natural and anthropogenic shifts in meiofauna populations.
Several marine studies employed both microscopic and metabarcoding approaches, and highlighted discrepancies in results between the two methods.Mostly, morphological determinations revealed a lower diversity than genetic analyses (e.g., Cowart et al. 2015;Lejzerowicz et al. 2015;Ammon et al. 2018;Cahill et al. 2018).Even juvenile stages, damaged specimens, and taxonomic groups susceptible to preservation and extraction processes are still detectable by metabarcoding and may be determined to species level.However, sometimes morphological identification of meiofauna organisms can reveal more taxa than metabarcoding (Pantó et al. 2021).Issues such as the inefficiency of universal primers, PCR failure rates, and the lack of reference barcodes for meiofauna species in publicly available data repositories have been identified as key challenges (e.g., Lejzerowicz et al. 2015;Pantó et al. 2021).Some studies strongly recommend metabarcoding as the preferred monitoring method (Lejzerowicz et al. 2015;Martínez et al. 2020), but the general consensus is to combine both methods to complement each other (Ammon et al. 2018;Cahill et al. 2018).
In our investigation, we took a step back and compared metabarcoding at the OTU level with morphological determinations at the higher taxon level.Both approaches showed similar results: distinct clustering patterns associated with the four sampling dates and significant differences in community composition between pre-and post-impact (Figs 4, 7, 8; Suppl.material 1: table S2).The number of OTUs (Fig. 5) and the number of higher taxa observed microscopically were not significantly affected by the impact and remained the same during the recovery phase.It is known that the morphological approach at the higher taxon level can adequately explain the ecological aspects of meiofaunal communities.A large number of studies conducted in different marine environments have shown that the structuring effects of biotic and abiotic environmental parameters can be addressed by morphological determination and counts at higher taxonomic levels (e.g., Kuhnert et al. 2010;Sajan et al. 2010;Sedano et al. 2014;Säring et al. 2022).However, even at this level, meticulous work is required to detect and accurately identify the smallest and rarest meiofauna (e.g., Tardigrada), which are often crucial for differentiating communities (Säring et al. 2022).
The metabarcoding approach is particularly suited to cases where the time-consuming microscopic detection and determination of soft-bodied taxa such as Platyhelminthes, small Annelida (e.g., members of the genus Diurodrilus) and rare and small groups such as Gastrotricha and Tardigrada are required.These taxonomic groups are also often sensitive to commonly used preservatives, such as formalin, which can lead to degradation of specimens (Higgins and Thiel 1988).Especially in Platyhelminthes, this preservation effect is a frequently cited source of strong discrepancies between morphological and metabarcoding results (Lejzerowicz et al. 2015;Rzeznik-Orignac et al. 2017).We took this factor into account in our study, fixed all samples in ethanol and found that Platyhelminthes were well preserved and recognisable.
Given the faster and more cost-effective nature of metabarcoding compared to traditional morphological methods (Rossel et al. 2019;Martínez et al. 2020), our study findings indicate that metabarcoding is a viable option for future monitoring of Baltic Sea sand-nourishment projects.However, it is crucial to analyse a significant number of metabarcoding samples.While prioritising metabarcoding, it is important to retain a subset of morphological samples for validation purposes, establishment of a reference barcode library through single specimen barcoding, taxonomic classification, and future quantification of metabarcoding outcomes.We recommended that appropriate sampling techniques, such as the use of ethanol instead of formalin, consistent sieve sizes and standardised extraction methods, are employed to ensure optimal alignment of results.

Conclusions
-High-throughput methods are useful to study impacts on meiofaunal communities and to track their recovery.Our study was carried out on a sandy beach with a meiofaunal community of low diversity.This may have facilitated the fact that the metabarcoding results followed the same trend as those obtained by traditional morphological determination.-Using both traditional morphological and high throughput genetic methods, our results were consistent and provide a broad basis for interpretation.
The results demonstrate the negative impact of sand nourishment on meiofaunal communities and a long recovery time.-After initial studies using both methods, metabarcoding alone may be considered for comparable sites and situations.

Figure 1 .
Figure 1.Meiofauna organisms collected at the beach-sea interface of Ahrenshoop (Baltic Sea).Scans taken by confocal laser scanning microscopy (CLSM), organisms stained with acid fuchsin and Congo red.Top left to right: one individual each from Tardigrada and Platyhelminthes (upper left scale bar), Copepoda Harpacticoida and Nematoda (upper right scale bar).Bottom left to right: Halacaridae, Annelida (bottom left scale bar).All scale bars = 100 µm.Scans taken by Iryna Kapshyna, Olena Uzun, Tobias Fischer.

Figure 2 .
Figure 2. Monitoring the impact of sand nourishment on interstitial meiofauna communities along the sandy beach of Ahrenshoop (2021-2023): Five sampling stations (red circles) lay within the impacted zone (AH01-AH05), one unaffected reference station (PAP) was located south of Ahrenshoop.Inlay map: Overview of the western Baltic Sea with arrow indicating the location of the beach on the Fischland-Darß-Zingst Peninsula.

Figure 4 .
Figure 4. Metabarcoding of meiofauna communities collected at the water line before and after sand nourishment at the Ahrenshoop beach (2021-2023).Bar plot of the relative numbers of reads per higher taxon per core on the sampling dates T0, T1, T2, and T3 (0-15 cm).Samples of the five stations AH01-AH05 lay within the impacted zone, samples of station PAP south of the nourished beach.

Figure 5 .
Figure 5. Number of meiofauna Operational Taxonomic Units (OTUs) per core sample (0-15 cm) before and after sand nourishment in Ahrenshoop (Baltic Sea, 2021-2023).Per sampling date, the values of the cores from stations AH01-AH05 are used to create the boxes (including median and upper and lower percentiles: 75%, 25% (box frame), 95% and 5% (whiskers).For each sampling date, the values from the control station PAP are shown as circles.

Figure 6 .
Figure 6.Rarefaction curves of meiofauna OTUs per core before and after sand nourishment in Ahrenshoop (Baltic Sea, 2021-2023).The values of cores from stations AH01-AH05 and PAP are used for each sampling date (T0-T3).Species accumulation curves and confidence polygons (unconditional standard deviation) were created using the R package vegan (specaccum).

Figure 7 .
Figure 7. Monitoring meiofauna before and after sand nourishment at the Ahrenshoop beach (2021-2023).Multi-dimensional scaling of similarity analyses based on results from single cores.A. Metabarcoding: Bray-Curtis Similarity of the number of reads per higher taxon (square-root transformed; ANOSIM global R = 0.712, p = 0.00001); B. Morphological identification: Bray-Curtis similarity of the individual numbers per higher taxon (square-root transformed abundance data, 2D stress 0.09; ANOSIM global R = 0.516, p = 0.0001).

Figure 8 .
Figure 8. Morphological counts of meiofauna communities of the water line (beach-water interface) before and after sand nourishment in Ahrenshoop (Baltic Sea, 2021-2023).Non-metric multidimensional scaling (nMDS; 2D stress = 0.09) of communities of single cores collected at the sampling dates T0, T1, T2, and T3 based on the Bray-Curtis similarity of square-root transformed abundance data (core depth 0-15 cm).Superimposed bubble plots show the number of individuals per 10 cm 2 for the meiofauna taxa Nematoda, Copepoda, Gastrotricha and Acari.Results for the combined group "Platyhelmin-thes+Diurodrilus sp." and copepod nauplii are also shown.