Integrating disparate datasets to model the functional response of a marine predator: A case study of harbour porpoises in the southern North Sea

Abstract Quantifying consumption and prey choice for marine predator species is key to understanding their interaction with prey species, fisheries, and the ecosystem as a whole. However, parameterizing a functional response for large predators can be challenging because of the difficulty in obtaining the required data on predator diet and on the availability of multiple prey species. This study modeled a multi‐species functional response (MSFR) to describe the relationship between consumption by harbour porpoises (Phocoena phocoena) and the availability of multiple prey species in the southern North Sea. Bayesian methodology was employed to estimate MSFR parameters and to incorporate uncertainties in diet and prey availability estimates. Prey consumption was estimated from stomach content data from stranded harbour porpoises. Prey availability to harbour porpoises was estimated based on the spatial overlap between prey distributions, estimated from fish survey data, and porpoise foraging range in the days prior to stranding predicted from telemetry data. Results indicated a preference for sandeels in the study area. Prey switching behavior (change in preference dependent on prey abundance) was confirmed by the favored type III functional response model. Variation in the size of the foraging range (estimated area where harbour porpoises could have foraged prior to stranding) did not alter the overall pattern of the results or conclusions. Integrating datasets on prey consumption from strandings, predator foraging distribution using telemetry, and prey availability from fish surveys into the modeling approach provides a methodological framework that may be appropriate for fitting MSFRs for other predators.


