1. Eco-geographic rules describe spatial patterns in biological trait variation and shed light on the drivers of such variation. In animals, a consensus is emerging that ‘pioneering’ traits may facilitate range shifts via a set of bold, aggressive, and stress-resilient traits. Many of these same traits are associated with more northern latitudes, and most range shifts in the northern hemisphere indicate northward movement. As a consequence, it is unclear whether pioneering traits are simply corollaries of existing latitudinal variation, or whether they override other well-trodden latitudinal patterning as a unique eco-geographic rule of phenotypic variation.
2. The tree swallow (Tachycineta bicolor) is a songbird undergoing a southward range shift in the eastern United States, in direct opposition of the poleward movement seen in most other native species’ range shifts. Because this organic range shift countervails the typical direction of movement, this case study provides for unique ecological insights on organisms and their ability to thrive in our changing world.
3. We sampled female birds across seven populations, quantifying behavioral, physiological, and morphological traits. We also used GIS and field data to quantify a core set of ecological factors with strong ties to these traits as well as female performance.
4. Females at more southern expansion sites displayed higher maternal aggression, higher baseline corticosterone, and more pronounced elevation of corticosterone following a standardized stressor, contrary to otherwise largely conserved latitudinal patterning in these traits. Microhabitat variation explained some quantitative phenotypic variation, but the expansion and historic ranges did not differ in openness, distance to water, or breeding density.
5. This countervailing range shift therefore suggests that pioneering traits are not simply corollaries of existing latitudinal variation, but rather, they may override other well-trodden latitudinal patterning as a unique eco-geographic rule of phenotypic variation.

Study sites and environmental data

Methods were approved by Appalachian State University IACUC #13-15 and US Master Banding Permit #23563. All animals were handled in such a way to reduce stress and avoid physical harm. All adults were released in their home territory. We conducted fieldwork during May-June 2015. Historical sites included Saukville, Wisconsin (43.382 N, 88.023 W), Long Point, Ontario (42.623 N, 80.465 W), and Wolfville, Nova Scotia (45.107 N, 64.378 W). Expansion sites included Bloomington, Indiana (39.142 N, 86.602 W); Ames, Iowa (42.073 N, 93.635 W); Davidson, North Carolina (35.438 N, 80.697 W); and Boone, North Carolina (36.196 N, 81.783 W). We recorded GPS coordinates at each nestbox (Garmin GPSmap 78s).

Sites were categorized as either historical or expansion based on prior publications (Lee, 1993; Shutler et al., 2012), bolstered by personal communications with local researchers and data from the Bird Breeding Survey (BBS; 1967 to 2017)(Sauer et al., 2017). Historical sites have abundant tree swallow breeding for >100 years (Winkler et al., 2020), whereas expansion sites have increased abundance since the 1960s. In the first 10 years of the BBS, Nova Scotia, Ontario, and Wisconsin (historic sites) reported an annual average (± SE) of 458 ± 87 breeding tree swallows, whereas Indiana, Iowa and North Carolina (expansion sites) reported only 2 ± 1. In the most recent 10 years, breeding numbers have increased by 40-fold in expansion sites (Figure 2b; visualized at the state-/province-level in Figure S1a). Although BBS data surely underestimate abundance, these trends reflect site-specific data (Lee, 1993; Shutler et al., 2012). Notably, not all expansion sites are at the southernmost range edge (i.e. Tennessee, the Carolinas), and some expansion sites (Davidson and Boone, North Carolina) have more recent histories than others (Iowa, Indiana). All expansion sites are nonetheless beyond the historic core of the species distribution.

We quantified three key ecological factors at each site, using measurements in the field as well as satellite data. The purpose of these data was two-fold: First, we sought to test whether the expansion and historic range differed in a core set of parameters with strong ties to tree swallow success. Second, we wanted to assess whether habitat variation might confound any range-related differences in traits. To achieve these goals, we characterized land use/land cover (LULC) using ArcGIS 10.2 (ESRI, Redlands, California). Analyses focused on a typical foraging range of 300m around each nestbox (McCarty and Winkler, 1999). LULC data were obtained from USGS National Land Cover Dataset (2011) or Gouvernement du Canada Land Use (2010), at 30x30m resolution. We selected LULC categories based on their relevance to tree swallows, which require open, wet habitat for aerial foraging of insect prey (Winkler et al., 2020). As such, we focused on the percent of land with open habitat (pasture, open water, or barren), and distance to water (streams, rivers, ponds, wetlands). In the field, we estimated conspecific density as the number of nestboxes within a 50m radius that were defended by other tree swallows; previous studies have linked density with aspects of aggressive behavior (Bentz et al., 2013). 

Adult capture, morphology, and blood sampling

Females were captured between 09:00–16:30h, during incubation or chick-rearing (expansion: 44 incubating and 48 provisioning females; historic: 44 incubating and 36 provisioning). For provisioning females, average chick age was 6.3 days post-hatch (95%CI 5.3–7.3 days). We started a timer as soon as we captured each bird and collected an initial blood sample (~40–80uL) from the wing vein. Samples that took >3min were not used because Cort can elevate rapidly after handling (Schoech et al., 2013). We measured wing length, tail length, and body mass, but excluded wing data due to inadvertent methodological differences among sites. Tail length is correlated with structural size (Bourret and Garant, 2017; Hainstock et al., 2010) and affects flying ability (Norberg, 1990), an important feature for aerial insectivores. We recorded age as “second year of life” (SY, or yearling) or “after second year of life” (ASY, or >1 year) using plumage coloration (Hussell, 1983). All birds were banded with an aluminum leg band, and we marked each female across the breast with a colored marker to facilitate identification in during behavioral assays (Whittingham and Dunn, 2001). After processing, we placed each bird in an opaque paper bag until 30 min post-capture, a standardized restraint protocol to measure Cort elevation, or ΔCort. We collected a second blood sample (~40–80uL) and released the bird shortly thereafter.

