Modeling Unveils How Kleptoplastidy Affects Mixotrophy Boosting Algal Blooms

Abstract

Kleptoplastidy is a nutrition mode in which cells of protists and some multicellular organisms acquire, maintain, and exploit chloroplasts of prey algae cells as photosynthesis reactors. It is an important aspect of the mixotrophic feeding strategy, which plays a role in the formation of harmful algae blooms (HABs). We developed a new mathematical model, in which kleptoplastidy is regarded as a mechanism of enhancing mixotrophy of protists. The model is constructed using three thought (theoretical) experiments and the concept of biological time. We propose to measure the contribution of kleptoplastidy to mixotrophy using a new ecological indicator: the kleptoplastidy index. This index is a function of two dimensionless variables, one representing the ratio of photosynthetic production of acquired chloroplasts versus native chloroplasts, and the other representing the balance between autotrophic and heterotrophic feeding modes. The index is tested by data for the globally distributed, bloom-forming potentially toxic mixotrophic dinoflagellates Prorocentrum cordatum. The model supports our hypothesis that kleptoplastidy can increase the division rate of algae significantly (by 40%), thus boosting their population growth and promoting blooms. The proposed model can contribute to advancements in ecological modeling aimed at forecasting and management of HABs that deteriorate marine coastal environments worldwide.

1. Introduction

Kleptoplastidy is one of the most intriguing though poorly studied aspects of mixotrophic feeding mode, which theoretically can boost reproduction of planktonic algae capable of harmful blooms. As a composite nutrition strategy, mixotrophy is exploited by many protists, including globally distributed dinoflagellates, some of which are toxic or potentially toxic and form harmful algae blooms (HABs), representing a serious threat to sea coastal environments worldwide [1,2,3,4].

Mixotrophy is a combination of autotrophy and heterotrophy. The former is based on photosynthesis and dissolved inorganic nutrients. The latter is realized in two forms. First, osmotrophy, or consumption of dissolved organic nutrients by osmosis. Second, phagotrophy, or consumption of prey cells of algae and/or bacteria by engulfing them. At present, photo-osmo-mixotrophy is viewed as a ubiquitous feeding strategy [4], whereas photo-phago-mixotrophy is less common and its realization varies significantly across the protist spectrum [5]. Once combined, the two strategies give rise to remarkable adaptability of mixotrophic protists to fluctuating environments.

Plankton capable of mixotrophic feeding is referred to as mixoplankton [6]. It is characterized by high environmental tolerance [7], as it can advantageously exploit broad ecological niches including diverse food spectra [8,9,10]. As a result, mixoplankton is often responsible for HABs in aquatic ecosystems [11,12,13,14]. Due to mixoplankton toxicity, or potential toxicity, HABs pose a major global threat, as they negatively affect aquatic communities and ecosystems, fisheries, aquaculture and wildlife, coastal tourism, and human health [15,16,17]. Therefore, the knowledge of kleptoplastidy role in boosting mixotrophy and thus promoting algal blooms is crucial. Meanwhile, for the majority of bloom-forming species, such information is very scarce or lacking.

This paper is concerned with kleptoplastidy (or kleptoplasty) of the protistan mixoplankton, a feeding mode in which a mixotrophic protist cell acquires, maintains, and exploits chloroplasts of prey cells of other algae as photosynthesis reactors [18,19]. Here, we hypothesize that kleptoplastidy can significantly increase the division rate of mixotrophic bloom-forming protists.

Following our perspective on the characterization of cellular and molecular drivers and environmental triggers leading to HABs [10,20,21,22,23], in this paper, we develop a mathematical (as opposed to simulation-based) model of kleptoplastidy. We rely on basic cellular mechanisms for developing the model and use this model for quantifying the contribution of kleptoplastidy to reproductive growth of mixotrophic protists. Finally, we focus on the development of a new ecological indicator, a kleptoplastidy index, as a function of basic eco-physiological variables (cells division rate and consumption rates of different resources). To the best of our knowledge, the proposed model has no counterparts in the pertinent literature, and neither such indicators nor models and experimental data suitable for model calibration and validation are widely available so far [24,25,26,27,28,29,30]. We believe that, as a broader implication, our model may be incorporated in the future ecosystems’ assessment strategies, and the proposed approach can significantly improve existing mixotrophy models, HABs forecasting, and environmental management.

2. Materials and Methods

The mathematical structure of the model is based on extending basic Monod’s equation [31] to a biological system consisting of a single protist, one inorganic resource, and several identical prey cells containing chloroplasts. This approach results in governing equations with two advantageous properties. First, those equations are sufficiently simple to be solved analytically. Second, the equations account for kleptoplastidy via measurable parameters representing essential morpho-physiological features of a mixotrophic protist: cells’ division rate and consumption rates of different resources. A kleptoplastidy model was developed by introducing the concept of biological time and three thought experiments involving controlled environments, which can be realized in laboratory settings only. The first experiment was designed so that only the autotrophic mode was active; in the second, only heterotrophic mode was active; and in the third experiment, all three modes, the autotrophic, heterotrophic, and kleptoplastidic, were active. Consequently, in processing experimental results, it was possible to isolate the kleptoplastidic feeding mode and evaluate its contribution to cells’ growth.

For demonstration purposes, the model is applied to the dinoflagellate Prorocentrum cordatum (Ostenfeld) J.D. Dodge, 1975. This mixotrophic protist has been investigated by several research teams [3,32,33,34,35,36] as it is potentially toxic and globally spread. Recently, P. cordatum was singled out as the most common, globally distributed mixotrophic protist responsible for HABs in 30 countries [3]. A suitable organic resource for studying kleptoplastidy in P. cordatum is the cryptophyte microalga Teleaulax sp. [3] as this organism contains chloroplasts.

Our research aims at the development of a kleptoplastidy model by combining generality and simplicity of models in physics with specificity of biological descriptions. For this reason, the model developed in Section 3.1 and Section 3.2 is for a generic mixotrophic protist. In contrast, in Section 3.3, the model is applied specifically to Prorocentrum cordatum and Teleaulax sp. The basic symbols used in the model and their explanations and units are given in Table 1. Details that are more specific along with the model description are provided in Section 3.1, Section 3.2 and Section 3.3.

3. Results

3.1. Model

3.1.1. Background of the Model

The objective of this section is to develop a basic model of kleptoplastidy for a conventional reference mixotrophic protist. To this end, the protist is viewed simply as a vessel containing native and acquired chloroplasts. Native chloroplasts form a discontinuous belt adjacent to the cell wall interior and function as photosynthesis reactors, using light and inorganic resources. Acquired chloroplasts function as additional photosynthesis reactors, and we refer to them as kleptoplasts. In addition, we identify a free volume available for accumulating kleptoplasts (Figure 1).

A prey cell is viewed as a source of easily digestible organic nutrients and laboriously digestible chloroplasts. Once a prey cell is captured, easily digestible components are rapidly metabolized in food vacuoles, while laboriously digestible chloroplasts become kleptoplasts. This cycle is repeated for subsequently captured prey cells. The ensuing accumulation of kleptoplasts and their use as additional photosynthesis reactors for reproductive growth of the protist is the essence of kleptoplastidy and the subject of our model.

The entire modeling setup is presented in Figure 1. There, the protist is shown as the large irregular centrally located shape, inorganic resources as small peripheral gray circles, and prey cells as large peripheral circles. Native chloroplasts of the protist are shown in blue, while native chloroplasts of prey cells are shown in red, and kleptoplasts are shown in red with dark outline, which represents the additional membrane of the former food vacuole. The free volume available for accumulating kleptoplasts is shown in green (Figure 1).

