Isotopic niche partitioning and individual specialization in an Arctic raptor guild

Intra- and inter-specific resource partitioning within predator communities is a fundamental component of trophic ecology, and one proposed mechanism for how populations partition resources is through individual niche variation. The Niche Variation Hypothesis (NVH) predicts that inter-individual trait variation leads to functional trade-offs in foraging efficiency, resulting in populations composed of individual dietary specialists. The degree to which niche specialization persists within a population is plastic and responsive to fluctuating resource availability. We quantified niche overlap and tested the NVH within an Arctic raptor guild, focusing on three species that employ different foraging strategies: golden eagles (generalists); gyrfalcons (facultative specialists); and rough-legged hawks (specialists). Tundra ecosystems exhibit cyclic populations of arvicoline rodents (lemmings and voles), providing a unique system in which to examine predator diet in response to interannual fluctuations in resource availability. Using blood δ13C and δ15N values from 189 raptor nestlings on Alaska’s Seward Peninsula (2014–2019), we calculated isotopic niche width and used Bayesian stable isotope mixing models (BSIMMs) to characterize individual specialization and test the NVH. Nest-level specialization estimated from stable isotopes was strongly correlated with indices of specialization based on camera trap data. We observed a high degree of isotopic niche overlap between the three species and gyrfalcons displayed a positive relationship between individual specialization and population niche width on an interannual basis consistent with the NVH. Our findings suggest plasticity in niche specialization may reduce intra- and inter-specific resource competition under dynamic ecological conditions.


