Spatial models of carbon, nitrogen and sulphur stable isotope distributions (isoscapes) across a shelf sea: An INLA approach

Spatial models of variation in the isotopic composition of structural nutrients across habitats (isoscapes) offer information on physical, biogeochemical and anthropogenic processes occurring across space, and provide a tool for retrospective assignment of animals or animal products to their foraging area or geographic origin. The isotopic differences among reference samples used to construct isoscapes may vary spatially and according to non‐spatial terms (e.g. sampling date, or among individual or species effects). Partitioning variance between spatially dependent and spatially independent terms is a critical but overlooked aspect of isoscape creation with important consequences for the design of studies collecting reference data for isoscape creation and the accuracy and precision of isoscape models. We introduce the use of integrated nested Laplace approximation (INLA) to construct isoscape models. Integrated nested Laplace approximation provides a computationally efficient framework to construct spatial models of isotopic variability explicitly addressing additional variation introduced by including multiple reference species (or other recognized sources of variance). We present carbon, nitrogen and sulphur isoscape models extending over c. 1 million km2 of the UK shelf seas. Models were built using seven different species of jellyfish as spatial reference data and a suite of environmental correlates. Compared to alternative isoscape prediction methods, INLA‐spatial isotope models show high spatial precision and reduced variance. We briefly discuss the likely biogeochemical explanations for the observed spatial isotope distributions. We show for the first time that sulphur isotopes display systematic spatial variation across open marine shelf seas and may therefore be a useful additional tool for marine spatial ecology. The INLA technique provides a promising tool for generating isoscape models and associated uncertainty surfaces where reference data are accompanied by multiple, quantifiable sources of uncertainty.


