Connecting foraging and roosting areas reveals how food stocks explain shorebird numbers

Shorebird populations, especially those feeding on shellfish, have strongly declined in recent decades and identifying the drivers of these declines is important for conservation. Changing food stocks are thought to be a key driver of these declines and may also explain why trends have not been uniform across Europe’s largest estuary. We therefore investigated how winter population trends of Eurasian oystercatchers (Haematopus ostralegus) were linked to food availability in the Dutch Wadden Sea. Our analysis incorporated two spatial scales, a smaller scale focused on roost counting areas and food available to birds in these areas and a larger spatial scale of tidal basins. A novelty in our study is that we quantify the connectivity between roosting and foraging areas, identified from GPS tracking data. This allowed us to estimate food available to roosting birds and thus how food availability may explain local population trends. At the smaller spatial scale of roost counting areas, there was no clear relationship between available food and the number of roosting oystercatchers, indicating that other factors may drive population fluctuations at finer spatial scales. At the scale of tidal basins, however, there was a significant relationship between population trends and available food, especially cockle Cerastoderma edule,. Mortality and recruitment alone could not account for the large fluctuations in bird counts, suggesting that the site choice of wintering migratory oystercatchers may primarily drive these large fluctuations. Furthermore, the relationship between oystercatcher abundance and benthic food stocks, suggests winter shorebird counts could act as ecological indicators of ecosystem health, informing about the winter status of food stocks at a spatial scale of tidal basins.


