Within-Site Variation in Feather Stable Hydrogen Isotope (δ2Hf) Values of Boreal Songbirds: Implications for Assignment to Molt Origin

Understanding bird migration and dispersal is important to inform full life-cycle conservation planning. Stable hydrogen isotope ratios from feathers (δ2Hf) can be linked to amount-weighted long-term, growing season precipitation δ2H (δ2Hp) surfaces to create δ2Hf isoscapes for assignment to molt origin. However, transfer functions linking δ2Hp with δ2Hf are influenced by physiological and environmental processes. A better understanding of the causes and consequences of variation in δ2Hf values among individuals and species will improve the predictive ability of geographic assignment tests. We tested for effects of species, land cover, forage substrate, nest substrate, diet composition, body mass, sex, and phylogenetic relatedness on δ2Hf from individuals at least two years old of 21 songbird species captured during the same breeding season at a site in northeastern Alberta, Canada. For four species, we also tested for a year × species interaction effect on δ2Hf. A model including species as single predictor received the most support (AIC weight = 0.74) in explaining variation in δ2Hf. A species-specific variance parameter was part of all best-ranked models, suggesting variation in δ2Hf was not consistent among species. The second best-ranked model included a forage substrate × diet interaction term (AIC weight = 0.16). There was a significant year × species interaction effect on δ2Hf suggesting that interspecific differences in δ2Hf can differ among years. Our results suggest that within- and among-year interspecific variation in δ2Hf is the most important source of variance typically not being explicitly quantified in geographic assignment tests using non-specific transfer functions to convert δ2Hp into δ2Hf. However, this source of variation is consistent with the range of variation from the transfer functions most commonly being propagated in assignment tests of geographic origins for passerines breeding in North America.


