Small-Scale and Long-Term Variability in Population Dynamics of the Cockle Cerastoderma edule in a Southern North Sea Tidal Flat System

The cockle Cerastoderma edule is one of the most common macrofauna species in the Wadden Sea areas of the North Sea. Cockle population dynamics are influenced by various abiotic and biotic factors such as temperature, food availability, and inter- and intraspecific competition. Cockles play an important role in the food web of the Wadden Sea, for instance, large shellfish-eating birds, such as oystercatchers and common eiders, use the cockle C. edule and the blue mussel Mytilus edulis as a main diet component. However, the populations of shellfish-eating bird species have been declining dramatically across the Wadden Sea since the beginning of the 21st century. While there are detailed monitoring programs in blue mussels due to commercial interests, little information is known about the stocks and long-term dynamics of cockles in the German Wadden Sea. To fill this gap, in 2005 a local conservation society (“Der Mellumrat e.V.”) initiated a study to sample cockles at one transect per year south of the island of Mellum, which was extended by 5 more transects in 2011. In addition to the spatial analysis, we analyzed the long-term variability in cockle population dynamics. Min/max autocorrelation factor analysis (MAFA) revealed a decline in cockle abundance, while no clear length trends were found. Canonical and spearman correlation analyses exposed significant correlations between cockle abundance and length and chlorophyll a, mussel bank area as well as oystercatcher and common eider populations. This study clearly shows that there is an urgent need for comprehensive time series of cockle data to analyze and explain ecological long-term changes in cockle population dynamics in relation to environmental changes and to point out how parts of the Wadden Sea food web, such as shellfish-eating birds are affected by these changes.