| INTRODUC TI ON
Prey populations are directly and indirectly affected by predation and their dynamics are influenced by long-term and short-term responses of predators (Holling, 1959;Murdoch & Oaten, 1975). The functional response helps to assess the potential impact that predators could have on their prey by describing the response of predator consumption rates to varying prey densities, providing insight into prey preference and general predator-prey interactions (Dale et al., 1994). High consumption rates indicate strong interactions between predators and prey, resulting from high encounter rates and/or active predator choice. Switching between prey species may occur if predator preference changes with prey density, for example, when predators avoid scarce prey (Holling, 1959).
Although the functional response has been subject to extensive empirical research, most studies have been conducted within a laboratory setting or have described relationships among a small number of species (Morozov & Petrovskii, 2013). Modeling the multi-species functional responses (MSFR) for wild animals is challenging because observing both consumption and prey availability outside a controlled environment is difficult. Parametrizing a MSFR requires substantial datasets on predator diet and distribution, and on the availability of multiple prey species covering a range of prey densities. It is not surprising, therefore, that the ecological role of most large predators has not been quantified and that we have an incomplete picture of their impacts in many ecosystems (Estes et al., 2011).
However, the use of Bayesian methods can overcome the problem of data sparsity, allowing MSFR models to be fitted for top predators (Smout et al., 2014;Suryawanshi et al., 2017).
The aim of this study is twofold. Firstly, we develop a framework to integrate long-term datasets on predator consumption, predator distribution, and prey abundance to model the MSFR of a marine high trophic level predator. The framework consists of a number of methodological steps for modeling changes in diet in relation to prey abundance, which are appropriate for mobile marine predators for which diet is sampled at specific locations. Our intention is that this methodological framework can serve as a model for other similar studies and thus help improve understanding of the ecological role of high trophic level marine predators. We develop and apply this framework using the harbour porpoise (Phocoena phocoena) in the southern North Sea as a case study to examine the methodology, model performance, model output, and the sensitivity of the results to variation in assumptions. We choose the harbour porpoise partly because there are data on prey consumption from the stomach contents of stranded porpoises in the Netherlands (Leopold, 2015), data on the distribution and movements of individual porpoises in the North Sea from satellite-linked telemetry (Sveegaard et al., 2011), and data on prey abundance from the ICES International Bottom Trawl Surveys (ICES, 2018).
Secondly, in choosing the harbour porpoise as a case study, we aim to improve ecological understanding of an important marine predator in European Atlantic waters. The harbour porpoise is the most abundant large marine predator in the North Sea (Hammond et al., 2013) and its diet includes species that are also targeted by commercial fisheries (Santos & Pierce, 2003), such as whiting (Merlangius merlangus), Atlantic herring (Clupea harengus), and sandeels (Ammodytidae).
Harbour porpoises have a high metabolic rate and only a limited energy storage capacity, which limits their ability to buffer against diminished food availability/quality and makes them more susceptible to starvation if they fail to meet their high metabolic demands (Rojano-Doñate et al., 2018;Spitz et al., 2012). They have high ingestion rates and probably must consume prey on a daily basis (Kastelein et al., 2019;Wisniewska et al., 2016), unlike other cetaceans or pinnipeds that might move through certain areas while not foraging, and thus are particularly appropriate for this study. North Sea wide surveys showed a major north to south shift in the summer distribution of harbour porpoise from 1994 to 2005, maintained through 2016, which was likely linked to changes in prey distribution (Hammond et al., 2002(Hammond et al., , 2013. Information on the dynamic relationship between harbour porpoises and their prey is largely lacking but improving understanding of harbour porpoise predator-prey relationships may help to explain the observed shift in distribution.

| MATERIAL AND ME THODS
As a framework for analysis, the following sequence of steps (described in detail below) was followed to parameterize the functional response: (1) estimation of diet composition; (2) estimation of foraging range; (3) estimation of prey availability; and (4) fitting the multispecies functional response. All data processing and modeling were performed in software R (R Development Core Team, 2018), and MSFR fitting was carried out in WinBUGS (Lunn et al., 2000).

| Diet composition
To obtain information on harbour porpoise prey consumption, diet composition was estimated from the hard remains of prey (fish otoliths) recovered from the stomachs of individual animals stranded along the Dutch coastline between 2006 and 2015. To match diet composition to availability of prey data strandings that occurred in November-April were assigned to quarter 1 and those in May-October to quarter 3 (see Section 2.1.3). Sample collection and analysis are described in Leopold (2015). Postmortem examinations were carried out on stranded animals documenting standard measurements (e.g., body length). Prey species were identified to the lowest possible taxon. Otoliths were measured, paired when possible, and graded for wear. Grade-specific correction factors were used to estimate undigested otolith size and prey weight was estimated by applying otolith size-fish mass relationships. Prey species that contributed ≥5% of the total estimated prey weight were selected as main prey species.
Uncertainty in diet composition arises from measurement (estimation of prey weight) and sampling error (Hammond & Rothery, 1996). Sampling error was estimated by nonparametric bootstrapping using individual stranded porpoises as the sampling unit, stratified by season. To balance carcass freshness and retain an adequate sample size, only individuals with decomposition codes less than 4 were included in analysis (see Leopold, 2015). Measurement error was not estimated.

| Foraging range
Foraging range was defined as the geographical range in which a stranded porpoise could have foraged. Note that this is different from the realized foraging area, which includes a component of predator "choice" regarding prey availability, which we want to avoid.
Estimating the foraging range of porpoises prior to stranding is difficult due to the unknown location of death. It is possible that a stranded porpoise was alive and swimming until just before it stranded, or carcasses could have drifted at sea for a considerable period of time (Peltier et al., 2013). This introduces uncertainty in defining the area where porpoises likely foraged. We used information on the rate at which porpoises could have moved prior to stranding to obtain informed estimates of their potential foraging range.
The foraging range was estimated using telemetry data from satellite-linked tags deployed on harbour porpoises in the Kattegat, Belt Seas, and Western Baltic between 1997 and 2015 (see Teilmann et al. (2007) and Sveegaard et al. (2011) for tagging procedures, tag settings, and data filtering). The movements of harbour porpoises in the Kattegat and Belt Seas differ from those further north in the Skagerrak and in the North Sea (Sveegaard et al., 2011). To ensure the data were as representative as possible for porpoises that stranded in the southern North Sea, data from the southern Kattegat and further south (south of latitude 57.30°N and east of longitude 9.37°E) were excluded.
The use of stomach content data to estimate prey consumption depends on knowledge of the temporal window within which porpoises could have obtained their last meal, which is dependent on how long prey remains stay in the stomach. In the absence of information on passage rates of hard prey remains for harbour porpoises, information for similar sized gray seals Halichoerus grypus and harbour seals Phoca vitulina, which consume similar prey species, was used. Two days after consumption >50% of all otoliths were recovered in gray seal (Grellier & Hammond, 2006) and >85% in harbour seals scats (Wilson et al., 2017). To estimate harbour porpoise foraging range, a minimum timeframe of 2 days was chosen. Additionally, timeframes of 4, 6, and 8 days were applied to explore how resilient the results were to variation in the likely foraging area, including to accommodate any drifting of carcasses postmortem.
Prior to modeling the telemetry data, the track line of each tagged porpoise was processed to create positions at regular intervals. These positions were used to generate minimal enclosing circles (MECs) from sets of consecutive points for timeframes of 2, 4, 6, or 8 days (Figure 1). Using a generalized linear model (GLM), the MEC diameter (response variable assumed to follow a gamma distribution with log link) was modeled as a function of timeframe and age, sex, season (quarter of the year), and all two-way interactions.
Model selection was based on AIC scores. The variance inflation factor (VIF) was used to detect multicollinearity using a threshold of 4 (Hair et al., 2010).
Tagged individuals are measured repeatedly, so a generalized linear mixed model (GLMM) including a random effect for individual was also investigated. However, the GLM was better supported than the GLMM according to AIC scores and log-likelihoods.
Stranded porpoises are located on the coast, so the diameter of the MEC estimated from the GLM was used to predict the radius of a circular buffer, centered on stranding location, to approximate the foraging range (at sea) prior to stranding for each stranded individual ( Figure 2). Uncertainty about foraging range was explored by fitting separate MSFR models (see Section 2.2) for each timeframe (2, 4, 6, 8 days).

| Prey availability
Relative fish abundances were estimated using data from the North Sea International Bottom Trawl Survey (NS-IBTS), available from the International Council for the Exploration of the Sea (ICES) (datras. ices.dk). NS-IBTS data were available for quarter 1 (January-March) and quarter 3 (July-September). where L is length class (in mm), indicated by the lower limit of that class, e is the resolution of the length, either 5 or 10 mm (depending on species), CPUE L is the catch per unit effort for length class L, and α and b are length-weight conversion parameters, the values of which were derived from Wilhelms (2013).
Generalized additive models (GAMs) were used to predict distribution for each prey species over the entire southern North Sea (south of 56°N latitude (Figure 3)). The response variable BPUE, logtransformed to reduce the effects of relatively high/low catches, was assumed to have a Gaussian error distribution. Covariates considered were longitude, latitude, depth, and year. Smoothing parameter selection was performed by restricted maximum likelihood (REML) (Wood, 2011). The model allowed the spatial pattern to change with time, by including a three-dimensional tensor product smooth for geographical space and year: For a given haul, the biomass per unit effort is represented by BPUE i,t having space coordinates i and a date/time t.
To avoid smoothing being adversely impacted by land boundaries we applied a soap film smoother (Wood et al., 2008). In generating the soap film, knots were placed over the data and land was set to zero which ensured smoothing toward data points and avoided predicting over the boundary. Comparing the soap film smoother with a conventional thin-plate regression spline smoother showed that the soap film improved the prediction of fish densities in areas with closely adjacent land boundaries (e.g., the Strait of Dover).
The predictions of the fitted model represent expected BPUE values. To estimate the true underlying fish biomass, predictions would need to be scaled using gear efficiency and catchability estimates. However, absolute estimates of fish biomass are not required to fit a MSFR (see Section 2.2).
Sandeels are not well represented in the NS-IBTS due to catchability issues. Therefore, for this species we used ICES estimates of sandeel spawning-stock biomass from other data sources, including commercial catches and dredge surveys (ICES, 2016). Gobies (Gobiidae) had to be excluded because they are almost absent in the NS-IBTS data due to their small sizes (Knijn et al., 1993), and there is no other source of data.
The relative availability, and associated uncertainty, of each main prey species to each porpoise prior to stranding was estimated as the relative amount of prey present within the area of sea within the estimated circular buffer (see Section 2.1.2). For each buffer, the SD of the availability of each prey species was obtained by parametric resampling the estimated coefficients from the fitted GAMs.

| Model development
A general equation allowing a single species functional response to take the form of a type I, II, or III is (Holling, 1959): where c is the predator consumption rate, is the attack rate, N is prey availability, t is the consumption/handling time, and m is a shape parameter.
The equation can be rewritten to include multiple prey species: where c i is the consumption of prey species i; i and t i are the attack rate and handling time of species i. There is a total of Z prey species in the system. We do not have information on harbour porpoise consumption rates, but the equation can be revised in terms of diet composition: Predicted foraging range prior to stranding assuming 2 days foraging using telemetry data. Green buffer denotes the estimated minimum enclosing circle (MEC) from Danish telemetry data. White circle is the stranding location of one porpoise. Purple buffer represents the predicted foraging range prior to stranding. The diameter of the green buffer was used as the radius of the purple buffer The spatial distribution of porpoise prey species for the southern North Sea illustrated for the year 2007, for Quarter one (January-March) and Quarter three (July-September). Density surfaces were produced from generalised additive models that were fitted to biomass per unit effort data (kg per half hour trawl, based on NS-IBTS catches  (Real, 1977). We compared two model types: a hyperbolic type II functional response with shape parameter m = 1 (model 1) and a sigmoidal type III functional response with m = 1.5 (model 2).
The relationship between relative prey availability and consumption was estimated for each main prey species in turn by setting the availability of all other prey to one of three specific constant levels (minimum, mean, and maximum

| Model prediction
To illustrate the model's ability to predict consumption under different regimes of prey availability, the estimated parameters of the best fitting model were used to predict diet composition in 2011 (a year of high sandeel spawning-stock biomass, SSB, in the southern North Sea) and 2020 (a year of low sandeel SSB, and the most current advice from ICES for all prey species considered), assuming similar prey distribution and porpoise stranding locations. Prey availability of 2011 was rescaled using estimates of SSB from the ICES Stock Assessment Database for 2020, following Smout et al. (2014). Changes in diet composition were estimated relative to the daily biomass or energy consumption of an average adult male porpoise (i.e., 1.7 kg or 6.7 MJ per day (Gallagher et al., 2018)). Estimates from the literature were used to convert biomass to energetic content (Table 5). Energy values for gobies were used for the "other" prey category because they were the most prevalent species in that group.

| Diet composition
Stomach content data were available from 455 harbour porpoises.
In models to estimate the diameter of minimum enclosing circles (MECs), all covariates had a VIF score lower than 1.4; therefore, multicollinearity could be disregarded. Model results are summarized in Table 1. Age, quarter, sex, and timeframe were all found to be significant predictors (p < .01) for the foraging range (MEC diameter) and explained 24.5% of the variation. Predicted foraging range was smaller for males than for females and for juveniles in comparison with adults. Foraging range was significantly smaller in spring in comparison with the other seasons.

| Prey availability
Correlograms of the final models of relative prey abundance indicated very weak autocorrelation, and deviance residuals were evenly spread. BPUE predictions in all grid cells, including unsurveyed cells, are shown in Figure 3. The final models explained between approximately one third to two thirds of the total observed variation in the BPUE values (Table 2).
As described above, the availability of each prey species was predicted for each individual porpoise, within the circular buffer that represented the foraging range for each timeframe (see Section 2.1.3). For illustration, Figure 4 displays the prediction of whiting availability for one porpoise for different timeframes.

| Multi-species functional response
The best MSFR model in terms of timeframe according to DIC scores ( Model predictions of diet composition captured the overall pattern in the observed diet (Table 4)

| D ISCUSS I ON
Integrating disparate datasets to model the MSFR for harbour porpoises in the southern North Sea provides a methodological framework that may be appropriate for other predators. Results from our case study show that sandeels are an important and possibly a preferred prey for harbour porpoise, thus increasing knowledge of the foraging ecology of this important marine predator.

| Method evaluation and sensitivity
Setting suitable spatial scales can be a major challenge in ecological studies and the accuracy of any modeled relationship between prey consumption and availability is strongly dependent on  achieving realistic spatio-temporal overlap. In this study, the foraging distributions of porpoises prior to stranding (our source of diet information) are unknown, so it is crucial to explore whether assumptions made about the foraging range of these animals are reasonable. Our novel approach was to find the most likely foraging area prior to stranding by predicting the range used as a function of time period based on telemetry data, and using the MSFR model fit to determine the appropriate time period of 4 days.
There is little relative difference in modeled prey distribution for each prey species in the areas where porpoises could have been foraging in the vicinity of the Dutch coast (Figures 3 and 4), indicating that the overall pattern of results is unlikely to vary much over the range of time periods modeled (2-8 days). This is confirmed by lack of variation in the emerging patterns of estimated attack rates or the shape of the functional response (Appendix S1).
In this case study, our methodology thus appears rather robust to this aspect of uncertainty.

| Ecological inference
Different shapes of the predator functional response have different implications for prey populations, especially at low prey densities. In our best fitting model with a sigmoidal type III functional response, predation mortality decreases when a prey species becomes rare and is indicative of prey switching when prey is at low abundance (that is, there is a change in preference dependent on prey abundance). This may result in persistence and/or stabilizing effects on predator-prey dynamics (Murdoch & Oaten, 1975) because it may prevent one prey species from outcompeting others (Roughgarden & Feldman, 1975). A type III response may result from a number of ecological mechanisms, including prey refuge (McNair, 1986), and learning time (Tinbergen, 1960).

F I G U R E 5
Relationship between prey availability and consumption by harbour porpoises for the 4 days MSFR model. Relationships are shown as a single-species plot at three different levels of alternative prey (all other prey) availability TA B L E 5 Predicted change in harbour porpoise diet composition in terms of percentages of total prey biomass (%M) and energy (%E) for 2011 and 2020, and prey energy density  Lawson et al. (1997).
Classically, the attack rate parameter a in the functional response equation can be interpreted as a form of relative preference of the predator for a certain prey type. Here, we interpret these values cautiously because of the nature of the prey abundance estimates we used. These were indices, scaled in proportion to maximum values, and they were not estimates of overall total biomass (which is difficult to calculate). Thus, for example the "maximum" value of sandeel abundance was 100 and so was the maximum value for whiting.
In this study, porpoises consumed a disproportionately larger proportion of the most abundant prey. Sandeel consumption remained high even when other prey were abundant and was considerably higher than the consumption of other prey at equal availability index values. At prey abundances similar to those available to our study animals, harbour porpoise diets often have a high proportion of sandeels (Jansen, 2013;Santos et al., 2004), and it also implies that sandeel availability might have a particularly strong effect on the consumption by porpoises of other prey species in this area. Habitat models of harbour porpoise in the North Sea have found that harbour porpoise density increases with decreasing distance to sandeel grounds (Gilles et al., 2016), suggesting that porpoises could be attracted to those areas.
Harbour porpoises in better body condition have been found to be more likely to have higher amounts of fatty fish, such as sandeels, in their diet (Leopold, 2015). Our results add to the body of evidence that sandeels are important to porpoises. Sandeels have high energy content and are abundant in the southern North Sea, forming an important forage fish resource that supplies a number of predator species including harbour porpoises, seabirds (Rindorf et al., 2000), and gray and harbour seals (Wilson & Hammond, 2019).
Despite considerable differences in predicted diet composition in 2011 and 2020 (Table 5), differences in predicted consumption were relatively small (~2%), illustrating little variation in overall biomass or energy intake. However, if energy-rich prey species (i.e., clupeids and sandeels) were reduced to low levels, this could result in porpoises needing to increase biomass consumed to avoid failing to meet their energetic requirements.
Sandeel abundance in the northwestern North Sea has declined since 2000 (MacDonald et al., 2019). Poor seabird breeding success in the northwest North Sea has been linked to a reduction in the availability and quality of sandeels (Wanless et al., 2005(Wanless et al., , 2018. Our results confirming the importance of sandeels to harbour porpoise and indicating their possible preference for this prey are consistent with the reduction of sandeel biomass in the northern North Sea being a driver of the distributional shift of porpoises from the northern to the southern North Sea between 1994 and 2005 (Hammond et al., 2013).
However, this needs further exploration of the impact of other potential drivers such as competition with other sandeel predators (i.e., sea birds, other marine mammals, foraging fish as well as fisheries).

| Data limitations
Foraging range was estimated from telemetry data collected in areas of the North Sea outside the study area. By including data only from the area believed to be most similar to the study area, we sought to minimize error in estimated foraging range. Estimates of foraging range using movement data are uncertain and conservative. Active swimming is faster than drifting, so true foraging range will be larger than that estimated from drifting alone. The fitted MSFRs gave similar results for different assumptions about the foraging area available to porpoises before stranding, so we conclude that the lack of telemetry data from the study area should not affect our conclusions appreciably. ARGOS data from telemetry tags are subject to location error which was not quantified in this study but is believed to be negligible in this context.
Our prey availability estimates assume that relative prey abundance is proportional to true abundance. We also assume that prey abundance reflects prey availability to predators; however, the relationship between prey abundance and prey availability is largely unknown, not least because differences in prey behavior (e.g., diurnal and seasonal variation in schooling and burying behavior) may affect this. Our methodology is appropriate if spatio-temporal trends in relative abundance reflect those in absolute abundance of prey, which seems a reasonable assumption.
Most fisheries surveys, including the North Sea IBTS, sample at a coarse spatial and temporal resolution. Some species, especially sandeels and gobies in this study, are poorly sampled. Given the importance of sandeels for many marine predators (Engelhard et al., 2013;Gilles et al., 2016;Wanless et al., 2005;Wilson & Hammond, 2019) and the lack of knowledge regarding spatio-temporal variability in their distribution and abundance, improving effective sampling and modeling of sandeel distribution would improve the quality of the inferences made from future studies. The inability to model sandeel and goby distributions spatially could have led to error in availability estimates, especially because sandeel distribution is extremely patchy (Wright et al., 2000) and largely unknown for gobies.
The importance of gobies could have been underestimated because they were excluded from the prey availability analysis. Although information on goby distribution and abundance is largely lacking, gobies are extremely abundant within Dutch coastal waters (Tulp et al., 2008). Therefore, it might be reasonable to assume that these species have a relatively consistent availability.
Care should be taken in making inferences from stranding data because they do not represent an unbiased sample of the population; there is likely an over-representation of individuals that are inexperienced, old, and/or in poor health  Ross & Wilson, 1996), would be an alternative way to look at harbour porpoise diet. However, the diet of by-caught porpoises diet could be biased toward target species of the fishery and "net selection" of inexperienced individuals (Santos & Pierce, 2003), and porpoises killed by seals or dolphins may be more vulnerable to predation and not representative of the population.
The predicted diet for 2011 and 2020 for an adult male porpoise assumed sampling the same porpoise stranding locations, fish distributions having the same pattern scaled by North Sea wide stock assessment estimates, and that the relationship between fish biomass and energy was constant for each species. However, these assumptions could be violated in several ways. For example, the distribution of porpoises and/or fish could have changed, there are differences in porpoise prey consumption according to sex, age (Booth, 2020;Leopold, 2015), and prey energy densities vary by size class, season, etc. (Pedersen & Hislop, 2001). Therefore, a more elaborate analysis is required to explore the impact of changing these factors on the predictions.

| Context and applications
Applying a Bayesian approach to model the MSFR appears to work well, allowing incorporation of uncertainty in prey availability and consumption estimates. These features, together with the resilience of the results, suggest that the modeled MSFR provides a strong methodological framework that can be applied (generalized) to a range of other species and might aid in quantifying the ecological role of other predators that consume a variety of prey. For example, similar data exist for seabirds (Wanless et al., 2005), gray seals, and harbour seals in the North Sea (Carter et al., 2020;Wilson & Hammond, 2019) and applying this framework could provide valuable new insights into their population dynamics, especially in the context of possible competition for prey between these two seal species (Wilson & Hammond, 2019). To take this further, the MSFR could be integrated into ecosystem models to predict and test how prey and predator populations are expected to change under different fisheries management and climatic scenarios that impact prey availability. This could also shed light on the extent of direct and indirect competition between marine mammals, seabirds, and fisheries and possibly on the outcomes of fisheries management and stock recovery programs.

ACK N OWLED G M ENTS
We would like to thank Iosu Paradinas, James Grecian, Esther Jones, and David Miller, who provided valuable support in the data analysis. Thanks to Jonas Teilmann and Rune Dietz, who were part of the harbour porpoise tagging team.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest. writing-review & editing (equal).

DATA AVA I L A B I L I T Y S TAT E M E N T
Data to estimate the relative availability of prey species are openly available from ICES webpage (datras.ices.dk). Data and code for modelling the MSFR are deposited in the Dryad Digital Repository: https://doi.org/10.5061/dryad.z8w9g hxdg.