Combining Models of Environment, Behavior, and Physiology to Predict Tissue Hydrogen and Oxygen Isotope Variance Among Individual Terrestrial Animals

Variations in stable hydrogen and oxygen isotope ratios in terrestrial animal tissues are used to reconstruct origin and movement. An underlying assumption of these applications is that tissues grown at the same site share a similar isotopic signal, representative of the location of their origin. However, large variations in tissue isotopic compositions often exist even among conspecific individuals within local populations, which complicates origin and migration inferences. Field-data and correlation analyses have provided hints about the underlying mechanisms of within-site among-individual isotopic variance, but a theory explaining the causes and magnitude of such variance has not been established. Here we develop a mechanistic modeling framework that provides explicit predictions of the magnitude, patterns, and drivers of isotopic variation among individuals living in a common but environmentally heterogeneous habitat. The model toolbox includes isoscape models of environmental isotopic variability, an agent-based model of behavior and movement, and a physiology-biochemistry model of isotopic incorporation into tissues. We compare model predictions against observed variation in hatch-year individuals of the songbird Spotted Towhee (Pipilo maculatus) in Red Butte Canyon, Utah, and evaluate the ability of the model to reproduce this variation under different sets of assumptions. Only models that account for environmental isotopic variability predict a similar magnitude of isotopic variation as observed. Within the modeling framework, behavioral rules and properties govern how animals nesting in different locations acquire resources from different habitats, and birds nesting in or near riparian habitat preferentially access isotopically lighter resources than those associated with the meadow and slope habitats, which results in more negative body water and tissue isotope values. Riparian nesters also have faster body water turnover and acquire more water from drinking (vs. from food), which exerts a secondary influence on their isotope ratios. Thus, the model predicts that local among-individual isotopic variance is linked first to isotopic heterogeneity in the local habitat, and second to how animals sample this habitat during foraging. Model predictions provide insight into the fundamental mechanisms of small-scale isotopic variance and can be used to predict the utility of isotope-based methods for specific groups or environments in ecological and forensic research.

Variations in stable hydrogen and oxygen isotope ratios in terrestrial animal tissues are used to reconstruct origin and movement. An underlying assumption of these applications is that tissues grown at the same site share a similar isotopic signal, representative of the location of their origin. However, large variations in tissue isotopic compositions often exist even among conspecific individuals within local populations, which complicates origin and migration inferences. Field-data and correlation analyses have provided hints about the underlying mechanisms of within-site among-individual isotopic variance, but a theory explaining the causes and magnitude of such variance has not been established. Here we develop a mechanistic modeling framework that provides explicit predictions of the magnitude, patterns, and drivers of isotopic variation among individuals living in a common but environmentally heterogeneous habitat. The model toolbox includes isoscape models of environmental isotopic variability, an agent-based model of behavior and movement, and a physiology-biochemistry model of isotopic incorporation into tissues. We compare model predictions against observed variation in hatch-year individuals of the songbird Spotted Towhee (Pipilo maculatus) in Red Butte Canyon, Utah, and evaluate the ability of the model to reproduce this variation under different sets of assumptions. Only models that account for environmental isotopic variability predict a similar magnitude of isotopic variation as observed. Within the modeling framework, behavioral rules and properties govern how animals nesting in different locations acquire resources from different habitats, and birds nesting in or near riparian habitat preferentially access isotopically lighter resources than those associated with the meadow and slope habitats, which results in more negative body water and tissue isotope values. Riparian nesters also have faster body water turnover and acquire more water from drinking (vs. from food), which exerts a secondary influence on their isotope ratios. Thus, the model predicts that local among-individual isotopic variance is linked first to isotopic heterogeneity in the local habitat, and second to how animals sample this habitat during foraging. Model predictions provide insight into the fundamental mechanisms of small-scale isotopic variance and can be used to predict the utility of isotope-based methods for specific groups or environments in ecological and forensic research.

