Natural grazing by horses and cattle promotes bird diversity in a restored European alluvial grassland

Study site

We carried out the study in a 32-ha test-enclosure that is part of an ecosystem restoration project of the nature reserve Petite Camargue Alsacienne, located in Alsace, France, north of the city of Basel, Switzerland. The area is situated on the Rhine Island (see the map of the study area and further details in Lovász, Korner-Nievergelt & Amrhein (2021)), where approximately 100 ha of former agricultural fields had been transformed into an alluvial environment. The ecosystem restoration process started in 2014 and involved establishing a heterogenous grassland scattered with shrubs (hawthorn, Crataegus sp; dog rose, Rosa canina), young and old trees (willow, Salix sp; oak, Quercus sp) and bare-ground gravel sites, surrounded by patches of old forests (oak; ash, Fraxinus sp). Water from the Rhine river is diverted through the area as a small stream, and ground-water ponds have also been created. The aim has been to maintain the heterogeneity of the area and to promote a species-rich alluvial ecosystem.

The study area has been managed using the approach of natural grazing, a low-intensity year-round mixed grazing regime (<0.5 animals per ha), with the aim of substituting extinct wild herbivores, such as the wild horse Equus ferus or the aurochs Bos primigenius, with domestic breeds kept in semi-wild conditions (Lovász, Korner-Nievergelt & Amrhein, 2021). This natural grazing approach includes minimal human intervention and no systematic winter feeding (Linnell, Kaczensky & Lescureux, 2016; Vermeulen, 2015).

Land-cover characteristics

We carried out surveys to characterize the landscape on the grazed study site and the non-grazed control area. We used satellite imagery (Google Maps) in QGIS (version 3.4.4-Madeira) to divide the study area into 50 × 50 m grid cells. With the help of a printed map of these grid cells (using the QGIS printout layer) and a handheld GPS (application: Gaia GPS), the observer (L.L.) visited each grid cell. The observer visually estimated the percentage of types of land cover inside the cell by marking the corner points with flags and standing in the midpoint of the cell (which was already indicated with permanent markings used in a different study). When it was necessary for adequate estimation of the cover types (for example, when visibility decreased due to trees), the observer walked along one diagonal of the grid cell to determine the cover attributes. Additionally, the observer compared the observations to the QGIS map image of the area while estimating the percentages. The comprehensive land-cover survey took place once during the study period in summer 2020. The observed large-scale habitat characteristics did not noticeably change during the 1.5-year study period, as indicated by weekly visual inspections.

We followed Cunningham & Johnson (2019) when defining a small number of categories for land-cover types in our study area: tree-cover (all trees of ca. ≥ 3 years old), sapling-cover (all growth stages of young trees of ca. <3 years old), shrub-cover, cover of meadow, surface of bare ground, and surface of water on the respective grid cell. Based on Hildén (1965), these characteristics were suggested as important determinants for habitat selection of birds in our study site.

We identified six different land-cover types, indicating a heterogenous mosaic-like environment. The cover by meadow (grassland) was the most abundant (61%). Bare ground, tree cover and shrub cover represented only a small proportion of all cover types (6.9%, 3.4%, and 1.7%, respectively), as the area was restored with the aim of recreating an open alluvial grassland. Single trees and shrubs were patchily distributed without providing extensive cover. The relatively high percentage of sapling cover (17.9%) consisted mainly of poplar (Populus sp) trees. Such sapling patches were often relatively large (up to more than 1 ha), with either bare ground or species-poor meadow patches as understory. Extensive areas with bare ground were only found near water bodies. Water bodies represented 9.2% of all cover types, and all of them were groundwater ponds with relatively stable water levels throughout the year.

Cattle and horse data

To counteract the natural succession by willow and poplar saplings encroaching upon the meadows and to foster a heterogenous environment, large grazers were progressively introduced into the test enclosure. The introduction began with horses in September 2018, followed by cattle in January 2019. The initial number of five Highland cattle and five Konik horses increased to seven cattle (five adults and additionally two calves born on-site) in the summer of 2019, and to seven horses (five adults and additionally two young mares brought from other areas) in the summer of 2020. Over the duration of our data collection, the overall stocking rate increased from ∼0.3 animals per ha to ∼0.4 animals per ha.

To determine the density distribution of the horses and cattle across the study area, we recorded hourly GPS positions of the grazers and considered these as a proxy for the space use of the animals. All grazers (14 in total) were equipped with GPS collars (Followit, type Pellego) upon their arrival to the study area by the management of the nature reserve. The collars registered the position of the animals once per hour. The first round of data collection was from January to July 2019 (the GPS data from this first data collection were used in Lovász, Korner-Nievergelt & Amrhein, 2021). The GPS collars recorded for 5 to 6 months with one charge, and it was only possible to subsequently recharge the devices of the horses; the cattle were not tame enough to handle the collars without the need of capturing and blocking the animal (which was only possible in winter during the regular roundup of cattle). Thus, more durable GPS collars (Followit, type Tellus) were used for two cattle starting from March 2020, instead of the five formerly used collars (the structure of the collected data of the two collars was identical). Data were then continuously collected until March 2021. In total, data collection spanned one year and six months. In our analysis, we accounted for different numbers of collars during the different time periods.

