# Geochemical landscapes as drivers of wildlife reproductive success: insights from a high-Arctic ecosystem

Floris M. van Beest, Niels Martin Schmidt, Lærke Stewart, Lars H. Hansen, Anders Michelsen, Jesper Mosbacher, Hugo Gilbert, Gaël Le Roux, Sophia V. Hansson

► To cite this version:

Floris M. van Beest, Niels Martin Schmidt, Lærke Stewart, Lars H. Hansen, Anders Michelsen, et al.. Geochemical landscapes as drivers of wildlife reproductive success: insights from a high-Arctic ecosystem. Science of the Total Environment, 2023, 903, pp.166567. 10.1016/j.scitotenv.2023.166567. hal-04246451

HAL Id: hal-04246451

https://hal.science/hal-04246451v1

Submitted on 17 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
1 **Geochemical landscapes as drivers of wildlife reproductive success:**

2 **insights from a high-Arctic ecosystem**

3

4 Floris M. van Beest<sup>a,b,*</sup>, Niels Martin Schmidt<sup>a,b</sup>, Lærke Stewart<sup>c</sup>, Lars H. Hansen<sup>a,b</sup>, Anders

5 Michelsen<sup>d</sup>, Jesper B. Mosbacher<sup>e</sup>, Hugo Gilbert<sup>f</sup>, Gaël Le Roux<sup>f</sup>, Sophia V. Hansson,<sup>a,f</sup>

6

7 <sup>a</sup> Department of Ecoscience, Aarhus University, Frederiksborgvej 399, 4000 Roskilde, Denmark

8 <sup>b</sup> Arctic Research Centre, Aarhus University, Ny Munkegade 116, 8000 Aarhus C, Denmark

9 <sup>c</sup> Department of Natural Sciences and Environmental Health, University of South-Eastern

10 Norway, Gullbringvegen 36, 3800 Bø, Norway

11 <sup>d</sup> Department of Biology, University of Copenhagen, Universitetsparken 15, 2100

12 Copenhagen, Denmark

13 <sup>e</sup> Norwegian Polar Institute, Framsenteret, 9296 Tromsø, Norway

14 <sup>f</sup> Laboratoire Ecologie Fonctionnelle et Environnement (UMR-5245), CNRS, Université de

15 Toulouse, Avenue de l'Agrobiopole, 31326 Castanet Tolosan, France

16
17 Abstract

The bioavailability of essential and non-essential elements in vegetation is expected to
influence the performance of free-ranging terrestrial herbivores. However, attempts to relate
the use of geochemical landscapes by animal populations directly to reproductive output are
currently lacking. Here we measured concentrations of 14 essential and non-essential
elements in soil and vegetation samples collected in the Zackenberg valley, northeast
Greenland, and linked these to environmental conditions to spatially predict and map
geochemical landscapes. We then used long-term (1996-2021) survey data of muskoxen
(*Ovibos moschatus*) to quantify annual variation in the relative use of essential and non-
essential elements in vegetated sites and their relationship to calf recruitment the following
year. Results showed that the relative use of the geochemical landscape by muskoxen varied
substantially between years and differed among elements. Selection for vegetated sites with
higher levels of the essential elements N, Cu, Se, and Mo was positively linked to annual calf
recruitment. In contrast, selection for vegetated sites with higher concentrations of the non-
essential elements As and Pb was negatively correlated to annual calf recruitment. Based on
the concentrations measured in our study, we found no apparent associations between annual
calf recruitment and levels of C, Mn, Co, Zn, Cd, Ba, Hg, and C:N ratio in the vegetation. We
conclude that the spatial distribution and access to essential and non-essential (trace)
elements are important drivers of reproductive output in muskoxen, which may also apply to
other wildlife populations. The value of geochemical landscapes to assess habitat-
performance relationships is likely to increase under future environmental change.

38

39 **Keywords:** Arctic, Ecogeochemistry, Herbivores, Plant quality, Zoogeochemistry

40

41
## 1. Introduction

How free-ranging animals exploit the spatial and temporal variability in habitat and resources available to them can have important implications for their survival and reproduction (e.g., McLoughlin et al., 2007). Quantifying and understanding so called habitat-performance relationships (Gaillard et al., 2010; Matthiopoulos et al., 2015) is critical for a range of topics in environmental management including assessments of climate and land-use change impacts (Verniest et al., 2022). Habitat-performance relationships can be estimated by linking differences in animal performance to variation in use of resource abundance or quality (Felton et al., 2020; Ofstad et al., 2020). The quality of forage resources is, however, often loosely characterized e.g., by assuming a direct relationship between consumer density and prey abundance (in the case of carnivores: Mosser et al., 2009) or by measuring digestibility of plant material (in the case of herbivores: Proffitt et al., 2016). The elemental composition of forage resources is rarely considered explicitly when estimating habitat-performance relationships. This is surprising given that the quality of forage resources is multidimensional and complex, and animals must balance and regulate the intake of multiple essential and non-essential elements to optimize maintenance, growth, and reproduction (Boersma and Elser, 2006; Sperfeld et al., 2017).

The fields of ecogeochemistry and ecological stoichiometry have made great progress in understanding how different biochemical markers relate to the health status of animals as well as ecosystem structure and stability (Elser, 2006; Trueman et al., 2016). Indeed, it is now well established that a balanced intake of macro elements such as nitrogen (N) and phosphorus (P) can positively affect growth (Plath and Boersma, 2001) and reproduction (Durso et al., 2020) of different animal species. Similarly, micro elements i.e., elements required in trace amounts for normal physiological function such as copper (Cu), selenium (Se), and zinc (Zn), are critical to many biochemical pathways and thus considered essential
to life in most organisms (Feinberg et al., 2021; Kaspari and Powers, 2016). Adverse effects may arise, however, if homeostasis is disrupted i.e., when reserves of essential elements are outside required physiological ranges (Bhattacharya et al., 2016; Kabata-Pendias and Szteke, 2015). In addition, some trace elements, such as arsenic (As), cadmium (Cd), lead (Pb), and mercury (Hg) hold no physiological function and can suppress individual health already at low concentration thresholds (Burger, 2008; Massányi et al., 2020).

A range of processes such as weathering and atmospheric deposition combined with gradients in climatic conditions and anthropogenic disturbance can cause substantial heterogeneity in the bioavailability of chemical elements in natural environments (Brown et al., 1999). For terrestrial herbivores, the intake of essential and non-essential elements is primarily influenced by their bioavailability in vegetation, the underlying soils, and plant-specific uptake rates. Thus, the way terrestrial herbivores use habitat and forage resources within the heterogeneous geochemical landscapes they occupy should influence the health and performance of individual animals and, subsequently, population dynamics (Flueck, 1994; French et al., 2017). While assessments of elemental status in free-ranging herbivores and the linkages to wildlife health are rapidly increasing (Jutha et al., 2022; Mosbacher et al., 2022; Rioux et al., 2022; Squadrone et al., 2022), the importance of the underlying mechanism, namely variability in consumer habitat or resource use, is rarely considered. In fact, only few studies have investigated how consumers' habitat use varies as a function of the bioavailability of chemical elements across natural landscapes (Balluffi-Fry et al., 2020; Leroux et al., 2017; Rizzuto et al., 2021; Sach et al., 2020) though the potential consequences for individual or population-level performance were not quantified.

Our objective was to assess if spatiotemporal variation in the use of forage resources within multidimensional geochemical landscapes is an important factor in habitat-performance relationships of wildlife. Specifically, and using the largest herbivore roaming
the Arctic tundra, the muskox (*Ovibos moschatus*) as a study organism, we quantified and mapped the distribution of 14 essential and non-essential elements in soils and vegetation across the Zackenberg valley in northeast Greenland. We then employed long-term location and demography time series of muskoxen to estimate the relative use of forage resources with divergent elemental concentrations and potential linkages to annual variation in calf recruitment. Given their importance in physiological function, we expected the relative use of vegetated sites with increased levels of essential elements to positively correlate with annual calf recruitment (Hypothesis 1), especially for trace elements Se and Cu as deficiencies of these elements are known to limit reproductive success of muskoxen (Mosbacher et al., 2022; Rombach et al., 2002) and other herbivores (Flueck et al., 2012; Michaluk and Kochman, 2007). In contrast, we expected the relative use of vegetated sites with increased concentrations of contaminants (e.g., As, Cd, Hg, Pb) to have a negative effect on annual calf recruitment (Hypothesis 2).

