1. Maintaining the biodiversity of agricultural ecosystems has become a global imperative. Across Europe, species that occupy agricultural grasslands, such as Black-tailed Godwits (Limosa limosa limosa), have undergone steep population declines. In this context, there is a significant need to both determine the root causes of these declines and identify actions that will promote biodiversity while supporting the livelihoods of farmers.
2. Food availability, and specifically earthworm abundance (Lumbricidae), during the pre-breeding period has often been suggested as a potential driver of godwit population declines. Previous studies have recommended increasing the application of nitrogen to agricultural grasslands to enhance earthworm populations and aid agricultural production. Here we test whether food availability during the pre-breeding period affects when and where godwits breed.
3. Using large-scale surveys of food availability, a long-term mark-recapture study, focal observations of foraging female godwits, and tracking devices that monitored godwit movements, we found little evidence of a relationship between earthworm abundance and the timing of godwit reproductive efforts or the density of breeding godwits. Furthermore, we found that the soils of intensively managed agricultural grasslands may frequently be too dry for godwits to forage for those earthworms that are present.
4. The increased application of nitrogen to agricultural grasslands will therefore likely have no positive effects on godwit populations. Instead, management efforts should focus on increasing the botanical diversity of agricultural grasslands, facilitating conditions that prevent hardening soils, and reducing the populations of generalist predators.

Study Area
We studied the relationship between the reproductive decisions of godwits and earthworm abundance in the Haanmeer (52.9222°, 5.4353°), a 200-ha polder in southwest Friesland, The Netherlands. The Haanmeer polder is one of the last remaining areas with high densities of breeding godwits (~0.5 pairs/ha) in southwest Friesland — which represents the core of the godwit breeding distribution (Kentie et al., 2016) — and is part of a long-term landscape-scale study of godwit demography and breeding biology (see Groen et al., 2012; Senner et al., 2015 for more details). The Haanmeer consists of fields that are maintained under two different management schemes: (1) high-intensity grasslands (70 ha) and (2) low-intensity grasslands (130 ha). The two management types differ in a number of ways, including water levels, botanical richness, the number of cattle grazed, amount and type of nitrogen application, and mowing regime (Table 1).

In particular, the application of cow manure — the main source of nitrogen input — averages ~ 70 tons/ha on the high-intensity grasslands, but only ranges from 0 to 25 tons/ha on the low-intensity grasslands (A. Stokman, S. Venema, and D. Postma pers. comm.). In addition, the water level in the ditches surrounding the fields is kept consistently lower in the high-intensity grasslands (unpubl. data Wetterskip Fryslân) resulting in a ~70 cm difference in spring groundwater level (unpubl. data MAV & NRS). These two management types and the differences between them are typical for our study area (see Groen et al., 2012) and the godwit breeding range in the Netherlands (Teunissen et al., 2012).

Data collection on godwit behaviour
During 2008–2012, we captured 87 adult and 392 juvenile godwits in the Haanmeer and surrounding polders. Upon capture, adults and chicks were individually marked, bled for molecular sexing, and their biometrics were measured (see Loonstra et al., 2019). The latter included the measurement of their exposed culmen to the nearest millimetre (hereafter ‘bill length’). The individual markings consisted of a unique colour-ring combination for adults and appropriately sized chicks, and a single flag with a unique alphanumeric code for smaller chicks.

            In 2013 and 2014, we searched daily for these individually-marked godwits as they returned to the Haanmeer. We started with the first arrival of adults and continued until the first egg was found (2013: 8 Mar–10 Apr; 2014: 8 Mar–9 Apr). When an individually-marked godwit was encountered, its colour code, location, and behaviour — foraging, preening, displaying, etc. — were recorded. During the same period, we made daily foraging observations on a subset of individually-marked females. For these females, we recorded the number of probes and successes during consecutive 3-min periods, where a probe was defined as a downward movement of the bill into the soil and a success as a swallowing motion (Senner and Coddington 2011). When possible, we also recorded the prey type consumed. In 2013, we observed a foraging individual for up to 10 3-min periods for a total daily foraging observation length of 30 min. In 2014, we did this for up to three, 3-min periods for a total daily foraging observation length of 9 min.

            In 2014, we also tracked the locations of three females throughout their pre-breeding period using 7.5 g solar-powered UvA-Bits GPS trackers (Bouten et al., 2013). We attached these trackers in 2013 with a leg-loop harness using 2 mm nylon rope (see Senner, Stager, Verhoeven et al., 2018 for more details). These trackers stored location estimates once every 5 min when the battery was fully charged and once every 15–30 min in all other instances. One of these females was still carrying its transmitter in 2015 and this provided us with an additional year of spatial distribution data during the pre-breeding period.

            After the first egg was found, we stopped our foraging observations and started intensively searching for godwit nests. We floated and measured the dimensions of eggs from nests with completed clutches to estimate their incubation stage and egg volume (Liebezeit et al., 2007, Schroeder et al., 2009). We monitored nests every 2–3 days to determine their precise hatching date. We made a particular effort to associate individually-marked birds with nests. For nests found in the laying phase, we determined their lay date based on the assumption that godwits lay one egg per day (Cramp & Simmons 1983). For nests that hatched but were not found in the laying phase, we back-calculated the lay date by subtracting 26 days (the combined average laying and incubation period; Verhoeven et al., 2020). In cases where we found a nest with a complete clutch, but did not know the hatching date due to predation or abandonment, we estimated the lay date by subtracting the incubation stage (as derived from the flotation method) from the date the nest was found.

