A multi-isotope (δ13C, δ15N, δ34S, δ2H) approach to establishing migratory connectivity in lesser snow geese: Tracking an overabundant species

Abstract

Expanding populations of North American midcontinent lesser snow geese (Anser caerulescens caerulescens) have potential to alter ecosystems throughout the Arctic and subarctic where they breed. Efforts to understand origins of harvested lesser snow geese to better inform management decisions have traditionally required mark-recapture approaches, while aerial photographic surveys have typically been used to identify breeding distributions. As a potential alternative, isotopic patterns that are metabolically fixed within newly grown flight feathers following summer molting could provide inferences regarding geographic breeding origin of individuals, without the need for prior capture. Our objective was to assess potential to use four stable isotopes (δ13C, δ15N, δ34S, δ2H) from feather material to determine breeding origins. We obtained newly grown flight feathers from individuals during summer banding at three Arctic and two subarctic breeding colonies in 2014 (n = 56) and 2016 (n = 45). We used linear discriminant analyses to predict breeding origins from models using combinations of stable isotopes as predictors and evaluated model accuracy when predicting colony, subregion, or subpopulation levels. We found a strong inverse relationship between δ2H values and increasing latitude (R2 = 0.83), resulting in differences (F4, 51 = 90.41, P < 0.0001) among sampled colonies. No differences in δ13C or δ15N were detected among colonies, although δ34S in Akimiski Island, Baffin Island, and Karrak Lake were more enriched (F4, 51 = 11.25, P < 0.0001). Using δ2H values as a predictor, discriminant analyses improved accuracy in classification level as precision decreased [model accuracy = 67% (colony), 88% (subregion), 94% (subpopulation)]. Application of the isotopic methods we describe could be used to provide an alternative monitoring method of population metrics, such as overall breeding population distribution, region-specific productivity and migratory connectivity that are informative to management decision makers and provide insight into cross-seasonal effects that may influence migratory behavior.

About the Authors:

Drew N. Fowler

Contributed equally to this work with: Drew N. Fowler, Elisabeth B. Webb, Frank B. Baldwin, Mark P. Vrtiska, Keith A. Hobson

Roles Conceptualization, Formal analysis, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Missouri Cooperative Fish and Wildlife Research Unit, School of Natural Resources, University of Missouri, Columbia, Missouri, United States of America

ORCID logo http://orcid.org/0000-0001-9347-4579

Elisabeth B. Webb

Contributed equally to this work with: Drew N. Fowler, Elisabeth B. Webb, Frank B. Baldwin, Mark P. Vrtiska, Keith A. Hobson

Roles Conceptualization, Project administration, Supervision, Writing – review & editing

Affiliation: U. S. Geological Survey, Missouri Cooperative Fisheries and Wildlife Research Unit, School of Natural Resources, University of Missouri, Columbia, Missouri, United States of America

Frank B. Baldwin

Contributed equally to this work with: Drew N. Fowler, Elisabeth B. Webb, Frank B. Baldwin, Mark P. Vrtiska, Keith A. Hobson

Roles Conceptualization, Writing – review & editing

Affiliation: Department of Sustainable Development, Province of Manitoba, Winnipeg, Manitoba, Canada

Mark P. Vrtiska

Contributed equally to this work with: Drew N. Fowler, Elisabeth B. Webb, Frank B. Baldwin, Mark P. Vrtiska, Keith A. Hobson

Roles Project administration, Supervision, Writing – review & editing

Affiliation: Nebraska Game and Parks Commission, Lincoln, Nebraska, United States of America

Keith A. Hobson

Contributed equally to this work with: Drew N. Fowler, Elisabeth B. Webb, Frank B. Baldwin, Mark P. Vrtiska, Keith A. Hobson

Roles Formal analysis, Methodology, Writing – review & editing

Affiliation: Department of Biology, Western University, London, Ontario, Canada

Introduction

North American midcontinent lesser snow geese (Anser caerulescens caerulescens) are long-distant migrants that predominately winter in temperate Arkansas and coastal Texas and Louisiana along the Gulf of Mexico [1–2]. Breeding distribution of midcontinent lesser snow geese ranges from 52° N to 85° N and spans an approximate 33° range in longitude. The extensive breeding range of lesser snow goose has resulted in numerous dispersed breeding colonies and the classification of northern and southern subpopulations (henceforth referred to as Arctic and subarctic subpopulations). Midcontinent lesser snow geese have been considered an overabundant species since 1999, based on explosive population growth facilitated by increasing agricultural production [3–4]. Efforts to understand the life history of lesser snow geese and their response to population reduction actions have been multifaceted; however, management actions taken to reduce the population have been deemed largely unsuccessful [5–6]. As the midcontinent population persists above desired levels, continued monitoring of metrics such as productivity, immigration / emigration rates among colonies, links between breeding and non-breeding areas (i.e., migratory connectivity) [7–8] and harvest demographics are important for tracking population change and informing future management decisions [5].

