Comparing trophic position estimates using bulk and compound specific stable isotope analyses: applying new approaches to mackerel icefish Champsocephalus gunnari

Sample collection

Sampling took place between January and March 2019 onboard the RV Kronprins Haakon and the FV Antarctic Endeavour in waters off 60°39′S and 45°00′W. Field sampling was approved by the Instituto Antartico Chileno (RT_68-18) and by the Office of the Commissioner from the Falkland Islands Government House (2018/046). For bulk stable isotope analysis, data were collected as previously described in Canseco et al. (2023). Specifically, we sampled muscle tissue from 284 individual C. gunnari, 75 pooled samples (combination of 3–5 similarly-sized individuals) of its main prey E. superba, muscle tissue from E. antarctica (n = 30), a secondary fish prey item, and pooled samples from invertebrate prey species A. maxima, T. gaudichaudii, and E. frigida (total n = 50). A total of 23 particulate organic matter (POM) samples were collected by filtering seawater (4 l) through pre-combusted 0.7 µm GF/F Whatman filters. A total of 32 samples were analyzed for nitrogen CSI-AA, distributed as follows: five samples per size-class of C. gunnari (n = 10), five samples per size-class of E. superba (n = 10), three samples of T. gaudichaudii, three samples of E. frigida, three samples of A. maxim a and three samples of E. antarctica. No samples of POM were analyzed for amino acid δ¹⁵N.

Stable isotope analysis

