A critical assessment of marine predator isoscapes within the southern Indian Ocean

Abstract

Background

Precise and accurate retrospective geolocation of marine predators via their tissues’ isotopic composition relies on quality reference maps of relevant isotopic gradients (“isoscapes”). Additionally, a good working knowledge of any discrimination factors that may offset a marine predator’s isotopic composition from baseline isotopic values, as well as tissue specific retention rates, are imperative. We provide a critical assessment of inter-specific differences among marine predator-level isoscapes within the Indian Sector of the Southern Ocean.

Methods

We combined fine-scale GPS tracking data and concurrent blood plasma δ¹³C and δ¹⁵N values of eight seabird species (three albatross, two giant petrel and three penguin species) breeding at Marion Island to produce species- and guild-specific isoscapes.

Results

Overall, our study revealed latitudinal spatial gradients in both δ¹³C and δ¹⁵N for far-ranging seabirds (albatrosses and giant petrels) as well as inshore-offshore gradients for near-ranging seabirds (penguins). However, at the species level, latitudinal spatial gradients were not reflected in the δ¹³C and δ¹⁵N isoscapes of two and three, respectively, of the five far-ranging species studied. It is therefore important when possible to estimate and apply species-specific isoscapes or have a good understanding of any factors and pathways affecting marine predators’ isotopic composition when estimating the foraging distribution of marine predators via their tissues’ stable isotope compositions.

Conclusions

Using a multi-species approach, we provide evidence of large and regional scale systematic spatial variability of δ¹³C and δ¹⁵N at the base of the marine food web that propagates through trophic levels and is reflected in the isotopic composition of top predators’ tissues.

Background

Some of the greatest threats faced by land-breeding marine predators are experienced at sea. These include bycatch-risk and changes in food availability as a result of competition with fisheries and climate change [1,2,3]. Therefore, to implement effective conservation-based marine spatial planning there is a growing need to better understand the at-sea distribution of marine predators [4,5,6]. This has led to an impressive growth in the number of tracking studies in recent years (reviewed in [7]), often with the general aim of providing policy-relevant information on important habitat for the respective study species [8]. However, dataloggers are still cumbersome for small species (e.g. some burrowing seabird species) and deployment of loggers on study animals requires significant amounts of time in the field, especially when instruments need to be retrieved. Stable isotope ecology as a tool for retrospective geolocation of predator foraging grounds has relatively recently emerged as an alternative and complimentary method to conventional tracking studies [9, 10]. Stable isotope analysis of body tissues is relatively cheap, less demanding in terms of field time and as a result allows for easy sampling of a greater number of individuals than would often be incorporated in tracking studies (e.g. [11]).

The precision and accuracy of retrospective geolocation of marine predators based on their isotopic composition is primarily reliant on two important factors. Firstly, the reliance on the availability of reference maps of the relevant isotopic gradients, known as “isoscapes” [12]. In marine systems, marine predator movement is commonly inferred by linking the ratios of the stable isotopes of carbon (¹³C/¹²C; δ¹³C), and to a lesser extent nitrogen (¹⁵N/¹⁴N; δ¹⁵N), of their tissues to known gradients of δ¹³C and δ¹⁵N values present at the base of their food webs (e.g. [13,14,15]). Across the global oceans, there is a strong negative latitudinal gradient in the δ¹³C values of phytoplankton, from the equator towards the poles [16], as well as from inshore benthic habitats to offshore pelagic habitats [17,18,19]. Whereas gradients of δ¹⁵N are not as strong or predictable, the δ¹⁵N values of phytoplankton tend to be lower or higher in areas of nitrogen fixation (e.g. pelagic oceans) or denitrification (e.g. upwelling regions around coastlines), respectively [20, 21]. Secondly, a good working knowledge of potential discrimination factors which may offset a consumer’s isotopic composition from baseline isotopic values is required. These discrimination factors may vary with diet composition, isotopic averaging as well as physiological fractionation through intermediate trophic levels, isotopic turnover rates and physiological transformation in the consumer [9, 12]. This includes tissue-specific retention times, as isotopic turnover of different tissues varies greatly [22].

Previously, studies which estimated oceanic δ¹³C and δ¹⁵N isoscapes have largely used organisms close to the base of the food web (e.g. [16, 23, 24]) or particulate organic matter [25,26,27]. The isotope ratios of organisms near the base of the food web (e.g. phytoplankton) or particulate organic matter are influenced by broad scale and localized chemical element circulation and physical oceanographic features [28]. Furthermore, due to the high turnover rate of these organisms, their isotope ratios may change daily [29]. However, due to temporal integration of isotopic ratios from the base of the web through to higher predators [30], the stable isotopic compositions of marine predator tissues are not expected to reflect these short-term changes in the baseline but rather reflect more consistent isotope gradients.

The estimation of δ¹³C and δ¹⁵N predator-specific isoscapes requires information on the movement of the predator during the time which the body tissue of interest incorporated its isotope composition (e.g. [14, 31]). While breeding, marine predators such as seabirds are central place foragers, regularly returning to their nest to provide care to their offspring [32]. Therefore, due to the ease of recapture of individuals after single foraging trips and collection of tissues for stable isotope analysis, seabirds represent ideal study species to investigate inter-specific differences in marine predator-level δ¹³C and δ¹⁵N isoscapes. In addition, the dichotomy in foraging mode presented in seabirds, and their associated foraging constraints, provides an opportunity to compare flying (e.g. albatrosses and giant petrels), and diving (penguins) seabirds. The vast inter-specific differences in the at-sea foraging distribution of seabirds furthermore allow for both the investigation of isoscapes over geographically extensive areas as well as inter-specific differences in isoscapes within locations utilised by multiple species.

