Stable Isotope Mixing Models Are Biased by the Choice of Sample Preservation and Pre-treatment: Implications for Studies of Aquatic Food Webs

Introduction

The development of stable isotope analysis (SIA) to study ratios of naturally occurring stable isotopes (SIs) of carbon (δ¹³C) and nitrogen (δ¹⁵N) has advanced the understanding of aquatic ecosystem functioning and in particular the knowledge about trophic interactions (Middelburg, 2014). Although SIA has been successfully applied in ecological studies for over three decades, many uncertainties about the effects of sample preservation and treatment still remain (Lau et al., 2012; Schlacher and Connolly, 2014).

In principal, SIA consists of the complete combustion of dehydrated and pulverized samples in an elemental analyzer and transferring the gaseous products to CO₂ and N₂, which are subsequently detected in a mass spectrometer allowing for distinguishing between molecules containing ¹²C and ¹³C, as well as ¹⁴N and ¹⁵N. Ship-based sampling and sampling in remote locations, however, seldom allow for immediate sample preparation and analysis. Accordingly, samples must be preserved until the time of analysis. In case of organisms and sediments it is commonly done by freezing, since fluid preservation in formalin and/or ethanol is known to affect SI ratios (Edwards et al., 2002). Such preservation effects are typically strong for δ¹³C (reported Δδ¹³C typically from −1 to −2‰), while they are small or insignificant for δ¹⁵N (reported Δδ¹⁵N typically from −0.5 to 0.5‰) (Lau et al., 2012; Rennie et al., 2012; Turner et al., 2015). These reports are in agreement with the chemical process of formalin preservation, which results in the formation of methylene bridges within or between protein molecules (Thavarajah et al., 2012). Since these bridges contain C originating from formalin, but no N, a change of δ¹³C toward the δ¹³C of formalin can be expected. Nonetheless, in an increasing number of studies formalin preserved organisms are measured for SIs, since it allows using community samples collected for other purposes or museum specimens (Chandra et al., 2005; Ozersky et al., 2012). Such studies typically use correction factors to adjust measured SI ratios (e.g., Chandra et al. (2005): addition of 1.28‰ to observed δ¹³C), although consensus about the appropriate correction has not been reached as the effect may vary with preservation duration, sample treatment, formalin brand, and studied taxa (Edwards et al., 2002; Rennie et al., 2012). Due to the general belief that freezing does not affect SI ratios (Bosley and Wainright, 1999), these correction factors have been established in most cases by comparison of formalin fixed and frozen samples without control for changes in SI ratios that are caused by freezing itself (Edwards et al., 2002; Lau et al., 2012). Some studies, however, found that freezing can have an effect on SI ratios of similar magnitude as preservation in formalin (Feuchtmayr and Grey, 2003; Syväranta et al., 2011). Furthermore, the reported effects of freezing on δ¹³C in those studies were opposing: according to Feuchtmayr and Grey (2003) freezing lowers δ¹³C values in zooplankton by approximately 1‰, while Syväranta et al. (2011) found higher δ¹³C (by approximately 2‰) values in the Asiatic clam (Corbicula fluminea) after freezing. For formalin preservation, however, both studies reported an enrichment in ¹³C in comparison to dried samples.

Another consideration in SIA is the use of sample pre-treatments to remove carbonates and/or lipids. When studying food webs, scientists are typically interested in organic (and not inorganic) carbon as the δ¹³C of an organism closely reflects the δ¹³C of its food (McCutchan et al., 2003). Inorganic carbon structures (e.g., carbonate shells) are strongly enriched in ¹³C and thus they have higher δ¹³C in comparison to organic tissues and are routinely removed from samples (Schlacher and Connolly, 2014). This is commonly achieved through treatment with hydrochloric acid (HCl). Such a treatment, however, has been reported to have variable effects on δ¹⁵N and therefore analyzing samples in duplicates for δ¹³C (acidified) and δ¹⁵N (non-acidified) is recommended (Schlacher and Connolly, 2014). However, since this increases the total number of analyzed samples many studies that investigate aquatic food webs use acidified samples for both SIs (e.g., Fredriksen, 2003; Kędra et al., 2012; Jankowska et al., 2018).

In contrast, the carbon in lipids originates from the diet of an animal. Nonetheless, chemical lipid removal is common when SIA is used to study aquatic food webs. The reasons behind this are the fast carbon turnover in lipids and their tendency to have lower δ¹³C in comparison to carbohydrates and proteins (Tieszen et al., 1983; Focken and Becker, 1998), which is likely introducing a strong bias in food web studies when samples with variable lipid content are analyzed. Consequently, chemical lipid removal is recommended when samples with presumably differing lipid content are being analyzed (Pinnegar and Polunin, 1999). Since chemical lipid removal can affect δ¹⁵N ratios (Søreide et al., 2006), this approach also requires the preparation of duplicate samples. This is, however, time consuming, increases the total number of analyzed samples, and it is often not possible in studies with small sample sizes (e.g., for benthic invertebrates). Consequently, many authors choose to use a mathematical lipid normalization of δ¹³C instead, based on the C/N ratio. Particularly, the normalization models of Post et al. (2007) and McConnaughey and McRoy (1979) are widely applied to aquatic invertebrates (Jardine et al., 2003; Coat et al., 2009; Arostegui et al., 2019). Lipid extraction and normalization in vertebrates, however, have shown that lipid normalization models should be tissue-, taxa-, and trophic group-specific (Cloyed et al., 2020) and accordingly it is questionable whether the abovementioned models should be applied so widely to all aquatic invertebrates.

While all above mentioned studies have addressed effects of sample preservation or pre-treatment separately, detailed knowledge about potential interaction effects of preservation and pre-treatment is still missing. Furthermore, our knowledge about preservation and pre-treatment effects on δ¹³C and δ¹⁵N ratios is based on univariate approaches that identify effects and potential corrections for each SI individually. This univariate approach to sample preservation and pre-treatment, however, is in stark contrast to the multivariate approaches that are commonly used to analyze SI data (Jackson et al., 2012; Parnell et al., 2013; Quezada-Romegialli et al., 2018). Particularly the use of Bayesian stable isotope mixing models (SIMMs) to quantify consumer diets have become widely applied (Parnell et al., 2013; Stock et al., 2018). Currently, Bayesian SIMMs provide a flexible modeling framework that allows to account for uncertainty in source data, covariance in tracer values in source and mixture data, and to incorporate categorical or continuous covariates (Parnell et al., 2013). This variety of mixing model parameterizations has been recently unified within the MixSIAR software package (Stock and Semmens, 2016a; Stock et al., 2018), making them easily available to a wide range of users.

So far, knowledge about whether guidelines regarding sample preservation and pre-treatments can introduce biases in SIMM results is missing. Furthermore, it is not known whether lipid normalization and formalin correction can correctly adjust SIMM results to compensate for the lack of chemical lipid removal or the use of formalin preserved samples. However, the outcomes of Bayesian SIMMs are strongly affected by mathematical lipid normalization when compared to lipid intact samples (Arostegui et al., 2019). Although Arostegui et al. (2019) demonstrated that the choice of lipid normalization can affect the result of mixing models, they did not assess whether chemical lipid removal would result in similar model outcomes as lipid normalization.