3.1.2. Monod’s Equation and the Maximum Cell Division Rate

The most fundamental model of reproductive growth of microorganisms is based on Monod’s equation [31]. In particular, if the system of interest consists of organisms (cells) of a single species and a single resource, the cells divide at the rate

(1)$\mu =\Omega {\displaystyle \frac{R}{K+R}}$

Here, $\Omega$ and $K$ are system-specific positive parameters, and $R$ is the environmental resource concentration. This equation implies

(2)$\mu <{\mu }_{\max }=\Omega$

and therefore $\Omega$ is referred to as the maximum cell’s division rate. The parameter $K$ is referred to as the half-saturation constant.

For our purposes, it is useful to introduce a normalized concentration, defined as

(3)$r:={\displaystyle \frac{R}{K}},$

so that (1) can be rewritten as

(4)$\mu =\Omega {\displaystyle \frac{r}{1+r}}.$

We define resource-rich environments in which the normalized resource concentration r >> 1. Then, (4) implies that resource-rich environments are necessary for realizing division rates close to ${\mu }_{\max }=\Omega$.

For a single species and two substitutable resources [37], (4) is generalized as

(5)$\mu ={\Omega }_1{\displaystyle \frac{r_1}{1+r_1}}+{\Omega }_2{\displaystyle \frac{r_2}{1+r_2}}.$

Here, $r_1$($r_2$) is the normalized concentration, defined in (3), and ${\Omega }_1$(${\Omega }_2$) is the maximum cell’s division rate due to the first (second) resource. According to (5), the maximum cell’s division rate is

(6)${\mu }_{\max }={\Omega }_1+{\Omega }_2.$

In this equation, ${\mu }_{\max }$ is associated with the optimal way of consuming the two resources.

In the context of mixotrophy, it is meaningful to regard the resources in (5) and (6) as inorganic and organic. Accordingly, we assign the subscripts 1 and 2 to the inorganic and organic resources, respectively. In what follows, it will be demonstrated that (6) underestimates the true maximum division rate for mixotrophic protists capable of kleptoplastidy. To this end, we will develop a kleptoplastidy model by introducing the concept of biological time and three theoretical thought experiments.

3.1.3. Biological Time

Biological time is a useful concept that allows one to analyze scenarios in which the cells’ division rate is time-dependent, as it turns out to be the case in the proposed kleptoplastidy model. It is defined as a dimensionless function of time, $\psi \left(t\right)$, which satisfies the differential equation

(7)${\displaystyle \frac{d\psi }{dt}}=\mu .$

Further, biological time is confined to the interval $0\le \psi \le 1$, with $\psi =0$ and $\psi =1$ corresponding to the beginning and end of a biological process of interest, respectively. For our purposes, the process of interest occurs between two consecutive divisions of the reference protist cell. Mathematical foundations of this approach have been firmly established in the literature concerned with predicting service life of structural components [38].

As a demonstration, let us consider the most basic situation, where $\mu =\Omega$ is constant:

(8)${\displaystyle \frac{d\psi }{dt}}=\Omega .$

Integration of this equation with respect to time yields

(9)${\displaystyle \underset{0}{\overset{1}{\int }}d\psi }={\displaystyle \underset{0}{\overset{T}{\int }}\Omega dt}.$

Here, $T$ is the time interval between two consecutive divisions. This equation yields

(10)$1=\Omega T$

and

(11)$T={\displaystyle \frac{1}{\Omega }}.$

Of course, the answer for $T$ could have been obtained immediately, without introducing biological time. But this is because, in this example, the right-hand side of (7) is independent of time, $t$. In general, one needs integration, and this will be the case with analysis of kleptoplastidy.

3.1.4. Thought Experiments

In this subsection, we present three thought experiments leading to the kleptoplastidy model. The first two experiments provide essential building blocks, while the third experiment reveals a synergistic contribution of kleptoplastidy.

First Experiment: Autotrophy

Let us consider a single reference mixotrophic protist in an environment characterized by the normalized resource concentrations r₁ >> 1 and $r_2=0$. That is, the environment is inorganic-resource-rich but it contains no organic resources. Suppose that the outcome of this experiment is the division rate ${\mu }^{\left(1\right)}$. Here, we use the superscript (1) to refer to the first thought experiment. For this experiment, (5) implies

(12)${\Omega }_1={\mu }^{\left(1\right)},$

so that the measurement ${\mu }^{\left(1\right)}$ yields the parameter ${\Omega }_1$. Note that, in the chosen environment, the protist behaves like an autotroph.

For future reference, let us define the maximum division rate due to a unit native chloroplast mass as

(13)${\omega }_1:={\displaystyle \frac{{\Omega }_1}{M_1}}.$

Here $M_1$ is the mass of native chloroplasts in the protist. We use the term maximum, as ${\omega }_1$ is directly related to the maximum division rate ${\Omega }_1$.

Second Experiment: Heterotrophy

Let us consider the protist in an environment characterized by $r_1=0$ and r₂ >> 1. That is, the environment contains no inorganic resources, but it is organic-resource-rich. Further, the organic resources are prey cells containing chloroplasts. Suppose that the outcome of this experiment is the division rate ${\mu }^{\left(2\right)}$. Here, we use the superscript (2) to refer to the second thought experiment. For this experiment, (5) implies

(14)${\Omega }_2={\mu }^{\left(2\right)},$

so that the measurement ${\mu }^{\left(2\right)}$ yields the parameter ${\Omega }_2$. Note that, in the chosen environment, the protist behaves like a heterotroph.

Since there are no inorganic resources, kleptoplasts do not contribute to reproductive growth. Furthermore, kleptoplasts may halt reproductive growth by fully occupying the free space (Figure 1), so that the protist can no longer consume prey cells. We assume that this does not happen, and the protist steadily consumes prey cells.

Let us define $\tau$ as the time it takes for the protist to digest one prey cell. During this time, two events occur. First, the protist acquires essentially indigestible native chloroplasts of the prey cell, and they become its kleptoplasts. Second, the protist metabolizes easily digestible remaining nutrients of the prey cell.

In the adopted scenario, acquisition of kleptoplasts is a cumulative process, whose rate is

(15)${\displaystyle \frac{\mathrm{d}{\widehat{M}}_1}{\mathrm{d}t}}={\displaystyle \frac{m_1}{\tau }}.$

Here, ${\widehat{M}}_1$ is the mass of kleptoplasts and $m_1$ is the mass of native chloroplasts in one prey cell. By assuming that the initial mass of kleptoplasts is zero, integration of (15) yields

(16)${\widehat{M}}_1={\displaystyle \frac{m_{1}t}{\tau }}.$

Third Experiment: Mixotrophy

Let us consider the protist in an environment characterized by the normalized resource concentrations r₁ >> 1 and r₂ >> 1. That is, the environment is rich with both inorganic and organic resources. As in the second thought experiment, the organic resource are prey cells containing chloroplasts. Suppose that the outcome of this experiment is the division rate ${\mu }^{\left(3\right)}$. Here, we use the superscript (3) to refer to the third thought experiment. For this environment, (5) implies

${\mu }^{\left(3\right)}={\Omega }_1+{\Omega }_2,$

but this would be incorrect, as (5) does not account for kleptoplastidy. Rather, we write

(17)${\mu }^{\left(3\right)}={\Omega }_1+{\Omega }_2+\widehat{\mu },$

where $\widehat{\mu }$ is the division rate representing photosynthetic activities of kleptoplasts, which were present but inactive in the second thought experiment.

Using results of the first two thought experiments, we can combine (13) and (16) to estimate