Several studies have attempted to reconstruct δ¹³C and δ¹⁵N predator-specific isoscapes by combining known location and distributional range [17] or movement of seabirds [11, 14, 31, 33, 34] to temporally matched tissue isotopic values. The blood plasma of wild birds has a half-life of a few days and it is generally assumed that δ¹³C and δ¹⁵N values of seabird blood plasma reflects approximately 7 days bioaccumulation of isotope values of prey prior to sampling [22, 35]. This makes it an ideal tissue to link to fine-scale tracking data as the movements of the birds can be known with high accuracy several days prior to sampling. However, integrating concurrently collected δ¹³C and δ¹⁵N blood plasma values with fine-scale tracking data has rarely been done and only in a single-species context [14, 31]. Here, we combine fine-scale GPS tracking data and concurrent blood plasma isotopic values of five Procellariiformes (albatrosses and giant petrels) and three Sphenisciformes (penguins) species to produce species and guild specific δ¹³C and δ¹⁵N isoscapes within the Southern Indian Ocean. We assigned an individual’s isotopic values to a time integrated mean foraging location and use these data to investigate spatial patterns of a species δ¹³C and δ¹⁵N plasma values. The study aims to investigate the influence of species on the determination of δ¹³C and δ¹⁵N isoscapes as well as differences between predators with different movement modes (flying vs. diving). We hypothesized that the δ¹³C and δ¹⁵N tissue values of most species exhibit a spatial gradient and the extent and strength of this gradient is dependent on the species foraging range with the following predictions:

• 1. Influence of seabird species on the detection of large-scale spatial variability in δ¹³C and δ¹⁵N isoscapes. We predict that due to the vast distances travelled by flying species during a single foraging trip, their δ¹³C and δ¹⁵N tissue values will reflect known baseline latitudinal δ¹³C and δ¹⁵N gradients (e.g. [16, 20]). However, due to sharp changes in baseline δ¹³C values at frontal zones [25], we predict that if individuals of a species do not cross different major fronts within the Indian sector of the Southern Ocean, the δ¹³C values of that species’ tissues will not reflect a spatial gradient.

• 2. δ¹³C and δ¹⁵N isotopic values of major fronts and water zones within seabird species isoscapes. We predict that within the isoscapes of species which crossed one or more major fronts, there will be consistent stepwise increases at major fronts from south to north.

• 3. Influence of seabird species on the detection of fine scale spatial variability in δ¹³C and δ¹⁵N isoscapes. Previously, it has been shown that an inshore/offshore and benthic/pelagic effect can be detected in the isotopic values of seabirds [17, 36, 37]. Thus, due to the more limited foraging range of penguins compared to that of the flying species in this study, we predict that the δ¹³C and δ¹⁵N plasma values of penguins will not reflect a latitudinal gradient, but will rather reflect the inshore/offshore and benthic/pelagic foraging habitat of the species.

Materials and methods

Study site and species

The Prince Edward Archipelago is located in the Indian sector of the Southern Ocean between the sub-Antarctic and Antarctic polar fronts [38]. Ascending from a depth of approximately 3000 m, the archipelago consists of two islands, Marion Island (~ 240 km²) and Prince Edward Island (~ 45 km²), which are located 19 km apart and separated by a shallow inter-island shelf that ranges from 40 to 400 m in depth [39]. These two islands provide breeding grounds for more than five million seabirds and seals [40].

This study includes three of the four albatross species (wandering Diomedea exulans, grey-headed Thalassarche chrysostoma and sooty Phoebetria fusca albatrosses), both giant petrel species (northern Macronectes halli and southern M. giganteus giant petrels) and three of the four penguin species (gentoo Pygoscelis papua, macaroni Eudyptes chrysolophus and eastern rockhopper E. filholi penguins) breeding at Marion Island.

Data collection

Field work was conducted between August to March 2015/16–2017/18, along the south-east coast of Marion Island (46°54′S; 37°45′E) during the breeding seasons of the respective study species (Supplementary Material S1). GPS data loggers (CatLog-S GPS loggers, Perthold Engineering LLC USA, 50 × 22 × 8 mm, ~ 24.5 g) were deployed on individuals that were either incubating or brooding small chicks. In total, loggers were deployed on 176 individuals, which included 32 wandering, 28 grey-headed and 22 sooty albatrosses, 27 northern and 23 southern giant petrels and 17 gentoo, 18 macaroni and 9 eastern rockhopper penguins. Loggers were set to record geographic locations at hourly intervals for flying birds (i.e. albatrosses and giant petrels) and two-minute intervals for penguins.

Loggers were retrieved after 1.2 ± 0.5 foraging trips for flying birds and 1.7 ± 0.9 foraging trips for penguins (Table 1). Upon retrieval of GPS loggers, ~ 1 ml of blood was collected from the tarsal vein of flying birds or the brachial vein of penguins using a sterile heparinised 25 gauge needle. Approximately 0.5 ml of the blood was centrifuged within 3 – 4 h after collection, separated into red blood cells and plasma, stored in 70% ethanol and frozen until preparation for stable isotope analysis [41].

Table 1 Delipidated or normalized plasma δ¹³C and raw plasma δ¹⁵N values of albatrosses, giant petrels, and penguins breeding at Marion Island, which were tracked simultaneously with GPS data loggers during 2015–2018. Number of individuals indicates number of birds with tracks and corresponding stable isotopic values. Number of tracks indicates number of tracks used to estimate mean foraging locations as multiple tracks were recorded for brooding individuals (See Materials and Methods). Values given as mean ± SD (range)

GPS analysis

A foraging trip was defined as the last location on land until the first location back on land. Land-based locations were removed for further analysis. Unrealistic GPS locations were identified and removed using a speed filter algorithm (R package trip; [42]). Unrealistic locations were identified as those requiring movement speeds greater than 135 km h^− 1 [32] for flying birds and greater than 10 km h^− 1 for penguins [43]. Post-filtering, positions were linearly interpolated at one-hour intervals for the flying birds’ trips (R package adehabitatLT; [44]). The filtered data for penguins were processed using a continuous-time correlated random walk (CRAWL) model to estimate the approximate movement track at regular intervals (R package crawl; [45]). The intervals set in the CRAWL models were the same as the intervals at which the loggers were set to record geographic locations (i.e. 2 min). When penguins are underwater, GPS signal is lost, resulting in irregular time intervals between locations; the CRAWL method fits a movement model to estimate locations at regular time intervals rather than assuming linear movement between irregular location estimates [46]. Prior to further analysis, all GPS locations within a 2 km (penguins) or 15 km (flying birds) buffer around the island were removed to avoid an upward bias created by birds leaving and returning from foraging trips. This resulted in removal of all giant petrels that foraged exclusively within penguin and seal rookeries at the island.