The data were downloaded through satellite processing from the interface of the GPS collar provider (Followit, Lindesberg, Sweden), and, thus, no contact with the animals was necessary to access the data. GPS collars comply with animal welfare requirements for horses (Collins et al., 2014; Hennig, Beck & Scasta, 2018) and have been widely used also on cattle for decades, without causing harm or disturbance (Turner et al., 2000; Ungar et al., 2005).

The GPS positions in our dataset may show some imprecision because GPS accuracy can vary depending on atmospheric conditions or due to satellite or receiver errors (Hurn, 1993), caused by satellite geometry (Dussault et al., 2001), topography, overhead canopies, or adjacent structures (Di Orio, Callas & Schaefer, 2003; Moen et al., 1996). The GPS collars did not record HDOP (horizontal dilution of precision), so we could not statistically account for inaccuracy (Langley, 1999). However, as the rate of GPS positions falling outside of the fenced area were negligible (see Lovász, Korner-Nievergelt & Amrhein, 2021), we assumed that imprecision would not strongly influence our results.

Bird abundance and species richness data

Permissions

The fieldwork was carried out with the permission of the national nature reserve Petite Camargue Alsacienne.

Statistical analysis

For the analysis, we used the 50 × 50 m grid cells of the Universal Transverse Mercator (UTM) coordinate system’s grid over the study area. We first calculated the mean numbers of species and of individuals in the main study site and the control site during the breeding season. Subsequently, to assess how bird abundance and species richness were related to grazer density, we used only data from our main study site (since the control site was not grazed).

To analyze how the abundance of birds of different guilds were related to grazer density, season, and habitat characteristics, we used generalized additive mixed models (GAMM) that we fitted to the data of each bird guild separately. Our response variable was the number of bird individuals counted per survey per 50 × 50 m grid cell, summed over all species belonging to a bird guild. First, we fitted a negative binomial model, but a posterior predictive check showed that the data contained too many zeros. We therefore assumed a zero-inflated negative binomial distribution of the bird counts with a constant proportion of zero values. We used the logarithm as link function for the count model. As predictors, we used grazer density, date (day of year), and habitat characteristics.

Grazer density was measured as the sum of the hourly horse and cattle GPS positions per grid cell over the 3 weeks prior to a bird survey. Our previous study (Lovász, Korner-Nievergelt & Amrhein, 2021) showed that the effects of horse vs cattle densities on bird count densities were rather similar, and, therefore, we pooled the data for horse and cattle GPS positions. We further assumed that grazer space-use patterns earlier than 3 weeks before the respective bird survey did not substantially influence space use of birds (based on Lovász, Korner-Nievergelt & Amrhein, 2021; this was an arbitrary cut-off time we chose for pooling grazer positions). Because we assumed that the detectability of birds was relatively homogeneous across the main site and control site areas, and our aim was not to estimate the total number of birds present, we did not take detectability into account in our analyses, following Buckland et al. (2001). As habitat characteristics, we used the land-cover variables “bare ground”, “sapling”, “shrub” and “tree”. To avoid redundancy in the predictor variables due to cover types summing up to 100%, we did not include the most abundant land cover type “meadow”. The cover type “water” was included as a binary variable indicating presence or absence of water bodies on a 50x50 m cell. All predictors were z-transformed so that their means were zero and their standard deviations were one.

We used a two-dimensional tensor product smooth for grazer density and day of year. To obtain the smooth terms, we used three cubic regression splines, as implemented in the function t2 of the package mgcv (Wood, Scheipl & Faraway, 2013). Because some grid cells at the edges of the study area had sizes smaller than 50 × 50 m, we used the logarithm of the grid cell size as an offset in the linear predictor of the counts.

The grid cell ID was first used as a random factor to account for repeated measures. However, Markov chains of these model fits did not converge well. We therefore omitted the random factor from the model. To assess the strength of the pseudoreplication that may have arisen by omitting the grid cell ID as a random factor, we estimated the proportion of residual variance that can be explained by among-grid cell variance using normal linear mixed models. This proportion was relatively low for all guilds (wetland birds: 8.1%, woodland birds: 2.6%, open-area birds: 0.7%, aerial birds: <0.1%). We thus considered that pseudoreplication due to repeated measures of grid cells was low in our study. However, we acknowledge in the interpretation of our results that our compatibility intervals (Amrhein & Greenland, 2022) may underestimate statistical uncertainty, particularly for wetland birds.

