Reproducing the virus‐to‐copepod link in Arctic mesocosms using host fitness optimization

By shunting material out of the predatory pathway toward detritus and dissolved material, viruses are believed to have an important impact on biogeochemical functions of the pelagic microbial food web. To include viruses as a single plankton functional type (PFT) in dynamic food web models is, however, not trivial since they will then compete with predators for the same host/prey community as a shared limiting resource. As recently shown, one can solve this problem by introducing adaptation in the defensive and competitive traits of the host (prey) community. We here show how this can reproduce central aspects of viral dynamics as observed in a set of Arctic mesocosm experiments. In these experiments, contrasting microbial trophodynamics have previously been linked to the trophic cascades generated by seasonal vertical migration of large Arctic copepods. This approach thus produces a quantitative theory for the mechanisms regulating virus‐to‐prokaryote and lysis‐to‐predation ratios, and integrates this with a central role of predator top‐down control in pelagic microbial food webs.

Lytic viruses are believed to play an important role in flux partitioning, biodiversity, and evolution in microbial ecosystems (Suttle 2007). However, even at the coarse level of variations in virus abundances, our understanding of how the set of underlying control mechanisms work in concert is limited. Such underlying mechanisms presumably span from the molecular mechanisms behind viral infection and host defense, via trophic interactions in the microbial food web, to the environmental drivers that act on this food web. A description integrating these levels would seem to involve a rather complex set of interactions and feedbacks, perhaps explaining why a generally accepted theory for viral abundance and activity in natural systems is still lacking, three decades after the general acknowledgement of high viral abundances in natural aquatic ecosystems (Bergh et al. 1989;Proctor and Fuhrman 1990). For the study of microbial food web dynamics in mesocosm experiments, a "minimum" food web model without viruses has proven useful (Thingstad et al. 2007;Larsen et al. 2015;Tsagaraki et al. 2018). It is in principle possible to combine this minimum model with a dynamic description of viruses that resolves the prokaryote community to species or even to strain level (Thingstad et al. 2014;Våge et al. 2016), but gives a substantial increase in the number of plankton functional types (PFTs) and model complexity, reducing the conceptual transparency.
Host-virus models have been constructed to simulate laboratory systems (Levin et al. 1977;Blackburn et al. 1996;Bohannan and Lenski 1999;Weitz et al. 2005) and for organisms of particular interest in mesocosm experiments (Ruardij et al. 2005). For viruses in natural environments, neural network models may have a considerable predictive potential (Winter et al. 2005), but gaining mechanistic insight from this type of models can be difficult. With a more theoretical focus, models of the dynamics (Weitz et al. 2005) and the steady state (Härter et al. 2014;Thingstad et al. 2014) of host-virus arms-races have been developed. It is possible to construct multi-trophic models with viruses that retain the simplicity of black-box descriptions (no internal community resolution). This, however, creates a classical coexistence problem (Hutchinson 1961) as predators and viruses then compete for the same prey/host community as their shared limiting resource. While resolved models have the possibility to reach steady state through within-community adjustments in the richness and evenness of the host and virus populations, black-box formulations need other feedback mechanisms (van Velzen 2000;Thingstad and Våge 2019) to allow such coexistence. In the multi-trophic model of Weitz et al. (2015), this feedback is obtained by closing the upper end of the food web with a top predator loss rate that is positively correlated with the community size of this top predator. An independent copepod abundance, decoupled from their immediate microbial food supply is an essential part of our interpretation of observed mesocosms dynamics. Models containing constraints on copepod abundances in order to allow for viruses were therefore not considered desirable.
As an intermediate alternative to either resolved or blackbox approaches, one can focus on the community-level consequences of internal shifts in community composition. Rather than attempting any explicit representation of all the phenotypic and/or genotypic processes that drive shifts in community composition, we assume that the net outcome of selection of species and strains with higher fitness is also a community with higher fitness. Technically, this is summarized in a dynamic community strategy index (S) that regulates the balance between the community's competitive and defensive traits. By assuming changes in S to be proportional to the local fitness gradient, this gives a host community that dynamically drifts toward fitness optimum. As this optimum depends on growth conditions and virus abundance, a feedback is created where high viral abundance favors the selection of defensive hosts, reducing viral production. As shown by Thingstad and Våge (2019), this can provide sufficient feedback for stable coexistence of viruses and predators.
We here explore the extent to which this type of description is able also to reproduce observed host-virus dynamics. We have combined the results from two mesocosm experiments with similarities in experimental design, but large differences in the observed prokaryote and viral responses (Larsen et al. 2015;Sandaa et al. 2017).

