Data & codes for "Changes in abundance and distribution of European forest bird populations depend on biome, ecological specialisation and traits"

1.    Selection of European forest bird species and classification of their biome preferences

We selected all species that are related to forest and woodland based on two data sources: Storchová & Hořák (2018) and Tobias et al. (2022), resulting in 107 bird species studied (Data S1). We defined forest bird species as those using environments ranging from closed-canopy forests to more open-canopy woodlands (A. Lehikoinen & Virkkala, 2018; Storchová & Hořák, 2018; Tobias et al., 2022). We determined their biome specialisation using breeding distribution centroids and the overall breeding distribution of each of the species, using the global map of terrestrial ecoregions from Olson et al. (2001) and range data from European Breeding Bird Atlas 1 and 2 (Hagemeijer & Blair, 1997; Keller et al., 2020). We categorised species as Mediterranean, temperate, or boreal based on their predominant biogeographic region. We considered species commonly occurring over several biomes as “generalists”. For instance, we reclassified the two typically boreal species Glaucidium passerinum Linnaeus and Strix uralensis Pallas as “generalists” due to significant range expansions into central and southern Europe in recent decades, therefore no longer restricted to the boreal region. For the complete list of species, biome specialisation, traits, and specialisation indices, refer to Data S1.

2.    Changes in abundance and distribution of European forest bird species

We assessed long-term changes in European forest bird populations through two approaches: (i) changes in estimated total European-level species abundance over a 40-year timeframe; and (ii) changes in species spatial distribution over a 30-year timeframe (Fig. 1).

We utilized the estimated trends in European-level population size (i.e., the total number of individuals) for each common native European bird species from 1980 to 2017, as reported by Burns et al. (2021). Three species out of the 107 studied forest species were missing in the original manuscript and we used data generated with the same method from 1980 to 2018 from the European assessment, Article 12 (https://nature-art12.eionet.europa.eu/article12/). These abundance trends were calculated by Burns et al. (2021) using multi-sourced annual times series. For each species, they gathered population estimates and trends from each European country as well as European Union (EU)-level population trends. They analysed these data with a Bayesian hierarchical model to reconstruct EU-level smoothed species population time series. The model outputs include an average annual rate of abundance change and an associated 95% credible interval (Burns et al., 2021). Therefore, we did not directly use the average annual rate of abundance change, as this would have led us to consider species with low uncertainty as similar to those with high uncertainty. To account for the uncertainty, we categorised species as (i) declining, i.e., annual rates below one, (ii) increasing, i.e., annual rates above one and (iii) stable, i.e., annual rate whose 95% CI overlap one, i.e., no significant change. To better acknowledge the magnitude of the abundance change, significant changes with rates below 0.98 were labelled as “strongly declining” (i.e., 6.5% of the 107 species), while those above 1.02 were labelled as “strongly increasing” (i.e., 11% of the 107 species). To evaluate the sensitivity of the decision to categorised abundance change data, we also analysed abundance trend as continuous variable (see Supporting Information Fig. S8).

To determine changes in species distributions, we used a comparison of species distributions between two periods (i.e., 1985-1988 and 2013-2017) using the European Breeding Bird Atlas 1 and 2 (EBBA 1 & 2; Hagemeijer & Blair, 1997; Howard et al., 2023; Keller et al., 2020). Howard et al. (2023) provided calculations of observed colonisation and extinction areas at a 50 x 50 km resolution across Europe. We measured changes in range as the difference between colonisations and extinctions of each species, with negative values indicating contracting ranges and positive values indicating expanding ranges. Additionally, we calculated the shift in the centre of gravity of the distribution range between the two periods, as a distance (km) along the south-north gradient for each species (Howard et al., 2023).

3.    Trait and specialisation data for European forest bird species

We extracted data for six functional traits from several sources (Table 1). (i) The species temperature index (STI)represents the long-term average temperature within the species’ breeding range (A. Lehikoinen et al., 2021). (ii) Diet data during the breeding season were obtained from Storchová & Hořák (2018), classifying species into binary variables as vertebrate carnivorous, invertebrate carnivorous, and herbivores (combining the leaf and seed eaters). Storchová & Hořák (2018) classified species into a diet category when the corresponding food resource represented at least 10% of the species diet throughout the breeding season. Therefore, one species can be in several categories (i.e., omnivores). (iii) We obtained nesting site data from Pearman et al. (2014), classifying species into binary variables as ground nesters, tree hole nesters, or elevated nesters (> 1 m in a tree or shrub). We also included data on (iv) species dependence on old-growth forests (Data S1; mostly from Fraixedas et al. (2015) and Mönkkönen et al. (2014), if present on both references, we classified them as “1” and if only in one reference as “0.5”), (v) migration distance (Howard et al., 2023), and (vi) body mass (Tobias et al., 2022).

Finally, we extracted and developed seven species specialisation indices. (i) We used an overall specialisation index based on multiple traits (i.e., temperature, diet, foraging behaviour and substrate, habitat, and nesting site), and (ii) a nesting specialisation index, both obtained from Morelli et al. (2019). Both indices represent species specialization based on the dispersion of trait preferences for each species: e.g., nesting specialism equal 0 for species that nest in all habitat type and equal 1 for species that nest in only one habitat type). They are both calculated using the Gini index of inequality, which measures overall dispersion across, e.g., all traits for the overall specialization, based on data from Pearman et al. (2014) and Storchová & Hořák (2018). For additional information, see Morelli et al. (2019). We also used (iii) the diet specialisation index, (iv) the species distribution range during the breeding season (hereafter “breeding range area”) and (v) the climatic niche breadth from Reif et al. (2016). The diet specialisation index was calculated as the coefficient of variation for diet preferences for each species, where high values denotes specialized species (Reif et al., 2016). The breeding range area was evaluated as the number of 50-km squares in the distribution maps in Europe occupied by each species during the reproduction period, and is based on EBBA 1 (Hagemeijer & Blair, 1997). The climatic niche breadth was calculated as the difference between the 5% hottest and the 5% coldest mean temperature between April and June in which each species occurs, using EBBA 1 (Hagemeijer & Blair, 1997; Reif et al., 2016).