To analyze how species richness was related to grazer density, day of year, and habitat characteristics, we used the number of species detected during one survey on each of the 50 × 50 m grid cells as the outcome variable in another GAMM. Similarly to the model for the abundance, here we also used a zero-inflated negative binomial model. We used the same predictors as those used for the analyses of abundance, except that the logarithm of the grid-cell size was used as a covariate instead of an offset. For analyzing species richness, we used the grid cell ID as a random factor, to account for repeated measures in the grid cells. Unlike for the abundance models, Markov chains for the species richness model including cell ID as random factor converged well.

All models were fitted using Hamiltonian Monte Carlo methods as implemented in Stan (mc-stan.org; Carpenter et al., 2017). We used the R-package brms as an interface to Stan (Bürkner, 2017). We used the default prior distributions that were optimized for z-transformed predictors and that did not markedly affect the results: the positive range of t (3, −2.3, 2.5) for the intercept of the negative-binomial model, a logistic (0,1) for the intercept of the zero-model, flat priors for all model coefficients, positive range of t (3,0,2.5) for the variance parameters (random effects and coefficients of the smoother), and gamma (0.01, 0.01) for the shape parameter of the negative-binomial model.

After fitting the models, we checked convergence graphically by looking at the Markov chains and using the R-hat value. We further inspected the goodness-of-fit via posterior predictive model checking using the package shinystan (Gabry, 2017). We compared density distributions of our data with density distributions of data replicated from the models, and we compared means, standard deviation, minimum and maximum of our data with replicated data based on the models. We also checked for spatial correlation by looking at the semi-variance and by visualizing the residuals on a map, and we did not detect conspicuous spatial autocorrelation.

Posterior predictive model checking revealed that a Poisson distribution or negative binomial distribution did not describe the data distribution well, but zero-inflated negative binomial distributions fitted well to all our response variables.

Bird counts on the main study site and the control site

In total, we made 5,447 observations of individuals of 87 bird species during the entire study (see Table S2 for the species list). On the main (grazed) study site, we surveyed 112 grid cells during a total of 73 surveys between 31 January 2019 and 22 March 2021. During these 73 surveys, we counted 5,104 individual birds.

For comparison, we conducted surveys at the (non-grazed) control site during a period of about two months in the breeding season, from the 89th to the 142nd day of the year. On this control site, we surveyed 121 grid cells during seven surveys and compared the number of individuals observed on the control site and the main site during the same period in 2019 (six surveys) and in 2020 (seven surveys). The mean (±SE) number of observed bird individuals per survey per grid cell was 0.93 (±0.14) on the main site, and 0.43 (±0.08) on the control site. Also, the mean species richness (number of observed bird species per survey per grid cell) was higher in the grazed area (0.27 ±0.01) than on the control site (0.18 ± 0.01) during the period of comparison.

Effects of grazer density, season, and habitat on bird abundance

In the following, we describe how the numbers of observations of aerial, open-area, wetland, and woodland bird individuals changed seasonally, in relation to grazer density (i.e., the numbers of GPS positions of grazers per grid cell summed over a 3-week period) and in relation to the habitat characteristics for our main study site.

Aerial-foraging birds

In total, we made 831 observations of aerial-foraging bird individuals. Almost all birds of this guild were swallows and swifts that do not breed at the study site (422 Barn swallows, 138 House martins, 269 Common swifts). During spring migration (April and early May), the numbers of individuals in the aerial foraging guild reached a peak, particularly in areas with low grazer density (Figs. 1, 2A and 2B). During the rest of the year, there was no obvious relationship between the numbers of aerial-foraging birds and grazer densities (Figs. 1, 2B; the apparent increase at high grazer densities on 1 September in Fig. 2A had high statistical uncertainty, as indicated by the wide compatibility interval). Aerial birds were observed in similar numbers (average 0.2 individuals per grid cell per survey: see the heatmap contour lines in Fig. 1) throughout the whole range of grazer density in the breeding season (April–July).

There were slightly more counts of aerial birds on plots with higher proportions of meadow cover (Fig. S2A). We expected to find more aerial-foraging birds over water bodies, since swallows often forage on insects above water. However, we did not find such a difference when comparing grid cells with and without water bodies. Also, the percentage of coverage in all other habitat categories showed almost no relationship with bird abundance in the aerial-foraging guild (Fig. S2A).

Effects of grazer density, season, and habitat on the species richness of birds

We found species richness to be relatively constant throughout the year, with an average of about 0.3 species counted per grid cell per survey (Fig. 3).

Independently of the land-cover variables included in the model, most species per grid cell were observed at average grazer densities (24.1 GPS positions per grid cell per 3 weeks) on the median day of the study period (1st June), and this pattern was similar for all seasons (Fig. 4 and see Fig. S3). However, the pattern at higher grazer densities was unclear due to the large compatibility intervals, at least for the median day of the study period (Fig. S3, right side).

The number of bird species decreased slightly with increasing cover of bare ground, meadow, and saplings (Figs. 5A/5C), while it increased with more extensive shrub and tree cover (Figs. 5D/5E) and was also higher in grid cells with water bodies (Fig. 5F).