香京julia种子在线播放

Both our N-limited and our P-limited models include heterotrophic bacteria (H) and their viruses (V_H), cyanobacteria (C) and their viruses (V_C), zooplankton (Z), dissolved nutrients in inorganic form (N_in or P_in), and dissolved nutrients in organic form (N_org or P_org). This simplified representation of the microbial loop is similar to the one introduced by Weitz et al. (2015) to study the N-limited case only, except that we omit eukaryotic autotrophs, which is justified by the typical dominance of prokaryotic phytoplankton in oligotrophic environments (Follows et al., 2007). See Appendix A for model equations.

The P-limited version of our model differs structurally from the N-limited version only with regard to which form of the nutrient can cells take up. Here, we consider nutrient limitation in terms of the Liebig's law of the minimum (Baar, 1994), which is typically implemented by calculating the growth rates associated with each nutrient independently, and subsequently choosing the smallest to set the growth of the population. This differs from experimental approaches, which typically identify a nutrient as limiting when it falls below its detection limit (e.g., 0.05 μmol/L Graziano et al., 1996; Maat and Brussaard, 2016 for nitrate and 0.01 μmol/L for phosphate Maat and Brussaard, 2016). Both approaches are connected, as the expectation is that (i) a population will draw down the most limiting nutrient as much as it possibly can to utilize it for growth, and (ii) the nutrient that shows the lowest concentration will produce the lowest growth rate.

In both versions of our model, bacteria take up and assimilate nutrients following a Monod (i.e., saturating) function with parameters r (maximum growth rate), K (half-saturation constant) and ε (efficiency) (Monod, 1950). However, in the N-limited version, heterotrophic bacteria only use organic N (Equations A1 and A7), whereas cyanobacteria can only use inorganic N (Equations A2 and A6). In the P-limited system, both heterotrophic bacteria and cyanobacteria take up only inorganic P for growth and secrete phosphatase in the water to degrade organic phosphorus into its inorganic form (Equations A8, A9, and A14). In both versions, a fraction of the nutrient assimilated by bacteria is lost to respiration at a rate λ. An additional basal exudation rate (σ) represents the loss of organic material.

Zooplankton grow by consuming bacteria following a linear functional response with predation rate γ (Equations A3 and A10). Only a fraction of nutrients is used for growth, encoded with the efficiency parameter p_g. A fraction p_ex is released as biomass and the rest is lost to the inorganic nutrient pool through respiration. In addition, zooplankton have a constant basal respiration rate, λ_Z.

Viruses produce a fixed amount of progeny (burst size, β), virions that are released at a constant lytic rate, ϕ (Equations A4, A5, A11, and A12). The number of infections is proportional to the bacterial and viral concentrations. Note that, in this model, we do not consider lysogeny, as it does not contribute directly to the viral shunt. Because we assume that all infections result in direct lysis, infected bacteria are not included as an explicit class in this model. Upon lysis, the intracellular nutrients that were not used for virion production are released in organic form (Equations A7 and A14). Free viruses decay at a constant rate θ.

Nutrient import into the system occurs through an exchange of inorganic nutrients with the subsurface at a rate ω (Equations A6 and A13). Nutrients are exported from the system through the consumption of zooplankton by higher trophic levels (larger zooplankton and fish), encoded following a density-dependent predatory term (i.e., predation rate proportional to the squared density of zooplankton biomass). Part of this consumption is exported as biomass (at a rate p_ex), and the rest is lost to the subsurface in the form of fecal pellets. In order to keep the model tractable, we assumed that all organic matter produced by the modeled zooplankton contribute to DOM (Equations A7 and A14) and can be ultimately reused by bacteria, whereas fecal pellets produced by zooplankton predators contribute to the carbon sink.

The variables and processes allowed us to define and monitor four fluxes, namely productivity, nutrient release, nutrient export to higher trophic levels, and carbon sink. We defined productivity as the total amount of biomass created by heterotrophic bacteria and cyanobacteria (Equation A15). Nutrient release represents the dissolved organic nutrients released by zooplankton and viruses through sloppy feeding and lysis, respectively (Equation A16). The predation term for zooplankton is used for nutrient export to higher trophic levels (Equation A17). Finally, carbon sink is the sum of the refractory carbon produced during lysis and of the sinking fraction of fecal pellets produced by predators of zooplankton (Equation A18). See Appendix A for their corresponding equations.