(18)$\widehat{\mu }={\omega }_1^{\ast }{\widehat{M}}_1={\omega }_1^{\ast }{\displaystyle \frac{m_{1}t}{\tau }}=\gamma {\omega }_1{\displaystyle \frac{m_{1}t}{\tau }}={\displaystyle \frac{\gamma {\Omega }_1{m}_{1}t}{M_1\tau }}.$

Here, ${\omega }_1^{\ast }$ is the maximum division rate due to a unit kleptoplasts mass, and $\gamma :={\omega }_1^{\ast }/{\omega }_1$. While (18) is useful as it relates kleptoplastidy to morpho-physiological parameters, it may be prohibitively demanding on a required experimental program. Therefore, it is prudent to replace (18) with

(19)$\widehat{\mu }=\beta t,$

and treat $\beta$ as a phenomenological parameter. This parameter can be determined from the measurements using (12), (14), (17), and (19):

(20)$\beta ={\displaystyle \frac{{\mu }^{\left(3\right)}-{\mu }^{\left(1\right)}-{\mu }^{\left(2\right)}}{t}}.$

Of course, if desired, $\beta$ can be related to the morpho-physiological parameters in (18):

(21)$\beta ={\displaystyle \frac{\gamma {\Omega }_1{m}_1}{M_1\tau }}.$

Equations (7), (17) and (19) imply that biological time in the third thought experiment is

(22)${\displaystyle \frac{\mathrm{d}\psi }{\mathrm{d}t}}={\mu }^{(3)}={\Omega }_1+{\Omega }_2+\beta t.$

Upon integration of this equation, as it was carried out in (9), we obtain

(23)$1=\left({\Omega }_1+{\Omega }_2\right)T+{\displaystyle \frac{1}{2}}\beta T^2.$

This quadratic equation for $T$ yields only one positive root, so that the maximum division rate for the protist is

(24)${\Omega }_{1+2}={\displaystyle \frac{1}{T}}={\displaystyle \frac{{\Omega }_1+{\Omega }_2}{2}}\left[\sqrt{1+{\displaystyle \frac{2\beta }{{\left({\Omega }_1+{\Omega }_2\right)}^2}}}+1\right].$

Here, we use the subscript 1+2 to emphasize the synergistic nature of the estimate for ${\Omega }_{1+2}$. Note that ${\Omega }_{1+2}$ is larger than ${\Omega }_1+{\Omega }_2$, and the difference is

(25)$\Delta \Omega ={\Omega }_{1+2}-\left({\Omega }_1+{\Omega }_2\right)={\displaystyle \frac{{\Omega }_1+{\Omega }_2}{2}}\left[\sqrt{1+{\displaystyle \frac{2\beta }{{\left({\Omega }_1+{\Omega }_2\right)}^2}}}-1\right].$

3.2. Kleptoplastidy Index

In this subsection, we consider environments rich in both resources and focus on evaluating the relative contribution of kleptoplastidy to the division rate. By combining (24) and (25), we introduce the kleptoplastidy index κ defined as

(26)$\kappa :={\displaystyle \frac{\Delta \Omega }{{\Omega }_{1+2}}}={\displaystyle \frac{\sqrt{1+\frac{2\beta }{{\left({\Omega }_1+{\Omega }_2\right)}^2}}-1}{\sqrt{1+\frac{2\beta }{{\left({\Omega }_1+{\Omega }_2\right)}^2}}+1}}.$

In what follows, we use this index as the measure of kleptoplastidy.

While it is clear that $0\le \kappa <1$, one can attain a better understanding for the range of $\kappa$ by introducing the dimensionless parameters

(27)$\alpha ={\displaystyle \frac{\beta }{{\Omega }_1{\Omega }_2}}=\gamma \left[\left({\displaystyle \frac{m_1}{\tau {\Omega }_2}}\right)÷M_1\right]$

and

(28)$\sigma ={\displaystyle \frac{4{\Omega }_1{\Omega }_2}{{\left({\Omega }_1+{\Omega }_2\right)}^2}},$

so that (26) takes the form

(29)$\kappa ={\displaystyle \frac{\sqrt{1+\frac{1}{2}\alpha \sigma }-1}{\sqrt{1+\frac{1}{2}\alpha \sigma }+1}}.$

This form is useful for two reasons. First, it involves only dimensionless parameters. Second, both $\alpha$ and $\sigma$ have sound physical and mathematical interpretations.

The physical meaning of $\alpha$ becomes clear once it is noted that in (27) the term in parentheses is the mass of kleptoplasts accumulated in the second thought experiment. Therefore, the ratio in the square brackets is the mass of kleptoplasts acquired in the second thought experiment to the mass of native chloroplasts. Since $\gamma :={\omega }_1^{\ast }/{\omega }_1$, the parameter $\alpha$ is the ratio of the photosynthetic contributions to the division rate of kleptoplasts and native chloroplasts.

The parameter $\sigma$ represents both ${\Omega }_1$ and ${\Omega }_2$. Further, it reaches its maximum, $\sigma =1$, for what one may refer to as perfectly balanced mixotrophy, characterized by ${\Omega }_1={\Omega }_2$. In contrast, for autotrophs $\left({\Omega }_2=0\right)$ or heterotrophs $\left({\Omega }_1=0\right)$, the parameter attains its minimum, $\sigma =0$. Note that the dependence of $\sigma$ on ${\Omega }_1$ and ${\Omega }_2$ is symmetric. This means that $\sigma$ does not favor one feeding mode over the other. Rather, it represents how balanced the two feeding modes are.

In Figure 2, the kleptoplastidy index $\kappa$ is plotted as a function of $\alpha$ and $\sigma$ for $0\le \alpha \le 5$ and $0\le \sigma \le 1$. Thus, the kleptoplastidy index $\kappa$ is an increasing function of $\alpha$ and $\sigma$, and $0\le \kappa <1$. It is clear from (29) and the plot that the maximum $\kappa$ is reached when each $\alpha$ and $\sigma$ attain their respective maximum. In particular, $\kappa \approx 0.3$ for $\alpha =5$ and $\sigma =1$. That is, kleptoplastidy contributes about 30% to the division rate. For small $\alpha$, $\kappa$ is small because only a small mass of kleptoplasts is involved. For small $\sigma$, $\kappa$ is small because balanced mixotrophy is the inherent basis for kleptoplastidy.

3.3. A Case Study

In this section, the developed model is applied to the dinoflagellate Prorocentrum cordatum, as the protist, and a cryptophyte microalga Teleaulax sp., as the prey cell. To emphasize the specificity, we refer to Prorocentrum cordatum as the P-cell rather than just a protist.

3.3.1. Background

Due to its significance to algal blooms in general and often to HABs, the P-cell has been studied by multiple research groups [3,32,33,34,35,36]. Therefore, it is an excellent candidate for a case study in terms of both available data and ecological importance. An electron-microscope image of a P-cell is shown in Figure 3. There, we identify native chloroplasts, a food vacuole, and the nucleus. Figure 3 confirms the schematic in Figure 1, as native chloroplasts form a discontinuous belt adjacent to the cell wall interior, kleptoplasts are stored in food vacuoles or released to the cell interior when vacuoles collapse, and the nucleus exterior is available for accumulating kleptoplasts. The presence of food vacuoles implies consumption of relatively large prey cells via phagotrophy. This is consistent with kleptoplastidy, as prey cells should be large enough to contain potential kleptoplasts.

Mixotrophy of P-cells has been quantified for microalga Teleaulax sp. as prey cells by H.J. Jeong and co-authors [3], who estimated that phagotrophy contributes 15.5 ± 4.7% to the division rate of P. cordatum. Further, the equivalent spherical diameters of the P-cell and prey cell are D = 12.1 µm and d = 5.6 µm, respectively [33]. This means that the P-cell and prey cell are not significantly different in size, and therefore phagotrophy should proceed one prey cell at a time.

