Spatial gradients in country‐level population trends of European birds

Population trends reflect influence of environmental drivers acting upon species' population dynamics. As the strength of this influence may change predictably in space, we test multiple hypotheses about spatial gradients in the effects of environmental drivers on bird population trends across the continent.


| INTRODUC TI ON
Population trends are informative for detection of factors limiting species' population size, as well as for identification of species of conservation concern (Mace et al., 2008). Recently, several studies related population trends to species' traits (Gamero et al., 2017;Jørgensen et al., 2016;Koschova, Rivas-Salvador, & Reif, 2018;Seoane & Carrascal, 2008). By that means, we can infer how species' populations respond to respective environmental factors (Webb, Hoeting, Ames, Pyne, & Poff, 2010). For instance, population decline of forest species may indicate deterioration of forest habitat and population increase of resident species may reflect weakening effect of harsh winter weather on their population limitation due to warming climate (Bowler, Heldbjerg, Fox, O'Hara, & Böhning-Gaese, 2018).
In European birds, various patterns in relationships between population trends and species' traits had been identified such as steep decline of farmland birds (Gregory et al., 2005), divergence of trends between increasing warm-adapted and declining cold-adapted species (Devictor et al., 2012;Seoane & Carrascal, 2008), and population increase of species protected by legislation (Sanderson et al., 2016). However, such patterns were always investigated at the continental and sub-continental level of the whole Europe or by country-specific studies. Despite considerable achievements revealed by these studies, they were unable to fully account for spatial gradients in the impact of underlying drivers varying across the whole continent (but see Diaz et al., 2015;Jiguet et al., 2010). Due to differences in environmental conditions, in land-use histories and in levels of socioeconomic development among European countries, we suggest that following major spatial gradients in the influence of factors shaping species' population trends could exist: (a) climate change velocity is increasing from Southern to Northern and from Western to Eastern Europe (Loarie et al., 2009); (b) intensity of land use (especially of agricultural management) is increasing from Eastern to Western Europe (Tryjanowski et al., 2011), but large-scale abandonment occurred in Southern Europe (Hatna & Bakker, 2011) resulting in possibility of south-north gradient; and (c) enforcement and effectiveness of conservation legislation is increasing from Southern to Northern Europe and from Western to Eastern Europe (Donald et al., 2007;Sanderson et al., 2016).
Up to now, only a handful of studies took this spatial variability in the impacts of respective drivers into account when investigating the relationships between bird population trends and species' traits (Donald, Green, & Heath, 2001;Donald, Sanderson, Burfield, & van Bommel, 2006). Indeed, we are not aware of any study focusing on these relationships using a comprehensive set of traits and studying all these drivers in concert. Here, we address this important knowledge gap by focusing on country-level population trends across almost the whole continent (including its often neglected eastern nations such as Ukraine and Belarus).
By that means, we aim to investigate the key spatial gradients in trend-trait relationships which should mirror the variability in the impact of underlying drivers. We focus on the country level being aware that countries are human concepts and not natural biological units. However, country-level population trends are biologically meaningful at the same time because countries formulate economic and environmental politics which, in turn, shape bird populations (Gamero et al., 2017;Pe'er et al., 2017;Reif, Böhning-Gaese, Flade, Schwarz, & Schwager, 2011).
Based on the current knowledge of spatial gradients in drivers of bird population trends, we formulate and test following predictions.
(a) Population trends of open-habitat species will become increasingly negative towards the West of Europe due to the increasing agricultural intensity along this gradient (Tryjanowski et al., 2011). (b) The trends of open-habitat species will be increasingly negative towards the South of Europe due to widespread land abandonment in this region. (c) The trends of cold-adapted species should be more negative along the south-north and (d) west-east gradients due to higher velocity of climate change at higher latitudes and longitudes, respectively (Loarie et al., 2009;Wauchope et al., 2017). Similar gradients driven by the climate change velocity can be predicted for population trends of the long-distance migrants because the negative impacts of climate change increase with migration distance (Both et al., 2010;Jones & Cresswell, 2010;Salido, Purse, Marrs, Chamberlain, & Shultz, 2012). Therefore, we predict that (e) the trends of long-distance migrants will be more negative along the south-north gradient and (f) west-east gradient. Finally, we predict that (g) the species listed for longer time under Annex I of the European Union's Birds Directive will show more positive trends in Western and (h) Northern Europe due to the more effective enforcement of conservation legislation in the west (Donald et al., 2007;Sanderson et al., 2016) and a generally low pressure on bird populations in the North due to low human population density (Araujo, 2003).
The aim of this study is to test these four latitudinal and three longitudinal gradients in the relationships between population trends and species' traits for European birds. From Koschova et al. (2018), we excerpted information on population trends for 249 bird species in 32 countries (Table S1). These data are based on estimates of overall change in breeding population size (in %) of each species in each country from 2001 to 2012 reported in the European Red List of Birds (BirdLife International, 2015). For purposes of the Red List, the trend estimates were provided by national BirdLife Partners and result from national bird monitoring schemes, atlas mapping, species-specific surveys and expert assessments (BirdLife International, 2015 for exact specification of data source for each species in every country). Koschova et al. (2018) transformed these estimates to the scale from −2 (steepest possible decline) to 2 (steepest possible increase) to reveal symmetric values for population declines and increases using formula of Lemoine, Bauer, Peintinger, and Böhning-Gaese (2007). For this study, we used these transformed trends.

| Bird population trends
Since species' geographic ranges do not cover all countries, we obtained a dataset with 4,282 species-country combinations. We did not consider Russia for our analysis because we were interested in latitudinal and longitudinal gradients in country-level population trends, and extremely large area of this country relative to the latitudinal and longitudinal extent of European continent (Russia covers 79% and 45% of Europe's latitudinal and longitudinal range extend, respectively) makes Russian data on population trends uninformative in this regard. For each country, we calculated latitude and longitude of its centroid and these geographic positions were used in further analysis of spatial gradients (Table S2). Our calculation of countries' centroids considered only European mainland area, that is, we did not include, for example, Azores to calculate the centroid of Portugal. The only exception was countries confined solely to islands such as Ireland or the United Kingdom.

