Thomas Durantel

Gabriel Girard

Emmanuel Caruyer

Olivier Commowick

Julie Coloigner
email: julie.coloigner@irisa.fr

A Riemannian framework for incorporating white matter bundle priors in ODF

Keywords: Anatomical prior, Diffusion MRI, Track Orientation Distribution, Tractography, Tractometry

I. INTRODUCTION

D IFFUSION magnetic resonance imaging (dMRI) is an MRI modality that allows to measure the thermal agitation of the water molecule in the brain [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF]. This agitation being constrained by the tissues micro-structures, typically the nervous system axons, dMRI enables the voxelwise estimation of the orientations of the white matter fibers [START_REF] Dell | Reconstructing fiber orientations with diffusion mri[END_REF], [START_REF] Jones | Diffusion MRI: Theory, Methods, and Applications[END_REF]. By randomly choosing seeds from within the brain white matter, then following, from one voxel to the next, the local fiber orientations, one can achieve to estimate the brain connectivity and characterize the physical connections that mediate information transfer between cortical regions [START_REF] Jones | Diffusion MRI: Theory, Methods, and Applications[END_REF]. This process is called fibers tracking, or fiber tractography, and the resulting set of white matter trajectories is called tractograms [START_REF] Mori | Fiber tracking: principles and strategiesa technical review[END_REF]. The ability of those approaches to delineate the white matter fiber pathways offers unprecedented insight into the structural connections of the human brain and enormous potential for the study of human brain anatomy, development and function [START_REF] Jeurissen | Diffusion mri fiber tractography of the brain[END_REF]. Moreover, tractography has been proven particularly useful to neurosurgeons for the planning of surgery, especially to preserve important white matter pathways during resections [START_REF] Mancini | Automated fiber tract reconstruction for surgery planning: extensive validation in language-related white matter tracts[END_REF].

A multitude of new tracking algorithms has been proposed to improve the quality of the tractograms [START_REF] Schilling | Tractography dissection variability: What happens when 42 groups dissect 14 white matter bundles on the same dataset?[END_REF], [START_REF] Girard | On the cortical connectivity in the macaque brain: A comparison of diffusion tractography and histological tracing data[END_REF]. The first methods were deterministic methods, where only the principal diffusion orientation of the diffusion tensor was used in each voxel [START_REF] Basser | In vivo fiber tractography using DT-MRI data[END_REF]. However, studies have shown the difficulties that these methods have to represent complex brain regions, like crossing or fanning fibers [START_REF] Gao | Validation of dti tractography-based measures of primary motor area connectivity in the squirrel monkey brain[END_REF]. In response, new diffusion model, like the orientation distribution function (ODF) which encodes continuously all the diffusion orientations within a voxel and thus can characterize the complex fiber structure, have been developed [START_REF] Descoteaux | Deterministic and probabilistic tractography based on complex fibre orientation distributions[END_REF], [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion mri: Non-negativity constrained super-resolved spherical deconvolution[END_REF]. Modern methods have then sought to estimate the fiber dispersion estimation using probabilistic methods [START_REF] Jeurissen | Diffusion mri fiber tractography of the brain[END_REF], [START_REF] Behrens | Probabilistic diffusion tractography with multiple fibre orientations: What can we gain?[END_REF], [START_REF] Tournier | Mrtrix3: A fast, flexible and open software framework for medical image processing and visualisation[END_REF]. Among them, existing local tractography approaches estimates the local fiber orientation at each voxel independently by simultaneously fitting the local model and propagating in the most consistent direction, such as Kalman filtering method [START_REF] Gössl | Fiber tracking from dti using linear state space models: detectability of the pyramidal tract[END_REF], particle filtering method [START_REF] Pontabry | Probabilistic tractography using q-ball imaging and particle filtering: application to adult and in-utero fetal brain studies[END_REF], [START_REF] Stamm | Adaptive multimodal particle filtering for probabilistic white matter tractography[END_REF], random walks methods [START_REF] Descoteaux | Deterministic and probabilistic tractography based on complex fibre orientation distributions[END_REF] and graph theory method [START_REF] Sotiropoulos | Brain tractography using q-ball imaging and graph theory: Improved connectivities through fibre crossings via a modelbased approach[END_REF]. Other approaches have proposed more global approach which compute multiple fibers pathways and select the best ones based on the diffusion characteristics, like using Bayesian model [START_REF] Jbabdi | A bayesian framework for global tractography[END_REF] or the Hough transform [START_REF] Aganj | A hough transform global probabilistic approach to multiple-subject diffusion mri tractography[END_REF] and machine learning more recently [START_REF] Poulin | Tractography and machine learning: Current state and open challenges[END_REF]. However, despite the improvement of the dMRI acquisition, the modeling and the tracking, white matter fibers tractography still can't overcome some limitations [START_REF] Jbabdi | Tractography: where do we go from here?[END_REF], [START_REF] Maier-Hein | The challenge of mapping the human connectome based on diffusion tractography[END_REF] which prevents those methods to be usable in the clinical context. Maier-Hein et al (2017). have shown that most modern algorithms can recover up to 90% of the ground truth white matter bundles, but these are poorly recovered in terms of spatial extent [START_REF] Maier-Hein | The challenge of mapping the human connectome based on diffusion tractography[END_REF]. Moreover, those tractograms also contained large amounts of invalid bundles [START_REF] Maier-Hein | The challenge of mapping the human connectome based on diffusion tractography[END_REF]. Overall, the iFOD2 [START_REF] Tournier | Improved probabilistic streamlines tractography by 2nd order integration over fibre orientation distributions[END_REF] and PTT [START_REF] Aydogan | Parallel transport tractography[END_REF] algorithms had good performances in recent international challenges [START_REF] Maier-Hein | The challenge of mapping the human connectome based on diffusion tractography[END_REF], [START_REF] Rafael-Patino | The diffusion-simulated connectivity (disco) dataset[END_REF].

To overpass the poor spatial extent of white matter bundles, one proposed solution is to increase the numbers of the total generated streamlines leading to an over representation of easy tracks and an under-representation of difficult tracks, inducing a density bias [START_REF] Rheault | Visualization, interaction and tractometry: dealing with millions of streamlines from diffusion mri tractography[END_REF]. For example, the local tractography methods aim to recover a global connectivity by inferring from locals directional information that causes the algorithm to follow, in complex regions, the easiest path available and sometimes the non-realistic decisions [START_REF] Maier-Hein | The challenge of mapping the human connectome based on diffusion tractography[END_REF]. Various studies have suggested that advanced diffusion microstructure modeling [START_REF] Daducci | Microstructure informed tractography: pitfalls and open challenges[END_REF], [START_REF] Smith | Sift2: Enabling dense quantitative assessment of brain white matter connectivity using streamlines tractography[END_REF], streamline filtering techniques [START_REF] Daducci | Accelerated microstructure imaging via convex optimization (amico) from diffusion mri data[END_REF] or advances in machinelearning-driven tractography [START_REF] Neher | A machine learning based approach to fiber tractography using classifier voting[END_REF] could make improvement on the false positive rate.

Other approaches have proposed to compensate the lack of global information by adding them in the form of anatomical priors whose purpose would be to guide the algorithm in complex regions. As for segmentation and label fusion approach, recent methods are based on the use of an anatomical atlas such as the TRACULA method [START_REF] Yendiki | Automated probabilistic reconstruction of white-matter pathways in health and disease using an atlas of the underlying anatomy[END_REF], with promising results combining deterministic tractography and anatomical prior. Then, diffusion priors during the tracking process were proposed in [START_REF] Cook | Atlas-guided probabilistic diffusion-tensor fiber tractography[END_REF]- [START_REF] Rheault | Bundle-specific tractography with incorporated anatomical and orientational priors[END_REF], improving the delineation of white matter bundles. In [START_REF] Rheault | Bundle-specific tractography with incorporated anatomical and orientational priors[END_REF], a bundle-specific method incorporates anatomical and orientational priors based on a template, to improve the reconstruction of long fibers and increase of reproducibility, sensitivity and specificity. On the other hand, in [START_REF] Wasserthal | Tractseg-fast and accurate white matter tract segmentation[END_REF], another well-known approach based on a machine learning method was proposed to automatically segment, with high precision, the overall shape of a bundle, being a potentially great source of prior.

Based on these promising results, we developed a method of anatomical prior creation and combination, which is usable with any tractography algorithm based on orientation distribution function (ODF) [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF]. Our method use pre-segmented fiber bundles by [START_REF] Wasserthal | Tractseg-fast and accurate white matter tract segmentation[END_REF], to agglomerate global information from several different brains and capture the orientational variability in complex brain region. In this work, prior are computed on this anatomical atlas and expressed in the form of voxel-wise TOD [START_REF] Dhollander | Track orientation density imaging (todi) and track orientation distribution (tod) based tractography[END_REF] and then combined with the orientation distribution function (ODF) [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF] data using a Riemannian framework [START_REF] Goh | A nonparametric riemannian framework for processing high angular resolution diffusion images and its applications to odf-based morphometry[END_REF]. We decided to incorporate those priors in two state-of-the-art algorithms using traditional ODF data, MRtrix iFOD2 [START_REF] Tournier | Improved probabilistic streamlines tractography by 2nd order integration over fibre orientation distributions[END_REF] and Trekker PTT [START_REF] Aydogan | Parallel transport tractography[END_REF]. We then evaluated those methods on the diffusion-simulated connectivity (DiSCo) dataset [START_REF] Rafael-Patino | The diffusion-simulated connectivity (disco) dataset[END_REF] and on the Human Connectome project (HCP) 1 data to show the increase quality of the tractogram and more specifically on the spatial extent of the reconstructed bundles.

II. METHODS

Our method can be separated into 3 distinct parts: i) the construction of an anatomical atlas from segmented fibers, considered as gold standard, ii) the extraction and estimation 1 https://www.humanconnectome.org/ of the TOD anatomical priors, from the atlas and iii) the combination between the priors and the subject data. The entire framework is illustrated in Figure 1.

