Bayesian estimation of predator diet composition from fatty acids and stable isotopes

A Bayesian mixing model for fatty acid profiles

Bayesian models for SI data are commonly based on the assumption that SI ratios are normally distributed. This assumption cannot be made for FA profiles, since for most methods of analysis, the concentration of individual FAs is normalized to the total lipid content of the sample. Thus, the FA profiles are a collection of proportions (referred to as a composition), which lie between 0 and 1, and are constrained to sum to 1. A common solution to this problem, however, is to consider transformations that make the data approximately normal (Budge, Iverson & Koopman, 2006). To construct our model, we considered the additive log ratio transformation (Aitchison & Bacon-Shone, 1999), also called alr transformation, such that (1)${\displaystyle y_{i,s}=alr\left({\varphi }_{i,s}\right)=log\left(\frac{{\varphi }_{i,s,1\dots p-1}}{{\varphi }_{i,s,p}}\right)}$ where ϕ_i,s is the p-variate fatty acid composition of individual i of prey species s, with a total of n potential prey species considered. Note that in the following we often drop the subscript for FAs, e.g., ϕ_i,s and y_i,s are thus p and p − 1 dimensional vectors, respectively. We assumed that the distribution of y is multivariate normal, with species specific mean μ_s and covariance matrix Σ_s, or y_i,s ∼ N(μ_s, Σ_s). A vaguely informative prior on μ_s and Σ_s allows for uncertainty in prey distributions to propagate to estimates of diet proportions (Ward, Semmens & Schindler, 2010).

Each prey source represents a proportion π_j of the diet of predator j, and analogous to stable isotope mixing models, predator FA profiles are then a linear combination of prey FA profiles, normalized to sum to one. Since predators may selectively assimilate or metabolize FAs (Iverson et al., 2004; Budge, Iverson & Koopman, 2006; Rosen & Tollit, 2012), we specified prey-specific conversion coefficients κ_s = κ_s,1…κ_s,p for each of the p FAs (Rosen & Tollit, 2012). Furthermore, the n prey species likely have different fat content Φ that will affect the total amount of FAs assimilated from each prey species by the predator. The expected FA profile of predator τ_j is then a linear combination of the prey FA profiles, modified by conversion coefficients for each fatty acid p and fat content for each prey i: (2)${\displaystyle t_j\sim N\left(alr\left({\tau }_j\right),{\Sigma }_{\tau }\right)}$ (3)${\displaystyle {\tau }_j=C\left\{\sum _s^n\left({\pi }_{j,s}{\Phi }_s\right)\left({\kappa }_s\otimes {\varphi }_{j,s}\right)\right\}.}$

Here, C is the closure operation which normalizes the FA profiles to sum to one and ⊗ is the outer (element wise) product. ϕ_s,j is the FA profile of prey items of species s consumed by predator j. Similarly to Parnell et al. (2013), we thus assumed that individual predators do not necessarily feed on ‘average’ prey items, but rather consume prey items with signatures drawn from the estimated prey distribution. We formulated predator signatures t as draws from a normal distribution after transformation. We further assumed that Φ and κ are log-normally and gamma distributed, respectively, around known mean and variance values (estimated or calculated from controlled diet experiments, see below). The closure operation in Eq. (3) (i.e., the sum-to-one constraint on the FA profiles) leads to κ being determined in terms of relative uptake of FAs (i.e., up to a multiplicative constant), and implicitly makes the multivariate prior distribution over all κ a Dirichlet distribution. The same logic applies to Φ, and in both cases we opted for formulations that can be readily parametrised from sample means and variances from controlled diet experiments.

The diet proportions π of predators are the main focus of investigation in diet studies. These may be modeled at the population level (thus dropping the subscript j in expressions (2) and (3)) or at the individual level, as suggested in expressions (2) and (3). In the latter case, individual predator FA profiles can be modeled as random samples from a population level distribution of predator diet proportions. Recent approaches to stable isotope mixing have focused on transformations of the diet proportion vector π to get around the problems associated with the compositional nature of the diet proportions in such a hierarchical setup, and we followed this approach in our model. The diet proportions were transformed using centered log-ratio (clr) transformations (Semmens et al., 2009), such that the support is the real line rather than the interval [0;1], and we then assumed that clr(π_j) ∼ N(Π, Σ_Π), where Π is the vector of mean (population level) diet proportions. It is then possible to model diet proportions as a function of covariates, such as size, sex, or region (i.e., in a regression formulation, Parnell et al., 2013). While this approach is appealing, it adds to computation time needed to estimate model parameters, and often results in slower convergence. We therefore use a vague Dirichlet prior on the proportions when convenient (i.e., when we estimate only population level parameters).