Assay of maternal aggression

Nest defense is a key component of adaptive maternal investment. We measured clutch and brood size and found no difference between expansion and historic sites (linear mixed model, LMM, with population as random effect, clutch: β = 0.38 ±0.32 SE, F_1,85=1.46, p=0.23; brood: β = -0.09 ±0.43 SE, F_1,77= 0.05, p=0.84), and we therefore focused on a behavioral aspect of maternal investment. Specifically, we assayed maternal aggression against a simulated nest predator ~24h after blood sampling, using an assay modified from (Winkler, 1992). Decoys were commercially-manufactured models of the American crow (Corvus brachyrhynchos), a widespread nest predator. We randomly rotated among six exemplars to limit psuedoreplication. Ahead of time, we affixed a decoy via wire, dangling from a ~0.7m pole. Trials began by visually identifying the female and then deploying the decoy alongside crow calls from the Cornell Lab of Ornithology. The observer quickly slid the pole onto the existing nestbox and pole hardware, rapidly suspending the decoy above the box in a semi-natural flight position. The observer retreated to ~40m and began the behavioral assay. For 5 min, we measured the number of dives towards the model, within 1m of the nestbox. In tree swallows, nest defense is a highly repeatable trait (Betini and Norris, 2012).

Hormone assays

We quantified plasma Cort levels using an enzyme immunoassay kit that has high accuracy and assay parallelism (Cayman ELISA #500655, see Rosvall et al., 2012). Briefly, we added 10uL plasma to 200uL ultrapure H₂O, extracted 3x with ether, dried with N2, and reconstituted with 600uL assay buffer. Each plate contained up to 33 samples in duplicate, an 8-point standard curve, blank, maximum binding, non-specific binding, and total activity controls, as well as 3 additional plasma pools used to calculate variability. Plates were balanced by site, breeding stage, and time point (baseline or 30-min). We read absorbance at 412 nm and interpolated concentration using Gen5 software (v2.09, BioTek, Winooski, VT, USA). Inter-plate variability was 10.2%, and intra-plate ranged from 3.7–12.4% (mean: 7.7%). We did not calculate extraction efficiency at a sample-level; instead, we spiked a separate set of samples (n=10) with 20uL H3-Cort (~2500 CPM), extracted 3x, and found average efficiencies to be 98.1%. Values may therefore underestimate true Cort concentrations, though our data are typical for this species (e.g. Zimmer et al., 2019), suggesting this effect is minimal.

Data Analyses

Statistical analyses were performed in R (v. 3.4.3, R Core Team, 2017). We report results as mean ± one standard error, and α=0.05. This study included data on 172 birds and the habitat surrounding each bird’s nestbox; however, we did not collect all data from every bird, e.g. because not all birds could be captured or identified during behavioral observations, some could not be bled quickly enough, some did not have full LULC data. Thus, sample sizes vary from 97 to 147 birds (Table 1).

To test whether ranges (expansion vs historical) differed in habitat parameters, we used a permutational multivariate analysis of variance (PERMANOVA) using the adonis2 function in vegan (Oksanen et al., 2016). We selected PERMANOVA because this non-parametric test is insensitive to multicollearity, and we found moderate correlations among parameters (0.34<|r|<0.53). We visualized differences with a principal component analysis (PCA) using the prcomp function.

We next tested how five putatively pioneering traits differ between historic and expansion ranges. A random effect of populations was included in all models. We performed separate analyses for each trait. LMM were used to analyze all traits, except for anti-predator aggression, which was over-dispersed count data better suited to a negative binomial generalized LMM. Models were run in lme4 (Bates et al., 2015), and we tested residuals for normality with a Shapiro-Wilk test. Two variables required log-transformation to meet model assumptions: baseline Cort and ΔCort. There were 7 individuals (n=5 historic; n=2 expansion) that showed negative ΔCort values. We had no indication that they were disturbed by us or something else, but their responses are not the biological norm, and so we removed them from subsequent analyses under the assumption that they do not represent the typical stress-induced Cort elevation that we seek to understand. Two extreme outliers were detected for tail length (Grubb’s test, p<0.05), and these points were removed prior to model comparison.

For each trait, we ran a series of models, all of which contained range (historic vs. expansion) and all combinations of: stage (incubating vs. provisioning), a range*stage interaction, PC1 of the habitat parameters, age, and capture time, in addition to the intercept only (null) model. Capture time was unrelated to anti-predator aggression and was not included in those models. We used Akaike Information Criterion corrected for small sample sizes (AICc) to find the most appropriate model (Burnham and Anderson, 2002). Models within 2 ΔAICc of the top model are equally fit, and, when this occurred, we present the most parsimonious model. Model selection was conducted with maximum likelihood estimation, but variable significance of the top model was done using restricted maximum likelihood. Both marginal (R²_marg; variance explained by fixed effects) and conditional (R²_cond; variance explained by fixed and random effects) R² are estimated using MuMIn (Bartoń, 2019). See Supplementary Materials for greater detail on (a) how we accounted for spatial parameters, (b) how we evaluated range vs. latitude as predictors of trait variation, and (c) how we used a heterogeneity analysis to assess the degree of population-specificity of any range-related trait differences.

Excel