3.3.2. The Maximum Division Rate and Kleptoplastidy Index

In this subsection, (24) and (29) are used for evaluating the maximum division rate during mixotrophy, ${\Omega }_{1+2}$, and the kleptoplastidy index, $\kappa$, for the P-cell and cryptophyte microalga Teleaulax sp. as prey cells. To this end, we rely on the values of ${\Omega }_1=0.22$ and ${\Omega }_2=0.38$ divisions per day [32], so that (28) yields

(30)$\sigma ={\displaystyle \frac{4{\Omega }_1{\Omega }_2}{{\left({\Omega }_1+{\Omega }_2\right)}^2}}\approx 0.93.$

Unfortunately, we could not find experimental data for estimating $\alpha$. Nevertheless, one can estimate an upper bound on $\alpha$ directly from (27) by adopting three assumptions. First, we assume that the mass density of every constituent is the same and equal to the mass density of water, so that (27) can be rewritten as

(31)$\alpha =\gamma \left[\left({\displaystyle \frac{v_1}{\tau {\Omega }_2}}\right)÷V_1\right],$

where $v_1$ and $V_1$ are the volumes occupied by the chloroplasts of one prey cell and one P-cell, respectively. Second, the entire volume available for kleptoplasts is fully filled. Third, the kleptoplasts are as effective as native chloroplasts, so that $\gamma =1$. As a result, we obtain

(32)${\alpha }_{\max }={\displaystyle \frac{V_f}{V_1}}\approx 10,$

where $V_f$ is the volume of the P-cell available for accumulating kleptoplasts, shown in green in Figure 1. The estimate $V_f/V_1\approx 10$ is obtained upon visual inspection of Figure 3.

For $\alpha =10$ and $\sigma =0.93$, (24) and (29) yield

(33)${\Omega }_{1+2}={\displaystyle \frac{{\Omega }_1+{\Omega }_2}{2}}\left(\sqrt{1+{\displaystyle \frac{\alpha \sigma }{2}}}+1\right)\approx 1 \mathrm{division} \mathrm{per} \mathrm{day}$

and

(34)$\kappa \approx 0.4.$

These two estimates imply that kleptoplastidy can significantly accelerate reproductive growth of mixotrophic protists.

It is instructive to plot the history of concentration for various conditions. To this end, we assume exponential growth. In Figure 4, we plot the current concentration $c$ normalized by the initial concentration $c_0$ as a function of time for a one-week period; note that the log scale is used along the ordinate. We consider four situations:

(1). The environment is inorganic-resource-rich, but it contains no organic resources. In this case, the division rate is ${\Omega }_1=0.22$ divisions per day.

(2). The environment is organic-resource-rich, but it contains no inorganic resources. In this case, the division rate is ${\Omega }_2=0.38$ divisions per day.

(3). The environment is rich in both resources, but kleptoplastidy is suppressed. In this case, the division rate is ${\Omega }_1+{\Omega }_2=0.60$ divisions per day.

(4). The environment is rich in both resources, and kleptoplastidy is operational. In this case, the division rate is ${\Omega }_{1+2}=1$ division per day.

Thus, Figure 4 clearly demonstrates that kleptoplastidy is a potent mechanism, capable of increasing the cells concentration by three orders of magnitude in environments rich in both resources. In contrast, in environments rich only in one resource, where kleptoplastidy is suppressed, the gain is only one order of magnitude.

4. Discussion

4.1. Interpretation of the Main Findings

In this paper, guided by the topical issue of unveiling cellular and molecular drivers of algal blooms’ formation, we developed a simple mathematical model for mixotrophic protists capable of kleptoplastidy. Since kleptoplastidy has been observed in many aquatic species ranging from microalgae to slugs [39,40,41,42,43], we believe that the model or its variants can be applicable to other aquatic species.

The model was developed using a generalization of Monod’s equations and three thought experiments involving controlled environments, which can be realized in laboratory settings only. The first (second) experiment was designed so that only the autotrophic (heterotrophic) mode was active, and in the third experiment all three modes, the autotrophic, heterotrophic, and kleptoplastidic, were active. As a result, in processing experimental results, it was possible to isolate the kleptoplastidic mode; see (17), (19), and (20). The environments were either resource-rich (r₁, r₂ >> 1) or resource-poor $\left(r_1,r_2=0\right)$, which was sufficient for our purposes because we were focused on the maximum division rates, believed to be associated with HABs. Using controlled environments with moderate resource concentrations $\left(r_1,r_2=O\left(1\right)\right)$ makes modeling more complicated as (5) is no longer reduced to (12) or (14), but rather to equations involving $r_1$ and $r_2$. However, environments with moderate resource concentrations are also beneficial, as they allow one to determine the half-saturation constants $K_1$ and $K_2$ from the first and second thought experiments. Then, the third experiment can be modeled by generalizing (17) as

(35)${\mu }^{\left(3\right)}={\Omega }_1{\displaystyle \frac{R_1}{K_1+R_1}}+{\Omega }_2{\displaystyle \frac{R_2}{K_2+R_2}}+\beta t{\displaystyle \frac{R_1}{K_1+R_1}}{\displaystyle \frac{R_2}{K_2+R_2}}.$

The validity of this equation is impossible to establish without adequate experimental data. But the need for such an equation is apparent, especially for improving existing mixotrophy models and HABs forecasting [24,25,26,27,28,29,30,44].

The model was developed for one protist, one inorganic resource, and one organic resource. Formally, it is straightforward to extend to $N$ protist species, $N_1$ inorganic resource types, and $N_2$ organic resource types. To this end, we rewrite (17) as

(36)${\mu }_i={\displaystyle \sum _{i_1=1}^{N_1}{\Omega }_{1i_1}}+{\displaystyle \sum _{i_2=1}^{N_2}{\Omega }_{2i_2}}+{\displaystyle \sum _{i_2=1}^{N_2}{\displaystyle \sum _{i_1=1}^{N_1}{\beta }_{i_1{i}_2}t}} i=1,\dots ,N.$

For each $i$, this equation contains $N_1$ constants in the first sum for modeling autotrophy, $N_2$ constants in the second sum for modeling heterotrophy, and $N_1{N}_2$ constants in the double sum for modeling kleptoplastidy. Thus, in total, the model involves $N\left(N_1+N_2+N_1{N}_2\right)$ parameters and requires the same number of experiments, all of which are structured like the thought experiments in Section 3.1.4. If desired, one can combine (35) and (36) to derive a generalization of (36) valid for moderate resource concentrations. In this case, the number of parameters and required experiments increases to $N\left(2N_1+2N_2+N_1{N}_2\right)$.

In the second and third experiment, it was assumed that kleptoplastidy is not affected by saturation of the free volume by kleptoplasts. This assumption needs to be questioned, especially considering possibilities of protistan daughter cells inheriting kleptoplasts of their parent cells. If the saturation occurs, one ends up with interesting possibilities and much more significant modeling challenges. For example, a saturated protist may behave as an autotroph, as it has no room for accommodating prey cells. Alternatively, a saturated protist may perform as a mixotroph, but instead of relying on phagotrophy, it switches to osmotrophy and consumption of dissolved organic materials (DOM). Mathematically, such switches in the feeding strategy are associated with nonlinearities. In this regard, the proposed model is applicable to situations occurring well before the saturation point is reached.

