Integrating broad‐scale data to assess demographic and climatic contributions to population change in a declining songbird

Abstract Climate variation and trends affect species distribution and abundance across large spatial extents. However, most studies that predict species response to climate are implemented at small spatial scales or are based on occurrence‐environment relationships that lack mechanistic detail. Here, we develop an integrated population model (IPM) for multi‐site count and capture‐recapture data for a declining migratory songbird, Wilson's warbler (Cardellina pusilla), in three genetically distinct breeding populations in western North America. We include climate covariates of vital rates, including spring temperatures on the breeding grounds, drought on the wintering range in northwest Mexico, and wind conditions during spring migration. Spring temperatures were positively related to productivity in Sierra Nevada and Pacific Northwest genetic groups, and annual changes in productivity were important predictors of changes in growth rate in these populations. Drought condition on the wintering grounds was a strong predictor of adult survival for coastal California and Sierra Nevada populations; however, adult survival played a relatively minor role in explaining annual variation in population change. A latent parameter representing a mixture of first‐year survival and immigration was the largest contributor to variation in population change; however, this parameter was estimated imprecisely, and its importance likely reflects, in part, differences in spatio‐temporal distribution of samples between count and capture‐recapture data sets. Our modeling approach represents a novel and flexible framework for linking broad‐scale multi‐site monitoring data sets. Our results highlight both the potential of the approach for extension to additional species and systems, as well as needs for additional data and/or model development.