Therefore, the objectives of the present study were to: (i) quantify how sample preservation and pre-treatment affect carbon and nitrogen SI ratios (δ¹³C and δ¹⁵N), (ii) identify potential interaction effects among preservation and pre-treatment methods, (iii) study how preservation and pre-treatment affect the results of Bayesian mixing models, and (iv) assess whether lipid normalization and mathematical formalin correction should be used to adjust data for the use in such models. To address these objectives, we chose the Gulf of Gdansk (southern Baltic Sea) as study area, and two benthic invertebrate species (Limecola balthica and Crangon crangon) as study objects. We assumed that SI ratios are not affected by oven-drying, but that all other preservation and pre-treatment methods will change δ¹³C and/or δ¹⁵N. We aimed to determine whether an ecological interpretation of data is possible and how the uncertainties of the used methods should be addressed.

Materials and Methods

Study Location and Sample Collection

The sampling site in the Gulf of Gdansk (Figure 1) is characterized by sandy sediments, and year-round relatively stable bottom salinity around 7 (Witek et al., 2003). During times of high freshwater outflow from the close by Vistula river, however, the river plume can extend over the sampling site (Cyberska and Krzymiński, 1988). The vertical extension of the plume varies between 0.5 and 12 m (Grelowski and Wojewódzki, 1996). The bottom water temperature follows a typical seasonal pattern ranging from approximately 4°C in winter to 15°C in late summer. The onset of the spring phytoplankton bloom in the study region typically occurs between February and April, with the phytoplankton biomass peaking in April (Witek et al., 1997). Following the spring bloom, the primary production in the Gulf of Gdansk remains high until October (Witek et al., 1997). However, pelagic community respiration can exceed primary production for short periods in late summer, when water stratification is strong and water transparency is reduced (York et al., 2001). Such conditions have been reported to cause temporary oxygen depletion in bottom water of the Gulf of Gdansk (York et al., 2001).

Samples were collected from s/y Oceania on 24th April 2019 at the study location in the Gulf of Gdansk, southern Baltic Sea (54.40°N, 19.08°E, depth: 32 m; Figure 1). Two dominant benthic species, L. balthica and C. crangon, were chosen as study objects. Both species have a wide distribution in European waters and play an important role in the southern Baltic Sea ecosystem. The bivalve, L. balthica is a primary consumer that has the ability to switch between suspension- and deposit-feeding (Törnroos et al., 2015), while the shrimp C. crangon is a benthic predator (Oh et al., 2001). For each species enough individuals for 60 replicate samples (Figure 2A) were collected by use of a benthic dredge. Specimens were kept alive in filtered seawater for 24 h to allow for gut clearance. Altogether 60 individuals of L. balthica (size: 1–1.5 cm) and 120 individuals of C. crangon (length: 1–1.5 cm) were collected. For each species, the collected material was divided into 3 equal parts, which were subsequently preserved by oven-drying at 60°C on board, freezing at −80°C, and in 10% formalin buffered with borax (Figure 2A). Formalin (supplier Chempur) is an aqueous solution of formaldehyde (36–38%) stabilized with 10–15% methanol to slow down polymerization to paraformaldehyde. Consequently a 10% formalin solution contains approximately 4% formaldehyde and 1% methanol. L. balthica was removed from its shell prior to oven-drying and freezing. Intact individuals were, however, used for formalin preservation, since community samples and museum specimens are typically preserved in that way.

Furthermore, potential OM sources of the collected fauna were sampled. Twelve replicate samples of surface sediment (0–1 cm depth) and subsurface sediment (1–3 cm depth) were collected with box corer (sSOM – surface sediment OM; ssSOM – subsurface sediment OM). Half of the samples were dried at 60°C and the remainder frozen at −20°C (Figure 2C). A Niskin bottle was used to collect water from 10 m depth. The sampling depth was selected after a CTD profile showed a slightly elevated temperature (6°C) at the surface and a well-mixed water column with constant salinity (7.5) and temperature (4.5°C) below. A sampling depth well above the sea bottom within the well-mixed zone was chosen to reduce the risk of collecting resuspended sediment with the water. Water was then filtered through pre-combusted and pre-weighed Whatman GF/F filters to collect pelagic particulate OM (pPOM). Additionally, water samples were taken from Vistula River (54.25°N, 18.95°E, Figure 1) on 6th May 2019 to obtain riverine particulate OM (rPOM) samples following the same filtration procedures. Half of pPOM and rPOM samples (filters) were dried at 60°C and the other parts frozen at −20°C (Figure 2B).

Sample Preparation and Analysis

After 6 weeks of preservation, all frozen and formalin preserved samples were dehydrated. This was achieved by oven-drying at 60°C for pPOM and rPOM and freeze-drying for sediments and fauna samples. Dehydration by oven-drying and freeze-drying were considered interchangeable methods since multiple studies have found no effect of drying method on δ¹³C and δ¹⁵N of POM, SOM, and fauna samples (Bosley and Wainright, 1999; Carabel et al., 2006; Bessey and Vanderklift, 2014; Akamatsu et al., 2016), with sample types comparable to the sample types in this study. Accordingly, freeze-drying was used for sediments and fauna as it allows for easier homogenization of these sample types. Organisms preserved in formalin were soaked in water for 6 h. The water was exchanged twice during this time. L. balthica bodies were then removed from their shells prior to dehydration. All samples were then ground to a fine powder with mortar and pestle and kept in a desiccator until further processing.

In the next step lipids were removed from half of each of 3 sets (different preservation methods) of fauna samples (Figure 2A) by extraction in Folch solution [chloroform-methanol 2:1 (v/v)] for 24 h (Søreide et al., 2006). Consecutively extractant was exchanged twice to ensure the majority of lipid was removed. Samples were then air-dried under a fume hood at room temperature.

After weighing the samples into silver capsules, half of them were repeatedly (× 4) treated with 2 M HCl to obtain maximum of combinations for different preservation and pre-treatments methods (Figure 2). After each individual acid wash, samples were oven-dried at 60°C. Similarly, all pPOM and rPOM samples were subject to acidification. Samples for sediment OM were prepared in duplicates with and without acidification.

Samples were analyzed for carbon, nitrogen as well as their SIs (δ¹³C and δ¹⁵N) with an Elemental Analyzer Flash EA 1112 Series interfaced to an Isotopic Ratio Mass Spectrometer IRMS Delta V Advantage (Thermo Electron Corp., Germany). Isotopic ratios, δ¹³C and δ¹⁵N, were corrected against IAEA standards (CO-8 and USGS40 for δ¹³C and N-1 and USGS40 for δ¹⁵N). Laboratory working pure gases (CO₂ and N₂) of known isotopic composition were used as reference gases during the analysis. SIs measurements are reported in the commonly used delta (δ) notation in parts per thousand (‰) relative to the international standard Vienna Pee Dee Belemnite and atmospheric air for carbon and nitrogen, respectively.