Depending on the number of samples for prey and predators, it may be necessary to use informative priors for Σ_s and Σ_τ. Both were given inverse-Wishart priors, and since both are co-variances of transformed data, it is not straightforward to formulate default priors for these parameters. We have found that, in practice, manual adjustment of these priors is often needed to be able to achieve convergence and mixing (efficient exploration of the posterior distribution by the sampling algorithm) of the Markov Chain Monte Carlo (MCMC) routine employed to estimate model parameters. This is especially true when there are few source and/or predator samples. Our code (see below) allows for high level adjustment of these parameters through the specification of the order of magnitude of the diagonal of each covariance matrix.

Joint diet estimation from FA profiles and SI

There are at least three potential benefits of integrating FA profiles and SI data: (i) increased information to discriminate among sources, (ii) the potential of SI to resolve predator prey relationships due to trophic enrichment of SI, and (iii) the potential reduction in estimation error due to a larger body of research on fractionation coefficients for stable isotopes as opposed to conversion coefficients in FA profiles. Integrating the two complimentary types of data in a single model to estimate diet proportions may thus considerably improve estimates of diet proportions over estimation from either data-source alone.

Our model for FA profiles is conceptually similar to recent models proposed for SI data, and integration of FA profiles and SI data into a single model is straightforward in the present setting. We assumed that the vector of SI signatures of sampled prey items q follow a multivariate normal distribution, such that $y_{q,s}^{SI}\sim N\left({\mu }_s^{SI},{\Sigma }_s^{SI}\right)$, where the superscript SI denotes that these are stable isotope signatures. Predator SI signatures are again a linear combination of prey SI, this time modified by additive fractionation coefficients γ. Fractionation may, in turn, depend on prey isotope concentrations (Hussey et al., 2014; Caut, Angulo & Courchamp, 2009). In our model, we assumed additive fractionation, and suggest that concentration dependence is taken into account when specifying distributions for prey and SI specific fractionation coefficients γ_s (see examples below). The expected SI signature for predator r is then (4)${\displaystyle t_r^{SI}=\sum _s^n{\pi }_{r,s}\left(y_{q,r}+{\gamma }_s\right)}$ (5)${\displaystyle clr\left({\pi }_r\right)\sim N\left(\Pi ,{\Sigma }_{\Pi }\right)}$ (6)${\displaystyle {\gamma }_{s,SI}\sim N\left({\nu }_{SI},{\sigma }_{SI}\right).}$

Note that the different subscripts to the FA profile model imply that there is no need to have SI and FA profiles from the same prey or predator samples, as long as we can assume that the prey samples are drawn from the same statistical population as those for FA profiles, and that individual diet proportions of predators are drawn from the same population distribution of diet proportions.

The exact formulation of the integration of SI and FA profiles depends on the assumptions that one is comfortable with in a given setting: identical dietary proportions may be appropriate if diets (and hence SI and FA profiles) are thought to be stable, or if both chemical tracers are thought to integrate over similar time-scales. If the time scales of these two elements are thought to be different (e.g., for different tissue types), individual diet proportions may be more appropriate, and may be drawn from an overall population distribution of diet proportions.

An R (R Core Team, 2014) package (called fastinR) implementing methods outlined here, along with simulated examples and the analysis of experimental data described further below, is available on the open source repository github.com/philipp-neubauer/fastinR. Models implemented in the package include the above-mentioned formulations for population level diet estimates, individual diet estimates as well as linear model (regression and ANOVA) formulations for diet proportions, all available for SI and FA profiles individually or as combined models (see below). Model parameters were estimated using MCMC methods implemented in JAGS (Plummer, 2003), called from R through higher level functions in the fastinR package that allow for data input, inspection and manipulation.

Simulation studies

