Modeling nonbreeding distributions of shorebirds and waterfowl in response to climate change

Abstract To identify areas on the landscape that may contribute to a robust network of conservation areas, we modeled the probabilities of occurrence of several en route migratory shorebirds and wintering waterfowl in the southern Great Plains of North America, including responses to changing climate. We predominantly used data from the eBird citizen‐science project to model probabilities of occurrence relative to land‐use patterns, spatial distribution of wetlands, and climate. We projected models to potential future climate conditions using five representative general circulation models of the Coupled Model Intercomparison Project 5 (CMIP5). We used Random Forests to model probabilities of occurrence and compared the time periods 1981–2010 (hindcast) and 2041–2070 (forecast) in “model space.” Projected changes in shorebird probabilities of occurrence varied with species‐specific general distribution pattern, migration distance, and spatial extent. Species using the western and northern portion of the study area exhibited the greatest likelihoods of decline, whereas species with more easterly occurrences, mostly long‐distance migrants, had the greatest projected increases in probability of occurrence. At an ecoregional extent, differences in probabilities of shorebird occurrence ranged from −0.015 to 0.045 when averaged across climate models, with the largest increases occurring early in migration. Spatial shifts are predicted for several shorebird species. Probabilities of occurrence of wintering Mallards and Northern Pintail are predicted to increase by 0.046 and 0.061, respectively, with northward shifts projected for both species. When incorporated into partner land management decision tools, results at ecoregional extents can be used to identify wetland complexes with the greatest potential to support birds in the nonbreeding season under a wide range of future climate scenarios.

posed by the natural spatial and temporal variations in the condition of wetland habitats. Impending climate change may alter inundation patterns and the function of critical wetland habitats via, for example, increased evapotranspiration associated with higher temperatures (Johnson, Werner, & Guntenspergen, 2015;Johnson et al., 2010;Niemuth, Fleming, & Reynolds, 2014;Sofaer et al., 2016). In addition, humans continue to alter the landscape, including the conversion of land to uses that produce food and fuel (Hurlbert & Liang, 2012;Thomas, Lanctot, & Székely, 2006). The combined effects of changing climatic conditions and additional habitat loss could pose a substantial future threat to the persistence of waterfowl and migratory shorebirds.
Wetland availability across the North American Great Plains (hereafter Great Plains) is affected by both factors, as inundation state depends on rainfall and runoff in combination with evapotranspiration rates and land cover (Bartuszevige, Pavlacky, Burris, & Herbener, 2012;Cariveau, Pavlacky, Bishop, & LaGrange, 2011). Weather can vary considerably across the large regions crossed by long-distance migrants (Millett, Johnson, & Guntenspergen, 2009;Tøttrup et al., 2008). The latest global climate models forecast more frequent extreme precipitation events, yet the southern Great Plains are expected to become drier as a result of increasing temperatures and evapotranspiration rates, leading to longer droughts .
Such changes can have large negative effects on freshwater wetlands (Kundzewicz et al., 2008), increasing erosion and wetland sedimentation during extreme events, reducing hydroperiods and the diversity of water regimes (Johnson et al., 2010), and changing landscape connectivity (McIntyre et al., 2014).
Playas and other wetlands within the Great Plains provide essential habitat for many wetland-dependent vertebrate species and are especially important as migration and wintering areas for shorebirds and waterfowl (Cariveau & Pavlacky, 2009;Haukos & Smith, 1994;Skagen, Sharpe, Waltermire, & Dillon, 1999). The density of wetlands holding water positively affects the abundance and richness of waterfowl and shorebirds (Albanese & Davis, 2013;Webb, Smith, Vrtiska, & LaGrange, 2010), potentially representing the most important habitat features on the landscape. A large proportion of the wetlands in the Great Plains are ephemeral, resulting in a temporally dynamic hydrology that can have large implications to bird migration and wintering success (Albanese, Davis, & Compton, 2012;McIntyre et al., 2014;Naugle, Johnson, Estey, & Higgins, 2001).
The Great Plains are relatively flat, interspersed with wetlands, and fragmented with cropland and ranchland. A large proportion of the remaining wetland habitat is privately owned and therefore vulnerable to the societal dynamics that drive land-use and land-cover decisions, including agricultural activities (Detenbeck et al., 2002). Wetlands within the Great Plains are threatened by sediment accumulation (Burris & Skagen, 2013), increased summer temperatures, changing precipitation patterns, and declining hydroperiods due in part to agricultural intensification.
Among the many applications, species distribution models have been used to determine environmental relationships and potential impacts of climate change (Elith, Kearney, & Phillips, 2010;Heinanen & von Numbers, 2009). Our objective was to model current and future distributions of several en route migratory shorebirds and wintering waterfowl relative to land-use patterns, the spatial distribution and composition of wetlands, and data from a wide range of global climate models. In general, our models estimated contemporary probabilities of occurrence and responses to climate change. We assumed that bird encounters can indicate the presence of suitable habitat and, therefore, predicted probabilities of occurrence could inform future conservation efforts.