| INTRODUC TI ON
The isotopic ratio of elements such as hydrogen, oxygen, carbon, nitrogen and sulphur varies systematically across the natural environment. Modelling these spatial differences through mechanistic or statistical models (isoscapes [West, Bowen, Dawson, & Tu, 2010]) offers insight into the biogeochemical processes leading to spatial variation in isotopic expressions, and provides a tool for the retrospective assignment of animals or animal products back to their origin or foraging area (Hobson, Wassenaar, & Taylor, 1999), with broad uses in animal and human migration and tracking (Ehleringer et al., 2008;Hobson, 1999;Hobson, Barnett-Johnson, & Cerling, 2010), trophic ecology (Jennings & van der Molen, 2015;Olson et al., 2010) and traceability within consumer goods supply chains (Chesson, Valenzuela, O'Grady, Cerling, & Ehleringer, 2010;Kelly, Heaton, & Hoogewerff, 2005). Isoscapes have been used extensively in terrestrial ecological and forensic applications, particularly isoscapes describing spatial variations in hydrogen and oxygen isotopes of precipitation (Bowen, 2010). Spatial variation in isotopic compositions in marine systems has also been explored (Cherel & Hobson, 2007;Schell, Saupe, & Haubenstock, 1989). However, relatively few continuous surface isoscapes have been published in marine compared to terrestrial systems, probably due to the difficulty in obtaining sufficient reference samples over appropriate spatial and temporal scales.
To construct a continuous surface isoscape model, isotopic compositions of reference materials or organisms are typically projected across space using either spatial interpolation methods (Trueman, MacKenzie, & St John Glew, 2017;Vander Zanden et al., 2015), by statistical inference based on correspondence between measured data and environmental correlates (Bowen & Wilkinson, 2002;Courtiol & Rousset, 2017;Jennings & Warr, 2003) or a combination of both (MacKenzie, Longmore, Preece, Lucas, & Trueman, 2014;Wunder, 2010). Ideally, the isotopic composition of reference samples should vary only according to spatially dependent effects.
However, additional isotopic variance among reference samples is commonly introduced through processes such as collection of samples over different time-scales or differences in ecology or physiology among individuals or species. As the spatial scale of a study area increases, the ease of collecting uniform reference samples generally decreases, especially where isoscape reference data are compiled from opportunistically collected samples. Accounting for and quantifying spatially dependent and spatially independent variance is a key component of isoscape model creation. In simple interpolation models, the variance associated with the prediction of expected isotopic compositions at any point in space increases with distance from discrete sampling points; therefore, irregular spacing of reference samples produces spatial gradients in isoscape uncertainty, which may bias interpretations. When environmental correlates are introduced, estimating variance becomes more complex due to error associated with the relationship between measured data and environmental correlates, which is itself spatially varying, but rarely quantified across space (Bowen & Revenaugh, 2003;Courtiol & Rousset, 2017). To date, many isoscape models either assume spatially invariant uncertainty in the relationship between measured data and environmental correlates (Jennings & Warr, 2003), infer spatial variance by interpolating residuals from regression models (MacKenzie et al., 2014) or draw on resampling methods to estimate spatially varying uncertainty (Wunder & Norris, 2008). Courtiol and Rousset (2017) introduced a frequentist mixed modelling approach that enables spatially explicit variance surfaces to be calculated by including location as a random effect but at the cost of slow computational processing. Here, we introduce an alternative approach to isoscape generation based on integrated nested Laplace approximations (INLAs). We aim to address the common issue of limited sample availability by modelling isoscapes using multiple species and explicitly addressing spatial isotopic variation due to mixed sample sources. Many commonly used isoscape prediction methods are unable to incorporate multiple sample sources while quantifying associated spatial variance and including boundary effects (Table 1). We explore the use of recently developed Bayesian hierarchical modelling techniques using INLA. We firstly produce isoscape models for a restricted region, the North Sea, using a single jellyfish species as a spatial reference dataset, and compare the assignment accuracy and precision associated with INLA-produced and alternative North Sea isoscape models . Secondly, we predict isoscapes for carbon, nitrogen and sulphur across the wider UK shelf sea area using multiple reference jellyfish species.
The North Sea and wider UK shelf seas host some of the most globally productive fisheries, regionally significant oil, gas and renewable energy resources and infrastructure and intensive shipping activity. The UK shelf region has received extensive detailed investigation into spatial isotopic variability, with carbon and nitrogen isoscape models previously produced using purpose collected baseline samples, rather than commonly adopted opportunistic sampling. Barnes, Jennings, and Barry (2009b) and Jennings and Warr (2003) and Jennings and van der Molen (2015) used queen scallops Aequipecten opercularis from known catch locations as reference samples, coupled with environmental variables; however, variance surfaces were only calculated by Jennings and van der Molen (2015).
High resolution, in situ sample-based isoscapes have been modelled for the North Sea using lion's mane jellyfish Cyanea capillata as reference organisms through ordinary kriging of evenly spaced samples  and with additional environmental variables (MacKenzie et al., 2014). Spatially explicit variance surfaces were calculated in both examples and initial assignments of invertebrate, fish and seabird samples have proven successful (St John Glew et al., 2018;Trueman et al., 2017). However, this approach is constrained by the availability and distribution of a single reference species across the region of interest, limiting marine isoscape modelling capabilities across larger spatial scales, as no single jellyfish species is distributed across the entire range of the UK shelf seas. In addition, barrier effects (e.g. uneven coastlines) are particularly important in basin scale marine isoscape predictions, yet many existing modelling techniques do not enable easy incorporation of coastlines and boundaries (Table 1).

| Data collection and stable isotope analysis
To construct isoscape models of UK shelf seas, we collected 627 Jellyfish bell tissue δ 13 C values showed a significant negative linear relationship with C:N ratios (p < 0.005, slope = −2.22, adjusted R 2 = 0.06).
To correct for lipid-related variance in δ 13 C values, we applied an algebraic correction (Kiljunen et al. (2006)). Lipid-corrected carbon and nitrogen isotopic data from Queen scallops of known location were taken from Jennings and Warr (2003) and Barnes et al. (2009b), for scallops collected between 25 July and 29 September 2001 and from Barnes, Jennings, and Barry (2009a) for scallops collected in 2010.
We estimated within-species variation in jellyfish stable isotope compositions by averaging (mean) the among-individual standard deviation of the same species occurring at the same sampling location. We calculated among-species average isotopic differences by calculating the mean isotopic difference between species at the same location and then averaging across all locations.  Simple interpolation/kriging of sample data

| Environmental data
Linear regression models of sample data (with species as a random effect) and additional environmental variables, followed by interpolation General additive models of sample data (with species as a random effect) and additional environmental variables, followed by interpolation

| Model formation
We describe a Bayesian hierarchical spatial modelling framework for the shelf seas, predicting δ 13 C, δ 15 N and δ 34 S values using INLA via the r-inla package (http://www.r-inla.org) (Rue, Martino, & Chopin, 2009). This approach differs from the frequentist mixed modelling approach introduced by Courtiol and Rousset (2017) by adopting a Bayesian framework enabling uncertainty to be more easily interpretable, allowing the inclusion of boundary effects and solving the spatial dependency term in an alternative and faster way.
When modelling across a spatial range, ordinary linear regression ignores spatial dependency between sampling locations. Through the latent Gaussian field with Matérn correlation, r-inla provides a means to explicitly incorporate spatial dependency: where y (si) are the response values at all sampling locations which are assumed to be normally distributed with mean x (si) and variance σ 2 . u (si) is the spatial dependency random effect. r-inla includes a stochastic partial differential equation ( species were sampled at the same location. Put simply, we modelled each isotope value as a function of a set of covariates X i where species and the underlying spatial structure were included as random effects. Models were specified as: where Y i is the isotope value (δ 13 C, δ 15 N or δ 34 S) at location i, X i is a vector containing the environmental covariates as linear fixed effects, β i is a vector of parameters to be estimated, U i is the species random effect with assumed Gaussian distribution, W i represents the smooth spatial effect, linking each observation with a spatial location, with the elements of Ω estimated using Sea, we also compared models with and without a spatial effect.
When similar DIC values were observed (within 2), the simpler model was selected (Burnham & Anderson, 2003). We only compared spatial models within the shelf seas to capture and model isotopic spatial variation driven both by the larger spatial extent of the UK shelf study area, and additional variance introduced by multispecies reference samples. We used non-informative default priors for each model. We used the best model for each isotope to predict isotopic com-

| Comparing INLA and kriging isoscape models for single species isoscapes
To assess the differences between the traditional ordinary krig- where the true location falls within the assigned area Vander Zanden et al., 2015). We compared accuracy and precision of known-origin scallops assignments sampled in 2001 and 2010 between INLA-predicted North Sea isoscapes and ordinary kriging isoscapes produced by Trueman et al. (2017). We performed all analyses using r 3.4.2 (R Core Development Team, 2016).

| Within-and between-species variability
The variations in average stable isotope ratios within species sampled at the same locations were relatively consistent across species in both carbon (0.37-0.62‰) and sulphur (0.41-0.63‰) apart from crystal jellyfish where among-individual variation was higher for both δ 13 C (1.08‰) and δ 34 S (0.75‰) ( Table 2). Within-species differences in nitrogen were more variable ranging from 0.44‰ in blue jellyfish to 1.65‰ in crystal jellyfish (Table 2). Among-species differences ranged considerably between species and isotope ( Table 2).
δ 34 S differences were relatively constrained with differences ranging from 0.01 to 1.45‰ whereas δ 13 C differences ranged from 0.03‰ (between mauve stingers and compass jellyfish) to 3‰ (between barrel and compass jellyfish). Among-species differences in δ 15 N varied over the largest range, from 0.02‰ between lion's mane and blue jellyfish up to 7.1‰ between barrel and mauve jellyfish.

| The North Sea isoscape models
The best-fit carbon and nitrogen isoscape models for the North Sea were non-spatial models excluding interaction terms (Table 3): In both cases, interaction and spatial dependency terms did not improve model fit; therefore, the simplest models were selected (Table 3).
Broad spatial patterns in δ 13 C and δ 15 N ranges (Figure 4a

| UK shelf sea isoscape models
Global models, including first-order interaction terms, were the bestfit for carbon, nitrogen and sulphur isoscapes; Best-fit models for carbon and sulphur UK shelf sea isoscapes had moderate fit (R = 0.47, p < 0.05 and R = 0.50, p < 0.05 respectively) ( Table 3). The best fitting nitrogen isoscape model had a stronger fit (R = 0.80, p < 0.05) ( Table 3).
Minimal residual isotopic variability between species remained with the chosen carbon and sulphur isoscape prediction models, indicating that the majority of species isotopic variability was able to be explained by the combination of covariates, and interactions between these covariates, included within the models (Figure 7a,c).
Residual nitrogen isotopic variability has a larger range between species (c. 6‰), particularly between mauve stinger jellyfish depleted  TA B L E 2 Within (red)-and between (row 1 -column 1)-species isotopic differences (black). Calculated at locations where multiple individuals of the same species or multiple species occur and averaged across locations. Within-species isotopic difference is the among-individual standard deviation of the same species occurring at the same sampling locations and averaged across all locations. Between-species isotopic difference is the difference between different species sampled at the same

| D ISCUSS I ON
This study has two main aims to introduce INLA as a powerful tool for creating isoscape models incorporating environmental correlates as predictors and where reference samples contain a source of variance that is not spatially dependent and to describe the spatial variation in δ 13 C, δ 15 N and δ 34 S values across the shelf seas of the British Isles. We have shown that INLA-generated isoscape models have comparable accuracy and precision to simple kriging where reference samples are evenly distributed and common form (in our case the same species). We then extended the approach to draw isoscape models and uncertainty surfaces across a wide shelf sea area where collection of reference samples from a single species would be impossible.

| INLA as a tool for creating isoscape models
Creating isoscape models with associated uncertainty surfaces in regions where reference samples are either irregularly spaced and/or contain additional sources of isotopic variability is challenging (Courtiol & Rousset, 2017). Spatial modelling using the INLA approach addresses many common constraints. By using the INLA approach, we are able to incorporate environmental data and a species random effect into our isoscape prediction models. Although this can also be achieved by using a mixed effects model approach (Courtiol & Rousset, 2017), the INLA approach is unique in that it provides a computationally rapid technique to quantify the spatial variance due to the species random effect, which is essential for accurate measures of variance and subsequent isoscape assignments. In addition, INLA enables the incorporation of boundary effects, to model and predict around physical barriers, which is particularly useful in marine environments.
TA B L E 3 Model fit results for the carbon and nitrogen North Sea, and carbon, nitrogen, and sulphur UK shelf sea isoscape prediction models. Global models excluding and including first-order interaction terms were tested. Inclusion of the spatial term was also tested within the North Sea models. Model fit was tested using the deviance information criteria (DIC), and assessing the Pearson correlation between observed and fitted values. The t, R, 95% confidence intervals around R, and the degrees of freedom (df) are reported. The models displayed in red were the chosen models for isoscape predictions The INLA-predicted North Sea δ 13 C and δ 15 N isoscapes ( Figure 4) are broadly similar to isoscapes produced from ordinary kriging of identical lion's mane jellyfish data with similar low variance estimates within the spatial confines of the reference sample (Figure 5c,d). Accordingly, the accuracy and precision by which scallops of known origin could be assigned back to origin were also comparable to that demonstrated by Trueman et al.
(2017) ( Figure 6). However, it must be noted that accuracy and precision results cannot be extrapolated outside the North Sea range and do not reflect wider shelf sea isoscape accuracy and precision.
Given that the INLA approach draws on environmental correlates to predict isotopic compositions, one might expect simple kriging to produce more accurate isoscape models where reference sample collection is evenly spaced and dense compared to the spatial scale of isotopic gradients. The similarity in uncertainty between the two methods found here reflects the relatively strong statistical relationships between environmental correlates and reference isotope data.
Where evenly spaced reference samples cannot be recovered across the entire region of interest, or where spatial variation in isotope values is expected to occur at smaller spatial scales than the spacing between reference samples, isoscapes drawn from environmental predictors may produce more accurate and precise assignments.
One significant benefit of the spatial INLA approach is the ability to account for sources of isotopic variance in the reference data other F I G U R E 6 Accuracy (the proportion correctly assigned) and precision (proportion of the total surface area) of assignment to both the original North Sea kriging isoscape models  shown in black, and the new integrated nested Laplace approximation (INLA) modelled North Sea isoscapes shown in blue for the 2001 and 2010 scallop datasets. The red line represents the accuracy and precision values if assignments were no better than random F I G U R E 7 Marginal posterior distributions of the species random effect for the chosen carbon, nitrogen and sulphur isoscape prediction models. π is the species-level deviation from the overall mean isotope value, and D is the data. Distributions represent the probability density of a given isotopic difference, given the data and represents species differences that remain after the models have been applied. Differences between species represent isotopic differences unable to be explained by environmental variables than spatially varying terms. In our case, INLA allowed us to identify and account for large, among species isotopic differences, ranging between 0.03-3.0‰ in δ 13 C, 0.02-7.1‰ in δ 15 N and 0.01-1.45‰ in δ 34 S (Table 2) into the spatial model. The 'species effects' are quantified as residual differences unaccounted for by the environmental predictors within the final models and displayed as marginal distributions in Figure 7. Both carbon and sulphur 'best-fit' models were able to explain all isotopic differences between species, whereas residual nitrogen isotopic differences were still observed. Mauve stinger and crystal jellyfish had markedly different δ 15 N values with mauve stingers displaying consistently low and crystal jellyfish consistently high δ 15 N values. Isotopic variation among different species is expected, likely due to different diets, habitat uses and metabolic processes. Deciphering the reasons behind these species isotopic differences is beyond the scope of this study, but we emphasize the importance of treating gelatinous zooplankton as separate species in any isotopic study.
In this example, we use INLA to incorporate isotopic differences between species; however, the same concept applies whenever data with known, or assumed, differences must be combined. For example, in isoscape models where plankton or zooplankton are sampled and grouped (McMahon, Hamady, & Thorrold, 2013;Schell, Barnett, & Vinette, 1998); where data have been collected from multiple sources (Bataille et al., 2018); or where different sampling techniques have been adopted. The same approach could also be used to incorporate temporal variability in sample collection (Bowen & Revenaugh, 2003;Flockhart et al., 2013). While samples in the current study were collected over 2 years, sampling locations did not overlap across different times, so temporal effects could not be explicitly quantified.
F I G U R E 8 UK shelf sea carbon, nitrogen and sulphur isoscape models (a, c, e) and associated variance of the posterior predicted distribution, after species random effects have been accounted for (b, d, f

| Isotopic variability across the UK shelf seas
Stratification and mixing extent are strong drivers of spatial isotopic variability, with front locations closely matching isotope ratio boundaries in carbon, nitrogen and sulphur (Miller & Christodoulou, 2014).
Isotopic ratios are also strongly influenced by freshwater and terrestrial inputs. Freshwater has a lower δ 34 S ratio compared to seawater (Fry, 2002), reducing δ 34 S values in regions with high freshwater input (e.g. eastern Irish Sea, southern North Sea and English Channel and areas off West Scotland) (Painting et al., 2013) ( Figure 8). Anthropogenic nutrient sources enter the marine environment through estuaries (Howarth, 1998) and influence productivity causing increased δ 13 C and δ 15 N values in coastal and estuarine environments (e.g. southern North Sea, eastern English Channel, eastern Irish Sea) (Painting et al., 2013).
Production source also influences isotopic variability, with phytoplankton community structure differing between the northern and southern North Sea (Ford et al., 2016), the presence of cyanobacteria within the western English Channel (Rees, Gilbert, & Kelly-Gerreyn, 2009) decreasing δ 15 N but increasing δ 13 C values due to nitrogen fixation (Levitan et al., 2007) and influence of microalgae increasing δ 13 C and δ 15 N values around the East Anglian coast and into the southern North Sea (Bristow et al., 2013) (Figure 8). Variance surfaces are similar for each isoscape, with uniform variance across the majority of the UK shelf, but greater values found within more dynamic regions such as the eastern English Channel and eastern Irish Sea.

| CON CLUS ION
The principle reason for adopting an INLA (or mixed model) approach to generate an isoscape is to account for variance in reference samples that is not explicitly spatial in origin. Where reference datasets can be assembled from the same species, collected at the same time and processed in the same way, simpler spatial modelling or kriging interpolation approaches may be favoured; however, in many cases, some extra non-spatially dependent variance terms are introduced because of the difficulty of obtaining uniform reference samples. In our example, jellyfish species have varying distributions across the spatial range, so to generate a single isoscape model required use of multiple species and therefore the introduction of random effect of species. The INLA approach is a promising method for accounting for additional non-spatially dependent isotopic variance within reference samples. Although our study focuses on marine carbon and nitrogen and the newly introduced sulphur isotopes, the same methods and benefits and limitations are applicable across all environments and isotope systems.