Statistical Analysis

A two-way ANOVA was used to test for differences in isotopic ratios between OM sources (4 levels) and to identify potential preservation effects (2 levels) on the OM sources. To identify differences among fauna samples, a factorial ANOVA with four factors was used: (i) species identity (2 levels), (ii) preservation method (3 levels), (iii) lipid removal (2 levels), and (iv) carbonate removal (2 levels). Shapiro-Wilk test in combination with QQ-plots of the residuals and Levene’s test were used to test for normality and homogeneity of variance, respectively. Three OM samples (1 × rPOM_frozen; 1 × pPOM_dried; 1 × ssSOM_dried) were identified as extreme outliers that caused a violation of ANOVA assumptions for one of the isotopes. These samples were excluded from the analysis. Post hoc tests of significant factors and interactions of factors were done by Pairwise Comparisons of Estimated Marginal Means with Bonferroni correction applied.

Bayesian SIMMs were run in MixSIAR (Stock et al., 2018) to identify impacts of all combinations of preservation and pre-treatment on the estimated contribution of different OM sources to the diet of the collected fauna. Individual models were run per species-preservation combination, with pre-treatment as covariate (fixed factor with four levels, i.e., all combinations of lipid removal and acidification; Figure 2A). We included three distinct potential OM sources in the model: rPOM, pPOM, and SOM (i.e., surface + subsurface). However, as our results showed, SOM δ¹⁵N was significantly affected by freezing and therefore we ran every model twice with two different source data sets: (i) source data set A with dried SOM and (ii) source data set B with frozen SOM.

One key assumption that has to be made for Bayesian SIMMs is the selection of trophic enrichment factors (TEF) (Bond and Diamond, 2011). To ensure the best possible modeling results, TEFs should be selected as predator-prey pair specific as possible (Caut et al., 2009; Planas et al., 2020). Furthermore, a faithful incorporation of TEF variance is essential for model performance when literature-based TEFs are used (Stock and Semmens, 2016b). We assumed a trophic fractionation of 3.4 ± 1‰ for δ¹⁵N for each trophic step (Post, 2002). This commonly applied TEF was chosen in our study as it was developed largely on filter-feeding bivalves (pelagic baseline) and grazing snails (littoral baseline) and accordingly it is suitable for L. balthica and other primary consumers in the studied ecosystem. The same trophic fractionation is assumed for the second trophic step from primary consumers to predatory C. crangon, since it compares well with the fractionation of 3.6‰ that has been assumed by Fry (1988) in a community food web study that also included Crangon sp. Furthermore, it is close to an experimentally determined TEF of 3.55‰ for a mysid shrimp, Mysis mixta, that was fed on Artemia (Gorokhova and Hansson, 1999).

For δ¹³C, we used a scaled trophic fractionation, assuming a TEF of 4 ± 1.3‰ for the first trophic transfer from OM source to primary consumer, and of 0.8 ± 0.5‰ from primary consumers to secondary consumer, i.e., C. crangon. This scaled approach is necessary due to the particularly strong trophic fractionation of δ¹³C from OM sources to benthic primary consumers, which is consistently observed in marine systems (e.g., Fry and Sherr, 1984; Iken et al., 2010; Nerot et al., 2012). We selected a trophic fractionation of 0.8 ± 0.5‰ as the best estimate for the trophic step from primary consumers to predatory C. crangon, the same as reported for Crangon uritai in comparison to a variety of potential prey species (Antonio et al., 2011). The total TEFs for C. crangon in comparison to OM sources were calculated by summing of the TEFs for the trophic steps from OM source to primary consumer and from primary consumer to C. crangon. Since the studied species required the use of different TEFs (L. balthica: TEF_δ¹³C = 4 ± 1.3‰; TEF_δ¹⁵N = 3.4 ± 1.0‰; C. crangon: TEF_δ¹³C = 4.8 ± 1.8‰; TEF_δ¹⁵N = 6.8 ± 2.0‰), the analysis in MixSIAR was performed separately for each species.

All models were run on “extreme” setting (i.e., chain length = 3,000,000; burn = 1,500,000; thin = 500; chains = 3) with an uninformative prior and a multiplicative error term (Stock and Semmens, 2016b). Gelman-Rubin diagnostics and Geweke diagnostics were used to verify if models have converged. Furthermore, posterior distributions of the three chains were visually compared with each other to ensure that all chains produced almost identical results. For pairwise comparisons of posterior distributions within and between models Bhattacharyya’s coefficients (BCs) were calculated (Bhattacharyya, 1943; Rauber et al., 2008). Silverman’s rule of thumb was used to determine the bandwidths for the calculations of BCs (Silverman, 1986). BCs range from 0 (no overlap) to 1 (complete overlap).

The C:N ratio (weight/weight) was used to correct δ¹³C values of all fauna samples with lipids intact by subjecting them to two separate mathematical lipid normalization models (Figure 2A). For lipid normalization according to Post et al. (2007) all δ¹³C values from samples with C:N > 3.5 were normalized by the following formula:

For lipid normalization according to McConnaughey and McRoy (1979) all δ¹³C values from samples with C:N > 4 were subjected to:

Accordingly, samples with C:N < 3.5 and C:N < 4 were not affected by normalization according to Post et al. (2007) and McConnaughey and McRoy (1979), respectively. Two separate factorial ANOVAs (one for each normalization formula) with four factors were conducted: (i) species identity (2 levels), (ii) preservation method (3 levels), (iii) lipid approach (2 levels), and (iv) carbonate removal (2 levels). This allowed testing for differences in δ¹³C between lipid normalization and samples with chemically removed lipids, while also controlling for potential interaction effects with species identity, preservation method, and acidification.

Furthermore, δ¹³C of all formalin preserved samples was mathematically corrected by adding 1.28‰ to the original data according to Vander Zanden et al. (2003), which was previously used to study OM cycling in benthic food webs (e.g., Chandra et al., 2005). MixSIAR was then used to assess how the two lipid normalization formulas as well as the formalin correction affect SIMM results.

All statistical analyses were run in R version 3.6.2 (R Development Core Team, 2019), making use of emmeans (Lenth, 2020), MixSIAR (Stock and Semmens, 2016a), and dispRity (Guillerme, 2018) packages.

Since dehydration of samples is an essential step in SIA for all preservation methods, sample preservation through immediate drying was considered the least invasive approach. Consequently, we followed previous studies that considered oven-dried fresh samples as reference (e.g., Feuchtmayr and Grey, 2003; Syväranta et al., 2011). Accordingly, modeling results acquired from dried fauna samples with source data set A (dried SOM) were considered to reflect the diet composition most accurately and strong deviations of other preservation methods from this result were interpreted as poor performance. Pre-treatment effects on SIMMs are described in comparison to each other for each preservation method.

The effect of lipid normalization and formalin correction on SIMM results were evaluated by comparison to the results for original data and the desired outcome. Accordingly, lipid normalized SIMM results (for each preservation and HCl treatment combination) were compared to not normalized as well as lipid removed samples of the same preservation and HCl treatment combination. For formalin correction, model results for each pre-treatment combination of formalin corrected data was compared to results for each preservation method for the same pre-treatment combination.

