Svalbard reindeer winter diets: Long‐term dietary shifts to graminoids in response to a changing climate

Abstract Arctic ecosystems are changing dramatically with warmer and wetter conditions resulting in complex interactions between herbivores and their forage. We investigated how Svalbard reindeer (Rangifer tarandus platyrhynchus) modify their late winter diets in response to long‐term trends and interannual variation in forage availability and accessibility. By reconstructing their diets and foraging niches over a 17‐year period (1995–2012) using serum δ13C and δ15N values, we found strong support for a temporal increase in the proportions of graminoids in the diets with a concurrent decline in the contributions of mosses. This dietary shift corresponds with graminoid abundance increases in the region and was associated with increases in population density, warmer summer temperatures and more frequent rain‐on‐snow (ROS) in winter. In addition, the variance in isotopic niche positions, breadths, and overlaps also supported a temporal shift in the foraging niche and a dietary response to extreme ROS events. Our long‐term study highlights the mechanisms by which winter and summer climate changes cascade through vegetation shifts and herbivore population dynamics to alter the foraging niche of Svalbard reindeer. Although it has been anticipated that climate changes in the Svalbard region of the Arctic would be detrimental to this unique ungulate, our study suggests that environmental change is in a phase where conditions are improving for this subspecies at the northernmost edge of the Rangifer distribution.


