Spatiotemporal variation in mechanisms driving regional‐scale population dynamics of a Threatened grassland bird

Abstract To achieve national population targets for migratory birds, landscape‐level conservation approaches are increasingly encouraged. However, knowledge of the mechanisms that drive spatiotemporal patterns in population dynamics are needed to inform scale‐variant policy development. Using hierarchical Bayesian models and variable selection, we determined by which mechanism(s), and to what extent, changes in quantity and quality of surrogate grassland habitats contributed to regional variation in population trends of an obligatory grassland bird, Bobolink (Dolichonyx oryzivorous). We used North American Breeding Bird Survey data to develop spatially explicit models of regional population trends over 25 years across 35 agricultural census divisions in Ontario, Canada. We measured the strength of evidence for effects of land‐use change on population trends over the entire study period and over five subperiods. Over the entire study period, one region (Perth) displayed strong evidence of population decline (95% CI is entirely below 0); four regions displayed strong evidence of population increase (Bruce, Simcoe, Peterborough, and Northumberland). Population trends shifted spatially among subperiods, with more extreme declines later in time (1986–1990: 28% of 35 census divisions, 1991–1995: 46%, 1996–2000: 40%, 2001–2005: 66%, 2006–2010: 82%). Important predictors of spatial patterns in Bobolink population trends over the entire study period were human development and fragmentation. However, factors inferred to drive patterns in population trends were not consistent over space and time. This result underscores that effective threat identification (both spatially and temporally) and implementation of flexible, regionally tailored policies will be critical to realize efficient conservation of Bobolink and similar at‐risk species.

