Exploring Overlap of Feather Molting and Migration in Tundra Swans Using δ2H Analysis

Abstract Determining the processes that shape the relative timing of energetically-costly events in the annual cycle of migrating birds is important to our understanding of avian phenology and ecology. We paired satellite tracking and hydrogen stable isotope analysis (δ2H) to examine the relative timing of two such events – migration and feather molting – in tundra swans from four breeding areas in Alaska, USA. Our results show a trend of increasing intra-individual variability in breast feather δ2H values with increasing migration distance, suggesting the overlap of breast feather molting and migration. However, when individual samples were pooled by breeding area, the δ2H values of breast and head feathers showed no trend with migration distance, presumably resulting from high levels of inter-individual variability in δ2H values within each breeding area. We explore potential reasons for this variability, propose potential mechanisms influencing feather δ2H values of tundra swans, and recommend further research into methods for exploring the temporal configuration of events in the annual cycle of migrating birds.


Introduction
Feather molting is an energetically expensive process, and molt timing in some species may have evolved to minimize overlap with other metabolically demanding events, such as reproduction and migration [1][2][3][4]. While this 'staggered cost' hypothesis has received considerable attention, several species of waterfowl have been identified with some degree of overlap between molting and other energetically costly events [1][2][3], implying that energetic constraints may not be the sole determinant in structuring the timing of molting [1]. By identifying potential non-energetic determinants of molt timing, we can develop a more thorough understanding of avian ecology and the processes that shape the annual cycles of birds.
Tundra swans (Cygnus columbianus) are a species for which overlap of feather molting with other energetically costly events has been documented. Tundra swan migrations are generally reported to occur in the autumn between September and December and again in the spring between February and May. Body feathers are molted between June and December, and primary feathers are molted during a much narrower period between July and September [5]. Imbrication in the reported general timing of migration and body feather molt should not; however, be interpreted as absolute overlap of these two events as substantial variability in the timing, durations, and patterns of both molt and migration have been observed among individuals and breeding populations [2,3,6,7]. In one of the earliest reports of the co-occurrence of these two events, Earnst [2] observed overlap between primary feather molt and breeding/pre-migration fattening. Earnst [2] suggested that the short breeding season and the spatially and temporally extensive migrations undertaken by tundra swans in which individuals can spend many weeks and travel hundreds of miles between breeding and wintering areas [7] may have selected for the overlap of these events. Similarly, Craigie and Petrie [3] observed evidence of contour feather molt during both autumn and, to a lesser extent, spring migrations in tundra swans. The birds examined in Craigie and Petrie's [3] work spend over half of the annual cycle in migratory staging areas [3,5], leading the authors to a similar conclusion as Earnst [2], that temporal constraints in the annual cycle may result in overlap of molt and migration. This concept of temporally-influenced overlap between molting and migration leads to the hypothesis that tundra swans with shorter migrations will exhibit correspondingly less overlap of feather molting and migration than populations with longer migration distances. However, testing this hypothesis using traditional methods, such as feather observations on migrating birds, presents logistical challenges to generating statistically representative sample numbers.
Hydrogen stable isotope analysis (δ 2 H) provides a potential method by which we can investigate the connection between feather molting and migration in tundra swans in a relatively practicable fashion. Because feathers incorporate the δ 2 H values of resources consumed during growth [8][9][10][11][12][13], the geographic locations of feather growth can be explored by comparing the δ 2 H values of feathers to known patterns in the δ 2 H values of precipitation over continents [14]. For the North American continent, the δ 2 H values of precipitation range from ~ -30‰ in the southeastern United States to ~ -200‰ in northern Alaska [15][16][17]. We would expect the δ 2 H values of feathers that are molted and regrown using resources acquired in a single location (e.g. a breeding area) to cluster around the δ 2 H values of precipitation for that area. Conversely, we would expect the δ 2 H values of feathers that are molted and regrown using resources acquired in multiple locations with different precipitation δ 2 H values (e.g. multiple stopover sites during migration) to demonstrate greater overall variability (Figure 1).
North American tundra swans are income migrators [18] that embark on migrations with limited nutrient reserves [6,18,19] and fuel migration by consuming resources at feeding sites en route [20]. Consequently, we expect that: a) contour feathers grown using resources consumed during migration will have δ 2 H values that reflect those of the precipitation at migration stopover sites, b) the δ 2 H of contour feathers grown during migration will differ from those grown prior to the initiation of migration, and c) the δ 2 H values of all contour feathers from tundra swans that engage in migration and molt simultaneously will show more variability than those from tundra swans that complete molt before migrating. If time constraints in the annual cycle are a driver of the overlap between molting and migration; we would (assuming minimal variability in dietary ecology among individuals) expect a positive correlation between migration distance/ time spent migrating and variability in the δ 2 H values of contour feathers.
Here, we present the results of a study that paired satellite tracking with δ 2 H analysis of head, breast and wing feathers from tundra swans from four distinct breeding areas within Alaska, USA to evaluate potential overlap in feather molt and migration. For each breeding area, we assessed the absolute distances of fall and spring migrations, feather δ 2 H values, and potential relationships between migration distance and intra-and inter-population variability in feather δ 2 H values towards understanding the potential for non-energetic influences of molt timing in tundra swans.