A. Atlas construction

The first step of our framework is the creation of a reference fiber atlas from a set of high-resolution diffusion images (see Figure 1.a). In this paper, we chose high-resolution diffusion images from the Human Connectome Project (HCP) young adult study [START_REF] Alexander | Diffusion tensor imaging of the brain[END_REF], which were acquired at high spatial resolution and b-values with 90 directions and thus contains high quality information on the fibers orientations.

The proposed atlasing method follows a modified version of [START_REF] Guimond | Average brain models: A convergence study[END_REF]. [START_REF] Guimond | Average brain models: A convergence study[END_REF] which was adapted for diffusion data, to compute an atlas of Tensor images from a set of control subjects. This approach is based on a procedure which compute iteratively the atlas by registering the tensor images of HCP images onto a current reference. After each iteration, it performs an average of all the previous unbiased atlases to compute the next one at the following iteration. The main modifications from the original [START_REF] Guimond | Average brain models: A convergence study[END_REF] method is to use diffeomorphisms encoded as Stable Vector Fields (SVF) and the log-Euclidean framework [START_REF] Arsigny | A log-euclidean framework for statistics on diffeomorphisms[END_REF] to compute the average transformations and the approach was adapted to diffusion data. The entire procedure is detailed in [START_REF] Suarez | Automated delineation of white matter fiber tracts with a multiple region-ofinterest approach[END_REF].

For each HCP data, 72 bundles were obtained by Wasserthal et al. (2018), using the TractSeg algorithm [START_REF] Wasserthal | Tractseg-fast and accurate white matter tract segmentation[END_REF] which perform fiber tracking and fiber filtering, first by regions of interest (ROIs), then by an expert. The overall process used to segment these bundles is detailed in [START_REF] Wasserthal | Tractseg-fast and accurate white matter tract segmentation[END_REF] and the data is available online [START_REF] Maier-Hein | High quality white matter reference tracts (1.2.0)[END_REF].

The transformation field obtained by the atlasing method from each individual space to the atlas space are then applied to 72 segmented fiber bundles of interest of each subject in order to all align them to the fiber atlas.

B. Priors Extraction

Once aligned to a common space, the bundles of interest are combined in order to represent the general shape of bundles. In order to capture complex orientation, a local fibre orientation prior was estimated using the track orientation distribution (TOD) [START_REF] Dhollander | Track orientation density imaging (todi) and track orientation distribution (tod) based tractography[END_REF], in each voxel. The tractogram is represented as a voxel-wise probability distribution function (PDF) in the image domain, rather than a set of individual tracks (samples from this distribution). The TOD, thus, captures the expected fibre directions [START_REF] Dhollander | Track orientation density imaging (todi) and track orientation distribution (tod) based tractography[END_REF].

Our TOD imaging method operates as follows : first, all the fibers directions of the voxel are extracted, then clustered using a k-means algorithm to define the main directions within the voxel and to correct for the density bias. This way, up to 4 principal directions can be extracted, allowing us to correct for the density bias as well as to characterize complex crossing regions. Then, the TOD in each voxel is represented using a set of modified spherical harmonics (SH) basis functions (see [ [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF]]) and constructed by projecting one spherical point spread function (PSF) per extracted main direction along the z direction. The PSF, along the z direction at position ϵ, can be easily obtained in SH basis:

δ z(ϵ) (θ, ϕ) = ∞ l=0 l m=-l c m l Y m l (θ, ϕ), (1)
where 0 ≤ θ ≤ π and -π ≤ ϕ ≤ π are the spherical coordinates, Y m l the spherical harmonic of degree l and order m and c m l , the coefficients, given as:

c m l = π -π π 0 δ z(ϵ) (θ, ϕ)Y m l (θ, ϕ)dθdϕ (2)
We found that resolving the equation 2 with the PSF pointing towards the direction θ = π 2 and then rotating them to match the direction extracted from the atlas greatly simplifies the calculations. In [START_REF] Geng | Diffusion mri registration using orientation distribution functions[END_REF], authors showed a simple method to rotate functions expressed in SH basis. Since 4 different directions can be represented in a voxel, the final step of the prior extraction is to average the projected PSF. The averaging of distributions defined on the sphere is performed through a Riemannian framework [START_REF] Goh | A nonparametric riemannian framework for processing high angular resolution diffusion images and its applications to odf-based morphometry[END_REF]. To do that, we define ψ(θ, ϕ) = δ z(ϵ) (θ, ϕ), the square-root density function of the PDF δ(θ, ϕ). The square-root is used to ensure that the logarithm maps are available in closed form. Thus, δ(θ, ϕ) has to be strictly positive. For this reason, in order to find the right PSF, several different distributions were tested (i.e. a Dirac distribution, a Watson distribution and a Gaussian distribution) and then were expressed on the SH basis. We found out that a 2D Gaussian distribution, defined on the unit sphere, is the only one whose number of negative values do not diverge with the degree l.

Then, in order to calculate the average distribution, we used the weighted Karcher mean, ψ, of a set of n points {ψ i } n i=0 in a Riemannian manifold defined by:

ψ = argmin 1 2 n i=0 ω i dist(ψ, ψ i ) 2 (3)
with ω i ≤ 0 and n i=0 ω i = 1 As described in [START_REF] Goh | A nonparametric riemannian framework for processing high angular resolution diffusion images and its applications to odf-based morphometry[END_REF], the unique solution is ψ, such as:

n i=0 ω i log ψ (ψ i ) = 0 (4)
where log ψ is the logarithm map from ψ to ψ i , given by:

log ψ (ψ i ) = ψ i -⟨ψ i , ψ⟩ψ 1 -⟨ψ i , ψ⟩ cos -1 ⟨ψ i , ψ⟩ (5)
with ⟨., .⟩, the normal dot product.

At the end of this step, using the equation 5 and ω i = 1/n, ∀i ∈ {1, ..., n}, we obtained ψ, which is the square-root of the prior for each voxel. An example of the prior with the SH coefficients, truncated at a degree of l = 8, of the PSF described previously are shown in Figure 4. The prior creation process in order to average all extracted directions has to be calculated only once.

C. Data and prior combination

As displayed in Figure 1, the ODF image of the individual diffusion dataset is calculated to estimate the next direction. During tractography, in order to inform the ODF with the anatomical prior, as previously explained in the previous section, a weighted Karcher mean is computed between the ODF and the prior using equation 4, to obtain the enhancedodf (EODF). However, this method involves calculating, at each step of the tracking, the dot product between distribution on the sphere, that implies an integral over the sphere which would complicate the process. To simplify the computations, Goh et al. (2009) proposed to discretize the PDF and to work with the square root of histograms [START_REF] Goh | A nonparametric riemannian framework for processing high angular resolution diffusion images and its applications to odf-based morphometry[END_REF]. In this case, the dot product become summations.

Since the overall purpose of our method is only to guide the tractography algorithms and not to reflect the structure of the priors, the weighting factor ω i must be well studied. We want to be able to distinguish between the simple linear region, where not much guidance is needed, especially in easier-totrack bundles with only one fiber direction, and the complex regions with crossing fibers, where the use of the prior is more needed. In order to address this specification, we choose to use 2 measures to weight the prior: the generalised fractional anisotropy (GFA), for the prior and the Akaike information criterion (AIC) [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] for the data. Left: DiSCo dataset pathways. Right: Prior weighting map example for the DiSCo dataset. Higher value denotes higher usage of the prior and thus more complex regions.

The GFA, is given by:

GF A = 1 -c 2 0 N i=0 c 2 i (6)
where the c i are the SH coefficients. The GFA is a measure of the anisotropy within the considered voxel, that is, in the case of the prior, this value can reflect the complexity of a region.

The AIC, is described by;

AIC = 2k -2 ln(L) (7)
where k is the numbers of parameters in the model and L the maximized likelihood. The AIC is a measurement of the estimation quality of a model. Empirically, to well differentiate linear and complex regions, we defined the prior weight given by:

ω prior = α(1 -GF A) + βexp( AIC min -AIC 2 ) (8)
where α represents the amount of information extracted from the priors, via the GFA, and β the amount of information extracted from the data, via the AIC. Those parameters needs to be adapted to the data. For the scenario presented in this work, we choose α = 0.35 and β = 0.65 (see Figure 3). 8) for the tracking of the CST and for the iFOD and PTT algorithms. The method performs most effectively with an alpha value of 0.35.

An example of a weighting map for the DiSCo dataset can be seen on Figure 2 with higher ω prior in crossing fiber areas.

The result of this averaging process is then an EODF, expressed in SH basis, usable in any ODF-based tractography algorithm. Figure 4 shows examples of data and priors combination, as well as the effect of weighting on the results. The first row represents a voxel where the ODF and the prior have the same direction, corresponding to an easy-to-track regions with one-way crossing. On the other hand, in the second and third rows, the orientation of the ODF and the prior are different. In this case, according to the value of ω prior , the EODF is almost equal to the ODF or a mixture of the prior and ODF.

III. EXPERIMENTS

In order to validate our method using priors, we conducted our first experiments on the DiSCo challenge dataset [START_REF] Rafael-Patino | The diffusion-simulated connectivity (disco) dataset[END_REF]. Then, we tested our algorithms on HCP data. For each experiment, we implemented our prior approach on two state-ofthe-art algorithms, MRtrix iFOD2 [START_REF] Tournier | Improved probabilistic streamlines tractography by 2nd order integration over fibre orientation distributions[END_REF] and Trekker PTT [START_REF] Aydogan | Parallel transport tractography[END_REF]. Thus, on the two dataset, 4 tractography reconstruction were performed, both with and without using priors. For iFOD2, we used a step size of 0.2 voxel and a maximum angle between successive steps of 20 degrees. For PTT, we used the default parameters.

A. DiSCo challenge dataset

The aim of this experiment is to test the ability of our novel method to improve the connectivity estimation and the effect of noise. In this context, we do not have access to several datasets to build an anatomical atlas, instead, we computed the TOD: (i) using the ground-truth fiber pathways as anatomical priors and (ii) using fiber tracked on high quality data as anatomical priors.

