Copyright © 2017 Everett, Baird, Buchanan, Bulman, Davies, Downie, Griffiths, Heneghan, Kloser, Laiolo, Lara-Lopez, Lozano-Montes, Matear, McEnnulty, Robson, Rochester, Skerratt, Smith, Strzelecki, Suthers, Swadling, van Ruth and Richardson. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Zooplankton Are Not Fish: Improving Zooplankton Realism in Size-Spectrum Models Mediates Energy Transfer in Food Webs

Ryan F. Heneghan1 \*, Jason D. Everett 2, 3, Julia L. Blanchard4 and Anthony J. Richardson1, 5

1 Centre for Applications in Natural Resource Mathematics, School of Mathematics and Physics, University of Queensland, Brisbane, QLD, Australia, 2 School of Biological, Earth and Environmental Sciences, University of New South Wales, Sydney, NSW, Australia, 3 Sydney Institute of Marine Science, Sydney, NSW, Australia, 4 Institute of Marine and Antarctic Studies, University of Tasmania, Hobart, TAS, Australia, 5 Ecosciences Precinct, CSIRO Ocean and Atmosphere, Brisbane, QLD, Australia

#### Edited by:

Christian Lindemann, University of Bergen, Norway

#### Reviewed by:

Jan Marcin Weslawski, Institute of Oceanology (PAN), Poland Francis Poulin, University of Waterloo, Canada

\*Correspondence:

Ryan F. Heneghan<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Marine Ecosystem Ecology, a section of the journal Frontiers in Marine Science

Received: 03 August 2016 Accepted: 29 September 2016 Published: 19 October 2016

#### Citation:

Heneghan RF, Everett JD, Blanchard JL and Richardson AJ (2016) Zooplankton Are Not Fish: Improving Zooplankton Realism in Size-Spectrum Models Mediates Energy Transfer in Food Webs. Front. Mar. Sci. 3:201. doi: 10.3389/fmars.2016.00201 The evidence for an equal distribution of biomass from bacteria to whales has led to development of size-spectrum models that represent the dynamics of the marine ecosystem using size rather than species identity. Recent advances have improved the realism of the fish component of the size-spectrum, but these often assume that small fish feed on an aggregated plankton size-spectrum, without any explicit representation of zooplankton dynamics. In these models, small zooplankton are grouped with phytoplankton as a resource for larval fish, and large zooplankton are parameterized as small fish. Here, we investigate the impact of resolving zooplankton and their feeding traits in a dynamic size-spectrum model. First, we compare a base model, where zooplankton are parameterized as smaller fish, to a model that includes zooplankton-specific feeding parameters. Second, we evaluate how the parameterization of zooplankton feeding characteristics, specifically the predator–prey mass ratio (PPMR), assimilation efficiency and feeding kernel width, affects the productivity and stability of the fish community. Finally, we compare how feeding characteristics of different zooplankton functional groups mediate increases in primary production and fishing pressure. Incorporating zooplankton-specific feeding parameters increased productivity of the fish community, but also changed the dynamics of the entire system from a stable to an oscillating steady-state. The inclusion of zooplankton feeding characteristics mediated a trade-off between the productivity and resilience of the fish community, and its stability. Fish communities with increased productivity and lower stability were supported by zooplankton with a larger PPMR and a narrower feeding kernel–specialized herbivores. In contrast, fish communities that were stable had lower productivity, and were supported by zooplankton with a lower PPMR and a wider feeding kernel—generalist carnivores. Herbivorous zooplankton communities were more efficient at mediating increases in primary production, and supported fish communities more resilient to fishing. Our results illustrate that zooplankton are not just a static food source for larger organisms, nor can they be resolved as very small fish. The unique feeding characteristics of zooplankton have enormous implications for the dynamics of marine ecosystems, and their representation is of critical importance in size-spectrum models, and end-to-end ecosystem models more broadly.

Keywords: zooplankton dynamics, marine size-spectrum, end-to-end modeling, fish productivity, ecosystem stability

### INTRODUCTION

In the 50 years since Sheldon et al. (1967) first hypothesized an equal concentration of biomass from bacteria to whales, a range of size-spectrum models have been developed to explain this remarkable consistency (Andersen et al., 2015; Guiet et al., 2016b). Size-spectrum models represent the entire marine community as a size distribution, and traditionally do not resolve species identity. Their simplicity and parsimonious parameterization makes it possible for them to be used to investigate human impacts at the community level, including fishing (e.g., Andersen and Pedersen, 2010; Jacobsen et al., 2014; Law et al., 2016), climate change (e.g., Blanchard et al., 2012; Woodworth-Jefcoats et al., 2013; Barange et al., 2014; Dueri et al., 2014), and habitat loss (Rogers et al., 2014).

The focus of these models has been on higher trophic levels—primarily fish and fishing—and in recent years there has been considerable effort in improving their parameterization (Andersen et al., 2015; Guiet et al., 2016b). For example, recent theoretical developments now allow size-spectrum models to resolve different functional groups and even species by their traits, and this has been implemented for various fish (e.g., Blanchard et al., 2014; Dueri et al., 2014; Zhang et al., 2015). The focus on fish has meant that the dynamics of the plankton-dominated lower trophic levels has been neglected in model formulations. Zooplankton, as the main consumers of phytoplankton and prey of small fish are the chief intermediaries between primary production and higher trophic levels, and thus play a critical role in marine food web dynamics (Carlotti and Poggiale, 2010; Mitra and Davis, 2010).

In current dynamic size-spectrum models, the minimum size of the dynamic consumer spectrum extends to mesozooplankton. For smaller zooplankton, there are three common representations. First, phytoplankton and small zooplankton are represented as a fixed resource spectrum (with a varying intercept and a slope held at −1), and are considered only as a food source for the smallest fish size classes (Maury et al., 2007; Blanchard et al., 2009, 2011, 2012; Law et al., 2009; Datta et al., 2010; Guiet et al., 2016a). Second, the phytoplankton and small zooplankton spectrum is determined by an external nutrient–phytoplankton–zooplankton (NPZ) model, with no predation feedback from the larger dynamic size classes (Woodworth-Jefcoats et al., 2013; Lefort et al., 2015; Le Mézo et al., 2016). Third, phytoplankton and small zooplankton are modeled as a semi-chemostat system, with a fixed carrying capacity and predation feedback from higher trophic levels (Hartvig et al., 2011; Blanchard et al., 2014; Scott et al., 2014; Zhang et al., 2015, 2016). The latter approach is the only one in which the size-spectrum of fish dynamically interacts with phytoplankton and small zooplankton. These current representations essentially group smaller zooplankton and phytoplankton together as food for the smallest dynamic size classes, and resolve larger zooplankton as small fish.

Assuming zooplankton have the same dynamics as phytoplankton or small fish is not only incorrect, but could have considerable effects on energy transfer in food webs. Zooplankton have feeding characteristics distinctly different from fish. For instance, the average predator–prey mass ratio (PPMR; in grams of wet weight) for fish is typically around 100 (Jennings et al., 2001) but for zooplankton it is >1000 (Kiørboe, 2008; Wirtz, 2012). Additionally, zooplankton exhibit vast phylogenetic biodiversity, with at least eight phyla commonly present (crustaceans, chordates, chaetognaths, molluscs, cnidarians, echinoderms, ctenophores, and annelids), each with considerable differences in their feeding ecology, from passive suspension grazing of the water column to active ambushing and carnivory (Kiørboe, 2011). Further complicating their feeding, various species of jellyfish, copepods, and microzooplankton can switch between suspension and ambush feeding modes, and this is reflected in different optimal prey sizes realized by the same species (Landry, 1981; Goldman and Dennett, 1990; Saiz and Kiørboe, 1995).

Size-based predation is the key driver of dynamics in sizebased ecosystems (Jennings et al., 2001; Woodward et al., 2005; Andersen et al., 2016) and is broadly defined by five key parameters: (1) PPMR; (2) search rate coefficient; (3) body-size exponent, which determines how the search rate scales with body-size; (4) assimilation efficiency; and (5) the width of the feeding kernel (the diet breadth around the preferred PPMR), and modeling studies of the size-spectra of fish have shown that these parameters have a large effect on food web dynamics (Law et al., 2009; Datta et al., 2011; Zhang et al., 2013). For instance, a wider feeding kernel and lower PPMR dampens traveling waves through the fish community size-spectrum (Blanchard, 2008; Law et al., 2009; Zhang et al., 2013). Further, there is evidence that higher assimilation efficiency has a similar effect on the steady state of the size-spectrum (Datta et al., 2011). The sensitivity of ecosystem dynamics to parameterization of fish feeding characteristics strongly suggests that zooplankton feeding characteristics could be important to energy transfer through the food web. Therefore, the first step toward including zooplankton in end-to-end size-spectrum models is to include an accurate representation of their feeding characteristics.

The extensive experimental work elucidating zooplankton feeding characteristics has formed the basis of several recent syntheses of size-based feeding (Fuchs and Franks, 2010; Kiørboe, 2011; Wirtz, 2012, 2014) and provides an opportunity for improving zooplankton parameterization in size-spectrum models. Wirtz (2012) used the data collected by Hansen et al. (1994) and Fuchs and Franks (2010) to develop a mechanistic model that links zooplankton PPMR with their feeding characteristics. In another paper, Wirtz (2014) derived an ideal feeding kernel width for zooplankton from simple biomechanical laws, which agrees well with empirical data. Fuchs and Franks (2010) synthesized data from previous studies to investigate the relationship between zooplankton PPMR and the width of the feeding kernel. They found that the feeding kernel width decreased with decreasing PPMR, suggesting increasing selectivity amongst individuals who prey on individuals closer to their own size. Kiørboe (2011) found that the size-specific zooplankton search rate is independent of body size across seven different functional groups.

Here, we evaluate how the size-dependent feeding characteristics of zooplankton affect the dynamics of higher trophic levels in size-structured pelagic ecosystems. We extract feeding characteristics from a range of syntheses of size-based feeding (Fuchs and Franks, 2010; Kiørboe, 2011, 2013; Wirtz, 2012, 2014) and implement them in a dynamic size-spectrum model framework (Datta et al., 2010; Andersen et al., 2015; Guiet et al., 2016b). To our knowledge this is the first dynamic size-spectra model to resolve predation-based growth and mortality of zooplankton. The model has three components—a static phytoplankton resource spectrum and two dynamic spectra representing a general zooplankton and fish community, respectively. In our model, biomass flows from smaller to larger size classes as a consequence of larger organisms consuming smaller organisms, and growth at one size is balanced by mortality in smaller size classes. We separate our findings in three parts. In Zooplankton Are Not Fish, we provide a size-spectrum model using the best parameter estimates from the literature, and establish the individual effect each of the five key zooplankton feeding parameters has on the community size-spectrum, by comparing against a base model where zooplankton are parameterized as just another fish community. In Sensitivity Analysis, we assess how varying the feeding characteristics of the zooplankton community impact the stability and productivity of the fish community size-spectrum. Finally in Mediating Primary Production and Fishing, we evaluate how the feeding characteristics of different zooplankton functional groups—salps, chaetognaths, herbivorous copepods, flagellates, and carnivorous copepods—mediate changes in phytoplankton abundance and increased fishing mortality on the fish community size spectrum. The purpose of this study is not to give a quantitative evaluation of zooplankton or fish abundance, rather we wish to illustrate how incorporating zooplankton-specific feeding characteristics could affect the dynamics of sizestructured ecosystems. Our ultimate aim is to investigate how zooplankton feeding characteristics influence energy transfer from phytoplankton and fish, and thus move toward a more realistic and consistent parameterization for the zooplankton component of size-spectrum models.

### METHODS

### The Model

We developed a size-spectrum model that consists of a sizespectrum comprised of three communities: phytoplankton, zooplankton, and fish (**Figure 1**; **Tables 1**,**2**). The phytoplankton component covers the smallest size classes [wp,wz] and is held constant as a background resource spectrum for zooplankton. Size-dependent processes of growth and mortality drive the zooplankton and fish components. These two components are delineated by different size ranges and feeding characteristics. The zooplankton community covers the size range between phytoplankton and fish [wz, Wz], and the fish community covers the largest size classes [wf , Wf ], although some of the smallest fish size classes extend into the zooplankton range (from wf = 0.1 g to Wz = 1 g). Fish community size classes that extend into the zooplankton range represent larvae and very small fish that are smaller than the largest zooplankton. Predation is sizedependent, with big things eating smaller ones, so depending on their size, zooplankton can feed on phytoplankton, smaller zooplankton, and the smallest fish size classes. Similarly, fish feed on zooplankton and smaller fish. The temporal dynamics of the zooplankton and fish communities are governed by separate McKendrick–von Foerster equations with second-order diffusion terms (Datta et al., 2010),

$$\begin{split} \frac{\delta}{\delta t} N\_i(\boldsymbol{\omega}, t) &= -\frac{\delta}{\delta \boldsymbol{w}} \left( \mathcal{g}\_i(\boldsymbol{\omega}, t) N\_i(\boldsymbol{\omega}, t) \right) - \mu\_i(\boldsymbol{\omega}, t) N\_i(\boldsymbol{\omega}, t) \\ &\quad + \frac{1}{2} \frac{\partial^2}{\partial \boldsymbol{w}^2} \left( f\_i(\boldsymbol{\omega}, t) N\_i(\boldsymbol{\omega}, t) \right). \end{split} \tag{1}$$

The density of individuals in community i (where i is either zooplankton or fish) per unit mass per unit volume (g−1m−3 ) is denoted by Ni(w, t). Growth, mortality, and diffusion rates of individuals of group i at size w and time t, are denoted by gi (w, t), µi (w, t), and fi (w, t), respectively. In this context, the diffusion