| INTRODUC TI ON
Widespread population declines, range retractions, and extinctions highlight an urgent need to better understand drivers of wildlife population dynamics (Ceballos, Ehrlich, & Dirzo, 2017;Tittensor et al., 2014). Climate variation and trends can play a crucial role in determining population trajectories (Stephens et al., 2016). Most studies that relate climate to populations across large spatial extents are based on occurrence or count data (Dawson, Jackson, House, Prentice, & Mace, 2011;Pacifici et al., 2015). These data types have become relatively common and available across large spatial extents (Sullivan et al., 2009) and can be used to model both population state and demographic rate parameters (Dail & Madsen, 2011;Royle, 2004). However, estimates of demographic rates from such models may not always be reliable (Barker, Schofield, Link, & Sauer, 2018;Dennis, Morgan, & Ridout, 2015;Zipkin et al., 2014). Capturemark-recapture (CMR) data provide additional information for modeling demographic rates, providing a more mechanistic link between population dynamics and climate than models based on occurrence or count data alone (Amburgey et al., 2018;Buckley et al., 2010;McMahon et al., 2011;Selwood, McGeoch, & Mac Nally, 2015).
However, CMR data are less available and are relatively costly to obtain across large spatial extents. Analyses that incorporate strengths of different data types and models have the potential to improve the inferences about population dynamics.
Integrated population models (IPMs) provide a formal framework for jointly modeling independent count and CMR data (Besbeas, Freeman, Morgan, & Catchpole, 2002;Hostetler, Sillett, & Marra, 2015;Schaub & Abadi, 2011) along with climate predictors of demographic rates (Zipkin & Saunders, 2018). With one recent exception (Zhao, Boomer, & Royle, 2019), IPMs that have included climate covariates of demographic rates have been limited to population studies across relatively small spatial extents (Woodworth, Wheelwright, Newman, Schaub, & Norris, 2017) or have modeled multiple local populations independently (Weegman, Arnold, Dawson, Winkler, & Clark, 2017). Data from long-running national and continental scale avian monitoring programs (Dunn et al., 2005;Gregory et al., 2005;Robinson, Julliard, & Saracco, 2009;) present a unique opportunity for extending IPMs to broad-scale applications. However, development and application of IPMs for broadscale bird-monitoring data are challenging because of a variety of sampling issues, including mismatches in sizes and spatio-temporal distribution of sampling areas between monitoring data sets and properly accounting for observation error associated with multi-site, multi-observer studies (Ahrestani, Saracco, Sauer, Royle, & Pardieck, 2017;Robinson, Morrison, & Baillie, 2014). Additionally, IPMs developed for local scale studies based on binomial and Poisson population processes may not be appropriate at regional scales where population responses represent averages of local studies (Zhao et al., 2019).
Analyses of broad-scale IPMs require spatial stratification at ecologically relevant scales. Stratification decisions may involve geopolitical boundaries that are meaningful in conservation applications (e.g., state × bird conservation region; . One approach would be to stratify on as fine of a resolution as possible and model spatial structure explicitly (Bled, Sauer, Pardieck, Doherty, & Royle, 2013;Saracco, Royle, DeSante, & Gardner, 2010). This approach offers the advantage of allowing poststratification summaries across any larger spatial resolution of interest but can be computationally expensive. Alternatively, stratification may be based on natural spatial structuring of population dynamics (Rushing, Ryder, Scarpignato, Saracco, & Marra, 2016) or genetics (Ruegg, Harrigan, Saracco, Smith, & Taylor, unpublished data;Ruegg et al., 2014).
Spatial structuring of migratory species may be maintained (strong migratory connectivity) or dissolved (weak connectivity) between breeding and nonbreeding seasons. Effectively linking environmental covariates to demography in these species requires understanding of spatial structuring throughout the annual cycle (Cohen et al., 2018;Ruegg et al., 2014) and consideration of potential "carry-over effects," whereby environmental conditions experienced in one season affect demographic rates and population changes in a subsequent season (Norris & Marra, 2007 (Ruegg et al., 2014). Wilson's warbler ( Figure 1) is a good candidate for development of a climate-informed IPM, as it is well-represented in BBS and MAPS data sets in the western United States; its patterns of migratory connectivity are relatively well understood (Ruegg et al., 2014); its population has declined in recent decades ; and it has been designated as an "at-risk" species due to climate change (http://clima te.audub on.org/birds ). We included remote sensed and modeled climate covariates in models of demographic rates. We hypothesized that F I G U R E 1 Adult male Wilson's warbler in California's Sierra Nevada. Photography credit: Gabriel Gonzalez breeding productivity would depend on drought conditions on the wintering grounds (a carry-over effect) and on spring temperature (Saracco, Siegel, Helton, Stock, & DeSante, 2019;Socolar, Epanchin, Beissinger, & Tingley, 2017) and that adult survival would depend on winter drought and wind conditions during spring migration (Drake, Rock, Quinlan, Martin, & Green, 2014;Huang, Bishop, McKibbin, Drake, & Green, 2017;LaManna, George, Saracco, Nott, & DeSante, 2012 to decompose variation in population growth rates among vital rate and demographic structure components and to examine how these demographic contributions depended on climate covariates (Koons, Arnold, & Schaub, 2017;Koons, Iles, Schaub, & Caswell, 2016).

| Bird-monitoring data and focal species
Our analysis incorporates annual counts of adult birds from the BBS and capture-recapture data on adult birds and age-specific capture data from the MAPS program. The BBS is a roadside bird survey established in 1966 that provides data on the status and population trends of >420 bird species Sauer et al., 2017); it is a core component of North American bird conservation efforts (Rosenberg et al., 2016). The MAPS program, established in 1989 andstandardized in 1992, uses data from a cooperative network of mist-netting and bird-banding stations to provide information on demographic rates of >100 landbird species (DeSante & Kaschube, 2009). Here, we analyze MAPS and BBS data for Wilson's warbler stratified by three genetically distinct breeding regions that overwinter in northwestern Mexico (Ruegg et al., 2014; Figure 2; Table 1). Although MAPS and BBS data were available from additional genetic breeding regions, we limited our analysis to just these three regions because other breeding populations regularly utilize wintering ranges farther south than our climate covariate data set extended (Wang, Hamann, Spittlehouse, & Carroll, 2016;Wang, Hamann, Spittlehouse, & Murdock, 2012). For illustration of our model, we limit the time window of our analysis to the 17 years 1992-2008 based on the earliest year and latest years of vetted MAPS data available when the analysis was undertaken (updated verified MAPS data base expected in 2020).

| Climate data
We calculated overwintering and breeding season climate covariates from the ClimateNA database (https ://sites.ualbe rta.ca/~ahama nn/data/clima tena.html). ClimateNA uses bilinear interpolation of monthly gridded climate data (Daly et al., 2008;Hutchinson, 1989) and local elevation adjustments to provide climate metrics for individual points (Wang et al., 2016(Wang et al., , 2012). An advantage of using this data set is the availability of interpolated values for fine spatial resolution and projected climate estimates for future time periods to assess population viability under climate change scenarios.

| Overwintering season
Drought associated with the dry season may limit vital rates of birds overwintering in western Mexico (LaManna et al., 2012;Nott, Desante, Siegel, & Pyle, 2002); and winter drought conditions in this region are expected to become more severe in the coming decades (IPCC, 2014). To characterize annual winter drought conditions, we used Hargreave's Climate Moisture Deficit (cmd), a derived variable calculated as the monthly summed difference between atmospheric evaporative demand (Hargreaves & Samani, 1982) and precipitation.
We extracted winter (December-February) cmd values for 10,000 random points across the winter range for 1992-2008 and for "normal" (i.e., mean) values for . We then calculated the deviation of the annual value from the normal value to derive a cmd anomaly (cmd) for inclusion in the analysis. For birds that breed in the coastal California region (cca), we used cmd values from both the F I G U R E 2 Breeding Bird Survey (BBS) routes (squares) and MAPS stations (circles) sampled between 1992 and 2008 where Wilson's warbler was detected or captured. The three genetically distinct breeding regions (Ruegg et al., 2014) included the Pacific Northwest (pnw; purple), coastal California (cca; blue), and the Sierra Nevada (sne; orange). Birds of all three breeding regions winter in northwest Mexico (blue), although migratory connectivity data suggest that only the cca breeding region includes the southern Baja California portion of the wintering range. Points for which spring migration wind data were used are shown for each breeding region (purple circles for pne; green ×'s for cca; orange +'s for sne) Baja California and western Mexico portions of the overwintering regions, based on evidence that birds from this breeding region use both of these areas nearly equally during the nonbreeding season (Ruegg et al., unpublished data). For the remaining two breeding regions, we used only cmd values from western Mexico, as birds from these two regions do not appear to use (or use minimally) the Baja California portion of the overwintering range.

| Breeding season
We calculated mean spring (February-May) temperature (temp) within each of the breeding regions from averaged ClimateNA values from 1,000 random points sampled within each breeding region. Temperature covariates were among the top explanatory variables in models of projected changes in the distribution of Wilson's warbler under climate change (C. Wilsey and N. Michel, pers. comm.), and spring temperature in particular may affect snowmelt, green-up, and food availability during the nesting season, which can affect timing of breeding and productivity (Saracco et al., 2019). As with the drought covariate, we subtracted the annual spring temperature values from 1961 to 1990 normal values to derive the covariate used in the analysis described below.

| Spring migration
Spring wind conditions may be an important factor affecting survival rates of western Neotropical migrant songbirds (Drake et al., TA B L E 1 Summary of the BBS and MAPS data sets used in the integrated population model  (Kalnay et al., 1996). We averaged 6

| Model development and implementation
Our overall model was comprised of three basic components

| Count model for BBS data
At the core of the IPM is a state-space model for the count data (solid blue box in Figure 3). We modeled the counts, y, from i in 1, …, Nro We modeled n m,t at each time step as the sum of surviving adults from the previous year, s m,t and the number of recruits (local recruits + immigrants) entering the population between years, g m,t . For Interpretation of m,t should also be cautioned by acknowledgement that this parameter may also absorb unexplained variation representing discrepancies in the sampling process between data sets. We decomposed m,t into a productivity component derived from age-specific MAPS data, RI m,t , and a latent parameter, ι m,t , which reflects variation in first-year survival and immigration (i.e., local and external "recruits" of any age; DeSante, Kaschube, & Saracco, 2015), based on m,t = RI m,t × m,t . We derived regional-scale route-level abundance indices as: where the w m are weights representing the proportion of BBS routes on which Wilson's warblers were encountered in the region and the 2 i and 2 i,m,t are variance components of route × observer and overdispersion effects . We calculated regional population trends as geometric means of the annual realized population growth rates (N m,t+1 N m,t ). For composite abundance and trend estimates, we weighted regional abundances by proportions of area encompassed by regions (Link & Sauer, 2002).
to the product of the individual's residency state, R (0 = transient; 1 = resident), its alive state in time t − 1 (0 = dead or per- We defined a logit-linear model for m,t that allowed survival to vary as a function of a stratum-specific mean on logit scale (logit( 0[m] )), the winter drought index, cmd m,t , the tailwind covariate, tw m,t , and a zero-mean random stratum-specific year effect, m,t : We defined an analogous logit-linear model for m,t with the exception that we did not include the climate covariates. For the stratum-specific year effects in both models, we allowed precision hyperparameters to be stratum-specific.
We also developed models for the observation components of the MAPS CJS model. We modeled captures, c j,k,m,t at k = 1,…, K stations as a function of the true alive state and recapture probability, p k : c j,k,m,t ∼ Bern z j,m,t p k ; and observations of predetermined residency status, r j,k,m,t, as a function of the true residency state, R j,k,m,first(j) , and the predetermined residency probability, ρ k : r j,k,m,first(j) ∼ Bern R j,m,t k .
Predetermined residency was based on multiple within-season captures ≥10 days apart in the year of marking (r = 1 indicating observed residents; r = 0 for unknown residency status). We defined logit-linear models for p k and ρ k that included intercepts and zero-mean random station effects.
The reproductive index included in the recruitment component of the population process model was based on a binomial model applied to age-specific MAPS capture data (dotted yellow box in Figure 3).
Specifically, we modeled the number of young (hatching year) birds captured, HY k,m,t , according to HY k,m,t ∼ Binom(pHY k,m,t ,Nind k,m,t ), where pHY k,m,t represents the probability of capturing a young bird and Nind k,m,t represents the total number of individual birds captured at station k in region m and year t. We defined a logit-linear model for pHY j,m,t as: where the 0[m] are fixed stratum-specific intercepts; ef is the coefficient for an effort covariate, ef m,t ; cmd is the coefficient for the region-specific carry-over effect of winter drought (cmd t ) on productivity; the temp[m] are the region-specific effects of spring temperature (temp m,t ) on productivity; and the yr m,t and sta k are zero mean region × year and station effects, respectively. The ef s,t covariate was calculated as a ratio representing summed effort across MAPS banding sessions when young birds were captured relative to summed effort during sessions when adults were captured. Effort totals for young and adult capture sessions were weighted by regional totals for each age class across years, such that effort completed during MAPS sessions of peak captures counted more than MAPS effort during periods of fewer captures. We then derived the index of region-and year-specific postfledging productivity, RI m,t , on a per-capita scale (young/adult) as: We implemented the model with JAGS 3.3.0 (Plummer, 2003) using the jags function of the jagsUI package (Kellner, 2015) in the R statistical computing environment (R Core Team, 2018). We assigned vague uniform U(0, 1) prior distributions for inverse-logit transformed intercepts of models for parameters on 0-1 probability scales. Regression coefficients for fixed effects of linear models were modeled with Norm(0, 10 -3 ) priors, and standard deviation hyperparameters were modeled with U(0, 10) priors.
We inferred support for vital rate-covariate relationships for regression coefficients with 95% credible intervals that did not overlap zero. The first-year survival/immigration parameter, ι m,t , was determined based on a weakly informative prior distribution, ters (Gelman & Hill, 2006). We present all parameter estimates as means ± 95% credible intervals.

| Demographic contributions to variation in population change
We used transient life table response experiments (LTREs) to decompose temporal variation in population growth rates among vital rate and demographic structure components (Koons et al., 2017(Koons et al., , 2016. Specifically, we considered contributions of adult apparent survival, ϕ; productivity (of both new recruits/immigrants and survivors from previous time step), RI; first-year survival/recruitment, ι.
Following Koons et al. (2017), we also considered contributions of demographic age structure; however, as in that study, we found that temp m,t + yr m,t + sta k age structure contributed virtually nothing to explaining variation in population growth. Thus, we do not report those results here.
Finally, we examined annual changes in population growth rate in relation to changes in climate covariates to better understand how climate variation drives population dynamics.

| Population status, trends, and vital rate dynamics
Area-weighted estimates of our population size index, suggest that

| Climate covariate relationships
Adult apparent survival was negatively related to winter drought for the two California genetic groups for (Table 2, Figure 5a), but not for the pnw group (Table 2). We found no relationship between tailwind and adult apparent survival (Table 2). We found no evidence of a carry-over effect of winter drought effect on productivity [̂c md = 0.01 (−0.12, 0.13)]. Spring temperature was only weakly related to productivity for the cca group (Table 2); this relationship was much stronger for the sne and pnw groups (Table 2; Figure 5b).

| Demographic and climatic contributions to population growth
Annual changes in population growth rates were driven principally by recruitment components (RI and ), rather than by changes in adult apparent survival ( Figure 6). Our index of first-year (HY) survival, , had the largest effect on annual population changes for all three populations.
We found relatively little evidence of a relationship between winter drought conditions and population change for the two California populations (Figure 7a,c). Nevertheless, differences in drought conditions between years did appear to be related to differences in population growth between years (Figure 7b

| D ISCUSS I ON
Integrated population models have gained wide usage in population ecology because they provide a cohesive framework for understanding demographic and environmental drivers of population change (Koons et al., 2017;Schaub & Abadi, 2011). However, these models have received little attention in applications that combine multiple independent surveys across broad spatial and temporal scales (Ahrestani et al., 2017;Robinson et al., 2014;Zhao, Boomer, & Kendall, 2018). The IPM presented here provides a flexible framework for these broad-scale multi-site applications by modeling the survival and recruitment processes as functions of continuous random variables, rather than as functions of binomial and poisson processes typical of most IPM applications. We follow Alternative models for the capture-recapture data, such as Jolly Seber (JS) models (Link & Barker, 2005;Royle & Dorazio, 2008) could be incorporated to inform recruitment directly. However, we have found these models impractical with large data sets due to long computation times, especially whenever capture probabilities are low (as in our example) because of the necessary addition of large numbers of additional all-zero capture histories. Reverse-symmetry models based on site-level data summaries could also be used (Pradel, 1996;Tenan et al., 2014); however, we have encountered long implementation times with these models and difficulty achieving convergence We linked demographic parameters to climate covariates that will be impacted by climate change. Spring winds have been shown to affect survival in other passerine birds (Drake et al., 2014;Huang et al., 2017); however, we found no evidence of spring wind effects on these three Wilson's warbler groups. It is possible that our wind covariate did not properly capture migration conditions due to how migration regions were delineated or that mean tailwind conditions across the entire migration period was not an appropriate metric.  (Trenberth et al., 2013). Although our results suggested that adult apparent survival contributed less to overall population trends than recruitment, conserving and managing for drought-resilience habitats will nevertheless likely be an important component of any conservation plan for these populations.
Our analysis incorporated age-specific capture data to model postfledging productivity, which allows assessment of relative contributions of reproductive output and first-year survival and immigration to recruitment and population change. Although we found no evidence of winter drought carry-over effects on productivity, we did find population-specific effects of spring temperature on productivity: Productivity of populations in the montane (sne) and northerly (pnw) regions increased with increasing spring temperature. This finding is consistent with results of multi-species productivity models applied to MAPS data within the sne region ( S. Albert also provided administrative support. This work was supported by USGS grant no. G16AC00424 awarded to IBP. We thank K. Ruegg and R. Harrigan for shapefiles delineating genetic groups.
Any use of trade, product, or firm names in this publication is for descriptive purposes only and does not imply endorsement by the U.S.
Government. This is IBP Contribution Number 645.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
Both authors conceived the study. JS prepared the data, developed the model, analyzed the data, and led the writing of the manuscript.
Both authors contributed critically to drafts and approved the final version for publication.

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https ://doi.org/10.21429/ 04ma-p963.