Benthic species distribution linked to morphological features of a barred coast

: The composition of benthic species communities in the nearshore zone is closely related to the hydrodynamic and morphodynamic conditions. Sustainable management of the coastal ecosystem requires knowledge about the natural dynamics as well as human-induced changes on the ecosystem. To improve our knowledge of the benthic species distribution along a dissipative sandy shore with multiple breaker bars, an extensive dataset was collected in the nearshore zone of the barrier islands Ameland and Schiermonnikoog in the Dutch North Sea. From 2010 to 2014, every year, approximately 180 grab samples along 18 cross-shore transects were collected and analyzed for sediment characteristics and macrobenthic species composition. Mixed-e ﬀ ect-models and partial redundancy analysis were used to analyze the importance of morphological features (i.e., slopes, bar crests, and troughs) as an explanatory variable for the benthic species distribution. The results indicate that the morphological features in themselves explain three times more variation than the environmental parameters used. This demonstrates the importance of morphological features as a factor in explaining the distribution of benthic species communities in the nearshore. Detailed information on morphological features is easy to obtain from bathymetry maps or visual inspection. Incorporating morphological features in species distribution models will therefore help to improve sustainable management of our valuable sandy coastal systems.


Introduction
Sandy beaches are one of the world's dominant coastal types [1]. They constitute important habitats for benthic species, fish, and birds and provide a multitude of ecosystem services ranging from flood safety to food provision and recreation [2]. Beach and foreshore erosion is a common problem in the world, that may aggravate due to climate change and sea level rise [3]. In many places, sand nourishments are applied to mitigate the effects of erosion [4,5], and these efforts are expected to increase in the future. Nourishments alter the natural morphodynamics of the coast [6,7] and may thereby affect the quality of nearshore habitats. This is especially relevant for macrobenthic species that live close to and in the bed.
Macrobenthic species play a central role in the ecosystem of sandy coastal systems. They serve as food for fish and birds, forming an important link in the food web between the primary production distribution of the benthic species and species communities in the nearshore, above and besides the impact of the general slope towards deeper water.
We investigated this hypothesis with extensive field surveys on two dissipative barred shores along the North Sea coast of the barrier islands Ameland and Schiermonnikoog, the Netherlands. From 2010 to 2014, each year approximately 180 samples of the seabed were collected along 18 cross-shore transects and analyzed for macrobenthic species composition and sediment characteristics. This resulted in an extensive dataset giving us the opportunity to relate the benthic species distribution to the morphological features present on a dissipative barred shore.

Study Area and Data Collection
Data were collected in the Netherlands along the North Sea coasts of Ameland island (53 • 28 N, 5 • 51 E) and Schiermonnikoog island (53 • 31 N, 6 • 16 E) ( Figure 1). Both study sites have a dissipative beach with a barred foreshore with two to three straight continuous bars in water depths between 5 m and 8 m. This profile with two or three breaker bars is common along the entire coast of the Netherlands [34].

Study Sites
The study area at Ameland is located on the eastern part of the island between beach post 17.2 and 23.4. Here the nearshore bed has a general slope of 1:275 with a breaker bar around 4 m water depth and a second breaker bar around 6m water depth. The depth at the most seaward locations, 3 km offshore, ranged between 11 m to 15 m below mean sea level (MSL). The nearshore of Ameland was nourished at water depths between 5 m and 8 m in the years 1998, 2003, 2006, and 2010-2011. The beach of Ameland was nourished in 1992, 1996, 2006, and 2010-2011. The second study area is located at Schiermonnikoog between beach post 11.6 and 13.8. Here the nearshore has a very gentle slope (1:350) and a double bar-system at 2.5 m to 4 m water depth. Compared to Ameland, the entire bed profile and all morphological features present are shallower. The deepest location was at approximately 8 m below MSL.
Both sites have a similar environmental setting, such as a west-east orientation of the coastline and a semidiurnal tide propagating from west to east with a tidal amplitude of approximately 2 m. The dominant wind direction is from the southwest, but large storm events are frequently associated with northwesterly winds. The dominant wave direction is northwest as southwesterly wind and waves are shielded by the island. The mean annual wave height is approximately 1.1 m. The sediment in both areas consists of fine sand [36].