The estimate for $\alpha$, adopted for calculating (33) and (34), should be regarded as an upper bound, as it assumes that (i) the free space is saturated with kleptoplasts and (ii) $\gamma =1$. The latter assumption is particularly optimistic, as native chloroplasts are located near the surface, in positions advantageous for photosynthesis (Figure 1 and Figure 3). Thus, kleptoplasts are relegated to positions away from the surface, with inferior conditions for photosynthesis, and therefore one should expect $\gamma <1$.

The kleptoplastidy index $\kappa$ is rather insensitive to small variations in the input parameters $\alpha$ and $\sigma$. One way of establishing this is by showing that the absolute values of the first derivatives of $\kappa$ with respect to these parameters are sufficiently small [45]. In the current setting, this requirement translates into the inequalities $\left|\partial \kappa /\partial \alpha \right|<1$ and $\left|\partial \kappa /\partial \sigma \right|<1$. For simplicity, we show this for $\sigma =1$, as this restriction is the most important for kleptoplastidy. Upon differentiation of $\kappa$, we obtain

(37)${{\displaystyle \frac{\partial \kappa }{\partial \alpha }}}_{|\sigma =1}={\displaystyle \frac{1}{2\sqrt{1+\frac{\alpha }{2}}{\left(\sqrt{1+\frac{\alpha }{2}}+1\right)}^2}} \mathrm{and} {{\displaystyle \frac{\partial \kappa }{\partial \sigma }}}_{|\sigma =1}={\displaystyle \frac{\alpha }{2\sqrt{1+\frac{\alpha }{2}}{\left(\sqrt{1+\frac{\alpha }{2}}+1\right)}^2}}.$

The bounds for each derivative can be found by finding its minimum and maximum. Accordingly, using basic calculus, we determine that $0\le \partial \kappa /\partial {\alpha }_{|\sigma =1}\le 1/8$ and $0\le {\partial \kappa /\partial \sigma }_{|\sigma =1}<0.172$. Thus, the absolute values of both derivatives are less than one, and therefore the model is stable.

The estimates in Section 3.3 are limited to environments rich in both inorganic and organic resources. While their generalization to environments with moderate resource concentrations needs to be supported experimentally, we can evaluate situations when the environment changes from being rich only in one resource to being rich in both resources. In particular, if the environment changes from being inorganic-resource-rich to being rich in both resources, then the division rate switches from ${\Omega }_1=0.22$ to ${\Omega }_{1+2}=1$ divisions per day. This is more than 4.5 times increase. Similarly, if the environment changes from being organic-resource-rich to being rich in both resources, then the division rate switches from ${\Omega }_2=0.38$ to ${\Omega }_{1+2}=1$ divisions per day. This is more than 2.5 times increase. Both increases are significant, especially if the reproductive growth process is exponential. These increases in the cells’ division rate may result in dramatic mixoplankton population growth, as shown in Figure 4.

4.2. Advantages and Limitations of the Study

The results of this study contribute to the theoretical biology by mathematical modeling of kleptoplastidy, regarded as a key component of the mixotrophic feeding strategy in protists [39,40,41,42,43,44]. This reveals another aspect explaining high adaptability of such organisms to fluctuating environments [1,3,4,5,6], due to wider trophic niches of mixoplankton compared to photoautotrophic plankton [4,5,8,9,10,34,35,36,44]. In addition, this allows one to regard kleptoplastidy as a harmful bloom forming mechanism [30]. To the best of our knowledge, the proposed model is the only one, which clearly demonstrates how kleptoplastidy accelerates the growth of mixotrophic algae populations. The model focuses on two quantities of interest. First, the division rate ${\Omega }_{1+2}$ of bloom-forming mixotrophic algae capable of kleptoplastidy. This quantity is essential for evaluating and forecasting algal blooms in general, and HABs in particular. Second, the kleptoplastidy index $\kappa$, or the relative contribution of kleptoplastidy to ${\Omega }_{1+2}$, which allows one to assess how strongly bloom-forming mixotrophic algae react to changes in nutrients’ composition and abundance in the environment. Therefore, the first numerical assessment of the kleptoplastidy contribution to population growth of bloom-forming mixotrophic protists by calculating the kleptoplastidy index is a new and advantageous approach to comprehension, evaluation and forecasting algal blooms in general, and HABs in particular.

The limitations of the proposed model are due to scarcity of experimental data on growth rates of bloom-forming mixotrophic protists in environments suitable for mixotrophy [44]. For this reason, our case study (Section 3.3) was carried out using only one bloom-forming potentially toxic dinoflagellate species, Prorocentrum cordatum, as the protist, and the cryptophyte microalga Teleaulax sp., as the prey cell. These calculations are based on several literature sources [3,32,33] that, to our knowledge, contain the most complete and representative experimental data on the physiological characteristics and mixotrophic feeding of P. cordatum, our major species of interest [8,20,21,34,35,36,44]. Therefore, we believe that usage of these particular sources of data secured the most adequate results of our case study. Consequently, the validation of our analytical model is the task for future experimental investigations.

4.3. Future Perspectives

The issue of the productive capabilities of organisms is one of the key questions in modern ecology [2,3,4,5,6,11,12]. A number of different disciplines in natural sciences require information on fine mechanisms that drive life in aquatic ecosystems [44]. Specifically, the processes that ensure high feeding efficiency are among basic prerequisites of successful performance of organisms, their populations and communities. This can be proven by a large number of publications demonstrating the enhancement and worldwide spreading of HABs in marine coastal ecosystems and their dependence on trophic and physical-chemical conditions varying in space and time [2,10,11,12,13,14]. Algal blooms may result in significant shifts in community structure and functions, and as such changes can occur without a single dominant ecosystem driver, emergent ecosystem-level responses are difficult to predict [22,30]. Generally, there is a fair understanding of how different components interrelate in aquatic ecosystems in terms of who-eats-whom [46,47]. Meanwhile, the inter- and intracellular interactions within microbial communities largely remain enigmatic, as well as the question of how exactly some feeding strategies are implemented.

Our results offer a new perspective on the study of microbial mixotrophy as a HAB-forming mechanism [30] and a factor potentially influencing the productivity of phytoplankton and entire aquatic communities. In this study, kleptoplastidy as a primary though poorly studied aspect of mixotrophy was for the first time modeled, aiming at evaluation of its effectiveness. This aspect so far has never been taken into account in mathematical models of mixotrophic population growth and dynamics [19,24,26,28,44,47]. The proposed model and the new kleptoplastidy index may provide quantitative assessment of the functional role of mixotrophic protists in aquatic ecosystems. As a broader implication, the future incorporation of the kleptoplastidy model into larger ecosystem models will improve forecasts of algal blooms, including HABs, and contribute to the development of new monitoring strategies and effective environmental management of marine coastal ecosystems.

5. Conclusions

In this study, a basic mathematical model of kleptoplastidy was first developed. The model demonstrates how kleptoplastidy may significantly accelerate the reproductive growth of bloom-forming mixotrophic protists and thus lead to algal blooms and, in the cases of toxic or potentially toxic protists, to HABs. The model showed that kleptoplastidy in mixotrophic protists is a synergistic process, requiring both inorganic and organic resources. Further, the model is developed for resource-rich environments, under the provision that such environments are necessary for HABs.

The new model gives rise to the kleptoplastidy index κ, a new ecological indicator, which allows one to assess how kleptoplastidy contributes to proliferation of mixotrophic protists. This indicator is a function of the parameters α and σ. The former represents the photosynthetic contribution of kleptoplasts normalized by that of native chloroplasts. The latter represents the balance between the autotrophic and heterotrophic feeding modes.