Introduction
Intertidal areas support a rich benthic fauna across the world and are a vital habitat for staging, breeding and wintering (migratory) waterbirds (Hua et al., 2015;Blew et al., 2017). Conserving intertidal habitats is consequently vital for protecting migratory waterbirds since the loss of habitat in either of these life stages, i.e. staging, breeding or wintering, could have population-level consequences (Stralberg et al., 2011;van Roomen et al., 2012;Hua et al., 2015). The Wadden Sea is an example of an important intertidal area for waterbirds; it is the largest coherent intertidal area in the world, a World UNESCO heritage site in Denmark, Germany and The Netherlands, a Natura2000 site of international importance, and designated under the RAMSAR convention (Common Wadden Sea Secretariat (CWSS), 2017). The macrozoobenthos community of the Wadden Sea has been changing in recent decades due to factors including climate change, eutrophication, exotic invasive species, and human activities like fisheries and construction of sea defences (Beukema and Dekker, 2020a). Concurrently, many waterbirds that depend on the intertidal flats of the Wadden Sea have experienced strong declines (van Roomen et al., 2012;Blew et al., 2017), especially bivalve specialists (van Roomen et al., 2012).
Several potential causes for the decline of waterbirds dependent on intertidal areas like the Wadden Sea have been suggested (van Roomen et al., 2012). A prevailing hypothesis is that shorebird population sizes are linked to benthic food stocks, and especially for species specialising on bivalves (Atkinson et al., 2003(Atkinson et al., , 2010Verhulst et al., 2004;Kraan et al., 2009). Studies have shown that over-winter survival is strongly dependent on available food stocks, especially during periods of severe weather (Camphuysen et al., 1996;Atkinson et al., 2003). Research has also shown that the effects of low food availability and quality in combination with severe weather may particularly influence bivalve specialists (Camphuysen et al., 1996;Schwemmer et al., 2014). However, many waterbirds that utilise intertidal areas are highly mobile and local declines of benthic prey may result in predator redistribution patterns to other foraging areas (Atkinson et al., 2003), as exemplified by how the winter distribution patterns of Common Eider (Somateria mollissima) appear to be linked to available food stocks (Cervencl et al., 2015). Stillman et al. (2010) also highlighted how a decline in population at a local site may not necessarily be due to in situ declining benthic fauna, but instead due to an improving situation in neighbouring areas. Studies that relate shorebird populations to benthic fauna should therefore not only consider benthic fauna at local sites, but also benthic stocks available in the wider landscape.
Given that shorebird population sizes may be related to benthic food stocks, and that waterbirds may act as indicators of environmental health (Gregory and van Strien, 2010), understanding the relationship between benthic fauna and shorebird abundance has become a topic of importance (Folmer et al., 2010;Ponsero et al., 2016;Horn et al., 2017). Substantial investments are made every year to measure and predict benthic fauna of intertidal areas (Ysebaert and Herman, 2002;Bijleveld et al., 2012;Compton et al., 2013). A challenge that arises in relating avian population trends to benthic stocks is that bird counts are usually performed at high tide roosting sites (Hornman et al., 2020); from these aggregation sites shorebirds may disperse over a large area to forage (Schwemmer et al., 2016;Enners et al., 2019). Therefore, the connection between foraging areas and different roost sites needs to be understood before population sizes can be related to benthic food stocks. Furthermore, shorebirds may not use all suitable foraging areas with sufficient biomass of benthic fauna (Bijleveld et al., 2016). However, identifying foraging areas used over time, and subsequently estimating the connectivity between roosting and foraging areas has only become possible in recent years with the advent of GPS tracking technologies (Schwemmer and Garthe, 2011;Dokter et al., 2017).
The Eurasian oystercatcher (Haematopus ostralegus) is an example of a shorebird for which the global population is heavily reliant on the Wadden Sea during summer and especially during winter, when bivalves are an important food source (van de Pol et al., 2014). The oystercatcher has experienced some of the strongest declines of waterbirds in the Wadden Sea (van Roomen et al., 2012;Blew et al., 2017). It is currently listed as near-threatened by the IUCN (BirdLife International, 2021), and multiple threats have been described for the species both within intertidal areas and inland breeding areas (van de Pol et al., 2014). One of the key threats, especially during winter, has been the decline in food sources and associated reduction in over-winter survival (Camphuysen et al., 1996;Verhulst et al., 2004;Allen et al., 2019).
The oystercatcher is an example of an interference sensitive species (Ens and Crayford, 1996) meaning that their distribution is expected to be better predicted by benthic biomass than interference insensitive species that are not only attracted to food, but also to conspecifics (Folmer et al., 2010). Our aim is therefore to determine whether changes in counts of the Eurasian oystercatcher at roosting sites can be explained by changes in the stocks of available benthic prey. Instead of a general measure of benthic biomass in a region, which may include areas never used by oystercatchers, a novelty is that we specifically quantify the interconnectivity between areas used during roosting and foraging, as identified from GPS-tracked individuals,. Without site-specific knowledge of different areas, one may expect oystercatchers to be drawn to regions with high food stocks and not too many competitors (Sutherland, 1983;Van Der Meer and Ens, 1997), and to subsequently distribute evenly among high tide roosts to access these food resources. However, a number of factors other than benthic prey stocks may also influence the number of oystercatchers at a high tide roost, such as distance to foraging areas, predation risk, disturbance and shelter from weather (van Gils et al., 2006;Rosa et al., 2006;Folmer et al., 2010;Horn et al., 2020), meaning that trends at local spatial scales may be more driven by other processes and more difficult to link to food resources than those at regional scales. We therefore conducted our analysis at two spatial scales to determine if changes in prey biomass led to changes in the population size of oystercatchers at local and regional scales. The larger spatial scale focused on regions comparable in size to tidal basins (sensu van Roomen et al., 2012), which we hereon refer to as "tidal basins". The smaller spatial scale focused on specific high tide roost counting areas, which we hereon refer to as a "local" spatial scale.

Study system
The Eurasian oystercatcher is a long-lived shorebird with average adult annual survival of ~90% (Allen et al., 2019). Breeding site fidelity is considered high, and fidelity to overwintering sites is also high, especially for adults, although individuals may move to other regions during extreme cold spells (Camphuysen et al., 1996). Individual oystercatchers can specialize on different prey items, which in turn lead to a particular bill shape (Swennen et al., 1983;Hulscher, 1985). Therefore, although cockles and mussels are commonly considered as principal prey species during winter (Camphuysen et al., 1996;Zwarts et al., 1996), a broad array of benthic fauna can be considered as available.
The oystercatcher is distributed along the coastline of Europe and may also breed up to a few hundred kilometres inland . The Wadden Sea supports ~50% of the Eurasian oystercatcher global winter population (van de Pol et al., 2014;Blew et al., 2017). Maximum estimates of the Wadden Sea population have declined from ~570,000 oystercatchers in the early 1990s to ~410,000 oystercatchers in 2013/14 (Blew et al., 2017). This study was conducted in the intertidal Dutch Wadden Sea area (Fig. 1), which contains 61% of the Dutch wintering population and 13% of the global Eurasian oystercatcher population (van de Pol et al., 2014). The study focuses on four regions of the Dutch Wadden Sea, namely Balgzand, Vlieland, Ameland and Schiermonnikoog (Fig. 1). These four regions approximate the tidal basins previously used to quantify spatial differences in shorebird population trends with regard to size (van Roomen et al., 2012), and although technically different, we refer to them as tidal basins.

Methodological framework
Our analysis focused on the winter period, which is an energy demanding period when oystercatchers largely depend on the intertidal areas for foraging (Goss-Custard et al., 1996a;van der Kolk et al., 2020), and the oystercatcher population in the Wadden Sea drastically increases with wintering migrants (van de Pol et al., 2014), meaning that competition for food is greater and thus competition-sensitive species like the oystercatcher can be expected to distribute according to available resources (Goss-Custard et al., 1996b;Folmer et al., 2010). We first identified where oystercatchers roost and forage using GPS tracking data. We then used the individual-level movement data to quantify the connectivity between roosting and foraging areas. Food availability for birds at each roost was then estimated from the amount of benthic biomass in connected foraging areas. We subsequently investigated how food availability (i.e. benthic prey stocks) in these areas were correlated to counts of roosting oystercatchers at the two spatial scales described.

Roosting and foraging areas
Movement data from several tracking studies of oystercatchers equipped with UvA-BiTS Global Positioning System (GPS) transmitters (Bouten et al., 2013) were combined. The studies were conducted between 2008 and 2017 in four tidal basins, namely Balgzand (Dokter et al., 2017), Vlieland (van der Kolk et al., 2020), Ameland and Schiermonnikoog   (Fig. 1). GPS tracks were standardised to contain a minimum of one-hourly intervals and 1000 data points. The final dataset consisted of 73 individuals (Balgzand = 17, Vlieland = 34, Ameland = 12, Schiermonnikoog = 22). GPS locations of oystercatchers were assigned as roosting or foraging based on the timing in the tidal cycle (Roosting: 4.5 h after low tide to 4.5 h before the next low tide; Foraging: 1.5 h before low tide to 1.5 h after low tide; Appendix A). Although accelerometery data was available to estimate behaviour, our method of tidal classification had an 87% agreement with the accelerometer data in classifying roosting and foraging behaviours (Appendix A: Table A1). Using the accelerometer data would also result in 13 fewer birds for the analysis (Appendix A).
We used a kernel-based approach (Calenge, 2006) to obtain non-overlapping utilisation distributions for roosting areas and foraging areas for each individual (see details in Appendix A). Next, principal roosting and foraging areas were identified by overlapping the individual-level roosting and foraging polygons used by at least three individuals (Appendix A). We then identified which foraging polygons were used by oystercatchers from a particular roost and summarized the result in a connectivity table. In this calculation, foraging GPS points were connected to the roost that an oystercatcher visited longest between 10 h before and after that foraging GPS point (hereon roost-forage connection), thus covering both the preceding and upcoming high tide. Based on these connections, frequency tables (Ft) were constructed for all individuals i by summing all roost-foraging area connections ( ∑ Ft i ). The frequency tables were converted to individual-level proportion tables (Pt) so that all individuals summed to one ( The individual-level proportion tables were subsequently summed to obtain a population level proportion table (Pt = ∑ i Pt i ). The proportional connectivity in this table was thus unaffected by individual differences in the amount of GPS data available. The connectivity table was obtained by normalizing the connections such that the connectivity between a roost and all available foraging areas summed to one (C r = Pt r / ∑ Pt r ; r = roost). The resulting connectivity table was used in all further calculations with values ranging from 0 to 1, a 1 indicated that all birds using a particular roost foraged at one particular foraging area, whilst a 0 indicated that no birds from the roost used a particular foraging area.

High-tide roost counts and trends
High tide roost count data between 2009 and 2017 (Hornman et al., 2020) were used to estimate population fluctuations in roost counting areas. The high tide roost counts were conducted within a roost counting area (Fig. 1) and most of the areas were counted during four Wadden Sea-wide counts in September, November, January and one other shifting date. Additionally, local counts were often performed on other dates. The number of oystercatchers at a roost varied not only across years, but also within seasons. Since the counts were generally not at regular intervals and there were only a limited number of datapoints, smoothing splines were used to model population trends. First, the average seasonal winter pattern, α, was fitted by fitting a general additive model (GAM; with p-splines, the default value for k where this equalled the degrees of freedom, a quasi-Poisson distribution family and a log link function; Wood, 2017) disregarding yearly differences. We estimated α, the average number of birds at a roost during winter, as the area under the curve -α-divided by the number of days in the season. Roost counting areas with a low number of birds (α<1000, three of the 17 roost counting areas) were excluded from the local scale analysis. The actual roost count was divided by the average seasonal curve -αto obtain a normalized roost count (β). The annual trends (γ) at a roost counting area were estimated by fitting a second GAM (with p-splines, the default value for k where this equalled the degrees of freedom, a quasi-Poisson distribution and a log link function to the normalized values β; Wood, 2017). The population trend across years (δ) were obtained by multiplying the annual trend (γ) with the normalized number of birds (α).
We also estimated population trends per tidal basin. We summed the roost counts per tidal basin and smoothed the within-winter and amongyear variation as described for the roost counting areas.

Foraging landscape
To identify the prey available to oystercatchers, a foraging landscape was estimated based on two benthic surveys. The first, Synoptic Intertidal Benthic Survey (SIBES), is performed by NIOZ principally in June and July using a 500 × 500m grid with a broad coverage over the Wadden Sea (Bijleveld et al., 2012;Compton et al., 2013). The second survey is a shellfish survey mainly focussing on cockles (Cerastoderma edule) and mussels (Mytilus edulis), performed by Wageningen Marine Research (WMR) between April and June (Troost et al., 2021). The surveys were used to create annual estimates of available biomass per foraging area  from 2008 to 2013. Since the benthic surveys were carried out between April and June, the available prey biomass in the subsequent winter still changed due to prey growth during the preceding summer (Beukema and Dekker, 2020a)) and mortality in winter (Burdon et al., 2014;Dokter et al., 2017;Beukema and Dekker, 2020a). Given these processes, the benthic surveys may provide a better indication of the available prey biomass in the previous winter rather than the coming winter. Therefore, we compared bird counts in year t with available biomass in the previous summer (year t) and with the available biomass in the subsequent summer (year t + 1).
Available biomass included benthic fauna that were known to be prey species of oystercatchers (Arenicola marina, Cerastoderma edule, Ensis directus, Hediste diversicolor, Limecola balthica, Mya arenaria, Mytilus edulis and Scrobicularia plana; Ens et al., 2015). In our analyses, we consider the biomass of each prey species separately. Scrobicularia plana was excluded because it was rare (the rarest species) and often absent from foraging areas especially in the tidal basin Vlieland. Six percent of the total area of foraging areas had no survey data available as SIBES excluded some elevated locations from their sampling scheme. However, these areas presumably contained food since oystercatchers were foraging in these areas. We therefore imputed biomass in these areas with the mean biomass per foraging area.
Roost counting areas were generally larger than the individual roosting polygons identified from the GPS data (Fig. 1). Therefore, a counting area may contain several roosts. Roosting polygons were assigned to a roost counting area based on the amount of overlap, and if a roost counting area contained multiple roosting polygons, the connectivity values between the high tide roosts and the foraging polygons were averaged. The available biomass per roost counting area was defined as the sum of benthic biomass in foraging polygons that were connected to the roost counting area. We only considered foraging areas with a connectivity (x) greater than 0.01. Low connectivity values were indicative of foraging areas that were infrequently used from a particular roost counting area, for example due to distance or profitability. Available biomass was also calculated for each tidal basin rather than for individual roost counting areas (Fig. 1) in a similar way as the roost counting areas, except that if a tidal basin contained multiple roost counting areas, the connectivity values were averaged.
Our method contained a number of analytical steps and to determine whether these were necessary, we also considered a null model where we estimated the available biomass in a tidal basin, without taking into account where oystercatchers were foraging. We estimated the biomass within a buffer zone of the roosting areas (the same roosting areas were selected as in our analysis), the size of which was based on the 95th percentile of distances between the roosting and foraging area (7.57 km). The null model enabled us to validate the added value of identifying the available biomass in regions where oystercatchers forage, as opposed to general availability in the tidal basin.

Statistical analysis
We analysed whether spatiotemporal variation in the number of roosting oystercatchers could be explained by spatiotemporal variation in food availability. To be able to compare different sites, the available biomass and bird counts were normalized by dividing the annual values by the site mean during the study period. For both scales of analysis of 1) roost counting area and 2) tidal basin, a model selection approach was used to identify which benthic prey species explained most variation in roost counts. Year was also included as a categorical variable in the analysis of roost counting areas to be able to distinguish between birds redistributing within a tidal basin and birds leaving or entering the tidal basin, thus changing the bird numbers in all roost counting areas within the tidal basin. A global linear model was made that included known benthic prey species of the oystercatcher (described in 2.5) and model selection was performed using the MuMIn package in R (Barton, 2019). The top performing model was chosen based on the Bayesian Information Criterion (BIC; Burnham and Anderson, 2004). We checked final models for collinearity using variance inflation factors (VIFs; O'Brien, 2007) in the R package car (Fox and Weisberg, 2019), whereby VIFs >4 indicate collinearity in the model.

Roosting and foraging sites
Fifty-eight roosting polygons and 140 foraging polygons were identified from the GPS data (see foraging areas near Ameland in Fig. 2; for other tidal basins see Appendix B: Figure B1-B3). Mean flight distances between roosting and foraging polygons varied among the four tidal basins, ranging from 3.82 km (standard error = 0.50) (Balgzand) to 2.31 km (se = 0.17) (Ameland; Fig. 3). An average roosting polygon was connected to nine different foraging polygons. Foraging polygons had an average size of 0.5 km 2 with no clear differences among tidal basins. In most tidal basins, distinct patterns emerged regarding connectivity between foraging and roosting polygons. Connectivity was highest to nearby foraging polygons and declined non-linearly as distances increased (Fig. 3). The distance between roosting sites may not necessarily reflect the similarity in foraging sites used, since nearby roosts may exhibit quite different connectivity patterns (Fig. 2 and Appendix B: Figures B1-B3).

Oystercatcher counts and prey biomass
The finer scale analysis of roost counts within a tidal basin indicated that available biomass of benthic stocks only explained significant variation in roosts counts at Vlieland (Table 1). Benthic stocks explained little variation in roost counts at Balgzand and Schiermonnikoog (the analysis was not performed for Ameland as there is only one roost counting area on Ameland), and variables in the top model (based on BIC) were not significant, although p-values were close to 0.05 (Table 1). Whether roost counts were compared to benthic surveys before winter (t) or after winter (t +1) had a negligible effect on explained variance (Table 1). For Vlieland, for both t and t + 1, roost counts increased with decreasing levels of H. diversicolor and increasing M. edulis. The top model for t included Year as a categorical variable, suggesting large year-to-year differences. Year was not included for t + 1, and instead C. edule biomass surveyed after winter (t + 1) were important in explaining variation in roost counts, whereby roost counts increased with increasing C. edule (Table 1; Fig. 4d). In contrast, C. edule explained minimal to no variation in roost counts at Balgzand (Fig. 4b) and Schiermonnikoog (Fig. 4d).
For the analysis at the scale of tidal basins, the top model revealed a significant positive relationship between the normalized roost counts  Fig. 1.   Fig. 3. Travel distances from roosts to foraging areas in each tidal basin (A) and the relation between the connectivity and the distance (m) between roosts and foraging areas (B). Black stripes are the median connectivity on a distance interval (m) indicated by the stripe width. Connectivity is on a log scale in this graph. and C. edule biomass (Table 1; Fig. 4a). Benthic surveys performed after winter explained almost twice as much variation in roost counts as surveys performed before winter (Table 1). This result was independent of any particular tidal basin, since rerunning the analyses but excluding a tidal basin yielded slightly lower R 2 when excluding Balgzand (0.41), and higher R 2 when excluding Schiermonnikoog (0.65), Ameland (0.52) and Vlieland (0.65). When using benthic surveys performed before winter (t), our approach of only estimating biomass in foraging areas connected to roosting sites performed similarly to a buffer zone of 7.47 km. In both analyses C. edule was the most important benthic species. However, when using benthic surveys performed after winter (t +1), our approach explained more variation in roost counts than a simple buffer approach (Tidal basin = 0.49, Buffer = 0.19; Table 1).

Discussion
In this study, significant relationships were identified between shorebird counts and benthic prey stocks, whereby the effect of cockles on oystercatcher counts was the most pronounced. Relationships were stronger at larger spatial scales of tidal basins than at the smaller scale of roost counting areas. Furthermore, benthic prey stocks explained more variation in shorebird counts when only considering the food utilised by the birds at the roosts, i.e. by estimating connectivity between roosting and foraging areas rather than food availability in a general buffer zone around the roost. This may explain why some studies have had difficulty in relating benthic prey stocks to shorebird population size variation (e. g. Verhulst et al., 2004), and supports the approach of identifying the connectivity between roosting and foraging areas from GPS tracking data for oystercatchers. The relationship between shorebird counts and benthic prey stocks at the scale of tidal basins was significantly stronger when these were related to benthic surveys performed after winter (t + 1) rather than before (t). Our study also revealed intriguing foraging patterns with areas that were visited at different intensities and others that were not visited, with the distance to the roost as an important explaining factor. Furthermore, foraging patterns varied amongst tidal basins in terms of the distance travelled between roosting and foraging sites. We identified areas that may be important for foraging oystercatchers, and our method of establishing connectivity between roosting and foraging areas may provide a tool for prioritizing conservation of intertidal areas used by shorebirds.

Trends of roost counts and available biomass
The availability of benthic biomass explained significant variation in roost counts at the scale of tidal basins, but much less so at finer spatial scales of the roost counting areas. Oystercatchers may thus be drawn to a tidal basin given available food stocks but do not appear to distribute evenly amongst available roosting sites as per our expectation. How they distribute amongst available high tide roosts may be affected by other factors, such as the distance from foraging sites, shelter from weather, predation risk or disturbance (Rogers et al., 2006;Rosa et al., 2006;van Gils et al., 2006;Folmer et al., 2010;Horn et al., 2020). Many of these factors vary among days, contrary to the benthic food stock, and can Table 1 Results of model selection for linear models of how availability of benthic biomass explain the number of roosting oystercatchers at the scale of roost counting areas (Balgzand, Schiermonnikoog, Vlieland), and the scale of the Wadden Sea (Tidal Basin). Analyses are shown for benthic surveys performed before winter (t) and after (t + 1). Only the top performing models in terms of BIC are shown. A + orindicates the sign of the coefficient. The top five models for each scale are shown in Appendix B: Tables B1 -B4.  thus influence numbers on a given day (Rappoldt et al., 1985). In contrast, we only detected limited exchanges of birds travelling from one tidal basin to another, indicating that individuals generally stay within a basin for the entire winter. Given the study setup, it was possible that GPS data would not be retrieved if a bird permanently emigrated. However, a study of colour-ring data supports our findings that during winter there was very little movement between the eastern and western parts of the Dutch Wadden Sea (Allen et al., 2019). The contrasting dynamics of these fine-and large-scale processes may thus explain why benthic biomass may better explain large-scale patterns of shorebird populations but not those at finer spatial scales. Given the limited exchange of individuals among tidal basins, it raises the question regarding the mechanism driving the inter-annual variation in the number of oystercatchers wintering in each tidal basin. Local population dynamics will be important, but the rates of change in population size appear to be too high for the trends to be driven by annual variation in survival and reproduction. The Dutch Wadden Sea oystercatcher population increases drastically from summer to winter with three to four times as many oystercatchers in winter compared with summer (Hornman et al., 2020). The trends in the winter distribution of oystercatchers may therefore be driven by the site choices of wintering migratory individuals, and given our results, the site choice appears to be influenced by available food stocks. In addition, understanding the age distribution of local populations may also provide important insights, for instance juveniles are known to be less site-faithful Cayford, 1996, 2014) and their increased mobility may partly drive the observed population trends in the different parts of the Wadden Sea.
Especially the yearly variation in cockle biomass explained variation in oystercatcher abundance among tidal basins. Cockles also represented a large proportion of the available biomass, and since they do not burrow deeply, they are always within reach of foraging oystercatchers (Zwarts and Wanink, 1993;Zwarts et al., 1996). In contrast, mussels, also considered as an important food source of oystercatchers, were not important for explaining variation in oystercatcher abundance among tidal basins. Mussels made up a much smaller part of the benthic biomass and the variation in the biomass of mussels was also less compared with the variation in cockle biomass. Furthermore, mussel beds are increasingly mixed with Pacific oysters (van der Meer et al., 2019) and these mixed musseloyster beds are less attractive to oystercatchers than the pure mussel beds that existed prior to 2002 (Waser et al., 2016).
Our results show how GPS movement data can improve our knowledge of a species' foraging dynamics. Our approach of linking available benthic biomass from connected (i.e. used) foraging polygons outperformed models that simply used a buffer zone around a roost and improves our ability to link bird counts with benthic biomass. Buffer zones likely include areas not used by oystercatchers whilst knowing where oystercatchers forage, and from which roost they travel, encapsulates various processes that a buffer zone excludes, for example how disturbances, low average benthic biomass and other unknown factors may influence foraging distributions. The size of the buffer zone would also need to vary per tidal basin since the GPS tracking data revealed how oystercatchers commuted on average further in Balgzand compared to Schiermonnikoog and Ameland (Fig. 3a). These differences could be explained by general food availability or the differing tidal dynamics between the islands and mainland locations. That oystercatchers commute different distances in different regions is important to know for conservation planning, for example when prioritizing protection of high tide roost and foraging areas. Finally, we identified how the intensity of use of foraging areas varied substantially, with a general pattern that foraging areas that were closer to roosts were used more intensively than those further away (Fig. 3b). A driver of this pattern is likely to be exposure time, given the tidal dynamics of the mudflats. Oystercatchers may also balance travelling and foraging time and avoid commuting to distant foraging areas when nearby flats are exposed, thereby limiting flight costs. These results could potentially be used to simplify our method in the event that GPS movement data is not available, for example defining a buffer zone where available benthic biomass is weighted based on the distance to high tide roosts.