Aiming to simulate the response of different microbial communities, we used a range of values to parametrize our equations. A literature search provided the minimum and maximum values possible to define a realistic range for each parameter (Appendix C, Table C-1). While the models presented here are general, parameters were selected when possible to represent the conditions in the North Atlantic Ocean—where both N and P limitation can be observed (Ammerman et al., 2003; Mather et al., 2008). Within these ranges, we used a Latin Hypercube Sampling (LHS) method on the log-transformed parameter ranges to generate random parameter sets (McKay et al., 1979). To constrain the trait combinations to realistic ones, we imposed well-known trade-offs for both bacterial and viral traits. For bacteria, we distinguished between fast-growing bacteria (i.e., high maximum growth rate, high half-saturation constant, and low phosphatase production) and slow-growing bacteria (low maximum growth rate, low half-saturation constant, and high phosphatase production) (Klausmeier et al., 2004; Litchman et al., 2007). Viruses were split between large viruses (high adsorption rate, low burst size, and low decay rate) and small viruses (low adsorption rate, high burst size, and high decay rate) (Bongiorni et al., 2005; Weitz, 2015). The range of each constrained parameter was approximately divided in half in log space, with the upper and lower half representing “high” and “low” values, respectively. For each parameter set, we randomly selected one type of bacteria (fast- or slow-growing) and viruses (small or large) and sampled from the corresponding restricted ranges. See Analysis for more details. We did not constrain the host-virus pairings (i.e., assumed no specificity) but, as with the rest of realized communities, we assumed instead as a realistic outcome any equilibrium for the system with coexistence among all organisms.

In order to characterize the realized behavior of each community, we solved analytically the system of equations composing the model to obtain the symbolic equilibrium expressions of each variable (Appendix B). In addition, and using the constraints explained in the previous section, we generated 1,000 random parameter sets (listed in Table C-1), each one representing a community. For each of these parameter sets, we then evaluated the equilibrium expressions, which we used to find the percentage of coexisting communities. In parallel, we also generated 500 random parameter sets to use as initial points for an optimization procedure: following Weitz et al. (2015), we applied a constrained optimization algorithm to find parameter values (within the imposed ranges) that minimized the difference between the models' outputs and values found in the field (Table 1). This optimization algorithm calculates the gradient of a specified function (in this case, a sum of the weighted differences between variable target values and values from the model) at a point in the parameter space, and aims to find the combinations of parameters that minimize this function. To this end, the algorithm takes steps in the direction of the largest negative gradient that are of magnitude proportional to this largest gradient, until a local minimum is reached. See the MATLAB documentation for the fmincon function for further details. We only used the resulting parameter sets if the algorithm found a minimum within the maximum number of steps that were allocated. The optimized parameter distributions across communities are presented in Figure C-1. We then evaluated our analytical expressions using these optimized parameter sets to find the associated equilibrium values for the different variables, and conserved only the parameter sets that resulted in coexistence with and without viruses. Note that this optimization algorithm is computationally expensive, reason why we reduced the number of replicates with respect to the ones generated to study the percentage of coexistence. These two approaches to generate realistic communities allowed us to study the role of viruses for both the N- and P-limited versions of the model. Because the P-limited system without viruses did not allow for coexistence between heterotrophic bacteria and cyanobacteria (see below), we limited the P-limited virus-free version of the model to heterotrophic bacteria only.

In addition to each dynamic variable (bacterial, viral, zooplankton, and organic and inorganic nutrient densities), we also evaluated fluxes of interest: productivity, organic nutrient release, export to higher trophic level, and carbon sink. For each equilibrium with coexistence, we assessed stability by using MATLAB to find numerically the eigenvalues of the Jacobian matrix evaluated at that equilibrium. This exercise allowed us to distinguish between stable nodes (real negative eigenvalues), unstable nodes (at least one real positive eigenvalue), stable spirals (complex eigenvalues with negative real part) and unstable spirals (complex eigenvalues, at least one of which shows a positive real part), the last two being characterized by oscillations around the equilibrium value.

