Geochemical landscapes as drivers of wildlife reproductive success Insights from a high-Arctic ecosystem

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

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  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

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 heterogenous 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).

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 1200 m 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 noncalcareous 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 600 m 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 weekly every summer (Jun-Sep) in Zackenberg valley with densities ranging between 0.5 and 3.0 individuals km − 2 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.

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-600 m) 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 by mimicking muskox grazing i.e., grabbing and pulling a handful of plant material.Plant material (min.~2 g wet weight) was subsequently transferred into paper bags.Soil was collected in each plot by mixing and homogenizing three sub-samples from the Ahorizon (i.e., the upper 0-10 cm layer) using nitrile gloves and a pre-cleaned stainless-steel garden shovel.Soil material (~10 g wet weight) was transferred into polyethylene zip-lock bags.Fourteen of the 50 field plots were sampled more intensively by collecting three soil samples and three plant samples from within 1 m 2 to assess and incorporate plotbased variability in element concentrations.In total, 78 soil samples and 78 plant samples were collected from within the field plots.All soil and plant material were stored frozen until sample processing and chemical analyses.

Sample processing and chemical analyses
Soil samples were dried in an oven at low temperature (35 • C) for 5 days under constant air flow to limit the loss of volatile elements (Roos-Barraclough et al., 2002).While we acknowledge that oven drying is not the preferred method to process biological samples for trace element analyses, doing so was necessary due to constraints in transporting samples from the field to the laboratory.Once dried and transported, soil  samples were sieved through a 200 μm nylon mesh to separate soil particles from small gravel and organic matter such as root fragments.Compact samples were manually ground beforehand with an agate mortar.Between each sample, all sieving material (mesh, container, mortar) was rinsed twice, first with regular tap water, then with mQ H 2 O (18.2 MΩ cm − 1 ) and finally disinfected with 70 % ethanol.Vegetation samples were lyophilized using a Christ Alpha 1-2 freeze-drier and ground using a laboratory mill Retsch® Mixer MILL MM 400, at 30 Hz for 120 s.All material was cleaned between each sample, i.e., rinsed under tap water then rinsed with mQ H 2 O and finally disinfected with 70 % ethanol.
For wet chemistry digestion of both vegetation and soil samples, approximately 150 mg of dried and ground sample was placed in precleaned (2 % HNO 3 and mQ H 2 O) 50 mL digitubes to which 3 mL of ultrapure HNO 3 (Fisher Chemical, Nitric Acid HNO 3 , Optima Grade, 67-69 %) was added.The tubes were then sealed with airtight caps and placed on a hot plate at 90 • C overnight (~12 h).After allowing the digest to cool to room temperature, the digested solutions were transferred to 50 mL Falcon tubes and diluted with mQ H 2 O to reach a mother solution (MS) of 45 mL.Prior to ICP-MS analysis, the MS was further diluted by transferring a 5 mL aliquot to a 15 mL falcon tube to which mQ H 2 O was 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 >10× 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 − 1 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-Pyrenees, Toulouse, France.
Following the trace element analyses, soil and vegetation samples were also measured for C and N concentrations (%), which were 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 Isoprime isotope ratio mass spectrometer (Isoprime Ltd., Cheadle Hulme, Stockport, UK) coupled to a CN elemental analyser (Eurovector, Milan, Italy) with continuous flow, with a front furnace operating at 1050 • C and a reduction furnace operating at 700 • C. For quality assurance, certified reference material (CRM NIST 1547 Peach leaves material, National Institute of Standards & Technology) was included with every 12 samples and processed the same way as the vegetation and soil samples.Certified reference material had a mean value ± standard deviation of 2.965 ± 0.054 for %N and 46.01 ± 1.52 for %C and were measured at 2.939 ± 0.080 for %N and at 46.67 ± 1.89 for %C in this study.The detection limits of the elements were 0.015 % for N and 0.100 % for C.

