To extract time-resolved otolith oxygen stable isotope values (δ¹⁸O) for use as a proxy for seawater temperature, otoliths were prepared for SIMS analyses as reported in Hanson et al. (2013). Oxygen isotope ratios were measured with a Cameca IMS-1270 ion microprobe (#309) (Edinburgh Ion Microprobe Facility) using a ~ 5 nA primary 133Cs⁺ beam focused to 15 μm diameter. Step sizes between sample spots averaged 63 μm and spot depth was 2 - 3 μm. Secondary ions were extracted at -10 kV, and ¹⁶ O and ¹⁸ O were monitored simultaneously on dual Faraday cups (L’2 and H’2) at intensities of 3 x 10⁹ and 4 x 10⁶ cps, respectively. Mass resolution was 2600 for oxygen to remove ¹⁶ O ² H and ¹⁷ O ¹ H on ¹⁸ O. Each analysis involved a pre-sputtering period of 40 s, followed by automatic secondary beam and entrance slit centering, and finally data collection in two blocks of ten cycles. Mean external precision based on daily analysis of an inorganic carbonate standard (University of Wisconsin Calcite, UWC) was ± 0.19SD. Otolith oxygen isotope values are reported in per mil (‰) notation relative to Vienna Pee Dee Belemnite (VPDB) (Coplen, 1996), according to the equation δ ¹⁸ O_sample = [(R_sample /R_standard) – 1] × 1000, where R is equal to ​ 18 O ​ 16 O . Each sampling transect originated at the outer margin of the otolith and was terminated once a freshwater isotope signal was recorded (when values dropped to< -4 ‰). To compare isotope profiles between individual salmon, and to estimate a chronology for the sample spots, the transect was re-scaled to relative distance from the transition between freshwater and marine residency of the juvenile fish (i.e. smolt emigration from the river).

The isotope fractionation equation, 1000lnα = a + b(10³ T ^-1) (Kim et al., 2007) was used to estimate seawater temperature values at each SIMS sample spot, where α = δ 18 O o t o l i t h + 1000 δ 18 O w a t e r + 1000 and temperature T is in Kelvin. For the present analysis, the value of a was estimated from similar data from 1SW Atlantic salmon (-33.43, Hanson et al., 2013), and assuming a slope b = 17.88 ± 0.13 derived using synthetic aragonite (Kim et al., 2007). δ ¹⁸ O_water was set to 0, which is a reasonable average across the North Atlantic (LeGrande and Schmidt, 2006). VPDB values were converted to the Vienna standard mean ocean water (VSMOW) scale using the Coplen et al. (1983) equation for temperature calculations. Estimated temperatures were converted to degrees Celsius (°C) by subtracting 273.15.

Detailed investigation of scale growth was necessary for the present analysis for two reasons; first, in the absence of an otolith growth model for adult Atlantic salmon, scale growth was used as a proxy with which to anchor δ ¹⁸ O values in time (Hanson et al., 2013). As salmon grow their scales increase in size, but not number, and growth of the scale radius is directly linked to growth in body length of the fish. It is unlikely that the relative distances along the otolith and scale growth axes would be exactly equivalent but we considered this assumption to be justified given that growth of both structures is linked to somatic growth (Degens et al., 1969; Campana and Thorrold, 2001) and is correlated in this species (Wells et al., 2003). Our second objective here was to include, as observational data in the movement model, the body length of the individual (and thence an estimate of swim speed) at each δ ¹⁸ O spot back-calculated using the incremental measurements of the scale circuli (see below).

Details of Atlantic salmon scale increment analysis, and the method for assigning calendar dates to each growth ring for an individual fish, are reported by Todd et al. (2013). This method requires information on three recognizable scale features. Emigration of the juvenile smolt from freshwater was assumed to be April 27 2007, or May 7 2008 based on median smolt emigration dates for those calendar years taken from the only available Scottish monitoring archive (for River North Esk, E Scotland; Marine Scotland Science, unpublished data). The midpoint growth circulus of the winter annulus (1SW adult), or both winter annuli (2SW adult) was assumed to be set at the winter solstice (December 21), and the final capture date (= scale edge) was known. This method allowed the date of each growth ring to be expressed as the relative distance from smolt emigration.

We back-calculated body length increment throughout the marine phase from scale circuli measurements. Smolt length at emigration was back-calculated using a modification of the Dahl-Lea method (see Hanson et al., 2019 for details), and all subsequent growth in body length during the marine phase was assuming to be isometric with scale radius increment. Finally, back-calculated body lengths at each observed δ ¹⁸ O otolith value were assigned by matching to the nearest calendar date.

Given the series of otolith-derived water temperature time series (T^O ), we used a state space model to estimate progressive changes in the geographic location of the sampled salmon throughout their marine phase. Our modelling of Atlantic salmon movement assumed otolith T^O are predominantly reflective of ocean surface temperature which is justified by archival DST data having confirmed that these fish forage predominantly within the top 10m of the water column, whilst also undertaking occasional intermittent deep dives (Guðjónsson et al., 2015; Hedger et al., 2017).

We adapted the methodology of Thygesen et al. (2009) who present an HMM for the geolocation of fish, in discrete space, using DST data. In this approach, the true location of the fish at a given time-step was modelled as an unobserved state in a Markov process. The state transition probabilities were determined by a movement model, with the DST data providing observations on the hidden states.

We used a 0.25° regular grid between -75.125 and 25.125^o longitude and 45.125 and 80.125^o latitude to define the state space, Z, of the HMM. The resolution of the state space was aligned to available SST data (NOAA OISST v2; https://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.highres.html ), with an associated land-sea mask. States that were designated as land in the NOAA SST data had a likelihood of zero for all time-steps. Transitions between states occurred at weekly time intervals (modelled as the first day in each 7-day period since the date of river emigration).

Following Thygesen et al. (2009), let p(τ|τ) denote the discrete probability distribution of the location of the fish, i.e. over Z, given all observations taken at, or before, time τ and p(τ + 1|τ) be the predicted distribution at time τ + 1 given all observations strictly before time τ + 1. Applying the movement model then corresponds to the update:

where H defines the transition probabilities between states as determined by the movement model, and ∗ denotes convolution over the state space (i.e. all locations), such that for state s ∈ Z we have:

where H_s,t is the probability of moving to state s from state t and p _t is the value of p(τ|τ) at t.

Previous HMMs of fish behavior have modelled movement as a diffusion process resulting from an assumed unbiased random walk at a finer temporal resolution than that of the HMM (Thygesen et al., 2009; Guðjónsson et al., 2015; Strøm et al., 2017; Strøm et al., 2018). Unbiased random walk models do not appropriately capture persistence (i.e. directional correlation) in animal movement and therefore may not always be suitable for application to migratory data such as those of the present study (Codling et al., 2008); this is particularly relevant when the data are sparse and have high uncertainty. Random walk models, usually modelled with a Gaussian distribution, also have the disadvantage of allowing for a non-zero probability at large distances in finite time. As a simple alternative we used a uniform circular movement model, with radius dependent on the swim speed estimated from back-calculated body length of the fish in meters, l (Andersen et al., 2007). Thus the probability of moving to location s given starting location t and body length l_τ at time-step τ is

where Δx_s,t , Δy_s,t are the x and y components of the horizontal distance in meters between the two locations s and t, calculated as:

where λ_z , ϕ_z are the longitude and latitude of state z in radians, and an equirectangular projection is used to approximate distance. Given that our state space was defined over a regular longitudinal and latitudinal grid, the movement model results in an ellipse with its extent dependent on the latitude of the location at which it was applied. V is the progression rate in a time interval of one week,

assuming a migration speed, b, of 1.25 body lengths per second (BLs^-1; Drenner et al., 2012; Thorstad et al., 2012). As a consequence of using a body length-dependent progression rate, the number of non-zero values in H increases as l increases over time.

Next, the predicted distribution at time τ + 1, p(τ + 1|τ) is updated to the estimate p(τ + 1|τ + 1) by taking into account the data at time τ + 1. For the present application, in the absence of DST data, we treat the series of otolith-derived SST values (T^O ) as observations of the sea surface temperature (T^SST ) at the unobserved location.

Given an observation of otolith-derived SST, T τ + 1 O , the likelihood of the fish being at location, s, at time τ + 1 is

where T s , τ + 1 S S T is the OISST at location s at time τ + 1. The error in the observations, σ, was set to 1°C to account for variability in T^O - T^SST from various sources such as the precision of the SIMS analyses (Hanson et al., 2013), spatial variability in δ ¹⁸ O seawater (LeGrande and Schmidt, 2006), and the accuracy and precision of the applied thermometry equation.

The above defines a discrete likelihood function L(τ + 1) over all states, and the updated estimate of the location conditional on the new observation is given by:

By iterating through all time-steps, τ = 1 … N - 1, applying Equations 1 and 2, a series of probability distributions over the state space for the location of the fish was generated, with the distribution at each step conditional on all prior observations. The initial and final locations were known, so L (1) and L (N) were set to 1 at the location of the assigned river of origin (northeast Scotland) and north Scotland capture point respectively, and 0 at all other locations.

Weekly time-steps were used in the model, but time assignments of the available otolith sample spots did not result in a data point for each time-step for each fish. Accordingly, there were numerous time-steps for which there were no data with which to apply Equation 2. In this case no observation update was performed, or equivalently:

In the absence of conditionality on observations, the resulting distributions were more diffuse than time-steps with associated observation data ( Figure S1 ).

Finally, a backward smoothing recursion was performed such that the probability distributions at each time-step were conditional not only on all prior observations but also all subsequent observations. Specifically, the smoothed probability distribution over the state space at time-step τ, given all observations, was calculated by iterating from τ = N … 1 applying the calculation:

resulting in the final series of smoothed probability distributions, p(τ|∞) defined over Z.

For visualization, we also derived the mean latitude and longitude at each time-step as follows:

λ ¯ = ∑ λ λ ( ∑ ϕ p ( τ | ∞ ) λ , ϕ ) ,

where p(τ|∞)_λ,ϕ is the value of p(τ|∞) at the state with latitude and longitude λ,ϕ.

We simulated simplified migration trajectories, and associated SST observations, to better understand the limitations of the HMM. Body length, otolith sampling frequency and origin locations were all taken from the 1SW salmon used in the study to mimic the data obtained by means of SIMS sampling. Trajectories were constructed by simulating movement from the genetically assigned origin (NE Scotland) with fixed orientation, but with distance travelled between time-steps drawn from a uniform distribution between 0 and Vl_τ . Two orientation profiles were used, either northward to the Norwegian Sea or westward toward the mid-Atlantic ridge then moving longitudinally or latitudinally, respectively, before returning to the origin travelling backward along the same vector as the outward migration. These orientation profiles represented simplified migration trajectories which were not strictly biologically plausible (e.g. fish would not travel in a straight line) but which captured the major axes of movement for this species emigrating from southern Europe within the North Atlantic (e.g. Rikardsen et al., 2021) in a way that allowed some isolation of model performance with respect to movement trajectory. T s , τ S S T was extracted from simulated locations and uncertainty added by sampling from a normal distribution with μ = T s , τ S S T and σ = 1°C. Sampled T s , τ S S T was then thinned to reflect the sampling frequency achieved by means of SIMS.

The HMM was applied to the simulated T s , τ S S T values allowing a direct comparison of the inferred location of the HMM against true values. The performance of the HMM was qualitatively evaluated by visualizing the true and mean inferred locations and quantitatively evaluated by calculating the absolute error (in km) between true and inferred locations over the time series, using the equirectangular projection approximation. For the latter, values were weighted by the likelihood of each inferred location. Specifically, at time-step τ the error in the location estimate, e_τ was calculated by

where

and

where λ τ ' , ϕ τ ' are the true latitude and longitude of the simulated migration at time τ and R is an approximation of the Earth’s radius set to 6378.1 km.

As an additional comparison, we visualise the inferred locations from the HMM against SST thermal bands encompassed by T s , τ S S T ±   1 ∘ C. This provides an indication of the impact of the movement model in reducing the uncertainty of the location estimates relative to a naïve estimate using only T s , τ S S T to infer location (e.g. Hanson et al., 2010).

All data analyses and plotting was performed in the statistical package R (R Core Team, 2022).

In both simulation scenarios, the thermal profiles retrieved from the mean path closely mirrored the pattern of SST along the true simulated path ( Figure 1A ). However, the mean migratory path inferred from the northward simulation more closely represented the true path than that fitted to the westward migratory scenario ( Figure 1B ), indicating the paucity of information within thermal gradients in the westward (longitudinal) direction with which to inform the movement model.

A subset of individual time-steps of the posterior likelihood of locations for the two migratory scenarios are illustrated in Figure 2 together with locations within σ = 1°C of that measured along the true simulated path. Using the HMM to infer location, conditional both on irregular observations of SST and a movement model constrained by fish size, improved the precision of inferred locations relative to estimates based on SST alone – and most especially on the longitudinal axis. The true simulated location was contained within the 90% confidence intervals of the posterior distribution of inferred locations in 76% of the time-steps for the northward simulation and in 56% of time-steps for the westward simulation ( Figure 2 ; Figures S2 & S3 ).

The absolute error (in km) between the true location and all inferred locations (including those time-steps with no associated data), weighted by their probability, varied between 0 and 2070 km. The HMM was more accurate and precise in reconstructing predominantly northward movement than predominantly westward movement. Median error was 519 km [IQR = 326] for the northward simulation and 1321 km [IQR = 871] for the westward simulation.

Figure 3A shows the δ ¹⁸ O profiles from the two salmon examined in the present study. The SIMS sampling transects yielded 34 and 38 δ ¹⁸ O values over 59 (1SW fish) to 112 (2SW fish) sea weeks. Profiles indicated a pattern of peaks and troughs in δ ¹⁸ O values, associated with colder months ( Figure 3B ). SST values extracted from the mean inferred location were consistent with, but less variable than, the otolith-derived estimates.

Both fish were genetically allocated with high confidence (p = 1) to an assignment unit comprising the rivers Oykel/Cassley/Shin in north Scotland. Within the Scottish baseline, this assignment unit has high assignment accuracy (p > 0.89; Gilbey et al., 2016). These three rivers share a common catchment, and are confluent at the Kyle of Sutherland which itself drains to the northern North Sea (57.86  N 03.98 W). Metadata for each fish can be found in Table 1 . Length profiles and emigration/capture locations are shown in Supplementary Material Figures S4 & S5 .

The long-distance out-and-back migration patterns estimated by the HMM were broadly consistent with available information on the migration routes for this species ( Figure 3C ) (Gilbey et al., 2021; Maoiléidigh et al., 2018; Rikardsen et al., 2021).

Individual time-steps of the posterior likelihood of migratory locations for the two wild Atlantic salmon are illustrated in Figure 4 , together with locations within σ = 1°C of that measured in the otolith. The inferred migratory locations of the 1SW individual ( Figure 4A ) were characterized by northward travel past the Faroe Islands into the Norwegian Sea, over-wintering around the Iceland-Faroes Ridge and the Vøring Plateau (67˚N, 04˚E), before return migration to the Scottish coast. The inferred migratory locations for the 2SW salmon ( Figure 4B ) also were characterized by northward travel to the northeast of Iceland, with seasonal transitions between the polar front zone and the Norwegian Sea. By spring of its second year at the sea the positions were located near the Faroe Islands prior to return migration to the Scottish coast. Full outputs can be seen in Figures S6, S7 .

In the present analysis, we explore an approach to inferring movement of individual fish within the marine environment from stable isotope data, which can be applied to basin-scale marine migration. Sequential sampling of ecogeochemical markers derived from biogenic structures, coupled with movement modelling techniques, provides a potentially powerful new method to track migratory animals through time and space. This general approach has been successfully applied to, among others, sub-stock movement of plaice (Darnaude and Hunter, 2018), point-to-point movement of Japanese sardine (Sakamoto et al., 2019) and chum salmon (Matsubayashi et al., 2020). Importantly, our approach allows inference on complex out-and-back migratory patterns over the time period encapsulated in the otolith (in this case, a period of one to two years from juvenile emigrant to return adult stages).

At a basic level, the temperatures estimated from Atlantic salmon otolith δ ¹⁸ O can be used to infer broadscale movement patterns on the basis of seasonal temperature variation within the North Atlantic ( Figure 4 ). Darnaude et al. (2014) demonstrated strong concordance between otolith δ ¹⁸ O values and environmental conditions (temperature and salinity) for wild free-ranging plaice, while von Leesen et al. (2021) showed that SIMS-measured mean otolith δ ¹⁸ O values were well correlated with predicted otolith δ ¹⁸ O values based on temperature and salinity derived from DST-tagged free-swimming cod. However, such inferences would necessarily be poorly defined in the longitudinal direction because thermal gradients are predominantly latitudinal within the migratory range of this species in the eastern North Atlantic (see Reynolds et al., 2007); they would also be restricted to specific time points measured on the otolith. By using a state-space hidden Markov model we found it possible to account for the unobserved fish movement between otolith δ ¹⁸ O observations, and to constrain progression between observation steps using a priori assumptions about salmon movement thereby markedly increasing precision of location estimates in the longitudinal direction and providing inferences of location at all time-steps ( Figures 2 , 4 ).

The available information regarding at-sea distribution of wild Atlantic salmon with which to evaluate the plausibility of our modelled migration distributions is limited both in time and space. Recent analyses of available research and fisheries cruise data found Atlantic salmon postsmolt aggregations concentrated within the North Atlantic Current (April-June), around the Vøring Plateau (June-July), further north off the Barents Sea Plateau (August) and along the Iceland-Faroes Ridge (October) (Gilbey et al., 2021). From tag recaptures and, more recently, genetic screening of fisheries samples, it is apparent that Scottish salmon returning as 1SW and 2SW adults are present around the Faroe Islands at least during November – May. Some 2SW salmon of Scottish origin do latterly migrate westwards to spend their second summer in the western North Atlantic and have been captured in Greenlandic fisheries during August – October prior to their return migration (Bradbury et al., 2016; Ó Gilbey et al., 2017; Maoiléidigh et al., 2018). Our modelled locations in the eastern North Atlantic ( Figures 4 , S6, S7 & S8 ) therefore broadly conform with available fishery and survey-derived capture data. Marine migratory tracks of individual salmon tagged with archival DST are not available for Scottish salmon but such information is available for adult salmon originating from other European countries which had previously spawned, been tagged, and returned to sea. Whilst those data for spawned kelts are not directly comparable to the paths we infer here for maiden salmon undertaking marine migration for the first time, there are remarkable similarities between the tracks inferred here and those derived for salmon migrating especially from Denmark and southern Norway ( Figure S8 ; Rikardsen et al., 2021).

Notwithstanding the fact that migratory inferences derived from our HMM were broadly plausible, it was important to test the limitations of the underlying model and to be mindful of its assumptions. While Darnaude et al. (2014) were able to validate movement inferences from otolith δ ¹⁸ O for plaice using DST, such data are not available for many fish species due to size and/or recapture probability limitations. In the absence of such validation data, we explored simple simulations of movements to test the limitations of our approach.

The performance of the HMM on simulated data highlighted some of the difficulties in recovering location from SST observations alone. In general the inferred temperature profiles were close to the true simulated values ( Figure 1A ); however the HMM had high error when predicting the most extreme longitudinal locations visited by the simulation. Temperature gradients are primarily north-south rather than east-west in the eastern North Atlantic (Reynolds et al., 2007), and consequently, in the presence of the uncertainty in the observation step, there is a lower signal-to-noise ratio contained within a thermal profile extracted from east-west movement. This resulted in the HMM over-smoothing the location estimates, and remaining close to the center of mass during the extremes of the migration.

For simplicity, and to test extreme scenarios, the simulated trajectories included persistent straight line migration over multiple time-steps. However, the uniform circular movement assumption gives equal probability to continue in the same direction as to move in any other direction, making it a sub-optimal movement model for the simulated scenarios. This also impacts the ability of the HMM to infer extreme locations.

Overall, the precision of our method was comparable to other studies using temperature or stable isotope data to infer marine animal distribution (MacKenzie et al., 2011; Chittenden et al., 2013; Guðjónsson et al., 2015; Trueman et al., 2017; Sakamoto et al., 2019; Matsubayashi et al., 2020). We are not aware of other attempts to quantify accuracy of movement inferences based on indirect observations using simulated data. Methods that are GPS or DST-based will be more accurate given the amount and accuracy of data that can now be stored in tags (e.g. Rikardsen et al., 2021). However, it is important to recognize that inferring movement from these data sources also requires assumptions about fish behavior and physiology that are not frequently explored. Whilst a stable isotope-based method of geolocating aquatic animals may be less accurate than more traditional telemetry approaches, it does allow lower resolution data to be readily obtained from a larger number of individuals ( Table 2 ). In the specific case of otolith-derived information, a major advantage is that the method may be applied retrospectively to any fish.

The accuracy and precision of our modelled movement inferences for individual fish appear limited by the predominant direction of travel within the North Atlantic, with constraints on inference regarding longitudinal (east-west) movement. However, there are several avenues which have the potential to improve the accuracy of inferences. A ‘center of attraction’ within the migratory domain could be included as a prior within the model if there is good information on the likely distribution of animals at a particular time (e.g. McClintock et al., 2012). Information derived from, for example, capture-mark-release-recapture studies (Ó Maoiléidigh et al., 2018), research trawls (Gilbey et al., 2021), or even isoscape-based foraging location assignments from scale stable isotope analyses (MacKenzie et al., 2011) may help to parameterize a center of attraction. The failure of the HMM to recover extreme longitudinal locations during simulated migrations suggests that a more sophisticated movement model may improve performance. Specifically, directional persistence and/or variable migration speed in combination with assumptions about stage of migration (e.g. ‘traveling’, ‘searching’, or ‘foraging’) might provide better results (e.g. Michelot et al., 2017). However, parameterizing such behaviors in the HMM would require additional data or strong prior information; we purposefully did not take that approach here due to the datasets being sparse and a lack of robust priors. For some species analysis of otoliths retrieved from individuals fitted with archival DST for the duration of their migration would permit cross-validation of different model formulations, including sensitivity to error in smolt emigration timing. However, for Atlantic salmon, this validation step is not readily available because of the present lack of sufficiently small DST for most juvenile smolts at river emigration. None the less, the use of downstream-migrating smolts that had been uniquely tagged (or recorded telemetrically) near the head of tide would offer the model advantage of confidence in the individual’s date of river emigration.

Some assumptions implicit within our approach may warrant closer consideration. Perhaps the foremost of these is that the error in SST values derived from otolith δ ¹⁸ O values is spatially random. Spatial variation in seawater δ ¹⁸ O values (LeGrande and Schmidt, 2006) must be incorporated into future iterations of the approach, and experimental determination of an Atlantic salmon-specific δ ¹⁸ O fractionation equation may also improve predictions for this species.