| Bird species traits
Species' traits were represented by 12 variables (Table S3)  because previous studies showed their important relationships to bird population trends, and their statistical effects should be taken into account even if they do not appear in the hypothesized spatial gradients. According to these previous studies, we predicted that population trends should be more negative: (a) in species with narrower habitat, climatic and trophic niches as they are less tolerant to the recent human-induced environmental changes than generalists (Bowler et al., 2018;Jiguet et al., 2006;Seoane & Carrascal, 2008); (b) in species with closer association to wetlands due to destruction of humid habitats (Amano et al., 2018); (c) in species not associated with human settlements because they cannot benefit from urbanization (Møller, 2009); (d) in species with faster life history strategies because they are forced to reproduce under suboptimal conditions due to their short life span (Kolecek et al., 2014); (e) in species with smaller breeding ranges and low migratory dispersion because such traits preclude effective buffering of environmental perturbations (Gilroy et al., 2016;Kolecek, Prochazka, Ieronymidou, Burfield, & Reif, 2018).

| Statistical analysis
We aimed to test for spatial gradients in the relationships of four traits to population trends taking the non-spatial effects of eight more traits into account. Due to extremely demanding computation time, it was not feasible to assess all combinations of these 12 predictor variables within a single analysis. Instead, we proceeded in two steps to reveal the subset of the variables used for the inference.
We first applied the information-theoretic approach to select the most relevant predictor variables from the eight traits for which we did not hypothesize any spatial relationships to trends (i.e., breadth of habitat niche along the forest open-landscape habitat gradient, position along the habitat humidity gradient, association with human settlements, diet niche position, life history strategy, climatic niche breadth, breeding range size and migratory dispersion). For this purpose, we employed linear mixed models (LMMs) with population trend as a response variable, the eight traits as fixed-effects explanatory variables, and taxonomic levels (i.e., family nested in order) and country as random effects. In LMMs, we considered only the main effects of species' traits since their interactions were not hypothesized. We constructed a global LMM containing all predictors and then assessed all their combinations contained in candidate LMMs using Akaike information criterion corrected (AIC c ) within the R-package MuMIn (Bartoń, 2018). We revealed the best subset of predictors based on the models with ΔAIC c < 2. These predictors were taken as covariates for the next step.
In the second step, we focused on four traits for which we formulated hypotheses about the spatial gradients in their effects on population trends (i.e., position along the forest open-landscape habitat gradient, climatic niche position, migration distance and time since listing under the Annex I). For this purpose, we employed a generalized additive mixed model (GAMM) using R-package gamm4 (Wood & Scheipl, 2017). In the GAMM, population trend was related to these traits in a three-way interaction of a given trait with countries' latitude and longitude. Their main effects and lower order interactions were also considered in the model, but not shown in the model output (Wood & Scheipl, 2017). We modelled smooth terms by cubic regression splines (Wood, 2017)  was substantially higher than number of EDF in all interactions (see Results), implying an adequate smooth (Wood, 2017). We assessed the statistical significance of each interaction, but we also plotted three contour plots for each interaction term showing changes of bird trends in space related to minimum, mean and maximum value of a given predictor. Moreover, the plots depict zones where 95% confidence intervals for trend estimates overlap zero (non-significant trends) and do not overlap zero (significant trends), respectively. In addition to the interaction terms, we included the traits selected in LMMs as the main effect covariates into the GAMM without defining their smooth terms. The random effects applied in GAMM were the same as in LMMs, that is, the hierarchical structure of taxonomic levels and country.
Before the modelling, we checked for possible relationships among all explanatory variables by calculating their variance inflation factors (VIF) and found no indications of elevated VIF for any of the variables (Table S4). Normality of model residuals was assessed by a visual check by qq-plots and indicated that the assumptions for statistical modelling were met. We also did not find any influential outliers after plotting the variables.