Methods
Ethical approval: The research related to animal use has been complied with all the relevant national regulations and institutional policies for the care and use of animals.

Satellite Transmitters and Feather Sample Collection
Tundra swans were captured when undergoing wing molt using hand-held nets and implanted with satellite transmitters (Platform Transmitting Terminals (PTTs), Microwave Telemetry, Columbia, MD, USA) from mid-July through mid-August of 2008. The birds were captured at four breeding areas: 1) Northern Alaska Peninsula near King Salmon, 2) Yukon-Kuskokwim Delta near the town of Bethel and close to the coast, 3) Kotzebue Sound region, and 4) the Colville River Delta (Figure 2, coordinates for the breeding areas are provided in Table 1). After implanting satellite transmitters, 20-30 crown feathers and 6-10 breast feathers were collected from each swan, the age (SYsecond year, or ASY -after second year) of each bird was determined by plumage color, and the sex was assessed by cloacal examination. PTT implantation, transmission, and data curation methods are described in Ely and Meixell [7]. Briefly, during pre-migration (20 Jul-15 Sep) and winter (2 Dec-15 Mar) PTTs transmitted for 5 hours within each 72 -hour interval. During autumn (16 Sep-1 Dec) and spring (18 Mar-15 May) migration periods, PTTs transmitted for 5 hours within each 23 -hour interval. After 20 July 2010, PTTs transmitted for 5 hours within each 80 -hour period until transmitters failed. Position location, including estimates of latitude and longitude, date, time, and quality of location, was obtained from the Argos Data Collection and Location System (2007; CLS America, Woods Hole Group, Lanham, MD, USA). Unlikely locations were filtered based on rate and angle of movement [21]. The highest quality location estimate was used to represent the daily position of a swan when multiple daily locations were received during a transmission cycle. If a PTT failed or a swan died, data were censored for that migration period. Migration distance was calculated by summing the total distance moved between sequential locations during a migratory period [21,22].  Table 1: Breeding ground locations, number of individual birds included in statistical analyses by feather type, mean migration (spring + autumn / 2) distance of each population, and total ranges of δ 2 H values of precipitation encountered over each migration route.