Data collection on food availability
On four separate occasions, we sampled the entire Haanmeer for below-ground invertebrates along a 100 × 150 m grid (Fig. 1) as local soil invertebrate populations can vary considerably (Timmerman et al., 2006). In 2013, we collected 130 samples from 1–3 March — before the arrival of the first godwits — and 136 samples from 4–6 May, just after the last of our individually-marked females had begun incubation. We sampled six additional fields in May that had not been sampled in March. In 2014, we collected 135 samples from 8–10 March and 136 samples from 1–4 May. We inadvertently skipped one sample at the end of a transect in March 2014.

            We collected invertebrates in 20 × 20 × 20 cm soil samples. To estimate the depth at which invertebrates were present, we split each sample into four pieces along its horizontal axis immediately after removal from the ground, resulting in four equally sized slices 5 cm in height. Each slice was sealed separately in a plastic bag and stored between −5°C and +5°C. We processed the samples within two weeks by hand-sorting each slice, which is considered the most reliable method (Nordström and Rundgren 1972; Edwards and Lofty 1977). We then cleaned all invertebrates, weighed them to the nearest 0.01 g, and measured their (relaxed) length to the nearest 0.01 cm.

            Soil penetrability affects earthworm availability for godwits by influencing the ease with which godwits can probe the soil for invertebrates (Kleijn et al., 2011). We thus also measured the penetration resistance of the soil in the same fields in which we sampled soil invertebrates with a penetrometer with a 1-cm² cone (Eijkelkamp, Giesbeek, penetrometer 06.01.14). From late March to mid-May in 2013 we measured the penetration resistance along four 60-m long transects daily. Two transects each were located in high-intensity and low-intensity grasslands. We sampled every 4 m along each transect, resulting in 15 measurements per transect. To ensure that the results from our original two transects adequately captured variation in penetrability across the study area, we added eight additional transects in low-intensity grasslands in 2014 for a total of ten transects in low-intensity grasslands and two in high-intensity grasslands. For these transects, we started in the middle of March and sampled every 2 m along 10-m long transects (six measurements per transect).

Statistical analysis of annual and seasonal differences in invertebrate biomass
The recorded invertebrate biomass was zero-inflated due to the number of samples without earthworms (Table 2) and included only positive values in those cases when invertebrates were recorded. As a result, we used a two-part binomial-gamma hurdle model. The initial binomial model predicted the probability of recording earthworms in a sample, using a generalized linear mixed model (GLMM) with binomial error structure and a logit link function. The gamma model then predicted the invertebrate biomass when invertebrates were recorded using a GLMM with a gamma error structure and log link function. Both models had year, month and grassland management intensity as categorical predictor variables, and sampling location as a random intercept to account for pseudo‐replication. We obtained chi-squared values for the significance of the predictor variables from likelihood ratio tests of nested models with and without the variable of interest. In the final model, we removed grassland management intensity as a predictor variable because it did not significantly improve the fit in either the binomial or gamma model (see Results). By multiplying the predictions of the two parts of the hurdle model (i.e., weighting predicted biomass by the probability of observing an invertebrate), we obtained predictions of invertebrate biomass at both sampling time points in each year.
            To additionally test for potentially confounding variables that might influence underlying worm abundances and, therefore, godwit breeding biology, we also ran the same binomial and gamma models with soil type — (1) more clay, (2) more sand, or (3) more peat (Fig. 1) — replacing grassland management intensity as a categorical predictor. We excluded three earthworm samples for this analysis, either because the sampling area had previously been filled with an unknown soil type (n = 2) or because the sample was the sole representative of a soil type category (more loam; n = 1; Fig. 1). Soil type, however, was not a relevant predictor variable of worm abundance in either the binomial (χ² = 3.07, df = 2, p = 0.215) or gamma models (χ² = 1.76, df = 2, p = 0.414), and thus was not included in our final models (Table S1).