Results

Preservation and Pre-treatment Effects

Organic Matter Sources

Allochthonous (rPOM) and autochthonous (pPOM, sSOM, ssSOM) OM sources were distinct from each other for both studied isotopes. Two-way ANOVA confirmed a significant difference between OM sources for δ¹³C [F(3,37) = 2062.89, p < 0.001] and δ¹⁵N [F(3,37) = 412.28, p < 0.001]. The following pairwise comparisons identified a significant difference for all pairs of OM sources, except surface and subsurface sediments, for δ¹³C and δ¹⁵N (Figure 3).

No effect of preservation method on δ¹³C of the OM sources was observed. For δ¹⁵N, however, a significant interaction effect between OM source and preservation method was identified by two way ANOVA [F(3,37) = 6.92, p < 0.001]. Pairwise comparison revealed that δ¹⁵N of sSOM and ssSOM was significantly lower in frozen samples in comparison to samples that were dried (Figure 3). Across all SOM samples, δ¹⁵N was on average 0.6‰ lower in frozen samples than in samples that were preserved by drying.

Fauna

Both studied isotopes showed a clear distinction between collected species (Figure 4). Across all preservation and pre-treatments C. crangon was on average enriched by 1.5 and 2.4‰ in ¹³C and ¹⁵N, respectively, in comparison to L. balthica. Overall, δ¹³C was more variable than δ¹⁵N. The total range of δ¹³C for each species for all treatments exceeded the difference between species and the δ¹³C range of the two species overlapped (Figure 4), which was not the case for the total range of δ¹⁵N.

In addition to this species effect, a significant effect of preservation and both pre-treatments on δ¹³C was identified by ANOVA (Table 1). Each of these effects was independent and no interaction effect was identified. Pairwise comparisons showed significantly lower δ¹³C values after formalin preservation in comparison to drying (p.adj < 0.001) and freezing (p.adj < 0.001). No significant difference between δ¹³C values of dried and frozen samples was found (p.adj = 0.36).

For δ¹⁵N, ANOVA identified significant interaction effects between species identity and preservation method as well as between species identity and carbonate removal (Table 1). The effect of both interactions on individual δ¹⁵N measurements were low. However, post hoc testing showed that sample acidification significantly affected δ¹⁵N in both species (p.adj_{L. balthica} < 0.001; p.adj_{C. crangon} = 0.02), while no effect of preservation method was identified as significant. The direction of change after acidification was opposite for both species. Acidification caused an average enrichment of ¹⁵N in L. balthica samples by 0.4‰, but a depletion (−0.3‰) in C. crangon samples.

Lipid Normalization

All samples with lipids intact (i.e., samples not treated with Folch solution) had C:N > 3.5 (Table 2) and accordingly all samples were normalized by Post et al. (2007) formula. Most C. crangon samples (24 out of 30) and one L. balthica (out of 30) had a C:N ratio between 3.5 and 4.0 and accordingly McConnaughey and McRoy (1979) lipid normalization was applied only to the remaining 6 C. crangon and 29 L. balthica samples.

Factorial ANOVA showed that δ¹³C values after lipid normalization differed significantly between the studied species, preservation methods, and acidification for both normalization formulas (Table 1). Furthermore, neither formula achieved a non-significant result for normalized data vs. chemical lipid removal (i.e., lipid approach in Table 1). In addition to the independent treatment effects, both normalization formulas introduced a significant interaction effect between species identity and lipid approach. This interaction effect was the result of each normalization formula performing well for only one species. Lipid normalization according to Post et al. (2007) performed well for C. crangon, but poor for L. balthica (Figure 5). More accurate lipid normalization of L. balthica was achieved by the formula of McConnaughey and McRoy (1979), which in turn gave worse results for C. crangon than Post et al. (2007).

Bayesian Mixing Models

Preservation and Pre-treatment Effects

Bayesian mixing models for L. balthica indicated a strong preservation effect on the model outcomes, which was further modified by pre-treatment methods (Figure 6 and Supplementary Material A1). For dried fauna pPOM was identified as main OM source across all pre-treatments for both OM source data sets (median BC = 0.01 for comparisons between pPOM and rPOM/SOM across all pre-treatments). However, the estimated pPOM contribution differed considerably between pre-treatments. pPOM contribution differed strongest between L. balthica that was treated with HCl only in comparison to samples that were treated with Folch solution only (BC = 0.55 and 0.57 for source data A and B, respectively). Estimated pPOM contributions were intermediate for samples that received no pre-treatment and samples that were treated with both: Folch solution and HCl. For frozen L. balthica samples the model results suggest a mixed diet of pPOM and SOM (Figure 6). A large overlap (BC > 0.65) between the probability distributions of pPOM and SOM was observed for all pre-treatments when source data A was used in the model. For source data B, the overlap between pPOM and SOM was smaller for L. balthica that was treated with Folch solution (BC = 0.56) and Folch solution and HCl (BC = 0.16). The estimated contribution of SOM was slightly higher for formalin preservation than for freezing. Nonetheless, pPOM and SOM contribution to formalin preserved L. balthica was overall similar to the models for frozen samples for both OM source data sets (Figure 6; median BC = 0.82). The estimated contribution of rPOM, however, (median diet proportion ∼0.2 across all pre-treatments, Figure 6) was higher for formalin than for dried and frozen samples. This difference in rPOM contribution between formalin and dried/frozen samples was particularly strong for samples that were treated with Folch solution (with and without HCl treatment) (median BC = 0.06). In contrast, rPOM probability distributions showed an intermediate overlap (median BC = 0.53) between formalin and dried/frozen L. balthica when no pre-treatment or only HCl was used. For L. balthica, the overlap between model results was large when results based on OM source data A were compared with OM source data B (median BC = 0.92 across all preservation and pre-treatments; Supplementary Table A29).

In contrast to L. balthica, model results for C. crangon indicated a strong reliance on SOM across all preservation methods when OM source data A was used (Figure 6). The used OM source data, however, had a considerable impact on the model outcome. Across all preservation and pre-treatments, source data set B resulted in an increased contribution of pPOM to C. crangon diet in comparison to source data set A (median BC = 0.77 across all preservation and pre-treatments; Supplementary Table A29). While only a weak preservation effect on the median model results was observed, the width of the probability distributions and accordingly the overlap between pPOM and SOM decreased from dried (median BC = 0.77) to frozen (median BC = 0.43) and further to formalin preserved samples (median BC = 0.01). The effect of sample pre-treatment on C. crangon mixing model results was overall lower than for L. balthica. Only sample pre-treatment with Folch solution and no HCl caused an increased estimated contribution of pPOM (Figure 6). The difference for this increase was low among preservation methods but increased when source data set B was used in the model. HCl treatment caused only minor changes to model outcomes. For samples that received both pre-treatments, however, HCl treatment seems to cancel out the effect of prior treatment with Folch solution.