INTRODUCTION
Variations in stable hydrogen (δ 2 H) and oxygen (δ 18 O) isotopic ratios of animal body tissues are commonly used in terrestrial migration ecology to reconstruct origin and migration (reviewed by Hobson, 1999;Hobson and Wassenaar, 2019). The underlying premise of this method is that the isotopic compositions of body tissues are linked to those of local precipitation through predictable relationships (e.g., Chamberlain et al., 1996;Hobson and Wassenaar, 1997;Hobson et al., 2004; but see Hobson et al., 2012;Nordell et al., 2016;Magozzi et al., 2019). Because the isotopic composition of precipitation varies predictably across space and time (Craig, 1961;Bowen and Revenaugh, 2003; for a review see Bowen, 2010), variations in tissue isotopic compositions can be compared with spatial patterns in precipitation isotope ratios (isoscapes; e.g., Bowen et al., 2005b) to infer differences in tissue origin or migration. For example, bird feather δ 2 H values have been widely used to constrain the location of molting grounds (Chamberlain et al., 1996;Hobson and Wassenaar, 1997;Kelly et al., 2002;Hobson et al., 2004;Rubenstein and Hobson, 2004; reviewed by Hobson and Wassenaar, 2019). Oxygen isotope ratios are often measured simultaneously with δ 2 H values, and can also potentially provide useful information on molt location Hobson and Koehler, 2015).
A key premise underpinning the utilization of isotope-based methods to infer origin and migration is that tissue grown at the same site are isotopically similar. However, isotopic data from known-origin samples commonly show high levels of variation even among conspecific individuals living at a common site (e.g., Kelly et al., 2002;Meehan et al., 2003;Smith and Dufty, 2005;Wunder et al., 2005;Rocque et al., 2006;Langin et al., 2007;Gow et al., 2012;Haché et al., 2012;van Dijk et al., 2014;Reese et al., 2018). Such variability can often be quantified and incorporated in geographic assignment tests (Hobson et al., , 2014 but may limit the utility of such tests. Specifically, where the local isotopic variation is small relative to variation among locations, accurate assignments can be made; as the local variability increases, distributions for geographically distinct populations increasingly overlap and assignments become less certain (Wunder et al., 2005;Wunder, 2010;Hobson et al., 2012). Feather δ 2 H values of mountain plover (Charadrius montanus) chicks from 20 single breeding sites from Colorado to Montana, for example, spanned an average range of 41‰ across sites (average SD = 12‰; 5 ≤ n ≤ 164; minimum range = 16.3‰; maximum range = 109.2‰; Wunder et al., 2005;Wunder personal communication), severely limiting the use of isotopes to study migration in this species (Wunder et al., 2005). Values of feathers from adult American golden plovers (Pluvialis dominica) presumably grown over a small geographic area in Alaska varied over an even larger range (113‰; mean ± SD = −129 ± 69‰, n = 5; Rocque et al., 2006; but see Larson and Hobson, 2009). For adult American redstarts (Setophaga ruticilla) breeding in southeast Ontario, Canada, in contrast, the local range in feather δ 2 H values was only as high as 22‰ in 2004 (mean ± SD = −82 ± 4‰, n = 42; Langin et al., 2007) but even in this case assignment tests including local among-individual isotopic variation had reduced resolution (Langin et al., 2007). Over continental scales, within-site taxon-specific δ 2 H values spanned an average range of 25‰ across 13 sites and eight taxa (average SD = 8‰; 5 ≤ n ≤ 17; minimum range = 4‰ for Hylocichla mustelina; maximum range = 45.1‰ for Catharus ustulatus; Hobson et al., 2012).
Although we have a solid understanding of the processes controlling broad-scale isotope patterns (reviewed by Bowen, 2010), a more fundamental understanding of the underlying mechanisms of small-scale isotope variance is needed to improve estimates of such variance and assess the potential utility of isotope-based methods for specific groups or environments. Several studies have investigated the causes of within-site amongindividual isotopic variance through exploratory approaches such as correlation analysis and model selection (e.g., Clegg et al., 2003;Meehan et al., 2003;Smith and Dufty, 2005;Langin et al., 2007;Reese et al., 2018; see also Hobson et al., 2012;Nordell et al., 2016). Relationships between feather δ 2 H values and age have been observed for several taxa (e.g., Meehan et al., 2003;Smith and Dufty, 2005;Langin et al., 2007;Haché et al., 2012;Studds et al., 2012;Gomez et al., 2019) and attributed to age-related differences in diet, access to water, metabolism (Meehan et al., 2003;Langin et al., 2007;Storm-Suke et al., 2012), and/or evaporative water loss (Wolf and Martinez del Rio, 2000;Meehan et al., 2003, McKechnie et al., 2004Smith and Dufty, 2005; but see Kirkley and Gessaman, 1990;Marder et al., 2003). Differences in the magnitude and patterns of local among-individual variance observed for shore birds (Wunder et al., 2005;Rocque et al., 2006), raptors (Meehan et al., 2003;Smith and Dufty, 2005), and passerines (Langin et al., 2007) have been attributed to differences in molt phenology and evaporative water loss among these taxa (Smith and Dufty, 2005). Temporal variation in precipitation δ 2 H has been suggested as a driver of isotopic variation in accipiters (Meehan et al., 2003;Smith and Dufty, 2005). In general, these field-data point to two potential sources of variability: (1) differences in the isotopic composition of the diets of individuals, and (2) differences in individual physiology affecting the transfer of isotopes from diet to tissue. A critical comparison and evaluation of these two mechanisms in the context of empirical and theoretical constraints is lacking.
Discriminating among sources of variance potentially contributing to field-observed data is challenging, especially with incomplete sampling of the spatial and temporal variance in environmental isotope sources, and limited knowledge of the behavior and metabolic status of individuals. Simulation modeling offers an alternative approach, complementary to measurement of tagged individuals or manipulative fieldexperiments, to explore the potential contributions and interactions of isotope variance-inducing mechanisms to the variance in a sample population. The value of using such an approach is that it allows testing of multiple processes that could contribute to isotopic variation and their interactions, and offers the possibility of conducting experiments in the model environment where each model parameter can be controlled. Simulation models, however, involve assumptions and parameterization, and parameter values may be poorly constrained if data on the movement, behavior, and physiology of the study species are scarce or unavailable. In such cases, simulation using "idealized" models can still provide valuable information on expectations associated with different varianceproducing mechanisms or traits, and can be thought of as a framework for generating predictions associated with different hypothesized mechanisms rather than a tool to reproduce actual behavior or physiology. Simulation models of this type have previously been applied to explore how different migration hypotheses may contribute to variation in carbon isotopes in seabirds (Carpenter-Kling et al., 2019) and whales (Trueman and St John Glew, 2019; and oxygen isotopes in fish otoliths (Sakamoto et al., 2019).
Here we present a process-based modeling framework to estimate isotopic variation among local, conspecific individuals. The framework combines an agent-based model of animal behavior and movement, a dynamic version of the physiologybiogeochemistry model of Magozzi et al. (2019), and isoscape models of environmental isotopic variability. We use the resulting model toolbox to predict isotopic compositions of individual birds residing at a common but environmentally heterogeneous site, and explore how variation in behavioral, physiological, and environmental traits may interact to produce within-site among-individual isotopic variance. We compare model predictions against data from locally breeding, hatchyear Spotted Towhees (Pipilo maculatus) sampled prior to the start of fall migration in Red Butte Canyon (RBC), Utah, to evaluate which sets of model mechanisms appear most important in producing observed local among-individual isotopic variation.

Agent-Based Model
An agent-based (AB) model was developed to simulate the interaction of individual birds with the environment within RBC (Figure 1). The AB model is idealized and does not attempt to fully reproduce bird behavior (and drivers of behavior) but rather is intended to generate a set of ecologically and environmentally plausible hypotheses in order to experimentally manipulate potential sources of variance and explore the resultant effects on the simulated isotopic data. The environmental domain for the AB model was a representation of the spatial distribution of three major habitat types within the canyon: riparian, meadow, and slope (Figure 1). The domain was obtained through (1) classification of LiDAR data (National Ecological Observatory Network, 2017; http://data.neonscience. org) at 1 m by 1 m resolution into four vegetation cover type classes (Supplementary Figure A1), (2) aggregation of vegetation cover type class data to 10 m by 10 m resolution (Supplementary Figure A2), and (3) classification of these data into three habitat type classes (riparian, meadow, and slope; Appendix A). The riparian class comprised closed-canopy valley-bottom forest in the southwest portion of the canyon, including the area around the mist net sites (Figure 1). The meadow habitat represented open-canopy forest with an understory of grasses. Slope habitat included exposed soil, bunch grass, and shrub oak. Grid cell values of other environmental variables that can determine habitat quality (water and tree presence and prey density) were randomly drawn from distributions parameterized depending on habitat type. Water and tree presence were drawn from binomial distributions with probabilities of 50%, 1%, and 0.5% for water presence, and of 40%, 20%, and 5% for tree presence in riparian, meadow, and slope, respectively. Prey density values were obtained by drawing from a gamma distribution with shape and scale parameters giving mean prey densities of 250, 200, and 25 prey items per cell in riparian, meadow, and slope, respectively, and standard deviations of 1/10 of the mean in the former two habitats and of 1/2 of the mean in the latter habitat. Prey density was assumed to be higher in grid cells with trees than in cells without trees, and values drawn for cells with trees were multiplied by a random value drawn from a uniform distribution with min = 1 and max = 10.
The AB model simulated the movement and behavior of 105 individual birds over 60 days at hourly timesteps. Each bird was assigned a nest location to which it returned periodically over the simulation period. Two sets of simulations were run: in the first set, an equal number of birds (n = 35) were assigned random nesting locations within each habitat, allowing for among-habitat comparisons in the data analysis; in the second set, nesting locations were randomly selected without reference to habitat type, leading to a higher number of birds nesting in the more abundant slope habitat than in the meadow and riparian habitats. Because towhees nest predominantly on or near the ground in shrub-land habitat (https://www.audubon.org/field-guide/bird/ spotted-towhee), results from the second set of simulations were considered more representative of the local towhee population and were used in comparisons with the measured feather data.
At each model timestep, an individual makes probabilistic decisions about whether and where to move, whether to drink, or whether to eat. Probabilities of moving, drinking, and eating depend on: (1) time of day, (2) hunger and thirst, and (3) environmental conditions. Thirst and hunger are tracked using FIGURE 1 | Maps of location of Red Butte Canyon (RBC) within Utah (red star) and spatial distribution of three major habitat types in the canyon: riparian, meadow, and slope. Sampling sites for environmental water (black circles) and prey (blue circles), and locations of mist nets (red pluses) are also shown. two indices, which increase with time since the last drinking and feeding events. The rate of thirst increase during the day is twice that during nighttime, and hunger increase rate during the day is five times higher than at night ( Table 1). During the day, thirst increase rate is 7% higher than the "base" rate if in the slope habitat and 7% lower if in the riparian habitat (the base rate applies in meadow habitat). In addition, the rate of thirst increase is 7% higher than the habitat-specific rate if the bird moves, 7% lower if it does not move. Because the flux of vapor water out (F vw ; see below) is ∼1/3 of the total water loss for most modeled birds, a 7% shift in thirst growth rate depending on habitat and movement is proportional to a ∼20% shift in evaporative water loss (see section Physiology Model). Hunger increase rate (metabolic rate in the physiology model, see below) is 20% higher than the base rate if the animal moves and 20% lower if it does not move. Index values >0 indicate that the birds are thirsty and/or hungry.
When present, water is assumed to be available in excess, and therefore every time a bird drinks its thirst index resets to −1. In contrast, food availability is finite, and food intake is determined by food availability in the local environment and the individual's feeding success rate, which are multiplied to obtain a number of prey items eaten per hour; food intake is limited to a maximum of 20 prey items per hour. Eating success rate is randomly drawn from a normal distribution with mean equal to average eating success rate, which is a fixed property for each individual ( Table 1), and SD = 0.05; eating success rate has a minimum value of 0.01. When a bird eats, its hunger index decreases in proportion to the number of prey items eaten. As an individual feeds, prey density at the bird's location is concurrently reduced by the number of prey items eaten. Because this work does not focus on competition, the prey density map is re-initialized for each individual.
The probability that a bird will move from its current location depends on whether it is day or night, whether the bird is in its nest, whether it is thirsty and/or hungry, and whether there are adequate local water and/or food resources. At night, a bird that is not in its nest has a 99% probability to move to its nest; a bird that already is in its nest will not move (probability of moving = 0). During the day, a bird that is thirsty and/or hungry but lacks adequate local resources (i.e., water available if thirsty, [how much] food present if hungry) to meet its water and/or food demands has a 90% probability of moving; if resources are adequate, the probability of moving is low (10%). A bird that is neither thirsty nor hungry has a 50% probability to move to its nest (20% chance) or to any cell within the search area (i.e., random exploratory movement; 80% chance). When moving, a bird selects a destination cell within a square search area around its location (Table 1). Within the search area, optimal cells are defined as those that meet the individual's water and/or food requirements (i.e., if a bird is both thirsty and hungry, optimal cells are those with both water and adequate food), suboptimal cells are those that partially meet requirements (i.e., cells with either water or adequate food), and poor cells as those that do not  meet any requirement (i.e., cells with neither water nor adequate food). The probability that an individual selects any given cell is proportional to a random variable drawn from distributions shifted to lower values for less desirable cells ( Table 1).
The AB model was coupled offline with a physiologybiochemistry model that predicts the H and O isotopic compositions of animal body water and keratins from flux and isotopic inputs to the individual's body . Outputs from the AB model, including time of the day, whether movement and drinking occur, and the number of prey items eaten, were used to define H and O mass fluxes in the physiology model. The individual's location was tracked throughout the AB simulations and used to specify the H and O isotopic compositions of water and food consumed at each timestep based on environmental isoscape models for RBC.

Environmental Isoscape Models
Idealized high-resolution (10 m by 10 m) isoscapes were generated to represent spatial variation in H and O isotopic compositions of environmental water, liquid water within prey items, and prey biomass within RBC. Isoscape models simulated H and O isotope ratios of these resources in each grid cell as a function of habitat and were based on a combination of empirical data and theoretical considerations. Environmental water δ 2 H and δ 18 O grid cell values in the riparian habitat were drawn from a multivariate normal distribution the parameters of which were estimated from H and O isotopic data for 29 streamwater samples from Red Butte Creek downloaded from the Waterisotopes Database (Waterisotopes Database, 2018; http://waterisotopesDB.org) (Figure 1). Standing water is scarce outside the riparian corridor, but the AB model assumes that sparse water resources, for example small puddles or leafintercepted water, are available in a subset of grid cells in the meadow and slope habitats. Environmental water δ 2 H and δ 18 O values in these cells were simulated by drawing from the riparian distribution and adding 2 H-and 18 O-offsets representing evaporative enrichment. Values of 18 O-enrichment were drawn from a normal distribution with mean = 2‰ for meadow, and mean = 6‰ for slope, and SD = 1/2 of the mean for both habitats. 2 H-enrichments were estimated from 18 O-enrichments by assuming an evaporation line with a slope of 4-5 (uniform distribution with min = 4 and max = 5). Thus, both means and variances for environmental water increased from riparian to meadow and slope habitats (Supplementary Figure A7), reflecting the expectation of greater and more variable evaporation of sparse and transient water resources in the latter habitats.
Isoscape models for prey water and prey biomass δ 2 H and δ 18 O values were informed by data for insects and arthropods collected in 1 m by 1 m plots in each habitat at three times during the nesting season (Figure 1). Prey items were collected during the day, which is when the towhees forage, to obtain a relatively representative sample of potential dietary items. Diurnal or seasonal changes in prey activity, however, and the possibility that towhees might forage on some plant material (but see Smith and Greenlaw, 2015; https://birdsna.org/Species-Account/bna/species/spotow), could potentially result in larger isotopic ranges for prey items than those estimated here. Water was cryogenically extracted from individual prey items where possible, or from pooled samples including functionally and taxonomically related species where necessary, and analyzed for H and O isotope ratios using a Picarro L2130-i analyzer. The residual solid mass was homogenized using a mortar and pestle, lipid-extracted in a 2:1 chloroform:methanol solution, allowed to equilibrate with ambient water vapor alongside three keratin standards, vacuum desiccated, and analyzed for isotopic ratios using a Thermo Fisher TC/EA coupled with a Delta V isotope ratio mass spectrometer (for details, see analysis of feather samples). The 1σ precisions for δ 2 H and δ 18 O analyses for water samples were 0.8‰ and 0.25‰, respectively, and for solid mass samples were 2‰ and 0.3‰, respectively. Because sample numbers were relatively small (i.e., prey water: N = 9 in riparian, N = 4 in meadow, and N = 3 in slope), we chose to analyze these data in terms of isotopic offsets between substrates and habitats rather than using the raw data to characterize prey biomass and prey water isotope distributions directly (for plots of raw data, see Appendix A). For prey water, like for environmental water, evaporation is the primary process expected to produce isotopic differentiation among habitats. We estimated the 18 Oenrichment between prey water (δ 18 O pw ) and average stream water (δ 18 O streamw ) in each habitat, and thus calculated δ 18 O pw values as: Maximum possible 18 O-enrichment (ε evap,max ) was given by the difference between the most 18 O-enriched prey water sample (across habitats) and the average stream water δ 18 O value. Percent enrichment (ε evap,% ) -and thus prey water δ 2 H and δ 18 O values -in each grid cell were modeled using a beta distribution with α = 2 and β = 4 for riparian (most samples have small evaporative effects, few samples have large effects), α = β = 2 for meadow, and α = 4 and β = 2 for slope. Prey water δ 2 H values were estimated from δ 18 O values by using a regression line fit to the measured prey water H and O isotope data (Supplementary Figure A4) and adding dispersion by drawing from the residual distribution. The resulting distributions of simulated prey water δ 2 H and δ 18 O data reasonably approximated patterns in measured prey water data: prey water δ 2 H and δ 18 O values increased from riparian to meadow to slope, and among-sample isotopic variance was relatively similar across habitats (Supplementary Figures A3, A7). Data for field-collected prey samples were standardized for variation in H content, likely related to residual lipid content (details in Appendix A, section 2; Supplementary Figure A5 Figure A6), consistent with the expectation that differences in evaporative enrichment of prey body water should translate to differences in the isotope ratios of prey body tissues. Because the numbers of pooled samples and sampled groups varied substantially among habitats (i.e., ant: N = 1 in riparian, N = 1 in meadow, N = 1 in slope; Hemiptera: N = 2 in meadow, N = 2 in slope; fly: N = 1 in riparian, N = 2 in meadow, N = 1 in slope; salbug: N = 1 in riparian, N = 1 in meadow, N = 1 in slope; beetle: N = 3 in riparian, N = 2 in meadow, N = 1 in slope; spider: N = 1 in riparian, N = 2 in meadow), data from the meadow habitat (which was the most extensively sampled for prey) were used to characterize a multivariate normal reference distribution. Grid cell prey biomass δ 2 H and δ 18 O values for all habitats were drawn from this distribution, and values for cells in the riparian and slope habitats were then shifted by the mean within-functional group isotopic offsets for those habitats and the meadow. In this work, model isoscapes describing environmental isotopic variability are static with respect to time.

Physiology Model
We adapted the steady state physiology model of Magozzi et al. (2019) to develop a time-explicit model that predicts the H and O isotopic compositions of animal metabolic stores, body water, and keratins. To couple the AB and physiology models, we developed a "resource storage" sub-model representing the transient storage of metabolic resources between ingestion and metabolism. We modeled a single undifferentiated resource pool and did not attempt to track macronutrients separately with the exception that we tracked two O isotopic compositions: that of the carbohydrate + lipid component of the resource pool, and that of its protein component. The isotopic composition of the former is a function only of the isotope ratios of dietary inputs, whereas the latter is affected also by carbonyl-O exchange with water during digestion. The resource storage model is therefore characterized by four equations (Equations 2-5) that describe the change over time in: (1) O 2in and O 2out (mol O 2 h −1 ; Equations A1, A11) are fluxes of energy accumulated from eating and energy consumed via metabolism, respectively. The former is proportional to the number of prey eaten, and the latter equals the resting metabolic expenditure times hunger increase rate (therefore depends on time of the day and movement). Given this scaling, the size of O 2 is directly proportional to the individual's hunger index, coupling the AB and physiology models. The resting metabolic expenditure corresponds to the average energy expenditure during one nighttime hour and is estimated from the average energy expenditure during one daytime hour (based on the metabolic theory of ecology and rescaled to approximate measured metabolic rates in songbirds; Tieleman and Williams, 2000) by assuming that daytime metabolic rate is five times higher than nighttime metabolic rate (as for hunger growth rate in the AB model; Appendix A). RH p is the H isotopic composition of prey items. RO pcarb , RO pfat , and RO pprot are the O isotope ratios of dietary (prey) carbohydrate, fat, and protein, respectively (Equations A27-A29), and P carb , P fat , and P prot the relative proportions of each macronutrient in the diet (Supplementary Table A2). The O isotopic input for the protein component is calculated from isotope ratios of dietary protein (RO pprot ) and gut water (RO gw ; Equation A38) assuming partial exchange of carbonyl-O with water in the gut during amino acid hydrolysis (Schoeller et al., 1986), where 1 -P fO is the proportion of O exchanged (Supplementary Table A2). αOc-w is the fractionation factor associated with the interaction between carbonyl-O and water (Supplementary Table A2). Equations 4 and 5 implicitly assume that macronutrients are metabolized in proportion to their abundance in diet, a simplification adopted throughout our modeling. As a result, the proportion of macronutrient O in the resource pool equals that in diet, meaning RO res(carb+fat) = 18 O 2res(carb+fat) /(O 2 × (P carb + P fat )) and RO res(prot) = 18 O 2res(prot) /(O 2 × P prot ).
The change over time in (1) body water pool size (H 2 O, mol H 2 O), (2) 2 H content (mol 2 H 2 ), and (3) oxygen-18 content (mol 18 O) is described as:  Table A2) were estimated for a generic towhee (for model constants, see Supplementary Table A3). Differential equations for model state variables were solved using the R package "deSolve" (method used for solving ordinary differential equations: rk4; Soetaert et al., 2010). Variable values were returned at hourly timesteps. The initial resource pool size was set to 15% of body mass, assuming a mass to energy conversion appropriate for an insectivorous diet (Equation A36). Similarly, the initial body water pool size was set to 68% of body weight (Equation 37). Initial values for 2 H and oxygen-18 contents in the resource (carbohydrate + fat component and protein component in the case of O) and body water pools were obtained by running a steady state version of the model. Parameters for the "ode" algorithm are outputs from the AB model (i.e., time of the day, number of prey items eaten, whether drinking occurs, and habitat type) at hourly timesteps.
Following Ehleringer et al. (2008), the H isotopic composition of keratin synthesized at each timestep was calculated as: RH follw = P bw · RH bw + (1 − P bw ) · RH res (10) RH ker = P fH · RH res + (1 − P fH ) · αH prot · RH follw (11) RH follw is the H isotopic composition of follicle water, which is a mixture of body water (RH bw ) and metabolic H (RH res ; for P bw , see Supplementary Table A2). In turn, the H isotope composition of keratin (RH ker ) is a mixture of the resource pool and follicle water. Since no O is considered to be added at the site of keratin formation, the O isotopic composition of newly synthesized keratin (RO ker ) equals that of the protein component of the resource pool (RO res(prot) ) at the time of synthesis.

Model Experiments
Nine model experiments were conducted by running the model system with different parameter sets (Table 1). Behavior, hereafter "B", experiments included among-individual variation in traits influencing bird behavior and nutritional physiology (thirst and hunger increase rates during the day and at night, average eating success, probability to select optimal and suboptimal cells over poor cells, and search distance) but no environmental isotopic variability. Parameter values were sampled from uniform distributions spanning the mean ± 10%. Although no empirical constraints were available to inform the selected range of values for the behavioral parameters the ±10% range was adequate to assess the nature of any isotope effects associated with this model component; for the physiology parameters, the range adopted here is consistent with reported inter-individual ranges in metabolic rate for many songbirds (Ambrose et al., 1996;Smit and McKechnie, 2015;Cooper et al., 2019). In these experiments the isotopic compositions of water and food resources were the same in all grid cells. Environment, hereafter "E, " experiments included variability in environmental water, prey water, and prey biomass isotopic compositions but no variation in behavioral traits among individuals. At each drinking and feeding event, isotopic inputs were drawn from habitatand substrate-specific isoscape distributions described above.
In the case of feeding, prey water and prey biomass isotopic compositions were means for the number of prey items eaten. For each set of experiments, we ran a reference case in which behavioral parameters were centered on reference values for towhees (B1 and E1), and three tests in which the species-average values of different groups of behavioral parameters were varied to test the sensitivity of the results to differences in the "typical" behavior of individuals. These tests included increased (by 50%) demand for food and water (B2A and E2A), improved (by 100%) foraging success and habitat selection (B2B and E2B), and increased (by 100%) mobility (larger search distance; B2C and E2C). A final experiment, E1B1, included both among-individual behavioral variation and environmental isotopic variability, using the reference values for behavioral parameter means ( Table 1).

Feather Collection and Isotopic Analysis
Spotted Towhee was selected as a target species as hatch-year individuals are abundant in RBC over the breeding season.
Locally grown feathers from hatch-year towhees were collected in the canyon prior to the start of fall migration in September (Smith and Greenlaw, 2015; https://birdsna.org/Species-Account/bna/ species/spotow) in 2013-2016. After hatching in May through July, towhees undergo a first pre-basic molt in April-July, which is complete and occurs in the nest. Between July and October, they undergo a preformative molt on or near breeding grounds, which includes body feathers but (in ca. 60-65% of the individuals) no rectrices or primaries. By sampling the feathers prior to start of fall migration, we eliminated the possibility of collecting feathers from birds that could have hatched somewhere else and used RBC as a stop-over site during migration. Although no home range data exist for hatch-year towhees, it is likely that these birds are tied to their nest location until migration. Birds were captured with mist nets located at the bottom of the canyon, near the river, in riparian, meadow, and slope habitats (for net sites, see Figure 1). For each capture, various metadata and biological data were recorded, and one or two tail feathers were pulled. Tail feathers are a common choice for analysis in isotopic research, and no evidence exists for significant isotopic differences between tail and flight feathers (e.g., Hobson et al., 2012;Nordell et al., 2016).  = 11.19‰;Chesson, 2012) were allowed to equilibrate with water vapor in the laboratory atmosphere for 7 days, then were vacuum dehydrated prior to isotope analysis (Bowen et al., 2005a). Values used for the reference materials were anchored to the KHS and CBS standards (KHS: δ 2 H = −54.1‰, δ 18 O = 21.21‰; CBS: δ 2 H = −197.0‰, δ 18 O = 2.39‰; Qi et al., 2011;Soto et al., 2017) based on an inter-laboratory roundrobin test. Samples and reference materials were transferred to a Zero Blank Autosampler (Costech Analytical) interfaced with a Thermo-chemical Elemental Analyzer (Thermo Fisher Scientific). They were pyrolized at 1400 • C in an oxygenfree environment to form H 2 and CO gases, which were chromatographically separated and introduced sequentially into the source of a Delta Plus XL isotope ratio monitoring mass spectrometer (IRMS; Thermo Fisher Scientific). The data for the reference materials were used to calibrate raw measurements of the δ 18 O and non-exchangeable δ 2 H values for unknown samples relative to the VSMOW-SLAP standard scale using the two-point regression-based approach of Wassenaar and Hobson (2003). Analytical precision, based on the repeated analysis of a check reference material (POW), was better than 2‰ and 0.3‰ (1σ) for δ 2 H and δ 18 O values, respectively.

Data Analysis
In order to remove potential influences of initial conditions on state variable values, only outputs over last 30 days of simulation were considered for data analysis. Individual average keratin δ 2 H and δ 18 O values were calculated over this time, the length of which also likely approximates the isotope incorporation time for bird feathers (Hobson, 1999 and refs. there in). Linear regression models were used to examine relationships between keratin isotopic compositions and model variables including body water turnover time, proportion of water influx coming from drinking (vs. food), and mean input flux isotopic ratios. Body water turnover time (h) was calculated as the ratio of individual average body water pool size (H 2 O, mol H 2 O) to output flux (F vw + F wasteH + F lw ; mol H 2 O h −1 ). Mean input flux isotopic composition equals the mean of individual average isotope values for all inputs (e.g., δ f , δ fw , δ dw , and δ 18 O O2 for O), weighted by their relative contributions to the total water flux.
Relationships between keratin isotopic compositions and multiple predictor variables were described by multiple linear regression models of the form Y = β 0 + β 1 X 1 + β 2 X 2 + β 3 X 3 + . . . + β i X i + ε, and visualized with added variable (partial regression) plots, which show the bivariate relationships between the independent variable Y and each of the predictor variables (X i ), conditioned on the other predictors. To compare modeled vs. measured isotopic variances, mean squared deviances of keratin δ 2 H and δ 18 O values were compared among model experiments and data; means, sample variances, SDs, and ranges for individual average keratin δ 2 H, and δ 18 O values were also calculated for each population.
Complete code to run the models and experiments and generate simulated data is available online (Bowen and Magozzi, 2020).

AB Model Results
Nesting location and habitat strongly influenced the habitat(s) in which individuals spend most time (Figure 2). Across experiments and nesting habitats, individuals preferentially occupy the more resource-rich riparian habitat (relative to its spatial coverage). This effect is accentuated in the experiments with higher food and water demands (e.g., B2A), presumably because individuals are more strongly dependent on the high-resource grid cells to meet their requirements. In contrast, slope-nesting birds spend slightly less time in the riparian habitat in the experiments in which foraging success and habitat selectivity and search distance were increased (e.g., B2B and B2C, respectively), reflecting their greater ability to meet their needs without conducting long-range foraging bouts to the riparian zone. Regardless, the dominant pattern in model results is that individuals are strongly tied to their local (nesting) habitat, suggesting that, within the model framework, nest location drives the locations and conditions of foraging for individuals within a local population.
Different habitat utilization could affect the tissue isotopic compositions of individuals in several ways. In the B experiments, among-habitat differences in resource density translate into differing drinking and feeding frequencies and amounts among individuals associated with different habitats. In general, birds associated with the riparian habitat have greater access to water, and therefore drink more frequently than meadowand slope-nesting birds (Figure 3A), potentially affecting their body water and isotope balance. Riparian nesters also have greater access to adequate food resources to meet their food demands, and therefore eat less frequently but consume greater amounts of food at each feeding event (Figure 3B). In the less resource-rich habitats individuals spend a higher proportion of time actively foraging, particularly in the B2A experiment and slope-nesting individuals in the B2B experiment. Some individuals are unable to obtain sufficient nutrition and completely deplete their resource pool, which we treated here as lethal. Deaths are most common among slope-nesting birds in the B2A and E2A experiments (18 and 16 dead individuals, respectively), due to the increased food demand in these runs. Among-individual variation in drinking and feeding frequencies is, in general, slightly higher for slope-nesting birds than others across the B experiments (Figure 3).

Behavioral Variability
Even in the absence of environmental isotopic variation (the B experiments), keratin δ 2 H and δ 18 O values differ among modeled individuals and habitats, with the lowest values in individuals nesting in the riparian habitat and highest values in birds nesting in the slope habitat (Figure 4). These patterns are consistent across the experiments, with subtle differences. The lowest isotope values in all habitats occur when average food and water demand is high (the B2A experiment). Across all B experiments, among-individual isotopic variation is somewhat higher in slope-nesting birds than in other habitats. In general, among-individual isotopic variance in the B experiments is small (i.e., a few ‰ for H; Figure 4). Analysis of the model output suggests that 94-95% of H and O isotopic variation among B experiments, habitats, and individuals can be attributed to two variables: (1) body water turnover time, and (2) the proportion of water influx coming from drinking (p Fdw ; Figure 5, Supplementary Figure B1). Although neither variable was forced directly, both respond to differences in resource availability within the model framework, wherein individuals with greater access to water and higher drinking frequency ingest more drinking water, driving greater excretion of liquid waste and shorter turnover times.

Environmental Variability
Between-habitat isotopic variation in environmental substrates in the E experiments produces a much wider range of modeled keratin δ 2 H and δ 18 O values than the behavioral/physiological FIGURE 2 | Mean proportion of time spent in each habitat for each experiment (see Table 1) and nesting habitat. Proportions of the study area covered by the riparian habitat (blue dashed line) and by the riparian and meadow habitats (red dashed line) are also shown.  Figures 4, 6). The relative pattern among birds nesting in different habitats is the same as that observed in the B experiments, but within-and between-habitat variability in keratin isotopic compositions is substantially higher for all the E experiments. The range of tissue isotope values is somewhat reduced in the E2A and E2C experiments, in which greater demand for resources and/or a larger search distance induce more vigorous and spatially dispersed foraging. In contrast, higher feeding success rates in the E2B experiment allow slope-nesters to focus their foraging more within the slope, and slightly increase the difference in mean values between this habitat and the others. Among-individual variability within nesting habitats is also much greater than in the B experiments, with the largest ranges for birds nesting in the slope (Figure 6).
Among-individual variation in mean input flux isotopic ratios is the dominant driver of variability in keratin isotopic compositions in the E experiments, with secondary influence of variation in turnover time and the proportion of water influx  and body water turnover time (X 1 ), given proportion of water influx coming from drinking (p Fdw ; X 2 ) for each B experiment. Relationships do not differ among experiments, but keratin absolute isotope values in B2A birds are more negative than in other B birds due to lower turnover rates, despite higher p Fdw (Supplementary Figure B1). Multiple linear regression models are: for δ 2 H, Y = −97.41 + 0.12X 1 -20.35X 2 , R 2 = 0.95; for δ 18 O, Y = 6.78 + 0.04X 1 -5.63X 2 , R 2 = 0.94. coming from drinking (Figure 7; Supplementary Figures B2,  B3). Together, these variables explain 97% of H and O isotopic variation predicted across the E experiments. These results suggest that isotopic variability among habitats within a local environment, combined with selective utilization of habitats as discussed above, may produce substantial interindividual isotopic variation in consumer tissues. The extent of this effect is sensitive to the behavioral characteristics of species and individuals, and is accentuated at the population level when individuals forage locally within isotopically distinct habitats. The extent of environmental heterogeneity, too, is predicted to be a strong determinant of inter-individual isotopic variability, within the extreme case of the B experiments (isotopically uniform environment) suggesting an order of magnitude-reduction in variability among individuals.  and mean input flux isotope ratios (X 1 ), given body water turnover time (X 2 ) and proportion of water influx coming from drinking (X 3 ) for each E experiment. Relationships do not differ among experiments, and keratin absolute isotope values span similar ranges across behavioral properties due to factors producing offsets in isotopic relationships among experiments held constant here (but see Supplementary Figures B2, B3). Multiple linear regression models are: for δ 2 H, Y = 2.58 + 1.28X 1 + 0.10X 2 + 26.73X 3 , R 2 = 0.97; for δ 18 O, Y = −5.17 + 1.05X 1 + 0.03X 2 + 19.26X 3 , R 2 = 0.97.

Comparisons and Measured Data
Lastly, we compare simulated isotopic variation with variance observed in locally grown feathers from hatch-year Spotted Towhees in RBC. Measured feather δ 2 H values (mean ± SD) were −105.32 ± 12.97‰ (range = 59.32‰) and δ 18 O values were 13.1 ± 2.53‰ (range = 11.85‰; Supplementary Table B2). Results for the B experiments show much smaller variation in keratin δ 2 H and δ 18 O values than do either the E and E1B1 experiments or field-collected towhee data (Figure 8; Supplementary Table B2). Variation is more similar in the E and E1B1 experiments, implying that the among-individual physiological differences considered here accentuate the tissue isotopic variability imparted by environmental heterogeneity, but by a relatively small amount. These experiments approximate the variation in the observed samples well, meaning that the model toolbox reproduces observed within-site amongindividual isotopic variance only when environmental isotopic variability is accounted for. Among the E experiments, betweenindividual variation is smallest in the E2A and E2C experiments (Figure 8). In the former case, that is likely because 16 slopeassociated model individuals died, compressing variation toward ranges typical of the riparian and meadow habitats; in the latter, reduced variance is attributable to the larger foraging area of individuals.

DISCUSSION
Within-site among-individual variation in consumer tissue H and O isotope ratios has been documented across a range of biological systems. Despite its critical importance for the use of these isotope systems in ecology, the drivers and characteristics of this variation are generally poorly understood (but see, e.g., Hobson et al., 2012). On one hand, differences in behavior and physiology might produce differences in the isotope massbalance and biochemistry of individuals through mechanisms such as different metabolic rates, water loss pathways, or diet composition (Wolf and Martinez del Rio, 2000;Meehan et al., 2003;Smith and Dufty, 2005;Langin et al., 2007). On the other, individuals that forage within isotopically heterogeneous environments might access resources with different isotope ratios, contributing isotopic variation to their tissues (Meehan et al., 2003;Smith and Dufty, 2005;Hobson et al., 2012).
Here we present the first model-based assessment of some of these mechanisms, evaluating whether, within the context of our current understanding, these different processes can replicate observed variability within a natural system. We benchmark our model results against a new dataset of keratin H and O isotopic compositions for conspecific individual birds growing feathers at a common site. Datasets of this kind are valuable to quantify the magnitude of local among-individual isotopic variance and investigate the causes of such variance, but only a few datasets currently exist (e.g., for H: Meehan et al., 2003;Wunder et al., 2005;Langin et al., 2007;Hobson et al., 2012Hobson et al., , 2014for O: Hobson and Koehler, 2015). Our dataset provides further evidence of large within-site among-individual isotopic variance. Values of δ 2 H and δ 18 O for towhees span ∼60 and >10‰ ranges, respectively. The δ 2 H range is far larger than that observed in American Redstarts (12-20‰; Langin et al., 2007), and more similar to those measured for shorebirds (60 and 113‰; Wunder et al., 2005;Rocque et al., 2006; but see Larson and Hobson, 2009) and raptors (∼50‰; Meehan et al., 2003). Although ranges observed for towhees do not exceed variation in precipitation between breeding and non-breeding grounds for this species (Bowen et al., 2005b), the substantial local among-individual isotopic variance documented here would dramatically reduce the certainty and usefulness of assignment tests for this species.
The importance of heterogeneity in the isotopic compositions of environmental substrates as a control of within-population isotopic variance of animals has been proposed previously (Meehan et al., 2003;Smith and Dufty, 2005), but few studies have measured this environmental variability. Here we document substrate heterogeneity by sampling potential prey items. Even though our results are based on small numbers of pooled samples (see section Environmental Isoscape Models), ranges in prey water and organic matter isotope ratios within RBC are large, and among-habitat differences provide evidence for evaporation as the dominant underlying driver of observed patterns (Supplementary Figures A3, A6). Even though we did not measure spatial variability in drinking water directly, similar patterns are expected, assuming that any drinking water available in the non-riparian habitats comes from transient puddles and/or leaf-intercepted water subject to relatively large evaporative enrichments (Supplementary Figure A7). Given the environmental isotopic variability that we observed (or modeled, for drinking water), preferential use of resources from differing, isotopically distinct habitats within the study area could offer a plausible mechanism for the high variance observed in towhee feathers.
We test this idea using a mechanistic modeling framework and a range of different assumptions about environmental, behavioral, and physiological heterogeneity. The models do reproduce the magnitude of the observed isotopic variance in the towhee population, but are only able to do so when environmental isotopic variability is accounted for (E and E1B1 experiments; Figure 8; Supplementary Table B2). This suggests that variation in diet isotopic ratios may be a strong driver of variance within this population (Figure 7). Small differences exist between measured and predicted variance even in the E and E1B1 models, which could reflect the inappropriateness of some model assumptions. For example, we use data from hatch-year birds to ensure that our dataset only includes locally grown feathers.
Hatch-year birds are often more isotopically variable than adults, which could be attributed to their limited access to drinking water (Meehan et al., 2003;Smith and Dufty, 2005;Langin et al., 2007), variable nest microclimate, habitat segregation between molting adults and hatch-year birds, and/or difference in feather growth rate. These factors may not be adequately reflected in the idealized model used here.
We find that behavioral and physiological differences among individuals, in this case largely driven by differences in local habitat quality, contribute more subtly to amongindividual isotopic difference than does environmental isotopic heterogeneity. Habitat characteristics affect individual drinking and foraging behavior (Figure 3), body water balance, and tissue isotopic compositions (Figures 4, 5). In the model system this effect accentuates that of environmental isotopic variability in an isotopically heterogeneous environment. This prediction is consistent with the proposal that difference in body water balance might contribute to the inter-individual isotopic variability observed in local populations (Wolf and Martinez del Rio, 2000;Meehan et al., 2003;Smith and Dufty, 2005) but also with a previous model-based assessment suggesting that this contribution is likely to be relatively small .
These results are somewhat dependent on the range of parameter values tested. Although the ranges chosen here are generally consistent with observed data, where these are available (Goldstein, 1988), the values and variability of some parameters are poorly constrained. Comparison of the inter-individual variation expressed within our individual B experiments, where the prescribed range of trait values was relatively small, with the variation across experiments (where mean trait value differences were 5-10 times larger; Table 1), however, shows low values of behaviorally and physiologically driven isotopic variance relative to the E experiments in all cases.
Moreover, all of the physiological mechanisms that we have identified as potential drivers of significant isotopic variation, including inter-individual differences in foraging success and drinking and body water turnover, are in turn strongly forced by environmental heterogeneity. Thus, we suggest that environmental heterogeneity, whether via variation in substrate isotope ratios or behavioral and physiological correlates, may be the dominant, first-order control on inter-individual H and O isotope variance within consumer populations. For example, the high variance observed in our system relative to that observed by Langin et al. (2007) might in part reflect the influence of aridity and strongly contrasting evaporative effects among habitats in the RBC environment. Similarly, in coastal environments generalist species that use both terrestrial and marine resources might be more isotopically variable (Rocque et al., 2006; but see Larson and Hobson, 2009), particularly where individuals specialize their feeding within this context (Vander Zanden et al., 2010). High local among-individual isotopic variance might therefore be expected in habitats characterized by large environmental and isotopic heterogeneity, particularly those encompassing gradients (e.g., humidity, terrestrial to marine) along which the isotopic composition of environmental substrates and/or stressors affecting the body water balance of individuals vary extensively. In realworld settings, temporal environmental variability may act to enhance or reduce the effect of variation among habitats, depending on its scale and frequency. This first-order expectation represents a hypothesis to be tested with additional fieldcollected samples, and complexities associated with patterns of habitat use and specialization may strongly modulate or overwhelm differences in environmental variability for some systems. If our hypothesis is supported, however, it would suggest that uncertainty in isotope-based methods will generally be higher for species occupying environments characterized by larger isotopic and/or environmental heterogeneity, implying predictable differences in the utility of these methods for different groups of organisms, depending on the characteristics of their environment.
Our model, like any other, is clearly not comprehensive and other relevant mechanisms not considered here may also contribute to consumer H and O isotope variability. Potential factors may include changes in metabolic pathways (e.g., increased reliance on lipid metabolism; Estep and Hoering, 1980;Sessions et al., 1999) during exercise or starvation, and differences in dietary macronutrient composition (Meehan et al., 2003;Langin et al., 2007). Given this, we emphasize that our results represent predictions of isotopic patterns expected given a limited suite of environmental, physiological, and behavioral mechanisms and range of parameter space. The results appear to be reasonable and meaningful when tested against the towhee dataset, but future models should incorporate and explore additional mechanisms, evaluating their effects on isotopic variance and testing whether they are consistent across systems or system-specific.

CONCLUSIONS
In this work we introduce a modeling framework providing explicit predictions of the magnitude and patterns of isotopic variance in H and O among individual terrestrial endotherms (birds) living at a common site under different sets of assumptions. This set of tools offers new opportunities to explore and constrain the potential mechanisms giving rise to such variance. Particularly, we demonstrate that simulation modelbased inference can successfully be used to understand the potential effects of multiple variance-generating processes on a measured variable (here isotopic composition of body tissues), offering a complementary approach to data-based inference. Within the model framework, isotopic heterogeneity in the local habitat emerges as the dominant driver of local among-individual isotopic variance, with its effect accentuated by differences in body water balance. If verified, this implies that amongindividual variance in body tissue H and O isotopic compositions may be strongly linked to the degree of isotopic heterogeneity in local environmental H and O sources.
Future studies should investigate the movement and behavior of birds in natural environments and data should be used to adapt the model mechanisms to specific groups or taxa. For example, insectivorous birds might acquire most water from diet, which would dampen the effect from variable contribution of drinking water to the body water pool seen here. Similarly, the post-fledging movements of hatch-year birds are poorly known but the habitat used by adult birds during molt to a prebasic plumage before migration is believed to greatly influence the variance measured in feathers from unknown locations and such a shortcoming in our understanding would also deserve further investigation. Future work should also test additional mechanisms that might potentially modify the magnitude and patterns of such variation such as nutritional status and macronutrient composition, and assess the consistency of each effect across systems. The models developed here may be used to predict the usefulness of isotope-based applications for specific groups of organisms or environments and develop optimal, system-specific sampling strategies in ecology and forensic research. As greater confidence in model mechanisms and predictions develops, these models could also be ultimately applied to produce variance estimates to be used in statistical estimation of origin where measurements are inadequate to do so.

DATA AVAILABILITY STATEMENT
Code to run the models and generate datasets used in this study can be found online (Bowen and Magozzi, 2020).

ETHICS STATEMENT
The animal study was reviewed and approved by University of Utah Institutional Animal Care and Use Committee (number 18-02002). Written informed consent was obtained from the owners for the participation of their animals in this study.

AUTHOR CONTRIBUTIONS
SM and GB conceived ideas and designed methodology, with support from CT, HVZ, and MW. DP and PD analyzed the LiDAR data. JH and CS provided the feather samples. KP collected samples of environmental substrates and prepared substrate and feather samples for isotopic analyses. SM and GB led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

FUNDING
This work was supported by U.S. National Science Foundation Grants EF-1241286 and DBI-1565128.