Impacts of sheep versus cattle livestock systems on birds of Mediterranean grasslands

Study area and parcel selection

The study was conducted in Castro Verde Special Protection Area (SPA), located in southern Portugal (Fig. 1). The climate is Mediterranean, with hot summers (30–35 °C on average in July) and mild winters (averaging 5–8 °C in January), and over 75% of annual rainfall (500–600 mm) concentrated in October–March. The landscape is flat or gently undulating (100–300 m), mainly dominated by open areas used for rainfed pastures (ca. 60%) and annual crops (ca. 25%), and to a less extent by open woodlands (ca. 7%)¹⁵.

Since 1995, part of the study area has benefited from a CAP agri-environment aiming to protect the traditional farming system¹⁶. This scheme provides financial support to farmers for agricultural practices considered favourable to conservation, including the traditional rotation of cereals and fallows, the maintenance of low stocking rates (usually related with sheep grazing systems), and sowing of crops benefiting grassland birds¹⁶. However, in recent years the traditional farming system has been declining, with many farmers converting to specialized livestock systems, mainly, cattle grazing systems, with an increase of stocking rates^7,15.

Parcel selection started by identifying grasslands grazed by either sheep or cattle, based on parcel-level statistical information from 2010 provided by the Portuguese Ministry of Agriculture⁷. To minimize potentially confounding effects of adjacent land uses (edge effects) and other non-crop elements within parcels on bird assemblages, we excluded parcels less than 100 m from shrubland or forested areas, with shrub and tree cover > 5% and with a minimum size of 10 ha. In January 2019 we visited 100 pre-selected parcels which were grazed by either sheep or cattle in 2010 in order to confirm the parcel land use in the agricultural year of 2018/2019, aiming to sample a balanced proportion of 50 sheep and cattle grazed parcels. Additional livestock information for the agricultural year of 2018/2019 was obtained during systematic visits to targeted parcels (see “Grazing Regime” section from Methods). We ended up with 23 cattle parcels and 27 sheep parcels (Fig. 1).

Bird and vegetation data

Breeding birds were sampled twice in each parcel during 7–16 April and 1–15 May 2019 respectively, always by the same observer (R.F.R). This was done to take into account species-specific breeding phenology in the area (early and late breeders)¹⁷ and minimize bias due to other factors (like weather or disturbance). Sampling was conducted using standardized 10 min point counts¹⁸ carried out at the central point of the parcel (Fig. 1). As the open terrain allowed for high visibility, a large detection radius was used, and all birds detected within 100 m of the central point were identified and counted. This radius is roughly similar to the one previously used for characterizing bird populations in the region¹⁹. All counts were carried out in the first four hours after sunrise and in the last two hours before sunset, with none in heavy or persistent rain, or in strong wind conditions. To estimate bird species richness and occurrences in each parcel, we pooled the data from the two counts. Species-level analyses focused on the six most common species, which occurred in > 30% of the parcels (see Supplementary Table S1). In addition to presence/absence, we also estimated population densities, using the count which yielded the highest estimate of density for each species (assuming this is the best indicator of population density, given the potential phenology and detectability biases above mentioned). Bird densities were based on the number of males simultaneously detected and expressed as breeding pairs/10 ha or males/10 ha (in the case of Little Bustard Tetrax tetrax and Common Quail Coturnix coturnix). Categorization to the genus level was made for the Crested and Thekla larks (Galerida cristata and G. theklae) due to difficulties in correctly identifying all individuals of these two very similar species in the field.

Vegetation height and cover were measured once in each parcel, between April 22 and May 6. Vegetation height was estimated in a set of ten 3 m radius plots defined inside the 100 m buffer (Fig. 1). In each plot, ten measurements of vegetation height were taken at random locations, for a total of 100 measurements per parcel. Vegetation height was measured using a 50 cm ruler and was defined as the highest point of vegetation projection within 3 cm of the ruler²⁰. All values were estimated to the nearest half centimeter. When no vegetation was present (bare soil, soil litter, rocks or animal dung) the height was set to zero (0) but these measurements were not considered to estimate the mean height of the sward. Vegetation cover was measured inside a 50 × 50 cm quadrat placed at each of the ten grid points, by visual estimation to the nearest 5% of the percentage of the quadrat area covered by vegetation²¹ (Fig. 1). Vegetation height and cover measurements were averaged within each parcel.

Grazing regime