Introduction
Each year, billions of birds migrate between breeding and wintering areas [1]. Migration and dispersal movements have been difficult to quantify for most of these species. Yet, this information is important to identify the spatial scale at which population dynamics are taking place and should be considered in conservation planning (e.g. [2,3]). Recent studies have highlighted the importance of quantifying the level of migratory connectivity and extent of natal dispersal, i.e. straightline distance moved by an individual from its natal area to its first breeding site, in full life-cycle assessments [4][5][6]. High migratory connectivity implies that the majority of individuals in a breeding population use the same wintering areas, whereas low migratory connectivity suggests sparsely distributed individuals from a breeding population across the species wintering grounds [7].
Intrinsic markers like naturally occurring stable isotopes of several elements can provide important information on where birds grow their feathers and have the advantage of not requiring individuals to be recaptured when assessing movement dynamics. Among stable isotopes that have been used in movement studies (i.e. C, N, H, O, and S; [8]), stable hydrogen isotope ratios ( 2 H: 1 H; depicted as δ 2 H) are now the most commonly used marker [9]. Predictable spatial variation resulting from amount-weighted, growing-season mean precipitation (δ 2 H p ) worldwide is well documented [10]. The resulting patterns of spatial variation in δ 2 H p (i.e. isoscapes) have been used in combination with δ 2 H from feathers (δ 2 H f ) to infer geographic origin of birds as the H in bird feathers is ultimately derived from environmental waters where the feather was grown [11,12]. In the Nearctic-Neotropical migration system, feathers that are grown prior to fall migration (pre-basic molt; [13]) may provide information on natal or previous breeding origin in North America, whereas feathers that are grown prior to spring migration (pre-alternate molt) can be used to identify the wintering grounds of breeding individuals [14][15][16]. However, before assigning birds to geographic origins, δ 2 H p values need to be adjusted to reflect the difference between local precipitation driving food webs (δ 2 H p ) and δ 2 H f (generally from regression models; hereafter "transfer functions"; [8]). Clark et al. [11] provided transfer functions for waterfowl and songbirds sampled in North America. More recently, Hobson et al. [17] found that foraging and migratory strategies were an important source of variation that influenced the transfer function applied to migratory songbirds. They suggested that δ 2 H f differences between ground and non-ground foragers may result from differences in ground-level evapotranspiration [17] or differences in trophic level [18].
While differences between groups of species in the transfer function between δ 2 H f and δ 2 H p have been documented, what causes this interspecific variation is not as well understood. Analytical approaches have been developed to deal with this uncertainty in Bayesian assignment tests to provide more accurate and less biased estimates of likely geographic origin [19]. However, a better understanding of variation in δ 2 H f among species with the same geographic origin should help in assessing those factors contributing to variation we see among various transfer functions. For example, evaporative processes in response to increased activity or ambient conditions may lead to higher δ 2 H in tissues [20]. Thus, birds nesting in open areas (e.g. unforested and forest canopy) could have higher δ 2 H f values relative to species associated with forest understory. Increased body water loss due to metabolic activity can result in enriched heavy isotopes in the body water pool, in turn leading to higher δ 2 H f values, and this effect might be more important in smaller individuals/species [21,22]. However, Betini et al. [23] reported the opposite pattern, where body mass of nestling Tree Swallows (Trachycineta bicolor) was positively correlated with δ 2 H from blood samples. Phylogenetically conserved life history traits [24,25] and, potentially, conserved biochemical pathways controlling H isotope discrimination in tissues [26,27] among closely related species suggest a potentially important phylogenetic component of δ 2 H f that could explain important interspecific variation. Inter-annual variation in δ 2 H f has been reported for a few songbird species [28,29] and, all things being equal, this variation should be consistent across species nesting in the same region but, to our knowledge, no studies have tested this assumption by quantifying the species × year interaction effect on variation in δ 2 H f of forest songbirds. Finally, different species may seek out different local microhabitats during the molt period and these might also account for differences in δ 2 H f from the same general area but little is known about habitat use during this typically "secretive" period of the annual cycle (e.g. [30]).
In this study, we used tail feathers collected from different songbird species during the same breeding season and geographic area (within a ca. 130 km radius) to test various hypotheses that have been proposed to explain intraspecific and interspecific variation in δ 2 H f . Specifically, we used a model selection approach to compare the relative importance of species, land cover type, forage substrate, nest substrate, diet composition, body mass, sex, and phylogenetic relatedness in explaining variation in δ 2 H f among experienced breeders (after second-year; ASY) of 21 songbird species (7 families and 16 genera). We predicted that ground foragers and insectivores (i.e. species feeding almost exclusively on invertebrates during the breeding season) would have higher δ 2 H f compared to non-ground foragers and omnivores (i.e. species feeding on seeds, fruits, and insects during the breeding season), respectively. Observed δ 2 H f for ground and non-ground foragers were used to validate predicted values from the transfer functions provided by Hobson et al. [17]. We also predicted that individuals nesting in open areas would have higher δ 2 H f than those nesting in forest understory. Finally, we predicted that closely related species (i.e. within the same genera) would have more similar δ 2 H f values than less related individuals. There were no a priori expectations regarding the direction of the effect (positive or negative) of body mass and sex (male or female having higher δ 2 H f ) on δ 2 H f . Females are more likely to experience breeding dispersal movements, i.e. straight-line distance moved by an individual from a breeding territory to another in subsequent years, and over larger distances than males [31]. Thus, we predicted that females would have similar mean δ 2 H f values than males, but larger variation, assuming no bias in direction of dispersal movements. For four of our most common species, we also tested for a year × species interaction effect on δ 2 H f to determine whether annual variations in δ 2 H f are consistent among species.