governed by Equation (1) and the Equations in Table 2.

Sources: 1. Zhou et al. (2010), 2. Blanchard et al. (2009), 3. Benoît and Rochet (2004), 4. Blanchard et al. (2011), 5. Andersen et al. (2015), 6. Wirtz (2012), 7. Wirtz (2014), 8. Kiørboe (2011), 9. Peters (1983), 10. Hall et al. (2006), 11. Brown et al. (2004), 12. Barnes et al. (2011), 13. Rousseaux and Gregg (2015).

term allows the model to incorporate demographic variation in the growth rates of each community. That is, within each community two individuals of the same weight eating the same food will not grow by the same amount (Datta et al., 2010). This not only increases model realism, but the stability of the system steady state, over the traditional first-order McKendrick–von Foerster equation (Datta et al., 2011).

Phytoplankton dynamics in the background resource spectrum are not explicitly modeled, with the density of individuals held constant through time:

$$N\_{\mathcal{P}}\left(\omega,t\right) = a\omega^{-b},\tag{E11}$$

Equations are also found in **Table 2**. We use an exponent of −1 for the background spectrum, implying equal biomass over logarithmically equal body–mass intervals in keeping with past dynamic size-spectrum models Benoît and Rochet, 2004; Law et al., 2009, 2016; Blanchard et al., 2012). The coefficient a for the background spectrum (Equation 11) was calculated using the empirical equation from Barnes et al. (2011) and annual median chlorophyll-a concentrations for different ocean basins between 2005 and 2010 (Rousseaux and Gregg, 2015). Unless specified otherwise, we use the global median chlorophyll-a value (0.16 mg m−3 ) to give a value of 0.017 for a.

From the predator's perspective, total consumption depends on the total biomass of suitable prey. For an individual of size w at time t from community i, this is determined by the predator's search rate:

$$V\_i(\omega) = \gamma\_i \mathcal{W}^{\alpha\_i},\tag{E5}$$

and the density of suitable prey:

$$D\_i\left(\boldsymbol{w},t\right) = \int\_{\boldsymbol{w}\_P}^{\boldsymbol{w}} \phi\_i\left(\boldsymbol{w},\boldsymbol{w}'\right) \sum\_j N\_j\left(\boldsymbol{w}',t\right) \boldsymbol{w}' d\boldsymbol{w}'.\tag{E6}$$

The growth rate of an individual of size w at time t is fuelled by consumption of prey from smaller size classes:

$$\mathcal{g}\_i(\boldsymbol{w}\_i, t) = K\_i V\_i \left( \boldsymbol{\omega} \right) D\_i \left( \boldsymbol{\omega}, t \right), \tag{E7}$$

where Ki is the assimilation efficiency of community i.

#### TABLE 2 | Model equations with their units.

An equation number is given that is used in the main text.

Kiørboe (2011) found that the search rate (Equation 5) for zooplankton, across a wide range of taxa is largely independent of organism size (αZ = 1.01). This stands in contrast to the scaling for fish (αF = 0.8; Peters, 1983) that implies the specific search-rate per unit mass declines with increasing body size. Further, the search rate coefficient is higher for zooplankton γZ = 875g−αZm−3 year−1 (Kiørboe, 2011), compared to fish γF = 640g−αFm−3 year−1 (Peters, 1983). The probability that a predator of size w will consume an individual of size w ′ is given by the log-normal function (Equation 4), where βi and σi are community i 's PPMR and feeding kernel width.

In previous size-spectrum models the PPMR is held constant across the entire size range of the community (Andersen et al., 2015). For zooplankton, the wide variation in observed PPMR across phyla suggests a constant value across all zooplankton size classes is inappropriate (Wirtz, 2012). We have thus used the mechanistic formulation from Wirtz (2012) who argues that for zooplankton, PPMR will increase non-linearly as predator size increases, due to the nonisometric scaling of feeding-related apparatus with body size (**Figure 2**):

$$\beta\_{\rm Z}(w) = \left( \exp \left( 0.02 \ln \left( D\_{\rm w} / D\_{0} \right)^{2} - m + 1.832 \right) \right)^{3}, \qquad \text{(E2)}$$

where Dw is the predator equivalent spherical diameter (ESD) in µm:

$$D\_{\le} = 2\sqrt[3]{3\pi \times 10^{12}/4\pi}. \tag{E1}$$

The mechanistic model from Wirtz (2012) also allows the range of feeding modes across different zooplankton functional groups to be quantitatively incorporated. The feeding mode of the zooplankton community is denoted by m, and ranges from −3 to 2. A larger, positive m-value (say m = 2) suggests a more raptorial, carnivorous feeding strategy with a lower PPMR (**Figure 1**). Alternatively, a negative m-value (say m = −3) implies a more passive, herbivorous feeding strategy. For the fish community, we set βF = 100 (Andersen et al., 2015).

A wider feeding kernel means an individual feeds over a wider range of size classes, implying a more generalist feeding strategy. Wirtz (2014) obtained a general feeding kernel width for zooplankton of 0.75 log10 grams body-size from simple biomechanical laws, and found this value agreed well with measured values from different species. Fuchs and Franks derived an empirical equation that links zooplankton PPMR βZ with the feeding kernel width:

$$
\sigma\_{\overline{Z}} = \mathfrak{Z} \left( 0.05 \log\_{10} \left( \beta\_{\overline{Z}} \right) + 0.11 \right). \tag{E3}
$$

This equation suggests a positive relationship between the width of the feeding kernel and the PPMR. In other words, more active, carnivorous groups (m > 0), have a narrower, more selective prey size range compared to passive, filter feeding groups (m < 0) that have a wider, more generalist prey size range. For the fish community, we set σF = 1 (Andersen et al., 2015).

For individuals in community i, a fraction Ki of consumed biomass, the assimilation efficiency, is assimilated into new biomass. Observational and experimental work across different zooplankton functional groups show that assimilation efficiency of ingested food into new biomass ranges from 0.3 to 0.9 (Landry et al., 1984; Kiørboe, 2008; Montagnes and Fenton, 2012; Abe et al., 2013). Assimilation efficiency of copepods (Landry et al., 1984), dinoflagellates and larval fish (Kiørboe, 2008) depends on whether they were acclimated to low or high food environments; those from low food environments have a higher assimilation efficiency compared to those from high food environments. Similarly, the density and type or prey available had a significant effect on zooplankton assimilation efficiency—higher density and higher carbon content of prey gave lower assimilation efficiencies (Montagnes and Fenton, 2012; Abe et al., 2013). In previous models, assimilation efficiency for zooplankton is usually held constant at 0.70 (e.g., Zhou, 2006; Fuchs and Franks, 2010; Ward et al., 2012, 2014). Unless specified otherwise, we keep KZ = 0.7 to follow previous size-based planktonfocused models (Baird and Suthers, 2006; Zhou, 2006; Stock et al., 2008; Fuchs and Franks, 2010; Banas, 2011). For the fish community, we set KF = 0.6, as this is a common values given in previous dynamic size-spectrum models (Andersen et al., 2015).

From the prey's perspective, the total predation pressure from all larger size classes gives the predation mortality rate:

$$\mu\_{\mathcal{P}}\left(\boldsymbol{w},t\right) = \sum\_{j} \mathbb{I}\_{\{\boldsymbol{w} < \overline{W\_{j}}\}} \int\_{\boldsymbol{w}}^{\overline{W\_{j}}} \phi\_{\boldsymbol{j}}\left(\boldsymbol{w}',\boldsymbol{w}\right) V\_{\boldsymbol{j}}(\boldsymbol{w}) \mathbf{N}\_{\boldsymbol{j}}\left(\boldsymbol{w}',t\right) d\boldsymbol{w}'.\tag{E8}$$

To account for other sources of mortality (e.g., disease), we include a U-shaped intrinsic mortality term:

$$
\mu\_{0\_l}(\boldsymbol{\omega}, t) = B\_0 \boldsymbol{\omega}^\varepsilon + \mathbf{S}\_{0\_l} \boldsymbol{\omega}^s \tag{E9}
$$

that covers non-predation sources of mortality such as disease and senescence (Brown et al., 2004; Hall et al., 2006; Blanchard et al., 2009, 2011). Since individuals grow through time, the background mortality term describes rapidly decreasing background mortality in the early stages of life, a constant mortality for middle-age individuals, and an increasing mortality with senescence. The increase in senescence mortality with body size acts as a closure term for the largest size classes, by preventing a buildup of very large individuals (Andersen et al., 2015).

#### Community Characteristics

To evaluate effects of feeding characteristics of the zooplankton community on the fish community, we calculated several community-level measures. The total biomass of community i was obtained by integrating the abundance in all size classes:

$$B\_i\left(t\right) = \int\_{\overline{W}\_i}^{\overline{W\_i}} \nu \, N\_i\left(\nu, t\right) \, d\nu. \tag{E15}$$

Similar to Blanchard et al. (2011), we defined the total throughput of community i as the total consumption rate:

$$T\_i\left(t\right) = \int\_{\mathcal{W}\_i}^{\overline{W\_i}} \nu \, V\_i\left(\nu\right) D\_i\left(\nu, t\right) N\_i\left(\nu, t\right) d\nu. \tag{E16}$$

The production to biomass ratio of a community i —where production was defined as the total flux out of the community from all sources of mortality (Brown et al., 2004)—was used to evaluate the total energy flux through a community:

$$PB\_i\left(t\right) = \left(\int\_{\mathcal{W}\_i}^{\overline{W\_i}} \nu \,\mu\_i\left(\nu, t\right) N\_i(\nu, t) d\nu\right) / \left(\int\_{\mathcal{W}\_i}^{\overline{W\_i}} \nu \, N\_i(\nu, t) d\nu\right). \tag{E17}$$

Total throughput is a measure of how energy moves internally through the system from predation processes, whereas the production to biomass ratio is an indicator of how much new biomass is produced to replace biomass lost to mortality, per unit of existing biomass. To evaluate the transfer efficiency from zooplankton to the fish community, we calculate the ratio of total fish biomass to total zooplankton biomass:

$$FZ\left(t\right) = B\_F\left(t\right) / B\_Z\left(t\right). \tag{\text{E.18}}$$

This is similar to the approach taken in previous studies evaluating the transfer efficiency of phytoplankton to zooplankton (Friedland et al., 2012; Havens and Beaver, 2013).

Heneghan et al. Zooplankton Are Not Fish

We use two measures to evaluate the stability and total variability of the system. First, the nature of the system steady state was determined using the Newton–Raphson multidimensional root-finding method (Press et al., 2007). For each configuration of zooplankton feeding characteristics in this study, the abundance of the zooplankton and fish communities was taken after 20 years. The stability of this abundance was determined by the maximum real part (λmax) of the eigenvalues of the Jacobian matrix calculated with the Newton–Raphson method. If λmax < 0, the entire system will settle into a stable equilibrium over time where the abundance in each size-class does not fluctuate. The more negative λmax is, the faster the system will recover from local perturbations to the steady state. Alternatively, if λmax > 0, over time the system will settle into a repeating, periodic traveling wave of abundance from smaller to larger size-classes. In this case, the more positive λmax is, the faster the traveling wave moves through the size classes. Second, we measured the total variability of the system by calculating the coefficient of variation (CV) of the time-series of biomass over the last 10 years of the simulation. The Newton–Raphson stability analysis and the CV work together; the first will identify if the steady state is stable or oscillating, and the CV gives a measure of the magnitude of these oscillations through time.

#### Numerical Implementation

Dynamics of the zooplankton and fish communities are modeled with Equation (1), which we solve numerically using a secondorder semi-implicit upwind finite difference scheme (Press et al., 2007). We present the results in log10 space for ease of interpretation, mathematical convenience and comparison with previous work. For the numerical implementation we discretize the dynamic size range [10−5 , 106 ] into equal 0.1 log10 size intervals (on a log10 gram scale), and use a daily-time step for the time interval. We chose these values to discretize the time and weight ranges to ensure convergence in our numerical implementation without requiring unnecessary computational effort, in keeping with past studies (Press et al., 2007; Plank and Law, 2012; Zhang et al., 2013; Law et al., 2016). For simplicity we are not explicitly modeling reproduction, thus the abundances of the smallest size classes in the zooplankton and fish communities are held constant. This implies that we are assuming constant recruitment for zooplankton and fish (Law et al., 2009; Blanchard et al., 2012). The assumption of constant recruitment permits a clearer evaluation of how the feeding characteristics of the zooplankton affect the dynamics of a fish community, in keeping with previous community size-spectrum models (e.g., Benoît and Rochet, 2004; Maury et al., 2007; Law et al., 2009; Zhang et al., 2013). For the zooplankton community, the density of individuals in the smallest size class is determined from the continuation of the phytoplankton size-spectrum:

$$\text{Nz } (\boldsymbol{w}\_{\boldsymbol{z}}) = a \boldsymbol{w}\_{\boldsymbol{z}}^{-b}, \tag{E12}$$

and the density of the smallest size class in the fish community is held equal to the equivalent zooplankton size class:

$$N\_{\rm F} \left( \boldsymbol{w}\_{\rm f}, t \right) = N\_{\rm Z} (\boldsymbol{w}\_{\rm f}, t). \tag{E13}$$