Identification of foraging behaviour along GPS tracks

Foraging activity along a seabird’s track is characterised by high sinuosity (i.e. frequent turning) and low flight speeds, and can be distinguished from direct and fast transit to and from the colony [47, 48].

To determine foraging locations along individual tracks of flying seabirds, Expectation Maximization binary Clustering (EMbC) was used (R package EMbC; [49]). This method uses an unsupervised clustering algorithm based on maximum-likelihood Gaussian mixture models that produces biologically interpretable behavioural classifications from flying seabird tracking data [50, 51]. Derived from the turning angle and speed between successive GPS locations, the EMbC determines four behavioural classification categories, obtained from the four combinations of high and low values of turning angle and speed. Locations categorized to have low speeds and high turning angles are considered as the birds ‘actively sitting’ and can be considered as the birds being in a behavioural foraging phase [50].

The EMbC did not perform well on the penguin data. However, for diving predators, speed is a useful proxy of foraging behaviour [52, 53], with slower speeds (reduced horizontal displacement) typically associated with increased foraging activity [54]. Speeds between successive GPS locations were calculated and speeds which were lower than a species average speed were identified as foraging locations.

Stable isotope analysis

The high lipid contents within blood plasma, detected by a C:N mass ratio > 3.5 in tissues for aquatic animals [55], lead to artificially low δ¹³C values [56, 57]. This may lead to misinterpretation when inferring foraging locations of consumers from their δ¹³C tissue values. Mathematical normalization equations developed for a particular species or tissue may not be appropriate for another [55, 58], while chemical extraction methods may artificially increase δ¹⁵N values [59]. Therefore, δ¹³C and δ¹⁵N values were obtained from lipid extracted and raw plasma, respectively.

Blood plasma was dried at 50 °C for 48 h before being powdered using a mortar and pestle. Where possible, each plasma sample was divided into two aliquots: lipids were extracted from one half of the sample while the other half was analysed without lipid extraction. Lipids were removed by immersing powdered plasma in a 2:1 chloroform: methanol solution with a solvent volume three to five times greater than sample volume. Samples were then vortexed for 10 s every 10 min for 1 hr before being centrifuged for 5 min. The supernatant containing lipids was discarded, and samples dried at 50 °C overnight.

For small plasma samples where lipid extraction was not possible (two wandering, six grey-headed and one sooty albatross samples and two gentoo and five macaroni penguin samples), raw material was analysed so that true δ¹⁵N values could be obtained. To calculate lipid corrected δ¹³C values, species-specific mathematical normalization equations were developed using delipidated and raw plasma samples used in this study as well as others collected during the same breeding seasons but not presented here (Pistorius unpub. data), so that sample size could be increased. A species-specific normalized δ¹³C value was calculated using the difference of δ¹³C values (Δδ¹³C) and C:N ratios (ΔC:N) between non-delipidated plasma (δ¹³C and C:N) and delipidated plasma (δ¹³C_del and C:N_del) for each individual (eq. 1 and eq. 2, respectively). A linear regression between Δδ¹³C and the ΔC:N ratio was then calculated, thus giving an intercept (c) and slope (m) for each species-specific equation (eq. 3). Eq. 3 was used to calculate Δδ¹³C for the samples that did not have a lipid free counterpart with ΔC:N for the raw samples calculated as the difference between their C:N and the average C:N of the delipidated samples. Finally, a corrected δ¹³C (δ¹³C_cor) was calculated for samples without a lipid free counterpart using eq. 4. Species-specific equations (eq. 3) and sample sizes are given in Supplementary Material S2.

$$ \Delta {\updelta}^{13}\mathrm{C}={\updelta}^{13}{\mathrm{C}}_{\mathrm{del}}-{\updelta}^{13}\mathrm{C} $$

(1)

$$ \Delta \mathrm{C}:\mathrm{N}=\mathrm{C}:{\mathrm{N}}_{\mathrm{del}}-\mathrm{C}:\mathrm{N} $$

(2)

$$ \Delta{\delta}^{13}{\mathrm{C}} = \frac{\delta \mathrm{C}:\mathrm{N}+\mathrm{c}}{\mathrm{m}}$$

(3)

$$ {\updelta}^{13}{\mathrm{C}}_{\mathrm{cor}}={\updelta}^{13}\mathrm{C}+\Delta {\updelta}^{13}\mathrm{C} $$

(4)

The isotopic values of carbon and nitrogen in aliquots (~ 0.4 mg) of homogenized delipidated and raw plasma samples were determined by combusting samples in a Flash 2000 organic elemental spectrometer via a Conflo IV gas control unit (Thermo Scientific, Germany). All samples were processed at the Stable Light Isotope Unit at the University of Cape Town, South Africa. Replicate measurements of internal laboratory standards indicated minimal standard deviations within and among runs (Merck gel: SD_δ13C = 0.2 ‰, SD_δ15N < 0.1 ‰; valine: SD_δ13C < 0.2 ‰, SD_δ15N = 0.1 ‰; seal bone: SD_δ13C < 0.2 ‰, SD_δ15N < 0.1 ‰). All in-house standards were calibrated against International Atomic Energy Agency standards. Results are presented in the usual δ notation relative to Vienna PeeDee Belemnite and atmospheric N₂ for δ¹³C and δ¹⁵N values, respectively.

Estimation of Isoscapes

Blood plasma has a half-life of a few days [22]. It was therefore assumed that the δ¹³C and δ¹⁵N values of the seabirds’ blood plasma reflected approximately 7 days bioaccumulation of isotope value of prey prior to sampling. The isotope value of a seabird’s tissue is a moving temporal window of what it has ingested [22, 35]. Thus, a geographic location that would represent an individual’s single δ¹³C and δ¹⁵N plasma values was estimated by identifying cells where the individual was likely foraging (see Identification of foraging behaviour along GPS tracks) and weighting each of these cells by the proportion of time spent foraging within a given cell. Calculating the mean of these weighted positions resulted in a time-integrated weighted mean foraging location for each individual. In the case of brooding albatrosses (grey-headed and sooty albatrosses only) and penguins, where more than one foraging trip was recorded, multiple trips were included in further analyses if they did not exceed 7 days prior to blood collection.