Quantifying geochemical landscapes
Plot-based measurements of %C, % N and C:N ratio as well as of trace element concentrations (μg g − 1 d.w.) (Fig. 2) were used to generate maps of their predicted distribution and potential bioavailability across the soils and vegetation types of the Zackenberg valley (i.e., geochemical landscapes).To do so, we used stoichiometric distribution models (StDMs) as proposed by Leroux et al. (2017).StDMs generate spatially explicit predictions of elemental patterns across the landscape based on statistical relationship between element concentrations/ratio and environmental parameters (e.g., elevation, vegetation cover, soil moisture, etc.).Such maps can then be used to link consumers' distribution or habitat use to element concentrations across natural landscapes (e.g., Sach et al., 2020) as described in Section 2.5.
StDMs were computed using generalized linear models (GLMs) and constructed for each element (n = 14) and for both soil and vegetation samples separately (i.e., a total of 28 models).To improve normality of the data, concentrations of trace elements were log-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 1 m (Fig. S1) including soil moisture (index), spring snow cover (%), ambient temperature ( • C) solar radiation (kWh m − 2 ), 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 in the soil.For an overview of how environmental variables contribute to the clustering and spatial distribution of element concentrations 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 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 mean squared error [MSE] (using the testing data) and the amount of variation explained R 2 (using the training data).We then compared the distribution of MSE and R 2 of the testing and training data to those generated by 100 models where the response variable (i.e., elemental concentration or ratio) of the original data was randomly resampled.We considered the performance of the StDMs to be satisfactory as the prediction error (MSE) was always lower and the amount of variation explained (R 2 ) was always higher than for the models based on randomized datasets (see Figs. S3 and S4).The output of the elementspecific 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.

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 (i.e., selection) 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 during weekly surveys every 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 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 (Spearman rank correlations (r s ) of >0.6 and <− 0.6 for some element combinations, Fig. S20), separate RSFs were fitted for each element.Values of the explanatory variables were scaled prior to analyses to ensure RSF coefficients across years and elements were comparable.In all RSFs, we also fitted Year as a random intercept to account for unbalanced data.
Because no evidence exists that muskoxen, or other herbivorous animals, can sense the elemental composition in plant material, they are unlikely to actively select forage resources based on the bioavailability of different chemical elements.Instead, plant biomass and composition and thus vegetation type is likely a more realistic driver of muskox habitat use in this area as plant productivity and biomass varies between vegetation classes (Tomassini et al., 2019).Therefore, and for comparative purposes, we fitted a classic RSF using the interaction between vegetation type and Year as the independent variable.
Finally, the predictive success of all RSFs was evaluated using the kfold cross-validation procedure as proposed by Boyce et al. (2002) and Bastille-Rousseau and Wittemyer (2019).For this, we calculated crossvalidated Spearman rank correlations (r s ) between ten RSF-bin ranks and five test-training sets.We repeated this procedure 10 times to determine whether the r s differed from random.

Temporal variation in muskox calf recruitment and the linkages to the use of geochemical landscapes
To examine whether inter-annual variation in the relative use of geochemical landscapes influences calf recruitment, we related the output of the RSFs (i.e., annual selection coefficients for essential and non-essential elements) to muskox demography time series [updated from Schmidt et al. (2015)].For each year during the period 1996-2021, we calculated calf recruitment rates as the percentage of calves (<6 months old) per 100 adult 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).

Essential and non-essential elements in soil and vegetation
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).

Relative use of geochemical landscapes and vegetation types
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 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 those 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 s of 0.69 (sd = 0.11) across all models, which was consistently higher (p < 0.001) compared to the randomized model versions (Fig. S23).
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 lower during 2007-2014 compared to 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. 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 s of 0.71 (sd = 0.16), which was significantly higher (p < 0.001) than for a randomized model with a mean r s of 0.03 (sd = 0.33).

Relating muskox calf recruitment to relative use of geochemical landscapes and vegetation types
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.
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 %).
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 was negatively correlated to the relative use of sites with increased concentrations of As (p = 0.013) and Pb (p = 0.045).

Discussion
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.