105

## 2. Material and Methods

107

### 2.1 Study area

All data presented here were collected in the Zackenberg valley, northeast Greenland (74°28 N, 21°33 W; Fig. 1). The climate in the region is high-Arctic with a mean annual ambient temperature of ca. -9°C, peaking in July with a mean of 6°C and lowest in February/March with a mean of -20°C. The elevational gradient in the valley ranges from sea-level to 1200m on surrounding hilltops with steep slopes (>20°) in some areas.

Geomorphological landforms in the valley include exposed bedrock, alluvial, fluvial and glacial landforms (Cable et al., 2018). Soil water content ranges from dry conditions at hilltops to wet conditions in the lower valley. Soils are dominated by non-calcareous sandy
sediments and are slightly acidic to neutral (pH range 5.0–7.2), with pH generally increasing
with soil depth (Elberling et al., 2008). Barren areas or poorly vegetated sites are common
above 600m in elevation, while the rest of the valley is mainly covered by a mosaic of tundra
vegetation types of varying productivity including *Cassiope* and *Dryas* heath, fen, grassland
and *Salix* snowbed (Elberling et al., 2008; Stewart et al., 2018). The growing season typically
extends from early June to late August and snow typically covers the ground from October to
May, though interannual variation in snow accumulation and duration is considerable
(Pedersen et al., 2016). Muskoxen are the only large herbivores in northeast Greenland and
are surveyed every summer (Jul-Aug) in Zackenberg valley with densities ranging between
0.5 and 3.0 individuals km⁻² since 1996 (Schmidt et al., 2015). The impacts of human
disturbance and predation by Arctic wolves (*Canis lupus arctos*) and polar bears (*Ursus
maritimus*) on muskox calf recruitment are considered negligible in this area (Schmidt et al.,
2015). Muskoxen at Zackenberg do exhibit anti-predator behaviors, such as forming
defensive circles, but this occurs at such low frequency and of such short duration that it is
unlikely to bias habitat use estimates.

132

## 2.2 *Soil and vegetation sampling*

Between 22–30 of July 2021, soil and plant material were collected from within 50
permanent field plots (Fig 1) that are part of a tundra vegetation monitoring study in the
Zackenberg valley (Stewart et al., 2018). Field plots were selected (from a total of 200 plots)
using a random stratified sampling approach to obtain balanced samples of material across
the elevation gradient (0–600m) and the five dominant vegetation types in the valley
(*Cassiope* and *Dryas* heath, fen, grassland, and *Salix* snowbed). In each field plot, samples of
the dominant plant species or, in case of fen and grassland, functional group (graminoids)
were collected (leaves and stems, no roots) using nitrile gloves. All vegetation was collected

134

135

136

137

138

139

140

141
142 by mimicking muskox grazing i.e., grabbing and pulling a handful of plant material. Plant
143 material (min. ~2g wet weight) was subsequently transferred into paper bags. Soil was
144 collected in each plot by mixing and homogenizing three sub-samples from the A-horizon
145 (i.e., the upper 0–10 cm layer) using nitrile gloves and a pre-cleaned stainless-steel garden
146 shovel. Soil material (~10g wet weight) was transferred into polyethylene zip-lock bags.
147 Fourteen of the 50 field plots were sampled more intensively by collecting three soil samples
148 and three plant samples from within 1 m² to assess and incorporate plot-based variability in
149 element concentrations. In total, 78 soil samples and 78 plant samples were collected from
150 within the field plots. All soil and plant material were stored frozen until sample processing
151 and chemical analyses.

152

153 2.3 Sample processing and chemical analyses

154 Soil samples were dried in an oven at 35°C for 5 days under constant air flow. While
155 we acknowledge that oven drying is not the ideal method, doing so was necessary due to
156 logistical constraints in the field. However, oven-drying soil samples at low temperatures
157 limits the loss of volatile elements (Roos-Barraclough et al., 2002). Once dried, soil samples
158 were sieved through a 200µm nylon mesh to separate soil particles from small gravel and
159 organic matter such as root fragments. Compact samples were manually ground beforehand
160 with an agate mortar. Between each sample, all sieving material (mesh, container, mortar)
161 was rinsed twice, first with regular tap water, then with mQ H₂O (18.2 MΩ cm⁻¹) and finally
162 disinfecting with 70% ethanol. Vegetation samples were lyophilized using a Christ Alpha 1–2
163 freeze-drier and ground using a laboratory mill Retsch® Mixer MILL MM 400, at 30Hz for
164 120s. All material was cleaned between each sample, i.e., rinsed under tap water then rinsed
165 with mQ H₂O and finally disinfecting with 70% ethanol.
For wet chemistry digestion of both vegetation and soil samples, approximately 150mg of dried and ground sample was placed in pre-cleaned (2% HNO₃ and mQ H₂O) 50mL digitubes to which 3mL of ultrapure HNO₃ (Fisher Chemical, Nitric Acid HNO₃, Optima Grade, 67-69%) were added. The tubes were then sealed with airtight caps and placed on a hot plate at 90°C overnight (~12 hours). After allowing the digest to cool to room temperature, the digested solutions were transferred to 50mL Falcon tubes and diluted with mQ H₂O to reach a mother solution (MS) of 45mL. Prior to ICP-MS analysis, the MS were further diluted by transferring a 5 mL aliquot to a 15mL Falcon tube to which mQ H₂O were added up to the 10 mL mark. For quality assurance and accuracy, blanks and certified reference materials (CRMs: IAEA 336, NIST-1515, BCR 142, DORM-4, and DOLT-5) were included with every 10 samples and processed the same way as the vegetation and soil samples (see Table S1 for details). All blanks were >10x lower than the average of the total measurement of each element, and the recovery of each element, for all five CRMs, fell within ± 20 % of the certified value, with the exception of Ba in IAEA 366 (77.3% recovery) and for Pb in DORM-4 (75.4% recovery). All sample preparations were performed under clean laboratory conditions (Clean room class 10000, ISO 7) using acid cleaned lab ware. All concentrations (µg g⁻¹ dry weight (d.w.)) for the elements Mn, Co, Cu, Zn, As, Se, Mo, Cd, Ba, Hg, Pb were above the detection limit (here set to 3 x standard deviation of blank for each element, Table S1) and obtained using a Triple Quadrupole Inductively Coupled Plasma Mass Spectrometer (Thermo scientific, iCap TQ-ICP-MS) at the ICP-MS platform of the Observatoire Midi-Pyrénées, Toulouse, France.

Following the trace element analyses, soil and vegetation samples were also measured for %C and %N and subsequently used to calculate C:N ratios. Subsamples were weighed (soil 10–12 mg dry weight, plants 3–4 mg dry weight) and packed into tin capsules. All C and N analyses were conducted at the University of Copenhagen, Denmark, using an
191 Isoprime isotope ratio mass spectrometer (Isoprime Ltd, Cheadle Hulme, Stockport, UK)
192 coupled to a CN elemental analyser (Eurovector, Milan, Italy) with continuous flow, with a
193 front furnace operating at 1050°C and a reduction furnace operating at 700°C. For quality
194 assurance, certified reference material (CRM NIST 1547 Peach leaves material, National
195 Institute of Standards & Technology) was included with every 12 samples and processed the
196 same way as the vegetation and soil samples. Certified values of the reference material had a
197 mean and standard deviation of N%: 2.965 ±0.054 and C% 46.01 ± 1.52 and were measured
198 at N%: 2.939± 0.080 and C% 46.67 ±1.89 in this study. The detection limits of the elements
199 were 0.015% for N and 0.100% for C.

200

201 *2.4 Quantifying geochemical landscapes*

202 Plot-based measurements of %C, % N and C:N ratio as well as of trace element
203 concentrations (µg g⁻¹ d.w.) (Fig. 2) were used to generate maps of their predicted
204 distribution and potential bioavailability across the soils and vegetation types of the
205 Zackenberg valley (i.e., geochemical landscapes). To do so, we used stoichiometric
206 distribution models (StDMs) as proposed by Leroux et al., (2017). StDMs generate spatially
207 explicit predictions of elemental patterns across the landscape based on statistical relationship
208 between element concentrations/ratio and environmental parameters (e.g., elevation,
209 vegetation cover, soil moisture, etc.). Such maps can then be used to link consumers’
210 distribution or habitat use to element concentrations across natural landscapes (e.g., Sach et
211 al., 2020) as described in section 2.5.