The proportion of time spent per cell was calculated (R package trip; [42]). Since the path length and duration of albatrosses and giant petrels was much greater than that of the penguins, two different grid sizes were used for flying and diving seabirds when calculating time spent per cell. The grid cell size was estimated by taking the average distance that the flying or diving species travelled in 2 h (as location interval for flying birds was 1 hr) and rounding up to the nearest 0.05°. This resulted in grid cell sizes of 0.5° for flying birds and 0.05° for diving birds.

Possible geographic gradients were investigated using Pearson’s correlation to identify relationships between δ¹³C and δ¹⁵N values of plasma and mean foraging latitude and longitude (in flying species) or distance to coastline (in diving species). As all penguin species only moved within ~ 0.5° of latitude, their blood plasma isotope values are unlikely to be affected by latitude and distance to coastline was thus used to investigate a possible inshore to offshore gradient. To visualize δ¹³C and δ¹⁵N isoscapes, species-specific geostatistical models were used to interpolate δ¹³C and δ¹⁵N values among the mean foraging locations. Ordinary kriging is a spatial interpolation method which relies on the notion of autocorrelation as a function of distance [60] and has been shown to perform well when estimating isoscapes from point locations [61, 62]. Thus, ordinary kriging (R package gstat; [63]) was used, with best fit variograms identified (R package automap; [64]). Interpolations were made onto the same grid used to calculate time-spent per cell per individual, with extent of the grid being a convex hull around a species mean foraging locations with a 0.5° or 0.05° buffer for flying species and penguins, respectively. Isoscapes estimated for each species were then subtracted from one another to allow for comparison among the isocapes produced for each species, following the methods used by St John Glew et al. [61].

The mean foraging location and plasma isotope values from all species within a guild (i.e. flying birds and penguins) were then combined to produce guild-specific δ¹³C and δ¹⁵N isoscapes. To account for the different number of individuals per species, a bootstrap approach was used. Whereby, 22 and 9 individuals (the lowest number of individuals per flying and penguin species, Table 1) were randomly sampled 1000 times from each of the flying and penguin species, respectively, and used to produce new δ¹³C and δ¹⁵N isoscapes. The mean of these isoscapes was then calculated to present guild-specific isoscapes which consider different number of individuals per species.

To investigate possible δ¹³C and δ¹⁵N values characteristic of main fronts and water masses within the Southern Ocean, isotopic values were extracted within one degree along the sub-tropical, sub-Antarctic and Antarctic polar fronts as well as for sub-tropical, sub-Antarctic, polar frontal and Antarctic water zones. Front positions were estimated from satellite derived sea surface height (Ssalto/Duacs produced and distributed by the Copernicus Marine Environment Monitoring Service http://marine.copernicus.eu) averaged over December of the 3 years of the study (2015–2017). Following Swart et al. [65], the fronts were identified from the following sea surface heights: sub-tropical: 0.92 m, sub-Antarctic: 0.03 m and Antarctic polar: − 0.48 m. The four water zones identified were considered as the following: the subtropical (north of the sub-tropical front), the sub-Antarctic (between the sub-tropical and sub-Antarctic fronts), the polar frontal zone (between the sub-Antarctic and Antarctic polar fronts) and the Antarctic (south of the Antarctic polar front). In addition, the δ¹³C and δ¹⁵N values of the Antarctic polar (51°S) and sub-tropical fronts (42°S) were back calculated using estimated regression equations following Jaeger et al. [14].

All data analyses were performed using R version 3.6.3 [66]. All values are reported as mean ± standard deviation and significance is specified as p ≤ 0.05.

Results

Over the duration of the study, 176 birds were tracked, and their blood plasma analysed for δ¹³C and δ¹⁵N values (Table 1).

Influence of seabird species on the detection of large-scale spatial variability in δ¹³C and δ¹⁵N isoscapes

The δ¹³C blood plasma values of all flying seabird species were positively correlated to their mean foraging latitude and significantly so with the exception of the sooty albatross (Table 2, Fig. 1). The correlation between the mean foraging latitude and δ¹³C plasma values of the grey-headed albatross was the strongest (R = 0.83) followed closely by the wandering albatross (R = 0.78) and southern (R = 0.37) and northern (R = 0.37) giant petrels. On the other hand, none of the species δ¹³C plasma values and mean foraging longitudes were significantly correlated (Table 2; Fig. 2). These findings were apparent in δ¹³C isoscapes interpolated for the wandering and grey-headed albatrosses where a clear north to south gradient of δ¹³C values can be seen (Fig. 2a and b). The comparison between the wandering and grey-headed albatross isoscapes showed that the grey-headed albatross δ¹³C isoscape was ~ 1‰ greater than that of wandering albatross within the area of overlap, that is north of the sub-tropical front (Fig. 3). This was reflected with the species regression equations (Fig. 1) with greater δ¹³C values being predicted for the grey-headed albatross than the wandering albatross at similar latitudes. Although the correlation between mean foraging latitude and δ¹³C values of northern and southern giant petrels were significant, the correlation was not as apparent in the interpolated surface (Fig. 2d and e). No pattern could be discerned from the interpolated surface resulting from the sooty albatross’s mean foraging locations and δ¹³C values (Fig. 2c).

Table 2 Pearson’s correlation coefficient (R), p-value (P) and formula resulting from a Pearson’s correlation between plasma δ¹³C and δ¹⁵N values of albatrosses and giant petrels versus latitude (lat) and longitude (lon) of their mean foraging locations calculated from GPS tracks. ^a indicates a significant correlation

Fig. 1

Plasma δ¹³C (left) and δ¹⁵N (right) values of wandering (WA), grey-headed (GHA) and sooty (SA) albatrosses, northern (NGP) and southern (SGP) giant petrels versus their mean foraging a) latitudes and b) longitudes calculated from GPS tracks. Only significant correlations are shown (see Table 2)

