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).