Essential elements and muskox reproduction
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 − 1 d.w. and below the deficiency risk  threshold of 0.05 μg g − 1 d.w.proposed for free-ranging herbivores (Flueck et al., 2012).Similarly, the mean Cu concentration across all plants sampled in our study was 7.73 μg g − 1 d.w. and below the deficiency risk threshold of 10 μg g − 1 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).However, the mean Mo concentration in vegetation samples collected in this study (0.98 μg g − 1 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.At the same time, 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).Most ecological stoichiometry studies focus on the C:N ratio 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 (relatively) high N levels has a direct positive effect on calf recruitment.
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., 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.2012).

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 − 1 d.w.),As (0.12 μg g − 1 d.w.), and Hg (0.02 μg g − 1 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 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).Thus, continued monitoring of contaminant levels in Arctic vegetation and their potential direct and indirect effects on resident wildlife is pertinent.

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 elements in plants and soils vary between years and throughout the growing season and to determine the impacts by and on resident wildlife.Long-term herbivore exclosures provide a suitable experimental set-up to address such outstanding questions.
Our study relied on relatively small-scale (47m 2 ), 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 landscapes in wildlife ecology, we see particular value in the application of high-resolution biologging data (e.g., GPS collar and 3-dimensional accelerometers).Doing so will allow for the estimation of spatial and temporal variation in feeding behavior of individual animals across geochemical landscapes and multiple seasons (Nunes et al., 2022).Potential consequences of individual variation in habitat use and behavior on Darwinian fitness (i.e., survival and lifetime reproductive success) can then be assessed via long-term monitoring data or through Dynamic Energy Budget models (Desforges et al., 2019;Sperfeld et al., 2017).
Free-ranging animals are confronted by a myriad of environmental conditions that they need to balance to optimize fitness (Gaillard et al., 2010;Stein et al., 2014;van Beest and Milner, 2013).In our analyses, we focused on vegetation types, snow cover, and multiple essential and nonessential elements but discounted other potential factors (e.g., ambient temperature, density dependence, predation risk) influencing habitat use.We argue that this analytical approach is reasonable given that habitat use of muskoxen in this low-predator system is primarily influenced by snow conditions, plant productivity and access to forage biomass during spring and summer (Beumer et al., 2020;Tomassini et al., 2019).Indeed, muskoxen use this period primarily to replenish energetic reserves to survive the upcoming winter and to be able to develop a fetus (Desforges et al., 2021).This is also why we ran separate RSF and PCA analyses for the chemical elements (Fig. 4) and the vegetation types (Fig. 5) to facilitate comparison.The 'classic' PCA and RSF based on vegetation types explained more variation in annual recruitment (78.1 %) than the element-based analyses (67.5 %), which supports the idea that vegetation types are a better predictor of muskox habitat use and performance.Nonetheless, the element-based RSF and PCA also showed high predictive power and, moreover, provide insight into the geochemical mechanism underlying 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 results in the release of 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.

Conclusion
Our study clearly shows that the spatial 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 northeast Greenland.The distribution and bioavailability of the nonessential 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 alter the composition and spatial distribution of essential and non-essential elements.Therefore, continued monitoring of chemical element concentrations in Arctic soils and vegetation as well as responses by resident wildlife are critical to assess habitatperformance relationships under contemporary and future environmental change.While our study is specific to muskoxen roaming the high-Arctic tundra, we posit that the influence of geochemical landscapes on animal performance also applies to wildlife species living in lower latitude systems.

•
Use of geochemical landscapes impacts reproductive success of muskoxen.•Essential elements Cu, Se, and Mo posi-

Fig. 1 .
Fig. 1.A) Map of Greenland with Zackenberg valley indicated by the black square.B) 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.

Fig. 2 .
Fig. 2. Boxplots showing concentrations of chemical elements in soil (yellow) and vegetation (green) samples collected during July 2021 in the Zackenberg valley, northeast Greenland.Concentrations of C and N are in percentages, C:N is a ratio, and the remaining elements are in μg g − 1 d.w.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).

Fig. 3 .
Fig. 3. Predicted spatial distribution of chemical elements in vegetation across the Zackenberg valley in northeast Greenland.Concentrations of C and N are in percentages, C:N is a ratio, and the remaining elements are in μg g − 1 d.w.

Fig. 4 .
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.

Fig. 5 .
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.