Sample Preparation and Stable Isotope Analyses
Feathers were cleaned with a 2:1 chloroform:methanol solution to remove surface contaminants, rinsed with deionized water, and air-dried. Feather barbs were removed from the rachis and minced into small sections using surgical scissors. Approximately 0.1 -0.2 mg of each minced feather barb sample was loaded into a silver capsule (Costech Analytical Technologies, Valencia, CA, USA) for δ 2 H analysis. Carbon ( 13 C) and nitrogen ( 15 N) stable isotope values of feathers were determined in order to examine dietary heterogeneity [23] as a potential source of variability in feather δ 2 H values. For δ 13 C and δ 15 N analysis, 1 mg of each sample was loaded into a tin capsule (Costech Analytical Technologies, Valencia, CA, USA). Because of their small mass, head feathers from the same individual had to be pooled in order to create samples with a sufficient mass for reliable stable isotope analyses. Breast feathers, which are substantially larger than head feathers, were not pooled for analysis. Prior to δ 2 H analysis, processed feather samples were benchtop equilibrated with ambient water vapor [24] for a minimum of two weeks to account for exchangeable hydrogen. Stable isotope analyses were performed at the University of Alaska Anchorage stable isotope facility using a Finnigan TCEA (Thermo Fisher Scientific Inc, Waltham, MA., USA) coupled to a ThermoFinnigan Delta Plus XP mass spectrometer (Thermo Fisher Scientific Inc, Waltham, MA., USA) and a Costech ECS 4010 elemental analyzer (Costech Analytical, Valencia, CA., USA) coupled to a Thermo-Finnigan Delta V Advantage mass spectrometer (Thermo Fisher Scientific Inc, Waltham, MA., USA) for δ 2 H and δ 15 N/δ 13 C analyses, respectively. All δ 2 H analyses were conducted using three keratin standards of known non-exchangeable hydrogen composition (CCHIX and TURK1 from the University of Georgia and BWBII from the University of Alaska Fairbanks; mean ± sd (‰) = -93.4 ± 1.3, -53.9 ± 1.4, and -108.0 ± 2.1, respectively). The standards were kept in a well plate exposed to laboratory air alongside the feather samples to be analyzed. A regression equation was created using raw versus known δ 2 H values of the keratin standards for each run, which allowed for the correction of any bias due to exchangeable hydrogen. Regression R 2 values were consistently 0.95 or higher. Long term analysis of quality control standards yielded precision of 2.5‰ for δ 2 H, 0.1‰ for δ 13 C, and 0.2‰ for δ 15 N.
During an independent set of collection trips in late July, 2010, actively-growing wing covert feathers from separate individuals were sampled at each of the five sites using the methods described above. Wing covert feathers were stored, processed, and analyzed using the methods described above.

Statistical Analysis
Migration routes and stopover sites for swans from the four breeding areas were determined using PTT data [7]. δ 2 H precipitation values for the endpoints and stopover sites for autumn (September, October, and November) and spring (March, April, and May) migration routes were estimated using the Online Isotopes in Precipitation Calculator [25].
Bayesian models were used to estimate mean δ 2 H, δ 15 N, and δ 13 C feather values among breeding areas and associated variability. By using a Bayesian approach, we were able to estimate the uncertainty of our variability estimates via posterior distributions. For head and wings feathers, which only had one sample per individual, isotope values were modelled using breeding area as a fixed effect with breeding area-specific standard deviations. For the breast feathers, where multiple samples were collected from each individual, we constructed hierarchical models to estimate breeding area-specific isotope values and to partition variance within breeding areas as well as within individuals. Non-informative priors were used for all parameters in all models. Summary statistics for posterior distributions are reported as the mode and highest posterior density 95% credibility intervals (CI). For each of the models, we ran three MCMC chains each consisting of 100,000 samples with the first 1,000 removed, and thinning by retaining every third sample. Convergence was confirmed by visual inspection of traceplots as well as ensuring the Rhat diagnostic was < 1.1 [26]. The models were run in JAGS (ver. 4.1.0) using R (v. 3.4.0) as an interface. Feathers that did not have all three stable isotope values (δ 2 H, δ 15 N, δ 13 C) were excluded from analyses, resulting in the removal of head feather samples from one bird from the Colville River Delta and one bird from the Yukon-Kuskokwim Delta (Table 1). In addition, analyses of head and breast feathers were restricted to ASY birds to avoid the potentially confounding effects of age class. Because wing covert feathers were sampled during active growth, we expected their isotopic values to represent resources consumed at the breeding area regardless of age class. Wing feather analyses were not; therefore, restricted by age.