We ran each simulation for a 20-year period. In each simulation, our initial condition starts the zooplankton and fish community spectra as a continuation of the resource spectrum (Equation 11). If the solution was stable (λmax < 0), there would initially be some oscillations around the steady state that would diminish over time. For a stable solution, the closer λmax was to zero the greater the initial variance and the longer it took for the system to stabilize. When the solution was a traveling wave (λmax > 0), the variance of the system and magnitude of the oscillations would increase over time until the steady state was achieved. In this situation, the closer λmax was to zero, the longer the system took to find the unstable steady state. In all simulations the system achieved steady state within the first 5 years, therefore we discarded the first 10 years as a burn-in period.

### Zooplankton Are Not Fish

To establish the individual effect each of the five zooplankton feeding parameters has on the fish community, we begin with a base model where zooplankton are parameterized as another general fish community. From the base model, we build up to a model where the zooplankton community feeding characteristics are parameterized to represent a general, mixed zooplankton community. To do this, we use m = 0 to represent the average PPMR of a zooplankton community characterized equally by herbivorous and carnivorous feeding behavior, and set σZ = 0.75, KZ = 0.7, γZ = 875g−αZm−3 year−1 , and αZ = 1.01 to reflect the average feeding characteristics of zooplankton across multiple functional groups.

We change each zooplankton feeding parameter one at a time, then all together, and evaluate their individual relative impact on fish community measures against the base model, by calculating the change in the measure against the base model. For example, the relative fish biomass (rFB) for a new parameterization of the zooplankton community is obtained by dividing the fish biomass from the new model by the fish biomass from the base model.

### Sensitivity Analysis

In this section, we assess how variation in the feeding characteristics of the zooplankton community affects the productivity and stability of the fish community. We focus on zooplankton feeding mode (m), feeding kernel width (σZ), and assimilation efficiency (Kz), since these parameters vary across different zooplankton functional groups and environmental conditions. We vary m between −3 and 2, σZ between 0.4 and 2.2, and KZ between 0.3 and 0.9.

### Mediating Primary Production and Fishing

In our final section, we assess how the feeding characteristics of different zooplankton functional groups affect the productivity and stability of the fish community, and mediate increased primary production and fishing pressure, by evaluating the effect of these changes on the average total biomass of the fish community. We use the m-values from Wirtz (2012) for five different zooplankton functional groups (salps, chaetognaths, herbivorous copepods, flagellates, and carnivorous copepods) and a general zooplankton community (**Table 1**). The width of the feeding kernel for each of the six groups was determined with Fuchs and Franks' (2010) empirical Equation (Equation 3), which links the average zooplankton community PPMR with the feeding kernel width. For all groups, we hold the search rate and assimilation efficiency constant (see **Table 1**).

We used chlorophyll-a concentrations from two ocean basins—the North Central Pacific (0.06 mg m−3 ) and the North Atlantic (high concentration, 0.28 mg m−3 )—to give a range of coefficient values (intercept of the spectrum; a) between 0.010 and 0.024, which corresponds to a total phytoplankton abundance in the background resource spectrum of between 0.23 and 0.55 g−1m−3 . To include fishing pressure, we incorporate an additive fishing mortality term with a value between 0 and 2 year−1 , for all individuals in the fish community > 200 g.

### RESULTS

### Zooplankton Are Not Fish

The base model (denoted as the dashed line in each of the sub-plots in **Figure 3**) was a stable spectrum (λmax of −0.58), with the dynamic zooplankton and fish communities essentially a continuation of the static background spectrum in the base model.

Individually changing the zooplankton assimilation efficiency KZ from 0.6 and 0.7 (**Figure 3A**) increased the total throughput and production to biomass ratio of the fish community, in comparison to the base model (**Table 3**), and increased resilience of the entire system to local perturbations, with λmax = −0.71. Increasing the zooplankton community search rate coefficient (γZ) from 640 to 875 (g−αZm−3 yr−1 ) (**Figure 3B**), had a negligible effect on the total biomass or productivity of the fish community, compared to the based model (**Table 3**), however it did increase the stability of the system, with λmax = −0.76. Changing the search rate exponent for the zooplankton community (αZ; **Figure 3C**) from 0.82 to 1.01 reduced the total fish biomass by almost 70%, and reduced the relative production to biomass ratio (45% decrease) and relative total throughput (87% decrease), against the base model. Updating γZ decreased the resilience of the system, with λmax = −0.04, however the steady state remained a stable spectrum.

Individually reducing the zooplankton feeding kernel (σZ; **Figure 3D**) from 1 to 0.75, and changing the PPMR (**Figure 3E**) of the zooplankton component increased the total biomass, throughput and production to biomass ratio of the fish community, in comparison to the base model (**Table 3**). Changing the PPMR of the zooplankton gave the most significant increase in relative production to biomass (75% increase) and relative total throughput (335% increase) of the fish community. Only changing the zooplankton PPMR affected the relative fish to zooplankton biomass significantly, with a 22% increase against the base model. Changing the feeding kernel width and the PPMR for the zooplankton community changed the steady state from a stable spectrum to an oscillating system. Between the two parameters, changing PPMR gave the fastest oscillations, with

FIGURE 3 | The zooplankton and fish community size-spectra when various parameters are updated (A–E) one at a time and (F) all together. The dashed lines in each plot represent the zooplankton and fish communities in the base model parameterization, and the solid lines denotes the average abundance of the fish and zooplankton communities over 10 years in the modified model. The shaded areas show the regions of the traveling wave solutions over 10 years if the steady state is unstable.

TABLE 3 | Fish community biomass (FB), fish to zooplankton biomass ratio (F:Z), fish community production to biomass ratio (P:B), and throughput (TP) relative to the base model (r), the variation in fish community biomass (coefficient of variation; CV) and the maxmium real part of the Jacobian (λmax) when the zooplankton community feeding parameters are updated one at a time, and all-together.

The system has a stable steady state when λmax < 0, and an unstable, oscillating steady state when λmax > 0.

λmax = 0.65 compared to λmax = 0.24 for σZ. Further, the magnitude of the oscillations was larger when the zooplankton community PPMR was updated (CV = 0.28), compared to σZ (CV = 0.07, **Table 3**, **Figures 3D,E**).

When all parameters were changed for the zooplankton community (**Figure 3F**) there were significant increases against the base model in total fish biomass (69%), the fish to zooplankton biomass ratio (44%), and the fish community production to biomass ratio and total throughput (44 and 140%, respectively). Except for the relative fish to zooplankton biomass ratio, the increase in the total fish biomass and productivity measures were lower when all the parameters were updated, compared to just updating the zooplankton PPMR. (**Table 3**) The overall system was not stable (λmax = 0.47), and the magnitude of the oscillations through the system were higher than any seen in a system with a single parameter updated, with CV = 0.62. However, oscillations were slower compared to the system where only zooplankton PPMR was updated.

### Sensitivity Analysis

The total biomass of the fish community increases exponentially as m decreases (**Figure 4A**). From m = 2 to −3 (corresponding to an average zooplankton community PPMR range of 1—7.5), total fish biomass increases over 3 orders of magnitude (0.3 g m−3 for m = 2 to 620 g m−3 for m = −3). The exponential increase in fish biomass and productivity measures with respect to zooplankton PPMR starts at around m = −0.5, which corresponds to an average zooplankton community PPMR of 4.5. Similarly, a smaller feeding kernel width (σZ)—indicating a predator that feeds on a narrower size range of prey—results in an almost exponential increase in total fish community biomass (**Figure 4B**). From σZ = 0.4 to 1.1, total fish biomass increases from 0.25 g m−3 to 2.9 g m−3 . There is a roughly linear, positive relationship between total fish biomass and the assimilation efficiency KZ of the zooplankton community (**Figure 4C**). As KZ increases from 0.3 to 0.9, total fish biomass increases from 0.20 to 0.65 g m−3 . Similar patterns can be seen in the relationship between zooplankton PPMR, feeding kernel width and assimilation efficiency and the fish community productivity measures (Supplementary Figure 1).

The fish to zooplankton biomass ratio peaks at around m = 0 (3.60), and stays around 2.9 for m < −1, and decreases for m > 0.5 to settle around 2.1 (**Figure 4D**). For σZ, the fish to zooplankton biomass ratio peaks at σZ = 0.5 around 4.8, before uniformly declining as σZ increases (**Figure 4E**). There is minimal change in fish to zooplankton biomass ratio with increasing KZ, which suggests zooplankton biomass and fish biomass increase at the same rate (**Figure 4F**). Except for m > 0.5, the fish to zooplankton biomass ratio was higher than the base model (dashed line in **Figures 4 D–F**) across the ranges of m, σZ and KZ.

With σZ = 0.75, the CV is zero for m-values above 1 (**Figure 5A**) which corresponds to a stable steady state region (**Figure 5D**). The CV increases as m decreases from 0.5 to −0.5, which implies increasing variability in total fish biomass as the zooplankton community shifts from carnivorous to herbivorous feeding behavior. The CV stabilizes between 1 and 1.5 for m < 0. This suggests that even though the total fish community biomass is still exponentially increasing as m becomes more negative, the relative variation in fish biomass through time does not increase. There is a negative relationship between increasing σZ and CV, indicating increasing stability with a wider feeding kernel (**Figure 5B**). A similar pattern is observed in **Figure 5D**; the range of m-values that enable a stable system is larger, as σZ increases. The CV of the fish community varies across the range of KZ-values but within a much smaller range than the other two parameters (**Figure 5C**). Increasing KZ slightly increases the minimum σZ, and decreases the minimum m required for a stable steady state (**Figures 5E,F**).

### Mediating Primary Production and Fishing

Our results suggest a trade-off between the stability of the overall system and the total average fish productivity and biomass for most zooplankton groups (**Figure 6** and **Table 4**). The herbivorous salp community (m = −2.68, PPMR ≈ 7) is an exception (**Figure 6A**). It supports the most abundant and productive fish community, yet it is a more stable system overall than the one dominated by herbivorous copepods and chaetognaths (**Figures 6B,C**; **Table 4**). The salp community has the widest feeding kernel (σZ = 0.70), which suggests a wider feeding kernel gives a more stable system without sacrificing the productivity of the fish community.

A lower, increasingly negative m-value results in a zooplankton community with a flatter abundance spectrum. In other words, increasing herbivory results in a higher abundance in the larger zooplankton size classes. For the fish community, a shallower zooplankton spectrum leads to a higher abundance in the smallest fish size classes. The overall average slope of the fish community spectrum is similar across the 6 plots (**Figure 6**). This suggests the average slope of the fish community spectrum depends more on the feeding characteristics of the fish, over the dynamics of the zooplankton community. The total average fish biomass increases with increasing phytoplankton

FIGURE 4 | (A–C) The total fish community biomass (g m−3 ) and (D–F) the fish to zooplankton biomass ratio (F:Z) for different values of zooplankton feeding mode (m) , feeding kernel width (σz) and assimilation efficiency (Kz). In each plot, the other feeding parameters not specified are held constant at m = 0, σz = 0.75 and Kz = 0.7. The dashed line in (D–F) indicates the F:Z in the base model, where the zooplankton community are parameterized as fish.

abundance, across all six systems (**Figure 7A**). The magnitude of the increase in fish biomass correlated with the F:Z and CV of the system (**Table 4**). More fish were associated with a higher F:Z, and lower CV. The general zooplankton community system had the highest fish to zooplankton biomass ratio (4.43) and had the largest increase in total fish abundance: an 800% increase in fish. In contrast, the herbivorous copepod and chaetognath systems had similar fish to zooplankton biomass ratios to the general community (4.40 and 3.78), but higher CV's (1.26 and 1.27). These systems' fish biomass

TABLE 4 | Fish community biomass (FB), fish to zooplankton biomass ratio (F:Z), fish community production to biomass ratio (P:B) and throughput (TP), the variation in fish community biomass (coefficient of variation; CV) and the maxmium real part of the Jacobian (λmax) when the zooplankton community is defined by the feeding characteristics of different functional groups.

The system has a stable steady state when λmax < 0, and an unstable, oscillating steady state when λmax > 0.

increased by 340 and 410%, respectively. The flagellate system had the lowest fish to zooplankton biomass ratio (1.73), the second lowest CV (0.11) and the lowest increase in total fish biomass (170%).

Fish communities supported by herbivorous zooplankton communities were more resilient to fishing pressure, compared to fish supported by more carnivorous zooplankton (**Figure 7B**). The salp system had a negligible decline in average fish biomass, and chaetognath, herbivorous copepod and general community systems declined by up to 1, 2, 5%, respectively, with increasing fishing pressure. The two systems with carnivorous zooplankton communities (flagellates and carnivorous copepods) had an almost identical relationship between total rFB and fishing pressure, with both losing up to 15% of their average unfished biomass.

### DISCUSSION

This study is the first qualitative assessment of how zooplankton feeding characteristics mediate the transfer of energy from phytoplankton to higher trophic levels with a dynamic sizespectrum model. Improving the realism of the zooplankton community with zooplankton-specific feeding parameters increased the transfer efficiency of the system and the total mean biomass of the fish community, but changed the steady state of the system from a stable linear spectrum, to a series of traveling waves of abundance from smaller to larger size classes (**Table 3**; **Figure 3**). The change in steady state came from updating the zooplankton community PPMR and feeding kernel width (σZ). The general zooplankton community had a m-value of 0, which corresponds to a log10 PPMR of between 3 and 5 across the size range of the zooplankton community, and σZ of 0.75. This is in contrast to the fish community log10 PPMR of 2, and feeding

kernel width of 1. This observed change in the steady state agrees with the observed effects of increasing PPMR and decreasing σZ for fish communities (Blanchard, 2008; Law et al., 2009; Datta et al., 2011; Zhang et al., 2013).