Fig. 2

Plasma δ¹³C (left) and δ¹⁵N (right) isoscapes of a) wandering, b) grey-headed and c) sooty albatrosses, and d) northern and e) southern giant petrels breeding at Marion Island, Prince Edward Archipelago (black triangle) interpolated using ordinary kriging from the isotopic values of the respective birds’ plasma which were simultaneously tracked with GPS data loggers. Black points represent an individual’s mean foraging location. Positions of the sub-tropical (STF), sub-Antarctic (SAF) and Antarctic polar (PF) fronts are indicated

Fig. 3

Difference between wandering (WA), grey-headed (GHA) and sooty (SA) albatross as well as northern (NGP) and southern (SGP) giant petrel plasma δ¹³C (top right diagonal) and δ¹⁵N (bottom left diagonal) isoscapes. Differences were calculated by subtracting the isoscape of the species indicated on the vertical axis from the species indicated on the horizontal axis. Due to the small overlap between isoscapes of the giant petrels their comparison is not shown

The δ¹⁵N plasma values of two of the five flying species were significantly and positively correlated to mean foraging latitude, namely the wandering albatross and southern giant petrel (Table 2, Fig. 1). Contrastingly, the δ¹⁵N plasma values of sooty albatrosses, which were not significantly correlated to mean foraging latitude (Table 2), were significantly correlated to mean foraging longitudes (Table 2, Fig. 1b). The resulting interpolated surface of mean foraging location and δ¹⁵N plasma values of the wandering albatross and southern giant petrel showed clear gradients of δ¹⁵N values from north to south (Fig. 2a and e). Although, the resulting wandering albatross δ¹⁵N isoscape was ~ 1‰ greater than that of the southern giant petrel across the interpolated surface (Fig. 3). On the contrary, the interpolated surfaces resulting from the mean foraging locations and δ¹⁵N plasma values of sooty albatrosses showed a gradient from west to east (Fig. 2c).

The δ¹³C and δ¹⁵N plasma values of all flying seabirds combined were significantly correlated with mean foraging latitude (R = 0.73, p < 0.01 and R = 0.44, p < 0.01, respectively; Table 2, Fig. 1a, b) but not to mean foraging longitude (R = 0.12, p = 0.17 and R = 0.24, p = 0.19, respectively). The resulting isoscapes, interpolated from mean foraging locations and plasma δ¹³C values of all flying seabirds showed a clear south to north gradient in δ¹³C values (Fig. 4a). As δ¹⁵N plasma values of the flying seabirds were not as strongly correlated to latitude (R = 0.44), the resulting δ¹⁵N isoscape did not result in a clear south – north gradient, relative to the δ¹³C isoscape (Fig. 4a). However, there was also an increase in δ¹⁵N values from south – north with the higher δ¹⁵N values in the north-west area of the isoscape compared to those of the north-east (Fig. 4a).

Fig. 4

a Plasma δ¹³C (left) and δ¹⁵N (right) isoscapes calculated as the mean of 1000 bootstrapped surfaces interpolated using ordinary kriging which randomly sampled twenty-two individuals from all five flying species to account for different numbers of individuals per species. b Isoscapes, using the same bootstrap technique, modelled without the species which did not exhibit a significant latitudinal gradient (i.e. sooty albatrosses were removed from the δ¹³C isoscape and sooty and grey-headed albatrosses and northern giant petrels were removed from the δ¹⁵N isoscape). Position of the sub-tropical (STF), sub-Antarctic (SAF) and Antarctic polar (PF) fronts and Marion Island (black triangle) are indicated

As a post-hoc analysis, all species for which δ¹³C (sooty albatrosses) and δ¹⁵N (grey-headed and sooty albatrosses and northern giant petrel) plasma values were not significantly related to their mean foraging latitude were removed and the remaining isoscapes were re-interpolated (Fig. 4b). This resulted in a similar δ¹³C isoscape to the previously estimated δ¹³C isoscape (Fig. 4a). However, the resulting δ¹⁵N isoscape now showed a clear north to south gradient (Fig. 4b).

δ¹³C and δ¹⁵N isotopic values of major fronts and water zones within seabird species isoscapes

From each species-specific isoscape derived from flying seabirds (Fig. 2), mean isotopic values of the sub-tropical, sub-Antarctic and Antarctic polar fronts as well as four geographical water zones (sub-tropical, sub-Antarctic, polar frontal and Antarctic water zones) were estimated (Fig. 5, Supplementary Material S3). There were noticeable increases in δ¹³C values at the Antarctic polar and sub-Antarctic fronts within isoscapes interpolated for wandering and grey-headed albatrosses as well as at the sub-tropical and sub-Antarctic fronts compared to surrounding water zones (Fig. 5, Supplementary Material S3). With regards to δ¹⁵N values, there were noticeable decreases for water zones from north to south within the interpolated areas for the wandering albatross and southern giant petrel, whereas the remaining species (sooty and grey-headed albatrosses and northern giant petrel) showed no clear trend (Fig. 5, Supplementary Material S3). Following Jaeger et al. [14], regression equations calculated here (Table 2) were also used to back calculate the δ¹³C and δ¹⁵N values in plasma of the Antarctic polar and sub-tropical fronts at 51°S and 42°S (Table 3). The resulting δ¹³C values varied much less at both the Antarctic polar (− 21.6 ± 0.6 ‰) and sub-tropical fronts (− 20.3 ± 0.9 ‰) than the resulting δ¹⁵N values (12.5 ± 1.4 ‰ and 13.1 ± 1.8 ‰; respectively).

Fig. 5

The means and standard deviations of (a) δ¹³C and (b) δ¹⁵N values estimated for seven water zones and fronts: Antarctic zone (AZ), Antarctic polar front (APF), polar frontal zone (PFZ), sub-Antarctic front (SAF), sub-Antarctic zone (SAZ), sub-tropical front (STF) and sub-tropical zone (STZ) from the species-specific isoscapes presented in Fig. 2: wandering (WA), grey-headed (GHA) and sooty (SA) albatrosses, northern (NGP) and southern (SGP) giant petrels