Results
PTT data demonstrated that swans from the four breeding areas had migration routes that differed in length ( Figure  2, Table 1; [7]). Comparison of migration distances to time spent in breeding areas suggested that a) swans with longer migration distances generally spent less time in breeding areas, and b) swans from all four breeding areas migrated at a comparable rate, regardless of migration distance [7]. Birds from the Colville River Delta encountered the largest range of δ 2 H precipitation values during migration. Birds from the Kotzebue Sound, Yukon-Kuskokwim Delta, and King Salmon breeding areas encountered similar ranges of δ 2 H precipitation values, although the exact values of δ 2 H precipitation encountered differed among the sample populations ( Figure 3, Table 1). Ranges of δ 2 H precipitation values encountered were consistently larger during spring migration than fall migration.
Model results showed variability in the estimated means and standard deviations of δ 2 H values among breeding areas, with mean values ranging between -97.8‰ to -124.1‰ for head feathers and -138.8‰ to -173.3‰ for breast feathers and mean standard deviations ranging from 17.2‰ to 28.3‰ for head feathers and 0.0‰ to 18.9‰ for breast feathers (Figure 4, Table 2). In many cases, within-breeding area standard deviations were greater than the differences in δ 2 H values among breeding areas. No substantial differences in the means or standard deviations of δ 2 H values of head feathers were detected among the breeding areas, as demonstrated by the overlap of all 95% credibility intervals ( Figure 4, Table 2). For breast feathers, there were no differences in mean δ 2 H values among breeding areas, but the King Salmon population had a lower standard deviation than the Yukon-Kuskokwim Delta, Kotzebue, or Colville populations ( Figure 4, Table 2). Neither the means nor standard deviations of δ 2 H values of either feather type demonstrated a trend with migration distance.
Variability in the estimated means of δ 2 H values for wing covert feathers was generally lower than contour feathers. Overlap of 95% credibility intervals suggests no substantial differences in the means or standard deviations of δ 2 H values of wing covert feathers among the breeding areas ( Figure 4, Table 2). Neither the means or standard deviations of δ 2 H values of wing covert feathers demonstrated a trend with migration distance.
Ranges of δ 15 N and δ 13 C values suggest moderate variability in dietary constituents, but this variability was not consistent among breeding areas. Mean δ 15 N values of breast feathers were lower for the King Salmon population, but did not differ among the other populations. For head feathers, mean δ 15 N values were modestly higher for the Colville population, but did not differ among the other populations. No differences in the standard deviations of δ 15 N values were detected among the breeding areas for either feather type. No differences in means or standard deviations of δ 13 C values were observed among breeding areas for either feather type, and there was no relationship between means or standard deviations of δ 15 N and δ 13 C values and population migration distance for either feather type (Figure 4, Table 2). and autumn (open circles) precipitation δ 2 H values at stopover sites used during migration (estimated using the Online Isotopes in Precipitation Calculator [23]). Breeding areas, wintering areas, and stopover site locations were determined from PTT data. There was noticeable variation in both the length and in the range of δ 2 H values in spring and autumn precipitation encountered among migration routes.  Intra-individual variability in breast feather δ 2 H values showed a trend with migration distance, with individuals sampled from breeding areas associated with longer migration distances exhibiting greater variability than those from breeding areas with shorted migration distances ( Figure 4, Table 2). Overlap of 95% credibility intervals suggests no substantial differences in within-individual variability in δ 13 C and δ 15 N values of breast feathers among the breeding areas ( Figure 4, Table 2), and there was no relationship between intra-individual variability in δ 13 C and δ 15 N values of breast feathers and migration distance.