Summary of previous interpretations of observed ecosystem responses
The two experiments denoted PAME-I and PAME-II were performed in Kongsfjorden at the Ny Ålesund research station, Svalbard, Atlantic Arctic (78 55 0 30 00 N 11 55 0 20 00 E) in late (PAME-I: July-August, 2007) and early (PAME-II: June-July, 2008) summer, respectively. Experimental design and results are discussed in detail by Larsen et al. (2015). Prokaryote-virus dynamics is treated in detail by (Sandaa et al. 2017). Aspects important in the present context are summarized below for ease of reference.
The two experiments were done before (PAME-II) and after (PAME-I) the period where the dominating copepod Calanus finmarchicus leaves surface water for deep-water winter hibernation (Falk-Petersen et al. 2009). Presumably because of this timing, PAME-II had an initial copepod population ca. 3× higher than PAME-I. Inferences on the top-down effects from copepods are based on comparison of the two experiments. All mesocosms received the same daily dose of mineral nutrients (orthophosphate and inorganic nitrogen in Redfield-ratio; 100 nM PO 4 -P, N : P = 16 : 1). In both experiments, mesocosms were arranged in two gradients with increasing daily additions of easily degradable dissolved organic carbon (DOC) as glucose, from 0× to 3× the Redfield ratio (C : P = 106 : 1, molar) in glucose-C relative to added orthophosphate-P. Both experiments thus had a design allowing detection of potential effects of degradable DOC. In PAME-I, one of the two gradients was kept silicate replete to allow the growth of diatoms. No silicate was added to the other. Differences between the two gradients in PAME-I are used to discuss differences in food web responses with and without diatoms. In PAME-II, the two glucose gradients both received silicate, but differed in the use of ammonium versus nitrate as N-source. No significant effects of this difference in nitrogen source were detected (Larsen et al. 2015) and differences between nitrate and ammonium as N-source are not discussed further here.
Observed food web responses have been interpreted using the "minimum" food web model in Fig. 1A (Thingstad 2020). Important in the present context is the model's combination of three linear food chains from mineral nutrients to copepods: (1) through bacteria via heterotrophic flagellates and ciliates (the microbial loop), modified by potential bacterial limitation by labile DOC, (2) through autotrophic flagellates and ciliates, and (3) through diatoms ("classical" food chain), modified by potential silicate limitation of diatoms. This model has an odd number of links in food chains 1 and 3, but an even number of links in food chain 2. Cascading effects from copepods to primary nutrient consumers therefore create two very different food webs at low (Fig. 1B) and high ( Fig. 1C) copepod abundances, subsequently termed the model's LZ and HZ state, respectively. The flagellatedominated HZ state (Fig. 1C) corresponds to the main observed characteristics of the PAME-II experiment (Larsen et al. 2015) with no significant observed effects of glucose addition. The diatom-bacteria balance in the LZ (Fig. 1B) situation is slightly more complicated. In an LZ situation with silicate present (Fig. 1B), most added mineral nutrients will in this model be immobilized in a slowly grazed diatom population. Silicate addition will then have an inhibiting effect on bacterial activity, mediated through strong mineral nutrient competition from the large diatom population. This was the response observed in an earlier mesocosm experiment (Isefjorden, Denmark; Thingstad et al. 2007) for which the model was originally developed. In PAME-I, the observed response was "opposite" with glucose addition leading to diatom disappearance. The suggested explanation (Larsen et al. 2015) is the species composition of the diatom community, in PAME-I totally dominated by a small (ca 8 μm diam.) Thalassiosira sp. Assuming small diatoms to be grazed by ciliates (Verity et al. 1988;Montagnes 1996;Hansen et al. 1997;Nejstgaard et al. 1997), the model response shifts to the "glucose-replete bacteria outcompete diatoms" situation observed in PAME-I.
In summary, the model predicts size of the bacterial community to respond to glucose and silicate manipulations in the LZ situation. In the HZ situation, where flagellates dominate the response, while both diatoms and bacteria are sandwiched between a high-predatory pressure and high-mineral nutrient competition from the large community of autotrophic flagellates, neither glucose nor silicate has much effect. A similar difference in responsivity was observed also in viral abundance and community composition between the PAME-II (LZ) and PAME-I (HZ) experiments (Sandaa et al. 2017). Such similarity in prokaryote and virus responsivity is what one would expect if changes in virus community size and composition follows as a response to food web induced changes in their prokaryote host community.