Table 3 Plasma δ¹³C and δ¹⁵N values calculated using species-specific regression equations (Table 2) for the Antarctic polar (PF, 51°S) and sub-tropical (STF, 42°S) fronts following the methods of Jaeger et al. [14]

Fine scale spatial variability in δ¹³C and δ¹⁵N isoscapes estimated from seabird tissues

No significant correlations were found between the distance to the coast and mean foraging location of the penguin species and their corresponding δ¹³C and δ¹⁵N plasma values (Table 4, Figs. 6 and 7). Except for the δ¹³C values of the eastern rockhopper penguin. However, the δ¹⁵N plasma values of all penguins combined was significantly correlated with the distance between their mean foraging locations and the coastline of Marion Island (Table 4; Figs. 6 and 8). Furthermore, both δ¹³C and δ¹⁵N penguin-specific isoscapes showed a clear gradient of higher values closer to the island, in waters shallower than 200 m, to lower values further away, in deeper waters (Fig. 8).

Table 4 Correlation coefficient (R), p-value (P) and formula resulting from a Pearson’s correlation between δ¹³C and δ¹⁵N plasma values of individual penguin species as well as all species combined versus the distance between their mean foraging locations and the coastline of Marion Island (distTOcoast). ^a indicates a significant correlation

Fig. 6

Plasma (a) δ¹³C and (b) δ¹⁵N values of gentoo (GP), macaroni (MP) and eastern rockhopper (ERP) penguins and all penguin species (ALL) combined versus the distance between their mean foraging locations and the coastline of Marion Island. Only significant correlations are shown

Fig. 7

Plasma δ¹³C (left) and δ¹⁵N (right) isoscapes of a) gentoo, b) eastern rockhopper and c) macaroni penguins interpolated using ordinary kriging from the isotopic composition of the respective birds’ plasmas which were simultaneously tracked with GPS data loggers. Points represent mean foraging locations. Isobaths (grey) are shown at 200 m intervals with 200 m isobath in black

Fig. 8

Plasma δ¹³C (left) and δ¹⁵N (right) isoscapes of three penguin species breeding at Marion Island, Prince Edward Archipelago calculated as the mean of 1000 bootstrapped surfaces interpolated using ordinary kriging which randomly sampled 9 individuals from each penguin species to account for different number of individuals per species. Isobaths (grey) are shown at 200 m intervals with 200 m isobath in black

Discussion

By combining fine-scale tracking data and stable isotope analysis, we present a critical assessment of concurrent δ¹³C and δ¹⁵N gradients across a marine predator assemblage. This study revealed latitudinal spatial gradients in both δ¹³C and δ¹⁵N plasma values for far-ranging seabirds (i.e. albatrosses and giant petrels) as well as inshore/offshore gradients for near-ranging seabirds (penguins). Despite demonstrating clear potential for using seabird tissue samples for retrospective geolocation, the inter-specific differences in δ¹³C and δ¹⁵N isoscapes cautions against using non-species-specific isoscapes for studying marine predator foraging distributions.

Influence of seabird species on the detection of large-scale spatial variability in δ¹³C isoscapes

Strong positive correlations were found between the mean foraging latitudes and δ¹³C plasma values of two out of the five flying seabird species, namely the wandering and grey-headed albatrosses. Differences between these two species isoscapes (Figs. 2, 3 and 4) as well as the lack of gradients found within the δ¹³C isoscapes of the remaining species is likely due to differences in diet [67] and the use of different oceanographic features as foraging grounds [68,69,70].

• 1. Influence of diet on a species-specific δ¹³C isoscape

It was difficult to discern any spatial gradient within the δ¹³C isoscapes of the northern and southern giant petrels and the sooty albatross. This lack of a strong latitudinal gradient in the δ¹³C plasma values of these species may be attributed to their more unpredictable diet. For example, at Marion Island the diet of all three species may contain large proportions of seabird and marine mammal carrion [71, 72]. The relatively narrow range of both the southern giant petrel and sooty albatross δ¹³C plasma values (Table 1) may suggest that either they preyed on the same seabird carrion or different carrion that had been foraging in similar areas. Without concurrent stomach content sample analysis this is largely speculation, however the lack of spatial gradient found for the δ¹³C and δ¹⁵N plasma values of the southern giant petrel and sooty albatross does indicate that caution is needed when using non-species specific isoscapes to infer foraging locations of seabirds from their δ¹³C tissue values.

• 2. Influence of foraging habitat on a species-specific δ¹³C isoscape

Grey-headed albatrosses from Marion Island remain largely south of the sub-tropical front (this study [68]), whereas wandering albatross (this study [69]) frequently move north of the sub-tropical front. However, both species are likely foraging within and around biologically productive eddies that result from interactions between the Agulhas Return Current, sub-tropical and sub-Antarctic fronts [73], or within the fronts themselves [68, 74, 75]. As δ¹³C values at the base of marine food webs have been shown to significantly correlate with primary productivity and nutrient availability [25, 27], birds foraging around or within mesoscale eddies originating from the same oceanographic feature are likely to have elevated δ¹³C tissue values compared to species that forage outside of fronts and eddies. Thus, the higher δ¹³C plasma values at lower latitudes for the grey-headed albatross compared to other species (Figs. 3 and 4) may result from the grey-headed albatross preferentially foraging in highly productive eddies [68]. This is further supported by the similar δ¹³C values found for the wandering and grey-headed albatrosses and northern giant petrel within the sub-tropical front as well as within the sub-Antarctic and sub-tropical water zones (Fig. 5; Supplementary Material S3).

• 3. Influence of baseline δ¹³C gradient on a species-specific isoscape

