Physics or biology? Persistent chlorophyll accumulation in a shallow coastal sea explained by pathogens and carnivorous grazing

One of the most striking patterns at the land–ocean interface is the massive increase of chlorophyll-a (CHL) from continental shelves towards the coast, a phenomenon that is classically linked to physical features. Here I propose that the coastal–offshore CHL gradient in a shallow sea has biological origins related to phytoplankton mortality that are neglected in state-of-the-art biogeochemical models. I integrate a trait-based ecosystem model into a modular coupling framework that is applied to the southern North Sea (SNS). The coupled model very well reproduces daily, seasonal and inter-annual (2000-2014) dynamics and meso-scale patterns in macronutrients, zooplankton biomass, and CHL as observed in situ and by remote sensors. Numerical experiments reveal that coast–offshore CHL gradients may predominantly arise from a trophic effect as resolved by an increase in carnivorous grazing towards shallow waters. This carnivory gradient reflects higher near-coast abundance of juvenile fish and benthic filter feeders. Furthermore, the temporal evolution of CHL can be much affected by viral infection as a fast-responding loss process at intermediate to high phytoplankton concentrations. Viral control in the model also prevents excessive and unrealistic blooms during late spring. Herbivores as often only ecological factor considered for explaining the spatio-temporal phytoplankton distribution are in this study supplemented by pathogens as well as pelagic and benthic carnivores as powerful agents, which are barely represented in current modeling but can mediate physical drivers of coastal ecosystems.

Remote sensing chlorophyll-a (CHL) plotted against in situ data from Dutch and German marine monitoring stations in the southern North Sea. The circle diameter indicates the time difference between the two measurements from below 1 day (larger circles) to 3 days (smallest circles). A non-linear transformation of the ESACCI CHL seeks to better fit the in situ data while allocating equal weight to each station, thus trying to minimize the error to the cluster center of each station data set (broad grey line): CHL=0.35·10 1.2c 3.1 /(8+c 3 ) with c denoting the ESACCI CHL in mg-Chl m −3 .

A. Benthic biogeochemistry
A variant of OMEXDIA (Soetaert et al., 1996) brought into the FABM format simulates diagenetic processes in the sea floor. The original version of the diagenesis model calculates a simplified carbon and nitrogen cycle, primarily consisting of respiration of benthic particulate organic matter (bPOC), nitrification and denitrification as visualized in Fig 2. In addition, the P cycle is represented by the turnover of a dissolved species, bPO4, and a particulate pool bPP comprising both inorganic and organic forms. Adsorption of bPO4 (transformation into the particulate pool bPP) ceases when a critical capacity of electron donators (C PAds , see Table A in S1) is reached, which in turn is given by the difference between oxygen demand units, bODU (a pool of reduced substances, e.g., iron and hydrogen sulfide), and dissolved oxygen concentration, bO2 (benthic oxygen concentration): with the sigmoid response function σ s (x) defined in Eq.(3). The time derivative of bPP is equal to Eq.(1) with negative sign. Numerical instabilities in the original variant of OMEXDIA are prevented through an alternative management of the oxygen dynamics. Immediate oxygen demand in OMEXDIA reflects the concentration of NH3, high quality bPOC, and bODU. Oxygen consumption therefore contains aerobic respiration of high and low quality bPOC, nitrification, as well as oxidation of reduced substances, and depends on a saturation function of bO2 with process specific half-saturation constants. These constants are in the new model version raised by a small relaxation parameter times the concentration of the competing oxygen pathways (e.g., bPOC and bODU for NH3 oxidation). This scheme downscales consumption rates under high competition of electron donators at low levels of bO2, which in the numerical integration scheme prevents the tendency towards negative bO2 under such conditions.