Mathematically Corrected Data

Formalin correction resulted in a changed model outcome for L. balthica (Figure 6). Modeling results after formalin correction reflected almost exactly the results of dried samples for pre-treatments with Folch solution and Folch solution and HCl for both OM source data sets (median BC = 0.97). For no pre-treatment (median BC = 0.53) or pre-treatment with HCl only (median BC = 0.61), formalin correction made results more similar to dried samples, but the contribution of pPOM was underestimated and overestimated for no pre-treatment and HCl pre-treatment, respectively (Figure 6). Due to the high similarity between frozen and formalin posterior distributions (median BC = 0.68 across all pre-treatments for both OM sources), formalin correction caused an overall greater dissimilarity in comparison to frozen samples (median BC = 0.47 across all pre-treatments for both OM sources).

Formalin correction lead to an increased contribution of pPOM to the diet of C. crangon, making the results most similar to dried C. crangon (Figure 6). However, the median contribution of pPOM was overestimated on average by 0.15 and 0.17 for OM source data A and B, respectively. This strong effect was, however, accompanied by high BC values (median BC = 0.95 and 0.90 for OM source data A and B, respectively), due to the fact that formalin correction caused similarly wide probability distributions as dried samples. Since the median reliance of C. crangon on pPOM was lower for all preservation methods, results after formalin correction differed considerably from all preservation methods.

The effects of mathematical lipid normalization were quite variable (Figure 7 and Supplementary Material A). In contrast to the above described species-specific applicability of the two normalization methods, modeling results for both methods were similar for dried and frozen fauna (Supplementary Figures A5, A6, A12, A13). Only when lipid normalization was applied to formalin preserved samples, modeling results clearly differed between the two methods (Supplementary Figures A7, A14).

Modeling results of lipid normalized [according to McConnaughey and McRoy (1979)] dried L. balthica almost perfectly reflected model results of samples that had lipids removed with Folch solution (with and without HCl treatment) (Figure 7; median BC = 0.98). For frozen samples, however, this effect was only achieved for samples treated with HCl (median BC = 0.97). For frozen samples that were not treated with HCl, however, the median BC (0.81) was slightly lower than for the comparison between lipid intact and lipid removed samples (median BC = 0.84). Although these BCs indicate lipid normalized data and lipid intact modeling results being almost equally different from lipid removed samples, the model results between lipid intact and normalized data differed strongly (Figure 7), since the lipid intact model underestimated the contribution of pPOM in comparison to lipid removed samples, while lipid normalized data overestimated the contribution of pPOM. For lipid normalized formalin preserved L. balthica, normalization according to McConnaughey and McRoy (1979) performed better than normalization according to Post et al. (2007) (Supplementary Figures A7, A14). Posterior distributions of pPOM and SOM were improved after normalization according to McConnaughey and McRoy (1979). rPOM contribution estimates, however, were more different from lipid removed samples after lipid normalization.

Lipid normalization [according to Post et al. (2007)] had only a small effect on modeling results for C. crangon and results for lipid normalized data mirrored the results of the original, and not normalized, data (median BC = 0.95 across all preservation and pre-treatments for both OM source data sets). Accordingly, lipid normalization failed to reproduce the increased contribution of pPOM that was found for samples that were treated with Folch solution only (Figure 7).

Discussion

We report diverse effects of different preservation and pre-treatment methods on δ¹³C and δ¹⁵N ratios of the two studied species: L. balthica and C. crangon. The effects of preservation, lipid removal, and carbonate removal were relatively large for δ¹³C values, while the absolute effect on δ¹⁵N values was generally small. This is in agreement with previous studies that consequently highlighted that the methodological approach should be adjusted to optimize δ¹³C (Søreide et al., 2006; Lau et al., 2012). However, in contrast to previous studies, Bayesian SIMMs were used in this study to assess impacts on the results of data analyses. Modeling results showed that preservation and pre-treatment effects were variable and not directly linked to the strength of the absolute change of δ¹³C and δ¹⁵N signatures. This was highlighted particularly by the contrast between univariate analysis and MixSIAR for L. balthica. While the univariate analysis found no significant difference between dried and frozen L. balthica for both isotopes, the mixing model results suggested a different diet composition. In contrast, the univariate analysis found a significantly lower δ¹³C in formalin preserved samples in comparison to both, dried and frozen samples. Nonetheless, mixing model results of frozen and formalin preserved samples resembled each other closely. Consequently, small and not significant changes in SI ratios introduced by different preservation and pre-treatments (this study; Lau et al., 2012) display a so far unrecognized source of error in SIMMs. The observed preservation effects were strongest for formalin preserved samples, but also freezing, which is currently thought to be a less invasive preservation method, had a clear impact on modeling results. In addition to effects on fauna’s SI signatures, freezing also affected the δ¹⁵N values of SOM, which introduced another alteration of the mixing model result, but only for C. crangon. This effect of freezing is in strong contrast to the current believe that freezing should be used as the best method to preserve samples for SIA (Bosley and Wainright, 1999; Edwards et al., 2002).

Diet Composition of Limecola balthica and Crangon crangon

It is well documented that L. balthica is capable of suspension feeding as well as deposit feeding (Hummel, 1985; Törnroos et al., 2015) and accordingly a reliance on pPOM and SOM seem equally possible. It has, however, been shown that suspension feeding allows for deeper burial than deposit feeding, since the clam does not have to extend its inhalant siphon over the sediment surface (Hummel, 1985; Lin and Hines, 1994). Accordingly, suspension feeding provides better protection from predation and is the preferred feeding mode in situations with abundant food suspended in the water and epibenthic predators, like C. crangon, present (Kamermans and Huitema, 1994; Lin and Hines, 1994). Consequently, a high reliance of L. balthica on pPOM, as indicated by mixing model for dried samples, is likely, since sampling was conducted toward the end of the spring bloom when suspended OM was abundant in the sampling area.

Mixing model results suggest a strong reliance of C. crangon on SOM across all combinations of preservation and pre-treatment methods. The only exception were samples treated with Folch solution (but not with HCl) which showed a trend toward a considerable reliance on pPOM (Figure 6). Furthermore, the reliance of C. crangon on pPOM was around 0.15 higher when source data B (frozen sediment) was used as opposed to source data A (dried sediment). The overall high reliance of C. crangon on SOM is in agreement with previous studies that describe this species as a generalist benthic carnivore (Pihl and Rosenberg, 1984; Oh et al., 2001). For the smaller size classes of C. crangon (similar to the individuals in this study) it is likely that meiofauna constitutes a large proportion of their prey (Pihl and Rosenberg, 1984; Feller, 2006). Meiofaunal taxa are often specialized on consumption of detritus, bacteria or microphytobenthos (Carman and Fry, 2002) and consequently a high indirect reliance of C. crangon on SOM seems likely.

