The future distribution of wetland birds breeding in Europe validated against observed changes in distribution

Wetland bird species have been declining in population size worldwide as climate warming and land-use change affect their suitable habitats. We used species distribution models (SDMs) to predict changes in range dynamics for 64 non-passerine wetland birds breeding in Europe, including range size, position of centroid, and margins. We fitted the SDMs with data collected for the first European Breeding Bird Atlas and climate and land-use data to predict distributional changes over a century (the 1970s–2070s). The predicted annual changes were then compared to observed annual changes in range size and range centroid over a time period of 30 years using data from the second European Breeding Bird Atlas. Our models successfully predicted ca. 75% of the 64 bird species to contract their breeding range in the future, while the remaining species (mostly southerly breeding species) were predicted to expand their breeding ranges northward. The northern margins of southerly species and southern margins of northerly species, both, predicted to shift northward. Predicted changes in range size and shifts in range centroids were broadly positively associated with the observed changes, although some species deviated markedly from the predictions. The predicted average shift in core distributions was ca. 5 km yr−1 towards the north (5% northeast, 45% north, and 40% northwest), compared to a slower observed average shift of ca. 3.9 km yr−1. Predicted changes in range centroids were generally larger than observed changes, which suggests that bird distribution changes may lag behind environmental changes leading to ‘climate debt’. We suggest that predictions of SDMs should be viewed as qualitative rather than quantitative outcomes, indicating that care should be taken concerning single species. Still, our results highlight the urgent need for management actions such as wetland creation and restoration to improve wetland birds’ resilience to the expected environmental changes in the future.