The number and type of livestock in each parcel as well as the extent of the grazing period since the start of the year (2019) were gathered from interviews (Supplementary Information S1) to land managers during 1–15 May 2019. This information was further validated, and corrected in a few cases, through field checks during regular visits (made at two-week intervals) to the parcels (see “Bird and vegetation data” section from Methods). Three grazing regime indicators were estimated for the whole period (January–May 2019): livestock type (either sheep or cattle), animal density, and grazing pressure. The animal density in each parcel was calculated as the average density (animals per hectare) of any species (regardless of being sheep or cattle) that grazed the parcel during the 5-months period. Stocking Rate translated animal density into livestock unit (LU) per hectare (LU/ha), between January and May, according to the following criteria: one adult bovine = 1 LU; bovine aged < 6 months = 0.4 LU; one adult sheep = 0.15 LU²². Using LUs allows the comparison of densities across livestock types after correcting for their relative feeding requirements²³. Grazing Pressure was estimated as the Stocking Rate times the number of days a number of Livestock Units (LU) spent in a plot (LU/ha × number of days)²⁴. The area used for these estimations corresponded to the available area where animals could freely roam, which in many cases was larger than the sampled parcel area, which was often not delimited by fences. The number of days in the parcel was collected mainly from the interviews. However, in some cases the extent of grazing period was expressed qualitatively and thus had to be inferred, from common expressions according to the following criteria: ‘few’ = 5 days, ‘some’ = 10 days, ‘a fortnight’ = 15 days, ‘many’ = 20 days, ‘almost all month’ = 25 days²⁵.

Data analysis

Five explanatory variables describing grazing regime and sward structure (Table 1) were used as predictors of bird species richness, occurrence and abundance at parcel level. The correlation and multicollinearity between them were tested and all presented values of r < 0.70 and of variance inflation factor (VIF) smaller than 3²⁶.

Univariate differences between sheep and cattle parcels, both for bird response variables (species richness, density and occurrence) and for the explanatory variables, were tested using Generalized Linear Models (GLMs). We used a Gaussian error distribution and an identity link for quantitative variables, and a binomial error distribution and a logit-link function for occurrence data²⁶.

A Structural Equation Modelling (SEM) approach was then used to investigate how grazing regime directly and indirectly affects the occurrence and density of birds. SEM are probabilistic models that hypothesize a causal network with multiple variables that can appear as both predictor and response variables²⁷, allowing to look at both direct and indirect effects. We performed a confirmatory-exploratory path analysis²⁸ in the form of a piecewise SEM conducted in the R software²⁹, using the package “piecewiseSEM”²⁷. In piecewise SEM the network is translated as a set of linear equations which can be evaluated individually, using R²³⁰. The goodness-of-fit of the entire model was quantified by a directed separation test (“d-separation test”), which tests the assumption that all variables are conditionally independent, i.e. that there are no missing relationships among unconnected variables^27,30.

We started by building a theoretical model of our system (Fig. 2) based on previous literature and knowledge about birds and grasslands (See Supplementary Information S2 for more details on model construction). In short, the model states that the impacts of grazing regime on birds can occur: (A) indirectly, via the impacts of grazing pressure and potentially associated (non-measured) management decisions (e.g. fertilizer use or pasture improvement) on vegetation structure (vegetation height and cover); (B) directly, through the effect livestock-specific (sheep or cattle) behavior (trampling patterns, impacts of feeding mode on food resources for birds, potential egg predation) on birds; or, (C) directly through the disturbance impacts of animal density, expressed as number of herbivores spread over the area, irrespective of livestock type, on birds (Fig. 2). We considered all paths as significant if they had a p-value < 0.1. This threshold was used assuming it could indicate the existence of an effect, even if not significant at the traditional 0.05 level given the relatively low sample sizes. Other authors have used a similar approach in a SEM context (e.g. Sanz-Pérez et al.¹¹). Moreover, we used the d-separation test from piecewise SEM output to evaluate our theoretical model and identify eventual significant paths not considered initially.

After obtaining the final model for species richness and for the occurrence and density of each species, we estimated the standardized model parameters (expressed as mean ± standard error; SE) of causal effects. Effect estimates were used to calculate the strengths of direct and indirect effects between variables in the system. Indirect effects were described as a predictor variable (P1) having an effect on the response variable (R) through a simultaneous response and predictor variable (P2), P1 → P2 → R³¹. All statistical analyses were performed within “R” software environment, version 4.0.2²⁹.

Source: Ecology - nature.com