INTRODUCTION
Marine ecosystems, especially shallow coastal ecosystems such as the Wadden Sea have undergone drastic changes during the last century (e.g., Beukema, 1992;Kröncke et al., 2001;Meyer et al., 2016). Changes were caused by several factors including the anthropogenic induced sea surface temperature (SST) increase of 1.5-2.0 • C since 1950 (Beare et al., 2002;van Aken, 2008;Beukema and Dekker, 2020b), increasing anthropogenic land use due to commercial interests in coastal regions (Reise, 2005;Baer et al., 2017;Schultze and Nehls, 2017), or changes in food availability through re-or eutrophication van Beusekom et al., 2019). Abiotic changes are followed by secondary effects such as habitat loss or changes in food webs (Lotze et al., 2005;Reise, 2005;Harley et al., 2006). Due to their relatively sessile habit and thus, limited ability to move, benthic species are valuable indicators in long-term studies for environmental changes and disturbances (Schückel et al., 2015;Ghodrati Shojaei et al., 2016;Clare et al., 2017;Beukema and Dekker, 2020b).
One of the most common macrofauna species of the Wadden Sea is the cockle Cerastoderma edule (Riesen and Reise, 1982;Flach, 1996). As suspension feeder, cockles are living in the upper layers of the sediment. As Flach (1996) already stated, an adult cockle of around 3 cm length is an excellent bioturbator that can move around 20 cm 2 sediment to a depth of 3 cm per week. The biomass of cockles made up to 16% of the total benthic biomass of the Balgzand (Beukema et al., 1993;Beukema and Dekker, 2006;Compton et al., 2013). During the second half of the past century, in the Dutch Wadden Sea abundances of 100-500 adult cockles (∼2-3 cm) per 1 m 2 were found, while locally densities can reach several thousand adult cockles (Kristensen, 1958;Beukema, 1982;Jensen, 1992;Flach, 1996). Cockles are widespread around the coasts of Europe and tolerate a wide temperature range (Malham et al., 2012), but they are vulnerable to low winter temperatures and do not tolerate ice cover and frost in winter (Strasser et al., 2001). During cold or ice winters, cockles are strongly affected (Beukema et al., 2010), which led to a stock reduction of around 80%. The cold winters are usually followed by a strong recruitment and recovery until autumn (Jensen, 1992;Strasser et al., 2001Strasser et al., , 2002. However, Magalhães et al. (2016) stated, that temperature is only one of several factors, influencing cockle recruitment success, which is highly site-dependent.
There are clear signs that benthic long-term variability in the North Sea is influenced by climate change and increasing SST Beukema and Dekker, 2020b;Moore and Smale, 2020). For some bivalve species such as the Baltic clam (Limecola balthica), strong negative effects of increasing temperature on population dynamics were found, e.g., reduction of recruitment or enhancement of seasonal weight loss, resulting in decreasing abundances and biomass (Beukema et al., 2009). In Wadden Sea areas, bivalve populations and their biomass fluctuate strongly, however, the underlying processes are largely unknown (Beukema et al., 2010). Most recently, Beukema and Dekker (2020b) analyzed effects of increasing temperatures on cockle variability since the early 1970s. In adult populations, they found no clear long-term trends, however, in hot summers the survival rate of older cockles was slightly reduced. In addition, declining recruitment was compensated by an increasing winter survival rate.
It is well known, that benthic macrofauna, living at the sediment-water interface plays a key role in the food web of the Wadden Sea ecosystem (Kröncke, 2006;Lynam et al., 2017). Macrofauna species are located at a lower level of the marine food web, building a link between pelagic and benthic structures (Graf, 1992). Specified for cockles it means on the one hand, that cockle population dynamics depends on food supply through exploitable nutrients in the water column and consequently is bottom-up controlled by primary production (Flach, 1996;Strasser et al., 2001;Beukema and Dekker, 2005;Kröncke, 2006). Other bottom-up effects controlling cockle variability are hydrodynamics, temperature, or oxygen intake. Top-down effects also affect cockle dynamics in the Wadden Sea, for instance in form of predation through epibenthic species such as shrimps or juvenile fish (Beukema and Dekker, 2005), but also negative impacts through fisheries and shellfish-eating birds are possible (Beukema and Dekker, 2006).
On the other hand, changes in cockle populations may affect lower or higher trophic levels (Schückel et al., 2015). An increase in cockle density can increase bottom roughness, leading to increased erodibility of sediment (Ciutat et al., 2007) and thus might change the environmental conditions for the accompanying fauna. Moreover, large shellfish-eating birds, such as oystercatcher and common eider, use cockles, along with the blue mussel (Mytilus edulis), as a main diet component in the Wadden Sea Beukema et al., 2010). Changes in cockle population dynamics may therefore influence food availability for bird populations (Camphuysen et al., 1996(Camphuysen et al., , 2002Beukema et al., 2010). It is therefore with concern that the populations of large shellfish-eating bird species such as common eider (Somateria mollissima) and the oystercatcher (Haematopus ostralegus) have been declining dramatically across the Wadden Sea since the beginning of the 21st century (Kleefstra et al., 2019). The herring gull (Larus argentatus) is another important and common bird species in the Wadden Sea, which frequently feeds on cockles (Becker and Exo, 1991). However, herring gulls are generalists in terms of their diet, which, next to bivalves, includes fish, crustaceans, but also recyclable food from landfills (Becker and Exo, 1991;Bukacińska et al., 1996;Dierschke, 2006). But despite their flexible diet, from 1995 onward the population of herring gulls on Mellum decreased (Thorup and Koffijberg, 2016), with the exception of only 1,228 breeding pairs in 2019 (Scheiffarth et al., 2020).
Due to the commercial interest in blue mussels, monitoring programs have been established in all Wadden Sea countries, so that a good overview of the blue mussel population and, in some cases, its growth rates and body conditions, such as weight, length or shell density, is available (Nehls et al., 2009). Since there was no commercial exploitation of cockles outside the Netherlands since 1990, little information is known about the large-scale and longterm variability of cockles, their stocks and dynamics in most areas of the German Wadden Sea (van der Graaf et al., 2009).
To fill this gap in knowledge, in 2005 the local conservation society ("Der Mellumrat e.V.") initiated this study, which is unique in the central Wadden Sea outside the Netherlands regarding its size, length, and completeness. Therefore, the study area seems suitable as a model region for further sampling on other tidal flats in the German Wadden Sea. Thus, the aim of this study was to characterize the spatial and long-term variability of cockles in the study area. To give insights in the function of cockles in the food web of the Wadden Sea, the cockle long-term variability was analyzed in relation to abiotic environmental parameters such as SST, salinity, and chlorophyll a and biotic parameters such as occurrence of shellfish-eating birds and mussel bed area.

Study Area and Sampling
Cockle samples were taken in the south of the island of Mellum in the core zone of the Lower Saxon Wadden Sea National Park, in the center of the UNESCO Natural World Heritage site. Six transects (L-Q) were sampled during low tide annually in autumn (August to October) (Figure 1). Transect M was sampled over 15 years, between 2005 and 2019, transects N, O, and P were sampled between 2011 and 2019, and transects L and Q were sampled between 2013 and 2019.
Covering a gradient from the high (HWL) to the low water line (LWL), samples were taken in 250 m steps, starting at 250 m from land ending at 1,500 m from land. At each sampling distance, six replicates were taken (Figure 2). The upper 10 cm were sieved over 1 mm mesh size. Cores with a diameter of 15 cm were used, accounting for a surface area of 177 cm 2 . In 2005 a core with a diameter of 30 cm was used, accounting for a surface area of 707 cm 2 . All data were standardized to an area of 1 m 2 , abundances are given as individuals per 1 m 2 (ind./m 2 ). The length of each cockle was measured with a caliper up to 0.1 mm (Figure 2).

Statistical Analyses
Statistical analyses were performed with R 3.6.3 1 , except for min/max autocorrelation factor analysis (MAFA) and canonical correlation analysis, which were accomplished using the software package Brodgar 2 .

Spatial Variability of Cockle Data
For spatial analysis, all data from 2005 to 2019 were used. Abundance and length data were not normally distributed (Shapiro-Wilk test of normality; p-value < 0.05) and there was no significant difference in abundance and length between the six transects (Kruskal-Wallis rank sum test; p-value > 0.05). Thus, for spatial analysis, abundance and length data of all transects were pooled.

Visualization of spatial data length
The boxplots (Frigge et al., 1989) were created, based on mean abundance and length per sampling distance and year. Absolute outliers in mean abundances >4,000 ind./m 2 were excluded from the visualization. The visualization of the size classes depending on the sampling distance was conducted in the form of histograms. Length data were clustered in 2 mm steps, revealing 21 size categories. Data are shown as relative densities (%) per sampling distance.

Variance analysis
An overall Kruskal-Wallis rank sum test was used to examine significant differences in abundance and length between sampling distances. The non-parametric Kruskal-Wallis rank sum test analyses, whether independent samples originate from a common population. It can be used to compare more than two groups, starting from the null hypothesis, that there is no difference between the groups.
A non-parametric pairwise Wilcoxon rank sum test using the Bonferroni correction, was used to test significant differences in abundance and length between each sampling distances. A nonparametric Wilcoxon rank sum test is used to test if it is equally likely that a randomly selected value from one population is greater or less than a randomly selected value from the other population. If this hypothesis is rejected, it is assumed that the values from one population tend to be larger or smaller than those from the other population.

Environmental Data Sea surface temperature anomalies
The abiotic environmental variables sea surface temperature (SST in • C) and salinity were measured near the island of Mellum (53 • 42 50 • N, 08 • 06 45 • E), monthly since 2005. In detail, anomalies of mean annually SST and salinity, mean sampling SST and salinity (mean of months August to October), and mean winter SST and salinity (mean of months January to March) were used. Cold winters were defined as winters with a SST anomaly > -2 • C, warm summers were defined as summers with a SST anomaly > + 2 • C. Cold winters with a SST anomaly > -2

Chlorophyll a
Chlorophyll a (µg/l) concentrations in the water column were measured irregularly in winter, summer, and autumn at 53 • 43 13.362 N and 8 • 3 43.678 E in the entrance of the Jade as part of a standard monitoring program (NLWKN, 2013). To increase sample numbers and comparability, yearly means were used for the analysis.

Mussel beds
Data on eulittoral mussel beds were provided by NLPV (Lower Saxon Wadden Sea National Park Authority). The spatial distribution and size of eulittoral mussel beds was determined by aerial photographs (Folmer et al., 2014). The eulittoral mussel beds in the Wadden Sea area of Lower Saxon are characterized by blue mussels (Mytilus edulis) and Pacific oysters (Magallana gigas; formerly Crassostrea gigas).

Shellfish-eating birds
Bird data were provided by NLPV and NLWKN (Lower Saxon Department for Water, Coastal and Nature Conservation), oystercatcher and herring gull counts were conducted by a local conservation society ("Der Mellumrat e.V."). Common eiders were counted from a plane by the NLPV.  The number of common eiders (S. mollissima) in the area around Mellum are given for summer and winter. Common eiders are counted twice a year from the air: in winter (January/February) and in summer, during molting season (July/August).
Values for the oystercatcher (H. ostralegus) population were obtained from biweekly ground-based counts. They were analyzed by bird year, starting in July, and ending in June of the following year. Given is the average monthly population over the year, which is a measure of the average use of the roosting area.
Breeding numbers of European herring gulls (L. argentatus) were identified by sampling on random areas distributed over the island area used for breeding. For this purpose, the breeding area was divided into different density classes (low, medium, and high), a 50 × 50 m grid was laid over the breeding area and a total of 50 random plots were selected according to the area shares of the three density classes. Since 1993, the breeding numbers have been recorded using this method every two years from mid to late May (Wilkens and Exo, 1998).

Correlation Analysis
A correlation analysis, using the spearman correlation coefficient, between all environmental parameters and between mean abundance and length data of transect M between 2005 and 2019 was conducted. Data were visualized with a corrplot (R package "corrplot"). R-values of significant correlations (pvalue < 0.05) were plotted, higher R-values were plotted in a more intense color.

Long-Term Analysis
There was no significant difference in abundance and length between the transects, thus, only transect M was used for a longterm analysis between 2005 and 2019. Using only transect M prevents a sampling bias of the data due to a higher sample size since 2011.

Min/max autocorrelation factor analysis (MAFA)
A MAFA is a type of principal component analysis (PCA) for time series. The MAFA-axis represents the autocorrelation of a variable within a time lag k (k = 1, 2, . . .). Trends in data relate to highest autocorrelation within time lag 1. The 1st MAFA-axis reflects the most common pattern reflecting the trend of most of the variables in the time series. The 2nd MAFA-axis reflects the second most important trend in time series, reflecting a minor part of variables. One MAFA was performed for abundance data and one for length data. As longterm variables abundance or length data of the six sampling distances were used. Due to the limited number of variables, only the first two MAFA-axis are shown, but only the first MAFA was evaluated.

Canonical correlation analysis
Canonical correlation analysis was implemented in the MAFA. Canonical correlations were used to identify significant relationships between MAFA axis and response variables (abundance and length data), and further between trends and explanatory variables (environmental data). Standard deviations were used for an easier interpretation of the estimated regression parameters.

Spatial Variability in Cockle Abundance and Length per Sampling Distance
The Kruskal-Wallis rank sum test revealed significant differences between the sampling distances of the pooled transects for abundance (p-value < 0.0001) and length (p-value = 0.0069) (p-value < 0.0001). Overall, using data of all samples a mean abundance of 693 ind./m 2 was found.
Median cockle abundance increased from HWL to LWL until the sampling distance of 750 m, where highest median of mean abundances of 687 ind./m 2 , but also highest interquartile range were found. From sampling distance 750-1,500 m median abundance but also the interquartile range decreased (Figure 4). In detail the Wilcoxon rank sum test revealed significant differences for greater distance of the sampling distances, while e.g., between 250 and 500 m or between 500 and 750 m no significant differences were found ( Table 1).
Median length was nearly stable, around 10-15 mm, except of sampling distance 1,000 m, where median length lay at about 18 mm. However, the interquartile range in length increased from HWL to LWL (Figure 4). Between sampling distances 750 and 1,000 m, Wilcoxon rank sum test revealed a significant difference between sampling distance 1,000 m to all lower distances ( Table 1).

Spatial Variability of Cockle Size Classes
Overall cockle length ranged from 0.7 up to 42.0 mm. Relative densities of size classes of sampling distances 250 and 500 m were rather similar, with a maximum between 0 and 10 mm. While at sampling distance 250 m highest density of size class 4-6 mm was found, at sampling distance 500 m highest density of size class 6-8 mm was found. In contrast to the assumption that the size classes with higher relative density increased toward the LWL, at sampling distance 750 m a nearly uniform distribution of size classes <10 mm were found. And even at sampling distance 1,000 m highest densities of the smaller size class 4-6 mm were found. At sampling distance 1,250 m highest densities of size classes measured between 20 and 30 mm, at sampling distance 1,500 m between 6 and 8 mm and between 24 and 28 mm were found. Only at sampling distance 1,500 the full-size range between 1 and 42 mm was found (Figure 5).

Correlation Analysis
Within the environmental data negative R-values were found between chlorophyll a and common eider in winter (Rvalue −0.71), while between chlorophyll a and mussel bed area positive values were found (R-value 0.71). Besides, a highly negative correlation between herring gull breeding pair numbers and SST during sampling was found (R-value −0.73). Lastly, correlation analysis revealed a high correlation between oystercatcher and common eider in summer (Rvalue 0.63).
A negative correlation (R-value −0.58) was found for cockle abundance and winter SST, while cockle and oystercatcher (Rvalue −0.57) and common eider (R-value −0.6) abundance were highly negative correlated. Cockle length was negatively correlated with chlorophyll a (R-value −0.56), while significant positive correlations were found between common eider abundance in winter (R-value 0.73) and oystercatcher abundance (R-value 0.52). Lastly, a negative correlation between cockle abundance and length was found (R-value −0.65) (Figure 6).

Long-Term Trends in Mean Abundance and Length of Cockles
Considering boxplots for visualization of long-term data, no clear trends in mean abundance and mean length can be identified (Figure 7).

Long-Term Variability of Cockle Size Classes
From 2005 to 2019, cockle length ranged from 0.7 up to 40.0 mm at the transect M. No clear trend can be identified over the whole period. Over the years, the size class with the highest relative densities is that between 0 and 10 mm, followed by the size class 10-20 mm. However, while between 2005 and 2012 the size class 20-30 mm was found regularly in higher proportions, while between 2013 and 2019 this size class decreased considerably. In particular, the years 2011, 2013, and 2018 stand out, in which over 50% of the relative density was in the 0-10 mm size class (Figure 8).

Long-Term Trends and Explanatory Variables
A MAFA was conducted for abundance and length data, separately. For abundance data for the first MAFA a high autocorrelation coefficient (AC) of 0.831 and a p-value of 0.089 was found. For length data, a high AC of 0.707, but a p-value of 0.463 was found. Both long-term datasets showed high AC, however, both trends were not significant due to a p-value > 0.05. For both, length, and abundance the 2nd MAFA axis showed low AC and no significant p-values (Figure 9). Abundance AC lay > 0 from 2005 to 2011, while from 2012 to 2018 negative AC values were found, which resulted in an overall decreasing trend for the 1st MAFA axis from 2005 to 2019. Contrastingly, for length data, no clear trend is visible in the data. In 2010, a maximum AC value of 0.53 was found, followed by a minimum of −0.34 in 2014 (Figure 9). Canonical correlations revealed a significant positive correlation of the 1st MAFA axis and herring gull abundances and a significant negative correlation with chlorophyll a content and mussel bed area. Oystercatcher abundances and salinity were significant correlated with the 2nd MAFA axis (Figure 9). The 1st MAFA axis of length data was negatively correlated with chlorophyll a content and mussel bank area. More significant correlations were found between the 2nd MAFA axis sampling salinity, herring gull abundance as well as for chlorophyll a and mussel bank area (Figure 9).

Spatial Patterns of Cockle Populations South of the Island of Mellum
The spatial variability of macrofauna species is mostly controlled by the abiotic parameters sediment structure, hydrodynamics, and food availability (Rosenberg, 1995;Kröncke, 2006;. The sediment composition of the study area south of the island of Mellum is not completely homogeneous. From sampling distance 250-750 m, and at sampling distance 1,500 m, sediments are sand dominated with mud contents <12%, classified as slightly muddy sand mixed flats (Flemming, 2012) sandflat. Sampling distances 1,000 and 1,250 m are characterized by a higher a mud content between 12 and 16% (pers. comm NLPV). They are populated or affected by pacific oyster reefs. Since the sediment composition was recorded not regularly, it cannot be excluded that samples were occasionally taken in different sediment facies.
In their study Magalhães et al. (2016) pointed out a global average of 167 ind./m 2 , while they found between 1 and 2390 ind./m 2 in medium sands on the southwest coast of France. In the Jade Bay, close to the study area,  identified C. edule as a characteristic species for mixed sediments. In the 1970s, they found abundances up to 16 862 ind./m 2 on mixed sediments. On sandflats generally lower abundances around 100 ind./m 2 and below 500 ind./m 2 were found (Sanchez-Salazar et al., 1987;Ramón, 2003;. On mudflats, abundances of 1,041 ind./m 2 were found. Contrastingly in 1974, Reise (1978), described around 100 ind./m 2 on mudflats as dense populations in tidal flats of the island of Sylt. In this study, with a total mean abundance of 693 ind./m 2 at the island of Mellum, densities lay clearly below characteristic abundances of mixed flats, but above characteristic abundances of sandflats and also above the global average. However, in most of the samples cockle abundances lay clearly below 500 ind./m 2 , which is characteristic for sandflats and corresponds to results of studies on spatial cockle variability at different areas (Reise, 1978;Sanchez-Salazar et al., 1987;Ramón, 2003;. Sand flats are characterized by a hydrodynamic gradient, with higher hydrodynamic stress, coarser sediments, and low food supply at the HWL, changing gradually toward the LWL (Sanchez-Salazar et al., 1987;Kröncke et al., 2018). Suspension feeding species such as cockles, occurring in highest abundances at areas with intermediate submersion time, ensuring enough time to ingest food and low bottom hydrodynamic stress (Jensen, 1992;Herman et al., 2001;Van Colen et al., 2010). At the island of Mellum, highest mean cockle abundances were found at 750 m from the HWL, characterized by intermediate hydrodynamic stress and submersion time, while length increased from HWL to LWL. These gradients in mean cockle abundance and length from HWL to LWL, indicate a hydrodynamically regulated smallscale spatial cockle variability, which correspond to results of Sanchez-Salazar et al. (1987) or Van Colen et al. (2010).
In contrast to generally higher abundances at intermediate distances from the HWL, there are differences in spatial variability of juvenile and adult cockles (Sanchez-Salazar et al., 1987;Bouma et al., 2001;Malham et al., 2012). Adult cockles >20 mm occurred all over the spatial gradient from HWL to LWL, but in a higher relative density and in higher abundances at the more distant sampling distances (Sanchez-Salazar et al., 1987;Jensen, 1992;Malham et al., 2012). Minimum size at first maturity of cockles is between 12 and 14 mm (Malham et al., 2012). Thus, in the study area reproductive cockles occurred in higher relative densities at sampling distances 1,250 and 1,500 m. FIGURE 6 | Spearman correlations coefficients between environmental parameters and biological data. SST, mean annual sea surface temperature anomaly; SST_SAMP, mean annual sea surface temperature anomaly during sampling; SST_W, mean annual sea surface temperature anomaly in winter; SAL, mean annual salinity anomaly; SAL_SAMP, mean annual salinity anomaly during sampling; SAL_W, mean annual salinity anomaly in winter; Chl_a, chlorophyll a; Eider_W, abundances of common eider in winter; Eider_S, abundances of common eider in summer; Mussel bed, mussel bed area; Herring Gull, abundances of herring gulls at Mellum; Oystercatcher, oystercatcher abundances; Abundance, mean cockle abundance; Length, mean cockle length.
According to Sanchez-Salazar et al. (1987), juvenile cockles with shell length of 1-4 mm occurred in highest abundances at distances between 300 and 800 m from the high-water line. At the island of Mellum, we found highest abundances of smaller size classes <10 mm at sampling distances up to 1,000 m from HWL, while highest relative densities of the size class 0-2 mm were found at sampling distances 250 and 500 m, which corresponds to previous mentioned results. The highest relative density of the size class 4-6 mm can be attributed to a late sampling (end of October in 2013), which was probably connected with a settlement and dispersal of early cockle recruiters (Legendre et al., 1997;Norkko et al., 2001) after an early spatfall in October. Bouma et al. (2001) studied the spatial patterns of early cockle recruiters (<1 mm) and found highest early recruitment near the HWL, which corresponds to results of Jensen (1992), who found highest densities of recruits in the highest intertidal levels (50 and 100 m from HWL). However, even if early recruiters were not recorded and the first sampling point lay at 250 m from the HWL, the results of our study in spatial variability of cockle abundance and size classes are congruent with results of studies of Sanchez-Salazar et al. (1987);Jensen (1992), and Bouma et al. (2001).

Cold Winter Effects on Cockle Long-Term Population Dynamics
From 2005 to 2019, no cold or ice winters comparable to the ice winter in 1995/1996 (Strasser et al., 2001) occurred. However, the winter 2009/2010 (in the following winter 2010) was reported as a cold extreme in an ongoing warming climate (Cattiaux et al., 2010;Wang et al., 2010). SST data from Mellum revealed the winters 2006, 2013, and 2018 as similar cold as the winter 2010 with an SST anomaly >2 • C.
In the study area, a higher variability but also generally higher abundances and higher relative densities of smaller size classes in years with preceding cold winters, e.g., 2006 and 2018, were found. Our results correspond to results of previous studies and from other areas that cold winters are normally followed by a strong recruitment and recovery until autumn (Jensen, 1992;Honkoop and Van der Meer, 1998;  Strasser et al., 2001Strasser et al., , 2002. It is noticeable that in years with cold winters, the variability in size classes is much lower and that the highest relative densities lay distinctly below a median length of 10 mm. Beukema and Dekker (2020b) attributed an enhanced reproduction more to the reduction of predators, which feed on smaller size classes Dekker, 2014, 2020b). This is probably supported by a significant negative correlation between cockle abundance and winter SST, while length data are not significantly correlated with any SST data.
Contrary to studies on effects of the cold winter 2010 on the Wadden Sea ecosystem (Büttger et al., 2011), our results of 2010 did not correspond to above described results. In 2010, intermediate abundances, similar to years with regular winters, were found. In addition, highest median length data, with a uniform size distribution and a strong 20-30 mm size class were found. A comparable increase in the abundances and the relative densities of the size classes <10 mm became apparent in 2011. In contrast to the other years, the winter of 2010 was followed by a summer with comparably low SST. A combination of low winter SST but also low summer SST could have caused effects in the size classes only becoming visible in the following year 2011. In addition, Meyer et al. (2016) found constantly increasing abundances of the shore crab Carcinus maenas in the Jade area since the late 1990s. Shore crabs are next to brown shrimps and demersal fish species FIGURE 9 | Results of min/max autocorrelation factor analysis (MAFA) and canonical correlations. SST, mean annual sea surface temperature anomaly; SST_SAMP, mean annual sea surface temperature anomaly during sampling; SST_W, mean annual sea surface temperature anomaly in winter; SAL, mean annual salinity anomaly; SAL_SAMP, mean annual salinity anomaly during sampling; SAL_W, mean annual salinity anomaly in winter; Chl_a, chlorophyll a; Eider_W, abundances of common eider in winter; Eider_S, abundances of common eider in summer; Mussel bank, mussel bank area; Herring Gull, abundances of herring gulls at Mellum; Oystercatcher, oystercatcher abundances; Abundance, mean cockle abundance; Length, mean cockle length. main predators of juvenile cockles (Sanchez-Salazar et al., 1987;Malham et al., 2012;Whitton et al., 2012). Thus, cockles are exposed to a higher predatory risk. Meyer et al. (2016) found exceptionally high abundances of the predatory shore crab Carcinus maenas in the study area in 2011, which contradicts the before mentioned results by Dekker (2014, 2020b) and which may have caused a reduction in cockle stocks in 2010 in contrast to 2006 or 2013, when low abundances of C. maenas were found.

Cockle Long-Term Variability in Relation to Environmental Parameters
Looking at the boxplots and histograms, no clear longterm trends can be identified and the MAFA revealed no statistically significant trends ("p-value ≤ 0.05") for abundance and length data. At first, this seems not surprising taking the high annual variability described by several studies into account (Beukema et al., 2010;Beukema and Dekker, 2020b). However, considering the results of the MAFA, a clear long-term trend can be seen (Figure 9). Therefore the question arises, whether the lack of statistical significance also means the lack of trends or ecological interpretable results (Wasserstein and Lazar, 2016;Amrhein et al., 2019)? The lack of statistically significant trends may be attributable to the characteristics of environmental studies. In case of the present study due to incomplete or infrequent sampling, interfering values due to inconsistent sediment structure, as mentioned before, or a different sampling time in some years, where e.g., late samplings in October may lead to higher abundances.
There are some macrozoobenthic long-term sampling stations in the Wadden Sea (Drent et al., 2017), however, these stations are not very suitable for the evaluation of cockle population dynamics. This study is a first approach in analyzing cockle long-term variability in the German Wadden Sea, due to the lack in commercial interest nearly no long-term data are available and only little information is known (van der Graaf et al., 2009). A lack in evaluable cockle longterm and large-scale data means that information of an essential ecological parameter in the Wadden Sea is missing. Cockle long-term and large-scale data could not only be used to evaluate variability in cockle populations in relation to environmental parameters, such as SST and chlorophyll a, but also to fill gaps in the analysis of long-term variability of the persistently decreasing stocks of shellfish-eating birds such as oystercatchers and common eiders (Ens, 2006;Beukema et al., 2010).
During the last decades significant changes and shifts in macrofauna communities in the near-and offshore areas of the south-eastern North Sea Meyer and Kröncke, 2019) and the Wadden Sea (Schückel et al., 2015;Meyer et al., 2016) occurred. These changes and shifts were linked with changes in the climate system and an increase in SST since the beginning of the 21st century. Due to the increase in SST of 1.5-2 • C in the Wadden Sea during the last 50 years (van Aken, 2008; Beukema and Dekker, 2020b), we cannot make a statement for exceptional hot summers, such as described by Beukema and Dekker (2020b) because only in the years 2006 and 2014 the mean summer SST anomaly lay over 1.5 • C but under 2 • C, which is not exceptional hot in relation to the continuously increasing SST.
In our study area, the decreasing trend in AC values of the MAFA indicate overall decreasing trends in cockle abundance since 2005. For length data, maxima in 2010 and 2018 and a minimum in 2014 were found. Statistical long-term analyses revealed no significant correlation of cockle long-term trends in the study area with SST data. As mentioned before, intertidal sandflats in the Wadden Sea are highly variable systems, characterized by strong hydrodynamic stress e.g., through tidal currents, wave impact of local wind fields, and a high seasonal temperature variability. Further, several studies described the high temporal variability of cockle population dynamics, caused by recruitment failures, larval and juvenile dispersal, and feeding and foraging pressure (Norkko et al., 2001;Beukema et al., 2010;Beukema and Dekker, 2020b). Thus, indications of increasing SST and climate changes might be overlaid by a high environmental and biological variability. Further, a high spearman correlation between winter SST and cockle abundance was found, which could indicate, that the absence of cold winters comparable to those of 1995/1996, may have caused regularly recruitment failures (Beukema and Dekker, 2005). Simultaneously, an increasing predatory risk through shore crabs, an interspecific competition with mussel bed areas and a lack in usable food supply seem to have a negative impact on cockles and, thus, declining trends in cockle abundance emerge.
Canonical correlations between environmental parameters and MAFA axes revealed highly negative correlations of cockle abundance and length with chlorophyll a content in the area around Mellum, similar to the spearman correlation. This result seems contradictory to many studies, that found positive correlations between chlorophyll a concentrations and cockle biomass and showed a clear food dependence of the suspensionfeeding cockle (Prins and Smaal, 1989;Kamermans, 1993). Higher chlorophyll a concentrations led to a higher growth performance because chlorophyll a is an important indicator of food availability (Wijsman and Smaal, 2011). Not only the quantity of chlorophyll a, but especially the quality and also the concentration of all particles in the water column play a role for the usability. A higher amount of suspended matter possibly makes the chlorophyll a less usable (Kamermans, 1993). For other species, such as the deposit feeding sea urchin Echinocardium cordatum, similar relationships were found, where abundance depends on food quantity, whereas the food quality influenced the size of individuals (Wieking and Kröncke, 2003). Additionally, constantly decreasing nutrient input from rivers and thus decreasing chlorophyll a content in the Wadden Sea (Beukema and Dekker, 2020a) may increase the intraspecific food competition. Furthermore, there is a negative impact at the presence of blue mussels and pacific oysters for cockles by reducing the food supply due to a reduction in chlorophyll a. Comparing blue mussels and cockles, blue mussels are capable for sorting particle prior the ingestion (Prins and Smaal, 1989;Kamermans, 1993;Smaal and Haas, 1997). This is confirmed by a highly positive spearman correlation coefficient between mussel bed area and chlorophyll a in combination with a negative spearman correlation coefficient between mussel bed area and cockle length and between length and abundance MAFA axes and mussel bank area.
Shellfish-eating birds are important predators for cockle population in the Wadden Sea (Ens, 2006;Piersma, 2006;Beukema et al., 2010). We found no significant canonical correlations between trends of the 1st MAFA axes and shellfish-eating birds, which indicates that there is no predation pressure by shellfish-eating birds and long-term cockle variability. However, spearman correlation revealed significant negative correlations between common eider and oystercatcher populations and cockle abundance and significant positive correlations with cockle length. The high correlation confirms that shellfish-eating birds are important in structuring cockle populations (Sanchez-Salazar et al., 1987;Nehls et al., 1997;Piersma, 2006;Malham et al., 2012). Our results indicate that feeding pressure by common eiders and oystercatchers reduces the abundances, which seems to have positive effects on the length.
Our results indicate that cockles might be under pressure by low recruitment after warm winters, food competition with pacific oysters and blue mussels and predation by oystercatchers and common eiders.

DATA AVAILABILITY STATEMENT
Data for common eider number and mussel bed distribution are fully available at: http://mdi.niedersachsen.de/HeronKaDI/ JAVA_SCRIPT/37_Portal/. All other data are fully available on request. Requests to access the datasets should be directed to GS (Gregor.Scheiffarth@nlpv-wattenmeer.niedersachsen.de).