Estimation of age-specific annual survival and harvest rates among colony or subpopulation groups are important monitoring metrics and rely on long-term, annual banding efforts, in conjunction with band-recovery analysis [9]. To determine subpopulation affiliation, capture and marking of individuals must occur on the breeding grounds. Due to the large geographic breeding range of the population, general inaccessibility of the region, and relatively short molting period (i.e., when mass capture is possible), banding occurs at 5 distinct regions in the Arctic and subarctic, and indices of productivity (age ratios) at this time are assumed to be representative of the mid-continent population [10]. As an alternative, differentiation in stable isotope ratios across the wide geographic breeding range of lesser snow geese may provide a method to identify subpopulation or colony association of individuals of unknown origin, and provide annual indices of region-specific productivity that are likely less biased than those currently derived from the banded sample. Such a technique could also enable more frequent and cost-efficient assessments of breeding distribution, which have previously only been available through aerial photographic surveys.

Stable isotopes of light elements such as carbon (δ13C), nitrogen (δ15N), oxygen (δ18O), hydrogen (δ2H), and sulfur (δ34S) can serve as endogenous markers that allow inference to prior life history locations without requiring previous capture and recapture [11–12]. Briefly, food-webs differ regionally in isotopic signatures of light elements, and consumers incorporate and reflect tissue isotopic values specific to the food-web they inhabit [13]. These isotopic patterns can become fixed in metabolically inactive tissues such as feathers, claws, or hair [14] or exist in dynamic equilibrium in metabolically active tissues (e.g. blood, muscle) of an organism that moves among isotopically different food-webs [13, 15–16]. Therefore, isotopic analysis of select animal tissues has the capacity to reveal information on prior discrete events at both spatial and temporal scales.

Lesser snow geese undergo a complete flight feather molt in mid-summer, usually after breeding. Therefore, values of stable isotopes of newly grown flight feathers are representative of local source materials consumed at the time of flight feather growth (but see Fox et al. [17]). If individual breeding colonies produce unique feather isotopic signatures, then data conferring an individual’s previous year’s breeding / molt location (collected via harvest) could provide an alternative means of elucidating migratory connectivity and metrics of spatial population productivity. For example, Hobson et al. [18] used δ2H from flight feathers of hunter-harvested lesser scaup (Aytha affinis) to identify variation in breeding productivity of unique geographic regions. Similarly, isotopic analysis of muscle tissue has proven particularly useful in lesser snow geese to identify geographic wintering locations [19]. Consequently, joint application of these methods could provide insight into cross-seasonal effects that may exist within or among species that influence migratory behavior [20–21]. For instance, Paxton and Moore [22] assessed δ13C and δ2H values of red blood cells and feathers, respectively, from black-and-white warblers (Mniotilta varia) and found that quality of available winter habitats established differential body condition and proximate cues for spring migration and influenced individual migration movements.

To test the utility of stable isotope ratios in predicting geographic breeding origin of lesser snow geese, we collected newly grown flight feathers of adult female lesser snow geese at three Arctic and two subarctic breeding colonies in 2014. We evaluated stable isotope ratios expected to vary across latitude (δ2H) [23], diet (δ13C and δ15N) [24] and proximity to marine coastal environments (δ34S) [25] to provide for a wide potential variation in isotopic differentiation across colonies.

Materials and methods

Collection

