2. MATERIAL AND METHODS
2.1 Study area
Most of the studied Ortolan Bunting subpopulations, hereafter singing groups, were located along the coast of Finland and in the southern parts of the country, where agriculture in general is concentrated in Finland. Few study sites were in the eastern and Central Finland (Figure 1). The study area represents the northernmost distribution range limit of the species.
2.2. Bird data
The territories of the Ortolan buntings were mapped annually using a two-visit mapping method, which has been shown to be effective and reliable for mapping the species (Tiainen et al., 1985). The visits were made in May and the first half of June. Special attention was paid to simultaneous observations of singing males, to the territory locations, and to accurate estimation of territory numbers in territory concentrations. For each observed territory, its centre coordinates were recorded. In total 4,430 territories were recorded over the 19-year study period, from 2000 to 2018.
2.2.1 Forming singing group units
To form a study unit, we aggregated observations of territories from years 2000–2017 (observations from year 2018 were added in after the forming of singing groups) based on their proximity and surrounding landscape. Territories that were less than 500 metres away from each other, were considered to belong to the same singing group, i.e., subpopulation. Thus, we drew a 250-meter buffer around each territory observation within a singing group, merged the overlapping buffers and interpreted the outline as delineating the area occupied by that singing group (Figure 2). Ortolan Bunting territories have an average radius of about 100–300 meters (Vepsäläinen et al., 2007).
This kind of automatic forming of singing groups based on distance only produced sometimes very large singing groups which, in reality, might have consisted of several smaller groups separated by some natural boundary or other unfavourable matrix, such as a narrow stretch of forest or a busy road. Sometimes separate smaller singing groups were formed, when, in reality, they might have belonged to the same group despite the distance between the observations (or despite the failure to observe individuals between the observations). Therefore, in few cases a spatially very large singing group was divided into two or several smaller ones, or a few small and separate but adjacent singing groups were merged into one large group. When splitting a large group, the area along the split was proportionally divided to the newly formed singing groups. This ensured that a certain habitat area was not covered twice in separate singing groups. These exceptions (splits and merges) were always based on the observer’s perception on the field.
In total we formed 277 singing groups. Each singing group had an individual identification code (ID) and its centroid coordinates were recorded. Additionally, we recorded the region (REGION) and the municipality (MUN) within which a singing group was located (seven regions and 65 municipalities in total).
For each singing group we counted the annual number of territories (TERRI). For some years and some singing groups this information is missing, as not all sites have been visited annually, or because new locations (singing groups) have been added to the monitoring scheme during the study period. In total, the data consists of 1,474 annual territory mappings of singing groups.
2.3 Habitat and climate data
For each singing group area, we defined variables describing the habitat they occupied. From the Finnish topographic database (National Land Survey of Finland, 2023) we counted the length of roads, ditches, main drains and riverbanks (in metres) and the number of buildings within the singing group area. As there are little annual changes in these measures, we used data from years 2005, 2010 and 2016 to cover the years 2000–2007, 2008–2012 and 2013–2018, respectively. We then divided each of these measures by the area of the singing group to form a variable describing the density of each small-scale element (AROAD, ASTREAM, ARIVER, ABUILD).
By combining data from the Finnish land and crop parcel registers (Finnish Food Authority, 2023) we counted the total area of crop species grown within a singing group area. For years 2000–2016, in cases where more than one crop species was grown per field parcel, the data were unspecific about the precise location of different crop plants within a field parcel. Therefore, in cases where a field parcel with multiple crop plants was only partly within the singing group area, the areas of crop species were estimated based on their relative proportions within the field parcel. In cases where crop data were untraceable or not available, we completed the data from field notes, when possible, otherwise the crop was classified as “unknown species”.
Based on plant type or growth form we classified crop species into ten crop types: spring cereals (proportion of all observed crop types 0.54), grasses (0.26), winter cereals (0.06), oil plants (0.06), open-ground vegetables (0.03), legumes (0.02), herbs (0.01), unknown species (0.01), fruits and berries (0.005) and bioenergy crops (mainly reed canary grass, 0.003). For each singing group area, we then counted the proportions of each of these crop types and based on those, a Shannon–Wiener diversity index of crop types (SHANG).
Crop types were further classified into two groups based on how much vegetative cover the plants form during the start of the Ortolan Bunting breeding season in early May. Crop types which provide only minor cover and leave substantial amounts of bare ground on the field, mainly as a result of the associated ploughing practice, were classified as ‘bare’ (spring cereals, open-ground vegetables, oil plants and legumes), and crop types that provide substantial vegetative cover were classified as ‘cover’ (grasses and winter cereals, autumn-sown oil-seed rape). We then counted the proportion of ‘bare’ crop types within a singing group area (BARE).
To assess the interconnectedness of the agricultural landscape surrounding the singing group, we created a 5-km buffer around the singing group centroid and calculated the proportion of agricultural land cover (AGRI5) within the buffer using the Finnish Corine land cover data provided by the Finnish Environment Institute (Finnish Environment Institute, 2023). We used data from years 2000, 2006, 2012 and 2018 to cover the years 2000–2003, 2004–2009, 2010–2015 and 2016–2018, respectively.
From the Finnish meteorological data (Aalto et al., 2016) we calculated the mean daily temperature (TEMP) and precipitation (PREC) of the previous summer (21. May–15. July).
2.4 Statistical analyses
To focus on the multiplicative population rate of change of Ortolan Buntings, we first delimited our data to only those cases that had data from two consecutive years, i.e., with no gap between annual observations. Then we created a variable describing the number of territories in a singing group during the previous year (TERPRE), and further removed from the data those observations, where TERPRE was zero (no territories observed the previous year).
These limitations reduced the number of observations (annual territory counts of singing groups) to 678. The number of individual singing groups also reduced to 238. The number of consecutive visits (annual territory counts) per singing group varied from 1 to 13 (mean 2.8) and the number of territories per singing group varied from 0 to 26 (mean 4.25). Singing group areas ranged from 0.19 to 2.98 km2 (average 1 km2, SD = 0.73). The spatial and temporal extent of the data remained approximately the same after the data limitations. However, from Central Finland and North Carelia only three and one observations, respectively, remained and they were, in both cases, from only one singing group.
To model the change in number of Ortolan Bunting territories per singing group, we applied generalized linear mixed models with logarithmic link-functions. In the statistical analysis we first did model selection based on the model with the simplest set of fixed effects, to determine the most parsimonious approach for describing randomness, i.e., the random effect structure and error distribution – to be further used for studying the drivers of the decline as fixed effects. We compared 12 different models, each with only the categorical variable region (REGION) as a fixed effect and number of territories in previous year (TERPRE) as a log-transformed offset, and either i) no random effect (i.e., resulting in a GLM rather than a GLMM), ii) ID as random effect, iii) MUN as random effect, or iv) ID nested within MUN as random effect. These four model combinations were evaluated using three error distributions: i) Poisson, ii) COM-Poisson and iii) negative binomial; all of which show different dispersion patterns or relationships between the mean and the variance. Option ii) accommodates both over- and underdispersion, while option iii) can model overdispersion only. As singing groups were sampled multiple times, we wanted to test applying a mixed model in which singing group is used as random intercept, as this models a dependency structure among observations (variation in the average rates of change) from the same singing group. Including municipality as random intercept controls for possible spatial correlation in the average rate of change between singing groups which are located near to each other. In the fitting of the models, we used restricted maximum likelihood estimation (REML) in order to acquire estimators for the variance terms which are not biased.
We ranked the candidate models by the Akaike information criterion with a correction for small sample sizes (AICc) to evaluate their relative fit with data (Burnham and Anderson, 2002) and chose the model with the lowest AICc value as the best model structure for further analysis of fixed effects.
Based on the selected random structure (see Table S1 in Supplementary material), we studied which habitat variables explain the decrease in number of Ortolan Bunting territories within a singing group, by creating 16 nested models varying in their fixed effects only. All models included region (REGION) as fixed effect and number of previous year territories (TERPRE) as a log-transformed offset variable. In addition to these, the models contained all the possible combinations of the following groups of variables: i) weather (TEMP and PREC), ii) crop (SHANG and BARE), iii) small-scale structural element (AROAD, ASTREAM, ARIVER and ABUILD) and iv) landscape (AGRI5) variables (see Table 1 for full list of models). In the model selection for different fixed effects models, we used maximum likelihood estimation. Again, we ranked the candidate models by the Akaike information criterion with a correction for small sample sizes (AICc) to evaluate their relative fit with data (Burnham and Anderson, 2002) and chose the model with the lowest AICc value as the best model. Prior to running the model, all continuous variables AGRI5, PREC, TEMP, SHANG, BARE, ABUILD, ARIVER, ASTREAM and AROAD were scaled by subtracting the mean and dividing with the standard deviation.
All statistical analyses were done in the programming environment R program, version 4.3.1 (R Core Team, 2023). Models were fitted with glmmTMB-package (Brooks et al., 2017) and model selection was conducted with the MuMin-package (Bárton, 2023).