The δ¹³C combined species isoscape as well as that of the wandering albatross closely aligns with previously published isoscapes for this region, which used similar methods but only with a single species (wandering albatross; [14]), used whole blood of close-ranging seabirds [17] or baseline values (e.g. [16]). Previous studies, which have used a coupled physics-biochemistry model [16] or satellite tracking coupled with predator tissue values [14] to produce δ¹³C isoscapes for the Indian sector of the Southern Ocean, have shown uniform south-north latitudinal gradient in δ¹³C values. In agreement, presented here is a large scale south-north δ¹³C isoscape, albeit with slightly non-uniform latitudinal gradient in δ¹³C values. The isoscape presented by Magozzi et al. [16] had a greater range of δ¹³C values between southern Africa and Antarctic (− 32 to − 16 ‰) compared to our isoscapes (− 23 to − 18 ‰). Our predator-based isoscape may differ to that of Magozzi et al.’s [16] for two reasons. Firstly they modelled the distribution of δ¹³C values of phytoplankton and not seabirds. Secondly, Magozzi et al. [16] showed that in the Indian sector of the Southern Ocean, δ¹³C values of phytoplankton become increasingly more negative with the transition of summer to winter months. Therefore, the narrower range of δ¹³C values obtained in our δ¹³C isoscape can be related to the enrichment of δ¹³C values due to fractionation through multiple trophic levels from phytoplankton to seabirds [55, 76], as well as the temporal integration of spatially and temporally variable δ¹³C phytoplankton values [77].

• 4. Other influences on seabirds δ¹³C isoscapes

The isoscape map and regression equation of Jaeger et al. [14] for the wandering albatross tracked from Crozet Archipelago (using similar methods used here), revealed higher δ¹³C values with a steeper δ¹³C gradient (range: − 25.0 to − 19.1 ‰) than we found for Marion Island wandering albatrosses (range: − 23.6 to − 18.4 ‰) over similar latitudes (~ 38–58°S). This could be an artefact of differences in diet between the two populations of wandering albatrosses [67]. However, the mean δ¹⁵N plasma value of the wandering albatross in our study (mean: 14.0 ‰; range: 12.1 to 15.5‰) and Jaeger et al.’s [14] study (14.0 ‰; 11.2 to 15.8 ‰) was the same, indicating that individuals from the two populations were foraging at similar trophic levels. Rather, a difference in the methods used to remove lipids from plasma between the two studies (cyclohexane [14] versus a 2:1 chloroform to methanol mix (our study)) may explain the slightly more negative δ¹³C values reported by Jaeger et al. [14]. This is supported by the fact that Jaeger et al. [14] reported a higher C:N ratio (4.1) than we found for the wandering albatross plasma (3.5), which indicates that Jaeger et al. [14]‘s samples had a higher percentage of lipids than ours. Because lipids can lead to artificially low δ¹³C values [56, 57], this may explain their lower δ¹³C plasma values.

These finding once again highlight the importance of removing or accounting for lipids when inferring the retrospective foraging location of seabirds from their δ¹³C tissue values [56, 78]. Especially, since during different life-history stages, seabirds may undergo differential levels of nutritional stress such as periods of fasting, which may result in the increase of lipids circulating in the bloodstream [79]. As previously mentioned, this can lead to artificially low δ¹³C values [56, 57]. However, as δ¹³C blood values are not (or little) otherwise affected when seabirds undergo nutritional stress [79], we believe the δ¹³C isoscapes produced here are transferable across life-history stages of seabirds.

Influence of seabird species on the detection of large-scale spatial variability in δ¹⁵N isoscapes

This study revealed a south to north gradient of the wandering albatross and southern giant petrel δ¹⁵N plasma values whereas δ¹⁵N plasma values of the grey-headed and sooty albatrosses and northern giant petrel lacked a south to north gradient. However, the δ¹⁵N plasma values of the sooty albatross did present a relatively weak west to east gradient. The δ¹⁵N value of predators’ tissues is a combination of their trophic position and δ¹⁵N values of organisms at the base of the food chain [55, 76]. Therefore, as presented by Jaeger et al. [14], there are two possible and non-exclusive hypotheses to explain systematic spatial variability in seabird δ¹⁵N tissue values: firstly, a dietary shift of the predators throughout the predator’s distributional range and secondly a change in the baseline nitrogen values.

• 1. Influence of diet on a species-specific δ¹⁵N isoscape

Jaeger et al. [14] provided evidence for the first hypothesis by using previous studies of the wandering albatross diet to show that, although the wandering albatross breeding at Crozet Archipelago feed on squid across their range, Antarctic squid species feed at lower trophic levels than sub-Antarctic and tropical squid species [24], which may result in the latitudinal gradient of δ¹⁵N plasma values. The same data is not available for the wandering albatross breeding at Marion Island, but it is reasonable to assume that this is also a factor driving the spatial gradient of the wandering albatross δ¹⁵N plasma values within this study.

• 2. Influence of baseline δ¹⁵N gradient on a species-specific δ¹⁵N isoscape

Although diet is undoubtedly an important contributing factor to the δ¹⁵N plasma values of the birds within this study, this study provides evidence towards Jaeger et al. [14]‘s second hypothesis: baseline δ¹⁵N values influence the latitudinal gradient in δ¹⁵N plasma values of seabirds. Evidence supporting this argument is that the wandering albatross and southern giant petrel had similar regression slopes between mean foraging latitudes and δ¹⁵N plasma values. The difference between these two species δ¹⁵N isoscapes was homogenously 1 ‰ across overlapping areas (Fig. 3), which may reflect difference in their diets [67]. Furthermore, the wandering albatross and southern giant petrel δ¹⁵N isoscapes as well as the δ¹⁵N isoscape resulting from the combination of their data (Fig. 4b) resembled the latitudinal gradients previously found for δ¹⁵N values across the Indian sector of the Southern Ocean for seabirds (wandering albatross: [14]) and baseline values [20]. The values within these δ¹⁵N isoscapes ranged between ~ 12 ‰ in the south and ~ 15 ‰ in the north, closely matching the ~ 2 ‰ modelled gradient in baseline δ¹⁵N values between Antarctica and the coastline of southern Africa [20]. This indicates that these spatially variable baseline values propagate up the food chain and can be observed in the δ¹⁵N tissue values of the seabirds. These are important findings for the use of stable isotope analysis for the retrospective geolocation of marine consumers. Previously, δ¹³C values of predator tissues have been used for retrospective geolocation, however, agreement among trends found within δ¹⁵N isoscapes in this study and two others [14, 20] provide evidence toward the use of δ¹⁵N tissue values as an additional tool for retrospective geolocation. It further emphasizes the importance of taking the spatial gradients of δ¹⁵N values into account when inferring trophic level from δ¹⁵N tissue values. Of course, it is also important to consider the diet composition and corresponding trophic level of predators and how this may vary within the population (e.g. among age-classes or between sexes) since there is an enrichment of 0.5–1 ‰ for δ¹³C values and 3–5 ‰ for δ¹⁵N values per trophic level [55, 76].

