Opposing responses to drought shape spatial population dynamics of declining grassland birds

The joint threats of climate and land‐use change require an understanding of how environmental variation influences species abundance and distribution. However, most species distribution models use static data and methods without considering how species respond over multiple temporal and spatial scales. Using a novel analytical approach, we show how multiscalar environmental variation drives spatial population dynamics of mobile species.


| INTRODUC TI ON
Rapid global environmental change is affecting biological systems leading to the redistribution of species, shifts in community composition, altered phenology and the disruption of ecological interactions (Scheffers et al., 2016;Thackeray et al., 2016). Despite the attribution of climate change effects across broad taxonomic groups, fundamental challenges persist on how to improve our understanding and predictions of species responses to climate change (Ehrlén & Morris, 2015;Urban et al., 2016). Recent work has emphasized the need to move away from static modelling of species distributions in favour of integrating key dynamic ecological processes reflective of population dynamics (e.g., dispersal, migration and demographic indices) and to consider the multiscalar influence of environmental factors (e.g., short-and long-term climate and lag effects) that affect where species occur and persist on landscapes (Franklin, 2010;Guisan & Thuiller, 2005;Maclean, Suggitt, Wilson, Duffy, & Bennie, 2017).
Many mobile species select habitats in temporally dynamic and complex ways with movement typologies ranging from predictable migratory movements between seasons to irruptions and nomadic behaviours occurring within seasons (Newton, 2012;Somveille, Rodrigues, & Manica, 2015). Despite the varied mechanisms involved in these movement behaviours, they are unified by an ability to respond and adjust to climatic conditions and fluctuating resources over time-scales within and between years (Bateman, VanDerWal, & Johnson, 2012;Runge, Tulloch, Hammill, Possingham, & Fuller, 2015;Wilson, Anderson, Wilson, Bertram, & Arcese, 2013). To capture the complexity of mobile species distributions and to accurately assess the risk of decline under threats such as climate change, previous work has highlighted the need for dynamic modelling (spatial and temporal) of species distribution information (Franklin, 2010;Naujokaitis-Lewis & Fortin, 2016;Runge et al., 2015). Most research has focused on the influence of climate on species distributions across different spatial scales, but temporal dimensions of climate (e.g., resolution, extent and lags) are increasingly noted as critical elements that govern the ability to detect species responses (Adrian, Gerten, Huber, Wagner, & Schmidt, 2012;Bateman et al., 2016;Chan et al., 2016). Prior research suggests that there is no one-size-fits-all when considering temporal scale of variables as biologically relevant scales will be informed by the characteristics of the species and the dynamics of the ecosystem (Briga & Verhulst, 2015;van de Pol & Cockburn, 2011;van de Pol et al., 2016). In many cases, species may respond immediately to short-term weather events; however, lags between the onset of an environmental cue and ecological response constitute an important and often overlooked temporal dimension (van de Pol & Cockburn, 2011;van de Pol et al., 2016).
Climate change is predicted to have strong effects on the intensity and frequency of drought across geographical regions and ecosystems (Dai, 2011), and therefore, understanding the extent and mechanism by which mobile species respond to drought is key for predicting future distributions and population dynamics. As a climate variable, drought is intrinsically multiscalar and complex (Vicente-Serrano, Beguería, & Moreno, 2011) and indices capturing drought over multiple time-lags could reveal different effects on mobile species due to alternative mechanisms of action. For example, individuals could respond positively to periods of wetness via mechanisms such as higher prey availability and increased survival rates (Skagen & Adams, 2012;Yackel Adams, Skagen, & Savidge, 2006), while delayed effects of longer-term drought might influence vegetation resulting in a change in the availability of suitable breeding sites (Boulangeat et al., 2014;Vicente-Serrano et al., 2013).
Despite increasing acknowledgement of the importance of selecting an appropriate climate window for understanding species response to changing environmental conditions, few studies comprehensively consider response to climatic variation across temporal and spatial scales.
We examined the spatial population dynamics and the role of short-and long-term climatic conditions for two threatened grassland birds (Lark Bunting, Calamospiza melanocorys, Chestnutcollared Longspur, Calcarius ornatus) that breed in the Great Plains of North America. The Great Plains are a dynamic ecosystem that evolved under periodic drought (Clark et al., 2002;Woodhouse, Lukas, & Brown, 2002) combined with frequent disturbance due to mammalian grazing (e.g., Bison Bison bison, Prairie Dogs Cynomys spp) and fire (Askins et al., 2007). Spatial variation in precipitation and the intensity of disturbance creates a mosaic of vegetation communities with different faunal species favouring particular conditions (Askins et al., 2007;Augustine, 2010;Samson, Knopf, & Ostlie, 2004). As such, grassland endemics often have a spatially and temporally patchy distribution with vagility being a frequent characteristic (Bateman et al., 2015;Shane, 2000). This tendency for movement in response to environmental variation now occurs in the context of extensive habitat losses that occurred over the past two centuries and continue today (Brennan & Kuvlesky, 2005;Samson et al., 2004).
To estimate range-wide and regional abundance of the two species, we used a spatial analytical approach applied to 48 years  of citizen science data (North American Breeding Bird Survey (BBS), Environment Canada, 2017, Sauer et al., 2017. We examined long-term spatial trends and asked whether regional abundance was influenced by drought conditions at two temporal scales, the month of arrival on the breeding grounds and over the preceding 4 years. Modelling drought at these two scales represents mechanisms of short-term resource availability (e.g., wet springs lead to higher insect abundance) and long-term vegetation structure (e.g., successive wet years change the structure of prairie vegetation). We also tested whether environmental conditions experienced along the migration pathway influenced range dynamics. Recent studies have revealed how latitudinal migrants track shifting vegetation phenology during the non-breeding period (Thorup et al., 2017). Here, we show evidence for spatially autocorrelated abundance across the range of a latitudinal migrant due to the influence of environmental conditions along the migration pathway on breeding settlement patterns.

| Study species
Lark Bunting and Chestnut-collared Longspur are endemic to the North American Great Plains and are obligate migrants with distinct breeding and non-breeding areas (Bleho, Ellison, Hill, & Gould, 2015;Shane, 2000). The breeding range of the two species includes the Canadian prairie provinces and the northern Great Plains states with the Lark Bunting range extending further south into northern New Mexico and Texas (see Supporting information Figure S1). Both species overwinter in grasslands of the southern US states and northern Mexico. Breeding habitat for Lark Buntings includes short-grass prairie and sagebrush steppe, but they also use cultivated areas, such as hay fields, to a lesser extent (Shane, 2000). Chestnut-collared Longspur breeding habitat consists of arid short-to mixed-grass prairie with a preference for recently disturbed sites after grazing or burns.

| Population data
We analysed 48 years of data  from the BBS to examine spatial population dynamics, distribution shifts and the influence of drought on the regional abundance of our two focal species. The BBS is conducted annually during the breeding period in late May to early July. The survey is based on ~40-km roadside routes with a 3-minute point count conducted at stops every 0.8 km for a total of 50 counts per route. The number of individuals for a species is typically summed across all 50 counts to provide a single abundance per route as the sampling unit. We used the route-level counts to estimate the annual abundance for each species within geographical strata defined by terrestrial ecoregions (Wiken, Francisco, & Griffith, 2011). The larger ecoregions were further split by lines of latitude and longitude to allow for greater spatial variation in parameters (e.g., across the Canada-USA border), while retaining a sufficient number of BBS routes in each stratum to generate reasonably precise estimates (see Supporting information Figure S1). For each species, geographical strata were included if they contained at least five routes on which the species was observed and at least one route on which the species had been observed in >10 years.

| Environmental covariates
We used the standardized precipitation evapotranspiration index (SPEI) to test the effect of different types of drought on population dynamics (Beguería, Vicente-Serrano, Reig, & Latorre, 2014). The SPEI quantifies deviations from the monthly climatic water balance, based on the differences in precipitation and potential evapotranspiration (PET). The SPEI is a site-specific drought indicator that incorporates temperature and PET into the calculation, which improves on previous drought metrics that use precipitation alone Spring migration generally attenuates by the end of May (Bleho et al., 2015;Shane, 2000), and therefore, this month represents conditions individuals experience at the start of the breeding season including settlement and nest site selection. However, we acknowledge that conditions later in the breeding season might also influence movement (e.g., Skagen, Augustine, & Derner, 2018).
SPEI calculations used the Penman-Monteith method for the estimation of PET and were accessed from SPEIbase v.2.4 at a resolution of 0.5 degrees and were obtained from http://digital.csic. es/handle/10261/128892. To reflect the scale of our abundance models, SPEI variables were calculated as annual averages within each of the 14 geographical strata (Supporting information Figure   S1, Table S1).

| Statistical analysis
The BBS count data were modelled as hierarchical overdispersed Poisson variables and models were fit using Markov chain Monte Carlo methods in JAGS (Plummer, 2003) implemented through R with package rjags (Plummer, 2015;R Development Core Team, 2004). Hierarchical Bayesian approaches are well suited for population time series analyses because they provide a robust framework to account for sources of variation at multiple scales including observer effects, overdispersion and the spatially specific effects of drought across the species range. For the two species, we modelled each count C i,j,t (for stratum i, observer-route combination j and year t) as a Poisson random variable with mean λ i,j,t defined as a log-linear relation to the predictor variables, which include process and sampling components: The model includes stratum-specific estimates for average abundance per route (α i ) and log-linear population trends (β 1 i ). β 2 i represents the effect of stratum-specific drought (SPEI 1 i,t ) in the preceding month, and β 3,i represents the effect of stratum-specific drought (SPEI 48 i,t ) in the preceding 48 months. For strata in the northern part of the species range, β 4,i represents the effect of drought conditions in the preceding month in the southern portion of the species range (SPEI 1S t ; I 1 (i) = 1 if stratum i is in the northern part of the species range and 0 otherwise). All β . i parameters were estimated as random effects drawn from normal distributions with ). We used uninformative normal priors with mean 0 and variance 10 2 for B 1 , B 2 , B 3 , and B 4 and inverse-gamma priors with shape and scale parameters equal to 0.001 for all variances. We also modelled variation in the counts associated with routes and observers (Link & Sauer, 2007); observer-route combinations (ω j ) were treated as normal random variables with a mean of 0. BBS observers also differ in their experience at identifying species, with variation being particularly high between their first year and all subsequent years (Link & Sauer, 2007). We estimated these start-up effects as I 2 (j,t) , where I 2 = 1 when year t was an observer's first year on route j and 0 otherwise, and η was assigned a diffuse, normal prior with mean 0 and variance 10 6 . The component ε i,j,t is an observation-level random variable and helps accommodate for overdispersion, which is common in count data. Variance associated with observer-route combinations ( 2 ω ) and overdispersion ( 2 ε ) were assumed to be constant across strata and were assigned vague inverse-gamma prior distributions with shape and scale parameters = 0.001.
We ran three chains, using the first 500 iterations as an adaptive phase; the next 100,000 iterations were discarded as a burn-in, and we retained an additional 50,000 iterations, thinned by a factor of 5 to give 10,000 samples from the posterior distribution for inference.
We assessed model convergence through the parameter history plots and R-hat convergence diagnostics (Plummer, Best, Cowles, & Vines, 2006). We interpreted the strata population trends and effects of drought for each species by examining the posterior median and associated 95% credible intervals (CI) for the β . parameters and the overall effects of the drought measures using the posterior distribution of the hyperparameters B 2 , B 3 , and B 4 .
We calculated the annual distribution centroid for each species from the centroid coordinates of each stratum weighted by the estimated annual abundance of the species in each stratum (Huang, Sauer, Swatantran, & Dubayah, 2015): where z i is the proportion of routes within stratum i on which the species has been observed. We mapped the annual centroids for both species to examine the annual variation in the population