Birds dependent on agricultural habitats for breeding have exhibited more significant population declines than birds of any other habitat type in North America (Herkert, 1995;Knopf, 1994;Peterjohn & Sauer, 1999). Considerable geographic variation exists in these trends (Sauer, Hines, & Fallon, 2005), and its drivers remain largely unresolved (Corace, Flaspohler, & Shartell, 2009). Given that the majority of native grassland ecosystems have been converted to agriculture (Fletcher & Koford, 2002), regional variation in agricultural practices is plausible explanations for geographic variation in avian population dynamics (Hill, Egan, Stauffer, & Diefenbach, 2014).
Correlations among changing farming practices and avian population trends are often used to identify potential drivers of population dynamics and inform the growing interest in landscape-level conservation approaches (e.g., North American Waterfowl Management Plan, Partners in Flight, The Nature Conservancy's Migratory Bird Program). Such relationships are poorly understood for many species (Chamberlain, Fuller, Bunce, Duckworth, & Shrubb, 2000), and/or best resolved at broad geographic extents (i.e., nation, state/province, or Bird Conservation Region). For example, research into the causes of state-scale patterns in grassland bird population declines implicates changes in habitat amount (Hill et al., 2014;Murphy, 2003;Perlut, 2014), pesticide use (Mineau & Whiteside, 2013), factors correlated with increasing human populations, and declines in the number of dairy farms (Perlut, 2014). While these results suggest large-scale shifts in population trajectories and underlying environmental processes, they may not reflect patterns and processes at smaller scales (Herkert, 1995). Specifically, results from broad-scale analyses could lead to incorrect inferences about threats to population persistence if (1) aggregating data at broad scales homogenizes regional variability in trend estimates (Bled, Sauer, Pardieck, Doherty, & Royle, 2013), or (2) the factors that are influential vary over time (Morris, 1994). Modeling frameworks are therefore needed that take into account spatiotemporal variability in landscape pattern and process to provide robust insights into the mechanisms driving avian population dynamics, and better inform policy development to achieve landscape-scale conservation targets (Flather & Sauer, 1996).
Bobolinks (Dolichonyx oryzivorous) are among many species of grassland birds that have declined over much of their breeding range in the past several decades, resulting in their listing as a species of conservation concern in both Canada and the United States (COSEWIC, 2010;U.S. Fish and Wildlife Service, 2008). Like other grassland birds, Bobolinks are now dependent on surrogate grassland habitats in agroecosystems for breeding (Hunter, Buehler, Canterbury, Confer, & Hamel, 2001), which implicates changes to land-use practices during this season as potential drivers of population change. Murphy (2003) identified positive correlations between Bobolink population trends and hayfield amount and negative correlations with pasture area. At the regional-scale, Bobolink population trends vary spatially and temporally (Ethier & Nudds, 2015), suggesting that variation in land management practices at this scale, rather than more broadly, may contribute to population change. We therefore tested predictions under eight alternative hypotheses about factors driving regional variation in Bobolink population trends, to determine whether, where, and to what spatiotemporal extent changes in the quantity or quality of surrogate grassland habitats contributed to population declines.
Specifically, we assessed how changes in habitat amount (hayfields and pastures, separately), hayfield composition, cattle stocking density, pesticide use, human population growth, habitat fragmentation, and factors that covary with latitude drive spatiotemporal fluctuations in Bobolink population trends (Figure 1; detailed hypotheses and predictions in Appendix S1).

| Study area
The study area covered 97,136 km 2 of the Ontario portion of North American Bird Conservation Region 13 (43°69′N 79°45′W) within the mixedwoods plains ecozone (Crins, Gray, Uhlig, & Wester, 2009 (Crins et al., 2009). The dominant agricultural sectors include fruit, row and forage crops, poultry, hogs, and beef and dairy cattle (Government of Ontario, 2013). An estimated 1.5 million hectares are surrogate agricultural grasslands (hayfield and pasture), the management of which is highly influenced by economic market forces (Agriculture and Agri-Food Canada, 2011).
These surrogate grasslands support the greatest abundance of breeding Bobolink nationally, at 10%-12% of the global breeding population (COSSARO, 2010). Bobolink are listed as a threatened species within Ontario as they have exhibited more precipitous declines here than in other provinces in Canada or in the USA (ESA, 2007). As a result, there is significant interest in understanding the factors driving regional variation in abundance to tailor effective and efficient conservation practices across the province.

| Response variable
To facilitate the use of Agricultural Census data as predictor variables, we adapted our previous analysis of Bobolink population trends (Ethier & Nudds, 2015) to align spatially and temporally with available agricultural land-use data (Statistics Canada 2011). Specifically, we summarized Bobolink trends to the extent of the agricultural census division boundaries (area range: 964-7441 km 2 ) over 25 years  and over five subperiods corresponding to census periods: 1986-1990, 1991-1995, 1996-2000, 2001-2005, and 2006-2011 (Robbins, Bystrak, & Geissler, 1986). The annual sum of Bobolinks counted over the 50 stops were used as an index of counts.
We used hierarchical Bayesian models with small-area estimation to generate spatially explicit estimates of population trends. The primary advantage of small-area estimation is that it permits more precise trend estimates for areas with limited or no data. The adjacency of neighboring areas facilitates borrowing information from one area, effectively increasing the sample size on which trend estimates are based in another (Bled et al., 2013), helping to improve the power of the model to predict trends by accounting for underlying geographic variation in the data (Bahn, O'Connor, & Krohn, 2006).
Each route was assigned to the census division in which the start point for the route was located. Counts of Bobolink ϒ i,t were modeled on route i in year t by a Poisson distribution with a mean λ i,t, which depends on the year-specific intercept (α t ), an observer effect (ω K(i,t) ) for each individual (K) carrying out the survey, and a spatial effect (b c(i),t ) at the level of the census division (c): A customary vague prior was used for the year effect parameters: Observer effects were assumed to be normal distributed random variables with a mean 0 and variance σ 2 Obs : Spatial effects b c(i),t were modeled by applying a Gaussian conditional autoregressive (CAR) model (Besag, York, & Mollié, 1991) and were allowed to vary for each subperiod separately. Neighborhood weights in a CAR model, where the spatial process is modeled over irregular regions, are generally informed deterministically using adjacency-based methods (White & Ghosh, 2009) and are widely applied in small-area estimations for landscape-scale analyses in the fields of agriculture, epidemiology, ecology, and census surveys (e.g., Bled et al., 2013;He & Sun, 2000;Lawson, Browne, & Rodeiro, 2003;Roa, 2003;Thogmartin, Sauer, & Kuntson, 2004). We assigned 1 to areas that shared a common boundary, and 0 otherwise. This neighborhood structure assumes that the spatial structure random effect for a given region has a mean equal to the average of the random effects of bordering areas and does not depend on distance between areas.
Population trends (Δ c,a,b ) could then be defined as the yearly change from year t a to year t b for census division c expressed as a percent: where n ct would represent the counts at time t of a theoretical route in cell c and is specified as: The observer variance component (0.5 × σ 2 Obs ) is the back transformed mean of the observer random effect (Smith, Hudson, Downes, & Francis, 2014).
The model was fitted by iteration using WinBUGS 1.4.3 (Lunn, Thomas, Best, & Spiegelhalter, 2000). Three chains were run to draw 50,000 samples, discarding the first 10,000 iterations. Samples were thinned by 1 in 3 to reduce autocorrelation and obtain regional mean trends (Δ c,a,b ) from the remaining posterior distributions. Convergence was checked using the Gelman-Rubin diagnostic for assessing convergence (maximum R-hat = 1.2; Brooks & Gelman, 1998) and WinBUGS traceplots. Significant posterior mean trends were interpreted as those with 95% credible intervals (CRI) that did not include zero (Bled et al., 2013). In total, 64 survey routes distributed among 35 agricultural census divisions were used to estimate Bobolink population trends F I G U R E 2 Distribution map of North American Breeding Bird survey routes (BBS) across 35 agricultural census divisions in Ontario, Canada (insert top). Posterior mean population trend estimates for Bobolink (Dolichonyx oryzivorous) over the entire study period (a) and by subperiod (b-f) are mapped to agricultural census divisions using spatially explicit hierarchical models in a Bayesian statistical framework. The indicators specify whether the 95% credible interval was above (+) or below (−) the zero mean population trend line

| Predictor variables
Direct measures (i.e., raw values) were used to reflect mechanistic hypotheses of Bobolink population change when land-use statistics were available (e.g., change in habitat amount). In instances when direct measures of mechanism were not available (e.g., cattle stocking density, hayfield composition), or hypotheses reflected aggregated mechanisms (i.e., human population growth), we derived indices which best approximated mechanistic hypotheses (Appendix S1).

| Analysis
We used Bayesian variable selection to measure the strength of evidence that environmental covariates (percent change or static measures in the case of latitude and fragmentation) for each subperiod and the whole time period were good predictors of trends in Bobolink abundance (response) (Mitchell & Beauchamp, 1988 where yi corresponds to the index of counts, β the regression coefficients of the known j-th environmental covariate, x ij is the known j-th covariate for the i-th census division, and ε is a random error term.
The embedded indicator variable I j for the j-th regression coefficient is supported at two points, 1 and 0, and is the tool for variable selection. If I j takes the value 1, the j-th regression coefficient is considered a relevant predictor of the response, and the value 0 otherwise. The regression model is fitted with a Bayesian approach, which required explicit statements of prior distribution on all unknown quantities.
We assume ε i follows a normal distribution, N(0, σ 2 ), and independent priors for β, I, and σ; where β is assigned a multivariate normal distribution; I j for j = 1,…, p are chosen independently, each with a prior Bernoulli distribution of 0.25 (Chipman, Hamada, & Wu, 1997;Meyer & Box, 1992;Meyer & Wilkinson, 1998); and σ has an inverse gamma distribution (Kuo & Mallick, 1998). The regression model was also fit with a 0.5 prior on the Bernoulli indictor (Royale & Dorazio, 2008) to determine whether the assignment of this prior affected the selection of likely important predictor variables. Results were robust to the selection of Bernoulli prior (unpublished data). We therefore retained 0.25 as our prior, which reflects the assumption that a majority of predictor variables are likely unimportant to the outcome of the response. This is a common assumption in ecological studies (e.g., Mutshinda, Finkel, & Irwin, 2013), particularly when predictor variables are derived proxies rather than direct environmental measures (Mundry, 2011).
Bayesian variable selection is generally robust to correlated predictors (Mutshinda et al., 2013), although very high collinearity (>0.90) may influence results (Ghosh & Ghattas, 2015). We therefore checked all correlation coefficients prior to analysis. All predictors were centered and scaled using the mean and standard deviation to improve the performance of MCMC and allow for comparisons among predictors (Gilks & Richardson, 1995). Variable selection was implemented in WinBUGS 1.4.3 (Lunn et al., 2000) using BugsXLA graphical user interface (Woodward, 2011). For the entire study period, and each of the temporal subperiods, we ran three chains, thinned by 1 in 3, and based our inference on 50 000 samples from the posterior distribution of parameters, after 5 000 interactions were discarded. Convergence was checked using the Gelman-Rubin diagnostic in WinBUGS (Brooks & Gelman, 1998).

| Bayesian variable selection
Pearson's correlation coefficient between predictors ranged from −0.71 to 0.51, so all variables were included in the analysis. Important predictors of spatial-temporal patterns in Bobolink population trends can be inferred from the plots of posterior probabilities of variable activity ( Figure 4). Over the entire study period , human development was positively correlated with Bobolink population trends, and to a lesser extent, fragmentation was negatively correlated with Bobolink trends (Figure 4a). By subperiod (Figure 4b-f), the likely important predictors shifted from latitude and hayfield amount in 1991-1995, to latitude in 1996-2001; none of the predictors we included were selected as likely important in 1986-1990, 2001-2005, and 2006-2011. Hayfield amount and latitude were both positively correlated with Bobolink population trends in all instances.

| DISCUSSION
Our results demonstrate that analyses performed at broad geographic extents have the potential to homogenize regional variability in population trends, which may diminish our ability to identify regions and factors responsible for avian population declines (Bled et al., 2013;Corace et al., 2009;Ethier & Nudds, 2015 (Lusk, Guthery, George, Peterson, & DeMaso, 2002;Saether & Bakke, 2000). Notably, apparent drivers of regional population change over the entire study period (human population growth and fragmentation) did not emerge as important factors in any of our temporal subperiods. This may be due to nonlinear change or nonstationarity in factors affecting trends, which are otherwise commonly not taken into account among landscapescale trend analyses (Hill et al., 2014;Mineau & Whiteside, 2006;Murphy, 2003;Perlut, 2014). Violations of assumptions of linearity and stationarity could lead to imprecise and/or inaccurate estimates of the risks that various threats pose to a population (Hill et al., 2014;Lusk et al., 2002). By allowing the slope of the relationship between response and predictors to vary by temporal subperiod, these potentially confounding effects can be partially accommodated. Results from temporally divided analyses may therefore provide additional insight into the processes that drive regional patterns in avian population dynamics at finer-temporal scales, in addition to analyses over broader time spans.
The positive relationship between Bobolink and human population trends over the entire study period is inconsistent with the hypothesis that the multiple processes by which humans could negatively affect Bobolink act at the spatialtemporal scales selected here. Conversely, Perlut (2014)  High levels of habitat fragmentation were also linked to long-term Bobolink population declines, which suggest area sensitivity and edge effects may suppress density and nest success, as it does for other grassland birds (Helzer, 1996;Johnson & Temple, 1990). These effects may accumulate with time, since a Bobolink's decision to return to, and breed at, a site is influenced by its own reproductive success, and that of others, in the previous year (Bollinger & Gavin, 1989). Thus, while our measure of fragmentation is static, the demographic processes it influences (i.e., apparent survival) could be cumulative, creating the negative trends we observe. If a temporally varying measure of landscape fragmentation were available, we would expect this effects to be exacerbated. The long-term viability of Bobolink may therefore depend on the maintenance of large tracts of surrogate grassland habitat where they currently exist, and/or restoring open habitat types in regions where they have been lost.
When our analysis was temporally segregated over five subperiods, Bobolink population trends varied widely (Figure 2b-f), similar to patterns observed for Cerulean warbler (Setophaga cerulean) and Red-bellied woodpecker (Melanerpes carolinus) over similar spatial and temporal scales (Bled et al., 2013). Factors associated with these shifts was with hayfield amount. Given that Bobolinks are upwards of four times more likely to be found in hayfields than any other surrogate grassland habitat type (Bollinger, Bollinger, & Gavin, 1990), the maintenance of hayfields on the landscape, regardless of composition, will likely be important to Bobolink conservation.
Bobolink population growth rates were positively correlated with latitude during the 1990s, consistent with the hypothesis that climatedriven gradients in plant maturation (hence, forage harvest timing) may contribute to spatial patterns in Bobolink population trends. In an average year, optimal hay harvest timing in Ontario is more than 3 weeks earlier in the southwest than in the northeast (J. J. Nocera, pers. commun.), similar to patterns found in Michigan (Corace et al., 2009 Our derived indicator of hayfield composition was slightly negative over the entire study period, in contrast to results of some other studies (McCracken, 2005;Mussel, Schmidt, Ethier, & Doug, 2013). One reason for this discrepancy may be that we corrected our indicator for area. Previous studies that assessed hay harvest timing and frequency also failed to establish relationships between hayfield composition and Bobolink population trends (Corace et al., 2009;Herkert, 1997), suggesting that our result is not an artifact of how the predictor variable was derived.
Results presented here corroborate those of Hill et al. (2014) and Perlut (2014) in that overall pesticide use on the breeding grounds did not appear to be an important predictor of Bobolink population declines (but see Mineau & Whiteside, 2013). This could be due to, for examples, an overall decrease in the use of pesticides, or the general reduction in the toxicity of pesticides in the past several decades (Mineau & Whiteside, 2006). Future efforts to collect and model spatially and temporally explicit data on the amount, type, and application rates of various pesticides used per crop types, and their lethal, sublethal, and indirect effects on food resources and refuges, would greatly improve our understanding of pesticide risk to avian populations.
Our models primarily addressed how the loss and degradation of grassland habitat during the breeding season may affect Bobolink population trends, as these factors are within the management jurisdictions of local governments. However, Bobolinks live the majority of the year away from the breeding grounds, such that changes to migratory and wintering habitats may have important effects on population dynamics (Webster, Marra, Haig, Bensch, & Holmes, 2002). In systems with strong migratory connectivity, habitat degradation on the nonbreeding grounds should result in high spatial variation in the breeding season population trends. Connectivity is weak between breeding and nonbreeding areas for Bobolinks; Bobolinks from across the breeding range synchronously congregate on the nonbreeding ground over a relatively small geographic range (Renfrew et al., 2013). Thus, if the important drivers of Bobolink population change operate on nonbreeding areas, we might anticipate some degree of spatial trend homogenization (as occurred only during 2006-2011). Thus, the strong spatial variation in population trends we observed is inconsistent with the hypothesis that nonbreeding ground effects are driving Bobolink population trends in Ontario. Local conditions in the breeding range not accounted for among our predictor variables may more likely explain the observed spatial structure of trends (Bled et al., 2013).
We chose to use an irregular lattice for our CAR model because it enabled calculation of Bobolink population trends for agricultural census divisions, the spatial units at which land-use statistics were available for testing our alternative hypotheses. However, different neighborhood structures are available, the selection of which can be critical for spatially explicit parameter estimation (Thogmartin et al., 2004). We assumed that the spatial random effect depended only on the adjacent neighborhood structure, and not on the distance between areas. Thus, the neighborhood structure used here may not fully reflect underlying environmental autocorrelations. Alternative spatial models include regular grids (Bled et al., 2013), geostatistical interpolation (Diggle, Tawn, & Moyeed, 1998), continuous surface models (Kelsall & Wakefield, 2002), or the replacement of spatial autocorrelation structure with appropriate environmental covariates at a scale relevant to underlying biological processes (Thogmartin et al., 2004). Further, for CAR models, the standard error of the spatial effect can have lower precision for cells with fewer neighbors (i.e., those at the periphery of our study area; Law & Haining, 2004). As a result, better-connected census divisions are more likely to be significant due to narrower CRI. Using a logistic regression, we did a post hoc assessment to determine whether significant trend estimates were found more often than expected in better-connected neighborhoods. We did not find evidence to support this assertion (p = .88, unpublished data), suggesting our results are not influenced by variation in the standard error of the spatial effect.
With a growing interest in developing regionally tailored management plans for priority grassland species (Cooper, 2007), consideration ought to be given to finer-scale mechanisms to develop effective and efficient policy options with potential to arrest or reverse declines.
Regional variability in Bobolink population trends, and the predictors identified to drive them, suggests that a one-size-fits-all approach to grassland bird conservation may not adequately address threats to populations. Our results therefore speak to the need for active adaptive management to accommodate complex environmental systems that are not constant in space and time (Walters & Holling, 1990). This framework would allow resource managers to incorporate and address uncertainty regarding the factors that account for population declines, including nonlinearity and nonstationarity, while reducing these uncertainties by deliberate experimentation through policy implementation.

ACKNOWLEDGMENTS
We thank J.B. Sayers and two anonymous reviewers for comments and edits on earlier drafts of this manuscript. We also acknowledge the countless volunteers involved in collecting Breeding Bird Survey data.