Discussion
The staggered cost hypothesis suggests separation of energetically costly annual events, such as molting and migration [1][2][3][4]. Our work investigated the application of δ 2 H analysis, a commonly-applied geographic assignment technique to assess this separation in a species for which temporal limitations may be of equal importance to energetic considerations in influencing the relative timing of molt and migration. We expected that birds that engage in complete separation of these two events would regrow all feathers in relatively constrained geographic spaces with limited variability in δ 2 H precipation , thereby resulting in consistency in the δ 2 H values of contour feathers. Conversely, we would expect complete intersection of body feather molting and migration to result in comparatively large variability in feather δ 2 H values. Under this paradigm, birds with comparatively longer migration distances might experience more overlap in body feather molting and migration, resulting in a relationship between variability in feather δ 2 H values and migration distance. Indeed, our results show a trend of increasing intra-individual variability in breast feather δ 2 H values with increasing migration distance, suggesting both overlap of breast feather molting and migration and a potential temporal determinant in the relative timing of the two events. While previously recognized in waterfowl [1,2,3], this overlap of molting and migration may occur in other bird taxa and has, as such, received considerable focus in the recent peer-reviewed literature [27 and references therein]. Additionally, while we do see substantial variability in δ 13 C or δ 15 N values, this variability does not reflect the consistent pattern observed in δ 2 H values, supporting the influence of migration, rather than differences in diet, in shaping intra-individual δ 2 H variability. Unfortunately, without a specific diet-to-tissue discrimination value for tundra swan feathers (∆ 2 H feather-preciptiation ), we are unable to a) determine specific locations contributing to feather growth, or b) identify specific areas of particular dietary importance to migrating swans (e.g. areas that provide for the energetic needs of migrating tundra swans at a disproportionately high level). We recommend further studies investigating ∆ 2 H precipitation -feather values and employing intra-individual δ 2 H analysis. Further, we recommend focusing intra-individual δ 2 H analyses on one feather type, as our results show inconsistency in the differences between head and breast feather δ 2 H values among individuals.
The relationship between variability in feather δ 2 H values and migration distance we observed for individuals was not apparent at the breeding area scale. This is likely the result of inter-individual variability in feather δ 2 H values within breeding areas caused by differences in habitat use, dietary ecology, and/ or physiology. For example, Ely and Meixell [7] noted that birds from the King Salmon breeding area (termed "Bristol Bay Lowlands" in Ely and Meixell [7]) make use of two wintering locationsone in northern California (the same site as the Kotzebue birds) and the other in the upper Pacific Northwest ( Figure  1). Some birds even alternate between the two locations [7]. Indeed, when we examine the distributions of individual breast feather δ 2 H values within breeding areas, the Colville, Kotzebue, and Yukon-Kuskokwim Delta breeding areas display multimodality, whereas the distribution for King Salmon is decidedly unimodal ( Figure S1). While the modalities observed for Colville, Kotzebue, and Yukon-Kuskokwim Delta may be interpreted as suggesting breast feather molting and growth at multiple locations, the unimodality associated with King Salmon is a) likely the result of a homogenizing effect consequent to introducing a third wintering area, and b) demonstrative of the influence of inter-individual variability in habitat use when attempting to draw conclusions at the breeding area scale. In addition to inter-individual variability in habitat use, subtle variations in dietary ecology, and/ or physiology may also occlude patterns revealed by intra-individual analyses. Previous studies examining the relationship between the δ 2 H values of precipitation and feathers have noted large variability in ∆ 2 H feather-precipiation values of up to 34‰ among individuals, populations, and species [28] caused by both ecological influences, such as diet, season, altitude, etc. and physiological processes, such as differences in fractionation during feather growth or metabolic routing [29][30][31]. Wolf et al. [31] reported standard deviations in δ 2 H values of up to 28.9‰ in birds maintained on the same diet and drinking water sources simply as a result of individual differences in drinking water consumption rates, condition, and physiology. In Centre (ASC). This work was conducted under Federal Permit #MB789758 (2007-2010), and ASC IACUC # 2008- 15 (2009-2010). Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. The authors declare that they have no conflict of interest. This research was supported in part by the USGS and U.S. Fish and Wildlife Service. The analysis was supported by a NSF MRI award (0953271) to JMW and we greatly appreciate the assistance of N Bargmann and B Cohn (UAA) who assisted in the preparation of feather samples for analysis. Supporting data are available in Ely et al. [33] and Wolf and Ely [34] our study, inter-individual variation in δ 2 H values resulting from these phenomena resulted in high relative variability within breeding areas compared to the variability observed both among the breeding areas and in the total range in precipitation δ 2 H values encountered during migration, functionally washing out the patterns identified by the intra-individual analyses. This strong influence of inter-individual variation is apparent in the variability we observed in wing covert feathers. Given that the wing covert feathers we sampled were actively growing, and that wing covert feathers must be completely regrown before migration occurs, we would expect relatively minor variability in the δ 2 H values of these feathers within breeding areas. Indeed, wing covert feathers do display the lowest amount of overall variability of the three feather types investigated in our study, but we still observed variability of up to 14.5‰ within breeding areas. Because of the substantial variability in δ 15 N and δ 13 C values within and among feather types and breeding areas, we are unable to partition the influences of physiological or behavioral (e.g. drinking water consumption rates) factors and differences in diet among individuals.
Because tundra swan populations have long migration routes and must dedicate significant portions of the annual cycle to travelling [20], understanding the mechanisms that influence the relative timing of energetically-costly events, such as molting and migration, is of great ecological importance, especially when considering potential environmental change. Changes to bird migratory ecology as a result of anthropogenically-induced modifications to migration corridors, habitat degradation, or adjustments in forage quality could result in misalignment in the timing of arrival to breeding or wintering areas, leading to potential detrimental effects on fecundity [20,32]. Our study used hydrogen stable isotope analysis to demonstrate overlap in breast feather molting and migration in individual birds. Unfortunately, inter-individual variability in hydrogen stable isotope analysis confounds the use of this technique to study possible overlap in molting and migration on larger scales. Given the importance of this phenomena; however, we recommend further study into the mechanisms influencing the relative timing of these two events.