DOI: 10.5061/dryad.4j0zpc8hq
title: Winter range shifts and their associations with species traits are heterogeneous in eastern North American birds

description: Many species’ distributions are shifting in response to climate change.
 Many distributional shifts are predictably poleward or higher in
 elevation, but heterogeneity in the rate and direction of shifts both
 within and between species appears to be common. We found high
 heterogeneity in the trajectory of winter range shifts for 65 species of
 birds across eastern North America and in the different traits and trait
 interactions associated with these shifts across the spatial scales we
 examined. We used data from the Christmas Bird Count to quantify the
 trajectory of winter latitudinal center of abundance range shifts over
 four decades (1980–2019) for 65 species of songbirds and woodpeckers in
 North America, both across eastern North America (ENA) as a whole and for
 the Atlantic (ATL) and Mississippi (MISS) flyways separately. We then used
 linear models and AICc model selection to test whether species traits
 could explain variation in range shifts or flyway discrepancies. Across
 ENA, most species showed northward latitudinal range shifts, but some
 showed no latitudinal shift while others shifted southwards. Amongst ATL
 and MISS, we documented both within- and between-species differences in
 the rate and direction of latitudinal shifts, complicating the results
 from across ENA. No single trait emerged as a dominant driver of range
 shift differences at the ENA and flyway scales. Migration strategy
 interacted with insectivory to explain variation at the largest spatial
 scale (ENA), whereas frugivory and mean winter latitude explained much of
 the variation in ATL and MISS, respectively. Exploring heterogeneity in
 range shifts within and between species, and in the associations between
 range shifts and life history traits, will help us better understand the
 mechanisms that mediate differing responses to environmental change and
 predict which species will be better able to adapt to those changes.