The kleptoplastidy phenomenon and possibility of its ecological indication in terms of the kleptoplastidy index κ was illustrated by using existing experimental data for potentially toxic bloom-forming mixotrophic dinoflagellates Prorocentrum cordatum feeding on the cryptophyte microalga Teleaulax sp. For this system, we calculated and compared the cells’ division rates in various resource-rich environments. Our modeling results support the hypothesis that kleptoplastidy, as a feeding mode, can increase the division rate of mixotrophic protists significantly (by up to 40%), thus boosting their population growth and, consequently, promoting algal blooms.

The new developments can contribute to advancements in mathematical modeling of cellular mechanisms that underpin basic ecological processes of utmost importance—the formation of HABs that deteriorate marine and coastal environments worldwide. Incorporation of the kleptoplastidy model into comprehensive ecosystem models will improve HABs’ forecasts and future monitoring programs thus securing effective management of aquatic ecosystems.

Author Contributions

Conceptualization, I.V.T., G.J.R., H.S. and S.O.S.; methodology, I.V.T., G.J.R., H.S. and S.O.S.; software, G.J.R.; validation, I.V.T., G.J.R., H.S. and S.O.S.; formal analysis, G.J.R.; investigation, I.V.T., G.J.R., H.S. and S.O.S.; resources, H.S. and S.O.S.; data curation, H.S.; writing—original draft preparation, I.V.T. and G.J.R.; writing—review and editing, I.V.T., G.J.R., H.S. and S.O.S.; visualization, I.V.T. and G.J.R.; project administration, S.O.S.; funding acquisition, S.O.S. and H.S. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The study does not report any data.

Acknowledgments

I.V.T. acknowledges the facilities of the Zoological Institute RAS available within the government assignment topic # 125012800888-5.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviation

The following abbreviation is used in this manuscript:

HABs	Harmful algae blooms

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1 The conventional reference mixotrophic protist (large irregular shape), inorganic resources (small grey circles), and prey cells (large circles). Native chloroplasts of the protist are shown in blue, native chloroplasts of prey cells are shown in red, and kleptoplasts are shown in red with dark outline. The volume available for acquiring kleptoplasts is shown in green.

Figure 2 The kleptoplastidy index $\kappa$ as a function of the dimensionless parameters $\alpha$ and $\sigma$: (A) 3-D plot, (B) contour plot (color gradient and numbers at the lines indicate the $\kappa$ values: <0.05 in blue color, >0.25 in light yellow). The parameter $\alpha$ represents the photosynthetic contribution of kleptoplasts normalized by that of native chloroplasts. The parameter $\sigma$ represents the balance between the autotrophic and heterotrophic feeding modes.

Figure 3 An electron-microscopy image of Prorocentrum cordatum (photo courtesy of M.A. Berdieva). Abbreviations: FV—food vacuole, N—nucleus, ChP—native chloroplasts of the P-cell.

Figure 4 The normalized concentration of protist cells $c/c_0$ (log scale) as a function of time: (1) the environment is inorganic-resource-rich, but it contains no organic resources, (2) the environment is organic-resource-rich, but it contains no inorganic resources, (3) the environment is rich in both resources, but kleptoplastidy is suppressed, and (4) the environment is rich in both resources, and kleptoplastidy is operational.

Table 1

Nomenclature.

Symbol	Description	Units
$\mu$	A generic division rate of the protist	day−1
${\mu }^{\left(1\right)}$	The protist division rate measured in the first thought experiment	day−1
${\mu }^{\left(2\right)}$	The protist division rate measured in the second thought experiment	day−1
${\mu }^{\left(3\right)}$	The protist division rate measured in the third thought experiment	day−1
$\Omega$	A generic maximum division rate of the protist	day−1
${\Omega }_1$	The maximum protist division rate due to inorganic resources	day−1
${\Omega }_2$	The maximum protist division rate due to organic resources	day−1
${\Omega }_{1+2}$	The maximum protist division rate due to inorganic and organic resources in the presence of kleptoplastidy	day−1
$R$	The concentration of a generic resource	µg L–1
$R_1$	The concentration of an inorganic resource	µg L–1
$R_2$	The concentration of an organic resource	µg L–1
$K$	A generic half-saturation constant	µg L–1
$K_1$	The half-saturation constant for inorganic resources	µg L–1
$K_2$	The half-saturation constant for organic resources	µg L–1
$r$	The normalized concentration of a generic resource	none
$r_1$	The normalized concentration of an inorganic resource	none
$r_2$	The normalized concentration of an organic resource	none
$\psi$	Biological time	none
t	Time	days
T	The time interval between two divisions of the protist	days
τ	The time interval for digesting one prey cell by the protist	days
$m_1$	The mass of chloroplasts in one prey cell	µg
$M_1$	The mass of native chloroplasts in the protist	µg
${\omega }_1$	The maximum division rate of the protist due to a unit native chloroplast mass	µg–1 day–1
${\omega }_1^{\ast }$	The maximum division rate of the protist due to a unit kleptoplast mass	µg–1 day–1
$\gamma$	The maximum division rate of the protist due to a unit kleptoplast mass normalized by that of the native chloroplasts	none
$\beta$	The phenomenological parameter representing the acceleration of the maximum division rate due to kleptoplastidy	µg–1 day–2
k	Kleptoplastidy index	none
α	The photosynthetic contribution of kleptoplasts normalized by that of native chloroplasts	none
σ	The balance parameter between the autotrophic and heterotrophic feeding modes of the reference protist	none
$v_1$	The volume of the chloroplasts of one prey cell	µm3
$V_1$	The volume of the chloroplasts of the protist	µm3

References