Study area and feather samples
The study was conducted at bird banding stations in the Lesser Slave Lake Provincial Park (Lesser Slave Lake Bird Observatory, LSLBO; 55°20' N, 114°40' W) and the lower Athabasca oil sands region (Owl Moon Environmental Inc.; 56°43' N, 111°22' W) in northeastern Alberta, Canada. The region is characterized by conifer (black spruce, Picea mariana; white spruce, Picea glauca; jack pine, Pinus banksiana), deciduous (trembling aspen, Populus tremuloides), and mixedwood stands typical of the western boreal forest. The banding stations in the oil sands region are similarly forested, but also occur in riparian areas and reclaimed mine sites.
Breeding birds have been captured at banding stations as part of the Monitoring Avian Productivity and Survivorship (MAPS) program [32]. In 2013, 32 MAPS stations were monitored in the oil sands region (Table 1). Each banding station operated from 8 to 14 mist nets (each net was 12 m × 2.6 m) and captured birds passively (i.e. no attempts were made to attract the birds). Starting at sunrise, birds were captured at each station for 6 hours every 10 days throughout the breeding season (mid-June-mid-August 2013; periods 5-10; [32]). The same year, LSLBO monitored 4 MAPS stations within a forested area of approximately 3 ha that bordered the eastern shore of Lesser Slave Lake (Table 1). In 2011, active netting was conducted within ca. 1 km of the closest MAPS station where males of four species (Ovenbird; Seiurus aurocapilla, Swainson's Thrush; Catharus ustulatus, American Redstart; Setophaga ruticilla, and Yellow-rumped Warbler; Setophaga coronata) were attracted to mist nets using conspecific song playback. Additional samples were also collected from these 4 MAPS stations. Two 3 rd rectrices were plucked from as many captured individuals as possible and individual body mass was recorded at time of capture. All feathers were stored in paper envelopes at room temperature. In the oil sands region, birds were captured and feathers were collected by KF, CG (banders-in-charge), and their banding staff. At Slave Lake, birds were captured and feather samples were collected by RK (bander-in-charge), LSLBO banding staff, and volunteers. Only feathers from males and females of at least two years old (i.e. after-second year individuals; ASY) were used for isotope analysis. These individuals were selected based on the higher site fidelity reported in adult songbirds compared to first-year breeders [31,34]. Also, adult Neotropical migratory passerines generally undergo complete pre-basic molts, including the replacement of primaries and rectrices, at or near their breeding site [13]. Thus, for most species, our feather samples should have reflected the isotopic signature incorporated at our study site in the previous breeding season. Individuals from LSLBO were aged and sexed by the banders-in-charge, while those from the oil sands region were aged and sexed by the banders-in-charge and photos were reviewed by P. Pyle (Institute for Bird Populations; http://www.birdpop.org/). We further determined age of some S. aurocapilla by quantifying the wear pattern of the 3 rd rectrices [35]. Feather samples were required from a minimum of nine individuals for a species to be considered. For each location, we also identified the dominant land cover type at the point level and in a 100 m radius (exposed land, wetland-shrub, wetland-treed, shrub-tall, broadleaf, coniferous and water), and wet area (ha) within a 100 m radius using the Earth Observation for Sustainable Development of Forests (EOSD; [33]). For each species, we assigned dietary categories (i.e. insectivore vs. omnivore), forage substrate (i.e. upper canopy, lower canopy/shrub, and ground), and nest substrate (i.e. agriculture, bogs, tree / shrub swamp, coniferous, deciduous, early successional, marsh, mixedwood, and open stands) using the Avian Life History Information Database (hereafter "ALHD"; http://www.on.ec.gc.ca/wildlife/wildspace/project.cfm; Table 2).
Fieldwork was conducted in accordance with the Canadian Council for Animal Care and the permit for this study was approved by the University of Alberta Animal Care Committee (permit # AUP00000100). Federal and provincial bird banding and feather collection permits for this study were approved by the Canadian Wildlife Service and Alberta Environment and Sustainable Resource Development, respectively.

Stable isotope analysis
Surface oils were removed from all feathers by using a 2:1 chloroform:methanol solution. The central vane of each feather was weighed (350 ± 20 μg) and samples were secured in silver capsules. Isotopic analysis was conducted at the Colorado Plateau Stable Isotope Laboratory at the Northern Arizona University for 2011 samples and the Stable Isotope Hydrology and Ecology Laboratory of Environment Canada for 2013 samples. Samples were exposed to high temperature (1350°C) flash pyrolysis and the separated H 2 pulses were used to measure δ 2 H f by continuous-flow isotope-ratio mass spectrometry (CF-IRMS). For both laboratories, we used the comparative equilibration approach with the same in-house keratin working standards (KHS [-54.1‰], SPK [-121.6‰], CBS [-197‰]) to account for exchangeable hydrogen in keratins where δ 2 H of nonexchangeable H was established [36]. Thus we are confident that isotope results are comparable within measurement error between the two laboratories (see [8]). Results were expressed for non-exchangeable H delta notation (δ 2 H f ) in units of per mil (‰) and the analytical error, based on within-run replicates of keratin reference standards (n = 5 per run) was ±2‰. Results are reported relative to the Vienna Standard Mean Ocean Water-Standard Light Antarctic Precipitation (VSMOW-SLAP) scale.

Statistical analyses
Linear mixed models were used to explore how δ 2 H f values from individuals of the selected species captured in 2013 were influenced by species, individual, and land cover factors (fixed effects) and phylogenetic relationships across species (random effects). A Multivariate Normal distribution was specified for each model where fixed effects are part of the mean specification, while random effects are part of the covariance structure. The response variable (y ij ) was δ 2 H f for individual i (i = 1. . .n) of species j (j = 1. . .s). We examined 15 different fixed effects as possible explanatory factors for δ 2 H f and compared them to a (0) null (intercept only) model. The , and 8-latitude), individual (9-linear and 10-quadratic individual mass, and 11sex), and interactive effects (12-forage substrate × diet, 13-forage substrate × land cover, 14 -nest substrate × diet, and 15-forage substrate × diet × nest substrate; Table 3). We used the Akaike Information Criterion corrected for small sample sizes (AICc; [37]) to select the most parsimonious model describing variation in δ 2 H f among feather samples.
Residual variance (σ 2 ) was included as a random effect in our 15 models and we divided each model into two subsets where variance in δ 2 H f among species was estimated as homoscedastic (variation among species was assumed equal [σ 2 ]; Hom1-Hom15) and heteroskedastic (the assumption of variance equality among species was relaxed [σ 2 j ]; Het1 -Het15). To examine the phylogenetic correlation in δ 2 H f , we estimated a multiplier of the off-diagonal elements in the Multivariate Normal covariance matrix defined as λ [38,39]. A value of λ = 0 indicated evolution of traits independent of phylogeny, while value of λ = 1 indicated traits evolving according to Brownian motion. Values <1 indicated that phylogeny effect was weaker than expected under the Brownian model. The covariance between two observations (y ij , y ik ) and corresponding Mutivariate Normal means (μ ij , μ ik ) was defined as cov(y ij , y ij ) = σ j σ k t jk λ, where t jk is proportional to the shared evolutionary path length between species j and k. Shared evolutionary path lengths were estimated using 5000 pseudo-posterior trees based on genetic data [40]. A phylogenetic correlation matrix [41] was constructed for each of the 5000 trees. Elements of the correlation matrix were defined as lengths of branches shared between species based on mean node heights [42,43]. The different number of observations per species meant we could not apply existing phylogenetic mixed model software which uses only 1 observation per species. Therefore, we implemented a maximum likelihood estimating procedure for Multivariate Normal mixed models using Markov-chain Monte Carlo methods and a data cloning algorithm [44] using the 'dclone' R package [45] and JAGS [46].
A two-way ANOVA was used to test for a Year (2011 and 2013) × Species (Ovenbird, Swainson's Thrush, American Redstart and Yellow-rumped Warbler) interaction effect on δ 2 H f . All analyses were conducted using R version 3.1.0 [47]. Results are presented as mean δ 2 H f ± standard deviation (SD) unless specified otherwise. We also extracted δ 2 H p values for our study area from the isoscape provided by Bowen et al. [48] (see also IsoMap, isomap.org). Values were expressed as the annual mean growing season (i.e. months where mean monthly temp > 0°C) and used to estimate δ 2 H f for omnivores and insectivores following Hobson et al. [17].

Results
We analyzed δ 2 H f from 278 ASY individuals of 21 songbird species captured in 2013 ( Table 2). Mean values for most species (15 species) ranged from -127‰ (Swainson's Thrush) to -162‰  [38]). Each model also described processes acting at different levels (i.e. Individual, Species or Station). a Description of variables: ForSub = ground vs non-ground foragers, Diet = insectivore vs omnivore, NestSub = typical nesting habitat, LandStation, Land100m and WetArea = land cover at station, dominant land cover within 100 m radius, and wetland are within 100 m according to EOSD [33], respectively, IndMass and IndMass 2 = linear and quadratic individual mass at time of capture, respectively. (Savannah Sparrow), while six species (Chipping Sparrow, Tree Swallow, Red-eyed Vireo, Least Flycatcher, Cedar Waxwing, Alder Flycatcher) had higher mean δ 2 H f (-97‰ to -51‰; Fig 1). Combined with large variation around mean values for these six species, these results suggested that rectrices from these individuals were grown on the wintering ground or during migration. These results, at the exception of those for Cedar Waxwing, are consistent with previous accounts [13,49] and these six species were removed from further analyses. The remaining samples included 192 ASYs (43 females and 149 males) from 15 species. We also report δ 2 H f for 113 ASY of the 4 species captured in 2011 (Fig 2).
Standard deviations were smallest for Mourning Warbler (± 4‰) and largest for Song Sparrow (± 22‰ ; Fig 1). The 95% confidence intervals for species-specific mean δ 2 H f overlapped with the predicted range of δ 2 H f values from the transfer functions of Hobson et al. [16] for ground and non-ground foragers in our study area (Fig 1). Mean values of ground foragers and non-ground foragers (irrespective of species) were -143 ± 14‰ and -144 ± 15‰, respectively,  [17]'s δ 2 H f transfer functions for ground (Ground) and non-ground foraging (Nonground) songbirds, respectively. The phylogenetic tree was derived from 5000 pseudo-posterior trees; shorter branch lengths occur between species that are more closely related with each other. Also indicated are the sample size for each species and the six species removed from further analyses because they were believed to have molted their tail feathers on the wintering grounds or during migration. both 95% confidence intervals overlap with the respective δ 2 H f estimates from Hobson et al. [17] for our study area. Males and females had similar mean and standard deviation in δ 2 H f (-143 ± 14‰ and -146 ± 11‰, respectively).

Discussion
This study demonstrates important interspecific differences in δ 2 H f among songbird species breeding in the same region. Species was the strongest predictor of the variation in δ 2 H f among individuals followed by a nesting substrate × diet composition interaction effect. Limited support was found for single effect models that included life history traits, individual (sex and body mass), or land cover predictors. Some studies have found support for these predictors [17,23,50], but they did not compare the relative importance of many predictors for different levels of organization. We also showed that differences in δ 2 H f between years are not always consistent among species. Our results provide important insights regarding different sources of error being propagated in Bayesian assignment tests and suggest that the use of species-and year-specific δ 2 H f isoscapes could improve the accuracy of assigned geographic origin of birds (see also [51], but see [52]).
Based on results from Hobson et al. [17], we predicted that ground foraging songbirds would have higher δ 2 H f compared to those foraging in the upper-or lower-canopy. Mean δ 2 H f values or corresponding 95% confidence intervals for the 15 focal species did indeed overlap with predicted values from the transfer functions of Hobson et al. [17] applied to our study area. However, models with foraging substrate as a single predictor received very little support (ΔAICc > 26.6). There was also no evidence that songbirds using open habitat or forest canopy during the nesting period had higher δ 2 H f compared to those using forest understories. In temperate or boreal ecosystems, differences in ground temperature between closed and open habitats may not influence evapotranspiration rates and, ultimately, δ 2 H f . For example, Hache et al. [29] found no difference in δ 2 H f between nestlings from recent selection harvesting (open) vs unharvested stands (closed). However, most adult songbirds start molting during the postfledgling period and some mature forest specialists have been shown to aggregate in open- Species-Specific Feather Hydrogen Isotope Ratios canopy clearcuts for concealment from predators during this period [30,53]. Thus, many focal species might have used similar microhabitat and experienced similar temperatures and sun exposure during the molting period irrespective of their nesting substrate. Species-specific variation in timing and location of molt (e.g. [54]) could also result in important spatio-temporal variation in δ 2 H f of ASY songbirds from a given region, but it is unclear to what extent resources consumed during the nesting period can be integrated in tissues that have grown during the post-fledging period (e.g. [55,56]).
Birds with larger body sizes could have lower δ 2 H f than smaller birds because larger individuals or species lose proportionately less of their total body water to evaporative processes [21,22]. However, there was no evidence in our study to suggest body size was important for δ 2 H f in adult songbirds, but differences in body mass among individuals at the time of capture might not reflect differences during the molting period. Females also did not have different δ 2 H f values than males [57] and variation around mean δ 2 H f was only slightly lower for males than females. These results indicate that population level assignments can be generated irrespective of body size and sex of individuals.
Breeding dispersal distances appear to be species-specific [57,58] and would often occur over relatively short distances. Variation in mean δ 2 H f can be used as a proxy for breeding dispersal rate and distance, i.e. species with larger variation might have a higher proportion of individuals that experienced breeding dispersal movements and/or dispersed over larger spatial extents. The relatively small variation around mean δ 2 H f reported for all focal species in this study suggests that most dispersal events would have occurred within the isocline corresponding to predicted δ 2 H f for our study area. However, differences in mean δ 2 H f among species could also reflect population-level post-breeding movements (i.e. within year dispersal movements prior to fall migration) and molt occurring elsewhere in the species breeding range (e.g. [59]).
The better performance of heteroscedastic compared to homoscedastic variance models suggests an important species-specific component for future songbird assignment of geographic origins using δ 2 H f . Differences in variation around species-specific mean δ 2 H f in a region might reflect interspecific differences in breeding dispersal rates, but also differences in niche breath. Habitat generalists use a broader range of land cover types which could result in greater variation in δ 2 H f than habitat specialists. However, species-specific tolerance (i.e. the range of environmental conditions over which a species occurs; [60]) was not a good predictor of variation in δ 2 H f for 13 of our focal species (R 2 = 0.05; p = 0.46; S1 File). However, the interspecific variation in the range of environmental conditions used by forest songbirds might not be correlated with the range of δ 2 H p values encountered by these species.
Studies have shown morphological traits to be more similar in closely related passerine species while for other traits, such as behavior and habitat associations, this might not be the case [24]. In this study, closely related species did not have more similar δ 2 H f values than less related species. Our 15 focal species represented only 3 families (i.e. Passeridae, Parulidae, and Turdidae) and may have been too closely related to detect phylogenetic effects. We also expected a consistent difference in δ 2 H f between years (2011 vs 2013) for the 4 focal species, but this effect was not observed. If year effects on δ 2 H f were exclusively driven by differences in δ 2 H p across years, we would expect a proportionate year effect across species breeding in the same region. Future studies should aim to identify the underlying mechanisms explaining this interaction effect. This could potentially be achieved by quantifying differences in δ 2 H f for a broader range of species and years. Nonetheless, along with species-specific variation in δ 2 H f , this species × year interaction effect was another important factor and likely corresponds to a large proportion of the variance being propagated in Bayesian assignment tests [52].
We provided evidence suggesting that the most important source of error being propagated in assignment tests of geographic origin of birds is likely interspecific variation in δ 2 H f. However, we had limited success identifying the underlying mechanisms. Future studies should determine whether missing variables or interaction terms could explain interspecific variation in δ 2 H f or more precisely determine where birds undergo the prebasic molt relative to breeding territories. A better understanding of the variation in δ 2 H f within and across years (i.e. across a larger number of species and years) should provide important information about sources of error being propagated in assignments of geographic origin. Transfer functions to explicitly account for these sources of error would likely help generate more precise and accurate assignments and provide better information on migratory bird movements to inform full-cycle conservation strategies. However, our ability to integrate these sources of variance in assignment tests will depend on the objectives. For example, it might be unrealistic to expect that speciesand year-specific δ 2 H f isoscapes would be available over large spatial extent. Alternatively, it would be possible to document regional dispersal movements [3,4]. Results from Vander Zanden et al. [52] also suggest that the need to use year-specific δ 2 H f isoscapes would depend on the study region because the degree of inter-annual variation in δ 2 H p differs among regions. It is also likely that we are close to identifying the limit of accuracy of this single isotope approach, highlighting the importance of combining multiple markers in multivariate assignment tests [59,[61][62][63].
Finally, there are a number of caveats that need to be considered in our study and future research is encouraged to investigate further potential mechanisms contributing to within-site variance in δ 2 H f values. For example, it is still unclear how local hydrology interacts with biology by influencing relationships between δ 2 H values of foodwebs used by birds during the molt period and their subsequent δ 2 H f values. At interior continental sites, there are extreme seasonal effects in precipitation δ 2 H values [64] and northern sites will also have potential contributions from snowmelt as well as growing-season precipitation. These factors can strongly influence resulting assumed environmental δ 2 H values experienced by molting birds. In northern regions where H from snowmelt contributes to the overall foodweb leading to bird feathers, the use of an amount-weighted mean annual precipitation δ 2 H value or an average corresponding to those months where δ 2 H from precipitation has the strongest influence on the foodweb may be more appropriate. Nonetheless, our study provides encouraging evidence that once isotopic variance is understood and accounted for, current assignment isotopic models provide a valuable means of propagating such variance into the most parsimonious depiction of origins for North American passerines.