We initially explored the feasibility and performance of our model setup in a range of simulations, which are illustrated (including code) in Supplemental Information 1. The initial feasibility simulations used a set of three potantial prey species (30 samples per species) with two stable isotopes (i.e., an under-determined system) and 12 fatty acids. To clearly illustrate the method itself, prey-source separation was chosen to readily discriminate prey items for both markers, but with enough variability in prey profiles to illustrate uncertainty propagating through to diet estimates.

Simulations were also used to explore sensitivities of inferred diet proportions to the source configuration and diet evenness in a series of simulation experiments, using 10 simulated fatty acids and four potential prey items. We hypothesized that estimated diet proportions are sensitive to diet source separation in FA profile space, co-linearity in FA profile space (Blanchard, 2011) and diet makeup (e.g., specialist versus generalist diets). Further details and simulation results can be found in Supplemental Information 2.

Selecting fatty acids for analysis: an ordination approach

A potentially large number of FAs are available from analysis methods such as gas-chromatography. A common practice is to simply set a threshold and keep the most abundant FAs for analysis. This practice may, however, discard potential useful information, and a more judicious approach is to retain FAs based on the among diet source variability that they explain. Wang, Hollmen & Iverson (2010) used a method by which they tested the QFASA method on a series of subsets to determine the subset that gave the best accuracy. Although feasible, such a method can be time consuming with fully Bayesian models, which can take a long time to run with a large dataset.

Here, we propose a variable selection method based on constrained ordination, which considers the contribution of individual fatty acids to axes separating diet sources. Based on this contribution relative to the overall separation, the user can choose FAs that contribute most to source separation. This procedure is intended to reduce computation time (and dimensionality) of the models, while retaining as much accuracy in diet estimates as possible. Further details about the procedure are given in Supplemental Information 3.

Estimating predator diets in a controlled experiment

To illustrate the potential of the models presented above, we analysed data from an experimental study by Stowasser et al. (2006), which investigated changes in squid FA profiles and SI as a function of diet treatments. The treatments consisted of exclusive fish and crustacean diets, as well as switched and mixed diets, with the former switching diets from fish (henceforth SF, n = 4) to crustacean (SC, n = 5) after 15 days of the 30 day experiment.

In order to apply our model, we first estimated conversion coefficients of FA profiles and fractionation in SI, using squid from the 30 day diet treatments feeding exclusively crustacean and fish diets. The model for estimation of SI fractionation followed the model in Hussey et al. (2014), thus accounting for diet δ¹⁵N and δ¹³C, and used their results as priors for fractionation parameters for δ¹⁵N, and results from Caut, Angulo & Courchamp (2009) to construct priors for δ¹³C. Estimation of FA conversion coefficients used expressions (2) and (3) with proportions assumed known from feeding trials. Computational details on the estimation of conversion coefficients and fractionation are given in Supplemental Information 4.

In our diet analysis, we analyzed samples from the switched diet treatments, and used both SI and FA profiles to investigate the ability of our models to accurately estimate diet proportions in either treatments. We subset the data to use only squid from the switched diet experiment that were analysed for FA profiles and SI after at least 10 days under the respective treatment. We only had overlapping SI and FA profiles for the SC treatment squid, and we therefore started by analyzing this treatment in isolation to demonstrate that both SI and FA profiles can resolve diet proportions, and to demonstrate the benefit of using the two tracers in a joint model. We then analyzed the SF treatment squid, for which we only had 3 specimens with FA profiles and 1 specimen with SI. The markers available for this treatment did not overlap for any of the sampled squid.

Lastly, we estimated individual diet proportions in the SC treatment. To demonstrate how the model based approach to diet estimation can be use to answer ecologically relevant questions about predator diets, we also analyzed SF and SC treatment squid together in a linear model setup that investigated treatment differences explicitly. The linear model used treatment dummy variables to estimate individual intercepts for each treatment and prey combination, and allowed us to test for significant differences in diet between treatment groups, conditional on the data and priors.

