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.1. 欧洲森林鸟类物种遴选及其生物群系（biome）偏好分类 我们依托Storchová & Hořák（2018）与Tobias等人（2022）两份数据源，筛选出所有与森林及林地相关的鸟类物种，最终纳入107种研究鸟类（数据S1）。我们将森林鸟类定义为栖息环境涵盖郁闭林至开阔林地的物种（A. Lehikoinen & Virkkala, 2018; Storchová & Hořák, 2018; Tobias et al., 2022）。我们通过每个物种的繁殖分布质心与整体繁殖分布范围，结合Olson等人（2001）发布的全球陆地生态区地图，以及《欧洲繁殖鸟类图集1》和《欧洲繁殖鸟类图集2》（European Breeding Bird Atlas 1 and 2，简称EBBA 1 & 2；Hagemeijer & Blair, 1997; Keller et al., 2020）的分布范围数据，确定其生物群系特化程度。我们依据物种的主要生物地理区域，将其划分为地中海型、温带型或北方型；将广泛分布于多个生物群系的物种归类为“泛化种”。例如，我们将两种典型北方物种——花头鸺鹠（Glaucidium passerinum Linnaeus）和乌林鸮（Strix uralensis Pallas）重新归类为泛化种，因其近数十年已显著扩散至中欧与南欧，不再局限于北方生物群系。如需完整物种列表、生物群系特化信息、功能性状及特化指数，请参阅数据S1。 2. 欧洲森林鸟类的种群数量与分布变化 我们通过两种途径评估欧洲森林鸟类种群的长期变化：(i) 40年时间尺度内欧洲全域物种总种群数量的估算变化；(ii) 30年时间尺度内物种空间分布的变化（图1）。 我们采用Burns等人（2021）发布的1980年至2017年间欧洲本土常见鸟类的欧洲全域种群规模（即总个体数）估算趋势数据。本次研究的107种森林鸟类中，有3种在原始文献中缺失相关数据，我们采用欧洲评估报告第12条（https://nature-art12.eionet.europa.eu/article12/）中1980年至2018年采用相同方法生成的替代数据。Burns等人（2021）通过多源年度时间序列计算了这些种群数量趋势：他们收集了每个欧洲国家以及欧盟（EU）层面的种群估算值与趋势数据，并通过贝叶斯分层模型（Bayesian hierarchical model）重构了欧盟层面平滑后的物种种群时间序列。模型输出结果包含年均种群数量变化率及对应的95%可信区间（Burns et al., 2021）。因此，我们未直接使用年均变化率，否则会将不确定性较低的物种与不确定性较高的物种等同看待。为考量不确定性，我们将物种划分为三类：(i) 种群下降：年均变化率低于1；(ii) 种群上升：年均变化率高于1；(iii) 种群稳定：年均变化率的95% CI包含1，即无显著变化。为更清晰体现种群数量变化的幅度，我们将年均变化率低于0.98的显著下降归类为“剧烈下降”（占本次研究107种物种的6.5%），将年均变化率高于1.02的显著上升归类为“剧烈上升”（占11%）。为验证种群数量变化分类决策的敏感性，我们同时将种群数量趋势作为连续变量进行了分析（详见支撑信息图S8）。 为确定物种分布的变化，我们采用《欧洲繁殖鸟类图集1》和《欧洲繁殖鸟类图集2》（EBBA 1 & 2; Hagemeijer & Blair, 1997; Howard et al., 2023; Keller et al., 2020）中两个时期（即1985-1988年与2013-2017年）的物种分布对比数据。Howard等人（2023）提供了欧洲范围内分辨率为50×50 km的观测拓殖与灭绝区域计算结果。我们通过每个物种的拓殖区域与灭绝区域的差值衡量分布范围变化：负值代表分布范围收缩，正值代表分布范围扩张。此外，我们还计算了两个时期内物种分布范围重心的位移，以每个物种沿南北梯度的距离（km）表示（Howard et al., 2023）。 3. 欧洲森林鸟类的功能性状与特化指数数据 我们从多个数据源提取了6项功能性状数据（表1）。(i) 物种温度指数（Species Temperature Index, STI）：代表物种繁殖范围内的长期平均温度（A. Lehikoinen et al., 2021）。(ii) 繁殖季节饮食数据：源自Storchová & Hořák（2018），将物种划分为二元变量类别：脊椎动物食性、无脊椎动物食性、植食性（涵盖叶食与籽食类群）。Storchová & Hořák（2018）规定，当对应食物资源在繁殖季节占该物种饮食比例至少10%时，即可将其归入该饮食类别，因此一个物种可归属多个类别（即杂食性）。(iii) 筑巢位点数据：源自Pearman等人（2014），将物种划分为二元变量类别：地面筑巢者、树洞筑巢者、高架筑巢者（在树木或灌木中高度>1 m处筑巢）。我们还纳入了以下数据：(iv) 物种对原始林的依赖程度（数据S1；主要源自Fraixedas等人（2015）与Mönkkönen等人（2014）：若两份文献均记载该物种依赖原始林，则赋值为1；若仅单份文献记载，则赋值为0.5）；(v) 迁徙距离（Howard et al., 2023）；(vi) 体质量（Tobias et al., 2022）。 最终，我们提取并构建了7项物种特化指数：(i) 基于多性状（即温度、饮食、觅食行为与基质、栖息地及筑巢位点）的整体特化指数，以及(ii) 筑巢特化指数，二者均源自Morelli等人（2019）。两项指数均基于各物种的性状偏好离散程度衡量物种特化程度：例如，筑巢特化指数为0代表该物种可在所有栖息地类型筑巢，为1则代表仅能在单一栖息地类型筑巢。两项指数均采用吉尼不平等指数（Gini index of inequality）计算，该指数用于衡量整体性状离散程度（如整体特化指数基于所有性状），数据源自Pearman等人（2014）与Storchová & Hořák（2018）。更多细节请参阅Morelli等人（2019）。我们还使用了：(iii) 饮食特化指数；(iv) 繁殖季节的物种分布范围（以下简称“繁殖范围面积”）；(v) 源自Reif等人（2016）的气候生态位宽度。饮食特化指数通过各物种饮食偏好的变异系数计算，数值越高代表物种特化程度越强（Reif et al., 2016）。繁殖范围面积通过EBBA 1（Hagemeijer & Blair, 1997）的分布地图中，每个物种在繁殖期占据的50 km方格数量评估。气候生态位宽度通过EBBA 1的数据计算，为每个物种出现区域内4-6月平均温度的最热5%与最冷5%的差值（Hagemeijer & Blair, 1997; Reif et al., 2016）。 此外，(vi) 我们基于二元森林栖息地偏好构建了阔叶林特化指数（Storchová & Hořák, 2018）：仅分布于阔叶林的物种赋值为1，仅分布于针叶林的赋值为0，同时分布于两类森林的赋值为0.5。最后，(vii) 我们基于物种栖息地偏好构建了森林特化指数（Storchová & Hořák, 2018）：该指数通过物种对各栖息地的亲和度均值计算。我们采用沿树木优势度梯度递增的栖息地权重：开阔生境为1，灌丛为1.5，林地为2（即与树木结构相关但密度低于森林的生境），森林泛化种（同时分布于针叶林与阔叶林密林中）为3，森林特化种（仅分布于针叶林或仅分布于阔叶林密林）为4。例如，仅分布于灌丛、林地或同时分布于阔叶林与针叶林的物种，其指数值为2.167。 4. 数据分析 数据分析采用R软件版本4.4.1（R Core Team, 2024）。鉴于物种间因遗传亲缘关系导致的非独立性，我们在所有模型中纳入了物种间的系统发育距离（phylogenetic distance）作为校正。我们通过‘rotl’与‘ape’ R包（Michonneau et al., 2022; Paradis et al., 2023）构建了107种欧洲森林鸟类的系统发育树（phylogenetic tree）：以rotl作为“生命之树开放项目（Open Tree of Life）”的接口，使用tol_induced_subtree R函数生成系统发育树，并通过compute.brlen R函数采用Grafen算法设置分支长度。我们分别为北方型（17种）、温带型（15种）、地中海型（16种）及泛化种（59种）构建了系统发育树，以开展生物群系特异性分析（详见补充信息图S1与S2）。 为探究功能性状与特化指数对种群数量变化、分布范围变化及分布位移的影响，我们采用了两种回归方法：所有方法均基于变化量指标与功能性状或特化指数之间的关联关系。本次研究的样本单元为单个森林鸟类物种（即每个物种对应一个数值，无论是种群数量变化、分布范围变化还是分布位移）。种群数量变化为分类变量（即剧烈下降–下降–稳定–上升–剧烈上升），而分布范围变化（即拓殖区域与灭绝区域的差值）与分布位移（即南北向位移）均为连续变量。因此，为研究种群数量变化，我们采用‘ordinal’ R包中的（系统发育）clmm R函数进行比例优势线性混合效应模型（proportional-odds linear mixed effects model）分析（Christensen, 2022）：将物种间的系统发育亲缘关系作为随机效应，反映基于系统发育距离的物种相关性（亦可参阅Hagge等人（2021）与Seibold等人（2015））。为分析分布变化，我们采用‘nlme’ R包中的gls R函数进行系统发育广义最小二乘回归（Phylogenetic Generalised Least Squares, PGLS）分析：通过Pagel的λ参数（λ; Pagel, 1999）将系统发育相关结构纳入PGLS模型，该参数是衡量系统发育信号强度的通用指标（例如，可参阅Hagge等人（2021）与Triviño等人（2013））。 此外，我们将纬度（物种繁殖分布的质心纬度）作为固定协变量纳入所有全局模型（即涵盖所有生物群系的物种），但在物种温度指数（STI）模型中因存在强相关性未纳入该变量。纬度是大尺度鸟类群落的关键驱动因子（Luoto et al., 2007）。在生物群系特异性分析中，我们仅在北方型物种的分布范围变化与分布位移模型中纳入纬度，因其显著提升了模型拟合度（ΔAIC < -2）；而在温带型、地中海型及泛化种的专属模型中未纳入纬度，因其未提升模型拟合度（ΔAIC > -2）。此外，我们在分布范围变化与分布位移模型中纳入了繁殖范围面积，假设分布范围更大的物种会表现出更大的位移。我们将所有预测变量标准化为均值0、标准差1，以方便效应量比较。我们采用Holm法对p值进行校正（n=3），以校正针对3个响应变量的性状与特化指数多重检验带来的误差。