To analyze the stable isotope ratios of invertebrate species and particulate organic matter (POM), we homogenized 3-5 similarly-sized individuals (pooled samples) from the same species and tow, as previously performed by Canseco et al. (2023). In the case of C. gunnari and E. antarctica, muscle tissue was collected from individual animals. Muscle tissue and invertebrate pooled samples were rinsed with freshwater, oven dried for 48 h (60 °C) and homogenized, 1 mg of homogenized powder was encapsulated in tin capsules (Jardine et al., 2003). Stable isotope ratios δ¹³C and δ¹⁵N were analyzed using a Thermo Scientific Delta V mass spectrometer with a dual inlet and Conflo IV interface connected to a Costech 4010 elemental analyzer (EA), and a high-temperature conversion elemental analyzer (TCEA) at the Center for Stable Isotopes of the University of New Mexico (UNM). Replicated measures of reference materials (UNM CSI protein standard #1, UNM CSI protein standard #4) showed that precision (±SD) for δ¹³ C and δ¹⁵N values were 0.04 and 0.11 ‰ for carbon and nitrogen respectively, as described in Canseco et al. (2023). We followed Coplen’s (2011) guidelines for stable isotopes reporting, and our measurements are relative to international standards atmospheric-N₂ for δ¹⁵N and Vienna PeeDee Belemnite for δ¹³C.

Compound specific stable isotope analysis

Amino acids were made suitable for gas chromatography-combustion-isotope ratio mass spectrometry (GC-C-IRMS) by derivatization as N-acetyl methyl esters (NACME; Corr, Berstan & Evershed, 2007; Yarnes & Herszage, 2017). Prior to derivatization, amino acids were liberated from sample material proteins by acid hydrolysis. GC-C-IRMS was performed on a Thermo Trace GC 1310 gas chromatograph coupled to a Thermo Scientific Delta V Advantage isotope-ratio mass spectrometer via a GC IsoLink II combustion interface. Each sample was analyzed twice, but, if the initial measurements showed a very large measurement error, more replicates were analyzed. Replicates of the reference materials for the quality and assurance were measured every five samples, as performed by Athaillah, Yarnes & Wang (2023). The individual amino acids of high purity utilized in quality assurance mixtures were individually calibrated by elemental analysis-isotope ratio mass spectrometry (EA-IRMS) using the primary isotopic reference material for each element (i.e., VPDB for δ¹³C and air for δ¹⁵N).

EA-IRMS was performed using secondary reference materials calibrated against certified standard reference materials from the US Geological Service (USGS), the National Institute of Standard and Technology (NIST), and the International Atomic Energy Agency (IAEA) (i.e., IAEA-600, USGS40, USGS41, USGS42, USGS43, USGS61, USGS64, and USGS65). Calibration procedures for CSIA were performed as previously described in Athaillah, Yarnes & Wang (2023). They were applied identically across reference and sample materials, and initially, the isotopic values for all amino acids were adapted in order to obtain the already known isotopic composition of an internal reference material (e.g., Nor). Afterwards, the stable isotopes of individual amino acids were adjusted using the initial reference mixture for quality assurance, UCD AA1, as a basis, in order to obtain the already known isotopic composition of each amino acid within the mixture. Finally, all measurements were standardized to the primary reference materials for δ¹⁵N (i.e., Air) using the second quality assurance reference mixture, UCD AA2. The final evaluation of quality was founded on the precision and unbiasedness of the control materials, which contained a calibrated amino acid mixture, UCD AA3, as well as multiple natural substances.

For summarizing nitrogen SIA results, amino acids were divided into “trophic” and “source” ones, following Whiteman et al. (2019).

Trophic position estimation

Trophic position (TP) was determined following a Bayesian approach (see below), employing either bulk- or compound-specific stable isotope analysis methods. To compute the bulk trophic position for all sampled species, we utilized the equation proposed by Cabana & Rasmussen (1996) (Ec.1). ${\displaystyle {TP}_{\mathrm{bulk}}=\lambda +\left({\delta }^{\text{15}}{N}_{\mathrm{Consumer}}-{\delta }^{\text{15}}{N}_{\mathrm{Baseline}}\right)/TE{F}_{\mathrm{bulk}}}$ Where λ is the trophic position of our baseline organism, δ¹⁵N_Consumer is the nitrogen isotopic value of each consumer, δ¹⁵N_Baseline is the nitrogen isotopic value of our baseline, TEF is the nitrogen trophic enrichment factor (TEF) of 3.4‰ (SD =0.51) (Post, 2002). Our baseline for this calculation was the δ¹⁵N values of POM as estimated by Post (2002). POM δ¹⁵N values serve as the foundation of the pelagic ecosystem, playing a crucial role in enhancing productivity off the South Orkney Islands which is key for E. superba, the main prey of C. gunnari (Everson, 2008; Zhu & Zhu, 2022; Canseco et al., 2023). Additionally, to test the accuracy of TP estimates, we calculated C. gunnari TP using the equation shown above although the baseline (δ¹⁵N_Baseline) for this calculation was E. superba, as it has been mentioned that it is the well-known, main prey.

When estimating trophic position using compound-specific stable isotope analysis, we employed two different methods following the general formula proposed by Chikaraishi et al. (2009): ${\displaystyle {TP}_{\mathrm{Trc}-\mathrm{Phe}}=1+\left({\delta }^{\text{15}}{N}_{\mathrm{Trc}}-{\delta }^{\text{15}}{N}_{\mathrm{Phe}}-{\beta }_{\mathrm{Trc}-\mathrm{Phe}}\right)/{TEF}_{\mathrm{Trc}-\mathrm{Phe}}}$ Where δ¹⁵N_Trc are the δ¹⁵N values in the trophic AAs glutamic acid (Glx) or proline (Pro), δ¹⁵N_Phe is the value in phenylalanine (Phe), β_Trc/Phe is the isotopic difference (3.4‰ for Glx and 3.1‰ for Pro) between the corresponding trophic AA and Phe in marine phytoplankton (Chikaraishi et al., 2009), and TEFs are empirical trophic enrichment factors between diet and consumer for Glx (Δ¹⁵N_Glx–Δ¹⁵N_Phe = 7.6 ‰ ± 1.2; Chikaraishi et al., 2009) and Pro (Δ¹⁵N_Pro - Δ¹⁵N_Phe = 4.5 ‰ ± 0.9; McMahon & McCarthy, 2016).

Therefore, to estimate TP we implemented Bayesian generalized multivariate multilevel models and statistical routines proposed in this study in the R package tRophicPosition v0.8.5 (Quezada-Romegialli et al., 2018), using Stan (Stan Development Team, 2023), rstan (Stan Development Team, 2023) and brms (Bürkner, 2017). To parse, compile and fit the models, we used Stan (Stan Development Team, 2022), a state-of-the-art platform for statistical modelling and high-performance Bayesian inference with Markov chain Monte Carlo (MCMC) sampling, as implemented in the R packages stan (Stan Development Team, 2023) and brms (Bürkner, 2017). In both bulk isotope and CSIA analyses, a Gaussian linear model with four chains, each with 2,000 iterations (1,000 burned), was implemented to calculate trophic position with no informative priors. Given that significant correlation was expected within hauls, this grouping factor was treated as a random variable, following a mixed effects approach (Bürkner, 2017). Both mean δ¹⁵N baseline and TEF were treated as normally distributed parameters, using the brms function ‘me’ (Bürkner, 2017). These bulk and CSI-AA models were implemented in the R package tRophicPosition v 0.8.5 (Quezada-Romegialli et al., 2018) in the functions stanOneBaseline() and stanOneBaselineCSIA(), respectively. This set of functions propagate the error term of β and TEF. The code can be accessed at https://github.com/clquezada/tRophicPosition/tree/stanCSIA.

The three trophic position estimation approaches were compared within and across consumers using medians, credibility intervals and the probability of differences derived from the posterior distributions from each model fit (Kruschke, 2010).

Bulk and individual amino acid δ15N by species

Among the different species analyzed, particulate organic matter (POM) exhibited the lowest mean δ¹⁵N_bulk value of −0.5‰ (±1.4), followed by E. superba, whose mean δ¹⁵N_bulk was 3.5‰ (±0.6). The highest δ¹⁵N_bulk values were observed for C. gunnari and E. antarctica, with mean δ15N_bulk values of 8.5‰ (±0.2) and 8.6‰ (±0.7), respectively (Table 1, Fig. 1). Regarding individual amino acid, nitrogen isotopic values were highly variable among species which ranged from more than 23‰ in glutamic acid (Glx) to -22‰ in threonine (Table 1, Fig. 1). Values of δ¹⁵N close to zero were observed in glycine, and serine, although differences of >3‰ were observed between species. Negative values were observed among all species for glycine and threonine. Overall, E. superba had lower values than the other species in all amino acids, except phenylalanine. Champsocephalus gunnari and Electrona antarctica had higher values for all amino acids than invertebrates (Table 1, Fig. 1), especially isoleucine, leucine, lysine, and valine (>2‰). Antarctomysis maxima and Themisto gaudichaudii had higher values than Euphausia superba and Euphausia frigida.

Electrona antarctica displayed the largest difference in Phe and Glx (ΔGlx-Phe) with a mean value of 23.1‰ (±1.0) followed by C. gunnari with 21.8‰ (±0.7) while the lowest difference was observed for E. superba (13.2 ± 2.2‰, Fig. 1). Euphausia superba also displayed the lowest difference between Phe and Pro (ΔPro-Phe) with a mean value of 9.7 (±2.2‰).Themisto gaudichaudii, E. frigida and A. maxima showed smaller differences (∼1‰) when using ΔPro-Phe than when using ΔGlx-Phe. Champsocephalus gunnari and E. antarctica values were slightly smaller when ΔPro-Phe (∼20‰) was used instead of ΔGlx-Phe (∼22‰).

Trophic position

Estimated TP varied depending on the specific method employed (Table 2, Fig. 2). TPs computed using bulk-SIA were similar to those computed using the Glx-Phe CSI-AA approach but showed more variability. Median TP_Glx−Phe estimates resulted in 0.01 levels higher, p(TP_Glx−Phe ≥ TP_bulk) =0.54 (i.e., probability of TP_Glx−Phe being higher than TP_bulk), than median bulk-SIA estimates, across all consumers which indicates that there was no difference between methods. TP_Pro−Phe estimates ranged between 0.2 and 0.8 levels lower than TP_bulk estimates, depending on the consumer, with a p(TP_Pro-Phe<TP_Bulk) =0.92 across all consumers. Similar differences were observed for C. gunnari TP estimates, whose median values reached 3.6, 3.4 and 3.1 (Table 2), using bulk-SIA, Glx-Phe and Pro-Phe CSI-AA , respectively.

Overall, there was greater disagreement between methods for invertebrates than for fish, although fish were consistently 1 TP higher than invertebrates. E. antarctica had a trophic level of 3.6 using bulk SIA, 3.5 using Glx-Phe and a TP of 3.3 when using Pro-Phe (Table 2). Unexpectedly, A. maxima appeared at a higher trophic position than the rest of invertebrates when using Glx-Phe (TP =3.2; Table 2) and its estimated TP values were approximately 0.4 higher than the rest of invertebrates. The lowest TP (1.9–2.2) was observed for E. superba (Table 2), while the trophic positions of T. gaudichaudii (2.1–2.6) and E. frigida (2.0–2-8) were similar regardless of the method (Table 2, Fig. 2).

Assessing relative accuracy in C. gunnari TP

The use of E. superba as baseline when estimating TP using bulk-SIA (TP_bulk) allowed us to assess precision of our TP estimates. The ΔGlx-Phe method overestimated the difference between C. gunnari and E.superba TP by 0.1 (CI_95%: 0.0−0.4) TP compared to the theoretical difference of 1 TP. The use of ΔPro-Phe yielded a 0.2 (CI_95%: 0.1−0.4) difference between the expected difference and the one calculated, while the use of bulk SIA overestimated C. gunnari TP by 0.4 (CI_95%: 0.1−1.0) (Fig. 3).