| INTRODUC TI ON
Arctic ecosystems are changing in response to a warmer and wetter climate, leading to complex and varied responses by different species (Box et al., 2019;Chapin et al., 2008;Post & Stenseth, 1999). The warming and lengthening of the Arctic growing season has induced plant community shifts, increased primary production and in some cases increased the availability of palatable forage for herbivores (Elmendorf et al., 2012;Mallory et al., 2020;May et al., 2020;Myers-Smith et al., 2020;Post et al., 2009). In contrast, winter warming may have adverse effects on forage availability. Winter precipitation at above-zero temperatures may result in rain-on-snow (ROS) and ice layers that encapsulate vegetation (Andrews, 1996;Putkonen & Roe, 2003;Robinson et al., 1998).
Svalbard reindeer population densities exhibit interannual variation caused by density-dependent (i.e., competition for food resources) and density-independent factors (i.e., climate; Pacyna et al., 2018). In the most productive parts of Svalbard, populations of Svalbard reindeer have increased since the turn of the century (Hansen, Gamelon, et al., 2019;Le Moullec et al., 2019). The positive population trends have been attributed to increased summer temperatures that have led to increased plant productivity and higher forage availability (Hansen et al., 2013;van der Wal & Stien, 2014;Vickers et al., 2016).
Furthermore, a lengthening of the snow-free autumn season allows reindeer easier access to food for longer periods of the year and thus gain more fat reserves before the winter (Albon et al., 2017;Loe et al., 2020). Severe winters have been shown to alter the population structure (demography), with elevated mortality in calves, males, and old-age classes, and poor calf production the following summer (Albon et al., 2017;Peeters et al., 2017). However, the negative effects of difficult winters are reduced by behavioural responses to ROS, such as Svalbard reindeer migrating in search of nearby ice-free foraging sites (Loe et al., 2016;Stien et al., 2010).
Climatic impacts on forage availability is one of the main drivers of Svalbard reindeer population dynamics (Albon et al., 2017).
Optimal foraging theory (OFT; Pyke et al., 1977) predicts that generalist herbivores will shift from an optimal high-quality diet to one that maximizes the quantity of forage when forage resources are reduced or inaccessible. Reindeer winter diets are typically dominated by energy-rich fruticose lichens (Skogland, 1984), but on Svalbard lichens are scarce. Low-quality, energy-poor mosses contribute between 30% and 50% to the winter diets of Svalbard reindeer (Bjørkvoll et al., 2009;van der Wal, 2006;Węgrzyn et al., 2019). However, the low growing stature of mosses is likely to make them inaccessible when iced over during ROS events. Some studies have tested OFT using the winter diets of Svalbard reindeer and have captured a temporal "snapshot" of either the inter-annual dietary variability over a couple of years or the direct response to a rare extreme event (Beumer et al., 2017;Bjørkvoll et al., 2009;Hansen & Aanes, 2012;Hansen, Lorentzen, et al., 2019). For instance, Bjørkvoll et al. (2009) concluded that during winters with ROS (2000ROS ( -2002, Svalbard reindeer forage less on mosses and, instead, prefer plants with erect growth forms, such as graminoids, or plants growing in more exposed habitats, such as the dwarf shrubs, polar willow (Salix polaris) and mountain avens (Dryas octopetala). It is currently unclear what the long-term patterns of Svalbard reindeer winter diets are, in response to the cumulative impacts of increased population density, increased frequency of ROS events, increased summer temperatures and potentially increased forage abundance.
Thus, diets and dietary niches reconstructed from the variations in δ 13 C and δ 15 N may also indicate how inland populations of Svalbard reindeer are and may be capable of adapting their foraging strategy to interannual climate variation and climate change.
The overall objective of this study was to identify how Svalbard reindeer modify their winter diets in response to long-term trends and interannual variation in forage availability over a 17-year period (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012). Svalbard reindeer are the only ungulate in a relatively simple food web where they are subjected to insignificant levels of interspecific competition and predation. This makes them an excellent study system to evaluate herbivore winter foraging ecology in the context of climate change . We collected serum from Svalbard reindeer in late winter, which is the critical period for their winter survival and reconstructed their diets using serum δ 13 C and δ 15 N values. Specifically, we asked whether Svalbard reindeer's late winter (1) serum isotopic values vary over time, (2) dietary composition shows temporal variation, and (3) isotopic niche spaces vary over time in terms of area and/or overlap. In addition to Bayesian stable isotope mixing and isotopic niche space models (Jackson et al., 2011;Parnell et al., 2010), we used linear regressions to explore (4) potential mechanisms that may drive trends in serum isotopic values. Furthermore, we evaluated how serum isotopic values were associated with variation in forage production (summer warming), forage accessibility (ROS), intra-specific competition (population density), and intrinsic drivers like nutritional status (body mass) and status concerning pregnancy. We hypothesized, that graminoids and the deciduous shrub, S. polaris, would increase in the winter diet over time, as these plants were expected to increase in their cover in Svalbard (Hansen et al., 2010b;van der Wal & Brooker, 2004).
Furthermore, we hypothesized that the proportion of mosses in the diet would decline when covered by ice due to ROS, while plants (i.e., graminoids) with erect growth forms growing in exposed habitats would often be accessible above ice layers and increase.

| Study area and species
The study was carried out in Nordenskiöld Land, western Svalbard, in Annual precipitation varies across the archipelago and is 337 mm on the W coast and decreases to 190 mm in the fjord areas (Førland et al., 2011;Johansen et al., 2012;Norwegian Meteorological Institute, 2021a). Between 1971-2017, the average precipitation increased by 4.1% (~4.1 mm) per decade (Svalbard Airport; Hanssen-Bauer et al., 2019). In winters, the frequency of precipitation events during above-zero temperature periods, ROS events, has increased ( Figure S11; Peeters et al., 2019) and simultaneously, the duration of snow-covered days has decreased from an average of 253 days (1976-1997) to 221 (2006Norwegian Meteorological Institute, 2021b).
Svalbard reindeer are found across the non-glaciated parts of the archipelago (Pacyna et al., 2018;Pedersen, 2017), the highest population densities are in Nordenskiöld Land, Edgeøya, and Barentsøya,

| Svalbard reindeer sampling
Reindeer, 3-to 15-year-old females, were caught annually in late winter (late April to early May) using a net between two snowmobiles and restrained without immobilization drugs as part of a long-term capture-mark-recapture study (Milner et al., 2003;Omsjoe et al., 2009;Stien et al., 2002). During capture, ear tags and neck marker straps on each reindeer were used to record the individual's identity, then live body mass was measured and pregnancy was diagnosed (Ropstad et al., 1999). Blood was collected from the jugular vein into plain vacutainers and centrifuged (
To control for this, vegetation was re-sampled in Adventdalen in February 2019 from five randomly selected 4 m 2 plots as no feeding craters were observed due to extreme ground icing. The plant samples collected were graminoids (n = 2), mosses (n = 4), S. polaris (n = 2), and D. octopetala (n = 5). Forbs (n = 4) were collected as they were present in the plots but not present in or collected from the craters sampled in 2013. All plant samples were stored frozen (−20°C) until drying (60°C, 48 h).

| Stable isotope analysis
The dried serum and plant samples were homogenized in a TissueLyser II (Qiagen GmbH) at 25 Hz with 5-mm stainless steel beads (Qiagen GmbH) that were rinsed with 100% ethanol and dried between samples. Powdered serum (0.8-1.2 mg) and forage However, the discrimination of δ 13 C during photosynthesis in C 3 plants has been found to increase in response to increasing concentrations of atmospheric CO 2 (pCO 2 ) and decreasing atmospheric δ 13 C CO2 (Schubert & Jahren, 2012). Thus, we used the equation of Schubert and Jahren (2012) and annual pCO 2 (Tans & Keeling, 2022) to correct all serum and forage samples δ 13 C values to expected 2012 levels. The full data set compiled for the analyses that follow are available from the Dryad Digital Repository .

| Data analyses
All analyses were done in R version 4.05. (R Core Team, 2021) in R Studio (version 1.4.1106).

| Svalbard reindeer diet reconstruction
The proportions of the different forage sources in the diet over time were estimated using Bayesian stable isotope mixing models in simmr (R-package simmr version 0.4.5; Parnell, 2019). The performance of these mixing models is sensitive to the selection of priors and when available, informative priors are recommended (Phillips et al., 2014). Informative priors were available from forage species found in Svalbard Reindeer rumen contents collected in late winter (April/May) during 2000-2002 by Bjørkvoll et al. (2009). The mean percentage of forage species were as follows: graminoids (32%), mosses (26%) and S. polaris (24%), forbs (4%) and D. octopetala (4% ;   Table S1).
Tissue discrimination factors (TDFs) or the change in isotope ratios between source and consumer are often the weakest link in the reconstruction of diets (Bond & Diamond, 2011). In this study, TDFs were estimated using SIDER (version 1.0.0.0), an R package that incorporates reported TDFs, phylogenetic relatedness, foraging ecology, and tissue type into the Bayesian imputation of the most likely TDF values (Healy et al., 2017). The imputed TDFs were ∆ 13 C = 4.51 ± 1.94‰ and ∆ 15 N = 3.36 ± 1.24‰.
Bayesian stable isotope mixing models are also sensitive to the forage source data used. As forage samples were not available for the same years and locations as the reindeer serum samples, the available forage data were evaluated and grouped according to best practice principles (Phillips et al., 2014). That is all used plant sources are included, and they match the consumer's samples spatially, seasonally, and temporally as closely as possible (Table S2).
Four combinations of the source data (Table S3) were tested in the stable isotope mixing models. The simmr models were fitted to the reindeer serum data from each year separately using a burn-in period of 5000 MCMC iterations out of a total of 130,000 iterations of four chains that were thinned to every 50th iteration before analysis (Parnell, 2019;Parnell et al., 2013). Convergence was evaluated using the Brooks-Gelman-Rubin diagnostics (Fernándezi-Marín, 2016). Posterior distributions, posterior predictive values, credible intervals (CIs), correlation matrices and the deviance information criterion (DIC) were all used to check the performance of all the models. The models run with the source data from all three studies combined (Table 1) were found to be the best performing in terms of having the lowest DICs and correlation matrices within acceptable ranges (Figures S1-S6; Tables S4-S7).

| Svalbard reindeer isotopic niche space reconstruction
The isotopic niche spaces of the reindeer were reconstructed for all years using SIBER (Stable Isotope Bayesian Ellipses in R package; version 2.1.5; Jackson et al., 2011). Standard ellipse areas after small sample size correction (SEA C ) have a 95% probability of containing sampled parameters and were plotted with δ 13 C and δ 15 N values in bivariate δ-space. We used Bayesian standard ellipse area (SEA B ) estimates from 10 4 replicates to test for differences between the years through a comparison of the proportion of posterior ellipses (PP) that differed between groups. Differences in SEA B were considered significant if PP ≥0.95. In addition, the Bayesian CIs derived from the SEA B were compared to determine if there were interannual differences in niche breadths (Cheeseman et al., 2021)

| Drivers of the temporal trends in δ 13 C and δ 15 N values of Svalbard reindeer
As the isotopic signatures varied over the years, we, firstly, performed a post hoc test to inspect for significant differences between the years (Tukey's HSD; α = .05). Then, we ran a selection of linear mixed-effects models (LMMs) to explore the underlying drivers of the variation in the response variables δ 13 C and δ 15 N (lme; R-package lme4 version 1.1-27.1; Bates et al., 2015). This modelling approach is recommended as it can account for repeated observations from the same individual as well as spatial and temporal autocorrelation (Zuur et al., 2009). To account for the repeated measurements of individuals between years as well as the environmental variables only differing between years; reindeer identity (ID) and year were fitted as random factors using restricted maximum likelihood. The intrinsic factors of nutritional status and pregnancy that are known to affect the response variables were incorporated into the null models as fixed effect predictor variables. As a proxy for nutritional status, Population density and ROS were ln-transformed before being fitted as predictor variables in the model for δ 13 C. Six biologically relevant candidate models for variation in δ 13 C and δ 15 N were developed based on the abovementioned predictor variables (Table S12). Models were compared using the conditional Akaike

| RE SULTS
Isotopic ratios from serum exhibited interannual variability with δ 13 C predominantly showing a negative temporal trend (estimates slope,

| Reconstructed diets of Svalbard reindeer based on mixing models
Dietary compositions reconstructed by the selected stable isotope mixing model showed interannual variation ( Figure 3; Table S7).

| Svalbard reindeer isotopic niche widths
The isotopic niches reconstructed from δ 13 C and δ 15 N values of Svalbard reindeer showed isotopic drift over time with high interannual variation in overlaps ( Figure S7; Table S10). In addition, all year specific ellipses (SEA B ) differed significantly in size from at least one other year ( Figure S8; Year δ 15 N (‰)

(b)
between 2007 and 2012. Extreme ROS events of greater than 60 mm of precipitation were observed in 1996 and 2012, and the overlap between these two years was 57% (Figure 4; Table S10). In contrast, the overlap between these extreme ROS years with the years prior (1995 and 2011), when no ROS was observed, was 0% between 1995 and 1996, and 71% between 2011 and 2012.

| Drivers of the temporal trends in δ 13 C and δ 15 N values of Svalbard reindeer
Among the six candidate models describing the relationship between intrinsic and extrinsic factors and δ 13 C, we identified three models with ΔcAIC ≤2 (Table S11). Of these models, the full model that included ROS, July average temperature, population density and body mass was selected as a parsimonious model (∆cAIC = 1.75; Table S11). The total variance explained by the fixed effects (R 2

LMM(m)
) in the selected δ 13 C model was only 27%, while 64% of the variance was explained by both the fixed and random effects (R 2 LMM(c) ; Table S11). In the selected model, δ 13 C decreased as ROS, July average temperature, population density and body mass increased, but the effect of July average temperature was not significantly different from zero ( Table 2). The fixed effect predictor variables explained some of the temporal trend in average δ 13 C, but not all. The correlation between year and average δ 13 C was r = −.75 in the raw data and was reduced to r = −.41 when δ 13 C levels were estimated  Four of the six candidate models identified to explain the temporal variance in δ 15 N in response to the intrinsic and extrinsic factors were relevant as their ΔcAIC ≤2 (Table S11). As with the δ 13 C, the full model that included ROS, July average temperature, population density, body mass and pregnancy was among these parsimonious models. The total variance in δ 15 N explained by the full model (R 2 LMM(c) ) was 78% and by the fixed effects alone (R 2 LMM(m) ) alone was 34% (Table S11). Increases in ROS, July average temperature, population density and body mass were all associated with increased δ 15 N while pregnancy was associated with reduced δ 15 N ( Table 2). The parameter estimate for the effect of population density was not significantly different from zero ( Table 2). As for δ 13 C, the fixed-effect predictor variables explained some of the temporal trend in average δ 15 N, but not all. The correlation between year and average δ 15 N was r = .58 in the raw data and was reduced to r = .32 when δ 15 N levels were estimated by the annual average residual values, corrected for estimated fixed effects. Of all the predictors in the selected models for both δ 13 C and δ 15 N, only population density and July average temperature showed a marked trend over the study period and are therefore driving some of the temporal trends in δ 13 C and δ 15 N.

| DISCUSS ION
In this study, we sampled the serum of Svalbard reindeer between 1995 and 2012 to identify how the winter diets of a High Arctic herbivore have changed in response to long-term climate trends and interannual variation, mostly depicted by the "normalization" of ROS events, increased summer temperatures and increased population density. Preceding work has reported short-term variations in the summer and winter diets of Svalbard reindeer (Bjørkvoll et al., 2009;Hansen, Lorentzen, et al., 2019;Pacyna et al., 2018;Zhao et al., 2019). We demonstrate changes in the form of long-term trends in the winter diets.

| Temporal increase in graminoid consumption by Svalbard reindeer
We expected reindeer winter diets to follow OFT and thus the changes in the abundance and availability of forage species.
Consequently, dietary shifts would be expected to parallel the observed increases in greening (Vickers et al., 2016), biomass production, and abundance in graminoids (Hansen et al., 2010b;van der Wal & Brooker, 2004) on Svalbard and shrubification across the Arctic (Myers-Smith et al., 2011. Accordingly, we found strong support for increased proportions of graminoids (44% → 58%) and decreased proportions of mosses (26% → 18%) in the reindeer late winter diets over time, but contrary to our expectations, the proportion of S. polaris did not increase (Figure 3). We also found that the reconstructed isotopic niches of Svalbard reindeer shifted, expanded, and contracted over the years with different ROS conditions in a manner suggesting dietary adaptation to more frequently occurring ROS and an increasing abundance of graminoids. Furthermore, the data suggested that the isotopic niche metrics differed in response to ROS in the early years of our study period in comparison to the later years. Lastly, we found that population density together with climate, that is, summer air temperature and ROS, were important drivers of serum isotopic variability highlighting the important role of both density-dependent effects and the changing climate for the physiology and realized trophic niche of these herbivores.
Similar to our findings climatic patterns, such as the broad scale Arctic Oscillation, have been found to influence caribou and reindeer body condition and subsequently population dynamics through TA B L E 2 Parameter estimates (β ± SE) from the selected linear mixed-effects models assessing the effects of extrinsic-rain-on-snow (ROS), July average temperature and population density-and intrinsic-body mass and pregnancy (only for δ 15 N)-predictors on δ 13 C and δ 15 N values of female Svalbard reindeer. Female identity (ID) and year were included as random intercept effects and natural logarithmic (ln) transformation was used for ROS and population density in the δ 13 C model effects on the quality of the summer-ranges Svalbard (Aanes et al., 2002), in Alaska (Joly et al., 2011) and in northern Canada .
Svalbard reindeer are not the only Rangifer spp. whose winter diets have shifted toward more graminoids. A similar trend was observed as caribou in the Western Arctic Herd, Alaska consumed more graminoids and fewer lichens when the abundance of lichens in the foraging habitat declined (Joly et al., 2007). However, these findings may not be generalizable across the Arctic, as there is variation in how well the different Rangifer subspecies can adapt their diets to changing foraging ranges and forage availability . Graminoids are grazing tolerant and adapted to defoliation with basal meristems located close to the soil surface and indeterminate growth allowing continued biomass production after leaf removal by grazers (Barthelemy et al., 2015(Barthelemy et al., , 2019Welker & Briske, 1992). They are also relatively tolerant to trampling while mosses and dwarf shrubs are more sensitive (Egelkraut et al., 2020).
Consequently, graminoids benefit both from increased tundra soil temperature and increased soil nutrient availability generated by higher summer temperatures, a reduction in the insulating shrub and moss layer, and more nutrients released from reindeer faeces and urine (Aerts et al., 2006;Barthelemy et al., 2019;Egelkraut et al., 2020;Marino et al., 2014;Milner et al., 2018;Parsons et al., 1995;Sjögersten et al., 2012;van der Wal et al., 2004;van der Wal & Brooker, 2004 (Bjerke et al., 2017(Bjerke et al., , 2018) and due to their erect growth form (Parsons et al., 1995), they may protrude above the ice layers thus remaining accessible to herbivores. Indeed, on Svalbard, increased graminoid and decreased moss abundance correlates positively with Svalbard reindeer density and increased grazing pressure (Hansen et al., 2010b;van der Wal & Brooker, 2004).
Contrary to our prediction, the proportion of S. polaris was stable with slight interannual variations in the reindeer winter diets between 1995 and 2012 ( Figure 3). These results were surprising, as S. polaris is an abundant species in exposed habitats (ridges and heaths; Le van der Wal & Stien, 2014) and even when encased in ice or exposed to erratic winter temperatures remains undamaged (Le Moullec et al., 2020).  (Bjørkvoll et al., 2009;Gaare et al., 1977;Staaland, 1984). SIA only detects the contributions of δ 13 C and δ 15 N from digested and assimilated forage, and may, therefore, underestimate the use of the poorly digestible plants (Ben-David et al., 2001). Nonetheless, there appears to be no change in the contribution of S. polaris to Svalbard reindeer diets. Thus it may be that the prostrate woody stems and the senescent leaves on the ground are ice-encased and inaccessible while in the ice-free exposed habitats the reindeer are actively selecting against S. polaris in favour of higher-quality graminoids. These patterns in Svalbard reindeer diets follow the premises of OFT under all scenarios and point toward increasing proportions of graminoids in their winter diets.
This suggests that the winter foraging niche of Svalbard reindeer has been changing for decades and not only in response to an increased frequency in ROS events (Beumer et al., 2017

| Extrinsic drivers of the shifts in Svalbard reindeer serum isotope values and diets
Preceding work has reported contrasting changes in caribou δ 13 C and δ 15 N values during winters with body mass declines (Barboza et al., 2020;Gustine et al., 2014;Parker et al., 2005) and fat and protein reserves are mobilized in response to declining energy and N intake (Mänttäri et al., 2013;Reimers, 1984). In the current study, a low body mass was on average associated with enriched δ 13 C values and depleted δ 15 N values, although the effect was small (Table 2).
In contrast, in 1996 when an extreme ROS event resulted in severe starvation and relatively high mortalities (Milner et al., 2003) we observed low body masses, serum δ 13 C depletion and δ 15 N enrichment ( Figure S9; Figure 2). Still, our results support the findings of previous studies, in that they suggest that intrinsic processes, such as pregnancy, and nutritional status (mobilization of body stores) have small contrasting effects on serum δ 13 C and δ 15 N values. Dietary supplies of energy and N are the largest contributors to the stable isotope pool in serum (Ben-David & Flaherty, 2012;Fry, 2006;Rogers et al., 2015Rogers et al., , 2020Stanek et al., 2017Stanek et al., , 2019, consequently, shifts in the serum δ 13 C and δ 15 N values are most likely due to increased utilization of δ 15 N enriched and δ 13 C depleted graminoids in late winter. Ultimately, we propose that the main extrinsic environmental drivers of Svalbard reindeer serum δ 13 C and δ 15 N valuespopulation density, summer temperature, and ROS-have acted collectively to increase the utilization of graminoids in later winter.

ACK N OWLED G M ENTS
We thank the Governor of Svalbard for permission to undertake the research. We thank our collaborators and field assistants from the Svalbard reindeer study system. We are especially grateful to the logistical and technical staff at the University Centre in Svalbard (UNIS) for supporting the field campaigns. The data collection would not have been possible without the contribution of numerous field assistants, including veterinary students from the Norwegian School of Veterinary Science. Matt Odell prepared the blood samples for isotope analysis that were run by Matt Rogers in the UAA SIL. Bart Peeters from NTNU provided the ROS data.

FU N D I N G I N FO R M ATI O N
The work was supported mainly by grants from U.K. Natural

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data and R scripts that support the findings of this study are openly available in Dryad Digital Repository at https://doi. org/10.5061/dryad.ghx3f fbs7 . The R scripts are available in Zenodo at https://doi.org/10.5281/zenodo.7049837.

Reindeer sample collection met all guidelines of the Research
Council of Norway and were collected with the permission of the