Organic matter from the Vistula river (rPOM) was not identified as a considerable OM source in the diet of both species. Sampling during the spring bloom (rich in marine production) and a long period with low river inflow during winter (low input of rPOM) support such a high reliance of fauna on pPOM and SOM. Nonetheless, it is important to point out that the low contribution of rPOM to the diet cannot be interpreted as an indication that the studied species do not consume rPOM. The reason for this interpretation is related to the fact that both pPOM and SOM that were collected at the study locations are potentially mixtures of autochthonous (marine OM) and allochthonous (rPOM) production (Maksymowska et al., 2000). The study location is influenced by the outflow of Vistula river and accordingly the OM collected at the study location cannot be considered entirely free of rPOM. Consequently, the strong contribution of SOM and pPOM to the diet of C. crangon and L. balthica, respectively, demonstrates that primary consumers do not select for rPOM. However, the preference of L. balthica is for suspended and possibly fresh OM while C. crangon preys on organisms that rely on SOM.

Preservation and Pre-treatment Effects

Mixing model results highlighted that frozen and formalin preserved L. balthica samples performed poorly, in comparison to dried samples. Formalin preservation performed worse than freezing, but the results suggest that both preservation methods should be avoided. Particularly the impact of freezing on modeling outcomes was unexpected, since samples are commonly frozen in other studies that use Bayesian SIMMs (e.g., Jankowska et al., 2018). Our modeling results for frozen acidified L. balthica with source data B (median estimated contribution of SOM: 29%) was similar to results in a previous study in the close-by Puck Lagoon that estimated 29 and 37% contribution of SOM to L. balthica diet during spring (Ziółkowska et al., 2018). Ziółkowska et al. (2018) preserved their samples by freezing and also acidified samples prior to SIA. Consequently, we suggest that the contribution of SOM might be overestimated in their study and the development of a mathematical correction for SIs for frozen samples, similar to the existing corrections for formalin preservation, might seem advisable. However, since freezing causes only small, taxon specific and often not significant changes in SI ratios (this study; Bosley and Wainright, 1999; Wolf et al., 2016) a possible development of any mathematical correction for frozen samples seems unlikely.

The potential mechanisms behind isotopic fractionation due to freezing are not well understood. Most commonly, authors who observed isotopic changes after freezing explain it by a potential damage to cells and tissues that can cause a loss of fluids that are potentially isotopically enriched (Feuchtmayr and Grey, 2003; Syväranta et al., 2011; Wolf et al., 2016). However, no study has measured whether fluids released through freezing and thawing have a different isotopic composition than the remainder of the material (Wolf et al., 2016). Freezing and formalin preservation also affected the mixing model outcome for C. crangon. The magnitude of change, however, was rather low and the choice of preservation method appeared unimportant for this species in the present study.

Since the effect of freezing on SIA is only poorly understood, it is worth re-considering whether freezing could be considered the better “reference” preservation method (de Lecea et al., 2011). Based on studies that suggested no difference between oven-drying and freeze-drying (Bosley and Wainright, 1999; Carabel et al., 2006; Bessey and Vanderklift, 2014; Akamatsu et al., 2016), we chose to freeze-dry frozen and formalin preserved fauna and frozen sediment samples to make the homogenization of the samples easier. In an experiment on the effects of freeze/thaw cycles on SI ratios de Lecea et al. (2011), however, found that oven-drying caused a slight increase of δ¹³C and δ¹⁵N in 6 fish species in comparison to freeze-drying. The authors suggested that this difference can be the result of the loss of volatile and isotopically lighter compounds, such as lipids, urea, or amino acids, during oven-drying. Based on this assumption we would expect higher δ¹³C and δ¹⁵N in dried fauna samples in comparison to frozen samples that were dehydrated by freeze-drying. However, no significant difference between dried and frozen fauna was found and it seems likely that only the freezing prior to oven-drying makes this loss of volatile substances possible. However, we observed a depletion in ¹⁵N for freeze-dried (frozen preservation) in comparison to oven-dried SOM, but not for ¹³C. Whether the loss of volatile compounds could cause such an enrichment in¹⁵N depends first of all on the presence of these compounds in SOM. Amino acids can be found in marine sediments, but in much lower concentrations than in fauna or POM (Wang et al., 2018). Accordingly, it is unlikely that a loss of lighter amino acids during oven-drying caused the observed difference in δ¹⁵N of SOM. Furthermore, individual amino acids differ in both δ¹³C and δ¹⁵N ratios (e.g., Pollierer et al., 2019) and we assume that the loss of specific amino acids should have caused a change in both isotopes. Consequently, the higher δ¹⁵N in dried SOM samples could indicate the breakdown of urea to ammonia and the consecutive loss of the gaseous ammonia (de Lecea et al., 2011). Urea is excreted by fauna and bacteria in sediments (Pedersen et al., 1993). However, urea has a short turnover time (Pedersen et al., 1993) and excretion of ammonia, and not urea, is typically dominating the excretion of excess nitrogen in benthic fauna (Allen and Garrett, 1971). Accordingly, it seems possible, but unlikely that considerable amounts of volatile amino acids or urea could have been present in our SOM samples and that they were lost during the oven-drying process.

In contrast to SOM, we did not observe a significant difference between dried and frozen fauna for either element. This difference in the effect of the frozen treatment between SOM and fauna cannot be explained by the dehydration method, but it could be related to the freezing temperature, since OM and fauna samples were frozen at −20 and −80°C, respectively. Specific fatty acids, proteins, and pigments degrade faster at −20°C than at −80°C (Wasmund et al., 2006; Nazemroaya et al., 2009), even though samples were only stored for 6 weeks. However, such a degradation does not affect SI ratios by itself. Only if there is a fractional loss or addition of C or N compounds from the samples the SI ratios might be changed. The biggest effect of freezing is the breaking of cells due to the freezing water and a following leak of compounds, which can be volatile or be part of water draining out of the samples. Furthermore, if degradation of pigments or proteins at −20°C should cause an enrichment in ¹⁵N, the same should have happened to POM samples. POM contains typically a higher fraction of fresh OM than SOM. Accordingly, the potential of degradation of labile compounds during freezing at −20°C to affect SI ratios is higher for POM and not for SOM. However, no difference in SI ratios of dried and frozen POM samples was observed.

Based on these considerations, we think that a potential loss of volatile substances during oven-drying or different freezing temperatures cannot explain the observed enrichment of SOM in ¹⁵N and accordingly agree with previous studies (e.g., Feuchtmayr and Grey, 2003; Syväranta et al., 2011) that considered fresh, oven-dried samples as the appropriate preservation method for control treatments.

Limecola balthica tissues were removed from their carbonate shells and the effect of carbonate removal on C. crangon isotopic composition has previously been shown to be non-significant (Jaschinski et al., 2008). Accordingly, HCl treatment has not been considered necessary for the fauna in this study. The effect of HCl treatment on δ¹³C, δ¹⁵N, and modeling outcomes was weaker but more variable than that of lipid removal. Nonetheless, ANOVA showed that both pre-treatments applied in our study had significant effects on at least one isotope for both species. Furthermore, modeling outcomes were affected by both pre-treatments with an interaction between the two pre-treatments for C. crangon. For L. balthica, the effect of treatment with Folch solution and HCl on SIMM results was approximately the additive effect of both individual pre-treatments. In contrast, HCl treatment did not affect mixing model results for C. crangon individually. However, in combination with Folch solution, HCl canceled out the strong effect of Folch solution on modeling results. This interaction of the pre-treatments for C. crangon has not been identified by ANOVA and depicts a strong discrepancy between univariate analysis and mixing model.