Model considerations
Except for the viruses, the response pattern described above has been reproduced in simulations using the "minimum" model (Larsen et al. 2015). With a model able to capture central aspects of microbial food web trophodynamics, we wanted here to explore whether the observed variations in viral abundances (Sandaa et al. 2017) could also be reproduced using the technique with an adaptive bacterial community.
Assuming bacteriophages to numerically dominate the virus community, only bacteria-virus interactions are added to the minimum model (Fig. 1D). Virus and bacterial community sizes are represented by single state variables (V and B, respectively. Units in nM-P, converted to abundance by assuming fixed phosphorus per virus and phosphorous per bacterium conversion factors. Symbols are summarized in Table 1). The full set of differential equations is given in Box 1. Key equations for the virus-bacteria interactions are explained in the following. Changes in V are governed by the differential equation: where α V BV is the infection = lysis rate (in nM-P h −1 ). Viral community properties are thus described by three parameters: an effective community adsorption constant α V (L nmol-P −1 h −1 ), a specific decay rate δ V (h −1 ) and a yield Y V representing the fraction of host phosphorous incorporated into viruses. The remaining fraction (1 − Y V ) is assumed to be instantaneously remineralized to orthophosphate upon cell lysis. Description of the host community's interaction with viruses follows the formalism developed by Thingstad and Våge (2019). The effective community adsorption constant α V and an effective community bacterial nutrient affinity α B are The "minimum" model previously used with three food chains from dissolved inorganic phosphorous to copepods via bacteria, autotrophic flagellates, and diatoms. As these three linear food chains have 3, 2, and 1 step(s), respectively, the cascading effects from copepods are different, with (B) bacteria or diatoms stimulated in the low-copepod (LZ) situation, while (C) flagellates dominate in the highcopepod (HZ) situation. In the work discussed here (D), bacterial viruses are added and diatoms are assumed to be grazed by ciliates (dotted arrow). A. flagellates and H.flagellates are auto-and hetero-trophic flagellates, respectively. LDOC, DIP, and Si are labile dissolved organic carbon, dissolved inorganic phosphorous, and silicate, respectively.
With this formulation, S = 0 corresponds to the undefended situation α V = α m V À Á with maximum community nutrient competition ability α B = α m B À Á , while S=1 corresponds to total community immunity (α V = 0) at the price of no nutrient uptake (α B = 0). For trade-off τ < 1, the loss in relative com- for an increase in S will be less than the accompanying reduction in vulnerability αV α m V .
Fitness of the bacterial community is defined as the community's net specific growth rate: where μ B is bacterial specific (per unit B) growth rate, and g H is their specific (per unit B) grazing loss to heterotrophic flagellates. As we here focus on host adaptation to viral attack, we ignore grazing defense mechanisms, and g H is therefore assumed to be independent of S.
The rate of change in S is assumed to be proportional to the local fitness gradient ( ∂F ∂S ) (Abrams 2001): where r is a proportionality constant summarizing the community level effects of the evolutionary and phylogenetic processes changing internal community structure. Since bacterial growth rate μ B is a function of nutrient concentrations (labile DOC [L] and orthophosphate [P]), the rate of change in S will depend on growth rate limitation as well as on the abundance of viruses.
The ambition has been to add this description to the minimum model, altering other aspects of the model as little as possible. New and previously used (Larsen et al. 2015) parameter values are listed in Table 1.
To be able to reproduce the observed silicate effect on bacteria in PAME-I, the particular feature with small C-rich diatoms subject to ciliate predation is also added to the minimum model (Fig. 1D, see Supporting Information for more model details). While this does not affect the bacteria-virus interactions directly, it is crucial to explain the silicate effect on bacterial host dynamics in the PAME-I experiment (Larsen et al. 2015) and therefore required to simulate host-virus dynamics.

Specific growth rates for bacteria (B), autotrophic flagellates (A), and diatoms (D)
Specific ingestion rates for heterotrophic flagellates (H), ciliates (C), and copepods (Z) Differential equations: Where fitness F is defined as (using Eq. 7): and: Using net bacterial growth rate as the state-dependent fitness function (F), a generalized form of it is: where μ B is bacterial specific (per unit B) growth rate, and g H is their specific grazing loss to heterotrophic flagellates. As we here focus on host adaptation to viral attack, we ignored grazing defense mechanisms, and g H is therefore assumed to be independent of S.
To prevent predation from reducing bacterial abundance too much, and thereby creating unrealistically high virus-tobacteria ratio (VBR) values, a lower prey limit (B lim , Table 1) for predation on bacteria is added to the minimum model.
Simulations are otherwise based on the assumption that the microbial part of the food web was in an approximate steady state when the mesocosms were filled (see Supporting Information for details). This steady state depends on two drivers: total phosphate (P T ) available to the system, and the amount of copepods (Z). The calculated initial state is thus different for the LZ and HZ runs.
The Matlab ® script used to run the model is included in the Supporting Information.