To ensure certainty of isotope values originating from specific breeding locations, we collected newly grown flight feathers (P1 feather) of adult female lesser snow geese during routine mid-summer banding operations at three Arctic and two subarctic breeding colonies (Fig 1) in July 2014. Arctic breeding colonies included Southampton Island (63° 48' 30.67", -85° 41' 50.71") (n = 12), Baffin Island (66° 42' 46.04", -72° 33' 26.32") (n = 8) and Karrak Lake (67° 15' 37.30", -100° 16' 25.10") (n = 7). Subarctic breeding colonies included Akimiski Island (53° 06' 17.42", -80° 57' 28.73") (n = 11) and La Pérouse Bay (58° 43' 5.62", -93° 53' 21.54") (n = 18) (Fig 1). Upon removal, feathers were stored dry in paper envelopes until processing. All feathers from breeding colonies were collected in accordance with Environment Canada Animal Care Committee permit number #16JL01. Whole specimens during spring migration were collected under the University of Missouri’s Animal Care and Use Committee permit number 8191 and the United States Fish and Wildlife Service scientific collection permit number MB47969B-1.

[Figure omitted. See PDF.]

Fig 1. Location of sampled breeding colonies.

Newly grown flight feathers of adult female (n = 56) midcontinent lesser snow geese (Anser caerulescens caerulescens) were collected for quantification of isotope signatures in July 2014. A second year of collections containing flight feathers from both adult males and females (n = 45) occurred at all colonies except Karrak Lake in July 2016.

https://doi.org/10.1371/journal.pone.0203077.g001

Collected feathers were cleaned of surface oils in chloroform:methanol (2:1 v/v) solvent rinse. Stable hydrogen isotope measurements were performed on H2 gas derived from high-temperature flash pyrolysis of 350 ± 20 μg feather subsamples from the distal section of feather vanes (packed in silver capsules) using continuous-flow isotope-ratio mass spectrometry. Briefly, pyrolytic combustion (1350°C) was on a reactor with glassy carbon chips under helium flow in a Eurovector (Milan, Italy– www.eurovector.it) elemental analyzer interfaced with an Isoprime (Manchester, UK) continuous flow isotope ratio mass spectrometer (IRMS). Estimates of the nonexchangeable H in samples were derived from two keratin hydrogen-isotope reference materials (CNS: -197 ‰; KHS: -54.1 ‰) following the comparative equilibration approach of Wasssenaar and Hobson [26]. All δ2H results are reported for nonexchangeable H in delta notation, in units of per mil (‰), and normalized on the Vienna Standard Mean Ocean Water (VSMOW) standard scale. Based on within-run replicate (n = 8) measurements of keratin reference materials, we estimated measurement error (SD) to be ~ ± 2 ‰.

For δ13C and δ15N analyses, 0.5–1.0 mg of feather material (packed in tin capsules) was combusted online using a Eurovector 3000 elemental analyzer. The resulting CO2 was separated by gas chromatography (GC) and introduced into a Nu Horizon (Nu Instruments, Wrexham, UK– www.nu-ins.com) triple-collector isotope-ratio mass-spectrometer via an open split and compared to a pure CO2 or N2 reference gas. Stable nitrogen (15N/14N) and carbon (13C/12C) isotope ratios were expressed in δ notation, as parts per thousand (‰) deviation from the primary standards, atmospheric AIR and Vienna Pee Dee Belemnite (VPDB). Using previously calibrated internal laboratory standards (powdered keratin [BWB II: δ 13C = -20.0‰, δ15N = -14.1‰ and gelatin: δ 13C = -13.6‰, δ15N = -4.7‰]) within-run (n = 5), precision for δ15N and δ13C measurements was ~ ± 0.15‰.

For δ34S analyses, 3500 ± 100 ug of feather material (packed in tin capsules) was combusted in a Vario Pyro Cube (Elementar, Langenselbold, Germany– www.elementar.de) elemental analyzer and the resulting SO2 gas was introduced into Isoprime IRMS. Our laboratory standard was BWB-3 keratin (δ34S = 13.2‰) and δ values were reported relative to the Canyon Diablo Triolite (CDT) standard. Measurement precision, based on within-run replicate measurements of the lab standard was ± 0.3‰.

Isotopic segregation of breeding colonies

To determine whether feathers from different breeding areas were isotopically distinct, we calculated mean feather stable isotope values within each breeding colony and ran a multiple analysis of variance (MANOVA; α < 0.05). Feather values for all four isotopes were normally distributed (Shapiro test: all P’s > 0.05) and homoscedastic (Levene’s test: all P’s > 0.05). We used Tukey posthoc tests to determine differences among individual isotope means and specific breeding colonies. To describe spatial variation in feather stable isotope ratios across the geographic range of selected breeding colonies, we used simple linear regression analyses to relate stable isotope values to variation in latitude and longitude.

We used linear discriminant function analyses to develop a suite of candidate models to predict breeding colony association using combinations of δ2H, δ13C, δ15N, and δ34S values (Table 1). Based on MANOVA results, we developed a second set of linear discriminant models in which we grouped individuals from Southampton Island and Baffin Island together as a subregion based on their similar feather δ2H values (Table 1). Finally, we evaluated a set of linear discriminant models where colonies were compressed into either Arctic or subarctic subpopulation association (Table 1). We used K-fold cross validation (K = 5) to evaluate the accuracy of assignment for discriminant functions built to determine colony and subpopulation association using the “lda” function in the MASS package in Program R [27]. We determined the strength of individual models based on the percent accuracy of individuals from each breeding colony in 2014 whose isotopic value correctly represented their breeding colony [28].

[Figure omitted. See PDF.]

Table 1. Models selected for linear discriminant function analysis.

https://doi.org/10.1371/journal.pone.0203077.t001

As a proof of concept, we collected flight feathers (n = 40) from adult, female snow geese harvested in Arkansas, Missouri, Nebraska, and South Dakota during spring migration 2015. In a preliminary investigation, feathers from geese of unknown origins were prepared in identical manner described above and, based on development of discriminant modeling results, analyzed for δ2H values. Unknown individuals were then assigned origin based on the greatest posterior probability of belonging to one of the five colonies (and subpopulation association) predicted by a linear discriminant model based on δ2H using the “lda” “predict” function in the MASS package.

Inter-annual variability in δ2H

The use of feather δ2H values as a reliable marker of goose origin is dependent upon how well modeled food-web δ2H and corresponding feather δ2H (δ2Hf) matches predictions. Such models are ultimately based on the long-term dataset provided by the International Atomic Energy Association (IAEA) Global Network of Isotopes in Precipitation (GNIP). Predicted local precipitation values are derived from kriged data using actual amount-weighted monthly average precipitation δ2H (δ2Hp) values from a network of sampling stations. Those data are then typically transformed into expected δ2Hf values for spatially explicit assignment. Here we relied instead on describing expected feather isotope values based solely on year-specific, ground-truthed sampling of feathers. For δ2Hf, some interannual variation can be expected [12, 29–30] and such uncertainty is propagated in most spatially explicit probabilistic assignments [31]. To estimate potential interannual variation in lesser snow goose δ2Hf, we obtained additional lesser snow goose feathers from known colonies in 2016. Here, secondary covert flight feathers from adult male and female birds were collected at Akimiski Island (n = 9), La Pérouse Bay (n = 10), Baffin Island (n = 11) and Southampton Island (n = 15) in 2016 and analyzed for δ2Hf. We compared differences in δ2Hf by colony, year, and the interaction with a two-way analysis of variance (α = 0.05). We created a discriminant function model using δ2Hf values based on feathers collected in 2016 and used feathers derived from the same colonies in 2014 as a test group. Similarly, we used the δ2Hf discriminant function model based on 2014 feathers to obtain model predictions of the known 2016 feathers and compare differences in error rate between the two years. We used R 3.3.3 (R Core Development Team 2017) for all statistical analyses.

Results

Comparison of isotopes by breeding colonies

Simple linear regression analysis indicated a strong inverse relationship between δ2Hf values and increasing latitude (R2 = 0.83; Fig 2). No other isotope showed any detectable relationship between δ values and latitude. Minor relationships between δ34S and δ15N and longitude existed as depletion increased slightly in more westerly colonies, although variation in the data explained by the relationship was low (R2 = 0.12, R2 = 0.07, respectively; Fig 2).

[Figure omitted. See PDF.]

Fig 2. Relationships between isotopes (δ2H, δ34S, δ13C, and δ15N) analyzed from feather material to latitude and longitude for adult lesser snow geese collected July, 2014.

The top panel presents relationships between isotope values and latitude. The bottom panel presents relationships between isotope values and longitude. Breeding colonies are Akimiski Island (closed circle), La Pérouse Bay (open triangle), Southampton Island (closed triangles), Baffin Island (crosses), and Karrak Lake (closed squares). Linear regression lines are provided only for slopes that were significantly (P < 0.05) different from zero.

https://doi.org/10.1371/journal.pone.0203077.g002

Results from the MANOVA indicated overall differences (F16, 204 = 7.51, P < 0.0001) in measured isotope values among breeding colonies for feathers of known origin collected in 2014 (Table 2). However, only δ2Hf and δ34S differed among colonies (δ2Hf: F4, 51 = 90.41, P < 0.0001; δ34S: F4, 51 = 11.25, P < 0.0001) whereas δ13C and δ15N did not (δ13C: F4, 51 = 0.29, P = .8818; δ15N: F4, 51 = 2.02, P = 0.1047). Post hoc Tukey tests indicated δ2Hf differed separately among Akimiski Island, La Pérouse Bay, and Karrak Lake. Feather deuterium values did not differ between Southampton Island and Baffin Island but were different from remaining colonies (Table 2). Values of feather δ34S did not differ among Akimiski Island, Baffin Island, and Karrak Lake but as a group were more enriched than values found at La Pérouse Bay and Southampton Island (Table 2).

[Figure omitted. See PDF.]

Table 2. Mean and SE of stable hydrogen (δ2H), nitrogen (δ15N), carbon (δ13C), and sulfur (δ34S) isotope values in feathers collected from adult female lesser snow geese at Arctic and subarctic breeding colonies during summer banding in 2014a.

https://doi.org/10.1371/journal.pone.0203077.t002

Values of δ2Hf and δ34S most accurately predicted breeding colony association, when using K-fold cross-validation, compared to the comprehensive model using δ2Hf, δ13C, δ15N, and δ34S or δ2Hf alone (overall model accuracy of 80%, 76%, and 67%, respectively; Table 3). When predicting at the colony level, misclassification was greatest between Baffin Island and Southampton Island (Table 3). However, grouping Baffin Island and Southampton Island as a subregion improved overall model accuracies. Prediction for the model using δ2Hf alone was more accurate (88%) than the model collectively using δ2Hf and δ34S (86%; Table 3). Finally, greatest prediction accuracy was achieved when breeding colonies were classified according to subpopulation. The model using δ2Hf alone was sufficiently accurate to predict subpopulation association relative to the δ2Hf and δ34S model (94% accuracy compared with 93%; Table 3). Posterior classification probabilities using the δ2Hf model for 40 geese harvested during the 2015 spring migration increased as classification precision level decreased (S1 Table). Classifying individuals to either the subregion or subpopulation level incrementally reduced uncertainty in assignment probabilities without substantially changing the overall proportions of classifications to other groups (Fig 3).

[Figure omitted. See PDF.]

Fig 3. Breeding origin classification percentages from isotopic signatures of lesser snow geese (n = 40) harvested during spring migration, 2015.

https://doi.org/10.1371/journal.pone.0203077.g003

[Figure omitted. See PDF.]

Table 3. Prediction accuracies by discriminant function models, partitioned across the colony, subregion, and subpopulation level of adult lesser snow geese collected in 2014 using K-fold cross validation.

https://doi.org/10.1371/journal.pone.0203077.t003

Inter-annual variability in δ2Hf

Values of δ2Hf differed interactively among colonies and year (F3, 864 = 4.24, P = 0.0076) as only lesser snow geese collected at Akimiski Island indicated differences in δ2Hf values between 2014 and 2016 (Table 4). This annual difference in δ2Hf values at Akimiski Island was reflected by low Akimiski Island prediction accuracies (44%) when predictions of 2016 colony classification were made using the δ2Hf discriminant function derived from feathers collected in 2014 (Table 5). In contrast, δ2Hf discriminant functions derived from each respective year to predict subregion and subpopulation level increased overall prediction accuracies by 29% and 40%, respectively (Table 5).

[Figure omitted. See PDF.]

Table 4. Comparison of mean and SE of stable hydrogen (δ2H) isotope values in feathers collected from adult lesser snow geese in 2014 and 2016a.

https://doi.org/10.1371/journal.pone.0203077.t004

[Figure omitted. See PDF.]

Table 5. Comparison of prediction accuracies of annually derived discriminant function models, partitioned across the colony, subregion, and subpopulation level of adult lesser snow geese collected in 2014 and 2016, using K-fold cross-validation.

https://doi.org/10.1371/journal.pone.0203077.t005

Discussion

The primary goal of this study was to evaluate potential to determine breeding origins of lesser snow geese using a suite of stable isotopes analyzed from feather material and evaluate the combination of spatial precision and accuracy achievable under various origin scales. We found that among the four stable isotopes tested (δ2H, δ13C, δ15N, and δ34S), deuterium (δ2H) was the most distinct across the sampled breeding colonies. Generally, deuterium values were more depleted with increasing latitude. The positive relationship between depleted precipitation deuterium and increasing latitude, and the positive relationship between precipitation deuterium and feather deuterium is well established [23–31]. However, variance in the latter relationship can vary among taxa [31]. Our study quantified deuterium specifically in feather material of a large-bodied migratory bird whose breeding range extends across a wide Arctic and subarctic latitude. In contrast, δ13C and δ15N values did not vary across sampled breeding colonies and their inclusion in linear discriminant analyses reduced overall prediction accuracy slightly. Lack of difference in δ13C and δ15N among colonies may be related to similar plant foraging material across breeding areas [6].

We found differences in δ34S values among two groups of colonies. Values of δ34S were more enriched in colonies at Akimiski Island, Baffin Island, and Karrak Lake, in contrast to more depleted values at La Pérouse Bay and Southampton Island. In general, δ34S values are expected to be enriched in breeding areas in proximity to coastal areas where prevailing winds transfer sulfates from sea spray onto land [32–33]. Additionally, δ34S in feathers could be enriched based on the predominant consumption of marine foods over terrestrial resources. For lesser snow geese breeding at La Pérouse Bay, coastal marsh degradation from heavy grazing forced the colony to shift more inland towards freshwater marshes [34]. This inland shift and utilization of freshwater wetlands may explain the depleted δ34S values we observed in feathers from lesser snow geese at La Pérouse Bay in this study. However, lesser snow geese breeding at higher latitudes are thought to rely more heavily on freshwater wetlands than southern breeders [35], yet δ34S values from individuals from these areas in our study were more enriched. Nonetheless, δ34S could be useful as a secondary discriminant predictor of breeding colonies whose classification is confounded using only δ2Hf because of similar values among colonies. This was the case for classification of individuals from Baffin Island and Southampton Island. Using feather δ34S in addition to δ2Hf values improved overall model prediction accuracy from 67% to 80%. However, while improving prediction between Baffin Island and Southampton Island, using δ34S in conjunction with δ2Hf increased misclassification rates for individuals originating from Karrak Lake.

Selecting to use a discriminant function analysis derived at the colony, subregion, or subpopulation level results in a tradeoff between accuracy and precision. Our discriminant functions based on δ2Hf that classified individuals at the subregion or subpopulation level were 34% and 40% more accurate compared to a model classifying individuals to the colony level; however, the spatial scope of inference was reduced. Using a discriminate function based on δ2Hf that classifies an individual to the subregion level may currently be the most appropriate balance between accuracy of classification and inferential spatial scope.

The ability to isotopically distinguish among origin of individuals, even at the scale of subregion, has several useful applications, but the most important may be an improved understanding of breeding distribution. Overall population size can be estimated from age-specific harvest rate estimates using Lincoln’s method [36] obtained from banding data and age-specific harvest estimates, which are obtained from long-established hunter harvest surveys in Canada and the USA wherein a portion of hunters submit tail and wing feathers. Lincoln abundance estimates have become an important metric for monitoring population status for some populations of North American geese [10]. The components used in Lincoln estimates are available over long time periods and are thought to be more reliable than other survey methods (e.g., mid-winter aerial surveys, aerial photo surveys) that may be biased to an unknown extent by issues related to detection, speciation, and sampling coverage. Although they account for an unknown proportion of the total population and are conducted infrequently, photo-surveys of known nesting colonies have been used to improve the understanding of the breeding distribution of this population [37]. Lincoln estimates have become the preferred metric for monitoring changes in overall abundance to avoid the need for costly and logistically challenging photo-surveys [10]. However, Lincoln estimates do not provide information about the distribution of the breeding population, thus identifying an alternative approach to estimating breeding distribution would be an important contribution to population management and allow for monitoring migratory connectivity. Isotopic classification of clipped primary feathers submitted by hunters through hunter harvest surveys would permit partitioning of Lincoln estimates, at least to the level of subpopulation. Our classification of forty adult migrants collected during spring migration illustrates the potential to provide an alternative index that describes spatial distribution of the population. An improved understanding of breeding distribution is of management interest because the relationship between adult survival and harvest differs between northern and southern nesting segments of the population [5]. Thus, course-scale population distribution should be a consideration when designing monitoring programs, or when determining where to establish or expand research on population-level management questions.

Another important monitoring metric available from annual banding programs is the development of a productivity index (i.e., age ratio at banding). However, the degree to which age ratios at banding reflect overall productivity could be influenced by local weather or predation pressure [38], age of young at capture (e.g., exposure days) [39], and individual decisions made by banders (e.g., targeting smaller brood flocks or large groups of non-breeders, removal of incomplete captures from calculations, focusing on the best areas of production). Determining isotopic values of a random sample of primary feathers submitted by hunters through the annual harvest survey would likely yield a more spatially representative index of annual production than can be provided by annual banding programs.

Predicting breeding origins of lesser snow geese sampled away from breeding colonies based on a discriminant analysis derived from δ2Hf values of feathers grown at known colonies constrains classifications only to those colonies where feathers have been analyzed and included in the analysis. One drawback to our approach is the potential to misclassify an individual to one of the known colonies when true origin was another colony not assessed in this study. Similarly, we primarily analyzed feathers from individuals who successfully bred and were molting in July. Because of logistical constraints, banding efforts at Southampton Island only target failed breeders or those who forego breeding altogether. Therefore, our models are primarily derived from one cohort of individuals, breeders, and do not account for individuals who initiate a subsequent migration elsewhere to molt and regrow feathers. Extensive detail regarding the destination of molt migrants from specific lesser snow goose colonies is limited. In one study on lesser snow geese at La Pérouse Bay, molt migrants moved north to the McConnell River, about 250 km away [40]. Additional research identifying additional molt migrant destinations from other colonies would assist ability to interpret isotopic feather values. Additionally, we used only adult feathers to derive our linear discriminant functions. Previous works in other avian guilds (passerines and raptors) have identified more depleted deuterium values in juveniles compared with adults from the same known origin, potentially thought to be influenced by differences in diet, feather growth rate and metabolism [41–43]. While we do not have reasons to suspect differences in diet among juvenile and adult lesser snow geese, evaluating differences in δ2Hf among age classes will be important to ensure predictive accuracy. We also observed some inter-annual variation of δ2Hf values at Akimiski Island between sampling in 2014 and 2016 that resulted in misclassification when using one year as a training set for a model and the other as a test set. These limitations emphasize the importance of continued multi-year sampling (of newly grown feathers during molt) and an increase in the number of colonies sampled. Data from additional sampling could be used to calibrate a functional based map of δ2Hf across the Canadian Arctic and result in increased accuracy and precision of lesser snow goose breeding origin classification using stable isotope analysis.

Conclusions

Using isotopic markers found in metabolically active or inactive animal material provides a pathway to identify heterogeneous breeding origins within harvested populations without requiring previous capture, marking, and recapture. Our study provides an initial demonstration of the utility of identifying migratory connectivity through stable isotopes in lesser snow geese. Yet, our approach could be further developed through the calibration of the long-term GNIP dataset from lesser snow goose δ2Hf values, resulting in a species-specific breeding isoscape. The application of these data could be used to provide a more robust annual estimate of subregion or subpopulation size based on isotopic partitioning of harvest survey parts. As lesser snow geese have increased in abundance, juvenile-to-adult age ratios have declined [44] in response to increasing density-dependence [45]. Using previous techniques demonstrated to assess productivity [18] and migratory connectivity in harvested gamebirds [46] could be implemented in midcontinent lesser snow geese and facilitate comprehensive monitoring informative to future management decisions.

Supporting information

[Figure omitted. See PDF.]

S1 Dataset. Stable isotope values of hydrogen (δ2H), nitrogen (δ15N), carbon (δ13C), and sulfur (δ34S) in feathers collected from adult female lesser snow geese (n = 99) at Arctic and subarctic breeding colonies during summer banding in July, 2014 and 2016.

Stable isotope values are reported in per mil notation (‰). δ15N , δ13C, and δ34S were not analyzed in feathers collected in 2016.

https://doi.org/10.1371/journal.pone.0203077.s001

(DOCX)

S2 Dataset. Stable isotope values of hydrogen (δ2H) in feathers collected from adult female lesser snow geese (n = 40) collected during spring migration in 2015.

Stable isotope values are reported in per mil notation (‰).

https://doi.org/10.1371/journal.pone.0203077.s002

(DOCX)

S1 Table. Posterior probabilities of colony, subregion, and subpopulation association for lesser snow geese (n = 40) harvested during spring migration, 2015, based on linear discriminant analysis using δ2Hf from feathers collected on known breeding colonies in 2014 as a discriminating variable.

Percentages in bold identify the greatest probability of colony and subpopulation association for each collected specimen.

https://doi.org/10.1371/journal.pone.0203077.s003

(PDF)

Acknowledgments

Funding for this project was provided by the U.S. Fish and Wildlife Service’s Sportfish and Wildlife Restoration Program through the Nebraska Game and Parks Commission and by Ducks Unlimited Canada’s Institute for Wetland and Waterfowl Research. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. We thank R. Alisauskas, R. Brook, D. Kellett, J. Leafloor, R. Rockwell and others who collected feathers on breeding colonies for this project. Chantel Gryba prepared samples for isotope analysis and Geoff Koehler performed CNS isotope analyses at the Environment Canada Stable Isotope Laboratory in Saskatoon, SK. The Missouri Cooperative Fish and Wildlife Research Unit is jointly sponsored by the Missouri Department of Conservation, the University of Missouri, the U.S. Fish and Wildlife Service, the U.S. Geological Survey, and the Wildlife Management Institute. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Citation: Fowler DN, Webb EB, Baldwin FB, Vrtiska MP, Hobson KA (2018) A multi-isotope (δ13C, δ15N, δ34S, δ2H) approach to establishing migratory connectivity in lesser snow geese: Tracking an overabundant species. PLoS ONE 13(8): e0203077. https://doi.org/10.1371/journal.pone.0203077

References

1. Hobaugh WC. Body condition and nutrition of snow geese wintering in southeastern Texas. J Wildl Manage. 1985;49(4):1028–1037.

2. Jónsson JE, Afton AD. Are wintering areas shifting north? Learning from lesser snow geese banded in southwest Louisiana. Southeast Nat. 2015;14(2):293–307.

3. Batt B. Arctic Ecosystems in Peril: Report of the Arctic Goose Habitat Working Group. U.S Fish and Wildlife Service, Washington, D.C. and Canadian Wildlife Service, Ottawa, Ontario; 1997.

4. Abraham KF, Jefferies RL, Alisauskas RT. The dynamics of landscape change and snow geese in mid-continent North America. Glob Chang Biol. 2005 Jun;11(6):841–55.

5. Alisauskas RT, Rockwell RF, Dufour KW, Cooch EG, Zimmerman G, Drake KL, et al. Harvest, survival, and abundance of midcontinent lesser snow geese relative to population reduction efforts. Wildl Monogr. 2011;179(1):1–42.

6. Leafloor JO, Moser TJ, Batt B. Evaluation of special management measures for midcontinent lesser snow geese and Ross’s geese. Arctic Goose Joint Venture Special Publication. U.S. Fish and Wildlife Service, Washington, D.C. and Canadian Wildlife Service, Ottawa, Ontario; 2012.

7. Webster MS, Marra PP, Haig SM, Bensch S, Holmes RT. Links between worlds: unraveling migratory connectivity. Trends Ecol Evol. 2002;17(2):76–83.

8. Webster MS, Marra PP. The importance of understanding migratory connectivity and seasonal interactions. Birds of two worlds: the ecology and evolution of migration. Johns Hopkins University Press Baltimore, Maryland; 2005. 199–209 p.

9. Calvert AM, Alisauskas RT, White GC. Annual survival and seasonal hunting mortality of midcontinent snow geese. J Wildl Manage. 2017;81(6):1009–1020.

10. Central Flyway Council. Management guidelines for Central Flyway Arctic Nesting Canada geese. U.S. Fish and Wildlife Service, Lakewood, Colorado; 2013

11. Hobson KA, Wassenaar LI. Tracking animal migration with stable isotopes. Terrestial Ecology Series. Vol. 2. London, United Kingdom: Academic Press; 2008.

12. Hobson KA. Isotopic ornithology: A perspective. J Ornithol. 2011;152(1 SUPPL):S49–66.

13. Inger R, Bearhop S. Applications of stable isotope analyses to avian ecology. Ibis. 2008;150(3):447–461.

14. Hobson KA, Koehler G. On the use of stable oxygen isotope (δ18 O) measurements for tracking avian movements in North America. Ecol Evol. 2015;5(3):799–806. pmid:25691999

15. Hobson KA, Clark RG. Assessing avian diets using stable isotopes I: turnover of C in tissues. Condor. 1992;94(1):181–188.

16. Carleton SA, Del Rio CM. The effect of cold-induced increased metabolic rate on the rate of 13C and 15N incorporation in house sparrows (Passer domesticus). Oecologia. 2005;144(2):226–232. pmid:15891840

17. Fox AD, Hobson KA, Kahlert J. Isotopic evidence for endogenous protein contributions to greylag goose (Anser anser) flight feathers. J Avian Biol. 2009;40(2):108–112.

18. Hobson KA, Wunder MB, Van Wilgenburg SL, Clark RG, Wassenaar LI. A method for investigating population declines of migratory birds using stable isotopes: Origins of harvested lesser scaup in North America. PLoS One. 2009;4(11).

19. Hénaux V, Powell LA, Vrtiska MP, Hobson KA. Establishing winter origins of migrating lesser snow geese using stable isotopes. Avian Conserv Ecol. 2012;7(1).

20. Marra PP, Hobson KA, Holmes RT. Linking winter and summer events in a migratory bird by using stable-carbon isotopes. Science.1998;282(5395):1884–1886. pmid:9836637

21. Norris DR, Marra PP, Kyser TK, Sherry TW, Ratcliffe LM. Tropical winter habitat limits reproductive success on the temperate breeding grounds in a migratory bird. Proc R Soc B Biol Sci. 2004;271(1534):59–64.

22. Paxton KL, Moore FR. Carry-over effects of winter habitat quality on en route timing and condition of a migratory passerine during spring migration. J Avian Biol. 2015;46(3):495–506.

23. Hobson K A., Van Wilgenburg SL, Faaborg J, Toms JD, Rengifo C, Sosa AL, et al. Connecting breeding and wintering grounds of Neotropical migrant songbirds using stable hydrogen isotopes: a call for an isotopic atlas of migratory connectivity. J F Ornithol. 2014;85(3):237–257.

24. Alisauskas RT, Hobson KA. Determination of lesser snow goose diets and winter distribution using stable isotope analysis. J Wildl Manage. 1993;57(1):49–54.

25. Rubenstein DR, Hobson KA. From birds to butterflies: Animal movement patterns and stable isotopes. Trends Ecol Evol. 2004;19(5):256–263. pmid:16701265

26. Wassenaar LI, Hobson KA. Comparative equilibration and online technique for determination of non-exchangeable hydrogen of keratins for use in animal migration studies. Isotopes Environ Health Stud. 2003;39(3):211–217. pmid:14521282

27. Venables WN, Ripley BD. Random and mixed effects. In: Modern applied statistics with S. Fourth. New York, New York: Springer; 2002. p. 271–300.

28. Szymanski M, Afton AD, Hobson KA. Use of stable isotope methodology to determine natal origins of mallards at a fine scale within the Upper Midwest. J Wildl Manage. 2007;71(4):1317–1324.

29. Van Wilgenburg SL, Hobson KA, Brewster KR, Welker JM. Assessing dispersal in threatened migratory birds using stable hydrogen isotope (δD) analysis of feathers. Endanger Species Res. 2012;16(1):17–29.

30. Vander Zanden HB, Wunder MB, Hobson KA, Van Wilgenburg SL, Wassenaar LI, Welker JM, et al. Contrasting assignment of migratory organisms to geographic origins using long-term versus year-specific precipitation isotope maps. Methods Ecol Evol. 2014;5(9):891–900.

31. Hobson KA, Wilgenburg SL Van, Wassenaar LI, Larson K. Linking hydrogen (δ2H) isotopes in feathers and precipitation: Sources of variance and consequences for assignment to isoscapes. PLoS One. 2012;7(4):1–9.

32. Drake NA, Eckardt FD, White KH. Sources of sulphur in gypsiferous sediments and crusts and pathways of gypsum redistribution in southern Tunisia. Earth Surf Process Landforms. 2004;29(12):1459–1471.

33. Zazzo A, Monahan FJ, Moloney AP, Green S, Schmidt O. Sulphur isotopes in animal hair track distance to sea. Rapid Commun Mass Spectrom. 2011;25(17):2371–2378. pmid:21818798

34. Abraham KF, Jefferies RL, Alisauskas RT, Rockwell RF. Northern wetland ecosystems and their response to high densities of lesser snow geese and Ross’s geese. In Evaluation of special management measures for midcontinent lesser snow geese and Ross’s geese. Arctic Goose Joint Venture Special Publication. Washington, DC, Ottawa, Ontario; 2012.

35. Slattery SM, Alisauskas RT. Distribution and habitat use of Ross’s and lesser snow geese during late brood rearing. J Wildl Manage. 2007;71(7):2230–2237.

36. Alisauskas RT, Drake KL, Nichols JD. Filling a void: Abundance estimation of North American populations of Arctic geese using hunter recoveries. In: Thomson DL, Cooch EG, Conroy MJ, editors. Modeling Demographic Processes in Marked Populations. Boston, MA: Springer US; 2009. p. 463–489.

37. Kerbes RH, Meeres KM, Alisauskas RT, Caswell FD, Abraham KF, Ross RK. Surveys of nesting mid-continent lesser snow geese and Ross’s geese in eastern and central Arctic Canada, 1997–98. Technical Report Series No. 447. Saskatoon: Canadian Wildlife Service; 2006.

38. Morrissette M, Bêty J, Gauthier G, Reed A, Lefebvre J. Climate, trophic interactions, density dependence and carry‐over effects on the population productivity of a migratory Arctic herbivorous bird. Oikos. 2010;119(7):1181–91.

39. Hill RJ, Alisauskas T, Davison C. Influence of body size and condition on harvest and survival of juvenile Canada geese. 2003;67(3):530–41.

40. Abraham KF. Moult migration of lesser snow geese. Wildfowl. 1980;31(31):89–93.

41. Smith AD, Dufty AM Jr. Variation in the stable-hydrogen isotope composition of Northern Goshawk feathers: relevance to the study of migratory origins. Condor. 2005;107(3):547–58.

42. Langin KM, Reudink MW, Marra PP, Norris DR, Kyser TK, Ratcliffe LM. Hydrogen isotopic variation in migratory bird tissues of known origin: implications for geographic assignment. Oecologia. 2007;152(3):449–57. pmid:17370093

43. Haché S, Hobson KA, Villard M-A, Bayne EM. Assigning birds to geographic origin using feather hydrogen isotope ratios (δ2H): importance of year, age, and habitat. Can J Zool. 2012;90(6):722–8.

44. Alisauskas RT. Arctic climate, spring nutrition, and recruitment in midcontinent lesser snow geese. J Wildl Manage. 2002;66(1):181

45. Ross M V., Alisauskas RT, Douglas DC, Kellett DK. Decadal declines in avian herbivore reproduction: Density-dependent nutrition and phenological mismatch in the Arctic. Ecology. 2017;98(7):1869–83. pmid:28403519

46. Hobson KA, Lormée H, Van Wilgenburg SL, Wassenaar LI, Boutin JM. Stable isotopes (δD) delineate the origins and migratory connectivity of harvested animals: the case of European woodpigeons. J Appl Ecol. 2009;46(3):572–81