Introduction
Trophic ecology is largely centered on the concept of a realized trophic niche-described as all resources an organism actually consumes-which is inherently communitydependent and limited both by local resource availability and inter-and intra-specific competition (Hutchinson 1957). Analyzing trophic niche width can reveal spatial and temporal patterns of resource use, determine the degree of dietary generalism or specialization for a population, and serve as a metric of resource competition (Schoener 1982). Characterizing changes in an organism's trophic niche in response to dynamic ecological conditions can illuminate its role in greater ecosystem function and indicate its resilience to environmental change.
Two predominant paradigms exist to explain how a population may alter its trophic niche in response to variable environmental conditions: Optimal Foraging Theory (OFT; Macarthur and Pianka 1966), and the Niche Variation Hypothesis (NVH; Van Valen 1965). Classic OFT postulates that stabilizing selection leads to the evolution of behavioral strategies that maximize the net rate of energy intake-resulting in preferred diets that are consistent across all individuals within a population (Pyke et al. 1977). Conversely, the NVH posits that populations are composed of Communicated by Seth Newsome. individuals that specialize on different subsets of available diet items, and thus a generalized population's niche width reflects the combined niches of many individual specialists (Bolnick et al. 2002). Although the NVH was initially postulated in the context of evolved individual specialization in populations via destabilizing selection (Van Valen 1965), niche variation can also arise as a plastic response to variability in resource abundance, potentially acting as a mechanism for avoiding resource competition (Svanbäck and Bolnick 2005). Under the concept of ecological opportunity, populations aligning with the NVH paradigm are predicted to expand their niche via increased individual dietary specialization as a greater diversity of prey becomes available (Araujo et al. 2011;Bolnick et al. 2007).
Individual specialization is an important driver of food web dynamics in a variety of taxa (reviewed in: Bolnick et al. 2003). The degree to which populations exhibit individual specialization under the NVH paradigm is largely dependent upon that population's degree of intra-specific trait variation associated with different phenotypes, experience, and personalities (Bolnick et al. 2011;Toscano et al. 2016). Species that are long-lived, exhibit some degree of site and pair fidelity, occupy upper trophic levels, and express the ability to learn complex behavior are expected to have high degrees of individual specialization (Araújo et al. 2008). On an inter-specific basis, generalist species are more likely to fall under the NVH paradigm than obligate dietary specialists which may be more constrained in these phenotypic or behavioral traits (Bison et al. 2015;Maldonado et al. 2017). By maintaining heterogeneity in foraging strategies, populations composed of dietary specialists mitigate intraspecific competition which promotes population stability under resource fluctuations (Bolnick et al. 2003). As such, evaluating whether a population exhibits individual specialization in alignment with the NVH has meaningful implications for its ecological resilience, resource competition, and resultant food web dynamics (Bolnick et al. 2011).
The Arctic raptor guild comprises the upper trophic level in tundra ecosystems, wherein several species coexist through slight differences in their realized niche spacetheoretically by partitioning limiting resources via small differences in phenology, life-history strategy, morphology, and foraging behavior. Arctic tundra ecosystems are characterized by pulses in resource availability, often connected to the abundance of cyclic populations of arvicoline rodents (lemmings and voles; Myodes, Microtus, and Lemmus spp.). Highs and lows within these population cycles have effects that reverberate through the food web, directly affecting the predator guild that relies on them and indirectly influencing predation pressure on alternate prey species (Ims et al. 2008). Although raptors exhibit characteristics that we may expect to be consistent with the NVH (e.g., long lifespan, site fidelity, capacity to learn complex behaviors; Araújo et al. 2008), many raptor dietary studies have largely ignored the implications of individual specialization on their trophic ecology (but see L'Herault et al. 2013). Further, Arctic raptor dietary studies typically focus on a single species, thus the trophic dynamics of sympatric species remain largely understudied, particularly under varying degrees of resource availability.
To investigate the trophic dynamics of the Arctic raptor guild, we focus on three sympatric species with different foraging strategies: the golden eagle (Aquila chrysaetos; a dietary generalist; Katzner et al. 2020); the gyrfalcon (Falco rusticolus; a facultative specialist; Booms et al. 2020); and the rough-legged hawk (Buteo lagopus; a dietary specialist on arvicoline rodents; Bechard et al. 2020). Although these species nest sympatrically and their diets overlap to some extent, it remains unclear how inter-specific competition and niche partitioning influence their trophic ecology. By focusing on the trophic ecology of Arctic raptors under dynamic ecological conditions (i.e., cyclic arvicoline rodent abundance), we seek to assess their agreement with the NVH paradigm, thereby characterizing their behavioral plasticity and providing broader context for the tundra food web.
We quantified the niche space of the three focal raptor species and compared their response across years with fluctuating resource availability. We further produced proportional diet estimates for each species and used them to assess the degree of individual (nest-level) dietary specialization within each population. We anticipated the focal raptor populations would exhibit varying degrees of individual diet specialization, and based on the NVH, we predicted that population niche width expands via increases in individual specialization on an interannual basis. We further anticipated that these shifts in population niche width would alter dietary overlap among species, potentially contributing to interannual variability in competition within the raptor guild.

General approach
Methods for estimating trophic niche space have advanced over recent years with the integration of stable isotope approaches, resulting in the concept of an "isotopic niche" (Bearhop et al. 2004;Newsome et al. 2012;Shipley and Matich 2020). Animal tissue is largely composed of carbon and nitrogen that naturally occur in both heavy and light isotopes, and their ratios (e.g., 13 C/ 12 C [δ 13 C] and 15 N/ 14 N [δ 15 N]) provide a useful tool for investigating trophic ecology and estimating diets integrated over time. Variation in δ 13 C and δ 15 N exists among taxa, which ultimately stems from the basal carbon and nitrogen of their food web and is integrated by trophic level in a predictable manner (e.g., Phillips et al. 2014). Thus, individuals and populations will occupy different isotopic space based on the resources they consume, making an organism's isotopic niche a useful proxy for its trophic niche (e.g., Herrera and Rodríguez 2013), although the isotopic niche and trophic niche may not be directly equivalent in some circumstances (e.g., Hette-Tronquart 2019). Isotopic niche width can be directly estimated and compared across communities via the calculation of stable isotope Bayesian ellipses (Jackson et al. 2011) or via kernel density estimates, which account for holes in utilized isotopic space (Eckrich et al. 2020). Further, by iteratively analyzing a consumer's isotopic position relative to that of its resources, Bayesian stable isotope mixing models (BSIMMs) can generate proportional dietary estimates (e.g., MixSIAR, Stock and Semmens 2016). Together, these methods create a framework to quantitatively analyze isotopic niche dynamics and estimate which resources may be driving them. Resulting proportional diet estimates can be used to distinguish within-individual from between-individual components (sensu Roughgarden 1972) of the diet to evaluate broad trophic hypotheses associated with individual specialization (i.e., the NVH).

Study system
We studied three sympatric raptor populations on the Seward Peninsula on the west coast of Alaska (64° N-65° N, 164° W-166° W; Fig. S1). The study area (14,150 km 2 ) is a spatially heterogeneous Arctic tundra ecosystem characterized by low-lying vegetation (e.g., Ericaceae and Saxifragaceae spp.) at higher elevations, and low-medium shrubs (Salix and Alnus spp.) along low-elevation riparian corridors (Viereck et al. 1992). Rocky outcrops provide nesting substrate for members of the cliff-nesting raptor guild: golden eagles, gyrfalcons, and rough-legged hawks; all of which raise nestlings during the short climactic window of the Arctic summer (May-August). The three focal raptor species of this study have similar requirements for nesting substrate and prey type, but differ in their foraging strategy, reproductive/ migratory phenology, and life-history strategy.
Golden eagles are considered to be opportunistic dietary generalists across much of their range (Katzner et al. 2020), and on the Seward Peninsula their primary prey resources include Arctic ground squirrels (Urocitellus parryii), ptarmigan (Lagopus lagopus and L. muta), and waterfowl (Anas and Branta spp.; Herzog et al. 2019). They migrate to the study site in the early spring from the continental US and are the first of the three focal species to initiate egg-laying, requiring a longer interval over which to raise up to two young (Fig. S2). Gyrfalcons are considered to be ptarmigan specialists throughout much of their circumpolar range , although during the breeding season on the Seward Peninsula they consume primarily either ptarmigan or Arctic ground squirrels and generally exhibit higher degrees of dietary plasticity than has been observed in other systems (Robinson et al. 2019). The Alaskan population of gyrfalcons are thought to be year-round residents (Eisaguirre et al. 2016) and typically initiate egg-laying in late April, requiring a hatch-fledge interval of ca. 45 days to raise up to four young ; Fig. S2). Roughlegged hawks specialize on cyclic arvicoline rodents when they are available, and exhibit a strong numerical response (i.e., a change in density) to rodent density across the Arctic (e.g., Fufachev et al. 2019). As such, they exert differential competitive pressure on other members of the raptor guild depending on the stage of the arvicoline rodent cycle. During low arvicoline rodent years on the Seward Peninsula, roughlegged hawks have been observed consuming a variety of alternative prey ranging from shorebirds to hares (Springer 1975)-thus although most pairs are thought to use other breeding grounds or forego reproduction when arvicoline density is low, some pairs may be capable of shifting to other diet sources. Rough-legged hawks migrate to the study site from wintering grounds in the continental US, arriving on their breeding ground later than golden eagles but requiring a shorter interval to raise up to six young (Fig. S2).