Additionally, (vi) we calculated a broadleaf forest specialisation index based on binary forest habitat preferences (Storchová & Hořák, 2018), assigning values of one for species found only in broadleaf forests; zero for those in coniferous forests, and 0.5 for those found in both. Lastly, (vii) we created a forest specialisation index based on the species habitat preferences (Storchová & Hořák, 2018). The forest specialisation index was calculated as the mean of species affinity across habitats. We used increasing habitat weights along a gradient of tree dominance: open habitats as 1, shrubland as 1.5, woodland as 2 (i.e., species associated with habitats structured by trees in lower density than in forest), forest generalist (found in both coniferous and broadleaf dense forests) as 3, and forest specialist (found only either in coniferous or broadleaf dense forests) as 4. For instance, the index value for species occurring either in shrubland, woodland or both broadleaf and coniferous forests is 2.167.

4.    Data analysis

Data analyses were conducted with R software version 4.4.1. (R Core Team, 2024). Given the non-independence of species due to their genetic relatedness, we accounted for interspecific phylogenetic distance in all models. We constructed the phylogenetic tree for the 107 European forest bird species using ‘rotl’ and ‘ape’ R-packages (Michonneau et al., 2022; Paradis et al., 2023). We used rotl as an interface with the "Open Tree of Life", employing tol_induced_subtree R-function to generate the phylogenetic tree and compute.brlen R-function to set branch lengths using Grafen’s computation. We generated separate phylogenetic trees for boreal (17), temperate (15), Mediterranean (16) and “generalist” (59) species to perform biome-specific analysis (see Supplementary Information, Figs. S1 & S2).

To investigate the effects of functional traits and specialisation indices on abundance, range changes, and distribution shift, we used two regression methods. All methods were based on the relationships between a measure of change and a functional trait or specialisation index. Our sample unit is an individual forest bird species (i.e., one value for each species, either abundance or range change, or distribution shift). Abundance change was a categorical variable (i.e., strong decline – decline – stable – increase – strong increase), while range change (i.e., difference between colonisation and extinction) and distribution shift (i.e., south-north shift) were continuous variables. Therefore, to study abundance changes, we used proportional-odds linear mixed effects model using (Phylo)clmm R-function from the ‘ordinal’ R-package (Christensen, 2022). Interspecific phylogenetic relatedness was included as a random effect, reflecting the correlation between species based on phylogenetic distances (see also Hagge et al. (2021) and Seibold et al. (2015)). For distribution changes, we employed phylogenetic generalised least squares regression (PGLS) using the gls R-function from the ‘nlme’ R-package (Pinheiro et al., 2023). The phylogenetic correlation structure was integrated into PGLS using Pagel’s lambda parameter (λ; Pagel (1999)) a widely used measured of phylogenetic signal strength (see, e.g., Hagge et al., 2021; Triviño et al., 2013).

Furthermore, we included latitude, a key driver of bird communities at broad scales (Luoto et al., 2007), as a fixed covariable (centroid latitude of the species’ breeding distribution) in all global models (i.e., species from all biomes together), except for the STI model due to strong correlation. For biome-specific analysis, we included latitude only in boreal species models for range change and distribution shift, as it significantly improved model fit (ΔAIC < -2). We did not add latitude for models specific to temperate, Mediterranean, and generalist species since it did not improve model fits (ΔAIC > -2). Additionally, we included breeding range area in range change and distribution shift models, assuming that species with larger ranges would exhibit larger shifts. We scaled predictors to a mean of 0 and standard deviation of 1 to facilitate effect size comparisons. We adjusted p-values using the Holm method (for n=3) to account for multiple testing of traits and specialisation indices on three response variables.