B. Viral dynamics: infection, replication, and mortality
The dynamics in intracellular viral density v is here described by three aggregate processes, (1) infection-replication, (2) virus removal by (selective) host reactions, and (3) virus mortality.
(1) Infection comprises multiple stages which entail processes outside and inside the host organisms. Viral adsorption rate at the cell membrane r ads is proportional to a specific rate r * ads (likely influenced by water turbulence), and to the bulk virus concentration in the water, which is in the model proportional to the internal density v times the phytoplankton concentration PhyC, thus v · PhyC, and to the ratio of cross-sections describing the chance of hitting an autotrophic host. The cross-section ratio relates the probability for the virus to collide with a phytoplankton cell (increasing with PhyC) to the collision probability for other particles such as bacteria, aggregates, or lithogenic material.
Hence, high concentrations of detritus protect phytoplankton cells, in line with the observed reduction in infection rate with increasing abundance of bacteria and colloids (Murray, 1995). However, the real effect might be stronger as the lithogenic contribution is neglected in the POC calculation. After infection, intracellular viral replication using either RNA or DNA of the host cell depends on temperature (Brussaard, 2004) and the carrying virus capacity of the (decaying) cell, which analogue to Eq.(6) is a smooth step function σ s (v max − v) with capacity v max . The same non-linear form is here assumed to apply to the dependency on the physiological state of the host cell as expressed by its C:N:P stoichiometry (Hadas et al., 1997;Clasen and Elser, 2007;Birch et al., 2012) as proxy for protein and RNA/DNA turnover. The C:N:P stoichiometry is in MAECS given by relative N-and P-quotas (q N and q P , see Wirtz and Kerimoglu (2016)). Replication then should cease for low N-and P-quotas and reach a maximum for sufficiently high quotas. Together, we have with maximal burst size n * v . The latter can be combined with the specific adsorption rate r * ads into the viral maturation rate r * v .
(2) Viral density within a host population declines due to endogenous mortality. The corresponding rate r mort increases with temperature but ceases at very low density (v < v low ), with which the model imposes a persistent background or latent virus concentration in the environment, The specific rate r * mort can be estimated from viral residence times in the surface ocean reported to be on the order of 1 day (Suttle and Chen, 1992).
(3) Viral density in addition decreases because of antiviral defense such as preferential decline of more infected hosts. I here introduce a mathematical expression for preferential host mortality, which should also describe the more general case of host defense. Infected phytoplankton often displays apoptosis, i.e. programmed cell death (Bidle and Falkowski, 2004). As a consequence of the survival of healthy or less infected individuals or groups/species, average pathogen density v decreases. Preferential mortality of infected individuals and concomitant "cleaning" of the population can be understood as competitive process. Within the general formalism of competitive trait changes introduced by Wirtz and Eckhardt (1996) , the resulting decline rate, r defense , is proportional to the diversity of infection degrees within the host population (δ v ) and the marginal growth change of hosts due to viral disease. This marginal growth change is given by the (negative) derivative of the phytoplankton mortality m P,vir with respect to viral density v, thus -∂m P,vir /∂v, or with Eq.(6), m * P,vir s σ 2 s e s · (1 − v) =m P,vir s σ s e s · (1 − v) . Diversity of infection degrees has been derived for (trait) variables strictly bounded between zero and v max to be proportional to v · (v max − v) (Wirtz and Eckhardt, 1996). According to our definition of v, such a biomechanical upper bound will be around v max =2. In addition, the diversity in host infection degrees depends on the structure of the host community. A community dominated by few species with relatively homogeneous infection level, displays a small diversity, whereas coexistence of a broad variety of species and thus pathogenic affinities and histories will increase infection diversity. This specificity is here estimated refering to studies on biomass-diversity relationships. For example, the data of phytoplankton species diversity presented by Interlandi and Kilham (2001) can be well fitted by an inverse function of total phytoplankton biomass (here e −PhyC/C low ), such that both host and disease diversity will be low at high total biomass, which will render countermeasures such as apoptosis more effective. Finally, I assume any defense mechanism to cease at very low (but non-zero) viral density v low like for virus mortality above. Together, virus removal by preferential hosts mortality, is (5) Notably, Eq.(5) can also be thought to describe the expression of defense machinery or occurrence of virophages; anti-viral defense will increase with the marginal impact of the disease. This response will be most pronounced at lethal infection level v=v max /2=1 where the defense rate r defense is around 2m P,vir at high host biomass.      (left) is compared to remote sensing data from ESACCI (center) and the scenario of lacking viral infections (right). In contrast to other maps, CHL is plotted using a linear (non-logarithmic) scale.