1) Data: We performed two different variants of this experiment: in the first one, priors are estimated on ground-truth fibers, not available on in vivo data; and the second experiment is used to demonstrate the benefit of our method, without having a ground-truth but priors build on high resolution data. in clinical context, we used degraded data with a grid of 20×20×20 voxels and only using the 3000 s/mm 2 b-shell with 60 directions. The prior is computed with the original data. In addition, to study the effect of noise, we added Rician noise on the dMRI. The resulting images have an average SNR of 0dB, 10dB, 20dB, 30dB, 40dB and 50dB. The individual ODF were computed using the method described in [START_REF] Aganj | Reconstruction of the orientation distribution function in single-and multiple-shell q-ball imaging within constant solid angle[END_REF] and represented in a spherical harmonic basis truncated at the 8 th order. Then, EODF were calculated using the method described in this paper.

(ii) multi-shell ODF: In this second set of experiment, we computed the priors on fibers tracked with the multi-shell multi-tissue constrained spherical deconvolution (CSD) fiber ODFs. The idea is to simulate a clinical context by computing the ODF on clinical dMRI data and combined with anatomical priors calculated on high resolution data. The anatomical prior are calculated on 40 × 40 × 40 voxels DiSCo dMRI at SNR 20. From those images, fiber ODF were computing using the MRtrix implementation of the multi-shell multi-tissue constrained spherical (msmt-CSD) method described in [START_REF] Jeurissen | Multi-tissue constrained spherical deconvolution for improved analysis of multi-shell diffusion mri data[END_REF] using default parameters and using 8 th order SH basis. Then, fibers were tracked using the MRtrix iFOD2 algorithm with a step size of 0.2 voxel and default parameters. Only the streamlines that reach the endings ROIs were kept. Finally, the TODs were extracted from those fibers using the method detailed in the previous section. After prior estimation on high resolution data, we estimated the ODF of a subject data with 20 × 20 × 20 grid and SNR= 20. Then, EODFs were obtained by combining the TOD with the clinical data fiber ODF.

2) Metrics: To quantify the results, we computed the Pearson correlation coefficient (r) between the ground truth connectivity matrices of the three DiSCo dataset and the resulting tractograms connectivity matrix. These matrices are computed by counting the numbers of streamlines that reach both start and end ROIs for each DiSCo fibers bundles.

B. HCP data 1) Data:

To test our method on in vivo high quality data, we used the HCP young adult data, acquired with 90 gradients distributed on 3 shells of b=1000, 2000, and 3000 s/mm 2 with 6 b= 0 acquisitions and a final resolution of 1.25mm 3 . Among them, the 105 pre-segmented HCP images were used, as described in II-A, to build fiber atlases of the Arcuate fascicle (AF), the Cingulum left (CG), the Corticospinal tract (CST), the Optic radiation (OR) and the Superior longitudinal fascicle I (SLF I) to study the influence of the addition of prior on bundles of different degrees of complexity and in different region of the brain. A cross validation was performed on the 105 subjects in which 100 subjects are used to build the atlas and the 5 remaining to perform the tractography algorithm. This process is repeated 8 times, giving us a total of 40 subjects to test the method on. For each image belonging to training set, ODFs were computed on the shell b= 3000 s/mm 2 using the method described in [START_REF] Aganj | Reconstruction of the orientation distribution function in single-and multiple-shell q-ball imaging within constant solid angle[END_REF] and represented in a spherical harmonic basis truncated at the 8 th order. Then, the TOD images for each bundle were registered on the subject space. Those priors were incorporating during the tracking process for each subject. The tractography algorithms were tested on the 5 remaining subjects.

In order to compare our method with state of the art anatomical priors methods, we used the BST algorithm ([ [START_REF] Rheault | Bundle-specific tractography with incorporated anatomical and orientational priors[END_REF]], using the default parameters, to compute enhanced ODF for the CST and the AF bundles. We then tracked those bundles with the iFOD2 and PTT algorithms using the same parameters as before.

2) Metrics: To quantify the overall shape quality of the tractograms, we computed the generalized Dice score [START_REF] Crum | Generalized overlap measures for evaluation and validation in medical image analysis[END_REF] between the segmentation of fibers obtained with tractography algorithm and the reference fibers. In order to measure the improvement on the numbers of streamlines reaching both ROI endpoints, we computed the percentage of streamlines that correctly connects the beginning and the end regions of bundles over the total number of streamlines in the reference bundle (noted valid streamlines VS). It is worth noting that this measure does not mean anything in absolute, being biased by the numbers of streamlines in the reference tracks. However, since the same numbers of streamlines is generated for each bundle, it allows us to compare the overall bundles shape and quality with and without the addition of priors.

IV. RESULTS

A. DiSCo challenge dataset

Figure 6 displays the correlation with the ground-truth connectivity of the two experiments, for the three DiSCo dataset. First, on Figure 6.a, we observed that lower the noise is, better are the correlation improvement for the experiment using the fiber pathways as ground-truth. Next, both algorithms performed better with the addition of priors for all level of noise and the three datasets. We can also notice that the improvement decreases when the quality of the data increases. In average, the correlation increases of 0.23 for an SNR of 10 dB but only 0.11 for an SNR of 50 dB. Moreover, the mean correlation improvement is of 0.19 and 0.17 for the iFOD2 and PTT algorithms, respectively. At last, the average maximum reached is r = 0.83 for PTT and r = 0.82 for iFOD2. There are no major differences between the three dataset. For the second experiment using high resolution DW-MRI data, we showed an average correlation improvement of 0.15 for the PTT algorithm and 0.17 for the iFOD2 algorithm, in figure 6.b. Incorporating anatomical priors calculating on high resolution data during the tracking process of low resolution data increases the Pearson correlation coefficient with the ground-truth connectivity matrix. For the Valid streamlines score, the addition of priors appears to always increase, on average, the numbers of streamlines that connect both end regions. However, enhanced tractography seems to also increase the variability of the results, as we found a mean variance across all bundles of 7.19% without the prior and 19.84% with. We also noticed that incorporating the priors in the PTT algorithm improve more the Valid streamlines than with the iFOD2 algorithm. Indeed, the average gain for iFOD2 is of 25.08% and 31.37% for PTT. The same observation can be made for the dice measurement, but as opposed to the Valid streamlines measurement, the variability seems to be decreased with the addition of prior. We notice that incorporating the anatomical priors improve the spatial coverage and that a higher fraction of streamline reach the endpoints of the bundle

B. HCP data

V. DISCUSSIONS