Finally, we used non-optimized equilibrium values to study the relationship between viruses and bacteria, and between viruses and carbon sink. Including non-optimized equilibria allowed for a large diversity of equilibria. Because the expectation for the relationship between bacteria and viruses is a power law i.e., V = αB^γ (Wigington et al., 2016), we calculated the coefficient of the power law that best described this relationship in our simulations. The power law coefficient corresponds to the slope γ of the regression between bacteria and viruses in log space:

In order to obtain an ensemble of regression coefficients, we repeated this process a certain number of replicates. More specifically, we performed 50 linear regressions between bacteria and viruses that used the coexisting communities/equilibria resulting from 200 random parameter sets each (~70 coexisting communities/equilibria per replicate/linear regression). This protocol allowed us to get an average value for the slope and coefficient to ensure that the power law coefficients were representative.

We used a similar procedure to describe the relationship between viruses and carbon sink.

With the N-limited version of the model, coexistence of all populations occurred in 87 out of 1,000 communities (8.7%) for random parameter sets. This number only refers to parameter sets that led to coexistence when used both with and without viruses. As a comparison, coexistence occurred in only 0.45% of randomly generated communities in a similar model comprising an additional class of eukaryotic phytoplankton (Weitz et al., 2015). The omission of autotrophic eukaryotes in our model thus facilitated coexistence, as expected given the reduction in the competition experienced by cyanobacteria. When we evaluated the steady-state concentrations obtained for parameters resulting from the optimization algorithm, coexistence emerged in 131 out of 500 communities (26.2%; see resulting stationary values in Figure 2). Out of these 131 coexisting communities, 70 (53.4%) contained fast-growing heterotrophic bacteria, while the remaining 61 communities were composed of slow-growing heterotrophic bacteria. Approximately 45% of the 131 communities harbored fast-growing cyanobacteria. In addition, 58% of the 131 communities were composed of large viruses that infect heterotrophic bacteria, while 55% of them were composed of large cyanobacteria-infecting viruses. Viruses infecting heterotrophic bacteria were the only group to show a significantly different distribution after optimization and selection for coexistence, compared to the original random growth rate distribution (Chi-squared test, p = 0.0338). There were no significant correlations between the growth rates of heterotrophic bacteria and cyanobacteria (Pearson's correlation test, p = 0.0522), or between the sizes of the two types of viruses (Pearson's correlation test, p = 0.786). Similarly, there was no significant association between bacterial growth rate and the size of their infecting viruses (Pearson's correlation test, H: p = 0.400, C: p = 0.841).

In the virus-free P-limited system, coexistence between heterotrophic bacteria and cyanobacteria was not possible. Without viruses, both types of bacteria rely on the same resource and are preyed upon by the same class of zooplankton. As a consequence, and in agreement with classic resource competition theory, one bacterial class always outcompetes the other. The presence of host-specific viruses prevented competitive exclusion from happening, consistent with the hypothesis that viruses increase bacterial diversity (Thingstad, 2000). In order to make the model with and without viruses as comparable as possible, and noticing that when viruses are present heterotrophic bacteria dominate over cyanobacteria (i.e., have a higher median for abundance, see Table 2, and biomass, see Figure 4), we focused the P-limited case on heterotrophic bacteria only (Figure 3). Importantly, all results were qualitatively identical when using instead a cyanobacteria-only version of the virus-free model (Figure D-1). Out of 1,000 random parameter sets used for the P-limited case, 358 (35.8%) allowed for coexistence of all the populations originally present, both with and without viruses. Reducing the virus-free case to one bacterial class in the P-limited model removed the pressure of bacterial interspecific competition, and explained partially the high rate of coexistence in the P-limited case compared to the N-limited one (100 and 92.8%, respectively). However, the percentage of coexistence was also different between the N- and P-limited versions with viruses (8.7 and 35.8%, respectively), suggesting that intrinsic differences like the possibility to take up also the organic form of the nutrient play a role. Following optimization, we obtained coexistence in 165 out of 500 runs (33%; see resulting stationary values in Figure 3), a percentage that is not significantly different from that of the random sampling (Chi-square test, p = 0.283). This is in contrast with the ~3-fold boost obtained in the N-limited case when using the optimization algorithm, and suggests that the optimization procedure was not significantly more successful than the random sampling in finding feasible equilibria for the P-limited case. Moreover, for some of the 500 optimized parametrizations the optimization function failed to find an optimum in the allocated number of steps, which suggests the existence of regions with low gradient in the parameter space. Out of the 165 coexisting communities, 59 (35.8%) and 63 (38.2%) contained fast growing heterotrophic bacteria and cyanobacteria, respectively. Both of these proportions were significantly different from the random distributions (50.6 and 53.2%) prior to optimization (Chi-squared test, p < 0.001 for both). Coexistence thus occurred more frequently when bacteria grew slowly. On the other hand, the sizes of heterotrophic-bacteria-infecting viruses and cyanobacteria viruses in coexisting communities after optimization did not differ significantly from random distributions (Chi-squared test, p = 0.360 and p = 0.146, respectively). There were significant correlations between the sizes of viruses (Pearson's correlation test, p = 0.392) and between bacterial growth rate and the size of their infecting viruses (Pearson's correlation test, H: p = 0.453, C: p = 0.126). The growth rates of heterotrophic bacteria showed a significant negative correlation with the growth rates of cyanobacteria in coexisting communities (Pearson's correlation test, ρ = −0.456, p = 7.37e−10). In other words, coexistence occurred more frequently when one type of bacteria grew slowly and was efficient at harvesting nutrients, while the other grew faster and was less efficient at getting nutrients.

A majority of coexisting equilibria were also stable. A community that reaches a stable equilibrium will be able to recover from small perturbations (mathematically, all the variables representing the system go back to their equilibrium value following small perturbations). In contrast, if the system is at an unstable equilibrium, it will continue to move away from it when perturbed. Thus, the (ins)stability of an equilibrium is directly related to the resistance and resilience of the community against external perturbations. Regardless of whether the system moves toward or away from equilibrium after a perturbation, this movement can either follow a monotonic path or an oscillatory one, in which case the equilibrium is called a node or a spiral, respectively. Thus, the system can show different degrees of temporal variability when responding to the perturbation. With viruses, we found 100 and 85.5% of stable spirals in the N- and P- limited systems, respectively. In the P-limited system, the remaining equilibria behaved as unstable spirals. Without viruses, stable spirals still dominated (N-limited: 77.9%, P-limited: 95.8%), but stable nodes were also present (N-limited: 7.6%, P-limited: 4.2%). The remaining equilibria in the N-limited system were unstable spirals. There were no unstable equilibria in the P-limited system.

The virus-to-prokaryote ratio (VPR) spanned multiple orders of magnitude in both N- and P-limited versions of the model, and for both heterotrophic bacteria and cyanobacteria (Figure 5). We observed the most diversity in VPR across communities for heterotrophic bacteria in the N-limited system, where the ratio ranged from 0.001 to 1000 with a median close to 10. We applied a log transformation to VPRs to reduce the skew of the data. The standard deviation was 0.98 in log space, corresponding to nearly an order of magnitude. In contrast, the VPRs for cyanobacteria were confined between 10 and 1000 with a median just above 100 and a standard deviation of 0.41 in log space (equivalent to a factor of 2.6). In the P-limited system, the VPR distributions for heterotrophic and cyanobacteria were similar to each other, with most of the communities having a VPR between 1 and 300, a median close to 30 and 100, respectively, and a standard deviation of 0.44 for both (equivalent to a factor 2.8). Although the consensus until recently was that viruses outnumbered bacteria 10 to 1 (i.e., constant VPR = 10), with typical respective abundances of 10⁹ and 10⁸ particles per liter (Bergh et al., 1989; Chibani-chennoufi et al., 2004), recent VPRs gathered from hundreds of studies of marine pelagic environments ranged from 0.0075 to 2150, with a mean of 26.5 (Parikka et al., 2017). Thus, variation of parameters within realistic ranges—representing a diversity of microbial communities and environmental conditions—allowed for wide ranges of VPRs that are consistent with the most current field observations.

Moreover, we found that viral abundance followed a power law as a function of bacterial abundance with an exponent below 1 for both types of bacteria in the P-limited system (H: slope = 0.691, p < 0.001, R² = 0.642; C: slope = 0.817, p < 0.001, R² = 0.68, Figure 6), and for cyanobacteria in the N-limited model (slope = 0.702, p < 0.001, R² = 0.635). In other words, an increase in bacteria was associated with a less-than-proportional increase in viruses. A power law was not a meaningful representation of the relationship between the concentration of heterotrophic bacteria and their viruses in the N-limited case (F-statistic = 0.0246, p = 0.876, R² = 1.96 × 10⁻⁴). Interestingly, heterotrophic bacteria in the N-limited system are the only ones to be limited by an organic nutrient in our model. Both types of bacteria only take up inorganic P, but cyanobacteria only use inorganic N. It is unclear how reliance on inorganic nutrients would result in a significant power-law relationship between bacterial and viral abundance. After generating an additional 50 linear regressions with at least 50 different surviving communities each, the majority of the slopes were between 0.6 and 0.9 (Figure 7), which corroborated the results discussed above. These results are consistent with Wigington et al. (2016), who analyzed data from surveys around the world and found that the relationship between viruses and bacteria is better described by a sublinear power-law than by a linear model.

We observed a four-order-of-magnitude variation across communities in the values of carbon sink that emerged from the model. Variation over multiple orders of magnitude of carbon flux across the world has also been observed empirically, especially close to the surface (Mouw et al., 2016). Importantly, in our N-limited and P-limited systems total viral abundance alone explained 89 and 68% of the carbon sink variation observed with our model, respectively (Figure 8). The linear regressions for both systems were very similar, which suggests that local nutrient limitation does not play a relevant role in this relationship. The associated error differed between the N- and P-limited cases but stayed relatively small for both systems, compared to the variation of the carbon sink values. In the N-limited system, the root mean square error (RMSE) was equal to 0.283 in log-space. This number, which cannot be converted to a single number out of log space (the error is exponential and depends on the predicted value), corresponds to a mean error of 90% above the predicted value and 48% below the predicted value. For a carbon sink of 0.1 μM of carbon per day, this corresponds to a mean error of 0.09 above and 0.048 below 0.1 μM of carbon per day. In the P-limited system, the RMSE was equal to 0.429, corresponding to an error of 168% above and 62% below the predicted value.

For this work, we decided to focus on N and P limitation based on the stoichiometric mismatch between bacteria and viruses, which could affect the proportion of nutrients released during lysis (Jover et al., 2014). However, P limitation and N limitation are only common in subtropical regions. At higher latitudes, iron limitation becomes widespread (Moore et al., 2002; Mather et al., 2008). In addition, certain organisms have specific nutrient requirements, such as diatoms that can be limited by silica (Leynaert et al., 2001; Brzezinski et al., 2003). To understand the effect of viruses globally, we would then need to include these other nutrients in our model. In this context, a more realistic understanding of the system would require including both nutrients simultaneously, potential interactions between them, and the possibility for co-limitation (Bonachela et al., 2015). How to best model both nutrients simultaneously is still an open problem, and beyond the scope of this work. In addition to considering that a single nutrient is limiting at any given time, we made two important simplifications in order to reach analytical expressions for the steady state. First, we did not include a class of infected cells. Depending on the phage and its host, the period of infection can vary from minutes to days (Weitz, 2015), time during which both infected host and virus show a behavior that differs from their free counterparts. This is important because viral genes manipulate the metabolism and nutrient uptake of infected bacteria to facilitate the production of progeny viruses (Zeng and Chisholm, 2012; Puxty et al., 2016), and progeny synthesis in turn depends on the physiological state of the infected host (Hadas et al., 1997; You et al., 2002; Choua and Bonachela, 2019). Furthermore, we determined a fixed nutrient content for each organism, which does not represent accurately the wide variation of cellular elemental composition observed depending on nutrient availability and other environmental conditions (Bonachela et al., 2015). This would, of course, affect the amount of nutrient released during lysis but also, e.g., nutrient uptake dynamics.

Despite these limitations, our model yields variable values consistent with those observed in the field (Table 1) when realistic parameter values are used. In addition, we showed that at least two of its emergent properties (carbon sink and VPR) are within the range observed in empirical studies.

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.