DOI: 10.5061/dryad.280gb5msz
title: The role of tropical rainfall in driving range dynamics for a long-distance migratory bird

description: Predicting how the range dynamics of migratory species will respond to
 climate change requires a mechanistic understanding of the factors that
 operate across the annual cycle to control the distribution and abundance
 of a species. Here we use multiple lines of evidence to reveal that
 environmental conditions during the nonbreeding season influence range
 dynamics across the lifecycle of a migratory songbird, the American
 redstart (Setophaga ruticilla). Using long-term data from the nonbreeding
 grounds and breeding origin estimated from stable hydrogen isotopes in
 tail feathers, we found that the relationship between nonbreeding season
 survival and migration distance is mediated by precipitation, but only
 during dry years. A long-term drying trend throughout the Caribbean is
 associated with higher mortality for individuals from the northern portion
 of the species’ breeding range, resulting in an approximate 500 km
 southward shift in breeding origins of this Jamaican population over the
 past 30 years. This shift in connectivity is mirrored by changes in the
 redstarts breeding distribution of abundance. These results demonstrate
 that the climatic effects on demographic processes originating during the
 tropical nonbreeding season is actively shaping range dynamics in a
 migratory bird.
description: See the manuscript and supplementary materials for further
 details on model formulation. The data included here include BBS and
 environmental data for the abundance and occupancy analysis. Geolocator
 data can be found in https://doi.org/10.1002/ecy.3938. 