1. Stoecker, D.K.; Hansen, P.J.; Caron, D.A.; Mitra, A. Mixotrophy in the marine plankton. Ann. Rev. Mar. Sci.; 2017; 9, pp. 311-335. [DOI: https://dx.doi.org/10.1146/annurev-marine-010816-060617] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27483121]

2. Glibert, P.M. Harmful algae at the complex nexus of eutrophication and climate change. Harmful Algae; 2020; 91, 101583. [DOI: https://dx.doi.org/10.1016/j.hal.2019.03.001] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32057336]

3. Jeong, H.J.; Kang, H.C.; Lim, A.S.; Jang, S.H.; Lee, K.; Lee, S.Y.; Ok, J.H.; You, J.H.; Kim, J.H.; Lee, K.H. . Feeding diverse prey as an excellent strategy of mixotrophic dinoflagellates for global dominance. Sci. Adv.; 2021; 7, 4214. [DOI: https://dx.doi.org/10.1126/sciadv.abe4214] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33523999]

4. Mitra, A.; Caron, D.A.; Faure, E.; Flynn, K.J.; Leles, S.G.; Hansen, P.J.; McManus, G.B.; Not, F.; Gomes, H.D.R.; Santoferrara, L.F. . The Mixoplankton Database–diversity of photo-phago-trophic plankton in form, function and distribution across the global ocean. J. Eukar. Microbiol.; 2023; 70, e12972. [DOI: https://dx.doi.org/10.1111/jeu.12972] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36847544]

5. Mitra, A.; Caron, D.A.; Faure, E.; Flynn, K.J.; Leles, S.G.; Hansen, P.J.; McManus, G.B.; Not, F.; Gomes, H.D.R.; Santoferrara, L.F. . Defining planktonic protist functional groups on mechanisms for energy and nutrient acquisition; incorporation of diverse mixotrophic strategies. Protist; 2016; 167, pp. 106-120. [DOI: https://dx.doi.org/10.1016/j.protis.2016.01.003] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26927496]

6. Flynn, K.J.; Mitra, A.; Anestis, K.; Anschütz, A.A.; Calbet, A.; Ferreira, G.D.; Gypens, N.; Hansen, P.J.; John, U.; Martin, J.L. . Mixotrophic protists and a new paradigm for marine ecology: Where does plankton research go now?. J. Plankton Res.; 2019; 41, pp. 375-391. [DOI: https://dx.doi.org/10.1093/plankt/fbz026]

7. Glibert, P.M. More than propagule pressure: Successful invading algae have physiological adaptations suitable to anthropogenically changing nutrient environments. Aquat. Ecosyst. Health Manag.; 2015; 18, pp. 334-341. [DOI: https://dx.doi.org/10.1080/14634988.2015.1027137]

8. Telesh, I.V.; Schubert, H.; Skarlato, S.O. Ecological niche partitioning of the invasive dinoflagellate Prorocentrum minimum and its native congeners in the Baltic Sea. Harmful Algae; 2016; 59, pp. 100-111. [DOI: https://dx.doi.org/10.1016/j.hal.2016.09.006] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28073501]

9. Telesh, I.; Schubert, H.; Skarlato, S. Ecological niches of bloom-forming cyanobacteria in brackish Baltic Sea coastal waters. Estuar. Coast. Shelf Sci.; 2023; 295, 108571. [DOI: https://dx.doi.org/10.1016/j.ecss.2023.108571]

10. Telesh, I.; Schubert, H.; Skarlato, S. Wide ecological niches ensure frequent harmful dinoflagellate blooms. Heliyon; 2024; 10, e26495. [DOI: https://dx.doi.org/10.1016/j.heliyon.2024.e26495] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38404903]

11. Anderson, C.R.; Berdalet, E.; Kudela, R.M.; Cusack, C.K.; Silke, J.; O’Rourke, E.; Dugan, D.; McCammon, M.; Newton, J.A.; Moore, S.K. . Scaling up from regional case studies to a global Harmful Algal Bloom observing system. Front. Mar. Sci.; 2019; 6, 250. [DOI: https://dx.doi.org/10.3389/fmars.2019.00250]

12. Carstensen, J.; Conley, D.J.; Almroth-Rosell, E.; Asmala, E.; Bonsdorff, E.; Fleming-Lehtinen, V.; Gustafsson, B.G.; Gustafsson, C.; Heiskanen, A.-S.; Janas, U. . Factors regulating the coastal nutrient filter in the Baltic Sea. Ambio; 2020; 49, pp. 1194-1210. [DOI: https://dx.doi.org/10.1007/s13280-019-01282-y] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31707582]

13. Hallegraeff, G.M.; Anderson, D.M.; Belin, C.; Bottein, M.-Y.D.; Bresnan, E.; Chinain, M.; Enevoldsen, H.; Iwataki, M.; Karlson, B.; McKenzie, C.H. . Perceived global increase in algal blooms is attributable to intensified monitoring and emerging bloom impacts. Commun. Earth Environ.; 2021; 2, 117. [DOI: https://dx.doi.org/10.1038/s43247-021-00178-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37359131]

14. Karlson, B.; Andersen, P.; Arneborg, L.; Cembella, A.; Eikrem, W.; John, U.; West, J.J.; Klemm, K.; Kobos, J.; Lehtinen, S. . Harmful algal blooms and their effects in coastal seas of Northern Europe. Harmful Algae; 2021; 102, 101989. [DOI: https://dx.doi.org/10.1016/j.hal.2021.101989] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33875185]

15. Brown, A.R.; Lilley, M.; Shutler, J.; Lowe, C.; Artioli, Y.; Torres, R.; Berdalet, E.; Tyler, C.R. Assessing risks and mitigating impacts of harmful algal blooms on mariculture and marine fisheries. Rev. Aquacult.; 2020; 12, pp. 1663-1688. [DOI: https://dx.doi.org/10.1111/raq.12403]

16. Brown, A.R.; Lilley, M.; Shutler, J.; Lowe, C.; Artioli, Y.; Torres, R.; Berdalet, E.; Tyler, C.R. Marine harmful algal blooms (HABs) in the United States: History, current status and future trends. Harmful Algae; 2021; 102, 101975. [DOI: https://dx.doi.org/10.1016/j.hal.2021.101975] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33875183]

17. Rattner, B.A.; Wazniak, C.E.; Lankton, J.S.; McGowan, P.C.; Drovetski, S.V.; Egerton, T.A. Review of harmful algal bloom effects on birds with implications for avian wildlife in the Chesapeake Bay region. Harmful Algae; 2020; 120, 102319. [DOI: https://dx.doi.org/10.1016/j.hal.2022.102319] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36470599]

18. Rumpho, M.E.; Dastoor, F.P.; Manhart, J.R.; Lee, J. The Kleptoplast. The Structure and Function of Plastids; Wise, R.R.; Hoober, J.K. Advances in Photosynthesis and Respiration Springer: Dordrecht, The Netherlands, 2007; Volume 23, [DOI: https://dx.doi.org/10.1007/978-1-4020-4061-0_23]

19. Millette, N.C.; Gast, R.J.; Luo, J.Y.; Moeller, H.V.; Stamieszkin, K.; Andersen, K.H.; Brownlee, E.F.; Cohen, N.R.; Duhamel, S.; Dutkiewicz, S. . Mixoplankton and mixotrophy: Future research priorities. J. Plankton Res.; 2023; 45, pp. 576-596. [DOI: https://dx.doi.org/10.1093/plankt/fbad020] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37483910]

20. Skarlato, S.; Filatova, N.; Knyazev, N.; Berdieva, M.; Telesh, I. Salinity stress response of the invasive dinoflagellate Prorocentrum minimum. Estuar. Coast. Shelf Sci.; 2018; 211, pp. 199-207. [DOI: https://dx.doi.org/10.1016/j.ecss.2017.07.007]

21. Skarlato, S.O.; Telesh, I.V.; Matantseva, O.V.; Pozdnyakov, I.A.; Berdieva, M.A.; Schubert, H.; Filatova, N.A.; Knyazev, N.A.; Pechkovskaya, S.A. Studies of bloom-forming dinoflagellates Prorocentrum minimum in fluctuating environment: Contribution to aquatic ecology, cell biology and invasion theory. Protistology; 2018; 12, pp. 113-157. [DOI: https://dx.doi.org/10.21685/1680-0826-2018-12-3-1]

22. Telesh, I.V.; Schubert, H.; Joehnk, K.D.; Heerkloss, R.; Schumann, R.; Feike, M.; Schoor, A.; Skarlato, S.O. Chaos theory discloses triggers and drivers of plankton dynamics in stable environment. Sci. Rep.; 2019; 9, 20351. [DOI: https://dx.doi.org/10.1038/s41598-019-56851-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31889119]

23. Telesh, I.; Schubert, H.; Skarlato, S. Abiotic stability promotes dinoflagellate blooms in marine coastal ecosystems. Estuar. Coast. Shelf Sci.; 2021; 251, 107239. [DOI: https://dx.doi.org/10.1016/j.ecss.2021.107239]

24. Glibert, P.M.; Allen, J.I.; Bouwman, A.F.; Brown, C.W.; Flynn, K.J.; Lewitus, A.J.; Madden, C.J. Modeling of HABs and eutrophication: Status, advances, challenges. J. Mar. Syst.; 2010; 83, pp. 262-275. [DOI: https://dx.doi.org/10.1016/j.jmarsys.2010.05.004]

25. Mitra, A.; Flynn, K.J.; Burkholder, J.M.; Berge, T.; Calbet, A.; Raven, J.A.; Granéli, E.; Glibert, P.M.; Hansen, P.J.; Stoecker, D.K. . The role of mixotrophic protists in the biological carbon pump. Biogeosciences; 2014; 11, pp. 995-1005. [DOI: https://dx.doi.org/10.5194/bg-11-995-2014]

26. Berge, T.; Chakraborty, S.; Hansen, P.; Andersen, K. Modeling succession of key resource-harvesting traits of mixotrophic plankton. ISME J.; 2017; 11, pp. 212-223. [DOI: https://dx.doi.org/10.1038/ismej.2016.92] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27482925]

27. Fernandes-Salvador, J.A.; Davidson, K.; Sourisseau, M.; Revilla, M.; Schmidt, W.; Clarke, D.; Miller, P.I.; Arce, P.; Fernández, R.; Maman, L. . Current status of forecasting toxic harmful algae for the North-East Atlantic shellfish aquaculture industry. Front. Mar. Sci.; 2021; 8, 666583. [DOI: https://dx.doi.org/10.3389/fmars.2021.666583]

28. Horemans, D.M.L.; Friedrichs, M.A.M.; St-Laurent, P.; Hood, R.R.; Brown, C.W. Forecasting Prorocentrum minimum blooms in the Chesapeake Bay using empirical habitat models. Front. Mar. Sci.; 2023; 10, 1127649. [DOI: https://dx.doi.org/10.3389/fmars.2023.1127649]

29. Medvinsky, A.B.; Nurieva, N.I.; Adamovich, B.V.; Radchikova, N.P.; Rusakov, A.V. Direct input of monitoring data into a mechanistic ecological model as a way to identify the phytoplankton growth-rate response to temperature variations. Sci. Rep.; 2023; 13, 10124. [DOI: https://dx.doi.org/10.1038/s41598-023-36950-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37349488]

30. Cagle, S.E.; Roelke, D.L. Chaotic mixotroph dynamics arise with nutrient loading: Implications for mixotrophy as a harmful bloom forming mechanism. Ecol. Model.; 2024; 492, 110714. [DOI: https://dx.doi.org/10.1016/j.ecolmodel.2024.110714]

31. Monod, J. The growth of bacterial cultures. Annu. Rev. Microbiol.; 1949; 3, pp. 371-394. [DOI: https://dx.doi.org/10.1146/annurev.mi.03.100149.002103]

32. Heil, C.A. Influence of humic, fulvic and hydrophilic acids on the growth, photosynthesis and respiration of the dinoflagellate Prorocentrum minimum (Pavillard) Schiller. Harmful Algae; 2005; 4, pp. 603-618. [DOI: https://dx.doi.org/10.1016/j.hal.2004.08.010]

33. Jeong, H.J.; Yoo, Y.D.; Kim, J.S.; Seong, K.A.; Kang, N.S.; Kim, T.H. Growth, feeding and ecological roles of the mixotrophic and heterotrophic dinoflagellates in marine planktonic food webs. Ocean Sci. J.; 2010; 45, pp. 65-91. [DOI: https://dx.doi.org/10.1007/s12601-010-0007-2]

34. Matantseva, O.; Skarlato, S.; Vogts, A.; Pozdnyakov, I.; Liskow, I.; Schubert, H.; Voss, M. Superposition of individual activities: Urea-mediated suppression of nitrate uptake in the dinoflagellate Prorocentrum minimum revealed at the population and single-cell levels. Front. Microbiol.; 2016; 7, 1310. [DOI: https://dx.doi.org/10.3389/fmicb.2016.01310] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27610101]

35. Matantseva, O.; Pozdnyakov, I.; Voss, M.; Liskow, I.; Skarlato, S. The uncoupled assimilation of carbon and nitrogen from urea and glycine by the bloom-forming dinoflagellate Prorocentrum minimum. Protist; 2018; 169, pp. 603-614. [DOI: https://dx.doi.org/10.1016/j.protis.2018.05.006] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30096707]

36. Khanaychenko, A.N.; Telesh, I.V.; Skarlato, S.O. Bloom-forming potentially toxic dinoflagellates Prorocentrum cordatum in marine plankton food webs. Protistology; 2019; 13, pp. 95-125. [DOI: https://dx.doi.org/10.21685/1680-0826-2019-13-3-1]

37. León, J.A.; Tumpson, D.B. Competition between two species for two complementary or substitutable resources. J. Theor. Biol.; 1975; 50, pp. 185-201. [DOI: https://dx.doi.org/10.1016/0022-5193(75)90032-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/1127957]

38. Kachanov, L.M. Introduction to Continuum Damage Mechanics; 1st ed. Springer: Dordrecht, The Netherlands, 1986; [DOI: https://dx.doi.org/10.1007/978-94-017-1957-5]

39. Nishitani, G.; Nagai, S.; Hayakawa, S.; Kosaka, Y.; Sakurada, K.; Kamiyama, T.; Gojobori, T. Multiple plastids collected by the dinoflagellate Dinophysis mitra through kleptoplastidy. Appl. Environ. Microbiol.; 2011; 78, 813. [DOI: https://dx.doi.org/10.1128/AEM.06544-11] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22101051]

40. Stamatakis, K.; Vayenos, D.; Kotakis, C.; Gast, R.J.; Papageorgiou, G.C. The extraordinary longevity of kleptoplasts derived from the Ross Sea haptophyte Phaeocystis antarctica within dinoflagellate host cells relates to the diminished role of the oxygen-evolving Photosystem II and to supplementary light harvesting by mycosporine-like amino acid/s. Biochim. Biophys. Acta; 2017; 1858, pp. 189-195. [DOI: https://dx.doi.org/10.1016/j.bbabio.2016.12.002]

41. Cruz, S.; Cartaxana, P. Kleptoplasty: Getting away with stolen chloroplasts. PLoS Biol.; 2022; 20, e3001857. [DOI: https://dx.doi.org/10.1371/journal.pbio.3001857] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36346789]

42. Jesus, B.; Jauffrais, T.; Trampe, E.C.; Goessling, J.W.; Lekieffre, C.; Meibom, A.; Kuhl, M.; Geslin, E. Kleptoplast distribution, photosynthetic efficiency and sequestration mechanisms in intertidal benthic foraminifera. ISME J.; 2022; 16, pp. 822-832. [DOI: https://dx.doi.org/10.1038/s41396-021-01128-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34635793]

43. Sørensen, M.E.; Zlatogursky, V.V.; Onuţ-Brännström, I.; Walraven, A.; Foster, R.A.; Burki, F. A novel kleptoplastidic symbiosis revealed in the marine centrohelid Meringosphaera with evidence of genetic integration. Curr. Biol.; 2023; 33, pp. 3571-3584. [DOI: https://dx.doi.org/10.1016/j.cub.2023.07.017] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37536342]

44. Skarlato, S.O.; Telesh, I.V. Kleptoplastidy as a principal driver of mixotrophy in bloom-forming dinoflagellates. Protistology; 2024; 18, pp. 286-292. [DOI: https://dx.doi.org/10.21685/1680-0826-2024-18-4-4]

45. Cacuci, D.G. Sensitivity and Uncertainty Analysis: Theory; Chapman & Hall: New York, NY, USA, 2003; Volume 1, [DOI: https://dx.doi.org/10.1201/9780203498798]

46. García-Oliva, O.; Wirtz, K. The complex structure of aquatic food webs emerges from a few assembly rules. Nat. Ecol. Evol.; 2025; 9, pp. 576-588. [DOI: https://dx.doi.org/10.1038/s41559-025-02647-1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/40021903]

47. Frank, A.S.; Pinke, A.; Read, E.; Röblitz, S.; Sadhu, S.; Subbey, S. Exploring multistability in marine ecosystems: Insights from food web model dynamics. Ecol. Model.; 2025; 508, 111245. [DOI: https://dx.doi.org/10.1016/j.ecolmodel.2025.111245]