212 StDMs were computed using generalized linear models (GLMs) and constructed for
213 each element (n=14) and for both soil and vegetation samples separately (i.e., a total of 28
214 models). To improve normality of the data, concentrations of trace elements were log-
215 transformed while C, N, and C:N ratio (%) were logit-transformed, and all were fitted as
response variables using a Gaussian error distribution. A total of nine predictor covariates were available as raster layers with a spatial resolution of 1m (Fig. S1) including soil moisture (index), spring snow cover (%), ambient temperature (°C), solar radiation (KWH m⁻²), elevation (m), slope (°), hill shade (unitless), geomorphological landforms (7-class factor), and vegetation type (7-class factor). Because the vegetation type raster layer was produced using several of the candidate predictor variables (see: Stewart et al., 2018 for details), we were unable to use all predictor variables in the same model. Hence, to avoid issues of circularity, predictor variables for the soil analyses included: geomorphology landforms, elevation, slope, spring snow cover and soil moisture. To construct geochemical landscapes based on vegetation data, we considered the predictor variables: vegetation type, hill shade and the element concentration or percentage in the soil. For an overview of how environmental variables contribute to the clustering and spatial distribution of all elements considered see Fig. S2. For the vegetation models, the percent cover of each target plant species or, in cases of fen and grassland, functional group (graminoids), within the sampling plots was used as a weighting variable. Sampling date (Julian day) was included as a predictor variable in all StDMs to incorporate any natural changes in element concentration over time. All continuous variables were included as a second order polynomial to allow for non-linearity in the data.

To assess the predictive performance of the StDMs, we ran each element specific GLM 100 times and for each iteration split the data into a new training (75%) and testing (25%) set for which we calculated the MSE (using the testing data) and R² (using the training data). We then compared the distribution of MSE and R² of the testing and training data to those generated by 100 models where the response variable of the original data was randomly resampled. We considered the performance of the StDMs to be satisfactory as the prediction error was always lower and the amount of variation explained was always higher than for the
models based on randomized datasets (see Figs. S3 and S4). The output of the element-specific StDMs (Figs. S5-S18) were used to produce spatially explicit maps of each element in soil (Fig. S19) and vegetation (Fig. 3) across the study area.

## 2.5 Relative use of geochemical landscapes

To assess how muskoxen exploit the geochemical landscapes they inhabit, we used resource selection functions (RSFs), which are statistical models that compare environmental conditions encountered at used locations to those at available locations (Boyce and McDonald, 1999). The output of RSFs provides estimates of the relative probability of use of vegetation classes or sites along an environmental gradient and can be fitted with logistic regression with used (1) and available (0) locations (Manly et al., 2002). In our RSFs, used locations were geographical coordinates of muskox observations (n=3142) collected annually during summer (Jun-Sep; Fig. 1C) as part of a long-term monitoring program at the Zackenberg valley initiated in 1996 (Schmidt et al., 2015). Only vegetated areas were considered in our analyses; thus, we excluded locations of muskoxen on barren sites or exposed bedrock (n=723). Available locations were randomly generated for each year by sampling locations from within the spatial domain of the geochemical landscapes (Fig. 3).

We used a 1:10 used:available location sampling scheme as it yielded stable model estimates and adequately captured the distribution of available conditions within the landscape, which are important aspects for robust model inference (Northrup et al., 2013).

Predicted concentrations or percentages of elements at used and available locations were extracted from the vegetation-based geochemical landscape maps (Fig. 3) and fitted as fixed-effect explanatory variables in the RSFs using an interaction with the categorical variable Year to obtain annual coefficients of the relative probability of use of each element. Because of multicollinearity among the different chemical elements at used locations
266 (Spearman rank correlations ($r_s$) of >0.6 and <-0.6 for some element combinations, Fig. S20),
267 separate RSFs were fitted for each element. Values of the explanatory variables were scaled
268 prior to analyses to ensure RSF coefficients across years and elements were comparable. In
269 all RSFs, we also fitted Year as a random intercept to account for unbalanced data.

270 Because no evidence exists that muskoxen, or other herbivorous animals, can sense
271 the elemental composition in plant material, they are unlikely to actively select forage
272 resources based on the bioavailability of different chemical elements. Instead, plant biomass
273 and composition and thus vegetation type is likely a more realistic driver of muskox habitat
274 use in this area as plant productivity and biomass varies between vegetation classes
275 (Tomassini et al., 2019). Therefore, and for comparative purposes, we fitted a classic RSF
276 using the interaction between vegetation type and Year as the independent variable.

277 Finally, the predictive success of all RSFs was evaluated using the *k*-fold cross-
278 validation procedure as proposed by Boyce et al. (2002) and Bastille-Rousseau and
279 Wittemyer (2019). For this, we calculated cross-validated Spearman rank correlations ($r_s$)
280 between ten RSF-bin ranks and five test-training sets. We repeated this procedure 10 times to
281 determine whether the $r_s$ differed from random.

282

283 *2.6 Temporal variation in muskox calf recruitment and the linkages to the use of geochemical* 
284 landscapes

285 To examine whether inter-annual variation in the relative use of geochemical
286 landscapes influences calf recruitment, we related the output of the RSFs (i.e., annual
287 selection coefficients for essential and non-essential elements) to muskox demography time
288 series [updated from Schmidt et al., (2015)]. For each year during the period 1996–2021, we
289 calculated calf recruitment rates as the percentage of calves (<6 months old) per 100 adult
290 females as observed during weekly surveys in July and August (Fig. S21). Annual
recruitment rates were then categorized as above or below 40% as this threshold typically indicates an increasing or declining population size (Reynolds, 1998; Schmidt et al., 2015).

We then used the multivariate ordination technique Principal Component Analysis (PCA) to explore whether RSF coefficients could be used to separate calf recruitment years above and below 40%. We also included spring snow cover (%) for each year as a variable in the PCA as increased snow cover is known to negatively influence muskox calf recruitment (Fig. S21; Schmidt et al., 2015). Muskox cows give birth to a calf in late winter if sufficient energetic reserves are acquired the preceding summer (Desforges et al., 2019). Therefore, the RSF coefficients acquired for year *t* were matched to the calf recruitment rate observed the year after (year *t*+1) in all analyses. Lastly, we fitted a separate PCA using the annual coefficients of the classical RSF, based on vegetation types, to assess if these could be used to separate calf recruitment years above and below 40% and to compare findings to the element based PCA results.

Because PCA is primarily used to explore patterns in multivariate data rather than to test formal hypotheses (Jolliffe, 2002), which is part of our study objective (Hypotheses 1 & 2), we used beta regression to statistically relate annual variation in the relative use of the different chemical elements to calf recruitment. Thus, muskox calf recruitment in year *t*+1 was the response variable, while the relative use of the various elements in year *t* (i.e., annual RSF coefficients) as well as spring snow cover were explanatory variables. Again, separate beta models were fitted for each element due to multicollinearity.

All statistical analyses for this study were conducted using the software package R, version 4.2.3 (R Development Core Team, 2023).

## 3. Results
316

### 3.1. Essential and non-essential elements in soil and vegetation

317

Measurements of the elements in soil and vegetation samples collected at the Zackenberg valley in northeast Greenland are presented in Fig. 2 and Table S2 and form the numerical basis for the construction of the geochemical landscapes (Fig. 3).

318

319

320

321

### 3.2 Relative use of geochemical landscapes and vegetation types

322

The relative use of the geochemical landscape by muskoxen varied over time and differed between elements (Fig. S22). For example, the element-based RSFs revealed that the relative use of vegetated sites decreased with increasing C during the years 2014-2017, a pattern that was also evident for C:N ratio and the trace elements Co and Hg (Fig. S22). The relative use of vegetated sites with high levels of N, Mn, Zn, and Ba was variable across years. In fact, for most elements and years, the 95% confidence intervals of the RSF coefficients overlapped with 0 (Fig. S22), which indicates that the relative use of the different essential and non-essential elements did not differ significantly from their availability in the landscape. Nonetheless, the predictive capacity of the element-based RSFs was satisfactory with a mean *r*<sub>s</sub> of 0.69 (sd=0.11) across all models, which was consistently higher (p<0.001) compared to the randomized model versions (Fig. S23).