| Selection of traits without spatial patterns in population trends
The information-theoretic approach was used for selection of the best subset of the traits for which no spatial gradients in their relationships to bird population trends were hypothesized. By that means, we revealed that only a single LMM received a high support with ΔAIC c < 2 (Table S5; see Table S6 for random effects statistics). The top model contained two traits: habitat niche breadth and breeding range size (Table S5). Both trait variables had significant relationships with population trend: population trends became more positive with increasing niche breadth (coefficient = 0.03, SE = 0.01) indicating that species with broader habitat niche (i.e., habitat generalists) had more positive population trends than specialists, and population trends became more negative with increasing range size (coefficient = −0.04, SE = 0.01) indicating that species with larger breeding geographic ranges had more negative trends than species with smaller ranges. These two traits were added into the GAMM together with the traits for which spatial gradients in their effects on bird population trends were hypothesized (see the next section).
TA B L E 1 Results of a generalized additive mixed model relating bird population trends in European countries to species' traits in the interaction with countries' latitude and longitude: (a) estimates of the main effects for the two traits not included in the interaction terms, (b) significance of the three-way interactions testing for the spatial gradients in the effects of four more traits  Table S7 for random effect estimates. a EDFmax = maximum number of the effective degrees of freedom based on the dimensions of the basis used (k = 5). b EDF = number of the effective degrees of freedom.

| Testing for spatial gradients in population trends
By running the GAMM, we assessed the non-spatial relationships of population trends with the two traits selected by LMMs in the previous step, which remained qualitatively the same as in the top LMM (Table 1a), and spatial gradients in the relationships between population trends and four traits (Table 1b, see Table S7 for random effects statistics).
We found a significant three-way interaction between the time elapsed since listing the species under Annex I and countries' latitude and longitude (Table 1b) (Figure 1h). Open-habitat species showed a weak gradient of increasingly negative trends towards the Western Europe, but only the values for Ireland and the UK were significant (Figure 1i).
The interaction between the migration distance and the geographic position of European countries in the effect on bird population trends was insignificant (Table 1b). Accordingly, spatial gradients in the trend estimates were almost absent in species non-migrating or migrating for intermediate distances having the trends close to zero nearly everywhere in Europe (Figure 1j,k). In species with the longest migration distances, we found slightly, but significantly more negative trends in Western-Central Europe (Figure 1l).