description: We requested Christmas Bird Count (CBC) data from 1980–2019 from
 the National Audubon Society for 65 species (Supplementary Material Table
 S1) that met the above criteria. We chose 1980 as a starting point based
 on data from Rushing et al. (2020) and Vose et al. (2017) that indicated
 that breeding distributional shifts of North American birds and
 temperatures started to change most noticeably in the mid-1980s. The CBC
 is an annual survey during a two-week window centered on 25 December.
 Surveys occur within 24.14-km diameter circles centered on the same point
 each year, with locations across North America (figure 1). We truncated
 this dataset to only those circles located in the MISS and ATL flyway
 states and provinces (shaded regions, figure 1). To avoid potential
 spatial-temporal bias that may arise when new CBC circles are initiated
 each year, we limited our analysis to only those circles that participated
 in at least 36 out of the 40 years (90%) of data collection during the
 study period. Thus, of the 1,961 CBC circles in the two flyways that
 collected data during any years of the study period, we used data from 629
 circles. The spatial distribution of CBC circles is biased towards the
 north and east of the study area (figure 1) (Meehan et al. 2019). This
 sampling bias may influence the static location of a species’ Latitudinal
 Center of Abundance (LCA), but it affects species similarly, and temporal
 changes in LCA location (our principal concern) arise only from
 differential abundance changes across sites over time, not from the static
 spatial sampling distribution. Latitudinal center of
 abundance calculation Latitudinal Center of Abundance
 is a standard expression of mean geographic range location, and LCA
 changes over time have been used to quantify range shifts (La Sorte and
 Thompson 2007, Paprocki et al. 2014). To calculate the LCA of each species
 for each year and flyway, we calculated the geographic mean of all CBC
 circles with that species present that year, weighted by the relative
 abundance of that species in each circle, using the R package geosphere
 (Hijmans et al. 2019). Because survey effort varies between circles and
 years, we used the number of birds per party-hour as an index of relative
 abundance for the weighted average to account for this variation in effort
 (Koenig and Liebhold 2016, Curley et al. 2020). For
 each species, we used linear regression to determine the strength and
 direction of the shift in LCA. We used year as the explanatory variable
 and LCA as the response variable, and the resulting regression slope
 estimated the rate of shift over the 40-year time period (see Figure 2).
 For our initial descriptive summaries of range shifts across species and
 flyways, we used only those slopes that were significantly different from
 zero as evidence for a range shift. However, all slopes were used
 unaltered for further statistical analysis (see below). To better
 illustrate effect sizes, we converted the regression slope from degrees
 latitude yr<sup>-1</sup> to cumulative distance (in km) during
 the study period using the conversion of 1 degree = 111 km multiplied by
 40 years. We performed these regressions for each species in eastern North
 America (ENA) as a whole, and for each flyway (ATL and MISS) separately to
 quantify variation in the direction and magnitude of range shifts between
 the flyways. We quantified flyway synchrony for each species by
 calculating Pearson’s correlation coefficient r between the
 flyway-specific LCA time series. This metric provides a quantitative
 comparison of the relative rate and direction of range shift between
 flyways for each species (see Figure 2). Within-
 versus between-species variation To compare
 between-species variation in the rate and direction of range shifts to
 within-species variation between flyways, we created and compared two
 indices. Our index of within-species variation was the absolute value of
 the difference between flyway regression slopes, for each
 species: s_within<sub>i</sub> = |
 ATL_slope<sub>i</sub> - MISS_slope<sub>i</sub>
 | For a comparable index of between-species variation,
 for each species, we calculated the absolute value of the difference
 between the ENA regression slope of that species
 (ENA_slope<sub>i</sub>) to the mean of all other species’
 slopes (ENA_slope<sub>k</sub>) in the ENA dataset:
 s_between<sub>i</sub> = |
 ENA_slope<sub>i</sub> -
 [1/n∑<sup>63</sup><sub>k=1</sub>(ENA_slope<sub>k</sub>)] | We converted slopes to cumulative distance (km) across the four-decade study period, as above.  These two indices allowed us to more directly compare within- and between-species variation in range shifts for a given species, because they are on the same scale. In particular, the between-species variation provides a meaningful benchmark against which to judge the importance of the within-species variation. Species-level traits We tested if species-level traits could explain variation in the strength and direction of winter range shifts in the study area as a whole, within flyways, and between flyways. We included the following species-level traits in our analysis: <em>Winter latitude (continuous variable): </em>We calculated the mean LCA for each species during the first three years of the study period (1980–1982). The LCA during the first three years provides a baseline mean latitude for each species, prior to shifts occurring over the subsequent four decades. We included this variable to test whether more northerly species were shifting at different rates than more southerly species. <em>Migration strategy (categorical variable):</em> We placed each species into one of four migration strategy categories: residents, short-distance migrants, moderate-distance migrants, and irruptive migrants. We removed the three irruptive species, Cedar Waxwing (<em>Bombycilla cedrorum</em>), Pine Siskin (<em>Spinus</em> <em>pinus</em>), and Red-breasted Nuthatch (<em>Sitta</em> <em>canadensis</em>) from the trait-based models (see below) due to small sample size for the category (three species), but we included them in the initial descriptive statistics. <em>Winter diet (continuous variable): </em>Using data from Billerman et al. (2020), we quantified winter diet by estimating insects, fruit, and seeds taken in the winter for each species as a percentage of the total. To simplify this analysis, we did not include carnivorous (e.g. Loggerhead Shrike, <em>Lanius</em> <em>ludovicianus</em>) or piscivorous species (e.g. Belted Kingfisher, <em>Megaceryle</em> <em>alcyon</em>) in this study.  Trait-based range shift models and model selection We created generalized linear models (glms) using the cumulative magnitudes of LCA range shift (in km) in ENA, ATL, and MISS as the response variables and species traits as explanatory variables, with interaction terms. For the flyway discrepancy response variable <em>s_within</em>, we created glms with Gamma error distribution to account for the log-normal distribution of this variable. We used all LCA regression slopes for the response variables rather than converting non-significant slopes to zero, which would have created a non-normal response variable, violating the assumptions of the linear models. In practice, most non-significant slopes were near zero (Figure 3). We used Akaike’s Information Criterion with small-sample size correction (AIC<sub>c</sub>) to evaluate models with and without interactions (Burnham and Anderson 2002, Burnham et al. 2011). We performed model selection separately for ENA, ATL, MISS, and for <em>σ_within</em>, examining the same set of models in each case. The model set included all combinations of the trait variables winter latitude, migration strategy, and winter diet. Models that included latitude or % fruit in diet also included a quadratic component for those variables, based on unimodal (hump-shaped) relationships with range shifts suggested in preliminary scatter plot visual assessments (Supplementary Material Figure S1). Winter diet was represented by the three continuous variables % insects, % fruit, and % seeds, but only one of these was included in any given model to avoid collinearity. Additive and interaction models were included for each variable combination, resulting in a total of 26 models including an intercept-only model (Supplementary Material Table S2). Model selection and assessment of overall variable importance (i.e., the summed Akaike weights (<em>Σw<sub>i</sub></em>) of all models that included the variable) were performed using the MuMIn package in R (Barton 2022). The glms and all other numerical analyses described above were performed in R version 4.0.5 (R Core Team 2021).
description: To calculate the LCA of each species for each year and flyway, we
 calculated the geographic mean of all CBC circles with that species
 present that year, weighted by the relative abundance of that species in
 each circle, using the R package <em>geosphere</em> (Hijmans
 et al. 2019). Model selection and assessment of overall
 variable importance (i.e., the summed Akaike weights
 (<em>Σw<sub>i</sub></em>) of all models that
 included the variable) were performed using the MuMIn package in R (Barton
 2022). The glms and all other numerical analyses described above were
 performed in R version 4.0.5 (R Core Team 2021).