Data Collection Techniques
For five years, samples were collected once per year. Sampling took place in late summer or early fall (August 2010, September 2011, October 2012, September 2013, and September 2014) to reduce seasonal variation. Moreover, this enabled us to sample fully grown specimens that are able to survive the next winter as well as to reduce the bias in abundancy introduced by spatfall during spring and summer. In total 766 benthic samples and 294 sediment samples were collected. For the year 2011, locations under the influence of the nourishment were excluded from the analysis. Before each campaign, sample locations were chosen along transects perpendicular to the coast (Figure 1), guided by morphological features as the crests and troughs of the bars and the seaward and landward slopes present ( Figure 2). The precise position of the morphological features was obtained from recent annual bathymetry measurements (JARKUS) of the Dutch Ministry of Infrastructure and Water Management, or if not available, the bathymetry of the bed was recorded using a single beam echo sounder with a frequency of 210 kHz. During the campaign, the position and depth of the sample locations were logged by the ships' navigation system.
At each location, a sample of the seafloor was taken using a Van Veen grab with a sampling surface of 0.1 m 2 . Only samples with a minimum penetrating depth of 10 cm were accepted. Otherwise, the sample was re-taken. Samples were wet-sieved using a 1 mm mesh sieve. The residue was stored in jars and preserved with a seawater solution of 6% buffered formaldehyde. The macrobenthic samples were sorted in the laboratory and species were determined under a binocular microscope using the most recent literature and Dutch standard nomenclature (TWN) [37][38][39][40][41]. Specimens were identified up to species level when possible. For each species, the number of individuals per sample was recorded. Infauna ash-free dry weight biomass (g AFDW/m 2 ) was analyzed using the "Blotted Wet Weight" method and conversion factors [42]. Sediment samples were taken of the top 5 cm from the Van Veen grab before the sieving and kept frozen until analysis. Organic matter was analyzed for 273 locations by loss on ignition at 550 • C. After the removal of organic matter and carbonate, the sediment grain size distribution was obtained using a Malvern Mastersizer 2000.

Data Analysis and Methods
The relationship between benthic species composition and morphologic features (MF) in the nearshore was investigated using univariate and multivariate analysis. Univariate analyses tested for effects of morphometric properties on sediment composition and univariate community properties total density, total biomass, number of species per sample, and species diversity. With multivariate regression the relationship between the morphological features and the benthic species composition was investigated. All analyses were done in R version 3.5.2 [43].

Definition of Morphological Features
Based on the bathymetric information we divided the nearshore into bar crests (B), troughs (Tr), and slopes (S). We distinguished four seaward slopes (index S1-S4) and one landward slope (index L). This resulted in 11 morphological features, illustrated in Figure 2 and defined in Table 1. At times the beach slope (S B ) was bordered at its seaward end by a bar-like elevation creating the beach trough (Tr B ). The transition between slope and trough or slope and bar crest were discretized in such a way, that 25% of the length of the slope was allocated to the nearest trough, and 25% of the length of the slope was allocated to the nearest bar crest. These boundaries account for the spread in sampling locations within each morphological feature. In seaward direction: (S B ) beach slope, (Tr B ) beach trough, (Tr 1 ) first trough, (B 1 ) first bar crest, (S S1 ) seaward slope first bar crest, (Tr 2 ) second trough, (S L ) landward slope second trough, (B 2 ) second bar crest, (S S2 ) seaward slope second bar crest, (S S3 ) seaward slope, (S S4 ) gentle slope.

Characterization of the Morphological Features
The morphological features were characterized using five morphometric properties: (1) the aggregated bed level along a transect (profH), representing the overall increase of water depth with distance to the shoreline. The aggregated bed level was determined by the least-squares fit of the logarithmic relation between bed level and distance to the low water line. (2) The relative height (relH), defined as the difference between the actual bed level and the aggregated bed level. To analyze the relationship between morphological features and their morphometric characteristics and sediment properties, we used mixed effect models (ME-models). ME-models include both fixed and random effects of the predictor values on the predictand value and enable accounting for the nested structure of the dataset (samples along transects over multiple years) [44,45]. The lme-routine of the R-package nlme [46] was used to relate the morphological features and their characteristics (fixed effects X ijk ) to the response variables Y ijk (here D50, mud percentage and OM). Random effects (b i + b ij + w i ε ijk ) were used to explain variation associated with the temporal and spatial structure of the dataset. Spatial variation (Transect) was nested within sampling year (Year) allowing both parameters to have a different mean value. Heterogeneity over the years was present and included by varying the residual variance associated with the year of the survey. This resulted in the following model structure.
where b i ∼ N 0, σ 2 1 is the variation associated with the year of sampling, b ij ∼ N 0, σ 2 2 is the variation associated with transect within a year, and w i ε ijk ∼ N 0, σ 2 3 is the heterogeneity varying with each year. The years of sampling (2010-2014) are represented by i = 1...5, the cross-shore transects are indicated with j = 1...18 and the sampling position along the transect is indicated with k = 1...11. The fixed component exists of a vector (X ijk ), containing the selected explaining variables (predictors), and the estimated slope vector (β), which contains the estimated contribution of the predictors to the value of the response variable (predictand). For each ME-model, the categorical variable morphological feature (MF) and the continuous variables relative height (relH), aggregated bed level (profH), slope, orientation (aspect) and sedimentation-erosion (SE) were considered as possible fixed variables. Forward selection was used to select the optimal variables [47,48]. Before model estimation, normality of each variable was evaluated using QQ-plots and when necessary the data was transformed. For the ME-models with MF selected as an explanatory variable, a post-hoc pairwise comparison with least-square means was used (emmeans-package [49]) to identify which morphological features were significantly different from each other for the specific response variable under investigation.