| D ISCUSS I ON
Our analysis focused on population trends of 249 bird species in 32 European countries from 2001 to 2012. We uncovered spatial gradients in population trends for three of four species' traits considered: the time since listing under Annex I of the EU's Bird Directive showed strong and complex spatial patterns in trends and the same applied, albeit less significantly, for the species' climatic niche position; the habitat niche position showed weaker spatial gradients in population trends, and the spatial patterns in trends were almost absent for migration distance. Overall, the support found for majority of the focal traits indicates that the spatial variation in trends is indeed common at large spatial scales, even though the spatial resolution was relatively coarse in our dataset. This general finding corresponds to previous studies showing that the species' population trends often vary among regions (Böhning-Gaese, Taper, & Brown, 1994) and that some traits show contrasting relationships to population trends in different countries (Diaz et al., 2015;Donald et al., 2001Donald et al., , 2007Donald et al., , 2006Koschova et al., 2018;Sanderson et al., 2016).
In the case of the Annex I species, we found an interesting difference in spatial gradients in trends between the species listed for longer time (25 years) and shorter time (12 years). Whereas the species listed for longer time showed positive trends in majority of the continent and the positive values increased towards the north-west, the opposite, but much weaker pattern was found in species listed for shorter time. It seems that Annex I is generally successful, but that species particularly benefit from the effective enforcement of conservation legislation in Western Europe (Donald et al., 2007;Koschova et al., 2018;Sanderson et al., 2016) and perhaps also from generally lower human pressure upon natural habitats in the North due to low human population density (Sanchez-Fernandez, Abellan, Aragon, Varela, & Cabeza, 2018). However, delivering conservation benefits for protected species takes relatively long time (Male & Bean, 2005) and the initial population status of such species may be indeed worse in regions where conservation actions then take place (Donald et al., 2007). This may explain the difference in patterns between the shorter and longer listed species.
Even more conspicuous differences in spatial gradients in population trends were observed with respect to the species' climatic niche position. Cold-adapted species showed strong declines towards the north-west. This pattern corresponds with high ecosystem sensitivity to warming in North-Western Europe (Li, Wu, Liu, Zhang, & Li, 2018) and unfavourable population status of the upland and Arctic species caused by the climate change (Scridel et al., 2018). The rate of warming is particularly high in northern latitudes (Gilg et al., 2012) with adverse consequences for the whole ecosystem (Descamps et al., 2017). For bird populations, the most adverse impacts can be represented by the increased predation pressure (Kubelka et al., 2018) and depletion of food resources (Loboda, Savage, Buddle, Schmidt, & Hoye, 2018). However, the warm-adapted species showed increasingly positive trends in the north-western direction corresponding to widespread range shifts and population increases of such species (Devictor et al., 2012;Pecl et al., 2017;Stephens et al., 2016). It is in accord with much quicker population dynamics at the leading edge of species' ranges than at the trailing range edge which may not move at all (Gillings, Balmer, & Fuller, 2015). The mechanisms behind this pattern are still unclear but the different warming rates may be one of them (Virkkala & Lehikoinen, 2014).
In respect to the niche position along forest-open country gradient, spatial patterns were less significant, but still visible. Forest species showed more positive trends towards North-Eastern Europe with significant increases from Germany to Finland. This pattern corresponds to generally positive effects of forest management performed in these regions (e.g., forest maturation, heterogeneity in tree species composition and afforestation) on forest bird abundances (Ram, Axelsson, Green, Smith, & Lindstrom, 2017;Schulze et al., 2019).
Surprisingly, we did not find much evidence for the differences in trends of open-habitat birds between Western and Eastern Europe.
This pattern contrasts with more negative trends of farmland birds observed in Western than in Eastern Europe from 1970 to 2000 (Donald et al., 2001(Donald et al., , 2006, frequently explained by lower intensity of agricultural management in the east (Sutcliffe et al., 2015;Tryjanowski et al., 2011). Indeed, only in the westernmost regions Europe (Reif & Vermouzek, 2019), whereas the decline of Western European population can be already levelled out in the time period covered by the trend data used in our analysis (Inger et al., 2015). In turn, these processes can produce similarities in population trends of farmland bird species across Europe at the time scale we focused on.
However, we should also bear in mind that our descriptor of species' habitat niche using its position along the gradient from forest to open habitats may not represent only the farmland birds occurring in agricultural habitats, but also other open-habitat species largely independent on agricultural management including those breeding in rocks (e.g., Snowfinch Montifringilla nivalis) or in steppes and deserts (e.g., Pin-tailed Sandgrouse Pterocles alchata). Focusing solely on farmland species might reveal a stronger west-east gradient than the one detected by our analysis.
For migration distance, we did not confirm the hypothesized gradients of increasingly negative population trends towards areas of higher warming rates, that is in northern and eastern directions.
One reason may be that the effect of the climate change on migrants is not so strong at the breeding grounds and that their populations might be more limited by conditions at non-breeding grounds (Gilroy et al., 2016;Kolecek et al., 2018). Overall, the trends of species migrating for different distances did not show any spatial patterns with the exception of more negative trends of long-distance migrants observed in a part of Western Europe. If we accept that their populations are limited at non-breeding grounds, one speculative explanation offers the serial residency hypothesis (Cresswell, 2014).
According to this hypothesis, dispersion of juvenile birds over the non-breeding grounds buffers environmental perturbations and thus prevents population declines. Since western populations of these species migrate and winter closer to the Atlantic coast and have thus less space to spread than the eastern populations (Møller, Garamszegi, Peralta-Sánchez, & Soler, 2011), their dispersion is more limited and declines may be thus more frequent.
Besides the spatial gradients, we found two more traits related to bird population trends irrespective to countries' geographic location. First, generalists had more positive trends than specialists, most likely due to the sensitivity of specialized species to recent human-induced environmental changes, which can be, at the same time, favourable for generalists able to exploit novel opportunities (Bowler et al., 2018). Second, the species with smaller breeding ranges have more positive trends than the species with larger ranges. This pattern contrasted with our expectation of the buffering effect of the large range size preventing population declines.
However, it corresponds with generally more positive trends of the rare species reported in other pan-European studies (Inger et al., 2015). Since our models took the negative effect of the habitat specialization, which is typically high in the species with limited geographic distribution (Reif et al., 2016), into account (see above), the pattern is produced by the small range size per se. We suggest that this relationship is a consequence of recent range expansions of the rare species when such species colonize new areas undergoing exponential population growth and thus highly positive trends (Saether & Engen, 2002). These expansions may be driven by climate change as many species originally breeding at the southern edge of Europe or in Africa spread further north (Hagemeijer & Blair, 1997), such as Cattle Egret (Bubulcus ibis) and Cetti's Warbler (Cettia cetti).
In addition, successful conservation actions to rescue some species F I G U R E 1 Spatial gradients in bird population trends across Europe based on results of a generalized additive mixed model relating trends in European countries to species' traits in the interaction with countries' latitude and longitude (see Table 1b). Maps were produced for minimum, mean and maximum values of following traits: (a-c) time since listing under the Annex I of the EU's Birds Directive (0, 12.5 and 25 years); (d-f) climatic niche position expressed as the mean temperature in species' breeding range during the breeding season (3.83, 10.59 and 17.35°C); (g-i) habitat niche position along the gradient from forest to open habitats (forest species-1, mid-position species-4 and open habitat species-7); (j-l) migration distance (0, 3,950 and 7,899 km). Contour lines depict curves of constant values for trends and coloured areas show trend values from the more negative (in red) to the most positive (in green). Bright zones represent areas with trends significantly different from zero according to their 95% confidence intervals. Points are centroids of particular countries from regional extinction (e.g., Bearded Vulture, Gypaetus barbatus) and protection some other species from hunting (e.g., White-tailed Eagle, Haliaeetus albicilla) may have similar positive effects on the rare species trends (Kolecek et al., 2014).
Our models did not fully account for relatedness of the focal species by incorporating phylogeny into the random effects because it was not technically feasible due excessively long computation time.
We overcame this problem by considering hierarchical taxonomic levels in our models being aware that taxonomy does not fully reflect phylogeny. However, various studies showed that the species' abundances and population trends do not have any significant phylogenetic signal (McGill, 2008;Møller, 2008;Reif et al., 2010). Moreover, our study focuses on conservation issues arising mostly from contemporary modification of landscape by human activities, and not on evolutionary processes (Seoane & Carrascal, 2008). Under these circumstances, controlling statistical models for phylogeny does not seem necessary (Westoby, Leishman, & Lord, 1995;de Bello et al., 2015, for further discussions on this topic).
Another caveat of our approach results from using single species-specific values of particular trait variables for all countries considered. Although this approach is commonly used in studies focusing on population trends of European birds (Donald et al., 2007;Gregory et al., 2005;Jørgensen et al., 2016;Sanderson et al., 2016), expressing country-level trait values would likely reveal stronger relationships and perhaps also stronger spatial gradients. However, such data do not exist in comprehensive datasets (even the most up-to-date sources such as Storchova & Horak, 2018, worked at the pan-European level) and their collecting would be very difficult and prone to errors because different approaches to express the trait values are often applied in different country-level or populationlevel studies.
In summary, our study demonstrated that the relationships between bird population trends and species' traits change along various spatial gradients across Europe. One part of these changes can be interpreted in the light of spatially differentiated climate change impacts and associated species' responses. It would be interesting to test whether these gradients hold true also in other continents such as North America or Asia. The other factors include the differences in land-use (e.g., agricultural intensity and deterioration of open habitats) and conservation effectiveness. We suggest that taking these spatial gradients into account is crucial for formulation of European environmental policies. For example, the approach found to be effective in one region may fail to work in the opposite side of the continent (Baldi & Batary, 2011) due to the differences in the impacts of major driving forces described above. Our analyses identified the factors and groups of species where such risks exist, and we encourage policy makers to adopt these findings.

ACK N OWLED G EM ENTS
Manuscript benefited from the comments of three anonymous referees and the associate editor. The study was supported by Charles University, Prague (grant no. PRIMUS/17/SCI/16).

DATA ACCE SS I B I LIT Y
All data used in this study are provided in the supplementary online materials.