In this work, we proposed a novel method for creating and incorporating anatomical priors to any ODF-based tractography algorithms. We showed that incorporating anatomical priors on two state-of-the-art tractography algorithms improve the overall quality of tractograms when priors are extracted from a ground truth, but also from high quality tractograms. Indeed, first, on the DiSCo dataset, the proposed prior-based tractography methods obtained better correlation scores between the reconstructed bundles and the ground truth for every level of noise than the standards probabilistic algorithms, but this improvement decreases with the increasing SNR. Thus, the addition of priors to ODF images improve the quality of tractograms in presence of noise. However, due to the actual form of the extracted TOD it may not be useful on high SNR data. More particularly, in the DISCO experiment, we tried to mimic real clinical context with low spacial resolution data and an SNR of 20. In this case, our results showed a better correlation for two methods incorporating priors between the estimated fiber and the ground-truth, reducing false positive fibers. This could allow us to consider using this method on clinical data of average quality, informed by prior extracted on high quality data. On the HCP data, the same results were obtained with an increase of the fraction of valid streamlines numbers when anatomical priors were incorporating during tracking process. In the same time, due to the increase of the dice score when adding prior, the overall shape of the bundle seems to be also modified. We also reported an increase in variability across subjects of the tractography with priors compared with the ones without priors.

The effect of the priors seems to be more effective in complex fibers configurations, like fanning and crossing fibers or in large curvature radius. Indeed, our method is able to a better delineation and estimation of precise details. For the CST, it is well known that the upper part, that fans into the cortex, is hard to estimate. Without prior, the iFOD2 and PTT algorithms only manage to estimate a portion of the CST. While, with the addition of anatomical prior, they achieve a more complete reconstruction of the fanning portion. Another example is the OR bundle, where the Meyer's loop, with its highly curving fibers, is still a tractography open challenge [START_REF] Bertani | Optic radiation diffusion tensor imaging tractography: An alternative and simple technique for the accurate detection of meyer's loop[END_REF]. When guided by EODF, both tested algorithms show results that, while containing spurious streamlines, allow a better recovery of the anterior extent of the Meyer's loop. This result can be easily explained by the fact that in regions with crossing fibers, more information is taken from the prior. See the right column in Fig. 5 where the brightest correspond to the most complex fibers micro structure regions and denotes a higher usage of the prior.

It is important to note that the results presented in this paper use iFOD2 and PTT algorithms as proof of concept, but the main idea behind this work is to develop a method to improve the spatial coverage and the numbers of streamline reaching the endpoints of the bundle of any ODF-based tractography algorithm. In fact, future works should study the contribution of prior addition with other algorithm but also generalize this method to other diffusion model.

In [START_REF] Rheault | Bundle-specific tractography with incorporated anatomical and orientational priors[END_REF], authors generated sharp TOD by projecting Dirac in the SH basis. However, our experiences have shown that this also generate a number of negative values which increase with the truncated SH basis order, due to the Gibbs phenomenon. Thus, to overpass this problem and remove the negative values, they use apodized delta function [START_REF] Raffelt | Reorientation of fiber orientation distributions using apodized point spread functions[END_REF] proposed in [START_REF] Dhollander | Track orientation density imaging (todi) and track orientation distribution (tod) based tractography[END_REF], but at the cost of a loss of angular resolution. Using this PSF, we would have to truncated at a higher order the TOD SH basis to obtain the same angular resolution that we have with the 8 th order of the Gaussian PSF [START_REF] Dhollander | Track orientation density imaging (todi) and track orientation distribution (tod) based tractography[END_REF]. So we decided to keep the Gaussian PDF and the ODF estimated without CSD, at the expense of less sharp ODF peaks, our interest here being to compare the tracking connectivity with and without prior. In the same way, in [START_REF] Rheault | Bundle-specific tractography with incorporated anatomical and orientational priors[END_REF], priors and data are combined by an element-wise multiplication between the two sets of SH coefficients, followed by a normalization. However, our approach (i.e. using a Riemannian framework), allows us a better control over the weighting between the TOD and ODF and minimizes the potential angular smoothing. It also enables, for future works, more complex operation on TOD, such as interpolation between different set of priors or even TOD filtering to remove unnecessary information, in order to improve the prior quality, at the cost of more approximation during the TOD estimation and an increased complexity.

Other works that studies the idea of guiding tractography algorithms, like [START_REF] Yendiki | Automated probabilistic reconstruction of white-matter pathways in health and disease using an atlas of the underlying anatomy[END_REF], [START_REF] Cook | Atlas-guided probabilistic diffusion-tensor fiber tractography[END_REF], [START_REF] Yap | Poptract: populationbased tractography[END_REF] or [START_REF] Rheault | Bundle-specific tractography with incorporated anatomical and orientational priors[END_REF], use non-linear registration and atlasing techniques during the atlas, or template, creation process, preferring to do an average of the references images. This could smooth the variability across the atlas and could also be the source of errors that would be accumulated in the final atlas and would, in the end, inject a bias inherent to the atlas in the tractography process. In our approach, using a modified Guimond method [START_REF] Guimond | Average brain models: A convergence study[END_REF], the references image are, at each step, iteratively registered with a non-rigid transformation onto the current space that become the reference space in the next step ant it is only the average of the transformation that is computed at the end of each step. With this approach, all the variability information is accumulated in the atlas and the errors are not included in the final atlas. Thus, no bias is introduced in the tractography.

Also, using pre-segmented fibers to build the atlas may not be appropriate in order to study the global connectivity of the brain, as this approach only allows the enhanced tracking to be performed within a unique bundle. But, although not showed here, our method also permits, by concatenating the individual bundles priors, the construction of full brain atlases, and thus enable guidance of whole brain tracking. Therefore, it could be used for clinical studies where priors would be constructed from high resolution dMRI data in order to guide tractography algorithm on poor, clinical, resolution data.

Finally, our method has proved to be able to increase the quality of the estimated fibers, in term of spatial extant and number of valid connection when utilized on healthy subjects. These promising results could already enable advances in an academic context for the study of the healthy brain. In a more clinical context, if used on patients suffering from pathologies that only slightly or moderately modify the white matter, such as psychiatric pathologies, this method could also bring improvements in fiber estimation, perhaps through a more in-depth study of α and β parameters. Research along these lines should be pursued in future work. However, when working on patient suffering from severe brain alteration, such as strokes, this method could produce a tractogram solely guided by the prior in the affected region, thus removing the specific microstructural modification and therefore may not be the most appropriate approach. In such scenarios, a solution might be to import other types of information into the priors. Either in the way the data and priors are weighted, using for example Apparent fiber density [ [?]] of the data to potentially detect these brain modifications and thus prevent the fibers from being tracked solely on the priors, or directly in the estimation of the priors using other imaging modalities, such as myelin-sensitive relaxometry in order to obtain priors images describing these brain regions in greater detail.

VI. CONCLUSION

In this paper, we developed a method of anatomical prior creation and combination, which is usable with any ODFbased tractography algorithm. The prior are computed on fiber atlases and expresses in the form of TOD, in order to characterize the brain variability and multiple fiber directions. Then, our prior is incorporated to guide the ODF-based tractography algorithm. Based on our results on DISCO and HCP data, incorporating our anatomical priors improve the fibers reconstruction, in terms of spatial extent and valid streamlines, especially in crossing fiber regions. Furthermore, we have shown that our method outperforms previous methods using anatomical priors in the tracking of the CST and the AF. Moreover, our approach could also greatly enhance the tractography in the context of clinical data, by incorporating priors estimated on high quality data, which could help for the study of neurological diseases.