Our results suggest a trade-off mediated by the zooplankton community, between the stability of the overall system and the total biomass and productivity of the fish community. A zooplankton community with a more generalist, carnivorous feeding strategy—defined by a lower PPMR (larger, positive m) and a wider feeding kernel—stabilized the steady state of the system (**Figure 5**), but the fish community was less abundant and productive (**Figure 4**). In contrast, a zooplankton community characterized by specialized, herbivorous behavior—defined by a higher PPMR (larger, negative m) and a narrower feeding kernel—increased the total average biomass and productivity of the fish community (**Figure 4**), but destabilized the system steady state (**Figure 5**). Herbivorous and mixed communities (m ≤ 0) with a narrower σZ had a higher ratio of fish to zooplankton biomass (**Figures 4D,E**), indicating a more efficient transfer of biomass from zooplankton to fish. This positive relationship between the zooplankton community PPMR and transfer efficiency corroborates with previous theoretical (Andersen et al., 2009) and empirical work (Jennings et al., 2002; Barnes et al., 2010); a higher PPMR yields a higher transfer efficiency between trophic levels, and fewer trophic levels separating phytoplankton from fish.

Zooplankton communities with a higher PPMR and narrower σZ had more variance in their abundance (**Figure 5**). These results suggest the abundance of zooplankton communities characterized by specialized herbivorous feeding behavior could exhibit more variation in their abundance than carnivorous communities. A similar relationship for fish species was found by Blanchard (2008), who established a link between the variation in fisheries catch of certain species of fish with their PPMR and σ; species with a higher PPMR and narrower feeding kernel had greater variability in their fishing catch through time. Further, Jennings and Warr (2003) identified a link between environmental stability and a smaller ecosystem average PPMR, which in this context means increasing herbivory amongst zooplankton in unstable environments. Such a relationship has been observed in marine ecosystems; herbivorous zooplankton dominate in unstable coastal and upwelling regions, whereas more carnivorous zooplankton are abundant in the open ocean (Raymont, 1980).

The resilience of the fish community to fishing pressure increased, and ecosystems became more efficient in mediating energy from phytoplankton to fish, when zooplankton communities had a larger σZ and higher PPMR characteristic of more herbivorous functional groups (**Figure 7**). The relationship between zooplankton community feeding characteristics, and the resilience of the fish community and ecosystem transfer efficiency, has potential implications for the marine environment under climate change. The world ocean's oligotrophic regions are expected to expand as a result of climate change (Doney et al., 2012; Sarmiento et al., 2004; Polovina et al., 2008). Food chains in warmer, oligotrophic oceans are traditionally believed to be longer than other regions, as a result of the dominance of smaller phytoplankton (Sprules and Munawar, 1986; Irwin et al., 2006; Morán et al., 2010), which would result in lower rates of energy transfer from primary producers to higher trophic levels. Further, recent studies suggest possible climate-driven shifts in the dominance of certain zooplankton functional groups, such as salps or jellyfish (Atkinson et al., 2004; Richardson et al., 2009; Schofield et al., 2010). Our results indicate that, everything else being equal, an increase in the dominance of carnivorous zooplankton groups could further decrease the transfer efficiency of expanding oligotrophic regions. Conversely, an increase in the abundance of herbivorous groups with a large PPMR, such as salps or herbivorous copepods, could decrease the number of trophic levels between phytoplankton and fish and thereby increase the transfer efficiency of these future oligotrophic regions.

Overall, increasing the zooplankton community PPMR had the greatest effect on increasing the total abundance, productivity and resilience of the fish community (**Figures 4**, **7**), and increasing σZ had the greatest stabilizing effect on the steady state of the system (**Figure 5**). Zooplankton have a higher average PPMR and smaller σZ in comparison to average observed values for fish and this difference has enormous implications for ecosystem transfer efficiency and stability (Barnes et al., 2010). This means that zooplankton feeding characteristics in particular PPMR and feeding kernel width—are a critical component to consider moving forward in how the transfer of energy from primary production to higher trophic levels is resolved in marine ecosystem models. This agrees with Jennings and Collingridge (2015), who suggest that a poor understanding of energy transfer in lower trophic levels is a potential cause for the order of magnitude discrepancy between model predictions and observed mesopelagic fish biomass over large spatial scales (Davison et al., 2013; Irigoien et al., 2014).

The large changes in fish biomass and productivity as a result of changes in the zooplankton community lead us to assess the implications of assuming a constant phytoplankton abundance spectrum within the model. In this study, we assume no feedbacks on the phytoplankton community from zooplankton (i.e., predation), however we know from empirical studies that the slope of the phytoplankton spectrum does change. The phytoplankton spectrum is shallower in eutrophic, upwelling systems—indicating a higher abundance of larger individuals such as diatoms—and steeper in oligotrophic systems where small-celled phytoplankton dominate (Sprules and Munawar, 1986; Irwin et al., 2006). The effects of eutrophy or oligotrophy on higher trophic levels could be investigated by varying not only the intercept, but also the slope of the phytoplankton community and incorporating feedback from zooplankton predation.

Our model did not investigate how changes in the body composition of different zooplankton functional groups affects energy transfer from phytoplankton to fish. Gelatinous zooplankton have around one-tenth of the carbon content per unit of live mass compared to non-gelatinous plankton (Kiørboe, 2013) and carbon content as a proportion of weight scales isometrically with increasing body size for carnivorous zooplankton (e.g., ctenophores and cnidarians), but decreases for filter feeders such as salps (Molina-Ramírez et al., 2015). This would have implications for the nutritional value of different zooplankton groups for the fish community, and the fish community's resultant growth rates. Future work could investigate the effect zooplankton body composition might have on energy transfer, by varying the assimilation efficiency of the fish community for different zooplankton functional groups.

Looking forward, theoretical studies have shown that including more traits than just individual body size increases the stability of the size-spectrum (Datta et al., 2011; Zhang et al., 2013) and improves the realism of modeled predator-prey dynamics (Boukal, 2014). Recent developments in dynamic sizespectrum theory now allow multiple functional groups and even species to be resolved within the community spectrum (Maury, 2010; Hartvig et al., 2011; Scott et al., 2014) and have been used to represent actual fish communities with increasing realism (e.g., Blanchard et al., 2014; Dueri et al., 2014; Spence et al., 2016; Zhang et al., 2016). The instabilities in our single-spectrum zooplankton community indicate that more complexity is needed if we are to represent realistic zooplankton communities within a dynamic size-spectrum framework. We envision the next steps toward this goal would involve a functional group approach, where the unique size-based characteristics of multiple sizebased zooplankton communities are represented, and the model is calibrated and compared with real-world data. The growing literature on the size-based behavior of zooplankton functional groups—coupled with the recent theoretical developments in dynamic size-spectrum modeling—means size-spectrum models that realistically resolve both zooplankton and fish may now be within reach.

### CONCLUDING REMARKS

In present size-spectrum model formulations focused on fish, small zooplankton are lumped together with phytoplankton in a background resource spectrum, and large zooplankton are represented as fish. The results of this study clearly demonstrate what we already know to be true: zooplankton are not fish, and nor are they phytoplankton. Current formulations that do not resolve the unique feeding characteristics of zooplankton are neglecting a significant factor in how energy is transferred from phytoplankton to fish. The results of this study motivate further work toward increasing the realism of zooplankton processes in size-spectrum models, and end-to-end marine ecosystem models more broadly.

### AUTHOR CONTRIBUTIONS

All authors were involved in conceiving the original idea for this study. JB provided code from past size-spectrum modeling studies. RH undertook the literature review to obtain zooplankton-specific feeding parameters, constructed the model, conducted the analysis, and wrote the manuscript, with input from AR, JB, and JE.

### ACKNOWLEDGMENTS

We thank Kai Wirtz for being available to answer questions about his mechanistic zooplankton feeding equations. We also thank Iain Suthers for helpful comments about the assumptions underlying our model, from a biological perspective. This study was supported by an Australian Research Council Discovery Grant (DP150102656).

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmars. 2016.00201

### REFERENCES

**Conflict of Interest Statement:** 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.

Copyright © 2016 Heneghan, Everett, Blanchard and Richardson. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Impacts of Intraguild Predation on Arctic Copepod Communities

#### Karolane Dufour 1 \*, Frédéric Maps 1 , Stéphane Plourde2 , Pierre Joly 2 and Frédéric Cyr 3

1 Takuvik Joint International Laboratory, Université Laval (Canada) – Centre National de la Recherche Scientifique (France), Québec-Océan and Département de Biologie at Université Laval, Québec, QC, Canada, 2 Institut Maurice-Lamontagne, Department of Fisheries and Oceans Canada, Mont-Joli, QC, Canada, 3 Aix-Marseille Université, Université de Toulon, Centre National de la Recherche Scientifique/INSU, IRD, Mediterranean Institute of Oceanography, UM 110, Marseille, France

Communities of large copepods form an essential hub of matter and energy fluxes in Arctic marine food webs. Intraguild predation on eggs and early larval stages occurs among the different species of those communities and it has been hypothesized to impact its structure and function. In order to better understand the interactions between dominant copepod species in the Arctic, we conducted laboratory experiments that quantified intraguild predation between the conspicuous and omnivorous Metridia longa and the dominant Calanus hyperboreus. We recorded individual egg ingestion rates for several conditions of temperature, egg concentration, and alternative food presence. In each of these experiments, at least some females ingested eggs but individual ingestion rates were highly variable. The global mean ingestion rate of M. longa on C. hyperboreus eggs was 5.8 eggs ind−1 d −1 , or an estimated 37% of M. longa daily metabolic need. Among the different factors tested and the various individual traits considered (prosome length, condition index), only the egg concentration had a significant and positive effect on ingestion rates. We further explored the potential ecological impacts of intraguild predation in a simple 1D numerical model of C. hyperboreus eggs vertical distribution in the Amundsen Gulf. Our modeling results showed an asymmetric relationship in that M. longa has little potential impact on the recruitment of C. hyperboreus (<3% egg standing stock removed by IGP at most) whereas the eggs intercepted by the former can account for a significant portion of its metabolic requirement during winter (up to a third).

#### *Edited by:*

Susanne Menden-Deuer, University of Rhode Island, USA

#### *Reviewed by:*

Jose M. Riascos, Universidad del Valle, Colombia Øyvind Fiksen, University of Bergen, Norway

> *\*Correspondence:* Karolane Dufour<EMAIL_ADDRESS>
#### *Specialty section:*

This article was submitted to Marine Ecosystem Ecology, a section of the journal Frontiers in Marine Science

*Received:* 29 July 2016 *Accepted:* 09 September 2016 *Published:* 23 September 2016

#### *Citation:*

Dufour K, Maps F, Plourde S, Joly P and Cyr F (2016) Impacts of Intraguild Predation on Arctic Copepod Communities. Front. Mar. Sci. 3:185. doi: 10.3389/fmars.2016.00185 Keywords: intraguild predation, copepods, *Metridia longa*, *Calanus hyperboreus*, Arctic regions, numerical modeling

### INTRODUCTION

Arctic and subarctic marine food webs are characterized by the presence of large calanoid copepods that channel primary production toward secondary consumers. Copepods have developed life cycle strategies that allow them to thrive in these highly seasonal environments. In the Arctic Ocean, several species of Calanus dominate the mesozooplankton biomass (Head et al., 2003; Hopcroft et al., 2010). During the short productive period in spring and summer, copepods feed on ice algae (when available in ice-covered regions) and phytoplankton and concentrate this energy into lipids, mostly stored as wax ester in hypertrophied oil sacs (Lee et al., 2006). In winter, they survive thanks to those lipid reserves that fuel their reduced metabolism during an extended period of dormancy (the diapause). All year long, copepods are an essential food source for many predators such as bowhead whale, little auk, and above all Arctic cod, a cornerstone component of the Arctic food web (Fort et al., 2010; Falardeau et al., 2013; Pomerleau et al., 2014). Thus, communities of large copepods form a critical hub of matter and energy fluxes in Arctic and subarctic marine food webs.

In the Arctic, copepod community biomass is mainly composed of the large species Calanus hyperboreus and Calanus glacialis (adult female median prosome length of 6.7 and 4.1 mm, respectively), the medium-sized Metridia longa (2.8 mm) and the small Pseudocalanus spp. (1.1 mm) (Darnis et al., 2012). Following the current dynamics, the boreal species C. finmarchicus is also regularly found in marginal Arctic seas, especially in the eastern Greenland Sea and Barents Sea (Conover and Huntley, 1991). Moreover, C. finmarchicus' biogeographic distribution in the surface layer is projected to move even farther northward in response to surface circulation and temperature forcing induced by climate change (Reygondeau and Beaugrand, 2011). The structure and functions of copepod communities are critical from an energetic point of view for marine ecosystems and depend on the actual assemblage of species. From one species to the other, the oil sac size is different and therefore the energetic content differs. C. hyperboreus and C. glacialis are bigger and contain more lipids than the boreal C. finmarchicus. Consequently, many small visual predators such as fish larvae and juveniles, would reap a larger energetic reward for a similar harvesting effort, making these Arctic copepod species the preys of choice. Hence, changes in the assemblage of copepods communities could impact marine predators' recruitment, resulting in a form of bottomup control (Mills et al., 2013; Greene and Pershing, 2014). For example, along the West coast of the Spitsbergen island in the Svalbard archipelago, the increase of warm Atlantic water masses that bring along abundant C. finmarchicus may have a negative impact on the reproductive success of little auks by reducing the relative abundance of its preferred prey, C. glacialis (Kwasniewski et al., 2010). In order to understand and predict impacts of environmental changes on Arctic marine ecosystems, it is necessary to better understand the mechanisms responsible for the specific assemblages of copepod communities.

