Mollusc shell detritus affects benthic subtidal community dynamics in the Northern Wadden Sea

A shallow subtidal area in the northern Wadden Sea was monitored over 17 years (from 2003 to 2019) for sediment parameters and macrobenthic fauna. Due to the sheltered position of the study area, sediment composition remained rather stable with only minor annual and seasonal variations in sediment granulometry. An intermittend storm (‘Kyrill’) had no significant effect on sediment composition parameters; the construction of an artificial dune along the southern border of the study area had minor transient but no lasting effects on sediment composition. Faunal species richness and total abundance showed the typical seasonality with minimum abundance in late winter and a peak after recruitment in summer. Variations between years in autumn (post-recruitment) abundances were best explained by the number of days with a seawater temperature <1°C during the preceding winter. Temperature during other seasons, salinity, chlorophyll concentration and NAO showed no significant correlations with total abundance, nor did storm (‘Kyrill’) or construction of the artificial dune. Within-site faunal variability was best explained by water depth and velocity of the tidal currents while sediment granulometry was of minor importance. However, the amount of bivalve shell detritus mixed with the sandy sediment proved to be the strongest covariant of species numbers, total abundance, and species-specific abundances. At the sediment surface, shells provide the only anchorage for epibenthic species which in turn attract associated fauna. Shells inside the sediment hamper movement of infaunal predators and epibenthic predators are handicapped in rooting about for prey. Thus, shell material is a highly significant structural factor for the macrozoobenthos in these shallow waters. In a future with increasing ocean acidification, the availability of benthos as a food source for higher trophic levels will depend on the balance between pH-driven shell destruction and the compensatory power of shell-building species.