Fig. 1 .

 1 Fig. 1. Overall method pipeline. (a) Atlas creation. (b) TOD estimation. (c) Data and prior combination.

Fig. 2 .

 2 Fig. 2.Left: DiSCo dataset pathways. Right: Prior weighting map example for the DiSCo dataset. Higher value denotes higher usage of the prior and thus more complex regions.

Fig. 3 .

 3 Fig. 3. Evolution of the Dice score versus α (from the Equation8) for the tracking of the CST and for the iFOD and PTT algorithms. The method performs most effectively with an alpha value of 0.35.

Fig. 4 .

 4 Fig. 4. Examples of data and prior combination. On the first row: prior and data are aligned. The EODF main orientation are unchanged. On the second and third rows: Prior and data are crossing. The corresponding EODF peak is more or less attenuated using the weight of the prior.

  (i) Fiber pathways ground-truth: The DiSCo dMRI images have a grid of 40 × 40 × 40 with voxel size of 1 mm 2 and are acquired with 4 different b-shell of 1000, 2000, 3000 and 1300 s/mm 2 , 90 directions per b-shell and 4 b0. In order to obtain synthetic data with the same quality as those

Fig. 5 .

 5 Fig. 5. Visual comparison of the resulting tractograms, with and without prior. From left to right, in column : iFOD2, PTT, iFOD2 with priors and PTT with priors. Last column on the right: Prior weighting map, where brighter pixel denotes a higher usage of the prior. From top to bottom : The Arcuate fascicle (AF), the Cingulum (CG), the Corticospinal tract (CST), the Optic radiation (OR) and the Superior longitudinal fascicle I (SLF I).

Figure 5

 5 Figure 5 displays the segmentation of AF, CG, CST, OR and SLF I obtained with iFOD2 and PTT, with and without prior. The results for the Valid streamlines and the Dice score are presented, respectively, in Figure ?? and Figure ??.For the Valid streamlines score, the addition of priors appears to always increase, on average, the numbers of streamlines that connect both end regions. However, enhanced tractography seems to also increase the variability of the results, as we found a mean variance across all bundles of 7.19% without the prior and 19.84% with. We also noticed that incorporating the priors in the PTT algorithm improve more the Valid streamlines than with the iFOD2 algorithm. Indeed, the average gain for iFOD2 is of 25.08% and 31.37% for PTT. The same observation can be made for the dice measurement, but as opposed to the Valid streamlines measurement, the variability seems to be decreased with the addition of prior. We notice that incorporating the anatomical priors improve the spatial coverage and that a higher fraction of streamline reach the endpoints of the bundle

Fig. 6 .

 6 Fig. 6. Pearson correlation for the DiSCo experiments: (a) Fiber pathways ground-truth experiment with respect to SNR (b) Multi-Shell ODF experiment for the 3 DiSCo dataset. First column for PTT and second for iFOD2. In all cases, blue points denotes measure without prior and the oranges ones, with prior.

Fig. 7 .Fig. 8 .

 78 Fig. 7. On the left, Dice score result for the HCP experiment. On the right, VS number result for the HCP experiment. From top to bottom : AF, CG, CST, OR, SLF I


Floris M Van Beest

Niels Martin Schmidt

Laerke Stewart

Lars H Hansen

Anders Michelsen

Jesper B Mosbacher

Hugo Gilbert

Gaël Le Roux

Sophia V Hansson

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

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

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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., [START_REF] Mcloughlin | Lifetime reproductive success and composition of the home range in a large herbivore[END_REF]. Quantifying and understanding so called habitat-performance relationships [START_REF] Gaillard | Habitat-performance relationships: finding the right metric at a given spatial scale[END_REF][START_REF] Matthiopoulos | Establishing the link between habitat selection and animal population dynamics[END_REF] is critical for a range of topics in environmental management including assessments of climate and land-use change impacts [START_REF] Verniest | Coupling future climate and land-use projections reveals where to strengthen the protection of Mediterranean Key Biodiversity Areas[END_REF]. Habitat-performance relationships can be estimated by linking differences in animal performance to variation in use of resource abundance or quality [START_REF] Felton | Varied diets, including broadleaved forage, are important for a large herbivore species inhabiting highly modified landscapes[END_REF][START_REF] Ofstad | Opposing fitness consequences of habitat use in a harvested moose population[END_REF]. 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: [START_REF] Mosser | Serengeti real estate: density vs. fitness-based indicators of lion habitat quality[END_REF] or by measuring digestibility of plant material (in the case of herbivores: [START_REF] Proffitt | Linking landscape-scale differences in forage to ungulate nutritional ecology[END_REF]. 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 [START_REF] Boersma | Too Much of a Good Thing: On Stoichiometrically Balanced Diets and Maximal Growth[END_REF][START_REF] Sperfeld | Bridging Ecological Stoichiometry and Nutritional Geometry with homeostasis concepts and integrative models of organism nutrition[END_REF].

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 [START_REF] Elser | Biological stoichiometry: a chemical bridge between ecosystem ecology and evolutionary biology[END_REF][START_REF] Trueman | Ecogeochemistry potential in deep time biodiversity illustrated using a modern deep-water case study[END_REF]. Indeed, it is now well established that a balanced intake of macro elements such as nitrogen (N) and phosphorus (P) can positively affect growth [START_REF] Plath | Mineral Limitation of Zooplankton: Stoichiometric Constraints and Optimal Foraging[END_REF] and reproduction [START_REF] Durso | Stoichiometric and stable isotope ratios of wild lizards in an urban landscape vary with reproduction, physiology, space and time[END_REF] 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 [START_REF] Feinberg | Reductions in the deposition of sulfur and selenium to agricultural soils pose risk of future nutrient deficiencies[END_REF][START_REF] Kaspari | Biogeochemistry and geographical ecology: Embracing all twenty-five elements required to build organisms[END_REF]. Adverse effects may arise, however, if homeostasis is disrupted i.e., when reserves of essential elements are outside required physiological ranges [START_REF] Bhattacharya | Nutritional aspects of essential trace elements in oral health and disease: an extensive review[END_REF][START_REF] Kabata-Pendias | Trace Elements in Abiotic and Biotic Environments[END_REF]. 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 [START_REF] Burger | Assessment and management of risk to wildlife from cadmium[END_REF][START_REF] Massányi | Effects of Cadmium, Lead, and Mercury on the Structure and Function of Reproductive Organs[END_REF].

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 [START_REF] Brown | Mineral surfaces and bioavailability of heavy metals: A molecular-scale perspective[END_REF]. For terrestrial herbivores, the intake of essential and non-essential elements is primarily influenced by their bioavailability in vegetation, the underlying soils, and plantspecific 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 [START_REF] Flueck | Effect of trace elements on population dynamics: Selenium deficiency in free-ranging black-tailed deer[END_REF][START_REF] French | Geochemical landscapes as drivers of trace and toxic element profiles in wild red deer (Cervus elaphus)[END_REF]. While assessments of elemental status in free-ranging herbivores and the linkages to wildlife health are rapidly increasing [START_REF] Jutha | Evaluating the use of hair as a non-invasive indicator of trace mineral status in woodland caribou (Rangifer tarandus caribou)[END_REF][START_REF] Mosbacher | Hair mineral levels as indicator of wildlife demographics?-a pilot study of muskoxen[END_REF][START_REF] Rioux | Links between individual performance, trace elements and stable isotopes in an endangered caribou population[END_REF][START_REF] Squadrone | Wildlife hair as bioindicators of metal exposure[END_REF], 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 [START_REF] Balluffi-Fry | Quantity-quality trade-offs revealed using a multiscale test of herbivore resource selection on elemental landscapes[END_REF][START_REF] Leroux | Stoichiometric distribution models: ecological stoichiometry at the landscape extent[END_REF][START_REF] Rizzuto | Forage stoichiometry predicts the home range size of a small terrestrial herbivore[END_REF][START_REF] Sach | Spatial geochemistry influences the home range of elephants[END_REF] 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 habitatperformance 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 [START_REF] Mosbacher | Hair mineral levels as indicator of wildlife demographics?-a pilot study of muskoxen[END_REF][START_REF] Rombach | Utilization of copper during lactation and neonatal development in muskoxen[END_REF] and other herbivores [START_REF] Flueck | The implications of selenium deficiency for wild herbivore conservation: a review[END_REF][START_REF] Michaluk | Involvement of copper in female reproduction[END_REF]. 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).