Results
A main difference observed between the PAME-I and PAME-II experiments was the responsivity in prokaryote abundance to glucose and silicate additions (Fig. 2). This difference is reproduced by the model (Fig. 2) which has a much more dynamic pattern in the LZ, compared to the HZ runs. In the LZ runs, the difference in timing of treatment effects are also reproduced with the glucose effect dominating around days 2-5, while the main silicate effect is delayed until days 9-12; consistent with observations from PAME-I. The model's rapid decrease in bacterial abundance in LZ around day 5 is primarily predator-driven, as indicated by the rapid concomitant increase in heterotrophic flagellates (Fig. 3). This supports the previous suggestion (Larsen et al. 2015) that the observed oscillatory pattern in total prokaryote abundance in PAME-I reflects a Lotka-Volterra type prey-predator dynamics between the prokaryote community and heterotrophic flagellates, largely absent from both model (HZ runs Fig. 2) and observations (PAME-II Fig. 2) when copepod stock is high. The late effect of silicate on bacteria is, in the model, an indirect effect from the Si-stimulated small diatoms, mediated both by the trophic cascade propagating via diatom-consuming ciliates, the diatom production of organic substrates (L) for bacteria, and bacteria-phytoplankton competition for mineral nutrients.
Simulated viral abundances also follow the main observed PAME-I vs PAME-II contrast (Fig 2) with less viruses and little response to glucose and silicate in the HZ simulations (2Bb).
The model largely reproduces the observed positive effect of glucose on viral abundance. The indirect silicate effect on bacteria is translated to viruses toward the end of the simulations, but not to the same extent as in the observations (Fig. 2).
With most modeled values within a factor of 2 from the observed, not only the effect patterns, but also the level of viral and prokaryote abundances are reasonably reproduced.
Although the modeled VBR does not contain fundamentally new information, it is a sensitive indicator to the timing and level of changes in bacterial abundance and the subsequent response in viruses. With virus production following the increases in host abundance, one would expect that rapid increases in the host abundance are reflected in decreasing VBR (rapidly increasing denominator). This effect can be seen in the LZ runs, both in the early phase with rapidly increasing bacterial abundance in all treatments, and toward the end, where the rapid increase in bacterial abundance in silicateamended treatments is reflected in reduced model VBR values. The intermediate peak in LZ VBR values around days 5-10 is reasonably relative to the observed virus-to-prokaryote (VPR) values (Fig. 2). The model also reproduces the low VPR values observed in the PAME-II experiment (Fig. 2).
For illustration of how the zooplankton effect cascades through the food web, the modeled food-web dynamics is presented in Fig. 3. For corresponding observational food web Fig. 3. Simulated dynamics of the microbial food web. HZ (red) and LZ (blue) situations with (thick lines) and without (thin lines) glucose and with (broken lines) and without (solid lines) silicate added. For of observed food web results in these experiments, see Fig. S2 and Larsen et al. (2015) data and discussion of these, see Supporting Information Fig. S2 and Larsen et al. (2015).