• 3. Influence of foraging habitat on a species-specific δ¹⁵N isoscape

Even though the δ¹⁵N plasma values of the sooty albatross did not show any correlation with latitude, the correlation between their δ¹⁵N plasma values and longitude may still provide evidence towards the second hypothesis (i.e. baseline δ¹⁵N gradient influence on marine predator δ¹⁵N tissue values). Foraging of sooty albatross breeding at Marion Island within 20–25° longitude is largely concentrated along the South West Indian Ridge [80]. The interaction between the fast-flowing Antarctic Circumpolar Current and this ridge results in zones of upwelling as well as meandering of fronts and eddy formation [81]. Consequently, pockets of high productivity all along the ridge may cause periodically higher δ¹⁵N values at the base of the food web [82]. This may explain the longitudinal gradient in the δ¹⁵N values between individual sooty albatross feeding closer to the ridge (between 20 and 25°E) and further downstream of it. Without concurrent stomach content samples of the tracked birds or samples of phytoplankton in the birds’ foraging areas, no definite conclusion can be drawn. Compound specific stable isotope analysis could be used in future investigations to overcome these uncertainties [83].

Fine scale spatial variability in δ¹³C and δ¹⁵N isoscapes estimated from seabird tissues

This is the first study to combine the use of GPS dataloggers and the carbon and nitrogen stable isotope values of penguins to investigate fine-scale isotopic variability across space. It was found that even though all three penguin species stayed much closer to the island (0.5° latitudinal range) than the albatrosses and giant petrels (30° latitudinal range), they had lower δ¹³C plasma values than most of the flying seabirds (Table 1). This indicates that, while their isotopic values may be affected by latitudinal gradients in δ¹³C and δ¹⁵N stable isotopes (this study [14, 16, 20]), there are other factors impacting their isotopic values. These factors are most likely related to the major differences in foraging ecology between penguins and flying seabirds as well as regional physical features which contribute to fine scale spatial variability in baseline isotopic values at the Prince Edward Archipelago compared to elsewhere [18, 84, 85].

The inshore consumers at the Prince Edward Archipelago, like the gentoo penguin [86], derive their energy from external and internal sources. The external source is a result of allochthonous prey species being advected and concentrated around the archipelago by the Antarctic Circumpolar Current [87, 88]. Internally, due to the island mass effect, there is localized enhancement of biological productivity [87, 89]. Due to these two different energy sources, autochthonous species (e.g. microphytoplankton, benthic krill Nauticaris marionis, and fish), particularly benthic species, have higher values of both δ¹³C and δ¹⁵N compared to allochthonous species (e.g. sub-Antarctic krill Euphausia vallentini) found in the vicinity of the archipelago or in the surrounding pelagic waters [18, 84, 85].

The δ¹³C and δ¹⁵N isoscapes of penguins are therefore a combination of their diet and foraging habitat. The gentoo penguin breeding on Marion Island are strictly inshore foragers and are typically benthic foragers but they can also forage pelagically [86]. Their diet at Marion Island is highly variable, composed of mainly benthic fish and an autochthonous benthic crustacean (i.e. Nauticaris marionis) but can contain large proportions of pelagic crustaceans and juvenile fish if they are available [90, 91]. In contrast, the crested penguins (macaroni and eastern rockhopper penguins) are typically pelagic foragers with the eastern rockhopper penguin also exhibiting the ability to forage benthically [46, 53]. Although the brooding macaroni penguin reach further distances away from the island than the brooding eastern rockhopper penguin (this study [53]), the diet of both species is almost completely dominated by allochthonous pelagic krill species (i.e. Euphausia vallentini and Thysanoessa vicina), with small amounts of fish being eaten by both species ([92]; Pistorius unpub. data; Department of Environment, Forestry and Fisheries unpub. data). Thus, the inshore/offshore gradient in both δ¹³C and δ¹⁵N plasma values of penguins in this study is largely related to the benthic and autochthonous foraging ecology of the gentoo penguin and the pelagic, allochthonous foraging ecology of the crested penguins. Cherel and Hobson [17] similarly attributed the higher δ¹³C tissue values of shelf-feeding gentoo penguin at the Kerguelen Archipelago to local enhanced productivity and benthic feeding. In addition, although it has not been investigated for seabirds, there is the possibility of foraging depth having an impact on the isotopic values of penguins’ tissues, as seen in tunas [93].

Conclusion

Using a multi-species approach this study shows that large and regional scale systematic spatial variability of δ¹³C and δ¹⁵N values at the base of the marine food web propagates through trophic levels and is reflected in the isotopic values of top predators’ tissues. This study provides evidence for the effect of baseline δ¹⁵N values on predator δ¹⁵N tissues values. This result emphasizes the importance of considering δ¹⁵N oceanic isoscapes in studies that incorporate stable isotopic values of marine predators to investigate their trophic ecology. Although these baseline values propagate through the food chain and are reflected across multiple marine predator species, they are not reflected in all species. This may be as a result of the movement of prey species or foraging of individuals across major water zones and fronts. It is therefore important when possible to estimate and apply species-specific isoscapes or have a good understanding of any factors and pathways affecting marine predators’ isotopic values when studying marine predator stable isotope ecology. It is further important to note that this study was performed using the stable isotope values of seabirds’ plasma, which may not reflect stable isotope values of other tissues such as whole blood or feathers [14]. With seasonal variation in the gradient of baseline δ¹³C highlighted in Magozzi et al. [16], future investigation should include intra- and inter-annual variability in stable isotope compositions detected at the top predator level.