Introduction
Considerable effort has been invested in conserving biodiversity over recent decades. Yet, biodiversity losses continue at an unprecedented rate, as reflected by ongoing declines in population size and range contractions for many species worldwide (Pievani 2014, Tittensor et al 2014. The observed changes in the distribution of many species during recent decades have been primarily attributed to the ongoing rapid climate change, and to large-scale habitat loss (Brommer et al 2012, Reif and Flousek 2012, Gillings et al 2015, Hovick et al 2016, Pavón-Jordán et al 2019. Historical data clearly show that species may respond to climate and habitat changes by adjusting their spatial distributions (Parmesan et al 1999, Thomas and Lennon 1999, Brommer et al 2012, Littlefield et al 2017, Pavón-Jordán et al 2019. Therefore, it is recommended to consider climate as well as land-use variables to better describe drivers of species distribution changes (Newbold 2018).
Bird species that are ecologically dependent on wetlands are commonly used as indicators of wetland ecosystem health (Williamson et al 2013) and provide valuable ecosystem services such as food supply, pest control, seed dispersal, and cultural services such as recreation and hunting (Hamilton et al 1994, Teo 2001, Green and Elmberg 2014, Lehikoinen et al 2017. Still, many species of wetland birds have been declining worldwide and a subset has been classified as threatened species during the 20th century (Wang et al 2018, BirdLife International 2021. The expected changes in environmental conditions due to increases in global temperatures and changes in the land-use patterns that are likely to affect species distributions in the 21st century (IPCC 2014). Determining the expected change in range dynamics such as the direction and the magnitude of change in range margins and centroid allows for evaluating current networks and boundaries of protected areas with the possibility of moving from static to dynamic designs where the boundaries of protected areas change over time (Rayfield et al 2008, Cashion et al 2020. The range centroid is the centre of gravity of a distribution polygon and represents the core distribution of a species, where the abiotic conditions are assumed to be optimal for the species' biological and ecological functions (Sales et al 2020). Range dynamics may differ between the centroid and the margins but relatively few studies have considered the multiple changes in range characteristics (i.e. changes in range size, centroid, and margins). Huntley et al (2007) used a climatic-surface model on European birds to predict overall changes in range characteristics considering climate scenarios only and did not incorporate land-use scenarios. However, studies have shown that incorporating land-use information with climate information can significantly improve the predictive ability of species distribution models (SDMs) (Lee andJetz 2011, Sohl 2014).
Despite the surge in use of SDMs to predict future distributions during the last two decades (Newbold 2018, Soultan et al 2019, few studies have been able to use independent data to evaluate the predictive accuracy and temporal transferability of SDMs (Areias Guerreiro et al 2016, Barbet-Massin et al 2018). Nevertheless, the few studies available have reported interesting differences between the observed and predicted changes in species ranges, which provide new insights that will help to improve SDM methods (Virkkala and Lehikoinen 2014, Brun et al 2016. Here, we investigate the potential impacts of projected climate and land-use changes on the breeding distributions of 64 non-passerine wetland bird species in Europe, based on distribution data collected for the first European Breeding Bird Atlas (EBBA1, Hagemeijer and Blair 1997). We advance upon previous analyses for wetland birds in Europe (Huntley et al 2007(Huntley et al , 2008 by (a) incorporating landuse change scenarios together with climate change scenarios, (b) using ensemble SDMs, and (c) comparing predicted changes from the SDMs to the actual observed changes from the second European Breeding Bird Atlas (EBBA2, Keller et al 2020).

Methods
Species occurrences for 64 non-passerine wetland bird species that breed in Europe were obtained from the first Atlas of European Breeding Birds (Hagemeijer and Blair 1997), hereafter 'EBBA1' , which was compiled and published by the European Bird Census Council (EBCC). Appendix S1: Species data and study area (available online at stacks.iop.org/ ERL/17/024025/mmedia).
We modelled the breeding ranges of wetland birds by fitting ensemble SDMs using four commonly used presence-absence SDM algorithms (GLM, GAM, GBM, and RF) with default settings available within the 'biomod2' R package (R Core Team 2016, Thuiller et al 2016). SDM predictive performance was evaluated using the area under the curve (AUC; a threshold-independent metric) (Fielding and Bell 1997) (Mesgaran et al 2014) was used to assess the presence of non-analogue environmental conditions and to determine the degree of extrapolation (appendix S5: non-analogue environments and extrapolation). Last, the reference and future distribution ranges were classified into suitable and unsuitable ranges using a threshold that maximizes both model sensitivity and specificity (Liu et al 2013).
Three metrics were used to quantify the impact of environmental changes on the dynamics of breeding ranges for wetland birds: (a) per cent change in the area of the breeding range, (b) directionality and displacement shifts for the range centroid, and (c) latitudinal shifts of the northern and southern margins of the range (km yr −1 ). Changes in breeding range size were measured by calculating the range expansion (number of gained pixels; G) and range contraction (number of lost pixels; L), and relating them to the size of the reference range (total number of pixels; N) using 'BIOMOD_RangeSize' function in the 'biomod2' R package (Thuiller et al 2016). Directionality and displacement shifts of the geographic range centroid were quantified by delineating standard deviational ellipse (SDE) (Furfey 1927, Johnson andWilson 2009) over the reference and future ranges of a given species. As such, the centroid of SDE was used to represent species ranges' centroid. We quantified the directionality and displacement shifts in the range centroid by calculating the direction as a bearing relative to true north (0 • ) and the linear distances respectively, between the centroids of the reference and future ranges. SDE was calculated using 'calc_sde' function implemented in 'aspace' R package (Bui et al 2012), while both bearing and linear distance were calculated using 'bearing' and 'distGeo' functions, respectively, implemented in 'geosphere' R package (Hijmans 2019). It is expected that in case of expanding range size, the ranges of southerly species breeding in southern Europe might move northward, whereas the range of northerly species breeding in northern Europe might expand southward. Similarly, in the case of contracting range size, the ranges of southerly species would retract southward, whereas the ranges of northerly species would retract northward (Thomas and Lennon 1999, Kujala et al 2013). Therefore, based on the centroids of breeding ranges (i.e. the centroid of SDE), we classified our species into either northerly or southerly species if the breeding range's centroid was above or below the mean latitude of the study area of 5500 000 m (Thomas and Lennon 1999, Zuckerberg et al 2009. For the latitudinal shifts of southerly species, we measured the linear distance between the northern margin at the reference period and the predicted future periods for a given species (Ordonez andWilliams 2013, Carroll et al 2015). The northern margin was defined as the mean value of the upper 90% latitudes (90th percentile) of the pixels that were predicted suitable. For northerly species, we measured the linear distance between the southern margin at the reference period and the predicted future periods for a given species (Ordonez andWilliams 2013, Carroll et al 2015). The southern margin was defined as the mean value of the lower 10% latitudes (10th percentile) of the pixels that were predicted suitable.
Shifts in the latitudinal range margins are sensitive to original range size because small ranges can have larger potential shift (Williams and Blois 2018), and to natural barriers within the species' biogeographic regions such as Arctic ocean for northerly species. Therefore, to test whether a relationship exists between the predicted shifts in a range's latitudinal margins and the predicted changes in range size, we applied the approach developed by Thomas and Lennon (1999). We statistically estimated shifts in the southern margins of the northerly species and northern margins of the southerly species as the intercept of a regression line depicting the linear relationship between shifts in species range latitudinal margins and the changes in range size (Thomas and Lennon 1999, Taheri

Comparing predicted changes in species range with observed changes
The EBBA2 was recently published by EBCC (Keller et al 2020). EBBA2 is based on nationally collected data on breeding birds' distributions in Europe between 2013 and 2017 at a spatial resolution of 50 × 50 km grid cell and using the same methodological standards as for EBBA1. Comparisons of bird distributions collected during two time periods that were three decades apart (EBBA1: 1984, EBBA2: 2015) gave us a unique opportunity to evaluate and test predictions of SDMs. Our study objective was to compare the predicted changes in range size and range centroid from EBBA1 data with the observed changes (log10-transformed) calculated from EBBA2 data, assuming a constant rate of changes (linear) over the time.
We measured the displacement shifts of the range centroid by delineating SDE over EBBA1 and EBBA2 data of a given species. We measured the displacement shifts in the range centroid by calculating the linear distances between the centroids of the observed SDE of EBBA1 and SDE of EBBA2 data. The shifts in range centroids and changes in range sizes were calculated over different time scales, ∼30 years for the observed and ∼85 years for the predicted. Estimated shifts in range centroids were scaled to average annual shifts by dividing the observed and the predicted shifts in range centroids by the number of years, i.e. 30 and 85 respectively. We ran a linear regression to quantify the relation between the observed and predicted average annual shifts in range centroids. In the same way, we compared the predicted changes in range size with the

Changes in breeding range size
Our ensemble models predicted significant changes in the breeding ranges for most wetland birds under the projected future environmental conditions in Europe. Almost 75% of the species were predicted to contract their ranges, whereas ca. 20% of the species were predicted to expand their ranges (figures 1, 2 and table S2). The extent of the change in species range varied among species and according to the four different RCPs. For instance, ca. 25% and ca. 20% of the species were predicted to expand their breeding ranges by 2070 according to RCP 2.6 and 8.5, respectively. Four species with south-central distributions, little egret (Egretta garzetta), great white egret (Ardea alba), red-crested pochard (Netta rufina), and Kentish plover (Charadrius alexandrinus), were predicted to markedly expand their breeding ranges in the future (table S2). Other species such as common moorhen (Gallinula chloropus) and little grebe (Tachybaptus ruficollis) were predicted to maintain their reference breeding range in the future.
Range contractions were predicted for several northerly species, whereas almost all species that were predicted to expand their breeding range were southerly species (table S2). The pattern of change in range size was fairly consistent among the RCPs, with only a few species showing an inconsistent pattern of change such as black-crowned night heron (Nycticorax nycticorax) and red-throated diver (Gavia stellata) (table S2).

Shifts in centroids of breeding ranges
All species were predicted to shift their breeding range centroids, irrespective of the RCPs. A majority of species were predicted to shift their breeding range centroids in a northerly direction (ca. 5% NE, 45% N, and 40% NW) (figure 2 and table S3). The mean displacement shift in range centroid was predicted to be ca. 5 km yr −1 across 64 wetland birds. Appendix S6: Shifts in breeding ranges centroids.

Shifts in range margins
Both northern and southern range margins were predicted to shift northward. However, the magnitude of margin shifts was dependent on the species. For northerly species, shifts in their southern margins varied among the RCPs, with a mean displacement shift of ca. 2 km yr −1 ( figure 3 and table S4). For southerly species, the shifts in their northern margins varied from ca. 0-25 km yr −1 depending on the RCPs with a mean shift of ca. 6 km yr −1 (figure 3 and table 1).
Changes in range size were positively correlated with the predicted annual shifts in northern margins of southerly species, which suggest an increase in the number of suitable sites at northern boundaries (table 1). Changes in range size were negatively correlated with the predicted annual shifts in southern margins of northerly species, suggesting a decrease in the number of suitable sites at southern margins.

Comparing the predicted change in species range with the observed change
There was a significant positive association between the observed and the predicted annual changes in breeding range size and also the annual shifts in centroids ( figure 4). The predicted contractions of breeding range sizes were in general larger than what was observed in EBBA2 (intercept = −0.29 ± 0.15). However, some species were predicted to contract  their breeding ranges while they showed no change or a small increase in range such as the tufted duck (Aythya fuligula), pochard and whooper swan (Cygnus cygnus) (figure 4(a) and table S2). The predicted shifts in range centroids were on average greater (ca. 5 km yr −1 ) than the observed ones (ca. 3.9 km yr −1 ) (intercept = 3.06 ± 0.31; figure 4(b) and table S3). The differences in predicted vs observed shifts were largest for species with small observed shifts in distribution (figure 4(b) and table S3).

Discussion
The ensemble SDMs based on the expected changes in climate and land-use in the coming decades predicted significant contractions in the breeding ranges of many wetland birds, while only a few species were predicted to expand their breeding ranges. In general, most species distributions, as estimated by range centroids and range margins, were predicted to move northwards. The predicted shifts in range centroids were positively associated with the observed shifts in centroids over the 30 years (the 1980s-2010s) from EBBA2 data. Similarly, the predicted and observed changes in breeding distribution range size were positively related although some species displayed marked differences between predicted and observed changes. Our SDMs predicted: (a) considerable reductions in the size of the breeding ranges size (>50%) for many European wetland birds in the coming decades (figure 1 and table S2), (b) an average northward shift in breeding range centroids of ca. 5 km yr −1 (figure 2 and table S3), and (c) corresponding shifts in range margins with average displacement shifts of 2 and 6 km yr −1 for southern range margins of northerly In our study, species with wide southerly breeding distribution such as red-crested pochard, great white egret, and little egret were among those that were predicted to expand their breeding ranges in the future (table S2). The pattern of expansion for these species was also supported by the observed expansion reported by EBBA2 (Keller et al 2020). Species with broad distributions often encompass several sub-populations each with distinctive ecological characteristics and dynamics (Stockwell and Peterson 2002). Furthermore, such species are characterized by a wider environmental domain than they currently occupy, so they might benefit from new environmental conditions and, therefore, be able to expand their ranges (Stockwell and Peterson 2002, Koschová et al 2014). A second explanation for expansion of the southerly species could be that their ranges are not constrained by the continental border in the north (Koschová et al 2014).
About 75% of the modelled bird species were predicted to contract their breeding ranges in Europe in the future. For some species, such as longtailed duck (Clangula hyemalis) and common snipe (Gallinago gallinago), our SDMs predicted major contractions by 2070s. The magnitude of the predicted contractions (>50%) were consistent with results for many other birds at local (Andriamasimanana and Cameron 2013), regional (Harrison et al 2003, Virkkala et al 2008, and continental-scale (Barbet-Massin et al 2012, Langham et al 2015. The contractions were partly inconsistent with the observed changes from EBBA2 (Keller et al 2020) as many species including long-tailed duck and common snipe were observed to largely have almost the same range size in 2015 as 30 years earlier (table S2).
Some species also show a marked opposite pattern between predicted and observed range changes such as common merganser (Mergus merganser) and smew (Mergellus albellus) ( figure 4(a)). Large discrepancies may have been a result of some biotic factors not considered in our model. For instance, over the last decades, some species have strongly benefitted from the increased protection and conservation, intensified farming, and milder winters (Keller et  We focused on conditions during the breeding season but milder winters have benefitted the population sizes of several short-distance migrants that are wintering in central-north Europe (Musilová et al 2015(Musilová et al , 2018. Positive effects of wetland protection and mild winters could be possible explanations for predicted decreases but observed increases in range sizes for grey heron (Ardea cinerea), common goldeneye (Bucephala clangula), smew, and great cormorant (Phalacrocorax carbo) (table S2). Similarly, the divergence between the predicted expansion and the observed contraction in the breeding range of Kentish plover (table S2) could be attributed to the development in coastal breeding habitats (Montalvo and Figuerola 2006), and changed grazing pressure at coastal grasslands and increased predator populations (Keller et al 2020). Further, we assumed a constant linear rate of changes in breeding ranges over time due to the lack of data that can inform a better realistic assumption. For some species, the environmental predictors might not be able to capture the main niche dimensions of species. Examples are many fish-eating species such as goosander, smew and great cormorant that probably increased in numbers as a result of changed fish communities (Østnes andKroglund 2015, Frederiksen et al 2018), and large grazing birds such as whooper swan and common crane (Grus grus) that have increased due to changes in agricultural practices (Montràs-Janer et al 2020).
Why are most species predicted to contract their breeding range? First, the majority of the species that were predicted to contract their ranges are breeding in northern Europe, and thus constrained by the northern continental border (Gregory et al 2009, Koschová et al 2014. Second, the rate of climate change at northern latitudes could be faster as compared to that of the southern latitudes (Jetz et al 2007, Koschová et al 2014.
The northward shift of the southern margins was mainly driven by losing suitable sites at lower latitudes (significant negative range shift in table 1), while the northward shift of the northern margins was driven by gaining suitable sites at higher latitudes (significant positive range shift in table 1). A similar pattern has been found in several observational studies and has mainly been attributed to the latitudinal temperature changes (Thomas and Lennon 1999, Brommer 2004, Hitch and Leberg 2007, Kujala et al 2013, Ordonez and Williams 2013, Huang et al 2017. The predicted average displacement shift of breeding range centroids (5 km yr −1 ) is consistent with the average shift predicted in previous SDMs' studies (Huntley et al 2007, Russell et al 2015. Although most other SDMs' studies predicting a shift in range centroids suggest a shift towards the north, observational data from atlas inventories at country scale suggest these shifts to be smaller than predicted ( Our models probably overestimate the short-term impacts of environmental change because some of the inherent uncertainties associated with SDMs. A primary source of uncertainty in our study is the unaccounted factors such biotic interactions microclimatic conditions and species adaptability (Polaina et al 2021). A further source of uncertainty is the nature of EBBA2 data, which represent the transient distributions for many species including occurrences collected from old steady-state and newly colonized sites.
Our study calls for urgent intervention to preserve, manage, and restore the wetlands across Europe, which requires applying conservation measures at continental and national scales. We recommend to continue applying effective conservation measures such as wetland restoration and creation (Kačergytė et al 2021). Where the economic cost for restoring the natural wetlands is high, wetland creation is a potential alternative (Sebastián-González and Green 2016, Lehikoinen et al 2017). Additionally, previous studies showed that under effective governance including controlling bird hunting and restoring their potential habitats, wetlands can be refugia for wetland birds (Amano et al 2018, Kirby et al 2008). We recommend also applying spatial conservation planning, as it may inform the conservationists and decision-makers where to prioritize the conservation efforts.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// doi.org/10.15468/adtfvf.

Acknowledgments
We are very grateful to EBCC and all its contributors and GBIF for providing, curating, and making publicly available occurrence data from the EBBA1. We also acknowledge people and institutions that make public all the environmental data that we used in this research. Thanks to EBBA2 contributors for providing us with their recently published data. We also acknowledge the Scientific Project of the State Order of the Government of Russian Federation to Lomonosov Moscow State University No. 121032300105-0 for participating in EBBA2 data. Our research was funded through the 2017-2018 Belmont Forum and BiodivERsA joint call for research proposals, under the BiodivScen ERA-Net COFUND program, with the following funding organizations: the Academy of Finland (Univ. Turku: 326327, Univ. Helsinki: 326338), the Swedish Research Council (Swedish Univ. Agric. Sci: 2018-02440, Lund Univ.: 2018, the Research Council of Norway (Norwegian Instit. for Nature Res., 295767), and the National Science Foundation (Cornell Univ., ICER-1927646), and we also acknowledge the Swedish Environmental Protection Agency.