323

324

325

326

327

328

329

330

331

332

333

Output of the classic vegetation-type based RSF also revealed considerable differences and annual variation in the relative probability of use of the dominant vegetation types (Fig. S24). The relative use of sites with *Cassiope* heath was higher during 2007-2014, but lower before and after this period. The opposite pattern was found for sites with *Salix* snowbed as the relative use was low during 2007-2014 but higher in other years, which is likely influenced by annual variation in the distribution of snow cover. The relative use of fen and grassland sites by muskoxen was stochastic over time as muskoxen used these vegetated sites significantly more in some years but used them significantly less in other years (Fig. S23).

334

335

336

337

338

339

340
341

S24). Only sites with *Dryas* heath were used proportional to their availability in the landscape in all years (i.e., no statistical difference between use and availability; Fig. S24). Predictive capacity of the classical RSF was also satisfactory with a mean *r*<sub>s</sub> of 0.71 (sd=0.16), which was significantly higher (p<0.001) than for a randomized model with a mean *r*<sub>s</sub> of 0.03 (sd=0.33).

342

343

344

345

346

347

### 3.3 Relating muskox calf recruitment to relative use of geochemical landscapes and vegetation types

348

349

PCA ordination on the element-specific RSF coefficients and snow cover resulted in a clear separation of high and low calf recruitment years, with the first two axes explaining 67.5% of the between year variation (Fig. 4). C:N contributed most (13.6%) to PCA axis 1, followed by Ba (13.3%), C (12.8%), and N (11.0%). Se and Mo contributed most to PCA axis 2 (20.9% and 15.2% respectively), followed by As (14.5%), and Cu (11.8%). The remaining variables contributed <10% to the first two PCA axes.

350

351

352

353

354

355

In comparison, PCA ordination on the coefficients derived from the vegetation type based RSF and snow cover resulted in a slightly poorer separation of high and low calf recruitment years, although the first two axes explained 78.1% of the between year variation (Fig. 5). Here, the relative use of *Cassiope* heath contributed most (31.4%) to PCA axis 1, followed by *Salix* snowbed (28.8%), *Dryas* heath (27.1%), and Fen areas (12.5%). Snow cover contributed most (40.7%) to PCA axis 2, followed by Grasslands (39.4%), and Fen areas (16.5%).

356

357

358

359

360

361

362

Output of the univariate beta regression (Fig. 6) revealed that muskox calf recruitment was positively correlated to the relative use of sites with increased N content (p=0.048) and concentrations of Cu (p=0.0003), Se (p=0.031), and Mo (p<0.0001). Muskox calf recruitment

363

364
365

was negatively correlated to the relative use of sites with increased concentrations of As (p=0.013) and Pb (p=0.045).

366

367

368

## 4. Discussion

369

By coupling the relative use of geochemical landscapes to long-term variation in annual calf recruitment of muskoxen, the largest herbivore roaming the Arctic tundra, this study provides empirical evidence that spatial variation of multiple essential and non-essential elements in vegetation and underlying soils can influence habitat-performance relationships.

370

371

372

373

374

### 4.1 Essential elements and muskox reproduction

375

We uncovered positive associations between annual calf recruitment of muskoxen and the relative use of vegetated sites with increased concentrations of the essential elements N, Cu, Se and Mo (supporting Hypothesis 1). It has long been established that deficiencies of the trace elements Cu and Se can impair reproductive success in both free-ranging and captive animals e.g., by reducing conception rates and increasing risks of embryonic and neonatal mortality (Flueck et al., 2012; Flynn et al., 1977; Robbins, 1983). The mean Se concentration across all plants sampled in our study was 0.04 µg g⁻¹ d.w. and below the deficiency risk threshold of 0.05 µg g⁻¹ d.w. proposed for free-ranging herbivores (Flueck et al., 2012).

376

377

378

379

380

381

382

383

Similarly, the mean Cu concentration across all plants sampled in our study was 7.73 µg g⁻¹ d.w. and below the deficiency risk threshold of 10 µg g⁻¹ d.w. proposed for deer (Grace and Wilson, 2002). It thus seems probable that poor muskox recruitment years were, at least partly, influenced by a shortage of Se and Cu in the landscape, as previously suggested for muskoxen in this region (Mosbacher et al., 2022). The positive association between Mo and reproductive output detected here aligns with experimental work on goats that showed a decline in conception and fetal survival rates reduced Mo intake (Anke et al., 2010).

384

385

386

387

388

389
However, the mean Mo concentration in vegetation samples collected in this study (0.98 µg g⁻¹d.w.) is within the range required for normal plant production (Kabata-Pendias and Szteke, 2015) suggesting that the bioavailability of Mo in the study landscape is not limited.

Moreover, the bioavailability of Mo is unlikely to exceed toxicity levels of muskoxen in this system, as has been reported for other herbivores (Reis et al., 2010).

Arctic and boreal ecosystems are characterized by strong N limitation and high C stocks in soils (Liu et al., 2018). As a consequence, most ecological stoichiometry studies focus on C:N ratios to assess vegetation quality and dynamics at northern latitudes (Mosbacher et al., 2019; Peng et al., 2011) or to estimate linkages with consumer distribution (Leroux et al., 2017; Rizzuto et al., 2021). In our study, however, we found no apparent association between C:N ratio in vegetated sites used by muskoxen and annual calf recruitment. Instead, muskoxen often avoided vegetated sites with high C content but selected for vegetated sites with increased content of N, with the latter having a strongly positive effect on annual calf recruitment. The importance of N in Arctic vegetation has previously been highlighted for caribou (*Rangifer tarandus*), and N limitation was suspected to have consequences for the distribution and growth of populations (Barboza et al., 2018). Our study corroborates this notion by revealing that N content in Arctic plants is a major driver of habitat use of muskoxen and that selection for vegetated sites with increased N content has a direct positive effect on recruitment rate.

Overall, we interpret these results as clear indications that the concentrations, distribution, and bioavailability of N and several essential trace elements in Arctic landscapes play critical roles in the reproductive success of muskoxen. Moreover, results of this study lend support to the notion that alternate year breeding in muskox populations in northeast Greenland is driven by depletion and slow recovery of elemental and nutritional reserves needed for a successful reproductive cycle (Desforges et al., 2019; Thing et al., 1987).
Finally, we suspect that constraints in the uptake of the essential elements N, Cu, Se, Mo in poor muskox recruitment years affects the reproductive cycle through reduced conception rates or increased embryonic and neonatal mortality as also observed in other herbivores (Anke et al., 2010; Flueck et al., 2012).

419

## 4.2 Non-essential elements and muskox reproduction

We detected a negative association between annual calf recruitment of muskoxen and the relative use of vegetated sites with higher concentrations of the non-essential elements As and Pb (supporting Hypothesis 2). At first glance, these findings suggest that the concentrations of these contaminants were sufficiently high to negatively affect the reproductive organs of muskoxen, as has been observed for other species occupying contaminated areas (Massányi et al., 2020). However, the average concentrations of Pb (0.5 µg g⁻¹ d.w.), As (0.12 µg g⁻¹ d.w.), and Hg (0.02 µg g⁻¹ d.w) measured in plant material collected in our study area were well below those reported for plants in areas considered to be (semi-) polluted with heavy metals (Kabata-Pendias and Szteke, 2015). Possible direct, degenerative effects of contaminants in the vegetation on the reproductive organs of muskoxen is therefore unlikely. Instead, these results may be a product of annual variation in the relative use of *Cassiope* and *Dryas* heath, i.e., the vegetation types with the highest concentrations of As and Pb in our study area (Figs. S12 & S18), and their negative correlations with muskox calf recruitment (Fig. 5). Muskoxen generally do not select for heath covered sites as forage plant productivity and biomass in those vegetation types is low compared to preferred and more productive fen and grassland areas (Tomassini et al., 2019). Thus, in years that muskoxen increase their use of heath vegetation, possibly because of constraints in forage accessibility elsewhere due to snow cover, recruitment rates decline and the relative use of sites with higher concentrations, and thus potential bioavailability, of As

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439
and Pb increases. It is, of course, possible that atmospheric deposition of contaminants will increase the concentrations and subsequent bioavailability of Cd, Pb, As or Hg over time, as seen elsewhere in the Arctic (Aslam et al., 2019; Couture et al., 2008; Dastoor et al., 2022; McConnell et al., 2019). Thus, continued monitoring of contaminant levels in Arctic vegetation and their potential direct and indirect effects on resident wildlife is pertinent.