Discussion
Recent theoretical work has demonstrated how coexistence on a single resource is possible, both between viruses and predators Våge 2019) andbetween predators in general (van Velzen 2000), in models where host/prey traits are made adaptive. We have here included such a mechanism in an established food-web model and shown how this can reproduce observed patterns in virus abundance and dynamics in two contrasting mesocosm experiments.
The technique allows "gray-box" host-virus descriptions where the internal structure of the host community is represented by a community level strategy S. S regulates the balance between competitive and defensive properties of the host community. With this, one can avoid explicit representations of poorly understood mechanisms such as how competitive and defensive traits are distributed among pelagic prokaryotes and how this relates to the selection of mutants and the arms races that presumably shape the internal structure of microbial communities (Thingstad et al. 2014). The price for this simplification is to loose the more explicit connection between host group diversity and nutrient cycles present in models which resolve the host community into subgroups. The potential for model comparison with the detailed diversity data obtainable with modern sequencing techniques is also limited.
The balance between competition and defense at community level seems crucial for how the prokaryote community behaves in a food-web context. Constructing the relationship between S, τ, competition, and defense (as in Fig. S1) may, however, be difficult; partly because of the somewhat abstract nature of the theoretical concept of a "community strategy." In the formulation used here, S will always move in the direction of increasing host community fitness. This optimization is clearly related to the concept of an evolutionarily stable community (ESC) (Edwards et al. 2018), but establishing an experimental fundament for ESC may also seem difficult.
In our model, lysis, and thus virus production, is directly coupled to food-web dynamics through host abundance (Eq. 13). To get variations in viral numbers correct, it is therefore crucial to get the variations in bacterial abundance correct. The minimum model, amended with the assumption of a small, DOC-producing diatom grazed by ciliates (Thingstad et al. 2008;Larsen et al. 2015), captures relatively well how the copepod stock modulates the effects of glucose and silicate additions (Figs. 2, 3) on bacterial abundance. The variations in virus abundance can essentially be understood as the result of high production in periods where host abundance is high, combined with a relatively slow viral loss rate (determined by the decay constant δ V ) when host abundance is low. With a reasonably correct pattern for host abundances, the model also reasonably reproduces the observed effects of copepods, glucose, and silicate on viral abundances (Fig. 2) With the set of parameters used here (Table 1), the values of S are within the narrow range 0.95 < S < 1 (Fig. 4), indicating a community dominated by defensive hosts. The range of S is sensitive to some of the other parameter choices. In particular, a reduction in the maximum viral adsorption constant α m V will reduce the value of S (less defensive community). The level and dynamics of S is different for the HZ and LZ runs (Fig. 4). So were also the fingerprint data obtained for prokaryote (DGGE) and virus (PFGE) community composition in the PAME-I and PAME-II experiments (Sandaa et al. 2017). Drawing a more direct and quantitative connection between the community strategy S and community composition data seems, however, difficult at this stage. Viral dynamics is sensitive to the host adaptation rate r and the decay rate δ V . We have little experimental evidence to constrain r. The value used here (0.001 h −1 ) corresponds to a characteristic time scale of nearly a month. This may be adequate for evolutionary processes such as the selection of new mutants and the dynamics of host-virus arms races. It seems too long, however, for shifts in dominance of existing species which, even in our Arctic experiments, can occur over time scales of days (Tsagaraki et al. 2018). Viruses produced will remain longer in the system if δ V is small. The relatively stable viral abundances through the period with low host abundances, and therefore low virus production (day 6 to day 9, LZ runs), is thus related to the choice of a relatively low decay rate (4.8% d −1 ).
The model uses a rather low temperature correction for the α parameters (Q 10 = 1.4), somewhat higher for I max (Q 10 = 1.9) (see discussion in Thingstad and Aksnes 2019). If this model is representative for the Arctic, environmental changes affecting external drivers such as water column stability (and therefore total available nutrients) or copepod abundance in surface waters likely have larger effects on structure and function of the microbial food web and its viruses than direct microbial responses to changes in water temperature.
While the parameter set of the minimum model without viruses has been found adequate for many experiments (Thingstad et al. 2007;Larsen et al. 2015;Tsagaraki et al. 2018), this is a first and exploratory attempt to include viruses. The new, virus-related, parameters should therefore not be considered well constrained.
Prokaryote defense mechanisms are not restricted to viral attack, but include also grazing resistance (Pernthaler 2005), shown to be enhanced in situations with a combination of high grazing pressure and nonlimiting concentrations of degradable organic-C (Matz and Jürgens 2003). Selection for grazing resistant forms in prokaryote community composition has been suggested to explain an observed strong dampening of the cascade from copepods at the level of prokaryotes, supported by a model where a predator-defensive bacterial community was added to the minimum model (Tsagaraki et al. 2018). It seems technically feasible to extend our model with an extra bacterial strategy index representing predator defense in the prokaryote community, but this has not been explored here.
In terms of "organism" size, the connection drawn here between viruses and copepods covers about half the pelagic food chain. It is thus a contribution to the discussion of viruses-towhales ecosystem modeling (deYoung et al. 2004). "Wasp-waist" control is defined as the situation where a key group of organisms act as a resource control for the food chain above and a predator control for the part below (Fauchald et al. 2011). The important role of Arctic copepods as a food resource for higher trophic levels is well recognized (Conover et al. 1995). In the model used here, the copepod-to-virus link is mediated through the cascading effects illustrated in Fig. 1 and how these control the bacterial host community. With this combined role in topdown and bottom-up control, copepods seem to have a key "wasp-waist" role in Arctic food chains.
As viral lysis shunts energy out of the predatory food chain(s) toward dissolved material and detritus, it affects the f-ratio (the fraction of total primary production fuelled by imported nutrients), illustrating the importance of incorporating viruses into large-scale ecological models. The technique used here suggests how this can be done at the relatively low cost of adding two new dynamic variables (V and S). As large data sets exist for viral abundance and virus to prokaryote ratios (Wigington et al. 2016;Parikka et al. 2017), this may open for a fitting to observational data, potentially providing some constraint on the "community trade-off" parameter τ.