Field methods
To locate occupied raptor nests and estimate the age of nestlings, we conducted aerial surveys of > 500 historically occupied raptor cliff-nest sites across the study site (Fig. S1; specific methods in Robinson et al. 2019). We installed motion-activated cameras at gyrfalcon nest sites (hereafter, "nest cameras") during mid-incubation, which captured prey deliveries to nestlings from when they hatched to when nestlings were measured and sampled (see details in Robinson and Prostor 2017). Dietary data from cameras placed in gyrfalcon nests were analyzed in previous studies (Johnson et al. 2020;Robinson et al. 2019) by visually identifying, quantifying, and estimating the biomass of prey deliveries. We visited nest sites midway through the brood-rearing period (golden eagles ca. 35 days, gyrfalcons ca. 25 days, rough-legged hawks ca. 20 days; Fig. S2) to collect blood samples. Red blood cells have been shown to have a half-life of ~ 2 weeks in similarly sized birds (Phillips and Eldridge 2006), thus we expect their stable isotope ratios to reflect this interval (i.e., from when nestlings hatched to when they were measured) with some bias toward more recently consumed prey. We elected to analyze red blood cells for stable isotopes because they have a longer turnover rate relative to plasma and the proximate composition of plasma varies with time since last meal (Adibi and Mercer 1973). We collected blood from each nestling at each targeted nest site and immediately deposited it into a heparinized vacutainer. Samples were stored on ice for < 6 h until they could be separated via centrifugation. After removing the plasma, red blood cells were immediately frozen at -20 °C until they could be processed for stable isotope analysis.
Raptor blood samples were collected between 2014 and 2019, although the sampling effort for each species varied across years (Tables S1 and S2). We use the term "nest" as a sampling unit for this study, which refers to all young within a nest within a year. Thus, diet estimates represent the prey each pair of adults fed to their nestlings during the brood-rearing period in a given year. Importantly, there is also the possibility that adults consume a different diet for self-maintenance-thus, our diet estimates more accurately reflect prey fed to young. Total raptor sample sizes are as follows: golden eagles, 44 nests (59 nestlings; 2014-2019); gyrfalcons, 40 nests (105 nestlings; 2016-2019); roughlegged hawks, 11 nests (25 nestlings; 2018-2019).
Arvicoline rodents exhibit complex population cycles on the Seward Peninsula, and we implemented an opportunistic trapping protocol to collect tissues for stable isotope analysis (see below) that also allowed us to estimate a population index during the 2018 and 2019 field seasons. Briefly, we set grids of 40 baited Museum Special snap-traps in 10 m × 10 m grids and checked them every 24 h. We set these grids in seven different locations across the study site constituting different tundra microhabitat types and revisited the same plots between years to reduce geographic variability in our sampling design. We calculated trap success rates and used them as a metric to generate a relative index of arvicoline rodent abundance across years.
In 2018 and 2019, we collected 219 muscle prey tissue samples representative of seven predetermined categories of the Seward Peninsula raptor diet (Arctic ground squirrels; migratory shorebird and passerine species; waterfowl; hares (Lepus americanus and L. othus); cyclically available arvicoline rodents (Microtus, Lemmus, and Myodes spp.); avian predators (gulls, jaegers, and other raptors); and ptarmigan (Table S3). We targeted > 5 individuals of each species and distributed our sampling efforts across space and time, focusing collection on a time interval contemporaneous with the raptor brood-rearing period. From each prey specimen, we dissected a ~ 1 g sample of muscle tissue (breast for birds and thigh for mammals) and stored samples at − 20 °C until they could be processed for stable isotope analysis. The carbon:nitrogen ratio of consumer and source groups was low (C:N < 4), and we found no significant relationship between d13C and the C:N ratio of prey (p > 0.15 for all prey categories). Thus, we did not extract lipids from red blood cells or muscle tissue samples (following Post et al. 2006).

Laboratory analysis
To prepare samples for stable isotope analysis, we freezedried prey muscle tissue and raptor red blood cell samples for 48 h, then ground them to a fine powder and subsampled 0.05 g into tins. We analyzed carbon and nitrogen stable isotope ratios using a Costech ECS 4010 (Costech Analytical Technologies, California, USA) connected to a Delta V Plus XP Mass Spectrometer (IRMS; Thermo Fischer Scientific) via a Finnigan Conflo III (Thermo Fischer Scientific). Stable isotope ratios are expressed using δ notation (parts per thousand-‰), relative to atmospheric N2 and Vienna PeeDee Belemnite for δ15N and δ13C, respectively. We calibrated our analyses with within-run peptone standards every 10 samples, estimating a precision of ± 0.13 (SD) for nitrogen and ± 0.10‰ for carbon across runs. Laboratory working standard peptone (No. P-7750 meat-based protein. Sigma Chemical Company) was calibrated to NIST standards (REF 8544,8550,8551,8549,8542,USGS 40 and 41) twice annually for quality assurance. Laboratory analyses were performed in the Alaska Stable Isotope Facility at the University of Alaska Fairbanks.

Isotopic niche width
Isotopic niche width in the context of this study was defined as the Standard Ellipse Area (SEA) of a stable isotopic Bayesian ellipse, calculated and visualized using the SIBER package in R (Jackson et al. 2011). Unlike other quantitative methods to characterize isotopic niche width (e.g., convex hull approaches), SIBER implements a Bayesian framework to generate credible intervals of niche width estimates (SEA B ), that are largely unbiased and robust to differences in sample size, although may have more uncertainty associated with sample sizes < 10 (Jackson et al. 2011). Thus, we used the SEA B as a metric for the cross-comparison of isotopic niche space among raptor species in different years of our study, as well as a metric by which to quantify niche overlap using maximum likelihood. Bayesian ellipses were calculated at the nest level using averaged isotopic values for nestlings within the same brood. To evaluate Bayesian ellipses as a niche width estimation method, we compared SIBER estimates to a kernel density approach (rKIN; Eckrich et al. 2020) for the multi-year dataset.

Bayesian stable isotope mixing models
We implemented BSIMMs using the MixSIAR package in R (Stock and Semmens 2016). Input parameters for the models were the raw δ 13 C and δ 15 N values of raptor nestling red blood cells, which were analyzed against the mean and standard deviations of δ 13 C and δ 15 N values for each of the seven prey categories (Fig. 1). We tested for the simultaneous overlap of δ 13 C and δ 15 N for each prey category using ANOVA with post hoc pairwise Tukey's tests for each isotope, after confirming data normality (Shapiro-Wilks p > 0.05) and homoskedasticity to meet the assumptions of the ANOVA.
Stable isotope mixing models require the inclusion of a trophic discrimination factor (TDF; expressed as ∆X = δX consumer -δX sources ) to account for differences in isotope ratios between source and consumer tissue that arise through physiological processes. We applied a TDF calculated for gyrfalcon nestlings in this study system based on nest camera data: ∆ 13 C = 0.89 ± 0.43 (SD); ∆ 15 N = 1.13 ± 0.53 (TDF CAM ; Johnson et al. 2020), which performed better than any TDF available from the literature. To assess the validity of applying Gyrfalcon TDF values to other species, we applied the same gyrfalcon TDF CAM to BSIMM studies of other nestling raptor species in different systems that had nest camera dietary data available for validation (Robinson et al. 2018;Swan et al. 2020) and achieved more accurate diet estimates than using any other available TDF (Johnson 2021). Thus, although this TDF CAM was developed for gyrfalcons specifically, we applied it to other raptors in our system because they are subject to similar environmental factors, and because samples were taken at a similar stage of nestling growth. We evaluated the mixing space using simulated polygons (following Smith et al. 2013) and found all consumer points to fall within the 95% mixing region (Fig. S3).
We applied separate informative priors to the mixing model for each raptor species using published values from other raptor diet studies on the Seward Peninsula (Table S4). We rescaled proportions into Dirichlet hyperparameters before their incorporation into BSIMMs (following Stock and Semmens 2016). We also constructed BSIMMs with uninformative priors (i.e., where each relevant prey category is weighted evenly) and compared diet estimates from informative vs uninformative BSIMMs for each raptor species and prey group. We used Bhattacharyya's Coefficient (BC) as a metric by which to compare posterior density distributions, considering a BC > 0.6 to indicate significant overlap between model estimates (following Bond and Diamond 2011). We tested the sensitivity of BSIMMs to informative priors for each raptor species (Figs S4-S8), concluded based on BC estimates that the priors were not disproportionately driving the models, and elected to use BSIMMs formulated with informative priors for all further inference.  Table S1). Shaded ellipses and boxes correspond to 50%, 75%, and 90% credible intervals of standard ellipse area for each raptor species across the years sampled. Silhouettes correspond to stable iso-tope values from each of seven prey categories, displayed as mean ± 1 SD (sample sizes; Table S1). Raptor isotope values are adjusted using TDF CAM mean values, and the TDF CAM SD is added to the error bar of each prey species. Inset (bottom left) represents numerical estimates of standard ellipse area for each raptor species All BSIMMs were formulated with 100,000 iterations thinned by 25 and a burn-in of 50,000 with three chains, that were considered to have converged when they passed the Geweke and Gelman-Rubin Diagnostics (following Stock and Semmens 2016). We generated proportional diet estimates (with credible intervals) for each raptor nest in each year, treating territory and year as random effects in MixSIAR because nestlings within a brood or cohort were expected to have more similar isotopic values to each other than the population as a whole. We tested the influence of each random effect in an LOOic framework following Stock et al. (2018) to verify this covariate structure represented the top model for each raptor species (Table S5). Mean estimates from mixing models were then used to estimate the degree of individual specialization for each nest.

Individual specialization and niche variation
Individual dietary specialization was calculated for each cohort (i.e., for each species in each year at the nest level) using the E index of individual specialization (Araújo et al. 2008) and produced through the RInSp package in R (Zaccarelli et al. 2013). Briefly, this index compares the diet of each nest to the diets of all other nests within the cohort to determine the average pairwise dissimilarity; values range from 0 (all nests have identical diets) to 1 (all nests have unique diets). We calculated variance around the E index via a jackknifing protocol (Araújo et al. 2008) and compared individual specialization across years for each raptor species using a 95% confidence interval.
We analyzed the relationship between individual specialization (E index) and total isotopic niche width (SEA B ) for each raptor species using generalized linear models with a beta distribution (using betareg; Cribari-Neto and Zeileis 2010). Due to a low number of sample years, rough-legged hawks were excluded from this portion of the analysis. We tested whether the degree of individual variation was greater than expected by chance by comparing the slopes of empirical models to those of null models generated under the assumption that each individual drew randomly from the population-level diet (following Bolnick et al. 2002). To generate null models, we calculated the E index of specialization on proportional diet estimates (n = the number of nests in that cohort) randomly drawn from posterior distributions of BSIMMs generated for each species within each year. We repeated this process 1000 times for each cohort to calculate a null degree of individual specialization (E Null), which was then used to generate null models. Null models may indicate a statistically significant correlation (e.g., in Bolnick et al. 2007). The slopes of the empirical and null models were compared for each species using GLMs to determine statistical support for the NVH.
To further assess support for the NVH in gyrfalcons, we repeated this process with proportional diet estimates from nest camera data (n = 58 nests 2014-2019: Table S6). Diet estimates reflect the first half of the brood-rearing period (nestlings 0-25 days of age) to match the temporal window of stable isotope data. From the camera data, we calculated the E index of individual specialization and TNW from proportional diet estimates of each cohort to create an empirical generalized linear model. To generate a null model, we calculated the E index on random proportional diet estimates drawn from bootstrapped distributions for each cohort.

Validation with nest camera data
There is disagreement regarding the substitution of the isotopic niche for the dietary niche to draw inference on trophic relationships (e.g., Hette-Tronquart 2019). Thus, we tested the validity of isotopic niche estimates in our study system by comparing isotopic niche width (SEA B ) to estimates of total niche width (TNW; Shannon Entropy Diversity Index; Roughgarden 1972) for the years of our gyrfalcon dataset which had contemporaneous stable isotope and nest camera dietary data (n = 20 nests 2016-2019). We calculated proportional diet estimates for each nest from nest camera data (see details in Johnson et al. 2020), and we applied a Bayesian bootstrapping protocol (bayesboot; Bååth 2018) before calculating TNW to make the output of both methods directly comparable. We also tested individual specialization metrics from stable isotopes and nest cameras by analyzing a correlation of the proportional similarity index (PSi; Bolnick et al. 2002) at the nest level.

Results
We analyzed 189 raptor nestling red blood cell samples (Table S1) for δ 13 C and δ 15 N and generated Bayesian ellipses to characterize their isotopic niche space and overlap ( Fig. 1; Table 1). Across years, population niche width (as measured by the SEA B ) differed among species consistent with our assumptions based on their foraging ecology: golden eagles had a significantly larger niche than both gyrfalcons and rough-legged hawks (90% CI-GOEA: 1.56-2.59, GYRF: 0.93-1.55, RLHA: 0.42-1.19; Fig. 1  inset). There was a high degree of isotopic niche overlap between the three species and across years (Table 1), with gyrfalcons and golden eagles overlapping more with each other (61.20%) than either species did with rough-legged hawks (41.37% and 44.16% respectively). SIBER estimates of isotopic niche width were comparable to estimates from a kernel density approach (rKIN; Fig. S9). We further characterized trophic niche dynamics in the raptor guild using proportional diet estimates for each nest from BSIMMs (sample sizes Table S1; Table S7; Fig. S8). The isotopic composition of muscle tissues differed among the seven prey categories: pairwise Tukey's tests indicated no simultaneous overlap in δ 13 C and δ 15 N between any two categories (ANOVA, Tukey's test, p < 0.01) with the exception of hares, which overlapped with arvicoline rodents (δ 13 C p = 0.799; δ 15 N p = 0.977) and Arctic ground squirrels (δ 13 C p = 0.196; δ 15 N p = 0.99). We concluded that each prey category was adequately isotopically differentiable for its use in a BSIMM framework, but that precision in mixing model output may be reduced for the overlapping prey categories. Golden eagle diet was relatively consistent throughout the course of the study, indicating that Arctic ground squirrels were their most important resource, followed by hares and ptarmigan. Modeled gyrfalcon diet indicated a heavy reliance on ptarmigan and a sharp increase in arvicoline rodents in 2019 (to 26% ± 11%). Rough-legged hawks similarly showed a dramatic increase in arvicoline rodent consumption from 2018 (34% ± 11%) to 2019 (73% ± 14%; Table S7). The increased use of arvicoline rodents by gyrfalcons and rough-legged hawks in 2019, relative to 2018, was consistent with arvicoline rodent trapping success data-arvicoline rodents increased 20-fold between 2018 and 2019 (Table S8), indicating that 2019 was a high year in the arvicoline cycle.
We found significant interannual variability in individual specialization for each raptor species. For golden eagles, individual specialization was significantly lowest in 2016 and highest in 2019; gyrfalcons also had the highest degree of individual specialization in 2019; and roughlegged hawks significantly decreased from 2018 to 2019 (Fig. 3). We examined whether individual specialization (E index) increased with population niche width (SEA B ) at a greater rate than expected by chance, consistent with predictions of the NVH. Golden eagles had a lower degree of individual specialization than would be expected by chance, and there was no significant difference between the slopes of the empirical and null models (p = 0.20). Gyrfalcons had a significant positive relationship between individual specialization and isotopic niche width (β = 0.66, R 2 = 0.94, p < 0.01), but the slope of the empirical model was not significantly different from that of the null model (p = 0.14; Fig. 4A). When we further tested this relationship for gyrfalcons using a larger sample of dietary data from nest cameras (n = 58 nests 2014-2019; Table S6), we similarly found a significant positive relationship between individual specialization (E Index) and population niche width (Shannon Index; β = 2.70, R 2 = 0.96, p < 0.01), and the slope was significantly higher than that of the null model (p < 0.01), consistent with the NVH (Fig. 4B).
The isotopic niche (SEA B ) reflected the trophic niche (Shannon Index; based on nest camera data) for the subset of gyrfalcon nestlings for which contemporaneous data existed (Fig. S10). Specifically, the 2019 gyrfalcon subset had a significantly larger niche relative to other years based on the nest camera data (2017: 1.02-1.14; 2018: 0.98-1.14; 2019: 1.46-1.59); stable isotope SEA B data for the same subset of gyrfalcons showed a similar pattern although the 90% credible intervals overlapped. At the nest level, we found a strong positive correlation between individual specialization metrics calculated from stable isotope data and those from nest camera data (Pearson's r = 0.68).

Discussion
We observed a large degree of isotopic niche overlap among golden eagles, gyrfalcons, and rough-legged hawks in the first comparative dietary study of its kind for a guild of raptors. Isotopic niche width generally reflected each species' expected foraging niche, and responses to interannual resource availability varied by species ( Fig. 1; Table 1). In gyrfalcons, indices of individual specialization estimated from BSIMMs were strongly correlated with indices based on camera trap data. Further, we also found support for the NVH in gyrfalcons-individual specialization within the population was higher than expected by chance and positively correlated with trophic niche width on an interannual basis (Fig. 4). Our findings point to niche variation as an important aspect of Arctic raptor trophic ecology, which may reduce niche overlap and promote population stability under dynamic ecological conditions. Characterizing the role of individual specialization in a population's trophic dynamics is important for our understanding of functional ecology because it can greatly affect ecological processes ranging from resource competition to climate change resilience (Bolnick et al. 2003). Under the traditional OFT paradigm, populations are relatively canalized in their foraging behaviors, which maximizes the net rate of energy intake at the cost of heightened intra-specific resource competition and reduced heterogeneity (MacArthur and Pianka 1966). For instance, Jesmer et al. (2020) found that although the dietary niche breadth of large ruminants increases with food limitation, it does so via increased diet breadth of individuals, consistent with OFT, rather than through increased individual specialization. They attribute this pattern to population-level homogeneity in the physiological and cognitive requirements of their study species (moose; Alces alces). Although recent models postulate  Table S1). On upper ellipse plots, different shapes and line types correspond to samples from different years, and darker shaded regions indicate isotopic niche over-lap between years. On lower SEA B density plots, points indicate the mean SEA B estimate and colored boxes correspond to 50%, 75%, and 90% credible intervals respectively. Asterisks indicate species/years that had sample sizes < 5 for SEA B estimates how individual specialists may operate within an OFT framework, they do not fully explain how they arise and are maintained within populations (Svanbäck and Bolnick 2005). Conversely, individual niche specialization (under the NVH paradigm) is predicted to persist and be most prevalent in populations which exhibit intra-specific trait variation, ultimately resulting in different individual foraging strategies/preferences based on phenotype (Bolnick et al. 2003), personality (Toscano et al. 2016), or learned foraging behavior (Araújo et al. 2008). We predicted that the NVH would be supported in the Arctic raptor guild because their high cognitive capacity, parental care, and social learning (e.g., Kitowski 2009) create conditions under which individual differences in foraging strategies may arise. The presence of prey population cycles may further promote the persistence of individual specialization within long-lived predator populations-one foraging strategy may confer high reproductive success in one phase of the cycle but lead to nest failure in the subsequent phase. Long-term studies that simultaneously measure prey abundance and predator diets are needed to establish whether this is a generalizable pattern in response to the arvicoline rodent cycle.
Several underlying traits may have contributed to the varying degrees of individual specialization observed between the three species: phenotypic variation in plumage coloration, territory/pair fidelity, and the cognitive capacity for learning new foraging behaviors. Gyrfalcons are highly phenotypically variable in their plumage coloration (ranging from dark brown to white; Booms et al. 2020) that could affect their detectability by certain prey species and drive individual variation in foraging strategies (e.g.,  Table S1). Error bars represent 95% confidence intervals Fig. 4 The relationship between individual specialization and niche width for gyrfalcons on the Seward Peninsula, Alaska. a Represents E index and SEA B estimates from stable isotope data (n = 40 nests; 2016-2019; Table S1). Panel (b) represents E index and niche width estimates (Shannon entropy diversity index) from nest camera data (n = 58 nests; 2014-2019; Table S5). Each point represents population-level data from each year. In b, the slope of the empirical linear model (green line) is significantly different from that of the null model (gray line), consistent with the NVH Bolnick et al. 2011). Whereas golden eagles and rough-legged hawks migrate south for the winter, gyrfalcon pairs stay on or near their home-ranges year-round (Eisaguirre et al. 2016), which may facilitate learning specialized foraging behaviors specific to their home-ranges. Use of gyrfalcons in the sport of falconry has also revealed they can be trained to hunt prey they would not encounter in the wild (Cade 1982), further highlighting their individual capacity for learning complex foraging behavior. Based on these same criteria, rough-legged hawks were expected to have the lowest degrees of inter-individual variation in foraging strategy: they are more phenotypically similar and as a migratory and semi-nomadic species, they likely have the lowest rates of site/pair fidelity (e.g., Bechard et al. 2020). We further expected rough-legged hawks to operate under the classic OFT paradigm because they may be more physiologically constrained to specialize on certain prey types (i.e., lemmings and voles) than their more generalized counterparts. Golden eagles are long-lived dietary generalists whose larger body size (3-4 × that of the other two focal species) allows them a larger prey base, they exhibit territory and pair fidelity in our study system and have clearly displayed the ability to learn complex foraging behavior on an individual basis (e.g., knocking goats from cliffs and fishing; reviewed in Katzner et al. 2020). We predicted that golden eagles would follow the NVH based on these traits, but instead found that they display lower degrees of individual specialization than would be expected by chance.
Our finding that golden eagles did not fall under the NVH may have an ecological explanation or could be attributed to methodological difficulties associated with constructing BSIMMs for a generalist species. In our study system, prey groups occupy different temporal niches [e.g., strictly diurnal ground squirrels and ptarmigan (Williams et al. 2017;Reierth and Stokkan 1998); nocturnal arvicoline rodents; (Peterson and Batzli 1975)]. Golden eagles may require a broad temporal foraging window to meet the nutritional needs of their growing young, potentially necessitating that each individual or pair within the population broaden their trophic niche as a result. Alternatively, this finding could have an analytical basis: BSIMM estimates for golden eagles were more sensitive to an informative prior (Fig. S5) and we had lower sample sizes in some years compared to gyrfalcons, which may have had a homogenizing effect on modeled inter-individual diets. Furthermore, golden eagles (as dietary opportunists) proved to be generally difficult to analyze in a BSIMM framework simply because they eat so many different prey types. Although we included source values from their main prey categories (Table S3; Herzog et al. 2019), underrepresented groups outside of the mixing space could have influenced fine-scale golden eagle diet estimates and resulting individual specialization metrics. Since gyrfalcons and rough-legged hawks were less sensitive to an informative prior and have prey groups that are better isotopically characterized for this study, we are more confident in our BSIMM diet estimates for those species than the diet estimated for golden eagles. Additionally, although physiological factors such as nutritional stress and variation in growth can influence TDFs and isotopic niche metrics (e.g., Gorokhova 2018), we found that diets (Johnson et al. 2020), niche metrics, and indices of specialization (this study) estimated from stable isotopes were strongly correlated with estimates based on nest camera data.
For the two years for which we have data on arvicoline rodent abundance, we found that inter-specific niche overlap decreased in the high arvicoline year and that individual specialization concurrently expanded in the species with greater dietary flexibility (i.e., gyrfalcons and golden eagles). This suggests that resource abundance may be influencing individual niche variation. The consequences of individual specialization are typically discussed in the context of intra-specific competition (e.g., Bolnick et al. 2005;Xia et al. 2020), but the same patterns may hold true for niche dynamics at an inter-specific level (e.g., Bison et al. 2015). A peak of cyclic arvicoline rodents in tundra ecosystems increases the diversity of the raptor prey base: both directly increasing rodent density but also indirectly increasing the density and reproductive success of their specialist predators which these raptors may also consume (Therrien et al. 2014; e.g., in our system; long-tailed jaegers [Stercorarius longicaudus], short-eared owls [Asio flammeus], short-tailed weasels [Mustela erminea]). BSIMM output indicates that some gyrfalcon pairs remained ptarmigan specialists (constituting up to 89 ± 5% of the diet), whereas others relied more heavily on arvicoline rodents directly (up to 62 ± 5%), and still others capitalized on arvicoline rodent predators (up to 19 ± 5%; Table S9). We thus suspect gyrfalcons in our study system may fall under a "competitive refuge" model (described by Svanbäck and Bolnick 2005), wherein individuals within the population share top-ranked prey (e.g., ptarmigan) but have different rankings for lower tier resources based on efficiency trade-offs associated with their phenotypes or learned foraging behaviors, which may vary on a spatial and temporal basis. Rough-legged hawks displayed a trend opposite to the other two raptor species, reducing inter-individual specialization in 2019 as all individuals were able to capitalize upon their preferred resource, consistent with a traditional OFT response. Collectively, these disparate foraging strategies appeared to result in reduced niche overlap between the three raptor species under higher resource diversity. This highlights the complex roles population cycles may play in tundra food webs, though more long-term studies are needed to determine the extent to which arvicoline cycles influence niche overlap among predator communities.
We provide evidence that individual specialization correlates with population niche width in an Arctic raptor population in agreement with the Niche Variation Hypothesis. Along with directly mitigating resource competition, this niche variation may promote population stability under dynamic ecological conditions by buffering against the loss (or periodic unavailability) of specific habitats or resources (Bolnick et al. 2003). This is particularly pertinent in Arctic tundra ecosystems, which are undergoing rapid change due to anthropogenic global warming that is causing unprecedented and unpredictable effects on ecosystem-level processes (review in Post et al. 2009). Specifically, arvicoline rodent population cycles are collapsing and dampening across Europe and the Arctic (Ims et al. 2008), influencing the population dynamics of their specialist predators (e.g., Fufachev et al. 2019) and illuminating the need to study food web dynamics in systems where these cycles are still intact. Our findings indicate that although arvicoline rodent abundance is directly important to its specialist predators (i.e., rough-legged hawks and some individual gyrfalcons) it also may indirectly widen the available prey base for other members of the predator guild, allowing for reduced resource competition via niche variation. Further work linking individual specialization to reproductive success and territory characteristics will be crucial in furthering our understanding of predator population resilience in a changing Arctic.