In addition to physical forcing, relationships occurring within copepod communities can influence their structures and functions. Species that follow each other and co-occur in a community not only share food and space resources, but also develop complex interactions between them. Intraguild predation (IGP) has been proposed as an ecological strategy that could structure copepod communities (Irigoien and Harris, 2006; Plourde et al., 2009; Darnis, 2013; Melle et al., 2014). This particular type of predation occurs between members of a group of species that share the same food resources (Polis et al., 1989). This phenomenon is widespread across a variety of marine and terrestrial ecosystems and through all taxonomic and trophic levels (Polis et al., 1989; Holt and Polis, 1997; Arim and Marquet, 2004). This complex interaction is particularly interesting since it results in an immediate energy gain for the predator as well as a long term reduction of its competition (Hiltunen et al., 2013).

In copepod populations, survival rate to adulthood is strongly influenced by egg and nauplii mortality (Davis, 1984; Plourde et al., 2009). Egg mortality is particularly high in broadcast spawning species that release their eggs in the water column (Ohman et al., 2004), such as Calanusspp and M. longa. Although predominantly considered to be herbivorous, most calanoid copepod species have a flexible diet and can be omnivorous or even cannibalistic (Landry, 1981; Ohman and Hirche, 2001; Bonnet, 2004; Basedow and Tande, 2006). These species usually generate a filtration current to obtain small prey, essentially phytoplankton and micro-zooplankton, or cruise through the water and attack when a prey is detected (Kiørboe, 2013). Thus, they can consume eggs and young nauplii stages with limited mobility.

Cannibalism on eggs and nauplii can control the phenology of C. finmarchicus recruitment (Ohman and Hirche, 2001). It is therefore likely that IGP can impact recruitment as well as influencing the temporal succession of dominant species (Irigoien and Harris, 2006). In environments as contrasted as the Beaufort Sea and the St-Lawrence estuary, reduced C. hyperboreus recruitment and abundance co-occurred with an increase in M. longa abundance (Plourde et al., 2002; Darnis, 2013). During the time of peak egg production by C. hyperboreus, the gut of individual M. longa is often observed to be orange, a particularity linked to the probable ingestion of lipid-rich C. hyperboreus eggs (Conover and Huntley, 1991). Hence it has been hypothesized that M. longa individuals that remain active at intermediate depths throughout winter (no diapause) could intercept and ingest buoyant C. hyperboreus eggs that are spawned in deep waters (Plourde et al., 2003; Darnis, 2013). Given that M. longa population dynamics resources are scarce in winter, lipid-rich eggs could represent an important energy source for such an omnivorous (opportunistic) species. Later in early spring, numerous C. hyperboreus and C. glacialis that emerge from diapause and ascend toward the surface ahead of the phytoplankton bloom could also be feeding on eggs and young nauplii stages. Consequently, IGP could impact the recruitment of the true Arctic C. hyperboreus.

Surprisingly, only few studies have been conducted on predation within calanoid copepod communities (Landry, 1981; Huntley and Escritor, 1992; Bonnet, 2004; Basedow and Tande, 2006; Vestheim et al., 2013). The majority of the studies targeted cannibalism and while none was focused on IGP between the dominant Arctic species M. longa and C. hyperboreus, Huntley and Escritor (1992) reported ingestion rates of the vicariant species M. gerlachei on eggs of the dominant Antarctic Calanoides acutus in the Austral Ocean. In all experiments, ingestion rates on eggs (and nauplii) vary according to their concentration. The influence of alternative food source on ingestion produced ambiguous experimental results. On one hand, experiments with female C. pacificus suggest a switch between herbivorous and carnivorous behavior that depends on the relative abundances of phytoplankton and their own nauplii (Landry, 1981). On the other hand, experiments with C. finmarchicus have shown that ingestion rates on its own nauplii are independent of the ambient algae concentration (Basedow and Tande, 2006).

Copepod community models have been developed with a focus on development, growth and competition for food resources (Record et al., 2012, 2013). However, community level processes involving interspecific interactions such as IGP have not yet been implemented into these models. In order to provide a better understanding of IGP and to provide a better parameterization for models of copepod communities, we conducted laboratory experiments that quantified the ingestion of copepod eggs by some of the dominant Arctic copepod species, M. longa and C. hyperboreus. We further conducted a simple numerical experiment in order to assess the potential ecological implications of our findings in Arctic marine ecosystems. Copepods were sampled in the Lower St-Lawrence estuary (LSLE) and feeding experiments on eggs were conducted under different conditions of temperature, egg concentration, and alternative food availability.

### METHODS

### Area of Study and Sampling

The LSLE is the southernmost sea directly influenced by Arctic water masses and that presents Arctic features in the North Atlantic (**Figure 1**). Arctic water masses enter in the Gulf of St-Lawrence (GSL) by the Strait of Belle Isle in the North, transporting Arctic copepods species such as C. hyperboreus, C. glacialis, and M. longa. Warmer North Atlantic water masses enter the system by Cabot Strait in the South and carry along boreal species such as C. finmarchicus. Both Arctic and temperate species thrive in the GSL system where cold and warm water masses are segregated vertically between the thin seasonal surface layer, the cold intermediate layer renewed locally each winter (core temperatures can be negative) and the deep Atlantic layer (between 4 and 5◦C). Colder water masses dominate the eastern and northern parts of the Gulf and the deep lower Estuary where upstream tidal pumping brings cold and nutrient-rich waters to the surface, whereas warm surface waters (up to 20◦C) develop throughout the summer over the shallow southern half of the Gulf (Saucier et al., 2003; Le Fouest et al., 2005). These particular water masses and species mix make the LSLE an exceptional experimental model to study IGP in the context of the rapid environmental changes that the Arctic is currently facing. Zooplankton was sampled in the LSLE (48◦ 40′N, 68◦ 35′W) in October 2014 with a vertical plankton net (158 µm mesh size, 1 m diameter) at 125 or 200 m from the surface at a towing speed of 0.5 m s−1 . The catch was diluted into 4 L jars filled with filtered seawater and maintained close to ambient temperature at 5 to 6 ◦C in coolers during transport to the laboratory (less than an hour).

In a subsequent numerical experiment, we applied our observations to a truly Arctic environment, the Amundsen Gulf. The Amundsen Gulf is located in the Canadian territory of Nunavut at about 71◦N and bordered by Banks Island and Victoria Island. It connects south-eastern Beaufort Sea to the Canadian archipelago (**Figure 1**). The Amundsen Gulf water masses are generally cold (around 0◦C) and are mainly discriminated by their different salinities: the Polar-Mixed Layer from the surface to c.a. 50 m (S ∼ 31.6), the Pacific Halocline below until 200 m (32.4 < S < 33.1) and the slightly warmer Atlantic Layer below 200 m (S ∼ 34) (Carmack and Macdonald, 2002). The Gulf covers about 60,000 km2 and the maximum

depth of this large channel is 630 m. In winter, it is entirely ice-covered except for sporadic polynyas and flaw leads and the spring ice-breakup his highly variable (Galley et al., 2008).

### Live Animal Sorting

Owing to the harsh conditions at sea during winter over the LSLE, copepods were sampled at the end of October 2014 prior to the beginning of C. hyperboreus reproduction and the formation of dense sea ice. Plourde et al. (2003) have previously shown, for a similar experimental setup, that capture and handling of C. hyperboreus females triggered gonad maturation, leading to egg production a few weeks after collection and approximately a month earlier than expected according to the in situ timing. A visual inspection of appendages and behavior under a binocular microscope allowed us to select adult female M. longa and C. hyperboreus in good condition from the live samples, within 48 h of the catch. Animals were kept in the dark between 3 and 5 ◦C in groups of 25 to 50 in 1 L beakers equipped with egg separators (mesh size = 333 µm) and filled with filtered seawater. Female M. longa were fed with solutions of concentrated diatoms Thalassiosira weissflogii (Instant Algae R TW 1200), whereas female C. hyperboreus where not fed since egg production is entirely fuelled by internal lipids during the dormant part of their life cycle. The beakers were inspected daily for egg production for about a month until female C. hyperboreus spawned enough eggs to start the experiments.

### Predation Experiments

Within 48 h of the start of the experiment, individuals used as predators were photographed laterally using a PixeLINK camera of 5 Mb (PL-E425CU) mounted on a LEICA MZ6 stereoscope. In order to minimize stress, individuals were kept in cooled seawater until the very moment that the picture was taken and gently manipulated with handling needles. Prosome length (distance between the tip of the cephalosome and the tip of the last thoracic segment), prosome area, oil sac length, and area were measured with the software ImageJ v. 1.49. The condition index was estimated as the oil sac area divided by the prosome area of the individual in order to give an indication of its lipid content. Predators' carbon content was estimated from species and stagespecific seasonal relationships between individual carbon mass and prosome length (Forest et al., 2010). The photographed animals were then placed individually in 45 mL Petri dishes filled with filtered seawater and equipped with egg separator (mesh size = 333 µm) for acclimation at the experimental temperatures.

Experiments were carried in an atmosphere-controlled chamber in November 2014 with female M. longa fed with C. hyperboreus eggs. Predators were placed individually in a bottle (1.35 or 2.37 L) filled with filtered seawater containing a precise number of eggs spawned within the previous 48 h. Bottles were turned upside down once per minute on a rotating wheel in order to maintain the eggs in suspension. Manipulations were conducted under a dimmed red light and the experiments were carried in the dark. Incubation time was kept relatively short (less than 24 h) in order to avoid a complete consumption of the eggs and hence a failure to accurately estimate ingestion rates. It also varied according to incubation conditions, with longer incubation times for lower egg concentration (lower encounter probability between eggs and the predator). After 4 to 22 h, the contents of the bottles were filtered with a 73 µm sieve and the remaining eggs were counted. A minimum of two controls without predator for each treatment was set up in order to check the accuracy of the egg recovery method. The status of each individual was verified at the end of the experiment. Ingestion rates were discarded for the few dead individuals, and the sluggish or unhealthy-looking ones as well.

We tested the effect of temperature on M. longa ingestion of C. hyperboreus eggs. Incubations were carried at 1, 4, and 8 ◦C to reflect the potential range of in situ water temperature encountered by this species between the Arctic and subarctic regions. Unfortunately, we could not test for negative water temperatures. In order to characterize functional responses, we tested the effect of C. hyperboreus egg concentration on M. longa ingestion rate. We were not aware of actual data about in situ C. hyperboreus egg concentrations in the water column, but we estimated it to be low (Huntley and Escritor, 1992). Hence we chose concentrations of 5, 10, 20, and 30 egg L−1 . Finally, we checked whether there was an influence of alternative food availability, by adding an additional food source in the form of concentrated diatoms T. weissflogii (Instant Algae R TW 1200) in half of our replicates. Algae were offered at about 50 µg C L −1 according to cell concentration determined with a Hausser Bright-Line Hemacytomer (couting chamber) and carbon to volume relationships for diatoms (Menden-Deuer and Lessard, 2000).

### Estimation of Ingestion Rates

The instantaneous feeding rate on eggs g (h−1 ) was derived from an exponential equation (Båmstedt et al., 2000):

$$\lg = \frac{\ln\left(\frac{E\_f}{E\_0}\right)}{t} \tag{1}$$

where E0 and Ef are respectively the egg concentration at beginning and the end of experiment (egg L−1 ) and t is incubation time (h). The number of C. hyperboreus eggs obviously did not increase during the experiments, and experiment duration was not long enough for hatching to occur (no nauplii was ever found in any of the control or experiment bottles).

We deduced the clearance rate F (L ind−1 h −1 ), which corresponds to the volume of water processed assuming 100% capture efficiency and a homogeneous food concentration, from both g and the volume of the incubation bottle V (L):

$$F = \text{g} \times V \tag{2}$$

Finally, we obtained the ingestion rate I (egg ind−1 h −1 ) with:

$$I = F \times \left[E\right] \tag{3}$$

where [E] is the average egg concentration as given by:

$$\mathbf{E}[E] = \mathbf{(E\_0 \times \frac{1 - e^{\left(-\mathbf{g} \times t\right)}}{\mathbf{g} \times t})} \tag{4}$$

Daily ingestion rates where assumed to be 24 times the hourly rates since M. longa is known to be active and swimming almost continuously (Hirche, 1987). Egg ingestion rates were then converted in carbon units using an egg carbon content of 0.84 (µg C egg−1 ; Plourde et al., 2003). The proportion of daily metabolic needs (%) satisfied by egg ingestion was estimated from the ratio of carbon ingestion rate IC = 0.84 ∗ I (µg C ind−1 d −1 ) and a mean and constant respiration rate (µg C ind−1 d −1 ) measured by Seuthe et al. (2006).

### Data Analysis

The number of ingested eggs followed a Poisson distribution (an asymmetric right-skewed distribution of discrete values). Hence, in order to minimize estimating errors, a generalized linear model (GLM) for Poisson distribution was used to predict the number of eggs ingested (EI) (the raw data) according to temperature (T), egg concentration (E), presence of additional food source (AC), prosome length (PL), and condition index (CI) of the predator. Bottle volume (V) and duration of experiment (D) were taken into account by using them in an offset term. Several models were tested (with and without interactions) and we computed Akaike's information criterion corrected for overdispersion (QAIC) as a decision-support metric. The GLMs formulae were of the form (here is the one with all the independent variables but no interactions):

$$\log\left(\frac{EI}{D \times V}\right) = \beta\_0 + (\beta\_T \times T) + (\beta\_E \times E) \tag{5}$$

$$+ (\beta\_{AC} \times AC) + (\beta\_{PL} \times PL) + (\beta\_{CI} \times CI)$$