Material and Methods

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 [START_REF] Cable | Geomorphological and cryostratigraphical analyses of the Zackenberg Valley, NE Greenland and significance of Holocene alluvial fans[END_REF]. 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 [START_REF] Elberling | Soil and Plant Community-Characteristics and Dynamics at Zackenberg[END_REF]. 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 [START_REF] Elberling | Soil and Plant Community-Characteristics and Dynamics at Zackenberg[END_REF][START_REF] Stewart | Forecasted homogenization of high Arctic vegetation communities under climate change[END_REF]. 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 [START_REF] Pedersen | Spatiotemporal characteristics of seasonal snow cover in northeast Greenland from in situ observations[END_REF]. 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 -2 since 1996 [START_REF] Schmidt | Long-term patterns of muskox (Ovibos moschatus) demographics in high arctic Greenland[END_REF]. 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 [START_REF] Schmidt | Long-term patterns of muskox (Ovibos moschatus) demographics in high arctic Greenland[END_REF]. 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 [START_REF] Stewart | Forecasted homogenization of high Arctic vegetation communities under climate change[END_REF]. 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 by mimicking muskox grazing i.e., grabbing and pulling a handful of plant material. Plant material (min. ~2g wet weight) was subsequently transferred into paper bags. Soil was collected in each plot by mixing and homogenizing three sub-samples from the A-horizon (i.e., the upper 0-10 cm layer) using nitrile gloves and a pre-cleaned stainless-steel garden shovel. Soil material (~10g 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 plot-based 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 35°C for 5 days under constant air flow. While we acknowledge that oven drying is not the ideal method, doing so was necessary due to logistical constraints in the field. However, oven-drying soil samples at low temperatures limits the loss of volatile elements [START_REF] Roos-Barraclough | An analytical protocol for the determination of total mercury concentrations in solid peat samples[END_REF]. Once dried, 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 H2O (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 30Hz for 120s. All material was cleaned between each sample, i.e., rinsed under tap water then rinsed with mQ H2O and finally disinfected 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% HNO3 and mQ H2O) 50mL digitubes to which 3mL of ultrapure HNO3 (Fisher Chemical, Nitric Acid HNO3, 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 H2O 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 H2O 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 -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 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 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℃ and a reduction furnace operating at 700℃. 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 values of the reference material had a mean and standard deviation of N%: 2.965 ±0.054 and C% 46.01 ± 1.52 and were measured at N%: 2.939± 0.080 and C% 46.67 ±1.89 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 [START_REF] Leroux | Stoichiometric distribution models: ecological stoichiometry at the landscape extent[END_REF]. 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., [START_REF] Sach | Spatial geochemistry influences the home range of elephants[END_REF] 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 logtransformed 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 - 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: [START_REF] Stewart | Forecasted homogenization of high Arctic vegetation communities under climate change[END_REF] 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 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 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 andS4). 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 [START_REF] Boyce | Relating populations to habitats using resource selection functions[END_REF]. 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 [START_REF] Manly | Resource selection by animals: Statistical analysis and design for field studies[END_REF]. 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 [START_REF] Schmidt | Long-term patterns of muskox (Ovibos moschatus) demographics in high arctic Greenland[END_REF]. 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 [START_REF] Northrup | Practical guidance on characterizing availability in resource selection functions under a use-availability design[END_REF].

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 (Spearman rank correlations (rs) 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 [START_REF] Tomassini | Density, snow and seasonality lead to variation in muskox habitat selection during summer[END_REF]. 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 k-fold crossvalidation procedure as proposed by [START_REF] Boyce | Evaluating resource selection functions[END_REF] and [START_REF] Bastille-Rousseau | Leveraging multidimensional heterogeneity in resource selection to define movement tactics of animals[END_REF]. For this, we calculated cross-validated Spearman rank correlations (rs) between ten RSF-bin ranks and five test-training sets. We repeated this procedure 10 times to determine whether the rs 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 [START_REF] Schmidt | Long-term patterns of muskox (Ovibos moschatus) demographics in high arctic Greenland[END_REF]]. 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 [START_REF] Reynolds | Dynamics and range expansion of a re-established muskox population[END_REF][START_REF] Schmidt | Long-term patterns of muskox (Ovibos moschatus) demographics in high arctic Greenland[END_REF].

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; [START_REF] Schmidt | Long-term patterns of muskox (Ovibos moschatus) demographics in high arctic Greenland[END_REF]. Muskox cows give birth to a calf in late winter if sufficient energetic reserves are acquired the preceding summer [START_REF] Desforges | Quantification of the full lifecycle bioenergetics of a large mammal in the high Arctic[END_REF]. 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 [START_REF] Jolliffe | Principal Component Analysis for Special Types of Data[END_REF], 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).