445

446 *4.3 Caveats and future prospects*

To link long-term (1996-2021) variation in reproductive output of muskoxen to the use of geochemical landscapes, we relied on chemical analyses of plant and soil samples collected during the peak-biomass period (July) of a single year (2021). An important assumption underlying our work is therefore that the spatial distribution, concentrations, and bioavailability of the elements considered here remained constant over time. Yet, multiple processes can alter the spatiotemporal dynamics in vegetation composition and productivity (Elmendorf et al., 2012), which may impact the bioavailability of chemical elements. In fact, through foraging, defecating, and trampling, wildlife species themselves can alter the distribution and cycling of essential and non-essential elements across landscapes (Abraham et al., 2023; Schmitz et al., 2018). Further work is, therefore, recommended to assess if and how the concentrations, and thus potential bioavailability, of (trace) elements in plants and soils vary between years and throughout the growing season and to determine the impacts on resident wildlife. Long-term herbivore exclosures provide a suitable experimental set-up to address such outstanding questions.

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

Our study relied on relatively small-scale (52m²), yet long-term, survey data of muskox observations and demography. However, muskoxen typically roam over larger areas (Schmidt et al., 2016) with substantial individual variation in habitat use (Beumer et al., 2023; van Beest et al., 2020). To further investigate the importance of geochemical

462

463

464
465 landscapes in wildlife ecology, we see particular value in the application of high-resolution
466 biologging data (e.g., GPS collar and 3-dimensional accelerometers). Doing so will allow for
467 the estimation of spatial and temporal variation in feeding behavior of individual animals
468 across geochemical landscapes and multiple seasons (Nunes et al., 2022). Potential
469 consequences of individual variation in habitat use and behavior on Darwinian fitness (i.e.,
470 survival and lifetime reproductive success) can then be assessed via long-term monitoring
471 data or through Dynamic Energy Budget models (Desforges et al., 2019; Sperfeld et al.,
472 2017).

473 Free-ranging animals are confronted by a myriad of environmental conditions that
474 they need to balance to optimize fitness (Gaillard et al., 2010; Stein et al., 2014; van Beest
475 and Milner, 2013). In our analyses, we focused on vegetation types, snow cover, and multiple
476 essential and non-essential elements but discounted other potential factors (e.g., ambient
477 temperature, density dependence, predation risk) influencing habitat use. We argue that this
478 analytical approach is reasonable given that habitat use of muskoxen in this low-predator
479 system is primarily influenced by snow cover, plant productivity and access to forage
480 biomass during spring and summer (Beumer et al., 2020; Tomassini et al., 2019). Indeed,
481 muskoxen use this period primarily to replenish energetic reserves to survive the upcoming
482 winter and to be able to develop a fetus (Desforges et al., 2021). This is also why we ran
483 separate RSF and PCA analyses for the chemical elements (Fig. 4) and the vegetation types
484 (Fig. 5) to facilitate comparison. The ‘classic’ PCA and RSF based on vegetation types
485 explained more variation in annual recruitment (78.1%) than the element-based analyses
486 (67.5%), which supports the idea that vegetation types are a better predictor of muskox
487 habitat use and performance. Nonetheless, the element-based RSF and PCA also showed high
488 predictive power and, moreover, provide insight into the geochemical mechanism underlying
489 linkages between habitat use and population performance.
The Arctic is now warming four times faster than the global average (Rantanen et al., 2022). Rapid thawing of permafrost threatens to release biological, chemical, and radioactive materials that have been sequestered for tens to hundreds of thousands of years. As these stored compounds re-enter the environment, they have the potential to disrupt ecosystem function endangering the health of humans and Arctic wildlife (Miner et al., 2021). While the potential risks of leaching of trace metals following permafrost degradation have been studied to some extent (Antcibor et al., 2014), little is known about possible changes in composition and spatial distribution of essential trace elements following permafrost thaw. We foresee great value in the application of ecogeochemical methods in combination with climate-forecast models to assess shifts in the bioavailability, concentration, and distribution of essential and non-essential elements over time and space. This will provide the necessary foundation to assess habitat-performance relationships under contemporary and future environmental change.

503

## 5. Conclusion

Our study clearly shows that the distribution and bioavailability of the essential elements N, Cu, Se, and Mo in Arctic vegetation are critical to the reproductive success of muskoxen and, moreover, possible contributing factors to alternate year breeding as observed in populations across northeast Greenland. The distribution and bioavailability of the non-essential elements As and Pb negatively correlated with calf recruitment of muskoxen, but the concentrations measured were unlikely to cause degenerative impacts on reproductive organs. Instead, this result likely reflects a constraint on forage accessibility in poor recruitment years due to extensive snow cover, leading to an increased use of less optimal habitat with higher concentrations of As and Pb and lower concentrations of the elements essential for reproduction.
The rapid warming of the Arctic and subsequent thawing of permafrost is suspected to severely alter the composition and spatial distribution of essential and non-essential elements. Therefore, continued monitoring and measurements of chemical elements in Arctic soils and vegetation as well as responses by resident wildlife are critical to assess habitat-performance relationships under contemporary and future environmental change. While our study is specific to high-Arctic muskoxen, we posit that the influence of geochemical landscapes on animal performance highlighted here also affect the ecology of wildlife species living in lower latitude systems.

523

## Funding

This study was financially supported by the Danish Environmental Protection Agency (J.nr. 2019-7208), by ANR JCJC (ANR-21-CE34-0001 ATCAF-project), the Carlsberg foundation (CF17-0097 and CF18-0279) and Independent Research Fund Denmark/ DFF Nature and Universe, grant No. 0135-00140B.

529

## Ethics approval

Not applicable.

## Data availability

Data on concentrations and percentages of the essential and non-essential elements are provided in Table S2 in the Supporting Information.

535

## Acknowledgements

We thank Greenland Ecosystem Monitoring for access to ecosystem data, and Aarhus University, Denmark, and Zackenberg Research Station, Greenland, for logistical support at the study site. We also thank Laetitia Leroy for assistance in the clean room, and Camille
540 Duquenoy and Aurélie Marquet (both at the ICP-MS platform at the Observatoire Midi-

541 Pyrenees, Toulouse, France) for assistance with the ICP-MS analyses. We also thank three

542 anonymous reviewers for their extremely constructive comments that helped to improve this
543 article.

544

545 **References**

546 Abraham, A.J., Duvall, E., Ferraro, K., Webster, A.B., Doughty, C.E., le Roux, E., Ellis-Soto,  
547 D., 2023. Understanding anthropogenic impacts on zoogeochemistry is essential for  
548 ecological restoration. Restor. Ecol. 31, e13778. https://doi.org/10.1111/rec.13778

549 Anke, M., Seifert, M., Arnhold, W., Anke, S., Schäfer, U., 2010. The biological and  
550 toxicological importance of molybdenum in the environment and in the nutrition of  
551 plants, animals and man: Part V: Essentiality and toxicity of molybdenum. Acta  
552 Aliment. 39, 12–26. https://doi.org/10.1556/aalim.39.2010.1.2

553 Antcibor, I., Eschenbach, A., Zubrzycki, S., Kutzbach, L., Bolshiyanov, D., Pfeiffer, E.-M.,  
554 2014. Trace metal distribution in pristine permafrost-affected soils of the Lena River  
555 delta and its hinterland, northern Siberia, Russia. Biogeosciences 11, 1–15.
556 https://doi.org/10.5194/bg-11-1-2014