Food stocks and shorebird counts
Benthic surveys are a composite of data sampled in spring (April/ May) and early summer (June/July), while our analyses are based on winter (September-February) roost counts. At the scale of tidal basins, an important finding was that benthic surveys performed after, rather than before winter, explained more variation in roost counts. It is likely that factors that influence benthic fauna before and during winter (i.e. after benthic surveys are performed) are important considerations when relating shorebird counts to available food stocks. For instance, cockle spatfall occurs in the summer and the growth of cockles primarily occurs during the period of the benthic surveys with peak growth in early summer (House and Farrow, 1968). Birds may also contribute to the depletion of the food sources (e.g. Bijleveld et al., 2015;Dokter et al., 2017), and harsh winters may cause mortality amongst the prey species (Burdon et al., 2014;Beukema and Dekker, 2020b). These patterns may explain why cockle biomass at t + 1 correlated most strongly with winter shorebirds counts.
Our results therefore have various implications when relating shorebird counts to benthic food stocks. The first is that benthic surveys at time t do not appear to explain as much variation in shorebird counts, and thus do not appear to be a suitable predictor for upcoming winter shorebird population sizes. Secondly, to understand interannual changes in oystercatcher population sizes at the scale of tidal basins, benthic surveys performed in the summer after the shorebird counts had the most explanatory power. Understanding the processes occurring between the two benthic sampling periods (i.e. our periods of t and t + 1) may further clarify the mechanisms driving oystercatcher winter abundances. Given these intervening processes, it also stands to reason that shorebird counts combined with our connectivity and home range analyses may provide additional information about over-winter food stocks at the scale of tidal basins, and thus serve as an indicator of environmental health for the winter period when benthic surveys are not performed. However, our analysis did not reveal strong correlations between the number of roosting birds and benthic stocks at local spatial scales, and further research is needed regarding how site-specific conditions determine roost choice and how these interact with available food resources.
We were able to explain variation in the size of wintering shorebird populations in tidal basins based on available food stocks, in particular available biomass in foraging areas that are connected to the high tide roost counting areas. Our results indicate how wintering shorebird populations may fluctuate among tidal basins interannually depending on available food stocks, and is an important consideration when analysing local population trends. Although available food stocks may explain variation in the distribution of oystercatchers among tidal basins, the most important prey items (cockles and mussels) have fluctuated without a clear trend in recent years (Beukema and Dekker, 2020b), meanwhile oystercatchers have continued to decline (van Roomen et al., 2012;van de Pol et al., 2014). An approach considering the full annual cycle would, therefore, be needed to determine the relative contribution of winter food stocks to the ongoing decline of the oystercatcher, along with other impacts that vary in space and time. The approach that we present may nonetheless benefit conservation planning, for example, the connectivity analyses enable predictions of how shorebirds like the oysteratcher may be impacted by changes to foraging or roosting sites, such as if foraging areas become affected by additional pressures from shellfisheries or if roost sites are lost or disturbed. Alternatively, the approach we describe could help identify locations where there are no or too few roosting sites based on available food stocks, meaning that new roosts could be created, or where actions could be prioritised to minimise disturbance. These conservation implications highlight the benefit of quantifying the connectivity between roosting and foraging areas in relation to the food landscape and emphasises how having a continuously tracked population of "sentinel" birds can inform about the birds themselves and their food distributions. These potential contributions align with current programs from the program 'Towards a rich Wadden Sea' (Programma naar een Rijke Waddenzee) and birdlife 'Rust voor vogels'.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.