where βi are the coefficient estimates for each variable. The number of eggs ingested per unit of time and volume can easily be obtained from this model results.

### Model of Egg Vertical Distribution

In order to assess the implications of our findings in the context of Arctic marine ecosystems, we developed a simple water column (1D) model of the vertical distribution of C. hyperboreus eggs in the Amundsen Gulf. The model computed the time evolution of egg concentration at a given depth according to advection and diffusion, gains by egg production and losses by development and predation (**Figure 4**). The rate of change of egg concentration followed the classical advection-diffusion-reaction formulation (Soetaert and Herman, 2009):

$$\frac{\partial E}{\partial t} = -\mathbf{w}\frac{\partial E}{\partial z} + K\frac{\partial^2 E}{\partial z^2} + P\_z - I \ast M + 1/H \ast E \tag{6}$$

where the first right-hand side term represents the effect of egg vertical velocity, the second the effect of diffusion and the others several biological reaction terms. More specifically, E was the egg concentration (egg m−3 ), t the time (h), w the egg velocity (m h−1 ) defined positive downward, z the depth (m), K the vertical eddy diffusivity coefficient (m2 h −1 ), PZ the depth-dependent egg production rate (egg m−3 h −1 ), H the egg hatching time (h), M the females M. longa concentration (ind m−3 ) and I the ingestion rate by other copepods (egg ind−1

h −1 ). The egg velocity w was given by Stokes' law (Visser and Jónasdóttir, 1999):

$$\omega = \Re 600 \times \frac{\text{g} \times d^2 \times (\rho\_{\text{egg}} - \rho\_{\text{water}})}{18 \times \mu} \tag{7}$$

where g is the gravitational constant (9.81 m s−2 ), d is the egg diameter (m), ρegg is the egg density (g m−3 ), ρwater is the water density (g m−3 ) and µ is the dynamic viscosity of the seawater (1.85 g m−1 s −1 ), here taken as a constant (**Table 1**).

The source term for eggs (PZ) came from the average daily production of C. hyperboreus population (30,000 eggs m−2 ) observed between February and April 2008 in the Amundsen Gulf (Darnis, 2013). C. hyperboreus females released more than 90% of their eggs during this 3 month-period, while remaining at depths between 200 and 300 m (Darnis, 2013). As a result, we computed the vertical profile of egg production rate PZ according to a normal distribution whose mean was centered at 250 m, its standard deviation 15 m and its integral equal to 30,000 eggs m−2 (99.9% of the eggs were released between 200 and 300 m).

The IGP rate I (h−1 ) exerted on C. hyperboreus eggs was simply the product of M. longa females' abundance (ind m−3 ) and the individual filtration rate (m3 ind−1 h −1 ) found in our grazing experiments.

Egg hatching time H (h) followed an empirical Belehrádek's ˇ function obtained from observed hatching times of C. hyperboreus eggs (Jung-Madsen et al., 2013):

$$H = a \times \left(T - \alpha\right)^{-b} \tag{8}$$

where T was water temperature (◦C), a (d ◦C −1 ), α ( ◦C) and b constants.

A simple ordinary differential equation framework was not optimal for the modeling of developing eggs because of the "numerical diffusivity" caused by the hatching rate (Gentleman et al., 2008). Simply put, with a development (hatching) rate, the progression through development stages is treated as a continuous process within the population, instead of a discrete event highly synchronized among many individuals. This can

#### TABLE 1 | Model parameters and references.

VJ1999, Visser and Jónasdóttir (1999); D2013, Darnis (2013); JM2013, Jung-Madsen et al. (2013)

lead to unrealistic and spurious results, such as a small fraction of the simulated egg population that has already hatched after the first time step! In order to prevent this, we used the simple approach of spreading the egg development over 20 numerical stages of equal length (Gentleman et al., 2008) and we further integrated the development throughout theses stages with a flux limiting numerical scheme (Record and Pershing, 2008).

### Simulations

The 1D water column model was split into 5 m vertical layers between 0 and 300 m, and the time step of integration was 12 h. For the model forcing, we used physical and biological datasets from the Circumpolar Flaw Lead System Study (CFL; Barber et al., 2011). We only selected profiles from stations that were at least 300 m deep and located within the Amundsen Gulf. As a result, the physical forcing came from mean vertical profiles of eddy diffusivity coefficient (K), temperature (T), and water density (ρwater) obtained at 20 different stations. These stations were sampled in November and December 2007 under dense sea-ice cover with a vertical microstructure profiler (VMP500, Rockland Scientific International). For two long-term stations the number of casts used in the average profile were respectively 24 and 25, whereas at least 5 profiles were used to build the 18 other mean profiles. Missing values in averaged profiles (e.g., near the surface or below 250 m) were dealt with according to Equation (1) from Bourgault et al. (2011). Further details about this dataset can be found in their study.

We also used 18 vertical profiles of M. longa female abundance obtained with a Hydrobios R multinet sampler. Details of the sampling procedure can be found in Darnis and Fortier (2014). We selected the stations sampled between February and April 2008, during the peak of C. hyperboreus reproduction. The layer thickness for the vertical sampling ranged from 10 to 144 m, with a median of 20 m.

In order to test the sensitivity of the model to both the physical properties of the water column and the vertical distribution of the predators (M. longa) we ran 360 simulations, one for each possible combination of the physical and biological forcing fields. In addition, we ran this ensemble of simulations for three different scenarios of egg density in order to verify the impacts of different egg velocities: the mean (scenario D0), minimum (Dmin) and maximum (Dmax) egg densities observed by Jung-Madsen et al. (2013). Simulations ran for 15 days in order to reach a quasi-steady state where the maximum local rate of change in egg concentration δE/δt was less than 10−6 .

### RESULTS

### Egg Ingestion Rates