557 Balluffi-Fry, J., Leroux, S.J., Wiersma, Y.F., Heckford, T.R., Rizzuto, M., Richmond, I.C.,  
558 Vander Wal, E., 2020. Quantity–quality trade-offs revealed using a multiscale test of  
559 herbivore resource selection on elemental landscapes. Ecol Evol 10, 13847–13859.
560 https://doi.org/10.1002/ECE3.6975

561 Barboza, P.S., Van Someren, L.L., Gustine, D.D., Bret-Harte, M.S., 2018. The nitrogen  
562 window for arctic herbivores: plant phenology and protein gain of migratory caribou  
563 (Rangifer tarandus). Ecosphere 9. https://doi.org/10.1002/ecs2.2073

564 Bastille-Rousseau, G., Wittemyer, G., 2019. Leveraging multidimensional heterogeneity in resource selection to define movement tactics of animals. Ecol. Lett. 22, 1417–1427. https://doi.org/10.1111/ele.13327

567 Beumer, L.T., Pohle, J., Schmidt, N.M., Chimienti, M., Desforges, J.-P., Hansen, L.H.,  
568 Langrock, R., Pedersen, S.H., Stelvig, M., van Beest, F.M., 2020. An application of  
569 upscaled optimal foraging theory using hidden Markov modelling: year-round  
570 behavioural variation in a large arctic herbivore. Mov. Ecol. 8, 25.
571 https://doi.org/10.1186/s40462-020-00213-x

572 Beumer, L.T., Schmidt, N.M., Pohle, J., Signer, J., Chimienti, M., Desforges, J.-P., Hansen, L.H., Højlund Pedersen, S., Rudd, D.A., Stelvig, M., van Beest, F.M., 2023. Accounting for behaviour in fine-scale habitat selection: A case study highlighting methodological intricacies. J Anim Ecol Early View. https://doi.org/10.1111/1365-2656.13984

577 Bhattacharya, P.T., Misra, S.R., Hussain, M., 2016. Nutritional aspects of essential trace elements in oral health and disease: an extensive review. Scientifica (Cairo) 2016, 5464373. https://doi.org/10.1155/2016/5464373

580 Boersma, M., Elser, J.J., 2006. Too Much of a Good Thing: On Stoichiometrically Balanced Diets and Maximal Growth. Ecology 87, 1325–1330. https://doi.org/10.1890/0012-9658(2006)87[1325:TMOAGT]2.0.CO;2
583 Boyce, M.S., McDonald, L.L., 1999. Relating populations to habitats using resource selection functions. *Trends Ecol. Evol.* 14, 268–272.

584 Boyce, M.S., Vernier, P.R., Nielsen, S.E., Schmiegelow, F.K.A., 2002. Evaluating resource selection functions. *Ecol Model* 157, 281–300.

587 Brown, G.E., Foster, A.L., Ostergren, J.D., 1999. Mineral surfaces and bioavailability of heavy metals: A molecular-scale perspective. *Proc. Natl. Acad. Sci. U.S.A.* 96, 3388–3395. https://doi.org/10.1073/pnas.96.7.3388

590 Burger, J., 2008. Assessment and management of risk to wildlife from cadmium. *Sci. Total Environ.* 389, 37–45. https://doi.org/10.1016/j.scitotenv.2007.08.037

592 Cable, S., Christiansen, H.H., Westergaard-Nielsen, A., Kroon, A., Elberling, B., 2018. Geomorphological and cryostratigraphical analyses of the Zackenberg Valley, NE Greenland and significance of Holocene alluvial fans. *Geomorphology* 303, 504–523. https://doi.org/10.1016/J.GEOMORPH.2017.11.003

596 Desforges, J.-P., Gonçalo, J., Marques, M., Beumer, L.T., Chimienti, M., Hansen, L.H., Stine, J., Pedersen, H., Schmidt, N.M., Van Beest, F.M., 2021. Environment and physiology shape Arctic ungulate population dynamics. *Glob Change Biol* 27, 1755–1771. https://doi.org/10.1111/gcb.15484

600 Desforges, J.P., Marques, G.M., Beumer, L.T., Chimienti, M., Blake, J., Rowell, J.E., Adamczewski, J., Schmidt, N.M., van Beest, F.M., 2019. Quantification of the full lifecycle bioenergetics of a large mammal in the high Arctic. *Ecol Model* 401, 27–39. https://doi.org/10.1016/J.ECOLMODEL.2019.03.013

604 Durso, A.M., Smith, G.D., Hudson, S.B., French, S.S., 2020. Stoichiometric and stable isotope ratios of wild lizards in an urban landscape vary with reproduction, physiology, space and time. *Conserv Physiol* 8, coaa001. https://doi.org/10.1093/conphys/coaa001

608 Elberling, B., Tamstorf, M.P., Michelsen, A., Arndal, M.F., Sigsgaard, C., Illeris, L., Bay, C., Hansen, B.U., Christensen, T.R., Hansen, E.S., Jakobsen, B.H., Beyens, L., 2008. Soil and Plant Community-Characteristics and Dynamics at Zackenberg. *Adv. Ecol. Res.* 40, 223–248. https://doi.org/10.1016/S0065-2504(07)00010-4

612 Elmendorf, S.C., Henry, G.H.R., Hollister, R.D., Björk, R.G., Bjorkman, A.D., Callaghan, T.V., Collier, L.S., Cooper, E.J., Cornelissen, J.H.C., Day, T.A., Fosaa, A.M., Gould, W.A., Grétarsdóttir, J., Harte, J., Hermanutz, L., Hik, D.S., Hofgaard, A., Jarrad, F., Jónsdóttir, I.S., Keuper, F., Klanderud, K., Klein, J.A., Koh, S., Kudo, G., Lang, S.I., Loewen, V., May, J.L., Mercado, J., Michelsen, A., Molau, U., Myers-Smith, I.H., Oberbauer, S.F., Pieper, S., Post, E., Rixen, C., Robinson, C.H., Schmidt, N.M., Shaver, G.R., Stenström, A., Tolvanen, A., Totland, Ø., Troxler, T., Wahren, C.H., Webber, P.J., Welker, J.M., Wookey, P.A., 2012. Global assessment of experimental climate warming on tundra vegetation: Heterogeneity over space and time. *Ecol Lett* 15, 164–175. https://doi.org/10.1111/j.1461-0248.2011.01716.x

622 Elser, J., 2006. Biological stoichiometry: a chemical bridge between ecosystem ecology and evolutionary biology. *Am. Nat.* 168 Suppl, S25-35. https://doi.org/10.1086/509048

624 Feinberg, A., Stenke, A., Peter, T., Hinckley, E.-L.S., Driscoll, C.T., Winkel, L.H.E., 2021. Reductions in the deposition of sulfur and selenium to agricultural soils pose risk of future nutrient deficiencies. *Commun Earth Environ* 2, 1–8. https://doi.org/10.1038/s43247-021-00172-0

628 Felton, A.M., Holmström, E., Malmsten, J., Felton, A., Cromsigt, J.P.G.M., Edenius, L., Ericsson, G., Widemo, F., Wam, H.K., 2020. Varied diets, including broadleaved forage, are important for a large herbivore species inhabiting highly modified landscapes. *Sci Rep* 10, 1904. https://doi.org/10.1038/s41598-020-58673-5
632 Flueck, W.T., 1994. Effect of trace elements on population dynamics: Selenium deficiency in free-ranging black-tailed deer. Ecology 75, 807-812. https://doi.org/10.2307/1941736

634 Flueck, W.T., Smith-Flueck, J.M., Mionczynski, J., Mincher, B.J., 2012. The implications of selenium deficiency for wild herbivore conservation: a review. Eur J Wildl Res 58, 761-780. https://doi.org/10.1007/S10344-012-0645-Z

637 Flynn, A., Franzmann, A.W., Arneson, P.D., Oldemeyer, J.L., 1977. Indications of Copper Deficiency in a Subpopulation of Alaskan Moose. J Nutr 107, 1182-1189.
https://doi.org/10.1093/jn/107.7.1182

640 French, A.S., Shaw, D., Gibb, S.W., Taggart, M.A., 2017. Geochemical landscapes as drivers of trace and toxic element profiles in wild red deer (*Cervus elaphus*). Sci. Total Environ. 601-602, 1606-1618. https://doi.org/10.1016/j.scitotenv.2017.05.210