| Spatial variation in abundance and temporal trends
Both species showed strong spatial variation in abundance and trends across their ranges (Figure 1

| Distribution shift over time
The

| Influence of short-and long-term drought on distribution
The ecoregion abundance of Lark Buntings showed contrasting responses to short-term (1-month) versus long-term (4-year) drought conditions ( Figure 3). There was an overall positive response to the Model coefficients for SPEI 1-month and SPEI 48-month broadly overlapped 0 with a negative response to the latter in only one ecoregion, NmgU ( Figure 3).
We hypothesized that wetter short-term conditions might influence the species distribution if individuals choose to settle in favourable areas as they pass through them during northward migration.
We found strong support for this hypothesis for Lark Buntings based on 1) strong, negative relationships between SPEI 1-month in the

| D ISCUSS I ON
Many species globally face the joint threat of land conversion combined with rapidly changing environments and ecosystems due to climate change (Mantyka-Pringle, Martin, & Rhodes, 2011;Travis, 2003). Amidst these threats, it is imperative to understand the patterns and mechanisms underlying spatial and temporal fluctuations in abundance and distribution of species. While species distribution models often employ a static approach, distributions are expected to fluctuate to varying extents as individuals move in accordance with environmental conditions (e.g., Bateman et al., 2012;Runge et al., 2015;Wilson et al., 2013). Using a novel analytical approach that allowed us to estimate annual ecoregion-level abundance over a 45-year period, we were able to examine how drought at different temporal scales influenced annual abundance across species ranges.
Our results revealed several important findings. First, we showed that species may respond both positively and negatively to the same environmental variable across time periods suggesting a temporal hierarchy in the mechanisms determining range-wide abundance. Our work provides strong support for the importance of incorporating greater temporal resolution and time-lags when examining how environmental variation affects species distributions (Adrian et al., 2012;Nadeau, Urban, & Bridle, 2017 The different responses to drought at the two scales almost certainly reflect two mechanisms by which environmental conditions influence movement and settlement behaviours. Our intention was not to identify the mechanisms influencing these patterns, and such an effort would be difficult at scales covering the distribution of a species. However, we propose that a short-term positive response to moisture likely reflects a mechanism related to resource availability (e.g., Skagen & Adams, 2012), while a long-term positive response to drier conditions may be due to the influence of drought on habitat structure (e.g., Vicente-Serrano et al., 2013). The latter response for Lark Buntings was most pronounced along the eastern periphery of the range where wetter periods may lead to a taller grassland habitat structure outside of what is typically preferred by this species (Niemuth, Solberg, & Shaffer, 2008;Shane, 2000). While core western areas likely provide a buffer when eastern and northern areas of the range are less suitable, these peripheral regions could allow the population to expand in some years when conditions are favourable (Bateman et al., 2012(Bateman et al., , 2015Magoulik & Kobza, 2003). Our results also revealed a long-term westward shift of the Lark Bunting distribution amidst the considerable interannual variability along a north-south axis (see further detail below). Changes in precipitation may have influenced this westward range shift, but the considerable conversion of grassland to cropland along the eastern edge of the range likely also had an effect (Wright & Wimberly, 2013).
Sympatric species may vary in the extent to which they respond to variation in environmental conditions (Gorzo et al., 2016;Greenville, Wardle, Nguyen, & Dickman, 2016;Martin et al., 2017;Skagen et al., 2018). In contrast to the strong response to environmental conditions displayed by Lark Buntings, we found no evidence that Chestnut-collared Longspur spatial dynamics were influenced by drought conditions at either temporal scale. Although it is possible that both species were influenced by other unmeasured climatic variables, our drought index incorporated temperature and precipitation in its derivation, and thus, we did not include these in our analyses. With little response to environmental variation, there was also little annual variability in the longspur distribution centroid. Historically, different types and extents of movement behaviours among species may have been equally suitable strategies to cope with the normal range of variability in climate experienced in the Great Plains. However, these strategies are now placed in the F I G U R E 4 Relationship between the latitude coordinate of Lark Bunting annual population centroids and spring moisture (SPEI 1 month -southern range) in the southern part of the species range. Points are scaled to reflect the relative precision of the estimated latitude of the population centroid in each year with larger points indicating years with more precise latitudinal estimates context of heavy losses of breeding and non-breeding grassland habitat (Askins et al., 2007;Pool, Panjabi, & Macías-Duerte, 2014) as well as climate change, where model predictions indicate an increase in drought severity over the next few decades (Cook, Ault, & Smerdon, 2015). Theoretical and empirical research suggests that more specialized and less mobile species may be most susceptible to the joint impacts of land-use and climate change (Travis, 2003;Warren et al., 2001).
However, our study highlighted another pattern that is much more rarely observed where the breeding distribution of an obligate directional migrant shifts annually along a latitudinal axis depending on the conditions experienced by individuals during migration as they first reach the southern portion of the breeding range. The more extensive settlement in these southern regions then leads to a decline in abundance in northern regions highlighting the influence of environmental conditions in driving spatially autocorrelated abundance across the range. This type of behaviour has been observed in other avian taxa including water birds and raptors but primarily with a shortening of the autumn migratory pathway and a northward shift of the winter distribution relative to the breeding range (reviewed in Elmberg, Hessel, Fox, & Dalby, 2014). There is little known on the extent to which environmental conditions drive these spatially autocorrelated patterns of abundance in migratory passerines, but these relationships may be underestimated because of the challenges in tracking small songbirds over large spatial scales. In one study in eastern North America, Rushing, Dudash, Studds, and Marra (2015) used isotopic markers to show that in years of earlier spring phenology, American Redstarts (Setophaga ruticilla) appear to shift their breeding distribution northward likely because optimal conditions to initiate breeding in the southern portion of the range had already passed by the time migrating individuals reached those locations.
There is some evidence suggesting that the breeding distributions of other prairie endemics may shift in response to precipitation (Green, Lowther, Jones, Davis, & Dale, 2002 most strongly in the north, they are declining consistently across the range. These declines are likely due primarily to the conversion of grassland to cropland on both the breeding (Askins et al., 2007;Vickery et al., 1999) and non-breeding grounds (Pool et al., 2014). An interesting result from the stronger declines in northern regions of the range is that the distribution centroid for Chestnutcollared Longspurs has shifted to the south. This southward shift is opposite to common expectations that species ranges will shift poleward with warmer temperatures (La Sorte & Jetz, 2012;Root et al., 2003), although relationships between climate change and range shifts are complex and species may not respond as expected (Coristine & Kerr, 2015;Currie & Venn, 2017 (Derner, Lauenroth, Stapp, & Augustine, 2009;Knopf, 1996). Our results highlight how variation in species responses to climate could further influence which sites are a priority at any one time. Efforts to meet conservation targets for such mobile species could incorporate dynamic conservation approaches (Reynolds et al., 2017) where the priority areas and strategies employed vary in time and space to more efficiently benefit such species.

ACK N OWLED G EM ENTS
We appreciate the contributions of all the volunteers who conducted Breeding Bird Survey routes and the work of C. Robbins and the United States Geological Survey for creating and maintaining the survey. We also thank two anonymous reviewers for valuable comments on an earlier version of this manuscript.

DATA ACCE SS I B I LIT Y
Avian count data used in this manuscript are publicly available from