Predation of C. hyperboreus eggs by female M. longa occurred in each of the incubation experiments. However, in each experiment there was high individual variability and several individual incubations showed no egg ingestion. After discarding incubations within which dead or unhealthy individuals were found at the end of the experiment, for each treatment approximately 8 individual replicates out of the initial 10 were used for further analyses. The frequency distribution of ingestion rates was positively skewed, i.e., the median value was lower than the mean. The global mean ingestion rate of C. hyperboreus eggs by M. longa was 5.8 eggs ind−1 d −1 (SE = 0.57, n = 141) and the median was 3.7 eggs ind−1 d −1 (**Figure 2**). In terms of carbon, the mean ingestion rate was 4.9 µg C ind−1 d −1 and the median 3.1 µg C ind−1 d −1 (SE = 0.48, n = 141; **Figure 2**). C. hyperboreus eggs constituted an energy-rich food source (0.84 µg C egg−1 ; (Plourde et al., 2003) and the average daily ration of M. longa females feeding on C. hyperboreus eggs was 37% of their estimated metabolic needs based on respiration rates (SE = 4, n = 141). Individual variability resulted in a contrasted pattern where about a quarter of M. longa did not ingest any eggs, while an equivalent proportion filled more than 50% of their daily energetic requirements through egg ingestion. Some individuals actually largely exceeded their daily metabolic needs (i.e., over 100%, **Figure 2**), even when eggs were offered at low concentrations.

### Influence of Incubation Conditions

According to the GLM analysis, presence of additional food source, prosome length, or condition index of the individuals had no discernable effect on egg ingestion (**Table 2**). The only predictor that improved the model and that significantly affected the number of eggs ingested by M. longa was egg concentration (**Figure 3**). In the full model with no interactions, temperature seemed to have a significant negative effect on ingestion rates (**Table 2**). However, in the model that only kept egg concentration and temperature as independent variables, the influence of temperature and the interaction term between them did not remain significant. Moreover, the size effect of the temperature and interaction term coefficients was small compared to the impact of egg concentration and the QAIC values were very close between the models. Hence, we decided to use the simplest model with only egg concentration as predictor of egg ingestion (**Figure 3**). We did not observe any feeding saturation for the range of egg concentrations offered.

### Simulated Egg Vertical Distribution

The velocity of C. hyperboreus eggs estimated with Equation (7) ranged between −9 and −4.1 m d−1 for the mean egg density scenario (D0) over the 20 physical forcing profiles. Negative velocities meant that eggs were positively buoyant, from the bottom of the water column up to the surface. Egg velocity was dependent on the water density profile and it decreased slowly with decreasing depth (**Figure 4**). This pattern was conserved among the 20 physical forcing profiles whose overall variability was low. Egg velocity was strongly influenced by egg density itself. C. hyperboreus egg density is highly variable, both between individual females and within clutches of the same female (Jung-Madsen et al., 2013). When we used the minimal egg density (scenario Dmin) the velocity tripled to range between −26.2 and −21.3 m d−1 , whereas for maximal egg density (scenario Dmax) the associated velocity was not negative throughout the water column. The egg velocity ranged between −2.3 and 2.6 m d−1 with a converging depth of neutral buoyancy around 100 m.

Egg velocity was critical for egg vertical distribution. In the ensemble of simulations for scenario D0, egg vertical

distributions showed higher concentrations between 150 and 275 m and peaked around 225 m to reach about 3000 eggs m−3 (or 3 eggs L−1 ) once the simulation reached its steady state (**Figure 4**). For scenario Dmax, the denser eggs ascend only slightly in the water column before hatching. They remained concentrated between 200 and 300 m with a maximum egg concentration a little less than 4000 eggs m−3 near 250 m (**Figure 4**). Contrarily, in the Dmin scenario, eggs moved rapidly upward and some even managed to reach the first 5 m of the water column to attain a concentration of a little less than 300 eggs m−3 (**Figure 4**). Eggs where spread over the entire column and concentration peaked above 160 m at about 1200 eggs m−3 .

### Impact of *M. longa* Predation on *C. hyperboreus* Eggs

For each egg density scenario tested, the proportion of egg biomass eaten by M. longa was more sensitive to the profiles of M. longa abundance than to the physical forcing (see **Figure 5** for scenario D0). Over the ensemble of 360 simulations of the D0 scenario, the percentage of C. hyperboreus egg standing stock ingested by M. longa ranged between 0.1 and 1.1% with most of the values being in the low end (**Figure 6**). The percentage of egg standing stock ingested by M. longa varied between 0.2 and 3.2% for scenario Dmin and between 0 and 0.8% for maximal egg density (scenario Dmax; **Figure 6**).

This modest impact of M. longa IGP is hardly noticeable between the egg concentration profiles simulated with and without egg predation, even for the combination of physical conditions and M. longa profiles that lead to the maximum difference (**Figure 7**). The corresponding daily egg ingestion of the whole population of M. longa ranged between a little more than 1 to almost 20 eggs m−3 d −1 (**Figure 7**). It is noticeable, though, that such egg ingestion values could allow M. longa individuals to satisfy almost 10% of their metabolic TABLE 2 | Results from the generalized linear models (GLM) fitted to predict the number of *C. hyperboreus* eggs ingested by female *M. longa* according to egg concentration (*E*), temperature (*T*), alternative food source (*AC*), condition index (*CI*), prosome length (*PL*), and interactions between the terms.

Estimated coefficients b, standard error SE, P-values (\*\*\*P < 0.001;\*\*P< 0.01; \*P < 0.1; ns not significant) and Akaike's information criterion corrected for overdispersion QAIC. Model entry in bold indicates the one selected for further numerical experiments.

needs according to respiration rates from Seuthe et al. (2006), or up to 37% if we consider the lower respiration rates reported by Hirche (1987).

### DISCUSSION

### Individual Variability

Our results show high individual variability of egg ingestion rates, with a quarter of all the individuals not ingesting any eggs and about the same proportion satisfying more than half of

FIGURE 3 | Ingestion rates of female *M. longa* on *C. hyperboreus* eggs according to the presence of alternative food source (left panel), temperature (middle panel), and egg concentration (right panel). Points are jittered to reduce overlap. Black and white line is the prediction of the mean from the selected generalized linear model (see text); dotted black lines are the corresponding 95% confidence intervals.

simulated for mean egg density (scenario D0) of 19.4 kg m−3 are presented in the right panel. Dashed line is the egg velocity (m d−1 ) computed according to Stokes' law. Dotted line is at the initial condition corresponding to the spawned eggs profile. Solid line is the egg concentration distribution after 15 days of simulation. Gray area is the female M. longa abundance (ind m−3 ).

their daily metabolic needs from egg grazing. This asymmetric and widely spread distribution has a geometric coefficient of variation of 149%. This pattern is consistent across the range of incubation conditions we tested, and this level of variability is common in any experimental setting measuring individual biological features (size, structural or storage weight, swimming behavior, etc.) and physiological rates (respiration, ingestion, growth, etc.) (e.g., Basedow and Tande, 2006). We do not have clear explanations for the level of variability we observed, but it is a useful observation to report (see Supplementary Material for a spreadsheet of individual observations). Individual variability has long been recognized as a key property of plankton ecology since population dynamics and trophic interactions (that are of primary interest for marine ecologists) are emerging properties of individual characteristics and behaviors (Båmstedt, 1988). Modern experimental and in situ observation methods

are providing increasingly detailed and abundant individuallevel data (e.g., Schmid et al., 2016), while current numerical approaches allow for testing how and how well individual-based models can effectively represent emerging properties at higher organizational levels (Neuheimer et al., 2010; Morozov et al., 2013).

### Dynamical Interactions in the Water Column

This modeling exercise provided insight on the physical and biological dynamical processes interacting in the water column and their relative importance for the vertical distribution of C. hyperboreus eggs.

We first performed a scale analysis of the two first terms on the right-hand side of Equation (7), i.e., egg buoyancy and vertical turbulence. For the simulation scenario illustrated in **Figure 4**, we can estimate the vertical gradient in egg concentration ∂E ∂z ≈ 15 egg m−4 near the maximum egg concentration at 250 m (an approximate increase of 750 egg m−3 over 50 m) and ∂ 2E ∂z 2 ≈ 1 egg m−5 ( ∂E ∂z varies between ± 15 egg m−3 over about 30 m from both sides of the maximum concentration), simple arithmetic suggests that

and

$$K \frac{\partial^2 \mathbf{E}}{\partial \mathbf{z}^2} \approx 10^{-5} \text{egg } \text{m}^{-3} \text{ s}^{-1}$$

egg m−3

s −1

≈ [1, 3] × 10−3

with w ε [4, 9.3] m d−1 (see Results) and K = 3.4×10 m−6 s −1 (background turbulent diffusivity) from (Bourgault et al., 2011). This suggest that from a physical point of view, the buoyant vertical displacement of eggs is dominant over the turbulent diffusion mechanism, and only turbulent events about a 100 times above the background value could effectively influence their distribution. The role of turbulent mixing in egg distribution (aggregation or spreading) is likely minimal and the use of a parameterization such as the one presented in Bourgault et al. (2011) could have been sufficient here. The effect of turbulence may only become important for denser eggs rising very slowly toward the surface or directly after the spawning if it occurs in a thin layer pattern, hence producing a high concentration gradient. This effect could be studied more efficiently with new in situ sampling devices such as the LOKI underwater imaging system that can provide highly resolved vertical distribution of adult females C. hyperboreus, their eggs, a whole suite of potential other intraguild predators beyond M. longa and the physical properties of the water column as well (Schmid et al., 2016).

From a biological point of view, egg density had an overwhelming impact on egg vertical distribution patterns (**Figure 7**). Egg density defined their vertical velocity w and as a result both the range of depth they could reach before hatching and their corresponding concentration. For two out of three egg density scenarios, eggs laid at depth did not manage to reach and accumulate within the surface layer. Even for the minimum density scenario Dmin, the amount of eggs reaching the surface remained marginal. This is coherent with the generally accepted idea that C. hyperboreus nauplii, rather than eggs, accumulate under the ice in advance of the phytoplankton bloom (Conover and Huntley, 1991). Meanwhile, the maximum egg concentration reached was about 4 eggs L−1 , close the minimal egg concentration used in our grazing experiments. Depth and concentration were crucial for the interaction with female M. longa whose vertical position and abundance vary a lot, and whose ingestion rate depends on the surrounding egg concentration. As a result, it appeared that the probability of encounter between a predator and an egg of C. hyperboreus was determined essentially by the density of the latter.

### Impact of Intraguild Predation on *C. hyperboreus* Recruitment

If we consider thin layer effects to remain marginal, our results suggest that impacts on C. hyperboreus population dynamics may remain limited in space and time since M. longa IGP was limited to a little more than 3% of C. hyperboreus egg biomass. However, young nauplii stages could also be preyed upon by M. longa and thus our figures could underestimate the actual impact of M. longa on C. hyperboreus recruitment. M. longa is a cruising feeder, i.e., it cruises through water searching for food and prey and captures them upon detection. Moreover, motile preys such as nauplii can generate a hydrodynamic trail while swimming that could render them be easier to detect by M. longa than nonmotile preys such as eggs (Kiørboe, 2011). This assumption is supported by feeding experiments with M. lucens and M. longa in which phytoplankton and much larger nauplii Artemia were offered together, illustrating selective feeding on Artemia nauplii (Haq, 1967).

In addition to females M. longa, other development stages of this species and other copepod species could also exert IGP on

w ∂E ∂z

(scenario Dmax) (right panel).

C. hyperboreus eggs. As already mentioned, adults and advanced copepodite stages of the large C. glacialis emerge from diapause and initiate ascent to surface layers prior the spring bloom to feed on ice algae (Daase et al., 2013). Therefore, they could also benefit from the energy rich C. hyperboreus eggs (and potentially the nauplii) and add to the predation pressure. Thus, we likely have underestimated IGP pressure on C. hyperboreus and its impact on its recruitment in several ways. However, it may remain very dependent on the seasonal timing of the feeding activity of these potential predators, as well as their vertical position in the water column as already demonstrated.

It seems that there is actually an ecological trade-off for C. hyperboreus females to lay low (lipid-rich) or high-density (lipid-poor) eggs. On one hand the obvious metabolic advantage for lipid-rich eggs is that the offspring can rely on abundant reserves to develop. The downside however could be that low-density eggs reach faster and "en masse" the layer where M. longa and other species are abundant, hence increasing their mortality risk. On the other hand, dense eggs ascend very slowly and actually hatch before having reached these dangerous depths. The wide range in egg density that has been observed among eggs from different C. hyperboreus females but also within the same egg clutch (Jung-Madsen et al., 2013) could actually be part of a strategy that mitigates predation risk by spreading the eggs across a range of ascent speeds. However, the delicate balance of both opposing effects on offspring fitness could only be assessed in a more detailed modeling study supported by finely resolved and concurrent observation of the zooplankton community vertical distribution.

### *C. hyperboreus* Importance in the Arctic System

Our results could offer an interesting contrast with the earlier study of Huntley and Escritor (1992) on a couple of homologous species from Antarctica, M. gerlachei and C. acutus. They showed from incubations of a group of individuals that M. gerlachei could reach daily rations ranging between 4 to 11% of its body weight when offered 1000 egg L−1 . The authors further estimated that this concentration was about three orders of magnitude higher than what it could be in situ (1 egg L−1 ) and thus concluded that this type of predation is likely insignificant. However, our experimental results showed that M. longa individuals could meet up to 75% of their daily metabolic requirements on eggs of C. hyperboreus offered at a concentration of 5 L−1 , based on respiration rates observed in Amundsen Gulf (Seuthe et al., 2006). We also simulated that the mean in situ ration of eggs for M. longa females should vary between 8 and 37% of their metabolic needs, depending on the respiration rate. The lower boundary is based on Seuthe et al. (2006), whereas the upper boundary is based on lower respiration rates (Hirche, 1987). Neither value is sufficient to ensure complete metabolic maintenance, but during wintertime it could be combined to some level of lipid stores to cope with the otherwise scarce food background available to these copepods. Lipid-rich C. hyperboreus eggs (and potentially nauplii) are a reliable and valuable food source that at least some M. longa individuals seem prone to take advantage of during several months of the year.

Our study stresses the importance of C. hyperboreus as a linchpin of Arctic marine trophic network. This key species has adapted to the extreme environment by evolving its ability to store considerable amounts of lipids. C. hyperboreus efficiently concentrate the energy from the short-lived primary production bloom and supports the entire marine trophic network. Higher trophic levels rely heavily on its large copepodite and adult development stages, while zooplankton species from a similar trophic level (from the same guild) also benefit from the smaller packaging of the bounty within its eggs. Moreover, since C. hyperboreus reproduction occurs for several months during a period of the polar year when primary production has shut down, it likely provides a precious and unparalleled resource for several species of planktonic predators that remain active in the dead of winter. Hence, IGP should not only be considered as an extra mortality source that could affect recruitment of the species that is preyed upon, but also as a crucial survival strategy that could shape the life cycle strategies of some opportunistic species.

IGP needs to be studied further as we move toward an integrated approach of marine ecology that recognizes the influence of both individual variability and community-level interactions. Other implications of IGP than those evoked in our study could be important. In the North Atlantic for example, the survival of C. finmarchicus early stages follows different seasonal patterns in areas where it is the dominant Calanus species than in areas where its larger congeners, C. glacialis and C. hyperboreus, co-occur. In such areas where all three species live together, the C. finmarchicus recruitment peak can occur several weeks after the spring bloom, much later than it usually does (Melle et al., 2014). In early spring, C. hyperboreus and C. glacialis, that are already active when C. finmarchicus initiates its reproduction fuelled by the phytoplankton bloom, could ingest C. finmarchicus eggs and affect its recruitment. Interestingly, IGP could contribute to the resistance of marginal Arctic marine ecosystems to the northward advance of the boreal C. finmarchicus under the pressure of climate change.

### AUTHOR CONTRIBUTIONS

KD designed and ran the laboratory experiments, participated to the in situ sampling, did the statistical analysis, ran the modeling experiment, and wrote most of the paper. FM contributed to the design of the laboratory experiment, contributed to the data analysis, designed, and executed most of the modeling experiment and contributed to the writing. SP contributed to the design of the laboratory experiment, provided access to the laboratory facility, and contributed to the writing. PJ contributed to the design and the execution of the laboratory experiment, and ran the in situ sampling. FC contributed to the design of the numerical experiment, provided data from Amundsen Gulf, and contributed to the writing.

### FUNDING

This work has been supported by an NSERC Discovery grant to FM, DFO, and Québec-Océan.

### ACKNOWLEDGMENTS

KD is grateful to PJ for his invaluable help in the field and in the laboratory, Geneviève Parent for her help with image analysis, Michel Starr and Liliane St-Amand for their help with laboratory analysis, David Levasseur for his assistance in the laboratory experiments and Daniel Bourgault for the Amundsen Gulf physical data. KD also thanks FM, SP, and Maurice Levasseur for the stimulating and constructive discussions and finally FC and Jeffrey Runge for their help on the manuscript. This is a contribution to the research programs of Québec-Océan, ArcticNet and UMI Takuvik.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmars. 2016.00185

The R code and forcing fields for the model can be found online at: https://github.com/NEOLab-Git/Chyp\_Egg\_1D.git

The detailed table of individual observations is provided.

### REFERENCES

**Conflict of Interest Statement:** 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.

Copyright © 2016 Dufour, Maps, Plourde, Joly and Cyr. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Modeling Plankton Mixotrophy: A Mechanistic Model Consistent with the Shuter-Type Biochemical Approach

Caroline Ghyoot 1 \*, Kevin J. Flynn2 , Aditee Mitra2 , Christiane Lancelot 1 and Nathalie Gypens 1

1 Ecologie des Systèmes Aquatiques, Université libre de Bruxelles, Brussels, Belgium, 2 Biosciences, Swansea University, Swansea, United Kingdom

Mixotrophy, i.e., the ability to combine phototrophy and phagotrophy in one organism, is now recognized to be widespread among photic-zone protists and to potentially modify the structure and functioning of planktonic ecosystems. However, few biogeochemical/ecological models explicitly include this mode of nutrition, owing to the large diversity of observed mixotrophic types, the few data allowing the parameterization of physiological processes, and the need to make the addition of mixotrophy into existing ecosystem models as simple as possible. We here propose and discuss a flexible model that depicts the main observed behaviors of mixotrophy in microplankton. A first model version describes constitutive mixotrophy (the organism photosynthesizes by use of its own chloroplasts). This model version offers two possible configurations, allowing the description of constitutive mixotrophs (CMs) that favor either phototrophy or heterotrophy. A second version describes non-constitutive mixotrophy (the organism performs phototrophy by use of chloroplasts acquired from its prey). The model variants were described so as to be consistent with a plankton conceptualization in which the biomass is divided into separate components on the basis of their biochemical function (Shuter-approach; Shuter, 1979). The two model variants of mixotrophy can easily be implemented in ecological models that adopt the Shuter-approach, such as the MIRO model (Lancelot et al., 2005), and address the challenges associated with modeling mixotrophy.

Keywords: constitutive mixotroph (CM), non-constitutive mixotroph (NCM), modeling, phytoplankton, zooplankton

## INTRODUCTION

Traditionally, planktonic protists are separated into two distinct functional groups: the phototrophic phytoplankton and the phago-heterotrophic microzooplankton. However, many protists assigned to these two groups are recognized as capable of combining phototrophy and phago-heterotrophy (Flynn et al., 2013). These mixotrophic protists have been reported in all planktonic functional groups (with the notable exception of the diatoms) and they include a large diversity of mixotrophic types (Flynn et al., 2013). In some mixotrophs, growth is mainly supported by phototrophy (e.g., Cryptomonas ovata; Tranvik et al., 1989) while in others, growth

#### Edited by:

Michael Arthur St. John, Institute for Aquatic Resources, Danish Technical University, Denmark

#### Reviewed by:

Peter Allan Thompson, Commonwealth Scientific Industrial Research Organisation, Australia Alberto Basset, University of Salento, Italy

> \*Correspondence: Caroline Ghyoot<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Marine Ecosystem Ecology, a section of the journal Frontiers in Ecology and Evolution

> Received: 21 June 2016 Accepted: 30 June 2017 Published: 18 July 2017

#### Citation:

Ghyoot C, Flynn KJ, Mitra A, Lancelot C and Gypens N (2017) Modeling Plankton Mixotrophy: A Mechanistic Model Consistent with the Shuter-Type Biochemical Approach. Front. Ecol. Evol. 5:78. doi: 10.3389/fevo.2017.00078 is mainly supported by phagotrophy (e.g., Poterioochromonas malhamensis; Sanders, 1991). Some mixotrophs are forced to use both trophic modes to grow and survive because some essential metabolites come specifically from one of the two metabolic pathways (e.g., the ciliate Laboea strobila; Stoecker et al., 1988); others are facultative mixotrophs (e.g., the dinoflagellate Fragilidium duplocampanaeforme; Park et al., 2015). Some mixotrophs use their second trophic mode to fulfill carbon requirements (e.g., the haptophyte Chrysochromulina brevifilum; Stoecker, 1998); for others, it is a means to fulfill nutrient requirements (e.g., the dinoglagellate Prorocentrum minimum; Stoecker, 1998); and for others yet, it is a mechanism to get specific metabolites (as phospholipids; Kimura and Ishida, 1989).

Among this large diversity of mixotrophic types, a major distinction can be made between mixotrophs depending on whether they photosynthesize using their own chloroplasts or using chloroplasts acquired from their phototrophic prey (Mitra et al., 2016). Mixotrophs that actively synthesize and maintain their own chloroplasts are called constitutive mixotrophs (CMs). Mixotrophs that do not constitutively synthesize chloroplasts but photosynthesize using chloroplasts acquired from their phototrophic prey are called non-constitutive mixotrophs (NCMs). The acquired chloroplasts in NCMs remain functional for periods ranging from hours to days, depending on the type of NCM (Mitra et al., 2016). Generalist NCMs, i.e., those that acquire chloroplasts from a broad range of phototrophic prey, have a poor ability to maintain functional chloroplasts (Dolan and Pérez, 2000). By contrast, specialist NCMs, i.e., those that need to acquire chloroplasts from specific phototrophic prey, can maintain a photosynthetic activity for long periods (Stoecker et al., 2009; Hansen et al., 2013). Specialist NCMs can further be divided into those that retain the entire photosynthetic prey (endosymbiosis) and those that only retain the chloroplasts (kleptochloroplasty; Mitra et al., 2016). In a functional capacity, the CMs are closer to the strict phototrophs while the NCMs are closer to the strict phago-heterotrophs.

While mixotrophy has been reported for a long time, it is now appreciated as being much more widespread in aquatic ecosystems than initially thought. Previously, the mixotrophic status was only accorded to some dinoflagellates, forams, radiolarian, and acantheria while now, mixotrophy has been acknowledged in all eukaryote planktonic microorganism groups, except diatoms (Sanders and Porter, 1988; Burkholder et al., 2008; Flynn et al., 2013). In addition, these species are not limited to a specific habitat: mixotrophs are observed in both freshwater and marine (Sanders, 1991; Stoecker et al., 2009), oligotrophic and eutrophic systems (Burkholder et al., 2008), and from polar to equatorial regions (Zubkov and Tarran, 2008; Stoecker et al., 2009; Sanders and Gast, 2012). In coastal areas, mixotrophic flagellates can account for some 50% of the pigmented biomass (Havskum and Riemann, 1996) and chloroplast-containing ciliates can account for 40–>60% of the planktonic ciliates in summer (Stoecker et al., 1987; Bernard and Rassoulzadegan, 1994). In oceanic waters, mixotrophs account for 40–95% of the bacterivory in the euphotic layer (Zubkov and Tarran, 2008).

Constitutive mixotrophy is suggested to play an important ecological role when inorganic nutrients are low (e.g., in oligotrophic systems; e.g., Arenovski et al., 1995) or unbalanced (e.g., in eutrophied systems; e.g., Nygaard and Tobiesen, 1993; Burkholder et al., 2008), and when light is limiting (e.g., during the polar night or occasionally in eutrophied systems; e.g., Bird and Kalff, 1986; Roberts and Laybourn-Parry, 1999; Jones et al., 2009). In such systems, prey ingestion provides nutrients and energy. Similarly, non-constitutive mixotrophy may be advantageous in "high light–low prey" ecosystems (e.g., in oligotrophic systems; e.g., Skovgaard, 1998; Hansen et al., 2013) because acquired phototrophy supplies carbon by photosynthesis. Finally, mixotrophy is worthy of interest in an environmental perspective because most marine harmful algae have been reported as mixotrophs (Stoecker et al., 2006; Burkholder et al., 2008; Jeong et al., 2010).

Several models have been developed with the specific aims of exploring planktonic mixotrophy from an ecophysiological or ecological (species competition) perspective. Most interest has been leveled at the potential impact of mixotrophs on the microbial food web structure and functioning, and the conditions under which mixotrophs may likely coexist with strict phototrophs and heterotrophs (Thingstad et al., 1996; Baretta-Bekker et al., 1998; Stickney et al., 2000; Jost et al., 2004; Hammer and Pitchford, 2005; Hood et al., 2006; Flynn and Mitra, 2009; Crane and Grover, 2010; Ward et al., 2011; Våge et al., 2013). By far the greater effort has been applied to CM organisms. The complexity of the model structure varies widely among such studies. The simplest models describe mixotrophy as the ability to combine both phototrophy and heterotrophy without any feedbacks or trade-offs between the two nutritional modes and organisms have a fixed stoichiometry (e.g., Hammer and Pitchford, 2005). At the other extreme, the most complex model explicitly describes the main regulative processes that occur between phototrophy and heterotrophy in the mixotroph and allows for a variation of the cellular stoichiometry (Flynn and Mitra, 2009).

Despite the potential significance of mixotrophy in aquatic systems and the existence of mathematical models describing this trophic mode, few ecological/biogeochemical models consider it explicitly (e.g., Mitra et al., 2014). Those models that do include mixotrophs demonstrate the potential for mixotrophy to significantly impact the flow of energy and nutrients in the system. In particular, primary production can potentially be increased by a factor 2 compared to a system in which mixotrophy is not considered, thanks to a shortened and more effective chain from nutrient regeneration to primary production (Mitra et al., 2014).

Considering mixotrophy in biogeochemical/ecological models is however challenging due to the large diversity of mixotrophic types, the scarcity of data allowing the parameterization of physiological processes, and the need to make the inclusion of mixotrophy into existing ecosystem models as simple as possible. Indeed, the addition of a mixotroph functional type description into ecosystem models will inevitably be tempered by the structure of the extant ecosystem model.

Regarding the diversity of mixotrophic types, this paper proposes a flexible mechanistic model featuring the two main types of mixotrophy: constitutive and non-constitutive. The two mixotrophic types have indeed different ecophysiology and different relation to their prey (the NCM being strictly dependent on the presence of prey); therefore, they affect differently the ecosystem dynamics (Mitra et al., 2016). In addition, the version describing constitutive mixotrophy can be configured to represent either CMs that favor either phototrophy or heterotrophy. Simulations have been performed to compare each mixotrophic type with the strict trophic type it most closely resembles from a functional perspective (the CM with the strict phototroph and the NCM with the strict heterotroph) in order to explore the competitive advantage of mixotrophy over strict forms.

As for the explicit inclusion of mixotrophy into existing biogeochemical/ecological models, we were mainly interested in introducing mixotrophy in a form consistent with the mechanistic structure of the AQUAPHY model (Lancelot et al., 1991); this model has been used to describe primary production via phytoplankton growth physiology within several biogeochemical models (Billen et al., 1994; Lancelot et al., 2000, 2005). A feature of AQUAPHY is that the total cellular carbon is divided into separate classes of components on the basis of their function: these comprise synthetic and structural material (i.e., photosynthetic apparatus, ribosomes, genetic material, membranes, etc.), stored carbon (i.e., carbohydrates, lipids), and carbon monomers (i.e., photosynthetic products and precursors of complex molecules). This mechanistic approach was suggested by Shuter (1979) as a means to successfully describe the observed behavior of a variety of unicellular algal species, qualitatively and quantitatively. The work is thus a first step in the implementation of mixotrophy in the biogeochemical models that adopt a "Shuter mechanistic approach" to describe the phytoplankton physiology.

### METHODS

### Constitutive Mixotrophy

The model of constitutive mixotrophy was constructed on the basis of the following hypotheses, summarizing the main qualitative observations related to CMs : (i) the CM is preferentially phototrophic but it can supplement its nutrient requirements (under inorganic nutrient limitation) or its carbon requirements (under light limitation) by ingesting prey (e.g., in Jones et al., 1995; Carvalho and Granéli, 2010; McKie-Krisberg et al., 2015); (ii) the maximum mixotrophic growth is not allowed to exceed the maximum phototrophic growth; (iii) the phototrophic growth has to account for at least 10% of the mixotrophic growth, such that we describe an obligate phototroph (Caron et al., 1993; Brutemark and Granéli, 2011); (iv) the nutrients regenerated by heterotrophic activity (PO3− 4 , NH+ 4 ) are re-assimilated to sustaining phototrophic growth with any surplus being excreted (Flynn and Mitra, 2009).

The model combines an adaptation of the AQUAPHY model (Lancelot et al., 1991; Ghyoot et al., 2015) for the phototrophic pathway, coupled to a simple zooplankton model (Lancelot et al., 2005) for the phago-heterotrophic path. It is important to note that, in keeping with the formulation of the AQUAPHY model, the description of mixotrophy given here describes biomass growth (i.e., molC m−3 time−1 ) rather than per capita growth (i.e., C C−1 time−1 ) as in some other models (e.g., Flynn and Mitra, 2009). **Figure 1A** shows the schematic representation of the CM growth physiology, linking the phototrophic (dark gray) and heterotrophic (light gray) components. The CM model contains six state variables (**Table 1**) describing intracellular components on the basis of their function: the functional and structural metabolites (e.g., nucleic acids, proteins) synthesized by phototrophic and heterotrophic activities (Fphot and Fhet , respectively), soluble carbon monomers (i.e., early products of photosynthesis; SC), carbon reserves (i.e., carbohydrates, fatty acids; RC), intracellular soluble phosphate (SP), and intracellular soluble inorganic nitrogen (SN). The total C-cell biomass (mmolC m−3 ) is given by the sum of Fphot, Fhet, SC, and RC. The Fphot and Fhet pools are assumed to have a fixed C:N:P stoichiometry, based on biochemical constraints (Geider and Laroche, 2002). Variable cellular stoichiometry is enabled by considering the additional C, P, and N accumulated as carbon monomers (SC), carbohydrates and fatty acids (RC), soluble inorganic phosphorus (SP), and soluble inorganic nitrogen (SN). The model of constitutive mixotrophy is linked to three state variables describing external inorganic nutrients—dissolved inorganic nitrogen (DIN = NO− 3 + NH+ 4 ) and phosphate (PO3− 4 )—and also the prey which in reality expresses a variable C:N:P stoichiometry (though here prey stoichiometry is fixed as the emphasis is on the description of the autecology of the mixotrophs, rather than system ecology).

**Tables 2**, **3** show, respectively, the nine conservation equations related to the state variables and the associated processes. Phototrophic growth (µ phot; Equation 17) is controlled by the concentration of the structural and functional metabolites related to phototrophic activity (Fphot), by the limitation in carbon monomers SC—either directly produced by photosynthesis (ϕ; Equation 20) or indirectly by RC catabolism (catRC; Equation 10)—and by the limitation in internal soluble inorganic nutrients (SN and SP). The SC limitation is formulated by a Michaelis–Menten equation in which the substrate concentration is expressed by XSC − kSC , with XSC = SC F phot and kSC , which is the minimum value for XSC (assumed to be equal to the half-saturation constant for SC assimilation). The inorganic nutrient limitation is formulated by the Liebig's minimum law in which the limitation for each nutrient is expressed by a hyperbolic function depending on XSN,P (the ratio between SN,P and the N,P contained in Fphot). The uptake of inorganic nutrients (uptDIN and uptPO4; Equations 26 and 27) depends on the external nutrient concentration and the status of the internal nutrient reserve. The phototrophic respiration (respphot; Equation 22) includes costs for cellular maintenance and for synthesis of new Fphot .

Prey ingestion (graz; Equation 14) is controlled by the concentration of the structural and functional metabolites related to heterotrophic activity (Fhet) and by prey availability; the latter

FIGURE 1 | Schematic representation of functioning of constitutive mixotroph (CM; A) and non-constitutive mixotroph (NCM; B). Total mixotroph biomass is divided into six components: structural and functional metabolites related to phototrophic activity (Fphot) and heterotrophic activity (Fhet), carbon monomers (SC), carbon reserves (RC), soluble inorganic nitrogen (SN), and soluble inorganic phosphorus (SP). The "dark gray part/solid lines" is the phototrophic part and the "light gray part/dashed lines" is the heterotrophic part.

is controlled by a sigmoid (type III) function. From the ingested prey, a non-assimilated fraction is egested (egest; Equation 11) as dissolved and particulate organic matter (loss terms), a fraction is respired to meet the heterotrophic metabolic costs and is released as carbon dioxyde (resphet; Equation 21), and the last fraction is assimilated and contributes directly to the heterotrophic growth (µ het; Equation 16). If the nutrient content of the prey is higher than the nutrient required by Fhet, the surplus is regenerated as NH+ 4 and PO3− 4 (regi with i = N, P; Equation 23). These regenerated nutrients can be retained up to a maximum level in the inorganic nutrient reserves SN and SP (reti ; Equation 24), contributing to the phototrophic growth. If the reserve capacity is full, the surplus is excreted to the environment (excri ; Equation 12).

TABLE 2 | Conservation equations for the constitutive mixotroph (CM).

Mixotrophic growth is computed as the sum of the phototrophic and the heterotrophic growth rate, but limited by the maximum mixotrophic growth. The latter is computed differently according to whether phagotrophy is used to offset a lack of carbon or a lack of nutrients. When the 24 h-average SC limitation is below a threshold value (threshSClim) set here as 0.15, we assumed that the mixotroph switches to perform increasing levels of phagotrophy to acquire carbon (in addition to nutrients if also nutrients limited) and, therefore, the maximum mixotrophic growth is equal to the phototrophic growth obtained with no limitation (µ mix max; Equation 18). When the 24 h-average SC limitation is above this threshold value, we assumed that the CM undertakes phagotrophy to acquire nutrients and, therefore, the maximum mixotrophic growth is equal to the phototrophic growth obtained with no nutrient limitation (µ mix max; Equation 19). If the sum of the phototrophic and the heterotrophic growth is higher than the maximum mixotrophic growth, either the grazing or the photosynthesis rate is regulated to TABLE 3 | Equations that describe the processes occurring in the constitutive mixotroph (CM).