Relating godwit behaviour to earthworm biomass
Because grassland management intensity and soil type were not related to earthworm biomass (see Results), and because prevailing management directives (Brandsma 1999; van der Weijden and Guldemond 2006; Kleijn et al., 2009a) posit a direct link between earthworm abundance and godwit behaviour and abundance — proposing that earthworms explain observed variation in godwit reproductive phenology and habitat selection — we focused our analyses on testing the relationship between earthworm biomass and godwit behaviour and abundance.

To test for an effect of earthworm abundance on godwit foraging rates, we used only the biomass of earthworms in our analyses because earthworms comprised ~90% of prey mass found in March (see Results). Furthermore, we used only the earthworm abundance in the top 10 cm because those earthworms found deeper in the soil are not within reach of female godwits (see Results). We plotted the prey intake rate of female godwits during the pre-breeding periods of 2013 and 2014 against the average earthworm biomass in the top 10 cm in March of the meadow in which each godwit was foraging (Fig. 2). We used only females that probed more than 220 times (25^th percentile) and less than 370 times (75^th percentile) to control for differences in motivation between individuals, as this is known to cause differences in observed intake rates (Duijns et al., 2015). We then used Holling’s type II functional response equation to explore the relationship between intake rate and earthworm density (Holling 1959). We had no empirical data for handling time (Th) or instantaneous area of discovery (a), and used the least squares method to calculate the values for both these parameters that yielded the best fit to our data.

            We next tested for an effect of earthworm abundance on where godwits foraged during the pre-breeding period. In a generalized linear mixed model with a Poisson error structure and log link function, we related the number of unique individuals seen on a meadow during the pre-breeding period to the average earthworm biomass in the top 10 cm of that meadow in March, while accounting for the size of the meadow (continuous covariate), the year (two-level factor), and the fact that samples from the same meadow are not independent (random intercept). We ran the same model to relate the number of nests found in a meadow to the average earthworm biomass in the top 10 cm of that meadow in March. For the latter analysis, we included only nests laid before 1 May in order to exclude replacement clutches as much as possible (see Verhoeven et al., 2020).

            Third, we related earthworm abundance to godwit breeding phenology. We used linear mixed models to relate the (1) arrival dates, (2) lay dates, and (3) duration of the arrival–laying interval of females to the average earthworm biomass in March in the top 10 cm within 200 m of their nests. We did the same, separately, for the biomass within 300 m of their nests. We included year as a two-level factor and individual as a random intercept. We used the earthworm biomass from within these distances because 70% of the foraging locations of the GPS-tagged females were within 200 m and 88% within 300 m (Fig. S1a), which is similar to the spatial distribution exhibited by the rest of the population (Fig. S1b). In these analyses we again included only nests laid before 1 May and used chi-squared values to assess the significance of each predictor variable.

Finally, we assessed whether earthworms influenced godwit reproductive investment. For this, we used linear models to relate the lay dates and average egg volume of all nests (including those of unknown females) laid before 1 May to the average earthworm biomass in the top 10 cm within 200 m and 300 m of these nests in March. We calculated egg volume using the formula: length × breadth² × 0.52 (Schroeder et al., 2009). Female size and other covariates were not included in these models, as we have previously shown that they have little effect on egg volume (Verhoeven et al., 2019).

Additionally, we tested the robustness of our worm sampling scheme and our choice to average worm measures across samples within 200 m and 300 m of a female’s nest or within the meadow in which they nested (i.e., within 500 m or less). Using both Moran’s I and Mantel tests, we found evidence of positive spatial autocorrelation or ‘clustering’ of earthworm biomass at these distances in March 2013 (Moran’s I for 200 m = 0.27, p = <0.001; 300 m = 0.15, p = <0.001; 500 m = 0.1 p = <0.001) and March 2014 (Moran’s I for 200 m = 0.20, p = 0.004; 300 m = 0.09, p = 0.019; 500 m = 0.05. p = 0.021). Our Mantel tests also indicated that sampling locations nearer to each other had more similar worm abundances, i.e., positive spatial autocorrelation (r_March2013 = 0.13, p = <0.001; r_March2014 = 0.11, p = <0.001).