| Study location
We modeled the distributions of several migratory shorebirds throughout spring migration, that is, stopover locations, as well as the winter, that is, December and January, distributions of two ducks (Table 1), within the boundary of the Great Plains Landscape Conservation Cooperative (GPLCC) (Figure 1), a public-private partnership that provides science assistance to natural resource managers within Bird Conservation Regions 18 and 19 (http://www.nabci-us.org/bcrs.htm). Relatively flat and once predominantly mixed-grass and shortgrass prairie, the GPLCC encompasses a large portion (>782,000 km 2 ) of the south-central Great Plains (Figure 1). Approximately one-half of the study area has been converted to agriculture, largely in the eastern regions.

| Observation data
We downloaded observation data from the eBird citizen-science database (Sullivan et al., 2009) for the states intersected by the GPLCC boundary for 2002-2010, years that coincide with the beginning of the eBird program and the most recent year in the empirical climate data that we used, respectively. We clipped data to the GPLCC boundary and applied several filters including one based on the scale at which shorebirds optimally respond to habitat patterns, that is, 1.25-2 km (Albanese et al., 2012;Cunningham & Johnson, 2006. We maximized the number of observations available by applying the 2-km threshold at several steps (see Fink et al., 2010;Sohl, 2014). We used all incidental observations, that is, from surveys with neither spatial nor temporal measures, and stationary counts, that is, from surveys with a known duration and across an area <30 m in diameter. For traveling and area counts, from surveys across known distances and areas, respectively, which also report to single point locations, we assumed a positive relationship between survey size and habitat heterogeneity and therefore only included traveling counts ≤4.023 km and area counts ≤1,257 ha (2-km circular radius).
An important component of eBird data is a field for birders to specify whether or not all species were recorded, information useful to the assumption of absences. We eliminated absences from incidental surveys, considering them less reliable than those from the other survey types. As with the presences, we also eliminated those from traveling counts >4.023 km and those from area counts >1,257 ha. Additionally, we eliminated duplicate absences, that is, surveys on the same date and at the same location as another absence or presence.
We added data from two surveys conducted between mid-March and mid-May in 2008, accounting for approximately 5% of the observations used to build models. One survey included roadsides within a random selection of townships first stratified on wetland and cropland areas. Roadsides were surveyed twice, and all shorebirds within 400 m were recorded. Surveys of 25 sites with known shorebird usage were also included. On the ground, these data were collected similar to those classified as stationary counts in eBird; data from all sources, eBird and supplemental, were converted to presence and absence categories.
We mapped observation locations by date, compared them to frequency histograms (Skagen et al., 1999), and selected March 1 to June 15 as the range for spring migration. The largest number of migrants occurs toward the middle of a season (Dunn, Hussell, & Adams, 1997).
To reduce sampling effects, we sorted data for each species by the survey date and clipped the smallest number of records that removed ≥5% from each tail, that is, the earliest and latest encounters within the spring migration window (see Farmer, Hussell, & Mizrahi, 2007).

| Predictor variables
Migration is an interaction between space and time; therefore, we modeled spring stopover locations in the GPLCC with variables of each type, including latitude and day-of-year, that is, 1-365; 60-166 thus represents spring migration, March 1 to June 15 (see Gowan & Ortega-Ortiz, 2014). Higher-order date variables have been important predictors elsewhere (see Dunn et al., 1997;Farmer et al., 2007), so we also included the square of the day-of-year term. Geographic coordinate variables can account for small-scale processes and address issues with spatial autocorrelation (Bailey & Gatrell, 1996;Sohl & Sayler, 2008).
Individual birds likely make decisions during migration reflecting both intrinsic and extrinsic factors such as body condition and proximal weather (Marra, Francis, Mulvihill, & Moore, 2005;Richardson, 1978). We elected to use weather data from Maurer, Wood, Adam, Lettenmaier, and Nijssen (2002) because they overlapped a large portion of the eBird data, 2002-2010, and were used to bias-correct the global climate model data. Avian migration studies have found important predictors in temperature (Huin & Sparks, 2000) and wind (Green & Piersma, 2003;Richardson, 1978). Additionally, precipitation is often key to wetland condition. For example, Bartuszevige et al. (2012) found that playa inundation is largely a result of the precipitation events in the preceding 2 weeks.
We time-matched observations to five month-long weather variables including total precipitation and averages of the daily wind speed and minimum, maximum, and average temperature (Table 2). We matched observations that occurred in the first 15 days of a month to weather data from the previous month and remaining observations T A B L E 1 Focal waterfowl and shorebird species-modeled, migration distance index, range of water depths used for foraging, and numbers of presences and absences in observation data Migration distance indices (× 1,000 km) for shorebirds, based on distances between breeding and wintering areas, are average of shortest distance, distance between midpoints, and distances between extreme edges of ranges (from Skagen & Knopf, 1993).
to the matching month. We also matched each record to the precipitation total from the previous calendar year, for example, all 2010 records, regardless of month, were matched to the total precipitation in 2009 for that location (see Johnson, Rice, Haukos, & Thorpe, 2011). In addition to phenological effects, weather variables can predict available wetland habitat because availability can vary as a function of precipitation, temperature, solar radiation, and wind speed (Albanese et al., 2012;Farmer & Wiens, 1998).
We created wetland variables from unique features, that is, nonoverlapping, with data from the Playa Lakes Joint Venture (PLJV) and the National Wetlands Inventory (NWI) ( Landscape context is also important to the state and condition of wetlands (Cariveau et al., 2011). For example, wetlands surrounded by native rangeland generally experience less sedimentation (Bartuszevige et al., 2012) and longer hydroperiods than those in cropland (Tsai, Venne, McMurry, & Smith, 2007). The inclusion of stationary land-cover variables improves the accuracy of predictions even when models are projected using future climate scenarios (Sohl, 2014); therefore, we used the cropland and grassland classes from the National Land Cover Database (Fry et al., 2011) (Table 2) (Table 2).
When comparing contemporary to future predictions, one should first project the model to contemporary GCM data (hindcast) because differences between the empirical and GCM data would otherwise contribute to the predicted changes in the probabilities of occurrence (Sofaer et al., 2016). We therefore projected models using GCM data averaged over 1981-2010 and 2041-2070, temporal ranges large enough to minimize the effects of natural variability. Specifically, we used the average annual total precipitation of each temporal range and the 30-year averages for the appropriate month.

| Modeling and evaluation
In response to the science needs of the GPLCC and the PLJV, we modeled nonbreeding distributions of two common waterfowl species and 14 shorebird species representing three families (Table 1). Collectively, these species range from short-to long-distance migrants that use a broad range of habitat types, including arid grasslands, unvegetated mudflats, and shallow-to-deep water wetlands. At each observation T A B L E 2 Resolutions and sources of the predictor variables and global climate models used to model shorebird and waterfowl probabilities of occurrence during spring and winter, respectively, across the region delineated by the Great Plains Landscape Conservation Cooperative location, we extracted a complete set of the empirical predictive variables from which we built models (Table 2). We also extracted variables across a 30 × 30 m lattice of points, including instead the GCM data, to which we projected the models.
We modeled probabilities of occurrence using Random Forests, a nonparametric method that combines the predictions from numerous trees. A tree is constructed from a bootstrap subsample of the observations and, for each branch on the tree, the data are split by selecting from a randomized subset of the predictor variables; this step minimizes modeling issues resulting from correlated variables (Breiman, 2001;Cutler et al., 2007). We were more interested in spatial predictions than inferences about the predictor variables, so reduction of correlation issues was a desirable property (see Dunn et al., 1997). Selected variables are used to hierarchically partition the data into increasing homogenous groups; thus, Random Forests inherently model variable interactions. We equalized the number of presences and absences at 80% of the number of presences for each species, specified that four randomly selected variables would be available for each branch, and built 4,500 trees for each forest.
We evaluated model performance with the out-of-bag observations, that is, observations not selected for a bootstrap sample, a standard Random Forests evaluation procedure. The class (presence or absence) of an out-of-bag observation is predicted using a majority vote across all of the trees where it was an out-of-bag observation.
The reported out-of-bag error rate is the average across all observations. Given our binary response variable, probabilities of occurrence are the proportion of trees that classified the site as a presence. For the ecoregional-level evaluations of probability of occurrence gain or loss, we used averages across the GPLCC.
We used the randomForest library (Breiman, Cutler, Liaw, & Wiener, 2014) with R version 3.1.1 (R Development Core Team, 2012) to build and test all models.

| Climate patterns
Across the GPLCC region, there is a north-south gradient in average annual temperature ranging from 5°C to 19°C  The changes predicted in spring temperature varied both spatially and across the GCMs (ranging from 2.5°C to 4.0°C) and with GFDL-CM3 predicting the highest temperatures, particularly in the western and southern portions of the region (Figure 3). An ensemble of the five GCMs revealed changes in precipitation throughout the region varying from −11 to +65 mm and changes in temperature from 2.2°C to 2.7°C (Figure 3).

| General shorebirds patterns
Models generally performed better without longitude; thus, we built models without it. Across all species, out-of-bag error rates averaged 5.70% and ranged from 1.19 to 10.51% (Table 3).  (Table 5, Appendix S1). Northward movement patterns across the spring migration window are apparent in the probabilities of occurrence of many of the species, notably Lesser Yellowlegs, Baird's Sandpiper, Least Sandpiper, and Long-billed Dowitcher (Appendix S1).
The five GCMs resulted in small differences in future probabilities of occurrence, as illustrated with predictions for the Long-billed Dowitcher in the middle of migration (Figure 4). Probabilities of occurrence were slightly larger with the two hottest GCMs, GFDL-CM3 and ACCESS1-0, in mid-April than the remaining models, and slightly smaller with the driest model, INM-CM4, than all other models. Probabilities of occurrence for all GCMs and species are presented in Appendix S2.

| Projected changes in probabilities of occurrence of shorebirds
Overall, changes in probabilities of occurrence differed among GCMs and among species relative to their spatial distribution across the region. Averaged across all locations, species, and models, probabilities T A B L E 3 Out-of-bag error rates (%) and the numbers of absences and presences correctly and incorrectly predicted. Latin names for shorebird species are provided in Table 1 Error rate    (Table 7).

| General waterfowl patterns and projected changes
T A B L E 5 General distribution pattern of shorebirds across the Great Plain Landscape Conservation Cooperative region relative to migration distance and foraging habitat Migration distance index and water depths from Skagen and Knopf (1993) and Skagen et al. (1999).

| DISCUSSION
Although we did not directly model hydrologic responses to climate in order to forecast future occurrence patterns of birds, we assumed The yellow-to-brown color ramp corresponds to small-to-large probability values. The diagrams indicate the position of GCM relative to changes in precipitation and temperature with wetter models appearing as an X above the horizontal line, warmer models as an X to the right of the vertical line, and the average model as an X in the center Latin names are provided in Table 1 2011; Sofaer et al., 2016) and, in turn, between wetland density and wetland bird abundance (Austin, 2002;Niemuth & Solberg, 2003;Steen, Skagen, & Noon, 2014 (Ambrosini et al., 2014;Skagen, Granfors, & Melcher, 2008;Skagen & Knopf, 1993;Warnock, Haig, & Oring, 1998). Such behavioral plasticity represents a natural adaptation to environmental change and contributes to the adaptive capacity of a species which allows it to cope with climate change "with minimal disruption" (Glick, Stein, & Edelson, 2011). Recent evidence reveals that long-distance migrating terrestrial birds have a greater capacity than short-distance migrants to adjust en route migration timing and trajectories (location during migration) in response to environmental variation (La Sorte & Fink, 2017). This is in part because long-distance migrants presumably experience the selective pressure imposed by greater degrees of environmental variation.
Consistent with this viewpoint, predicted increases in probabilities of occurrence of shorebirds in this study were positively related to migration distance (F 1,12 = 12.3, p = .04), suggesting greater flexibility of longer distance migrants.
The speed of migration could influence the degree to which bird species can respond to climate change (Hurlbert & Liang, 2012). Some of the predicted shifts in probability of occurrence in our study suggest that shorebirds may advance their migration calendars in response to weather. Most assessments of the effects of climate change on avian migration phenology have used data from arrival areas on the breeding grounds, and far fewer have used field data from passage areas (Gordo, 2007; but see La Sorte & Fink, 2017). Weather, and in particular temperature in the northern hemisphere, plays a large role in the speed and timing of migration (Chambers, Beaumont, & Hudson, 2014;Gordo, 2007). Although photoperiod is generally accepted as the primary trigger for onset of migration (Gwinner, 1996), environmental conditions in departure areas allow birds to fine-tune their departure dates. Furthermore, weather patterns encountered en route may influence progression speed in response to wind favorability and stopover duration in response to body condition and food supplies (Gordo, 2007). Rainfall can slow passage rates, but can also greatly influence the availability of stopover habitats in the dynamic wetlands across the Great Plains (Bartuszevige et al., 2012;Cariveau et al., 2011;Gordo, 2007;Sofaer et al., 2016). Wind, including speed and direction aloft, is also an important factor for migration (Thorup & Rabol, 2001).
Here, we deliberately specify weather, rather than climate, as it has proven to be a better predictor of vagile species (i.e., en route migrants, nomads) and is needed for modeling yearly variations. Predictions based on climate can overestimate the "availability of suitable habitat and species climatic tolerances, masking species potential vulnerability to climate change" (Reside, VanDerWal, Kutt, & Perkins, 2010).
Our empirical models based on historic weather and bird survey data captured the tight correspondence of bird occurrence and weather because of the small time frame (1 month) within which we associated birds with locations and weather variables.

(b) Baird's Sandpiper
There remains considerable uncertainty regarding the factors used as signals for speed of migration (Marra et al., 2005). Empirical models of habitat used during active migration appear to be a rare (but see Gowan & Ortega-Ortiz, 2014) and potentially difficult endeavor (Joseph & Stockwell, 2000). Not surprisingly, there has been considerably more modeling work completed on the relatively static periods of the annual cycle, that is, breeding and wintering habitat, as well as on the dates of first encounter at some geographic location, often based on data availability (e.g., Chambers et al., 2014;Murphy-Klassen, Underwood, Sealy, & Czyrnyj, 2005).
Distribution models are possibly the best available tool for predicting responses to changes in habitat and environmental conditions (but see Sinclair, White, & Newell, 2010 Migrant shorebirds and waterfowl require wetland complexes for food and protection (Farmer & Parent, 1997;Hoekman, Mills, Howerter, Devries, & Ball, 2002;Skagen, 2006), but losses of potentially usable wetlands have already exceeded 50% in the conterminous USA (Dahl, 2000;Keddy et al., 2009) due to a combination of pressures from agriculture, development, and sedimentation (Luo, Smith, Allen, & Haukos, 1997). In addition, nearly all freshwater wetland loss since the mid-1980s has occurred in the Great Plains (Dahl, 2000). Protections favor those wetlands that are relatively permanent, although the suitability of larger wetlands can be altered by the presence of smaller ones (Naugle et al., 2001). In addition, wetland complexes offer the variety of habitat preferred by some species.
Through the aggregation of newly emerging spatial tools, managers can address societal needs while basing ultimate decisions on sound science and landscape ecology foundations. The findings of our study have applicability to ongoing research and management efforts within the Great Plains Landscape Conservation Cooperative boundary. An advantage to conducting our analyses regionally, rather than from national or hemispheric perspectives, is that findings are sufficiently fine-tuned to incorporate into land management and acquisition decisions. The findings allow our partners, the Rainwater Basin Joint Venture (RWBJV) and PLJV, to evaluate their current habitat priorities and delivery actions for nonbreeding shorebirds and waterfowl and to identify new landscapes to target for current or future conservation action, thereby enhancing the adaptive capacity of target species in the face of climate change (Stein et al., 2013). These findings, when used in combination with the Decision Support Tools of the RWBJV and the landscape-design planning process undertaken by the PLJV (2014), can lead to site-and time-specific water management options that may help to counteract the detrimental effects of climate change on migrating and wintering wetland birds (Uden et al., 2015). In addition, because there is substantial uncertainty regarding future climate projections, the application of risk diversification tools such as Modern Portfolio Theory to conservation planning may serve to optimize spatial targeting of conservation actions (Ando & Mallory, 2012