Mechanisms behind isotopic fractionation due to acidification are diverse (Schlacher and Connolly, 2014) and accordingly the causes for the isotopic ratios and modeling results in this study are difficult to address. However, the species-specific effect of acidification on δ¹⁵N indicate a potential role of chitin in the exoskeleton of crustaceans. Serrano et al. (2008) argued that HCl treatment of crustaceans may cause a shift of δ¹³C and δ¹⁵N through the known C and N loss due to hydrolytic depolymerization and deacetylation of chitin caused by the acidification (Percot et al., 2003). Although this has not been tested (Schlacher and Connolly, 2014), it provides a possible explanation for the species-specific acidification effects observed in this study.

Mathematically Corrected Data

In comparison to dried samples, ¹³C was depleted after formalin preservation on average by 1.1 and 1.4‰ for L. balthica and C. crangon, respectively. The depletion in comparison to frozen samples was slightly lower with 0.7 and 1.1‰, respectively. Consequently, a greater similarity to dried sample modeling results was expected after formalin preserved samples were corrected by addition of 1.28‰, since average δ¹³C after correction was within the range of dried samples. This assumption was, however, only met for L. balthica, while C. crangon modeling results overestimated the pPOM contribution after formalin correction. In contrast to the low reliance of C. crangon on pPOM for all preservation methods, the estimated reliance on pPOM increased after formalin correction (Figure 6). Since δ¹³C values were much closer to dried samples after the correction, the small and insignificant (according to post hoc test) species-specific effect of formalin preservation on δ¹⁵N values is likely responsible for the contrasting effect of formalin correction on the model outcomes. δ¹⁵N values of formalin preserved L. balthica were on average by less than 0.1‰ higher than for dried samples, and those of C. crangon by 0.3‰. These enrichments were comparable to the majority of formalin induced changes to δ¹⁵N reported in the literature (Lau et al., 2012). Our results for C. crangon, however, demonstrate that the correction of δ¹³C alone cannot compensate for a preservation effect on both SIs, even if the effect on δ¹⁵N is small and not significant. For L. balthica, on the other hand, with the preservation effect on δ¹⁵N being close to 0, the correction of δ¹³C achieved an improvement of modeling results by making them similar to results for dried samples.

The main physico-chemical processes of formalin preservation are the diffusion of methylene glycol and small amounts of non-hydrated formaldehyde into the tissue, followed by the covalent binding of the formaldehyde to some amino acids (Thavarajah et al., 2012). This reaction consumes the formaldehyde, shifting the balance between methylene glycol and formaldehyde, and causes more methylene glycol to dehydrate, making new formaldehyde available to bind to amino acids. Following this primary reaction of formaldehyde with an amino acid, a methylene bridge within or between protein molecules is formed in a secondary reaction over a long time (Thavarajah et al., 2012). This reaction is to some extend reversible by transfer of samples into water before sample preparation for SI analysis. However, considerable amounts of formaldehyde remain within the sample and accordingly a shift of δ¹³C toward the δ¹³C of the used formalin must be expected. Since formalin does not contain N, no shift in δ¹⁵N should be expected due to the described chemical processes of formalin preservation. Indeed Ogawa et al. (2013) demonstrated that formalin preservation does not affect δ¹⁵N of individual amino acids in gobiid fish. This observation was recently confirmed for δ¹⁵N of the amino acids of 5 fish species that were preserved in formalin and then transferred to ethanol (Chua et al., 2020).

Based on this lack of an effect of formalin fixation on the δ¹⁵N values of the molecules involved in the fixation process, one might argue that the sole focus on correcting the δ¹³C value for formalin preserved samples is justified. However, we argue that preservation in formalin should not be reduced to the chemical processes involved in the formation of methylene bridges between amino acids. It is well known that formalin fixation causes a strong loss of dry weight (range 10 to > 50%) in marine invertebrates (zooplankton and benthos) and fish larvae (Hay, 1984; Johnston and Mathias, 1993; Evgeny, 2003). This dry weight loss is a clear indication that formalin fixation cannot be reduced to the above described chemical processes since this would cause an increased dry weight. The dry weight loss is likely a combination of various co-occurring phenomena. Formalin has a high osmolarity (>1300 mO under standard conditions for 10% buffered formalin as used in this study). While this high osmolarity ensures the rapid diffusion of methylene glycol and formaldehyde into the tissue (Thavarajah et al., 2012), it also causes water to leave cells together with water soluble compounds, like sugars or salts. Further, formalin preservation has been suggested to cause loss of lipids (Jones, 1976), proteins and free amino acids (Hopkins, 1968; Williams and Robins, 1982; Arrington and Winemiller, 2002), or the leaching or fragmentation of mucosal membranes (Hay, 1984). Accordingly, all these compounds can be lost in the formalin solution.

No matter the cause of the dry weight loss, it constitutes a considerable loss of C and N. Accordingly, unless the dry weight is lost with the isotopic ratios of the entire sample, the effect of formalin fixation on SI ratios should not be seen as the result of the addition of isotopically different C from formalin only.

ANOVA results suggest that different lipid normalization formulas for δ¹³C should be selected for the two studied species. Currently, the common approach to lipid normalization of benthic invertebrates, however, is to apply one formula across all taxa when whole food webs are studied (e.g., Silberberger et al., 2018). This approach is very likely to introduce a considerable bias in such studies, considering how differently L. balthica and C. crangon are affected. Based on our results it has to be concluded that normalization according to Post et al. (2007) should be avoided for bivalves. The data used to develop this normalization formula was largely dominated by fish, but also included multiple arthropod species. In addition, the data contained a single sample of a unionid freshwater mussel. This single sample, however, seems to be an outlier in the data and can be perfectly corrected by the formula of McConnaughey and McRoy (1979) and (Post et al., 2007).

Many bivalve species, including L. balthica, can store energy in the form of glycogen (Wenne and Styczyńska-Jurewicz, 1987; Patterson and Carmichael, 2016). Patterson and Carmichael (2016) showed that high glycogen content in oysters results in C:N ratios greater than 3.5, even when lipid content is low. In contrast to lipids, glycogen is not depleted in ¹³C (Patterson and Carmichael, 2016). Accordingly, lipid extraction is not necessary and lipid normalization has to be avoided when glycogen is the major form of energy storage. L. balthica (among other bivalve species) is known to store energy in the form of lipids as well as glycogen (Wenne and Styczyńska-Jurewicz, 1987). It is therefore impossible to know whether high C:N ratios are an indication of high lipid or glycogen content. Furthermore, glycogen content can exceed the lipid content in other benthic invertebrate taxa [e.g., the polychete Hediste diversicolor (Durou et al., 2005)]. Consequently, we suggest that lipid normalization of δ¹³C based on C:N ratios should be abandoned for bivalves and other taxa with potential energy storage other than lipids.