643 Gaillard, J.-M., Hebblewhite, M., Loison, A., Fuller, M., Powell, R., Basille, M., van Moorter, B., 2010. Habitat-performance relationships: finding the right metric at a given spatial scale. Philos Trans R Soc Lond B Biol Sci 365, 2255-2265.
https://doi.org/10.1098/rstb.2010.0085

647 Grace, N.D., Wilson, P.R., 2002. Trace element metabolism, dietary requirements, diagnosis and prevention of deficiencies in deer. New Zealand Veterinary Journal 50, 252-259.
https://doi.org/10.1080/00480169.2002.36321

650 Jolliffe, I.T. (Ed.), 2002. Principal Component Analysis for Special Types of Data, in: Principal Component Analysis, Springer Series in Statistics. Springer, New York, NY, pp. 338-372. https://doi.org/10.1007/0-387-22440-8_13

653 Jutha, N., Jardine, C., Schwantje, H., Mosbacher, J., Kinniburgh, D., Kutz, S., 2022. Evaluating the use of hair as a non-invasive indicator of trace mineral status in woodland caribou (Rangifer tarandus caribou). PLOS ONE 17, e0269441.
https://doi.org/10.1371/JOURNAL.PONE.0269441

657 Kabata-Pendias, A., Szteke, B., 2015. Trace Elements in Abiotic and Biotic Environments. CRC Press, Boca Raton. https://doi.org/10.1201/b18198

659 Kaspari, M., Powers, J.S., 2016. Biogeochemistry and geographical ecology: Embracing all twenty-five elements required to build organisms. Am. Nat. 188, S62-S73.
https://doi.org/10.1086/687576/ASSET/IMAGES/LARGE/FG4.JPEG

662 Leroux, S.J., Wal, E.V., Wiersma, Y.F., Charron, L., Ebel, J.D., Ellis, N.M., Hart, C., Kissler, E., Saunders, P.W., Moudrá, L., Tanner, A.L., Yalcin, S., 2017. Stoichiometric distribution models: ecological stoichiometry at the landscape extent. Ecol Lett 20, 1495-1506. https://doi.org/10.1111/ELE.12859

666 Liu, X.-Y., Koba, K., Koyama, L.A., Hobbie, S.E., Weiss, M.S., Inagaki, Y., Shaver, G.R., Giblin, A.E., Hobara, S., Nadelhoffer, K.J., Sommerkorn, M., Rastetter, E.B., Kling, G.W., Laundre, J.A., Yano, Y., Makabe, A., Yano, M., Liu, C.-Q., 2018. Nitrate is an important nitrogen source for Arctic tundra plants. Proc. Natl. Acad. Sci. U.S.A. 115, 3398-3403. https://doi.org/10.1073/pnas.1715382115

671 Manly, B.F.J., McDonald, L.L., Thomas, D.L., McDonald, T.L., Erickson, W.P., 2002. Resource selection by animals: Statistical analysis and design for field studies. Kluwer Academic Publishers, Dordrecht, The Netherlands.

674 Massányi, P., Massányi, M., Madeddu, R., Stawarz, R., Lukáč, N., 2020. Effects of Cadmium, Lead, and Mercury on the Structure and Function of Reproductive Organs. Toxics 8, 94. https://doi.org/10.3390/toxics8040094

677 Matthiopoulos, J., Fieberg, J., Aarts, G., Beyer, H.L., Morales, J.M., Haydon, D.T., 2015. Establishing the link between habitat selection and animal population dynamics. Ecol. Monogr. 85, 413-436. https://doi.org/10.1890/14-2244.1

680 McLoughlin, P.D., Gaillard, J.M., Boyce, M.S., Bonenfant, C., Messier, F., Duncan, P., Delorme, D., Van Moorter, B., Saïd, S., Klein, F., 2007. Lifetime reproductive
682
683

success and composition of the home range in a large herbivore. Ecology 88, 3192–3201.

684
685

Michaluk, A., Kochman, K., 2007. Involvement of copper in female reproduction. Reprod Biol 7, 193–205.

686
687
688
689

Miner, K.R., D'Andrilli, J., Mackelprang, R., Edwards, A., Malaska, M.J., Waldrop, M.P., Miller, C.E., 2021. Emergent biogeochemical risks from Arctic permafrost degradation. Nat. Clim. Chang. 11, 809–819. https://doi.org/10.1038/s41558-021-01162-y

690
691
692
693
694

Mosbacher, J.B., Desforges, J.-P., Michelsen, A., Hansson, S.V., Stelvig, M., Eulaers, I., Sonne, C., Dietz, R., Jenssen, B.M., Ciesielski, T.M., Lierhagen, S., Flaten, T.P., Roux, G.L., Aggerbeck, M.R., Schmidt, N.M., 2022. Hair mineral levels as indicator of wildlife demographics?—a pilot study of muskoxen. Polar Res. 41. https://doi.org/10.3326/polar.v41.8543

695
696
697
698
699
700

Mosbacher, J.B., Michelsen, A., Stelvig, M., Hjermstad-Sollerud, H., Schmidt, N.M., 2019. Muskoxen modify plant abundance, phenology, and nitrogen dynamics in a high Arctic fen. Ecosystems 22, 1095–1107. https://doi.org/10.1007/s10021-018-0323-4

Mosser, A., Fryxell, J.M., Eberly, L., Packer, C., 2009. Serengeti real estate: density vs. fitness-based indicators of lion habitat quality. Ecol. Lett. 12, 1050–1060. https://doi.org/10.1111/j.1461-0248.2009.01359.x

701
702
703

Northrup, J.M., Hooten, M.B., Anderson, C.R., Wittemyer, G., 2013. Practical guidance on characterizing availability in resource selection functions under a use-availability design. Ecology 94, 1456–1463. https://doi.org/10.1890/12-1688.1

704
705
706
707
708
709
710

Nunes, G.T., Efe, M.A., Barreto, C.T., Gaiotto, J.V., Silva, A.B., Vilela, F., Roy, A., Bertrand, S., Costa, P.G., Bianchini, A., Bugoni, L., 2022. Ecological trap for seabirds due to the contamination caused by the Fundão dam collapse, Brazil. Sci. Total Environ. 807, 151486. https://doi.org/10.1016/j.scitotenv.2021.151486

Ofstad, E.G., Markussen, S.S., Saether, B.-E., Solberg, E.J., Heim, M., Haanes, H., Røed, K.H., Herfindal, I., 2020. Opposing fitness consequences of habitat use in a harvested moose population. J Anim Ecol. https://doi.org/10.1111/1365-2656.13221

711
712
713
714
715

Pedersen, S.H., Tamstorf, M.P., Abermann, J., Westergaard-Nielsen, A., Lund, M., Skov, K., Sigsgaard, C., Mylius, M.R., Hansen, B.U., Liston, G.E., Schmidt, N.M., 2016. Spatiotemporal characteristics of seasonal snow cover in northeast Greenland from in situ observations. Arct. Antarct. Alp. Res 48, 653–671. https://doi.org/10.1657/AAAR0016-028

716
717
718

Peng, Y., Niklas, K.J., Sun, S., 2011. The relationship between relative growth rate and whole-plant C : N : P stoichiometry in plant seedlings grown under nutrient-enriched conditions. J Plant Ecol 4, 147–156. https://doi.org/10.1093/jpe/rtq026

719
720
721

Plath, K., Boersma, M., 2001. Mineral Limitation of Zooplankton: Stoichiometric Constraints and Optimal Foraging. Ecology 82, 1260–1269. https://doi.org/10.1890/0012-9658(2001)082[1260:MLOZSC]2.0.CO;2

722
723
724

Proffitt, K.M., Hebblewhite, M., Peters, W., Hupp, N., Shamhart, J., 2016. Linking landscape-scale differences in forage to ungulate nutritional ecology. Ecol. Appl. 26, 2156–2174.

725
726
727
728
729
730
731

R Development Core Team, 2023. R: A language and environment for statistical computing.

Rantanen, M., Karpechko, A.Y., Lipponen, A., Nordling, K., Hyvärinen, O., Ruosteenoja, K., Vihma, T., Laaksonen, A., 2022. The Arctic has warmed nearly four times faster than the globe since 1979. Commun Earth Environ 3, 1–10. https://doi.org/10.1038/s43247-022-00498-3