FA profile analyses used data obtained by analyzing digestive gland tissue, which is thought to rapidly assimilate dietary FAs in relatively unmodified proportions relative to the original diet (e.g., Phillips, Jackson & Nichols, 2001). SI were analyzed from muscle tissue since we had more individuals sampled for SI from this tissue, which may be more prone to fractionation and slower turnover than digestive glad tissue. In the original study, a total of 25 FAs were reported. Here, we selected FAs using ordination methods described above. For estimation of model parameters, priors for prey and predator specific variances were adjusted manually to give reasonable behaviour in the MCMC algorithm. The analyses are detailed in Supplemental Information 5.

Simulation studies

Simulated test cases suggested that our model can estimate diet proportions from both SI and FA profiles (Supplemental Information 1), with accuracy depending mainly on source separation and diet evenness (Supplemental Information 2). For very uneven diet proportions, such as in the feeding trials analyzed in the squid example, we found the choice of posterior means as point estimate for diet proportions inevitably introduced error at the margins of the 0–1 interval when compared to true simulated diet proportions.

Models with low accuracy conversion coefficients (with prior mean for all FA set to 1 and large prior variance) also performed substantially worse than models with accurately specified coefficients when comparing point estimates of diet proportions to simulated diet proportions (Supplemental Information 2), showing decreasing accuracy with increasing variance among simulated convergence coefficients.

Estimating predator diets in a controlled experiment

Dimension reduction by NMDS on FA profiles of squid and their potential prey suggested that crustacean diets were readily distinguishable from fish diets (Fig. 1A). For fish diet items, however, no single fish species could be clearly distinguished from any other fish species. Predator signatures of switched diet squid aligned with their respective diets after correcting by posterior means of estimated conversion coefficients. The latter were different from expected (1/p) for many FA in the analysis (Supplemental Information 4).

Selection of FAs using constrained ordination lead to four FAs, 22.6n.3, 20.5n.3, 20.4n.6 and 18.1n.9 being retained for analysis (Fig. 2), accounting for a total of 74% of total among source variation on ordination axes while maintaining a low prey matrix condition number (κ = 15.67), suggesting limited co-linearity. The matrix condition number nearly doubled when the next most important fatty acid (κ = 29.17) was added and increased exponentially thereafter with the addition of other FAs. The resulting NMDS plot suggested that the reduction from 22 to four FA did not significantly alter the configuration of predators and prey items in FA profile space, despite the drastically lowered number of input dimensions (Fig. 1B). Retaining a larger subset of FAs (8 FAs) did not qualitatively alter the results, but did lead to lower uncertainty in diet proportion estimates, suggesting that we lost some relevant information by retaining only four of 25 original FAs to reduce computational requirements.

SI also showed clear separation between crustacean and fish prey (Fig. 3), but showed two groups for fish prey items, both consisting of specimens from more than one fish species. Squid δ¹⁵N was also substantially lower than any of the prey species analysed even after correcting for estimated fractionation coefficients.

FA profiles were able to resolve population level SC treatment squid diets, correctly suggesting a diet predominantly based on crustaceans (Fig. 4). While uncertainty about the exact diet proportions remained for both crustaceans and fish, most of the posterior density for squid diet proportions was clearly concentrated towards high proportions of crustaceans. For fish, posteriors were peaked near zero, however, all fish species posteriors had long tails that spanned nearly the whole interval of possible diet contributions. An analyses based on SI alone gave very similar results, despite different tissue types examined (Fig. 4).

Combining the two markers lead to a substantial reduction in the uncertainty of estimated diet proportions (Fig. 4), and suggested a clear dominance of crustaceans in the diet. For the combined analysis, the spread of the posterior distribution for crustaceans in the squid diet was reduced by approximately 30%, and most of the probability density was shifted closer to one, and the reductions in the spread of posterior distributions for fish diet items were as high as 70%. Lastly, estimates of individual diet proportions closely mirrored population level estimates (Fig. 5).

Due to overlap of fish species in FA profile and SI space, similar models for SF treatment fish were unclear about the contribution of individual fish species (Fig. 6), but suggested that crustaceans were a small part in the diet of these squid. SI and FA profiles combined (i.e., adding one squid with SI but no FA profile data) did not provide much improvement for individual fish species, and the linear model setup was unable to identify significant differences between diet proportions of individual prey items in the two treatments due to uncertainty about individual fish species’ contributions (Supplemental Information 5). However, combining fish species post-hoc as the sum of individual posterior distributions clearly shows a fish based diet (Fig. 7).