description: # The role of tropical rainfall in driving range dynamics for a
 long-distance migratory bird
 *Hierarchical Breeding
 Bird Survey Abundance Model* To evaluate the effect of winter rainfall on
 population abundances of American redstarts on the breeding grounds we
 relied on 30 years (1990-2019) years of the long-term Breeding Bird Survey
 (BBS) data to coincide with our 30-year banding and isotope dataset. Data
 from 554 BBS routes included in this analysis were selected based on the
 likely breeding origins of this population (*SI Appendix,* Figure S4).
 Stable-hydrogen isotope analysis indicated that this wintering population
 originates from a broad latitudinal band (35 to 50 °N). To determine the
 longitudinal band of our population we calculated the 95% CI of breeding
 longitude estimated from a sample of geolocator tagged individuals (-93 to
 -83 °W). The model includes estimates of log-linear population trends
 (![]()).  ![]()represents the effect of week-of-year on abundance,
 *β*<sub>3 </sub>accounts for breeding season (May – July)
 Total Rainfall, *β*<sub>4 </sub>accounts for breeding season
 (May – July) average max temperature,
 and *β*<sub>5</sub> represents the effect of total winter
 rainfall (Jan – Mar) in Jamaica. ![]() accounts for the effect of Latitude
 on abundance and ![]() allows for a varying effect of winter rainfall by
 breeding latitude. We used uninformative normal priors with mean 0 and
 variance 10<sup>2</sup> for all beta parameters and diffuse
 uniform distributions for all variance parameters. To account for the
 non-independence of observations on a given route, we have included a
 random effect on the intercept for each route
 (*Route<sub>i</sub>*) and observer
 (*Observer<sub>i,t</sub>*), both with a Gaussian prior
 centered on 0. To compare the relative contributions of each factor on
 breeding bird abundance and to aid in computations, we centered and scaled
 all predictors so that values were in units of SD about the mean.     We
 ran three chains of 25000 iterations, discarded the first 5000, thinned by
 a factor of 10 to give 6,000 samples from the posterior distribution for
 inference. We assessed model convergence through the parameter histogram
 plots and Gelman-Rubin (R-hat) convergence diagnostics. We interpreted all
 effects by examining posterior mean and associated 95% credible intervals
 (CI) for all β<sub>i</sub> parameters.   *Hierarchical
 Multi-Season Occupancy Model* We modeled changes in occupancy probability
 with multi-season Bayesian hierarchical occupancy model implemented in the
 *spOccupancy* package in R (60). We utilized the BBS data from the
 abundance model above as occurrence data (presence/absence) by treating
 each of the 5 segments along a BBS route (pooled 10-stop summaries, 50
 stops in total per route) as a replicate and converting all counts &gt; 1
 as a detections and counts &lt; 1 as an absence. Conventional practice in
 occupancy modeling is to utilize temporal replicates to estimate detection
 probability. BBS routes are only sampled once per year and so occupancy
 models have been adapted to utilize spatial replicates (stops) as a
 substitute (61). However, it’s important to note that this space-for-time
 substitution introduces potential spatial autocorrelation that may bias
 occupancy if not properly accounted for. We attempt to account for the
 correlation by including a route-level random effect on detection
 probability which would allow for stops within a route to be correlated.
 Ideally, modeling the spatial correlation directly would provide a more
 exact approach towards accounting for the survey structure (61). We
 modeled detection as a function the week of the year the route was run
 (centered and standardized) while accounting for random observer affects
 and route effects. Occupancy was modeled as a function of proportion of
 forest, a year trend (centered and standardized), and random route
 effects. The proportion of forest represents the area within 40km of a BBS
 route with forest cover &gt; 10%. These data were extracted from the
 Global Forest Change dataset (v1.6; 62) using Google Earth Engine. We
 incorporated a AR(1) temporal covariance structure to formally account for
 its effect on occupancy dynamics. We fit the model with default
 non-informative on all ![]() and weakly informative priors on the AR(1)
 variance and correlation parameters following Doser et al. (2022). We ran
 3 chains for 10000 iterations with a burn in of 2000 iterations and
 thinning rate of 10 resulting in 2400 posterior samples. We assessed model
 convergence through the parameter histogram plots and Gelman-Rubin (R-hat)
 convergence diagnostics. The significance of parameter estimates was
 evaluated by whether or not the 95% credible interval overlapped 0. ##
 Description of the data and file structure Survival Data are available
 upon request; however, the code to run survival analysis is included here
 as Survival Analysis_Final.r. BBS_Abundance_Data.rds contains all
 necessary data to replicate the abundance analysis, including winter and
 breeding season environmental covariates scaled and centered. Below is a
 quick description of the variables: Count = # of Redstarts observed Route
 = Route Identification Number Observer = Observer Identification Number
 Year = Year of Route scaled to number of years since 1989. Year 1989 would
 be coded as 0. Week = Week of the Year with 0 indicating the first week of
 the year Latitude = scaled and centered latitude (degrees north converted
 to units of standard deviation) BreedingRain = Total Rainfall in mm (May -
 July) scaled and centered (units of standard deviation) BreedingTemp =
 Average Daily Max Temperature Celsius (May - July) scaled and centered
 (units of standard deviation) WinterRain = Total Winter Rain in Jamaica mm
 (Jan - Mar) scaled and centered (units of standard deviation) NAs in
 count, observer, and week represent instances where the route survey was
 not conducted. The model formulation is presented in the
 BBS_Abundance_Model.r script but can easily be implemented in the brms
 package with model description in manuscript methods.
 BBS_Occupancy_Data.rds contains all the necessary data in the proper
 format to replicate the occupancy analysis in BBS_Occupancy_Analysis.r.
 BBS Occupancy data are presented in list format. The first list item are
 the observed counts in an array. Rows represent the different routes,
 columns represent the different years, and matrix slices represent 10-stop
 groups (e.g., Count10, Count20, Count30, Count40, Count50). NA represent
 instances when the route was not surveyed. The next dataset within the
 list is the spatial coordinates of each route. This is given in latitude
 (degree North) and longitude format (degree East). The next dataset within
 the list includes all of the detection covariates. This includes Observer,
 site, and week of the year. Observer is stored as a matrix organized by
 route and year. The site includes a unit identifier for the route, and the
 week is organized like an observer but with a week of the year denoted as
 above (weeks since Jan 1) scaled and centered (units of standard
 deviation). The final dataset within the occupancy data list includes all
 of the occupancy covariates. These include year, proportion of forest, and
 site. The year is organized by route and year, the site is organized as it
 was in the detection covariates, and the proportion of forest is organized
 at the route level and represents the total area covered with forest
 (fraction from 0-1) within 40km of the route center point. ##
 Sharing/Access information Isotope data are available upon request.
 Geolocator data are available in Dossman et al. (2022)
 .