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.

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^-1 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_i = | ATL_slope_i - MISS_slope_i |

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_i) to the mean of all other species’ slopes (ENA_slope_k) in the ENA dataset:

s_between_i = | ENA_slope_i - [1/n∑⁶³_k=1(ENA_slope_k)] |

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:

Winter latitude (continuous variable): 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.

Migration strategy (categorical variable): 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 (Bombycilla cedrorum), Pine Siskin (Spinus pinus), and Red-breasted Nuthatch (Sitta canadensis) 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.

Winter diet (continuous variable): 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, Lanius ludovicianus) or piscivorous species (e.g. Belted Kingfisher, Megaceryle alcyon) 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 s_within, 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_c) 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 σ_within, 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 (Σw_i) 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).

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).

Model selection and assessment of overall variable importance (i.e., the summed Akaike weights (Σw_i) 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).