We obtained F-values for the significance of all predictor variables from F-tests of nested models with and without the variable of interest. We ran all models in the R programming environment (R Core Team 2018, version 3.5.1) using the package “lme4” (Bates et al., 2015).

All files are numbered, a new number for each analysis/topic. Related files are given the same number but a different letter, i.e. files 3a, 3b, 3c, 3d are related to each other. 

1a.March2013_soil invertebrate sampling data
This file contains the results from our invertebrate sampling in March 2013. It has a tab containing all the worms and a tab containing all the non-worms. The section refers to the depth of the sample; 1=0-5cm, 2=5-10cm, 3=10-15cm, and 4=14-20cm.

1b.May2013_soil invertebrate sampling data
This file contains the results from our invertebrate sampling in May 2013. It has a tab containing all the worms and a tab containing all the non-worms. The section refers to the depth of the sample; 1=0-5cm, 2=5-10cm, 3=10-15cm, and 4=14-20cm.

1c.March2014_soil invertebrate sampling data
This file contains the results from our invertebrate sampling in March 2014. It has a tab containing all the worms and a tab containing all the non-worms. The section refers to the depth of the sample; 1=0-5cm, 2=5-10cm, 3=10-15cm, and 4=14-20cm.

1d.May2014_soil invertebrate sampling data
This file contains the results from our invertebrate sampling in May 2014. It has a tab containing all the worms and a tab containing all the non-worms. The sections refer to the depth of the sample; 1=0-5cm, 2=5-10cm, 3=10-15cm, and 4=14-20cm.

2.Bill length of Female godwits caught in Haanmeer
This file contains the measured bill lengths of females in the Haanmeer

3a.2013+2014_Worms_top 10cm_with management intensity
This file contains the average worm biomass of worms in the top 10cm in March and May and has a column that identifies the management intensity.

3b2013+2014_.Worms_top 10cm_with soil type
This file contains the average worm biomass of worms in the top 10cm in March and May and has a column that identifies the soil type.

4a.2013 soil penetration data
This file contains the soil penetration resistance data collected in 2013.

4b.2014 soil penetration data
This file contains the soil penetration resistance data collected in 2014.

5.2013+2014_foraging observation data
This file contains the foraging observations made in 2013 and 2014 and the average number of worms in the field of foraging in order to fit intake rates.

6.2013+2014_sightings in field and worm biomass
This file contains the number of sightings in a field in 2013 and 2014, the size of the field and the field’s worm biomass in order to look for a relationship.

7.2013+2014_nesting density and worm biomass
This file contains the number of nests in a field in 2013 and 2014, the size of the field and the field’s worm biomass in order to look for a relationship.

8.2013+2014_Arrival-Laying Date of marked females and worm biomass
This file contains the arrival date, laying date and the arrival-laying interval of females as well as the worm biomass within 200 and 300 meters of the nest.

9a.2013+2014_all laydates and worm biomass within 200m
This file contains all of the known lay dates of godwits in 2013 and 2014 and the worm biomass within 200 meters of those nests in order to look for a relationship.

9b.2013+2014_all laydates and worm biomass within 300m
This file contains all of the known lay dates of godwits in 2013 and 2014 and the worm biomass within 300 meters of those nests in order to look for a relationship

10a.2013+2014_all egg volumes and worm biomass within 200m
This file contains all the average egg volumes of nests from 2013 and 2014 in which all 4 eggs were measured and includes the worm biomass within 200 meters of those nests in order to look for a relationship.

10b.2013+2014_all egg volumes and worm biomass within 300m
This file contains all the average egg volumes of nests from 2013 and 2014 in which all 4 eggs were measured and includes the worm biomass within 300 meters of those nests in order to look for a relationship.

11a.2013+2014+2015_Distance of tracking location from nesting location for FigS1
This file contains the distance from the nest of three females that were tracked via GPS-transmitters in 2014 and 2015.

11b.2013+2014_Locations of observation and field centers for FigS1
This file contains the distance from the nest to the center of the field of observation of 42 individually-marked females that were resighted during the 2014 pre-breeding period.

12a.2013_March_sampling location and worm biomass
This file contains the worm biomass in the top 10cm in March 2013 for each sampling point and the location (latitude + longitude) of that sampling location, used to look for the presence of spatial autocorrelation of worm biomass.

12b.2014_March_sampling location and worm biomass
This file contains the worm biomass in the top 10cm in March for each sampling point and the location (latitude + longitude) of that sampling location, used to look for the presence of spatial autocorrelation of worm biomass.