The SIMM results after lipid normalization, however, were improved for L. balthica (exception frozen L. balthica with no pre-treatment), while hardly any effect was observed for C. crangon. Accordingly the presented results only confirmed to a limited extend the results from a previous study that found lipid normalization to affect mixing model outcomes in comparison to lipid intact samples (Arostegui et al., 2019). In contrast to Arostegui et al. (2019), however, the present study was able to compare lipid normalized also with lipid extracted samples. This demonstrated that the lipid normalization of C. crangon failed to correct for the effect lipid removal had on mixing models. As a conclusion at this point, lipid normalization is not advisable for benthic invertebrates when used in mixing models. With regard to the above described uncertainties about potentially high glycogen storage in benthic species and the failure of the lipid normalization to resemble lipid removal effects for C. crangon, chemical lipid removal should be used for benthic invertebrates. This statement, however, should be revisited when progress with species- and tissue-specific lipid normalization models for invertebrates has been made, comparable to the current state for vertebrates (Cloyed et al., 2020).

Methodological Considerations

As described in the methods, the TEFs in the current study have been selected as species specific as possible. While it is never entirely possible to verify TEFs, a comparison of the used TEFs with the average observed differences in stable isotope ratios between OM sources and the two studied species (Table 3) suggests that the chosen values and variance were appropriate. Nonetheless, the possibility remains that a wrong TEF was applied in the current study and as a consequence modeling results would be biased. However, since the same TEFs were applied to all models of the same species, we argue that all models are equally right or wrong, and accordingly the described preservation and pre-treatment effects are true differences.

Most studies that were used as reference for the selected TEFs used drying (Fry, 1988; Gorokhova and Hansson, 1999; Post, 2002; Antonio et al., 2011) or freezing (Iken et al., 2010; Nerot et al., 2012) for sample preservation. This agrees well with the chosen TEFs being closer to the observed trophic fractionation for dried and frozen fauna than to the observed trophic fractionation for formalin preserved samples. While this observation might suggest that TEFs should be developed for different preservation methods, we consider this a wrong approach as it mixes a natural process with a methodological error. Trophic enrichment of stable isotopes is a natural process that happened prior to sample collection and is therefore independent of any preservation or pre-treatment that is applied after sample collection.

The Bhattacharyya’s coefficient was used for pairwise comparisons between posterior probability distributions of SIMMs. In contrast to previous studies that used a significance level of 0.6 (i.e., BC > 0.6 indicates significant overlap), BCs were only used as a relative measure of overlap between probability distributions and patterns in median BCs were described. The decision for this use was made because the significance level of 0.6 has no clear underlying statistical or biological meaning. It was first introduced by Catry et al. (2009) for the Horn-Morisita Index. As pointed out by the authors, the value of 0.6 was arbitrarily selected and no objective justification for selecting this level was given. Bond and Diamond (2011) then introduced the same significance level for the BC to compare mixing model probability distributions, arguing it is analogous to other indices of dietary overlap (like the Horn-Morisita Index). However, if the significance level was chosen arbitrarily for the Horn-Morisita Index, then it is also arbitrary for BCs. Furthermore, while it is correct that both indices range from 0 (no overlap) to 1 (complete overlap), a clear difference between them exists. The Horn-Morisita Index measures actually observed dietary overlap between two samples. The BC, however, measures the overlap of two probability distributions and not actual dietary overlap. Accordingly, we suggest that BCs should be interpreted for what they are, a measure of overlap between two probability distributions. Consequently, a BC < 0.05 indicates significant different probability distributions and a BC > 0.95 indicates significant similar distributions. All values in between do not allow a statement about significance and should just be interpreted as probability of overlap.

Furthermore, like the Horn-Morisita Index is depending on the categorization of prey items (Diamond, 1983), the BC is strongly dependent on the binning/bandwidth used for its calculation. Previous studies (e.g., Jardine et al., 2015), however, did not report what binning was used and accordingly a fixed significance level is meaningless. In the current study Silverman’s rule of thumb was used to determine the bandwidths for the calculations of BCs (Silverman, 1986).

Conclusion

Our study confirmed the results of previous studies that the absolute preservation and pre-treatment effects on δ¹³C are high, while effects on δ¹⁵N are rather low. The strength of these effects is, however, not directly linked to their influence on the performance of Bayesian mixing models and consequently previously developed guidelines about sample preservation and pre-treatment need to be re-evaluated for this purpose. Furthermore, many observed effects on SI ratios and model performance were species-specific in this study and consequently every preservation and pre-treatment choice (for fauna and OM sources) may bias ecological interpretation. We assume that this has a potentially strong impact in the context of multispecies studies and suggest that future research should address how this might bias studies of whole food webs.

Furthermore, our study has demonstrated that mathematical correction of δ¹³C for high lipid content or formalin preservation produces variable mixing model results and should be avoided or at least be applied with extreme caution if data is used in Bayesian mixing models. Even if δ¹³C is corrected adequately, the model outcomes are often not improved and sometimes even worsened. Consequently, we recommend that the overall small, but apparently species-specific preservation and pre-treatment effects on δ¹⁵N should not be neglected in ecological studies that use SI ratios to study diet composition or food-web structure.

For future studies, we recommend that researchers (i) preserve samples by oven-drying, (ii) avoid lipid normalization, and (iii) prepare non-acidified duplicate samples for δ¹⁵N when data will be used in Bayesian SIMMs. We advise against the use of bulk δ¹³C and δ¹⁵N from formalin preserved museum specimens to reconstruct dietary composition with SIMMs, and suggest that isotopic niche parameters of bulk SI ratios (Le Bourg et al., 2020) or compound specific δ¹⁵N of amino acids (Ogawa et al., 2013; Chua et al., 2020) should be considered as these methods seem to be not sensitive to formalin preservation.

Furthermore, we call for increased efforts to improve our understanding on how our methodological choices in SIA affect our interpretation of the data, since the previous focus on quantifying the absolute changes of individual SI ratios does not allow for the ideal method selection.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

MS and KK-M carried out the fieldwork and performed the laboratory analyses. MS preformed the statistical analysis and drafted the manuscript. MK acquired funding that made this study possible. All authors designed the research, interpreted the data, provided comments during writing, and proofread the final manuscript.

Funding

This study was funded by Polish National Science Centre, project no. 2017/26/E/NZ8/00496 (COmEBACk). Sample collection was additionally supported through statutory funds of IOPAN. KK-M’s participation in this study was further supported by the Foundation for Polish Science (FNP).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We want to thank Captain and crew of s/y Oceania. Thanks to P. Makuch, M. Szczepanek, and Z. Borawska for help during sampling. Furthermore, we want to thank B. Szymczycha for guidance for chemical lipid removal and E. Jordà Molina for drawings of species used in the figures. A previous version of this manuscript appears online as a preprint at Biogeosciences Discussions (Silberberger et al., 2020).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2020.621978/full#supplementary-material

References