Benthic Species Data
Species were identified at species level if possible. Ringworms (Nemertea) were not further specified into species. Rare species with minor abundances were grouped at the genus level. Species found in less than 20 samples (approximately 2.5% of the total number of samples) were excluded. Large or fast-moving epibenthic species (Decapoda and Mysida) are not properly sampled with a Van Veen grab and were excluded from further analysis. Finally, hard substrate species from the class Hydrozoa and Bryozoa were also excluded from the analysis as only their presence in a sample was recorded. At three locations the samples were not properly conserved and therefore removed from the dataset.
For each sample location, total abundance, total biomass, number of species (S), and species diversity expressed as the Shannon diversity index (H') were determined. Shannon's index accounts for both abundance and evenness of the species present at a location. Again, ME-models were used to relate the morphological features and their characteristics (X ijk ) to the response variables Y ijk (here total abundance, total biomass, S and H'). The selection procedure of the explanatory variables, the random structure of the ME-models and the final step using post-hoc pairwise comparison to identify which morphological features were significantly different from each other for the specific response variables under investigation, was kept the same as for the ME-models set up for the abiotic response variables (see Section 3.2).

Relationship between Benthic Species and Morphologic Features
To reveal the relationship of the benthic fauna with morphological features, sediment characteristics and the other morphometric parameters, partial redundancy analysis (pRDA) was used. Redundancy analysis is a multivariate method to relate species composition to several underlying environmental parameters [47], comparable to the univariate ME-models used to investigate the characteristics of the morphological features. As an RDA is expecting linear relationship preserving the Euclidean distance among sites, it is not naturally adapted to the analysis of species abundance data. First, a log-transformation was applied to weigh down the effect of large abundances, followed by a Chord transformation [50]. The latter was used to account for the fact that certain species are missing in both compared community samples, also known as the 'double-zero-problem'. Partial ordinations (pRDA) corresponds to partial regression. The influence of covariables is first removed and the explanatory power of the remaining environmental variables is tested [48,51].
We were interested in the overall effect of morphological features on benthic species distribution. Therefore, the longshore variation along the coast of Ameland and Schiermonnikoog represented by the transects and the temporal variation represented by sampling-year were included as covariables. In the nearshore, the overall increasing water depth represented by the aggregated bed profile has a very prominent effect on the species distribution and is largely entangled with the morphological features [23,35]. Therefore, the aggregated bed profile for each transect was also included as a covariable. A forward selection was used to derive the optimal combination of environmental variables and covariables according to the amount of variance in the species data they captured [52]. The statistical significance of the final model was tested using a Monte Carlo permutation test (1000 unrestricted permutations).
In a final step, variation decomposition [51] with adjusted R 2 values [53] was used to partition the variation explained represented by the pRDA into three compartments: (i) the effect of the morphological features, (ii) the effect of the selected environmental parameters, (iii) the effect of the covariables.
As the sediment characteristics were not determined for each grab sample taken in the field, the D50, mud content and OM were predicted based on the ME-models described under Section 3.2.
To identify the representative species for a specific morphological feature, an indicator value was calculated for each species using the 'IndVal-routine' of the package labdsv [54]. This indicator value was derived by taking the product of the relative frequency and relative abundance of a species in a morphological feature. A high value represents a combination of specificity, i.e., a large mean abundance of a species within a morphological feature compared to the other morphological features, and fidelity, i.e., the presence of the species in most locations of that morphological feature. Values over 0.05 were used to select the indicator species. There is no absolute (e.g., significance-based) criterion for setting a threshold. In our study, a threshold of 0.05 was high enough to provide a manageable number of indicator species in the presentation of the results, while it was low enough to ensure inclusion of most species that show some specificity and fidelity.

Shape Characteristics of the Morphological Features
The elevation of the bed in the study area varied between MSL and 11 m below MSL ( Figure 3A). A clear double-bar system was present with bar crests at median 3.6 m (IQR = 2.1-3.9) (B 1 ) and 5.5 m (IQR = 3.5-5.7) (B 2 ) water depth. The median bed level of the troughs was 3.8 m (IQR = 3.5-4.5) (beach trough: Tr B ), 5.2 m (IQR = 2.6-5.6) (first trough: Tr 1 ), and 6.3 m (IQR = 4.3-6.6) (second trough: Tr 2 ) below MSL. The local slope of the bed was gentle and varied in the cross-shore direction between 0.2% and 4.1% ( Figure 3B). The slopes were mostly oriented in northern direction (seaward slopes S B , S S1 , S S2 , and S S3 ) i.e., seaward dipping, except for the beach slope (S B ) with some locations directed in a southern direction. This possibly represents an early stage of new nearshore bar formation close to the shoreline. The landward slope (S L ) was directed in a southern direction. The slope orientation of the bar crests and troughs (B 1 , B 2 , Tr B , Tr 1 , and Tr 2 ) was a mixture of northward and southward orientations ( Figure 3C), in line with the fact that these features are centered around local minima and maxima in the topography.
The nearshore is a very dynamic area subject to the action of waves and tidal currents causing both sedimentation, erosion and bar migration ( Figure 3D). This is reflected in the net annual bed level changes up to one meter in the morphological features landward of the second breaker bar. Seaward of the second breaker bar the net yearly sedimentation-erosion was almost zero.

Sediment Characteristics Related to the Morphological Features
The median grain size, organic matter, and mud percentage varied in cross-shore direction within small ranges. The sediment is characterized as medium-fine sand (150-210 µm) with an overall D50 of 172 µm (IQR = 164-184). Despite this small range, there was a seaward gradient with coarser

Sediment Characteristics Related to the Morphological Features
The median grain size, organic matter, and mud percentage varied in cross-shore direction within small ranges. The sediment is characterized as medium-fine sand (150-210 µm) with an overall D50 of 172 µm (IQR = 164-184). Despite this small range, there was a seaward gradient with coarser sediment at the beach slope and finer sediments at the foot and downward slope of the outer beaker bar ( Figure 4A). Model selection resulted in a ME-model for the median grain size (ME D50 ) with, aggregated bed level (profH) and relative height (relH) in combination with morphological features (MF) as explanatory variables (Table 2). A likelihood ratio test between the ME D50 model and the ME D50 model without the morphological features showed that the ME D50 model with morphological features is the better model (ME D50 AIC = 2073, ME D50 no-MF : AIC = 2190, with p < 0.0001). The post-hoc pairwise comparison over the morphological features of the ME D50 model showed that the median grain size gradually changed over the morphological features (see letters in Figure 4A).  Table 2. Model estimates for selected explanatory variables (fixed effects) and variance (random effects) for the ME-models predicting median grain size (ME D50 ), percentage mud (ME Mud ), and organic matter (ME OM ). The mud percentage was very low throughout the study area, with an overall median of 0.8% (IQR = 0-1.5) ( Figure 4B). Over the morphological features the median mud percentage varied between 0.5% (IQR = 0-0.7) at the beach slope and 1.5% (IQR = 1.1-4.2) at the most seaward slope. The optimal ME-model for mud percentage (ME mud ) included relative bed level (relH), aggregated bed level (profH) and year of survey (Year) as explanatory variables (Table 2). Morphological features did not contribute to the ME mud model and were not included as explanatory variables (ME Mud AIC = 824, ME Mud with MF AIC = 825, p = 0.04) Therefore, no post hoc test over the morphological features was performed.

ME
The median percentage of organic matter found was 0.5% (IQR = 0.3-0.7). The minimum and maximum values ranged between 0.2% at the beach slope (S B ) and 1.9% in the first trough (Tr 1 ). There was an increase in organic matter in the seaward direction (0.3% (IQR = 0.3-0.4) -0.6% (IQR = 0.4-0.8) ( Figure 4C). The optimal ME-model for organic matter (ME OM ) had aggregated bed level (profH), relative bed level (relH) and year (Year) as explanatory variables (Table 2). No post hoc test for the morphological features was performed as the morphological features did not modulate the effect of the relative height and aggregated bed profile on organic matter (ME OM-MF AIC = −128.8, ME OM AIC = −127.2, p = 0.0176).

Macrobenthic Species in the Near Shore
In total 136 unique taxa were found in the nearshore of Ameland and Schiermonnikoog. After processing the raw data 50 taxa from 5 phyla were used in the analysis. Most taxa (n = 25) were polychaetes, followed by crustaceans (n = 13 species) and bivalves (n = 9 species). Furthermore, two echinoderm species were found and finally nemertean species which were not further specified into species. In total, 40 taxa were identified at species level and 10 at genus, family, or higher level.

Total Abundance, Biomass, Number of Species, and Diversity of the Benthic Fauna
Across the nearshore, abundance varied between 50 to 2110 individuals/m 2 at the beach slope locations (SB) up to 120 to 8370 individuals/m 2 at the first seaward slope (SS1). The median total

Total Abundance, Biomass, Number of Species, and Diversity of the Benthic Fauna
Across the nearshore, abundance varied between 50 to 2110 individuals/m 2 at the beach slope locations (S B ) up to 120 to 8370 individuals/m 2 at the first seaward slope (S S1 ). The median total abundance per morphological feature was highest at the first slope (S s1 = 1655 ind/m 2 (IQR = 517-3695) and lowest at the beach slope and first breaker bar (S B = 590 ind/m 2 (IQR = 335-1804), B 1 = 650 ind/m 2 (IQR = 365-1390) ( Figure 6A). Total abundance was best described with a ME-model (ME AB )with aggregated bed level (profH), relative height (relH), mud content, slope, and morphological features as explanatory variables. The morphological features contributed considerably to the ME AB model (ME AB AIC = 1692, ME AB no-MF AIC = 1756, p < 0.001) ( Table 3).
Median total biomass was highest at the foot of the second breaker bar and further offshore (S S3 = 29.9 g/m 2 (IQR = 8.7-51.3) and S S4 = 17.0 g/m 2 (IQR = 8. 8-30.8) and lowest at the beach slope and first breaker bar (S B = 0.4 g/m 2 (IQR = 0.2-0.8) and B 1 = 1.8 g/m 2 (IQR = 0.5-5.0)) ( Figure 6B). The ME-model with the explanatory variables aggregated bed profile (profH), relative height (relH), mud percentage and morphological features as fixed effects described the total biomass best (ME BM ) ( Table 3). Model comparison between the ME BM model with and without morphological features showed the significance of the morphological features (ME BM AIC = 2398, ME BM no-MF AIC = 2496, p < 0.0001). The pattern for biomass is almost equal to that of the total abundance with the seaward slope of the first breaker bar (S S1 ) and the most offshore location (S S4 ) considerably different from the other morphological features.  The number of species per sample increased from 1 to 12 species at the beach slope (S B ) to 3 to 21 at the first slope (S S1 ), increasing further to 5 to 21 species in offshore direction (S S4 ). The highest number of species were found at the seaward slope of the second breaker bar (S S2 = 3 to 22 species) (Figur 6C). Explanatory variables average bed profile (profH), sedimentation-erosion after one year (SE) and the morphological features were used as fixed effects in the ME-model for the species number (ME SN ) ( Table 3). Model comparison between the ME SN model with and without morphological features showed the significance of the morphological features (ME SN AIC = 3703, ME SN no-MF AIC = 3744, p < 0.0001).
Median diversity at the beach slope was the lowest (H' = 1.0 (IQR = 0.5-1.3)) and increased in the offshore direction up to 1.6 (IQR = 1.4-1.9) at the deepest location (S S4 ) ( Figure 6D). The environmental variables average bed profile (profH) mud content, slope and the morphological features described the diversity best and were used as fixed effects in the ME-model for the diversity (ME H ) ( Table 3). Model comparison between the ME H model with and without morphological features again showed the significance of the morphological features (ME H AIC = 958, ME H no-MF AIC = 972, p < 0.0001). The diversity index of the beach slope (S B ) and the seaward slope (S S3 ) were considerably different from the other morphological features. Table 3. Model estimates for the selected explanatory variables (fixed factors) and variance (random effects) for the models predicting, total abundance (ME AB ), total biomass (ME BM ), species number (ME SN ), and diversity (ME H ).

Species Composition Correlated with Environmental Parameters
Multivariate analysis using pRDA and the partitioning of variance revealed that the overall variance in the species abundance data explained by the selected variables is 38%, leaving 62% of the variance unexplained. The latter represents mostly sampling variance: the random effect of sampling (or not) of an individual in a relatively small sample. Such division between explained and unexplained variance is fairly typical of benthic samples. The parameters median grain size, mud percentage, organic matter, annual sedimentation-erosion, the slope of the bed, bed level, morphological features, and orientation of the slope (aspect) were selected for the model. The explained variance was partitioned into three classes (a) the morphological features, (b) the covariates Year, Transect and ProfH, used to remove the general effect of the temporal and longshore variation and the general decreasing water depth, and (c) the environmental parameters for the sediment characteristics (D50, mud content and organic matter) and the morphology (slope, aspect, bed level, and annual sedimentation-erosion).
In the multivariate analysis, the largest part of the variance explained was related to the covariates (75.5%) including their correlation with the environmental parameters and the morphological features. The explained fraction by environmental parameters (53.5%) and by morphological features (48.9%) were almost equal (Figure 7). The variation shared between the environmental parameters and the covariates is 39.9%. This is mostly caused by the fact that environmental parameters vary between years and in longshore direction. Therefore, the largest part (approximately two thirds) of the variance explained by the environmental parameters correlates with the covariates. As expected, the independent contribution of the covariates is also the largest (30.1%). The independent contribution of the morphological features (10.9%) is almost three times larger than the independent contribution of the environmental parameters (3.5%). Figure 7. Division of the explained variation in the species abundance over the environmental parameters (slope, sedimentation-erosion, bed level, aspect, D50, organic matter, and mud percentage), the covariables (alongshore variation, temporal variation, and average bed level) and the morphological features.
Permutation tests showed that the first three axes of the pRDA, with the covariables (Transect, Year, and profH) removed, were significant (p = 0.001). The pRDA given for the first two axes reveals the relationship between morphological features and benthic species (Figure 8). For visibility, only species with a contribution of more than 0.05 to the axes are shown. Additionally, species representative as indicator species and dominantly present at the specific morphological feature are shown with a color for each morphological feature.
In the plot the bar crests (B 1 , B 2 ) and the landward slope (S L ) are positioned close together with species such as the sand digger shrimp Bathyporeia elegans, a subsurface deposit-feeder preferring medium to coarse sediment, the bristle worm Spiohanes bombyx, mostly found in coarse sediment with low mud content and the free-living bristle scavenger worm Nephtys cirrosa found in muddy to coarse sediments. The troughs (Tr B , Tr 1 , and Tr 2 ) are positioned opposite to the bar tops. In the troughs, we found species favoring areas with less dynamic conditions such as Capitella capitata and Spio martinensis. The seaward slopes S S1 , S S2 , and S S3 are also opposite of the breaker bars but more to the left, indicating a stronger relation with mud and organic matter. Here we found species such as Limecola balthica, Nephtys hombergii, and Urothoe that may depend on the slightly finer fraction in the sediment. Finally, in the far left of the plot and at right angles with the bar-trough distinction (i.e., poorly correlated), we find the species related to the beach slope (S B ) and to some extent to the deepest location (S S4 ). These morphological features are related to the bristle worm Scolelepis (Scolelepis) squamata able to both suspension and deposit-feeding and preferring a large range of sediments, and the amphipods Bathyporeia pelagica, a subsurface deposit-feeders preferring medium to coarse sediment and Haustorius arenarius an intertidal species able to withstand very exposed habitats. The position of the species in the pRDA-plot corresponds with the general picture of the kite diagram ( Figure 5), showing that some species are mostly found at specific locations (very shallow, at the crest and troughs of the bars or in the deeper part) and also coincides with the largest part of the assemblage at these locations. Figure 8. Ordination diagram first two-axis pRDA analysis based on the transformed abundance data with the environmental parameters conditioned for the spatial and temporal variation and the average bed profile. Parameters with continuous values (D50, bedlevel, Slope, SE, OM, and mud) are given as arrows. Parameters described as a factor are given with triangles (black = morphological features, grey = aspect). Key species per morphological feature are indicated with a color (Green = S B , brown = S L , yellow = B 2 , light blue = S S1 , sky blue = S S3 , cyan = S S4 ).

Discussion
In our study area, Magelona johnstoni, Nephtys hombergii, Spio martinensis, Limecola balthica and Ensis were found most frequently. These are typical species of the dynamic shallow sandy coast [17,[55][56][57] and characterize the Magelona-Ensis community as defined for the very nearshore waters of the Belgian part of the North Sea [58].
The univariate analysis of the environmental parameters and the macrobenthic species in the near shore revealed a cross-shore gradient from the beach (S B , Tr B ) to the zone where the waves break (B 1 to S S2 ) with the bar crest and troughs standing out, and the slope of the outer breaker bar to deeper water (S S3 , S S4 ). In the offshore direction the number of species increased in the seaward direction ( Figure 6C). This contradicts with the general scheme given by MacLachlan and Brown [35] who indicate a decrease in species number from the low water line to the point where the waves break. The pattern of increasing species number in combination with slightly lower abundance and biomass at the bar crest is comparable with earlier findings in cross-shore profiles for Ameland and Schiermonnikoog [59] and along the Holland coast near Egmond [23,55], Terschelling [60], and Spiekeroog [61]. Nevertheless, spots of high diversity and abundance in the troughs, as described for Egmond [55], have not been recorded during our five-year survey. Possible reasons are the different physical conditions at Egmond (coarser sediment, steeper bars and troughs, and different hydrodynamics). Another likely reason is they sampled in a patch of Lanice conchilega, an ecosystem engineer capable of altering the sediment and attracting other species [62]. Despite our large sampling effort, we did not find patches with higher densities of Lanice conchilega in one morphological feature compared to the rest of the profile. Lanice conchilega was equally present over the whole cross-shore profile.
The near shore is a highly dynamic environment. In our study area, yearly net sedimentation-erosion rates were over half a meter between the breaker bars ( Figure 3D). Strong tidal currents and wave action reach the bottom and stir up fine particles of sediment and organic matter. Morphological features, such as bars and troughs, interact with the waves and tidal flows, causing wave breaking and more stirring up of the sediment [28,32]. The coarser particles end up at the beach and the finer particles are brought to slightly deeper and calmer water just outside the breaker-zone. In the nearshore of Ameland and Schiermonnikoog, the result of these processes was reflected in the median grain size, which gradually changed from 205 µm at the beach slope to 170 µm in offshore direction and the low content for mud and organic matter (Figure 4), which is common for exposed sandy shores [63,64].
Including morphological features as a factor improved the explanation of the species distribution in the nearshore of the Dutch barrier islands, Ameland and Schiermonnikoog. This indicates that species distribution in the nearshore relates not only to water depth and sedimentological parameters, often presented in studies relating benthic species to the environment [14,16,17,24,65,66]. This knowledge can be used to assess the effect of human impact on the benthic species distribution. For instance, in the early-design phases of a nourishment, as it helps to predict the ecological impact of morphological alterations due to the nourishment. Our results showed that the contribution of morphological features such as bar crests, troughs, and slopes was three times larger than the environmental parameters together, above and besides the general impact of the slope towards deeper water (Figure 7). The most likely reason is the fact that the morphological features are the expression of the mutually correlated cross-shore hydrodynamic regimes related to the bed topography in barred systems. As such, the morphological features summarize a combination of wave impact on the sediment, erosion-deposition dynamics, bottom shear stress, sediment grain size distribution, and possibly other characteristics that systematically vary between the different morphological features. Especially the hydrodynamic variables, such as the time-varying bottom shear stress, are very difficult to estimate from field sampling even-given the high spatial resolution needed-from dedicated models. The morphological features, in contrast, are easy to observe. It is shown in our study that they contain information on relevant environmental variables beyond the easily measured ones, but they function as black boxes: the exact causal mechanisms behind the correlations remain unclear. However, closer examination of the species typical for the different morphological features reveals some of the important contributing factors.
Beach slope-The beach slope is the transition between the intertidal beach and the subtidal beach with almost no water depth, waves rolling up and down the slope resulting in coarser sediment compared to the deeper part of the profile. We found a very distinct assemblage of species at the beach slope with Scolelepis (Scolelepis) squamata -Haustorius arenarius -Bathyporeia pelagica. These species represent the intertidal beach [23,55,[67][68][69] and the beach slope just below the low-waterline [55]. H. arenarius [70] and B. pelagica [71] prefer clean oxygen-rich sand, dig very fast into the sand, and search for food attached to the sand grains. S. squamata can be found in the North Sea at the slope of exposed beaches and in very dynamic sands offshore [61]. The species burrows in sediments with low mud content [72] and a high permeability allowing their ventilation current for respiration [70]. Mud-free, permeable and well-oxygenated sediment is a common requirement for these typical species. It is a result of frequent stirring and movement of the sediment top layers.
Bar crest-Waves approaching the shore increase in height, change shape, and eventually break at the second or first breaker bar creating turbulence and stirring up of the sediment. These circumstances relate to lower abundance, biomass, and species number compared to other parts of the profile ( Figure 6). Typical species we found here were Pontocrates and Bathyporeia elegans at the first breaker bar (B 1 ) and Spiophanes bombyx, Nephtys cirrosa, and Bathyporeia guiliamsioniana at the second breaker bar (B 2 ). B. elegans is a sand hopper preferring small cavities in the sand where it waits for sand grains to deposit [72]. The genus Pontocrates prefers medium sand and has a burrowing life mode. Both taxa are highly mobile and adapted to frequent movements of the sediment. The small slender bristle worm S. bombyx lives in a sandy tube and can withstand moderately strong tidal currents and some wave action. The polychaete N. cirrosa is a eurytope species with preference for rather dynamic conditions, it burrows through the sand and predates on small invertebrates and prefers clean sand. It is found at high densities in the Southern Bight, mostly at the Brown Bank and near the coast [56,72]. B. guiliamsoniana is the largest species of the genus Bathyporeia and resembles the other species within the genus. Both B. elegans and B. guiliamsioniana are predominantly subtidal species extending their distribution to the deeper part of the shore [71].
Bar trough-The organic matter content in the nearshore is low. However, in the troughs, median grain size is smaller compared to the bar crests and a little extra organic matter is present compared to the other morphological features (Figure 4). These conditions were reflected by finding Nephtys longosetosa and Capitella capitata in the trough areas. N. longosetosa burrow through the sediment, feed on small invertebrates, and prefer fine to medium sand. It occurs in the shallow waters along the coast, the Wadden Sea and the Dogger Bank [72]. C. capitata is an opportunistic species capable of building tubes at or just below the seabed. It is a deposit feeder ingesting microorganisms and detritus [72].
Landward slope-The landward slope is closely related to the bar crest (Figures 4 and 8) and this is also reflected in the key species found for this feature. The key species identified for the landward dipping slope were the bristle worms Nephtys cirrosa and Scolelepis bonnieri, and the worm Magelona mirabilis. M. mirabilis is a deposit feeder and mostly found in sandy sediments in the intertidal zone [73]. Like Scolelepis (scolelepis) squamata, S. bonnieri lives in mobile sand and builds loosely constructed burrows [72].
Seaward slope-The hydrodynamic conditions at the slope of the outer breaker bar changes rapidly between the upper part (S S2 ) and the lower part (S S3 and S S4 ). Wave action decreases from the upper to the lower part. In the latter, wave action may be absent during calm conditions [74]. Key species identified at the foot of the outer breaker bar (S S2 ) were Echinocardium cordatum, Spisula subtruncata, and Abra alba followed by Limecola balthica and Urothoe. All species are larger compared to the species of the breaker bars. The bivalves S. subtrucata, L. balthica, and A. alba prefer medium sand to muddy sand. They are suspension feeders but also capable of deposit-feeding [72]. Except for A. alba, these species have a lifespan of more than 2 years. Further offshore (S S4 ) the slope flattens, and the waves reach the bed only during moderate weather conditions. The influence of the long-shore current is stronger resulting in sediment with less mud. This is characterized by a large group of species consisting of bristle worms Malmgrenia, Magelona johnstoni, Scoloplos armiger, Eumida and tube building worms Lanice conchilega in combination with the bivalve Fabulina fabula. Figure 9 summarizes the key species and their traits along the continuum from beach to deeper parts. With decreasing wave and current energy exerted on the sediment, we observe a transition from highly mobile to sedentary species, from particle-centered deposit feeding, over bulk deposit feeding to suspension feeding, and from dependence on clear high-oxygen sands to sands with an admixture of mud. Also, the size and longevity of the animals increases along the general gradient. We note, however, that the bar crests and the portions of the slopes immediately surrounding them are characterized as focal points for wave energy dissipation. In accordance, the fauna in these morphological features is adapted to a higher level of hydrodynamic energy dissipation than elsewhere. Figure 9. Conceptual drawing of the cross-shore wave related processes (adapted from van der Zanden (2016) [27]), environmental parameters measured and the macrobenthic species identified as indicator species for the morphological features.

Conclusions
The distribution of benthic species in the near shore with a barred system at Ameland and Schiermonnikoog was studied taking morphological features (slopes, bar crests, and troughs) into account. Information was used of a detailed 5-year monitoring campaign measuring benthic species and sediment composition within each morphological feature. Due to this extensive dataset we were able to account for the alongshore and temporal variation within the benthic species distribution found at each location. Although depth is an important factor explaining species distribution, our study showed that the morphological features in themselves contribute three times more to the variance explained of the species variability than the environmental parameters describing the sediment characteristics and the morphology together.
The benthic species found were typical for the nearshore communities in the North Sea. Specific indicator species were found for the slopes, bar crests, and troughs, which form micro-habitats. These differed in their ecological traits, showing adaptations to the levels of stress in their preferred habitats.
Because of their correlated hydrodynamic, sediment texture, and depth characteristics, morphological features can be used as an effective proxy that incorporates most of the factors explaining the distribution of benthic species. In addition, morphological features can easily be extracted from bathymetry maps.
Incorporating explicit spatial information of morphological features can contribute to the prediction of benthic species distribution and community composition in the nearshore. This may facilitate the assessment of human impact (e.g., nourishment) or natural processes on the benthic species distribution, as this impact is often translated into changes of the morphological features.
The results presented are based on data from two locations that were similar in terms of morphological features. To assess the importance of the morphological features as a structuring factor for the benthic species community in general, additional studies at locations with explicit morphological features is required.