Reis, L.S.L.S., Pardo, P.E., Camargos, A.S., Oba, E., 2010. Mineral element and heavy metal poisoning in animals. J Medicine & Med. Sci. 1, 560–579.
732 Reynolds, P.E., 1998. Dynamics and range expansion of a re-established muskox population.
733 J Wildl Manage 62, 734-744. https://doi.org/10.2307/3802350
734 Rioux, È., Pelletier, F., Mosbacher, J.B., Lesmerises, F., St-Louis, R., Kutz, S., St-Laurent,
735 M.-H., 2022. Links between individual performance, trace elements and stable
736 isotopes in an endangered caribou population. Glob. Ecol. Conserv. 38, e02234.
737 https://doi.org/10.1016/j.gecco.2022.e02234
738 Rizzuto, M., Leroux, S.J., Wal, E.V., Richmond, I.C., Heckford, T.R., Balluffi-Fry, J.,
739 Wiersma, Y.F., 2021. Forage stoichiometry predicts the home range size of a small
740 terrestrial herbivore. Oecologia 1, 327-338. https://doi.org/10.1007/s00442-021-
741 04965-0
742 Robbins, C.T., 1983. Wildlife feeding and nutrition. Academic Press Inc, New York, USA.
743 Rombach, E.P., Barboza, P.S., Blake, J.E., 2002. Utilization of copper during lactation and
744 neonatal development in muskoxen. Can. J. Zool. 80, 1460-1469.
745 https://doi.org/10.1139/Z02-139
746 Roos-Barraclough, F., Givelet, N., Martinez-Cortizas, A., Goodsite, M.E., Biester, H.,
747 Shotyk, W., 2002. An analytical protocol for the determination of total mercury
748 concentrations in solid peat samples. Sci Total Environ 292, 129-139.
749 https://doi.org/10.1016/s0048-9697(02)00035-9
750 Sach, F., Yon, L., Henley, M.D., Bedetti, A., Buss, P., de Boer, W.F., Dierenfeld, E.S.,
751 Gardner, A., Langley-Evans, S.C., Hamilton, E., Lark, R.M., Prins, H.H.T.,
752 Swemmer, A.M., Watts, M.J., 2020. Spatial geochemistry influences the home range
753 of elephants. Sci Total Environ 729, 139066.
754 https://doi.org/10.1016/J.SCITOTENV.2020.139066
755 Schmidt, N.M., Pedersen, S.H., Mosbacher, J.B., Hansen, L.H., 2015. Long-term patterns of
756 muskox (Ovibos moschatus) demographics in high arctic Greenland. Polar Biol. 38,
757 1667-1675. https://doi.org/10.1007/s00300-015-1733-9
758 Schmidt, N.M., van Beest, F.M., Mosbacher, J.B., Stelvig, M., Hansen, L.H., Nabe-Nielsen,
759 J., Grøndahl, C., 2016. Ungulate movement in an extreme seasonal environment:
760 year-round movement patterns of high-arctic muskoxen. Wildlife Biol 22, 253-267.
761 https://doi.org/10.2981/wlb.00219
762 Schmitz, O.J., Wilmers, C.C., Leroux, S.J., Doughty, C.E., Atwood, T.B., Galetti, M.,
763 Davies, A.B., Goetz, S.J., 2018. Animals and the zoogeochemistry of the carbon
764 cycle. Science 362, eaar3213. https://doi.org/10.1126/science.aar3213
765 Sperfeld, E., Wagner, N.D., Halvorson, H.M., Malishev, M., Raubenheimer, D., 2017.
766 Bridging Ecological Stoichiometry and Nutritional Geometry with homeostasis
767 concepts and integrative models of organism nutrition. Funct. Ecol. 31, 286-296.
768 https://doi.org/10.1111/1365-2435.12707
769 Squadrone, S., Robetto, S., Orusa, R., Griglione, A., Falsetti, S., Paola, B., Abete, M.C.,
770 2022. Wildlife hair as bioindicators of metal exposure. Biol. Trace Elem. Res.
771 https://doi.org/10.1007/s12011-021-03074-6
772 Stein, A., Gerstner, K., Kreft, H., 2014. Environmental heterogeneity as a universal driver of
773 species richness across taxa, biomes and spatial scales. Ecol. Lett. 17, 866-880.
774 https://doi.org/10.1111/ele.12277
775 Stewart, L., Simonsen, C.E., Svenning, J.C., Schmidt, N.M., Pellissier, L., 2018. Forecasted
776 homogenization of high Arctic vegetation communities under climate change. J.
777 Biogeogr. 45, 2576-2587. https://doi.org/10.1111/JBI.13434
778 Thing, H., Klein, D.R., Jingfors, K., Holt, S., 1987. Ecology of Muskoxen in Jameson Land,
779 Northeast Greenland. Holarctic Ecol. 10, 95-103.
780 Tomassini, O., van Beest, F.M., Schmidt, N.M., 2019. Density, snow and seasonality lead to

781 variation in muskox habitat selection during summer. Can. J. Zool. 97, 997–1003.

782 https://doi.org/10.1139/cjz-2018-0292

783 Trueman, C.N., Chung, M.-T., Shores, D., 2016. Ecogeochemistry potential in deep time

784 biodiversity illustrated using a modern deep-water case study. Philos Trans R Soc

785 Lond B Biol Sci 371, 20150223. https://doi.org/10.1098/rstb.2015.0223

786 van Beest, F.M., Beumer, L.T., Chimienti, M., Desforges, J.-P., Huffeldt, N.P., Pedersen,

787 S.H., Schmidt, N.M., 2020. Environmental conditions alter behavioural organization

788 and rhythmicity of a large Arctic ruminant across the annual cycle. R. Soc. Open Sci.

789 7, 201614. https://doi.org/10.1098/rsos.201614

790 van Beest, F.M.M., Milner, J.M.M., 2013. Behavioural responses to thermal conditions affect

791 seasonal mass change in a heat-sensitive northern ungulate. PLOS ONE 8, e65972.

792 https://doi.org/10.1371/journal.pone.0065972

793 Verniest, F., Galewski, T., Julliard, R., Guelmami, A., Le Viol, I., 2022. Coupling future

794 climate and land-use projections reveals where to strengthen the protection of

795 Mediterranean Key Biodiversity Areas. Conserv. sci. pract 4, e12807.

796 https://doi.org/10.1111/csp2.12807

797
<figure><img src="image_2.png" /></figure>

800 Graphical Abstract
A)

<figure><img src="image_3.png" /></figure>

B)

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

C)

<figure><img src="image_6.png" /></figure>

803

<caption>Fig. 1: A) Map of Greenland with Zackenberg valley indicated by the black square. B)</caption>

804

Elevation map of the Zackenberg valley with the locations of the 50 field plots used in this study to sample plant and soil material. C) Elevation map of the Zackenberg valley with the locations of muskox observations (n=3142) during summer (Jun-Sep) for the years 1996–2021.

805

806

807

808

809
<figure><img src="image_7.png" /><figcaption>**Fig. 2:** Boxplots showing measurements of chemical elements in soil (yellow) and vegetation (green) samples collected during July 2021 in the Zackenberg valley, northeast Greenland.

The triangle plotted within the box shows the mean value. Note that elements are arranged based on their position in the periodic table (atomic number).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 3:** Predicted spatial distribution of chemical elements in vegetation across the Zackenberg valley in northeast Greenland. The unit for C, N and C:N are in percentages while the remaining elements are in µg g⁻¹ d.w.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 4:** Biplot showing how spring snow cover and the relative use of chemical elements in vegetated sites during summer (vectors) relate to muskox calf recruitment years in the Zackenberg valley, northeast Greenland.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 5:** Biplot showing how spring snow cover and the relative use of vegetation types during summer (vectors) relate to muskox calf recruitment years in the Zackenberg valley, northeast Greenland.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 6:** Effects of the relative probability of use of chemical elements in vegetated sites during summer (in year *t*) and calf recruitment (in year *t* +1) of muskoxen in the Zackenberg valley, northeast Greenland (1996-2021). Panels with solid black lines (predicted mean) and shaded areas (95% confidence intervals) indicate statistically significant relationships (p<0.05), while dashed lines indicate no statistical relationships.</figcaption></figure>