Results

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 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 rs 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 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. 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 rs of 0.71 (sd=0.16), which was significantly higher (p<0.001) than for a randomized model with a mean rs 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 [START_REF] Flueck | The implications of selenium deficiency for wild herbivore conservation: a review[END_REF][START_REF] Flynn | Indications of Copper Deficiency in a Subpopulation of Alaskan Moose[END_REF][START_REF] Robbins | Wildlife feeding and nutrition[END_REF]. 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 [START_REF] Flueck | The implications of selenium deficiency for wild herbivore conservation: a review[END_REF].

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 [START_REF] Grace | Trace element metabolism, dietary requirements, diagnosis and prevention of deficiencies in deer[END_REF]. 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 [START_REF] Mosbacher | Hair mineral levels as indicator of wildlife demographics?-a pilot study of muskoxen[END_REF]. 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 [START_REF] Anke | The biological and toxicological importance of molybdenum in the environment and in the nutrition of plants, animals and man: Part V: Essentiality and toxicity of molybdenum[END_REF].

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.

Moreover, the bioavailability of Mo is unlikely to exceed toxicity levels of muskoxen in this system, as has been reported for other herbivores [START_REF] Reis | Mineral element and heavy metal poisoning in animals[END_REF].

Arctic and boreal ecosystems are characterized by strong N limitation and high C stocks in soils [START_REF] Liu | Nitrate is an important nitrogen source for Arctic tundra plants[END_REF]. As a consequence, most ecological stoichiometry studies focus on C:N ratios to assess vegetation quality and dynamics at northern latitudes [START_REF] Mosbacher | Muskoxen modify plant abundance, phenology, and nitrogen dynamics in a high Arctic fen[END_REF][START_REF] Peng | The relationship between relative growth rate and whole-plant C : N : P stoichiometry in plant seedlings grown under nutrient-enriched conditions[END_REF] or to estimate linkages with consumer distribution [START_REF] Leroux | Stoichiometric distribution models: ecological stoichiometry at the landscape extent[END_REF][START_REF] Rizzuto | Forage stoichiometry predicts the home range size of a small terrestrial herbivore[END_REF]. 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 [START_REF] Barboza | The nitrogen window for arctic herbivores: plant phenology and protein gain of migratory caribou (Rangifer tarandus )[END_REF]. 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 [START_REF] Desforges | Quantification of the full lifecycle bioenergetics of a large mammal in the high Arctic[END_REF][START_REF] Thing | Ecology of Muskoxen in Jameson Land, Northeast Greenland[END_REF].

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 [START_REF] Anke | The biological and toxicological importance of molybdenum in the environment and in the nutrition of plants, animals and man: Part V: Essentiality and toxicity of molybdenum[END_REF][START_REF] Flueck | The implications of selenium deficiency for wild herbivore conservation: a review[END_REF].

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 [START_REF] Massányi | Effects of Cadmium, Lead, and Mercury on the Structure and Function of Reproductive Organs[END_REF]. 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 [START_REF] Kabata-Pendias | Trace Elements in Abiotic and Biotic Environments[END_REF]. 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 [START_REF] Tomassini | Density, snow and seasonality lead to variation in muskox habitat selection during summer[END_REF].

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

Caveats and future prospects

To link long-term 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 [START_REF] Elmendorf | Global assessment of experimental climate warming on tundra vegetation: Heterogeneity over space and time[END_REF], 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 [START_REF] Abraham | Understanding anthropogenic impacts on zoogeochemistry is essential for ecological restoration[END_REF][START_REF] Schmitz | Animals and the zoogeochemistry of the carbon cycle[END_REF]. 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.

Our study relied on relatively small-scale (52m 2 ), yet long-term, survey data of muskox observations and demography. However, muskoxen typically roam over larger areas [START_REF] Schmidt | Ungulate movement in an extreme seasonal environment: year-round movement patterns of high-arctic muskoxen[END_REF] with substantial individual variation in habitat use [START_REF] Beumer | Accounting for behaviour in fine-scale habitat selection: A case study highlighting methodological intricacies[END_REF][START_REF] Van Beest | Environmental conditions alter behavioural organization and rhythmicity of a large Arctic ruminant across the annual cycle[END_REF]. 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 [START_REF] Nunes | Ecological trap for seabirds due to the contamination caused by the Fundão dam collapse, Brazil[END_REF]. 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 [START_REF] Desforges | Quantification of the full lifecycle bioenergetics of a large mammal in the high Arctic[END_REF][START_REF] Sperfeld | Bridging Ecological Stoichiometry and Nutritional Geometry with homeostasis concepts and integrative models of organism nutrition[END_REF].

Free-ranging animals are confronted by a myriad of environmental conditions that they need to balance to optimize fitness [START_REF] Gaillard | Habitat-performance relationships: finding the right metric at a given spatial scale[END_REF][START_REF] Stein | Environmental heterogeneity as a universal driver of species richness across taxa, biomes and spatial scales[END_REF][START_REF] Van Beest | Behavioural responses to thermal conditions affect seasonal mass change in a heat-sensitive northern ungulate[END_REF]. In our analyses, we focused on vegetation types, snow cover, and multiple essential and non-essential 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 cover, plant productivity and access to forage biomass during spring and summer [START_REF] Beumer | An application of upscaled optimal foraging theory using hidden Markov modelling: year-round behavioural variation in a large arctic herbivore[END_REF][START_REF] Tomassini | Density, snow and seasonality lead to variation in muskox habitat selection during summer[END_REF]. Indeed, muskoxen use this period primarily to replenish energetic reserves to survive the upcoming winter and to be able to develop a fetus [START_REF] Desforges | Environment and physiology shape Arctic ungulate population dynamics[END_REF]. 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 [START_REF] Rantanen | The Arctic has warmed nearly four times faster than the globe since 1979[END_REF]. 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 [START_REF] Miner | Emergent biogeochemical risks from Arctic permafrost degradation[END_REF]. While the potential risks of leaching of trace metals following permafrost degradation have been studied to some extent [START_REF] Antcibor | Trace metal distribution in pristine permafrost-affected soils of the Lena River delta and its hinterland, northern Siberia, Russia[END_REF], 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 climateforecast 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 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 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 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.

Fig. 1 :

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

Fig. 2 :

 2 Fig. 2: Boxplots showing measurements of chemical elements in soil (yellow) and vegetation

Fig. 3 :

 3 Fig. 3: Predicted spatial distribution of chemical elements in vegetation across the Zackenberg valley in northeast Greenland. The unit for C, N

Fig. 4 :

 4 Fig. 4: Biplot showing how spring snow cover and the relative use of chemical elements in

Fig. 5 :

 5 Fig. 5: Biplot showing how spring snow cover and the relative use of vegetation types during

Fig. 6 :

 6 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

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

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.

Data availability

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

Ethics approval

Not applicable.

Duquenoy and Aurélie Marquet (both at the ICP-MS platform at the Observatoirie Midi-Pyrenees, Toulouse, France) for assistance with the ICP-MS analyses. We also thank three anonymous reviewers for their extremely constructive comments that helped to improve this article.

Graphical Abstract