Introduction
The northern Wadden Sea is a unique landscape renowned as a bird resting place during their annual migrations between breeding and hibernation areas. Its importance as a resting place mainly originates from the food supply offered by intertidal and shallow water benthos; hence, benthic organisms play a crucial role in bird population development. Consequently, the longterm development of benthic species in the Wadden Sea has been well studied for the tidal flats, in particular in the Dutch Wadden Sea (Beukema and Decker 2020a). The temporal variability of macrozoobenthos was mainly driven by winter temperatures, with gradual climate warming increasing total abundance and biomass over the past 50 years Beukema and Decker 2020a;Beukema and Dekker 2020b). Several factors contributed to this net effect of global warming, including higher over-winter survival of cold-sensitive species (Beukema 1989(Beukema , 1992, enhanced growth during warmer summers (Beukema et al. 2017a), and rather complex predator-prey interactions (e.g. Beukema et al. 2000;Strasser and Günther 2001). Few species such as the clam Macoma balthica suffered losses from warming (Beukema et al. 2017b), but their losses were over-compensated by introduced species such as razor clams Ensis leei and pacific oysters Communicated by U. Schückel This article is dedicated to the memory of Christian Hass ( †13.9.2020).
This article is a contribution to the Topical Collection Biodiversity and Ecology of the Wadden Sea under changing environments Crassostrea gigas. As a result, total biomass increased over time but since some of the introduced species grow outside the size ranges suitable for wading birds, the food supply available for these birds showed no long-term trend over the past decades (Beukema and Dekker 2019). Among other factors, eutrophication affected the macrozoobenthos only for a short period of time until legislative regulations reduced anthropogenic nutrient release again (van Beusekom 2005;Philippart et al. 2007) while sea level rise (in the Wadden Sea some 2 mm year− 1 ) was largely compensated by sedimentation, so far (Beukema 2002;. Contrary to this rather comprehensive knowledge on the long-term performance of intertidal benthos, far less is known about temporal variability of subtidal benthos although the areal extent of the subtidal by far exceeds the area of tidally emerging flats. So far, only the macro-scale changes that occurred in the subtidal during the past century have been analysed in more detail. These include the disappearance of the reefs of colonial polychaetes Sabellaria spinulosa due to dredging and trawling activities and the disappearance of natural oyster (Ostrea edulis) beds due to overexploitation (Riesen and Reise 1982;Reise and Schubert 1987). Mussel (Mytilus edulis) beds extended their range, partly in replacement of extinct oysters and partly due to anthropogenic mussel culture plots (Reise et al. 1989). Finally, subtidal seagrass (Zostera marina) beds decayed in 1933/34 (Wohlenberg 1935) due to wasting disease brought about by the slime mould Labyrinthula zosterae (Muehlstein et al. 1991) and never recovered (Reise et al. 1989).
Sabellaria, Ostrea, Mytilus and Zostera are all epibenthic ecosystem engineers, each of them creating a typical habitat type with many associated species. Therefore, changes in their presence are likely to change the infaunal communities. In addition, since epibenthic structures affect small-scale hydrodynamics there may be changes in sediment characters as well, and potential feedback of sediment composition on faunal composition. Unfortunately, for the subtidal Wadden Sea, there were no comprehensive long-term data, neither on infaunal performance nor on sediment composition. Therefore, we initiated a long-term study on sediments and fauna of the shallow subtidal off Königshafen near the island of Sylt to answer the following questions: (1) Are there long-term trends for change in faunal composition or sediment granulometry over years or just accidental variability between years? (2) What are the agents that rule (or at least correlate with) faunal variability? Can climate effects (such as cold and warm winters, e.g. Beukema and Dekker 2020b) known from intertidal flats be transferred to the adjoining shallow subtidal? Or is subtidal variability in the Wadden Sea more similar to the North Sea (e.g. Meyer and Kröncke 2019)? (3) Storms strongly structure species distributions in the adjoining coastal North Sea (Armonies et al. 2014)-may have storms the same effects in the Wadden Sea subtidal?

Material and methods
The study area (Sylt-Rømø bight) is located in the northern Wadden Sea, which is the fringe of the NE North Sea that is sheltered by barrier islands (Fig. 1 top). Inside this bight, we investigated the shallow subtidal area of some 4 km 2 between the tidal inlet 'Lister Ley' (the southward prolongation of the main tidal inlet 'List deep') and the intertidal flats of Königshafen from 2003 to 2019 (Fig. 1 bottom).
Using a rectangular grid, 50 sampling positions were fixed between mean low tide level (mltl) and the deep channel of Lister Ley (roughly along the 5 m depth contour, mean tidal range +1.8 m). During sample collection, the ship always piloted the fixed sampling positions (GPS navigation) but was free to drift during subsequent countersink of the corer; thus, the sampling design is 'stratified random'. For comparison, the southern part of the Königshafen intertidal (yellow dots in Fig. 1) was sampled in parallel during the first 10 years from 2003 to 2013.
For subtidal sampling, we used Research Catamaran 'Mya' until it was replaced by RV 'Mya II' in 2013. The latter ship has a deeper drought which enforced a reduction in the number of sampling sites to 33, i.e. the shallowest sites (orange dots in Fig. 1) were only studied during the first 10 years. Sampling was done with a 0.02-m 2 box corer of the Reineck-type, always around preceded high tide. A subsample of 10-cm 2 surface area was taken from each core for granulometric sediment analysis and the remaining sample was sieved through 1-mm meshes. The retained material was fixed in buffered formalin solution and later sorted and analysed for fauna. Granulometric analyses were done with a CILAS 1180L diffraction laser particle size analyser which provides grain size information between 0.04 μm and 2 mm. Particles (shell fragments and occasional stones) that had been retained by the sieve were quantified during faunal analyses as total wet weight (>1 mm) and weight of very coarse particles (>4 mm) retained on an additional sieve. For intertidal sampling, we used a hand-operated corer, with faunal and sediment analyses as above.
Sediment granulometric composition is characterised by the median grain diameter of the sediment fraction <2 mm (D 50 ) and by the 10 th (D 10 ) and 90 th (D 90 ) percentiles of the accumulative grain size curve; i.e. 10% (by weight) of the grains in a sediment core are finer than the value given by D 10 and another 10% are coarser than the value given by D 90 (details in Blott and Pye 2001). Since the values of these three sediment parameters all derive from the same graph, they are correlated and partly redundant. We nevertheless used all three of them, hoping to get an idea on which sediment component may vary the most over time, and which one might exert most influence to each of the species. A general description of the sediments of the study area was already published by Dolch and Hass (2008). Altogether, the subtidal grid was sampled 76 times but the sampling frequency varied. During the first 4 years, we aimed at monthly sampling (limited by ship availability and weather conditions) to record short-term changes, then (2007 to 2013) at seasonal sampling, and after 2013, we only sampled annually in September/October to record annual variability (Table 1; sampling dates and position details in Online Resources 1 and 2). Accordingly, statistical analyses are based on different subsets of the data collection: analyses for seasonality and comparisons between the subtidal and the adjoining intertidal are based on the 2003 to 2013 period with seasonal sampling of the initial 50 subtidal (and 12 intertidal as appropriate) sampling sites while analyses for long-term change are based on autumn data of the 33-site subset of deeper sampling positions that was sampled over the entire period (Table 1) An adaptation also occurred in sediment analyses; until 2006, the sediment subsamples comprised the top 4-5 cm of the sediment which is habitat for most of the fauna. But when we noted that the sediments were rather stable in the study area, we switched to analysing the top 1-cm sediment layer only from 2007 onwards, hoping to get a finer resolution for potential changes in surface sediment composition. This shift in sampling depth significantly affected the sediment parameters (D 10 , D 50 , D 90 ) used in this study. To eliminate these method-related differences, the series of 4-5-cm and 1-cm sub-sampling depths were separately normalised (xmean of the respective sampling position) to be comparable (Online Resource 1).
When storm 'Kyrill' passed the study area on 18 and 19 Jan 2007, we used this as an opportunity to study storm effects on the subtidal benthos by before-after comparisons of sediment granulometry and fauna, respectively. Another unique disturbance was the construction of an artificial sand dune for coastal protection along the southern border of the study area. This was done by dredging sediment in an offshore area, ship transport into the Wadden Sea, and flushing the suspended sediment onto the beach. Coarser sand grains remained on the beach but small grained particles and shells were partly washed ashore with the run-off. This sand replenishment was conducted in two campaigns in summer 2013 and summer 2014. Again, we used before-after comparisons to check for significant effects of this coastal protection measure on ambient sediments and fauna.

Sediment
Most of the subtidal study area showed a median grain size in the fine sand range, with a small spot of very fine sand and highest mud accumulation between mean low tide level and 1 m below, and coarser sand in the northwest close to the main tidal channel of Königshafen bight (Fig. 2a, b). The adjoining intertidal had similar median diameters and fine particle content but sharply deviated from the subtidal in the amount of coarse and very coarse sand fractions. With respect to coarse sediment compounds, neap low tide level (mean low tide level +10 cm in this area) formed a very distinct separation line (Fig. 2c).
Granulometric sediment composition showed significant variations both annually and seasonally (Table 2). Annually, years with coarser and finer sediment are not arranged in random order but sequences of years with decreasing particle sizes were followed by a number of years with increasing grain sizes (Fig. 3). Only the year 2016 made an exception with apparently out-of-order fine sediments.
Seasonally, the sediment was slightly coarser in winter (Fig. 3, right) but different seasonal patterns are seen on shorter time scales (Online Resource 1). Including water depth, latitude, or longitude of the sampling positions as covariants in ANCOVAs of sediment parameters vs. seasonality all resulted in non-significant (p>0.05) covariant effects (Online Resource 1). The amount of shells mixed in the sediment only varied over years but there was no significant seasonality (Table 2).
During the beginning of this study, most of the study area was covered by a superficial mud layer that vanished within half a year and never occurred afterwards. This transient event affects analysis of the long-term tendency of change in sediment granulometry: including the year 2003 results in a significant increase in median grain size (D 50 ), while omitting this year results in a significant decrease of D 10 (Table 3). However, in both cases, r 2 -values <0.03 indicate that granulometric sediment composition remained basically stable over the 16-year study period.
When storm 'Kyrill' passed the study area on 18 and 19 Jan 2007, it increased high tide water level in Sylt-Rømø bight by up to 2 m. However, before-after comparisons indicate that this did not significantly change any of the sediment composition parameters in the study area (Online Resource 1).
During construction of the artificial dune along the southern border of the study area, suspended sediments from the run-off of the flushing area significantly increased average diameter of  Error D 10 , D 50 , and D 90 were all position-wise normalised prior to analyses. SS sum of squares, DF degrees of freedom, MS mean square, F F-statistic, p probability the finest 10% of the sediment (D 10 ) in 2013 and 2014 but there were no lasting effects on granulometric sediment composition after the construction activities in the subtidal section. However, bivalve shells and shell fragments that had been included in the nourishment sediment were partly washed ashore with the runoff and remained in the subtidal study area. This significantly increased the amount of shell detritus in the subtidal and the elevated level of shell detritus content persisted until the end of this study (Online Resource 1). Once constructed, the artificial dune partially eroded again. While this had no significant effect in the subtidal, sand suspended from the dune during storms changed sediment composition in the intertidal and created a sandy hook in the intertidal over time (Fig. 4). In 2020 and 2021, this new hook approached the older 'List hook', separated only by a small tidal inlet. Since then, the area enclosed by both hooks turned increasingly muddy.

Differences between tidal levels
Some 150 benthos taxa were recorded during the 10-year period from 2003 to 2013 of joint intertidal and subtidal sampling (Online Resource 2). ANOVAs on the effect of tidal level (i.e. categorial factor intertidal or subtidal) revealed that 21 of the species were significantly more abundant in the intertidal and 33 the subtidal while two-thirds of the species were indifferent to tidal level ( Fig. 5a; Online Resource 3). In the intertidal, these three groups contributed rather equally to the species pool while abundance was strongly dominated by intertidal species (Fig. 5b). In the subtidal, the species pool was dominated by indifferent species while both intertidal and subtidal species contributed equally to abundance (Fig. 5c). Thus, in terms of abundance, intertidal species play a major role in the shallow subtidal while subtidal species are rather negligible in the intertidal.
Mean species number per sediment core constantly increased with water depth until limited along the physically unstable rim of the adjoining tidal channel (Fig. 6a). Abundance showed the same pattern in the subtidal but then strongly increased in the intertidal (Fig. 6b).

Seasonality
Both abundance and species density showed the typical seasonality with highest values in summer and lowest in late winter though development was not entirely synchronous: the spring increase of abundance started earlier (in April) in the intertidal against June in the subtidal (Fig. 7 bottom left; ANOVA results on the species level in Online Resource 3). As a result of common seasonality, there are positive correlations between inter-and subtidal abundances (r=0.5104, r 2 =0.2605, p=0.0003) and species numbers per core (r=0.6038, r 2 =0.3646, p<0.0001). But within seasons, significant correlations only occurred between inter-and subtidal winter abundances (r=0.7350, r 2 =0.5402, p=0.0378) and species numbers (r=0.7845, r 2 =0.6154, p=0.0212) while there were no significant correlations during the other seasons. The number of species recorded per sample set was remarkably constant throughout the year (one-way test of significance, p>0.9 for both tidal levels).

Long-term trends in the subtidal
In the subtidal, species density and total abundance significantly varied between years (Table 4, Fig. 8). Since the number of species was highly correlated with log(abundance) (r=0.8446, r 2 =0.7134, p<0.0001), both parameters changed largely in parallel over years. Both species density and abundance showed statistically significant long-term temporal trends towards higher numbers; however, very low r 2 -values indicate that the increases were quantitatively not important (Table 5). Species richness per sample set (33 cores each) did not significantly change over the 17-year period (p=0.1496).
Most of the higher taxa showed corresponding variability over years; only oligochaete abundance was remarkably stable (Table 6; in echinoderms and tunicates, the lack of significant effects is presumably due to low overall abundances). Polychaete, bivalve and amphipod abundances increased over time while gastropods, decapods and Pycnogonida decreased; however, again the very low r 2 -values indicate that these are rather spurious results.
On the level of single species, abundance of nearly all of the more abundant species correlated positively with total abundance of the other species; hence, abundance of most species followed the same temporal pattern (Online Resource 4). Among the potentially controlling agents, winter cold correlated the best with total abundance in the following year; the number of days with a seawater surface temperature (SST) <1°C gave the best fit (Table 7). Temperature during other seasons, chlorophyll concentration and NAO showed no significant correlations with total autumn abundance. Salinity correlated significantly positively with autumn abundance only in May (Table 7). Contrary to that, there are numerous significant correlations between these factors and abundance of single species (Online Resource 4), some referring to wellestablished physiological constraints or ecological chains of effects, others still unexplained or (because of the multitude of tests) just spurious.
Storm Kyrill (18 and 19 Jan 2007) increased mean low and mean high water levels in Sylt-Rømø bight by up to 2 m. This neither changed sediment composition (see above) nor reduced benthos abundance or species numbers. Instead, in autumn 2007, abundance and species numbers were even higher than those in the previous and following years (Table 8).

Linking faunal variability with sediment data
Abundance and species richness significantly varied over classes of sediment parameters (Table 9). On an average, both numbers increased with increasing grain size values (Fig. 9). However, ANCOVAs (Online Resource 5) revealed that the amount of shell gravel >4 mm by far had the strongest effect, with strongly increasing total abundance and species numbers with increasing shell detritus. With very few exceptions, this also applied on the species level (Online Resource 5). Closest correlations with the amount of shell gravel occurred in the more epibenthic species that need hard substrate for attachment (e.g. hydrozoans, anthozoans, blue mussels Mytilus edulis and slipper limpets Crepidula fornicata) and in species associated with them. But positive correlations also occurred in entirely infaunal species such as oligochaetes and the polychaete Scoloplos armiger while few species were negatively affected (Online Resource 5). In a few species, the effect of shell detritus changed seasonally. The abundance of Nephtys species, for example, correlated negatively with the amount of shell detritus in autumn but not on an annual base.

Sediment stability
The surface sediments proved to be highly stable over the studied 17-year period; obviously, for the study area, Sylt island is a highly effective barrier to North Sea waves. Accordingly, no significant change in sediment granulometry could be detected after storm Kyrill; more generally, we could not find any correlation between sediment changes and storm frequency or intensity (though this may partly be due to insufficient data), or NAO indices. Thus, higher hydrodynamics brought about by global change is unlikely to affect sediment composition in the study area unless it becomes high enough to decrease the barrier function of the island itself. This stability is in strong contrast to the tidal channels which have the highest morphological activity along the East Frisian coast (Winter 2011). However, morphological changes in these tidal inlets are likely to influence the Fig. 6 Macrozoobenthos species density (a mean number per core) and abundance (b logtransformed) along the tidal gradient. Vertical bars indicate the 95% confidence ranges; water depth refers to mean tidal level routes of the tidal waters and the position of small-scale turbulence areas in the following tidal channels. Thus, we suggest that the periods of change towards finer or coarser sediment we observed in the study area (Fig. 3) might be brought about by Fig. 7 Seasonality of macrozoobenthic species number per sample set, species density per sediment core, and abundance per 0.02m 2 (log-transferred) in the studied inter-and subtidal sites. Vertical bars indicate 95% confidence ranges  morphological changes in the tidal inlet 'List deep' and/or its southward branch 'Lister Ley'. This implies that the tendency of change towards finer or coarser sediment will depend on the position of the studied area within the catchment area of the tidal inlet: local change cannot spatially be upscaled. However, the observed long-term stability of sediments was only true for the subtidal section of the study area while the intertidal was subject to ongoing morphological changes accompanied by changes in sediment composition; sediment delivery from erosion of the artificial sand dune along the southern border of the study area still increased these changes (Armonies 2020).

Faunal variability
The overall pattern of macrofaunal distribution with highest species richness in the subtidal and highest abundance in the intertidal of the study area conforms well with previous findings in the southern North Sea . Within the subtidal, the spatial variation of macrofaunal abundance and species richness correlated positively with water depth and negatively with increasing current towards the tidal channel (Fig. 6); the same factors were found to best explain the spatial variations in macrofaunal community structure in the tidal channels of Jade Bay (Schückel et al. 2015). This highlights the importance of water depth (respectively, tidal level) and current velocities for macrofaunal community structure in the Wadden Sea. But despite these common structuring factors, there are major differences between the macrobenthos communities of the northern and southern Wadden Sea, likely attributable to differences in sedimentological characteristics (Schückel et al. 2015). Seasonally, abundances in the inter-and subtidal were expected to correlate to some degree because individuals of both tidal sections derive from a common larval pool and are widely subject to the same physical factors. Within single species, we only found minor differences in seasonality, presumably due to slower warming of the subtidal in spring. Nevertheless, total macrozoobenthic abundances in the inter-and subtidal were uncorrelated; hence, abundance estimates from one tidal level cannot be generalised to the other. Presumably, this is due to species-specific preferences for tidal levels and levelspecific predation pressure which combined result in different community compositions.
In the subtidal, abundances and species numbers significantly varied over years (Fig. 8). Just as in sediment parameters, long-term change occurred in 'waves', i.e. some years with increasing numbers were followed by a sequence of years with decreasing numbers. However, the timing and duration of faunal waves do not coincide with the waves of changing sediment parameters. Therefore, we assume different causes.
In the intertidal of the Dutch Wadden Sea, species numbers have been shown to increase over the past 50 years (Beukema and Decker 2020a); this was partly due to global warming favouring winter-sensitive species, and partly to the arrival of new species while no established species disappeared. In the shallow subtidal of the Sylt-Rømø Bight, we see the same trend (Fig. 8), and as in the Dutch Wadden Sea, this is partly due to the arrival of new species (Lackschewitz et al. 2015). However, with respect to winter-sensitive species, our data  either show no significant trend (e.g. Lanice conchilega and Nephtys spp.) or even a decreasing trend as in Carcinus maenas and Crangon crangon. For the first two species, the lack of a long-term temporal trend may be due to the virtual absence of cold winters with massive ice formation during the study period (but see below). But the decreasing trend in the latter two species cannot be explained by winter temperatures; possibly they suffered from competition for food with introduced species as Hemigrapsus takanoi (Cornelius et al. 2021) or from changing predation/fishery pressure (Temming and Hufnagl 2015). Anyway, since the latter two species are important predators on infauna (in particular on bivalve recruits, Beukema and Decker 2020b; Beukema and Dekker 2020a), their decrease in abundance may have contributed to the overall increasing trends of abundance and species number. Abundance of total macrozoobenthos in autumn correlated negatively with SST in the preceding winter and significantly positively with the number of cold days, with best fit for the number of days with an SST <1°C (Table 7). However, 'really' cold winters with massive ice cover of the Wadden Sea were frequent in the preceding decades while during the study period only the winter of 2010 was relatively cold; in the same year, total abundance and species density peaked (Fig. 8). Apart from this year, we assume that winter temperatures during the study period were high enough to avoid massive mortality of winter-sensitive species but still low enough to promote recruitment, either by lower survival of predators like C. maenas and C. crangon or by reduced recruitment success of the latter species by temporal mis-match effects (Strasser and Günther 2001). Thus, benthos abundance in autumn seemed to be fine-tuned by subtle differences in winter sensitivity between predators and their prey, though other factors like variations in fishery intensity and top-predator abundances (such as seal Phoca vitulina and harbour porpoise Phocoena phocoena) are likely to contribute to local changes of epibenthic predation pressure. This fine-tuning may be another reason for the lack of correlations between inter-and subtidal abundances: since birds mainly feed in the intertidal and fish in the subtidal, different predator spectra may give different fine-tuning results over tidal levels.
In summary, winter temperature effects in the northern Wadden Sea during the past two decades resemble the effects monitored in the Dutch Wadden Sea during the past five decades as summarised by Beukema and Decker (2020a). In the light of global warming, this is not astonishing: meanwhile, the rise in SST (since 1962 +1.76°C, i.e. +0.31°C per decade, Wiltshire et al. submitted for publication) may have levelled the latitudinal differences between the Dutch and the northern Wadden Sea. Thus, past temperature effects in the Dutch Wadden Sea may be well suited as a model for today's effects in the northern Wadden Sea. Accordingly, current changes in the Dutch Wadden Sea such as the ongoing replacement of the clam Macoma balthica by Abra tenuis on the high tidal flats (Dekker and Beukema 2021) are expected to forecast the future in the Northern Wadden Sea.
Variations of temperature during other seasons as well as variations in NAO, salinity and chlorophyll concentrations were all checked for correlations with abundance over month-wise increasing lag-times. There were a few significant correlations on the species level but not on the community level. Compared to the open North Sea, short-term variations in these factors are common in the Wadden Sea. Accordingly, species permanently living in the Wadden Sea need to be adapted to these variations. At the same time, these variations may be a major reason for the low abundances in the Wadden Sea of a number of species that are common in the adjoining North Sea. About a third of all species recorded belong to the latter group.   . 9 Variations of species richness (left) and total abundance (right) with sediment parameters. D 10 , D 50 and D 90 median diameters of the sediment particles at the 10 th , 50 th and 90 th percentiles by weight of the cumulative particle size graph; shell detritus is the total weight per sediment core of shells retained on a sieve with 4-mm mesh size and occasionally included some stones NAO was often used as a proxy for climate variability integrating several climatic variables such as SST, prevailing wind direction and wind speed. In the shallows of the southern North Sea, NAO was clearly associated with the development of benthos until the year 2001 but no more afterwards (Kröncke et al. 1998. For our 2003 to 2019 study period, we also failed to find a significant correlation. We assume NAO mainly indicated effects of the series of cold and warm winters that occurred before the year 2000 while there were no series of cold winters afterwards. In the southern North Sea, Kröncke et al. (2013) obeserved the importance of winter cold fading away around the turn of the century. This allowed more wintersensitive species to come in and finally resulted in a regime shift (Meyer and Kröncke 2019). Acoordingly, NAO partly lost its indicative value for benthos development (Dippner et al. 2014). In addition, the relationship between climate proxies and benthic communities is not consistent across regions since they may be influenced by the local environmental conditions and the local species composition (Birchenough et al. 2015). In our study area, it is the sheltered position behind a barrier Island which prevents westerly storms associated with strongly negative NAOs to exert influence on the benthic community while the same storms may strongly structure the shallow water benthos communities in the adjoining North Sea (Nehls and Thiel 1993;Armonies et al. 2014). On the other hand, in the twentieth century cold winters, effects on benthic communities were more pronounced in the nearshore than in the offshore regions of the North Sea (Reiss et al. 2006). With respect to global warming, we thus expect that benthic community dynamics in the near-and offshore areas may become more similar to each other.

Animal-sediment relationships
In the shallow subtidal, abundance of most species was positively correlated with increasing grain size; accordingly, total abundance showed the same tendency. However, since granulometric sediment composition had been highly stable over the study period and because the entire study area was rather homogenous fine sand, most of these correlations are rather marginal, i.e. they turned statistically significant because of the high sampling intensity but explained little of the variability as estimated by r 2 numbers in regressions, or F-values in ANOVAs. The amount of coarser shell material, in contrast, turned out to be both highly significant and quantitative important, for the majority of species as well as for total abundance. In an environment of unconsolidated soft sediments, shells provide the only anchorage for the more sessile fauna, and most species apparently do not differentiate between life bivalves and empty shells. And once established, epibenthic species attract further species. Thus, positive correlations between the amount of shell material and epibenthic species had been expected. However, most of the infaunal species correlated positively with shell detritus as well. We assume that shells buried in the sediment hamper predators; infaunal predators may be restricted in movement ability and epifaunal predators in the efficiency to ransack the sediment, and both types of predators may have problems to localise and catch their prey. These results corroborate the finding that shell debris constitutes small-scale habitat structure which promotes beta diversity (Hewitt et al. 2005;Casado-Coy et al. 2022).
However, on the level of species, the importance of shell debris may turn out to be more complicated because shell effects may vary over size classes. For the Nephtys species (free-moving infaunal predators, Schubert and Reise 1986), abundance correlated positively with the amount of shells in spring and summer when recruits were abundant; but as individuals had grown to size in autumn and winter, the correlation became negative. Presumably, the shelter-from-predation function of the shells surpassed the effects of restricted movement in juveniles but not in grown-ups.
The relationship between the amount of shells and faunal abundance became most evident when we ordinated abundances according to classes of shell weight while it was less visible in simple correlations. This is so because our sampling depth in the sediment (usually) surpassed the sediment layer that was actually populated while the amount of shells was evaluated from the entire sediment core. Thus, 'amount of shells' included material from deep sediment layers that had not been populated during sampling and, therefore, could not (directly) influence faunal abundance. For further studies, we therefore suggest to consider the vertical distribution of shells and fauna, e.g. by subdivision of the sediment cores into depth layers. Alternatively, high-resolution hydroacoustic seafloor mapping may be used to classify sediments according to the amount of shell material at the sediment surface (Mielck et al. 2014) which potentially gives a chance to estimate infaunal abundance from hydroacoustic data.
Our results corroborate previous findings that mollusc shells are a structural element in the shallow subtidal that may strongly influence faunal abundance and diversity (Gutiérrez et al. 2003). Shells are continuously produced by life organisms and decomposed by mechanical stress including shell-boring organisms and chemically. Since most of the shells are based on CaCO 3 , seawater chemistry in general and pH values in particular are crucial for the rates of chemical destruction. Thus, ocean acidification by still rising CO 2 levels is likely to accelerate shell decomposition; in Sylt-Rømø bight, mean seawater pH decreased by 0.23 units during the past 40 years (Rick et al. submitted for publication). At the same time, pHdependent decreases in calcification rate, shell dissolution and reduced growth have been observed in life bivalves and gastropods Guinotte and Fabry 2008). Thus, reduced shell growth and calcification in life animals meets with increased chemical decomposition. Over time, this might result in sediments depauperated in shell fragments and, accordingly, populated by less abundant and less diverse fauna, and this effect of reduced structural complexity is likely to enhance the direct effects of ocean acidification on the species' physiology (Kroeker et al. 2011). From this point of view, ocean acidification might turn shell debris a threatened resource. However, the effects of ocean acidification vary over habitats (Waldbusser and Salisbury 2014). In the Wadden Sea, bivalve species introduced during recent decades such as American razor clams Ensis leei and pacific oysters Crassostrea gigas meanwhile attained high abundance. Since both species grow rather large, they might compensate or even over-compensate potential losses of shell material due to ocean acidification. Thus, the future availability of benthos as a food resource for higher trophic levels depends on the balance between ocean acidification effects and the compensatory power of the new species. In cases of conflicting evidence like this, further attention to the problem is advisable.

Conclusions
In the twentieth century, the Wadden Sea benthos community was strongly shaped by cold winters with prolonged ice cover. With global warming, the importance of winter cold faded away. Today, benthos community structure seems to be mainly regulated by predation, with just fine-tuning of predator effects due to differential sensitivity to winter cold. Accordingly, other factors limiting or promoting predation pressure seem to have become more important. Our results indicate that in the fine sandy sediments of the Wadden Sea, the amount of coarse shell detritus is such a factor. This deserves further attention, in particular because the shell pool in the sediments is potentially threatened by ocean acidification.
Acknowledgements This study had been impossible without the tireless operational readiness of the crews of research vessels 'Mya' and 'Mya II' during some 80 cruises. Elisabeth Herre assisted in field sampling and participated in intertidal faunal analyses and Monika Hellwig-Armonies participated in subtidal faunal analyses. The reviewers critical comments helped to improve the manuscript. Thanks are due to them all.
Funding Information Open Access funding enabled and organized by Projekt DEAL.

Declarations
Conflict of interest The authors declare no competing interests.
Ethical approval No animal testing was performed during this study.
Sampling and field studies All necessary permits for sampling and observational field studies have been obtained by the authors from the competent authorities and are mentioned in the acknowledgements, if applicable. The study is compliant with CBD and Nagoya protocols.
Data availability All data generated or analysed during this study are included in the supplementary information files.
Author contribution This study was designed by WA and Christian Hass ( †13.9.2020). Field sampling was done by CB, FM and WA who also analysed the fauna. FM performed the sediment analyses. JR provided the oceanographic and plankton data. Data collection, statistical analyses and the first draft of the manuscript were performed by WA. All authors commented on previous versions of the manuscript and all authors read and approved the final manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/ .