Acquired Phototrophy and Its Implications for Bloom Dynamics of the Teleaulax-Mesodinium-Dinophysis-Complex

Introduction

Dinophysis is a harmful dinoflagellate of global distribution (Figure 1) that is the causative agent for diarrhetic shellfish poisoning (DSP) events leading to the closure of aquaculture facilities (review by Reguera et al., 2012, 2014). Significant effort has been expended on understanding the dynamics behind the occurrence of these, and other harmful algal bloom (HAB), events to find ways to forecast them (Raine et al., 2010; Velo-Suárez et al., 2010; Ajani et al., 2016; Moita et al., 2016; Pinto et al., 2016); many monitoring programs target toxigenic Dinophysis species (Campbell et al., 2010). Traditional methods of monitoring Dinophysis by counting cells and using satellite chlorophyll concentrations are not robust for forecasting blooms of these species (Ruiz-Villarreal et al., 2016). This is because the cell concentration threshold for Dinophysis, at which this organism can pose a threat, is only 200 cells L^–1 (Yasumoto et al., 1985), a cell abundance level that can be very difficult to detect (Berdalet et al., 2017). Satellite imagery is inappropriate because blooms often occur in subsurface water layers and are, therefore, not detectable with this method (e.g., Dinophysis acuta; Ruiz-Villarreal et al., 2016).

Modeling offers tools to forecast blooms and also to aid understanding of the dynamics behind their occurrences (Cusack et al., 2016; Ruiz-Villarreal et al., 2016; Flynn and McGillicuddy, 2018). Most models for Dinophysis bloom formation emphasize the importance of physical oceanographic drivers (Raine et al., 2010; Velo-Suárez et al., 2010; Moita et al., 2016; Pinto et al., 2016), with a particular focus on upwelling and downwelling processes and coastal advection that bring Dinophysis cells into coastal waters from offshore, or into surface waters from deeper water layers (Table 1). However, there is still a lack of predictive power in the existing models suggesting that there may be one or more key features missing from these models, for example, the biological characteristics (Ajani et al., 2016; Moita et al., 2016).

An important aspect of Dinophysis is that most of the toxigenic species (Anderson et al., 2002; Park et al., 2006; Jaén and Mamán, 2009; Nishitani et al., 2010; Rodríguez et al., 2012) are mixoplanktonic (Leles et al., 2017). Mixoplankton are protists that employ photosynthesis and phagotrophy for their nutritional needs (Flynn et al., 2013, 2019). There are important differences between mixoplankton functional types which greatly complicate the mechanistic modeling of these organisms (Mitra and Flynn, 2010; Mitra et al., 2016). Dinophysis sp. is plastidic specialist non-constitutive mixoplankton (pSNCM, defined by Mitra et al., 2016, modified by Flynn et al., 2019) that acquire their phototrophic ability from ingestion of ciliates from the Mesodinium genus (e.g., Mesodinium rubrum, Park et al., 2006; Mesodinium major, Rial et al., 2015; Mesodinium chamaeleon, Moeller et al., 2021; Mesodinium coatsi, Nam et al., 2015). Mesodinium is itself a pSNCM and acquires phototrophic potential from consumption of cryptophytes from the Teleaulax-Plagioselmis-Geminigera (TPG) clade (Gustafson et al., 2000; Johnson et al., 2006; Hernández-Urcera et al., 2018). However, there are records of Mesodinium and, as a result, Dinophysis (using Mesodinium as an intermediate), acquiring phototrophy using cryptophytes outside of the TPG clade (Díaz et al., 2020; Moeller et al., 2021). The cryptophyte Teleaulax are constitutive mixoplankton (CM; Mitra et al., 2016; Flynn et al., 2019), having their own innate photosystems, and can feed on bacterioplankton and the cyanobacterium Synechococcus (Yoo et al., 2017). While both the pSNCM (ciliate and dinoflagellate) are capable of maintaining their acquired phototrophy (kleptoplastids) for some time, ultimately they rely on this trophic linkage to proliferate (Hansen et al., 2013). Despite this strict linkage within the TMD-complex, while the global distribution and occurrence of Dinophysis are well-recorded, the distribution of Mesodinium is less well-documented, and this is even less so for Teleaulax (Figure 1). However, due to the species-specific plastid requirements of Dinophysis, the distribution of Mesodinium and also Teleaulax may be expected to be as widespread as Dinophysis.

Blooms of Mesodinium and Dinophysis can occur a couple of weeks apart, assuming that sufficient prey in the form of a Mesodinium bloom is available to support the initial Dinophysis growth (Moita et al., 2016). This is possible because Mesodinium retains the prey nuclei for acquired phototrophy with a half-life of 10 days or so and can maintain plastids potentially for months (Park et al., 2008). Accordingly, it has been suggested (Harred and Campbell, 2014) to use the presence of M. rubrum as an indicator for blooms of Dinophysis ovum in the Gulf of Mexico. There have been suggestions of including Mesodinium and their cryptophyte prey in the modeling of Dinophysis blooms (e.g., Glibert et al., 2010). While the conceptual model for the Rias in NW Spain combines physical oceanographic drivers with Mesodinium-Dinophysis occurrence to improve predictions of Dinophysis blooms (Velo-Suárez et al., 2014), first and foremost such blooms depend on the unique relationship the Dinophysis have with their prey. To further explore the dynamics of this biological complex formed of Teleaulax-Mesodinium-Dinophysis (hereafter TMD-complex) requires a simulation model that focuses on the unique physiology of the three organisms. Hitherto, however, there has been no attempt to mechanistically model these interactions.

In this study, we present the first system dynamics model describing the functioning of the TMD-complex. For the in silico experiments, the Dinophysis component was configured as Dinophysis acuminata because it is a well-studied common toxigenic species. There are relatively more physiological data available for this particular species to aid model configuration. Teleaulax amphioxeia and Mesodinium rubrum were chosen as species representatives for the cryptophyte and Mesodinium components, respectively, as these are associated with acquired phototrophy of D. acuminata. However, it must be noted that the TMD-complex model is designed on the ecophysiological understanding of these organisms; being a system biology-based platform, the model can be configured for different TMD species. We conducted in silico experiments using this configuration to explore the dynamics of the TMD-complex in response to predator-prey ratios, under different environmental conditions, such as nutrient concentrations, irradiance, and temperature. This food-web simulator focuses on the effect of the unique biology and interactions between the different mixoplankton types of the complex. From a conceptual angle, the model thus also proves a platform for improving our understanding of any future new-found mixoplankton interactions involving the exchange of phototrophic potential. Furthermore, the relatively simple construct of the TMD-complex model offers the potential for future coupling of this to a hydrodynamic model.

Materials and Methods

What follows is a general description of the simulation model. The complete equation set of the model and additional information are provided in the Supplementary Material.

Model Structure

The TMD food-web simulator is a development of the N-based system dynamics model of mixoplankton functional types described by Anschütz and Flynn (2020). For this study, the original model was extended to describe three protist functional types, namely, two pSNCM (configured for Mesodinium and Dinophysis) and one CM for Teleaulax. From here on, these organisms are referred to by their genus name only except when discussing specific species. Generic phytoplankton (termed “Alga”; equivalent spherical diameter, ESD = 6 μm) was included as an additional prey for Mesodinium and to provide competition to the Teleaulax for nutrients and light in phototrophic growth.

The full construct is shown in Figure 2. The food-web comprised the CM Teleaulax feeding only on the cyanobacteria Synechococcus (Yoo et al., 2017), the pSNCM Mesodinium feeding on the CM Teleaulax, Synechococcus, and the generic phytoplankton Alga (Myung et al., 2006; Jeong et al., 2015), while Dinophysis fed only on Mesodinium. Synechococcus is configured as a generic species representative of a picophytoplankton competitor and potential prey for Teleaulax. Dissolved inorganic nutrients (DIN) were described via two-state variables representing nitrate (NO₃^–) and ammonium (NH₄⁺). All the plankton functional types could uptake NH₄⁺. All except Dinophysis could also use NO₃^– as this genus seems unable to use this N-source (Hattenrath-Lehmann et al., 2015; García-Portela et al., 2020). However, the construct of the Dinophysis model is based on the generic construct of a pSNCM (Anschütz and Flynn, 2020) and, therefore, should evidence of Dinopysis using NO₃^– become available, then the model can be readily modified to enable NO₃^– uptake for Dinophysis just as it is enabled for the other pSNCM, Mesodinium. Nutrient regeneration occurs in the simulator by an implicit regeneration term. Bacteria and the microbial loop components are not explicitly described.

The operational maximum growth rate for each organism (UmT) was calculated with reference to the current temperature (T), the reference temperature (RT) at which the reference growth rate was declared (UmRT), and a Q₁₀ factor.

A constant loss rate through mortality (0.05 day^–1) was applied to all organisms. In the case of Mesodinium and Dinophysis, this constant mortality was in addition to the mortality rate, described in the original model, accounting for a failure in photosystem activity (Anschütz and Flynn, 2020).

In addition to the two-state variables describing DIN (i.e., nitrate and ammonium), two additional state variables were added describing labile and semi-labile dissolved organic nitrogen (lDON and slDON; Figure 2). These accounted for outflows of voided material from predation, releasing organic nitrogen due to mortality and DON leakage concurrent with photosynthesis (Alldredge et al., 1993; Bronk, 1999; Passow, 2002). Both lDON and slDON were decomposed to ammonium at a rate of 0.5 day^–1. The exact value of this regeneration did not affect the outcomes of simulations. Cyanobacteria, including Synechococcus (Wawrik et al., 2009), and most if not all protist plankton (Flynn and Butler, 1986; Flynn et al., 2019) are capable of exploiting lDON to support their growth. Thus, we enabled all organisms to use lDON with a modification of the routine used by Fasham et al. (1990). Details of model modifications and additions are provided in the Supplementary Material.

Model Parameterization

In this study, the model presented was configured specifically to represent D. acuminata, assuming a nitrogen content of 2,318 pg N cell^–1 (cellular carbon content after Menden-Deuer and Lessard, 2000; cellular N:C ratio of D. acuminata after Rao and Pan, 1993). D. acuminata was chosen for model configuration because it is a common toxigenic species (Figure 1), and relatively more physiological data are available to aid model configuration. T. amphioxeia and M. rubrum were chosen as representatives for the cryptophyte Teleaulax and the ciliate Mesodinium, respectively, for the same reasons of data availability for the configuration of the model and also as the species associated with acquired phototrophy in D. acuminata. If the model was to be used to describe other species, then values for constants controlling physiological features would need to align with traits of those species (Supplementary Data Sheet 1 Table S1).

All constant parameters defining the activity of Synechococcus sp., T. amphioxeia, M. rubrum, and D. acuminata were accorded values consistent with literature data (Table 2). Cellular nitrogen to carbon ratios (N:C) vary among species, with the availability of external nutrients and (for predators) with prey quantity and quality. In this N-based model, the assumed N:C value was constant with the value based on literature data of organisms grown in balanced nutrient conditions. In Anschütz and Flynn (2020), the prey-predator encounter model assumed swimming speeds calculated with reference to cell radius (Flynn and Mitra, 2016). As the motilities (swimming speed; v) of organisms in the TMD-complex organisms are well-studied (Jiang and Johnson, 2021; Supplementary Data Sheet 1 Table S1), and in the case of M. rubrum is also unique with a high escape swimming speed (Riisgård and Larsen, 2009), the auxiliary controlling motility (in Anschütz and Flynn, 2020) was changed to a constant (v, m s^–1) reflecting literature values (Table 2). The two pSNCM organisms, M. rubrum and D. acuminata, have different abilities to maintain optimal functionality of their photosystems which they acquire from their prey; accordingly, these were configured using literature values (critIR_Meso and critIR_Dino for M. rubrum and D. acuminata respectively, Table 3).

Dynamic Sensitivity Analysis

The dynamic sensitivity analysis was performed in the same manner as described by Anschütz and Flynn (2020). The conditions under which the food-web simulator was run are listed in Supplementary Data Sheet 1 Table S2, and the parameters tested for the sensitivity analysis are listed in Supplementary Data Sheet 1 Table S3.

In silico Experiment Scenarios

The TMD-complex food-web simulator was developed and run within Powersim Studio 10¹ using an Euler integration routine of step size 0.0625 day^–1. Simulations were run with default conditions (unless otherwise noted) of 15°C and surface PFD 500 μmol m^–2 s^–1. Other abiotic conditions were set as described in Tables 3, 4. Steady-state and dynamic sensitivity analyses were conducted as in the study by Anschütz and Flynn (2020).

For most simulations, the mixed layer depth (MLD) was set to 10 m. Dinophysis often forms high-density blooms in thin subsurface layers as a result of stratification (Gentien et al., 1995; Anderson et al., 2005; Farrell et al., 2012; Moita et al., 2016). To simulate the physical conditions of a thin layer for subsurface blooms of D. acuminata (e.g., Baldrich et al., 2021), the MLD was set to 1 m and irradiance (de facto at the top of that thin layer) to 50 μmol m^–2 s^–1 to account for the light attenuation in the overlying water.

In most of the scenarios, Dinophysis was not subjected to grazing by a higher predator (i.e., there was no closure function applied to this submodel; cf. Mitra, 2009). This was because this study explores the specific interactions within the TMD-complex rather than the complete ecosystem. To test the impact of a grazer of Dinophysis on the overall bloom dynamics of the TMD-complex, a simple Dinophysis specific grazing function (closure term) implicitly describing activity by a higher predator was implemented. In silico experiments at different nutrient loadings were conducted with and without this grazing function. All other simulations were run without this grazing function.

The control (default) scenario was run at eutrophic conditions (700 mgN m^–3), 15°C, 10 m MLD, and PFD 500 μmol photons m^–2 s^–1. In all other treatments, one of these factors was changed (Table 4). Relative changes in peak biomass and changes in bloom duration were calculated with reference to the control. The organisms, namely, Synechococcus sp., generic phytoplankton “Alga,” Teleaulax, Mesodinium, and Dinophysis, were inoculated at biomass abundances equivalent to 2, 2, 1.5, 0.1, and 0.01% of initial DIN concentrations, respectively. The biomasses of all organisms were described in nitrogen biomass units (mg N m^–3). The complete list of configurations for each in silico experiment is given in Table 4.

To aid contextual understanding, Dinophysis blooms are expressed as cells L^–1 in addition to biomass. For this transformation, cellular (i.e., biomass) carbon was computed from cellular nitrogen assuming cellular C:N as per the Redfield ratio. Then, cellular carbon was converted to ESD using the equations of Menden-Deuer and Lessard (2000). These values were then used to compute cells L^–1.

Results

Biotic Interactions

To provide a reference point for comparisons between different test scenarios, a summary of results is shown in Figure 3. All simulations shared a common pattern of organism succession where Synechococcus was grazed down early to low residual biomass, followed by growths of “Alga” and Teleaulax which were then followed by Mesodinium (Figures 4–8). The accumulating Mesodinium biomass not only outcompeted these species for nutrients but also exerted grazing pressure on them. The ciliate always generated the highest biomass. Dinophysis growth started during the active growth phase of Mesodinium and then exerted grazing pressure on Mesodinium accelerating its eventual decline. Ultimately, both Mesodinium and Dinophysis declined as their growth rates decreased due to the subsequential loss of their source of acquired phototrophy over time. The top-down control on Mesodinium by Dinophysis would alleviate the grazing pressure of the ciliates on the cryptophytes; this enables the cryptophytes to bloom again.

Ammonium was used preferentially over nitrate, along with lDON (Figures 4–8). Both lDON and slDON accumulated and ammonium began to increase due to the decomposition of these organic forms. High levels of ammonium later in the simulations could not be used by Mesodinium and Dinophysis because the ultimate source of their chloroplasts, Teleaulax, was depleted and, thus, the two pSNCM species lost their ability to photosynthesize, and therefore, could no longer exploit inorganic nutrients (DIN). Depending on the timing of events within the 90-day simulation period, in some simulations, this collapse of the pSNCM allowed Synechococcus, Alga, and Teleaulax to bloom again.

By the time Dinophysis blooms, nitrate concentrations were already much lower or depleted due to prior use by the other organisms (Figures 3, 4). However, as the (initial) nutrient load greatly affected the bloom development of those earlier organisms, it affected the timing and magnitude of the much later occurring Dinophysis bloom (Figure 4). Also, the earlier plankton activity caused ammonium to accumulate, which Dinophysis can use while it cannot use nitrate. Oligotrophic conditions (i.e., DIN < 5 μM) did not enable any significant growth of Dinophysis nor of any blooms (Figures 3, 4). Teleaulax and Mesodinium consumed the population of Synechococcus. Teleaulax and Alga bloom together and were followed by Mesodinium. Mesotrophic conditions in the same depth gave rise to a bloom of Mesodinium as well as Dinophysis (Figures 3, 4). Inorganic nutrient concentrations declined most rapidly when Mesodinium reached large biomasses. Eutrophic conditions resulted in the highest cell numbers of Dinophysis (Figure 4). More nutrients drove earlier, sharper, and higher peak blooms.

In this study, most of the simulations presented featured no grazing pressure on Dinophysis, in order to focus on the trophodynamics of the TMD complex. In reality, Dinophysis experiences grazing pressure from copepods and other protozooplankton, such as Fragilidium duplocampanaeforme (Kozlowsky-Suzuki et al., 2006; Lee and Park, 2017). Figure 4 shows the consequence of including a generic grazing function for Dinophysis. The bloom dynamics of the complex including biomass yield of Teleaulax and Mesodinium and the timing of the blooms were affected very little by predators grazing on Dinophysis. The Mesodinium bloom lasted slightly longer in eutrophic conditions. The Dinophysis blooms began at the same time but had a much lower biomass yield that is closer to the cell numbers observed in the field.

The proportions and concentrations of each organism at the start of simulations affected both the maximum cell numbers of all organisms and bloom dynamics in terms of bloom peak and duration to varying degrees (Figures 3, 5). Initial higher proportions of Synechococcus greatly increased the biomass of Teleaulax and to a small extent also that of Mesodinium while it had a negligible effect on Dinophysis growth (Figures 3, 5). Small initial concentrations of Mesodinium in proportion to Dinophysis resulted in the largest biomass of Mesodinium and the smallest biomass of Dinophysis (Figure 3), as Dinophysis could not develop and, in fact, initially declined (died) due to low prey availability. In consequence, Dinophysis could not exert significant grazing pressure on Mesodinium later on in the simulation. In contrast, a high initial concentration of the ciliate promoted a more rapid sequence of events that eventually (because both Mesodinium and Dinophysis blooms crashed) led to a later bloom of Synechococcus that used the available regenerated ammonium. High concentrations of Dinophysis inoculum promoted a quick bloom succession with smaller maximum biomasses of all organisms and no succeeding blooms (Figures 3, 5).

Abiotic Interactions

Higher irradiance increased the peak size of blooms of each organism type and also prolonged the bloom duration (Figure 3), but the timing of their blooms did not change (Figure 6). In contrast, temperature increases had a large effect on the outcome of the Dinophysis bloom with the potential to accelerate the appearance of its bloom as well as decrease its duration (Figures 3, 7). In addition, these factors promoted earlier blooms of Dinophysis and shorter growth cycles, but they could also lower the maximum cell numbers. The lower the temperature, the more the peak of Dinophysis was delayed (several days or weeks). The cell numbers of the Dinophysis bloom varied with the phase of the Mesodinium bloom in which Dinophysis began to grow.

In the comparison of the simulations at 10 m MLD vs. 1 m MLD, the maximum biomass of all organisms was lower at 10 m MLD, except Synechococcus which was grazed away early in all instances (Figure 8). The duration of the individual blooms remained broadly similar, but in 10 m MLD nitrate concentrations did not decline to the same low level as in the 1 m depth simulation. The predator-prey encounter rates were higher in 1 m depth. While the encounter rates between Teleaulax and Synechococcus are roughly the same in 1 m and 10 m depths, the encounter rate between Mesodinium and Teleaulax was almost double within the 1 m thin layer simulation. The encounter between Dinophysis and Mesodinium was also higher in 1 m depth and began to increase earlier than in 10 m.

Overall, Teleaulax peak biomass increased under any conditions that decreased the biomass of Mesodinium (Figure 3). While higher temperature increased the peak biomass of Teleaulax, it decreased that of Mesodinium and Dinophysis. Higher irradiance raised the peak biomass of all three organisms. Changes in the peak biomass of Teleaulax and Mesodinium were often coupled with the opposite effect in bloom duration. For example, while lower nutrient concentrations decreased and higher irradiances increased peak biomasses, bloom duration increased and decreased, respectively (Figure 3). In Dinophysis, there was much less variation in bloom duration compared with Teleaulax and Mesodinium.

Discussion

Understanding the bloom dynamics of the TMD complex is greatly complicated by the multitude of biotic and abiotic factors. These include water column stratification, salinity, advection of organisms and nutrients, changes in predator-prey interactions, and the mixoplanktonic nature of the organisms. In this study, we focused specifically on the interaction between the organisms of the TMD-complex itself. Thus, we assumed a simplified mesocosm-like situation to constrain impacts of hydrodynamics and higher predatory activity. We took this approach because, ultimately, blooms of Dinophysis can only develop if the other members of the TMD-complex have already proliferated and are present in the environment. Thus, Dinophysis depends mostly on the availability of Mesodinium at the right time. Likewise, Mesodinium depends on the cryptophyte or similar prey for its acquired phototrophy.

In our simulations, abiotic factors, directly and indirectly, shaped the timing and magnitude of Dinophysis blooms as they influenced the development of the Mesodinium bloom from which Dinophysis itself started to proliferate. The effects of the initial conditions accumulate over time and lead to, sometimes, counter-intuitive results. The collective interactions of the organisms of prey and predator and competitors for nutrients alongside nutrient regeneration all impacted the dynamics. Consistent with Gaedke et al. (2010) who found biotic interactions in phytoplankton blooms to be highly complex with the potential to overrule climate change effects, our results indicate that the biotic interactions between the mixoplankton species in the TMD-complex to be a key factor in their response to environmental changes and thus need to be considered in predictions about the effect on climate change on this plankton.

The organisms in this model were configured as Teleaulax amphioxeia, Mesodinium rubrum, and Dinophysis acuminata as those were the TMD-complex species with the most physiological data available in published literature. However, the model can be configured to represent any species of the genus Teleaulax, Mesodinium, or Dinophysis, provided there are available data. The dynamic risk analysis showed the model to be sensitive to variations in the maximum growth rate (which is sensitive to temperature) as well as the pSNCM maximum plastid retention time (crit_IR) and the optimal prey capture rates (Supplementary Data Sheet 1 Figure S1). Different species of Dinophysis are known to vary in their physiology and ecology and toxin production (Table 2 and Supplementary Data Sheet 1 Table S1; Reguera et al., 2012; García-Portela et al., 2018). Also, the plastid retention time may differ between species of Dinophysis. With more data becoming available on different species of the TMD-complex, the model could be used for exploring the potential implications of trait differences between the different species.

Apart from physiological and morphological variations, species of Dinophysis differ in seasonality, geographical distribution (Figure 1), and the hydrodynamics that support their proliferation (Reguera et al., 2012; Cusack et al., 2016; García-Portela et al., 2018). For example, in Australia, seasonal stratification and nutrient loads were found to be triggers for blooms of D. acuminata, while growths of Dinophysis caudata were linked to nutrients, salinity, and dissolved oxygen (Ajani et al., 2016). The importance of the TMD-complex is seen associated with a common trait of most species of Dinophysis in which they seem more dependent on the availability of their prey Mesodinium than on light availability (Reguera et al., 2012). Below, we discuss different facets of the interactions as driven by abiotic and biotic forcing as revealed in our simulations.

Predator-Prey Ratio Within the Teleaulax-Mesodinium-Dinophysis-Complex

Above all else, the ratio of encounters between the constituents of the TMD-complex affect the dynamics. Unlike CM (which may not need to feed regularly, if at all) or generalist non-constitutive mixoplankton (GNCM; which need to feed very frequently), pSNCM species, such as Mesodinium and Dinophysis, need to feed regularly but not necessarily frequently on their plastid donors. The initial ratio of the organisms at the start of simulations had a considerable impact on the subsequent bloom dynamics (Figures 3, 5). The concentrations of Mesodinium in comparison to the other organisms appeared to have the biggest impact on overall bloom dynamics as the ciliate is both an important grazer on Teleaulax as well as being the specialist prey for acquired phototrophy to Dinophysis. However, the results were sometimes counter-intuitive as low initial concentrations of Mesodinium could result in comparatively low biomasses of Dinophysis while the ciliate itself subsequently grew to larger numbers (Figure 3). If the Mesodinium started out in low numbers in comparison to Dinophysis, the dinoflagellate was starved of sufficient plastids to support its population early on. At the same time, Teleaulax could develop larger cell numbers due to a lower grazing pressure by Mesodinium. By the time the cell numbers of Mesodinium increased, the ciliate found itself in perfect conditions with high biomass of its prey Teleaulax and low grazing pressure by Dinophysis. While M. rubrum is considered non-toxic, it is known to form large blooms (red-colored tides) that can, on collapse, lead to oxygen depletion and thus suffocation of other organisms (Yoo et al., 1998). Dinophysis acuminata is thus dependent both on the concentration of M. rubrum as well as the timing of its appearance. Prey availability early on in the growth phase of Dinophysis in the simulations appears to be crucial as the pSNCM will be able to maintain its chloroplasts for a time. In this situation, the relative and absolute amounts of the TMD-complex species are affected by nutrient loading, irradiance, and temperature.

Nutrient Loading

Our results suggest that eutrophic conditions increase the potential for Dinophysis blooms because these conditions prepare the system for larger Mesodinium blooms which then increase the potential for Dinophysis blooms (Figures 3, 4). The importance of prey abundance for Dinophysis has been observed in the field (Moita et al., 2016; Smith et al., 2018), and DSP events have been related to high levels of dissolved organic and inorganic nutrients (Hattenrath-Lehmann et al., 2015; Hattenrath-Lehmann and Gobler, 2015). In the simulations, by the time Dinophysis starts to bloom, levels of nitrate could already be very low with ammonium starting to accumulate as a result of decomposition of lDON and slDON; this supported the growth of Dinophysis that had access to plenty of Mesodinium as a source for plastids. Dinophysis itself appears either unable or poorly equipped to use nitrate (Hattenrath-Lehmann et al., 2015; García-Portela et al., 2020); thus, in this study, we disabled the use of nitrate by Dinophysis (Figure 2). In any case, the intracellular regeneration of ammonium within Dinophysis on account of its feeding and then, later the increasing availability of external DON and ammonium would likely repress any significant need to use nitrate. Dinophysis has been found successfully growing in nitrogen-depleted waters (Seeyave et al., 2009). Thus, elevated nutrients have an indirect impact on Dinophysis blooms via the promotion of Mesodinium growth, plus a potential direct effect. This duality of direct and indirect effects of eutrophication on Dinophysis blooms also finds support by experimental studies on D. acuminata and M. rubrum (Hattenrath-Lehmann et al., 2015; Hattenrath-Lehmann and Gobler, 2015).

Anschütz and Flynn (2020) showed, through in silico experiments, that good growth conditions allow blooms of SNCM of such a magnitude that the essential provider of acquired phototrophy (e.g., Teleaulax for Mesodinium) is driven to extinction, and then, the SNCM bloom subsequently also becomes extinct (or the organism is forced to encyst). Lower nutrient systems prevent boom-and-bust dynamics and thus can support lower (but perhaps still actionable) levels of Dinophysis. The implication is that decreasing sources of eutrophication need not necessarily remove the potential for harmful blooms of Dinophysis.

With the exception of the simulations at different nutrient loads shown in Figure 4, the simulations were run without grazing on Dinophysis to focus on the effect of the interactions of Teleaulax, Mesodinium, and Dinophysis on their respective bloom dynamics. In the field, Dinophysis is subjected to grazing by protists (e.g., Fragilidium; Park and Kim, 2010) and metazoan zooplankton (e.g., copepods; Kozlowsky-Suzuki et al., 2006), which keeps their cell numbers lower. The activity of such grazers can be a major factor leading to the decline of Dinophysis blooms, together with changes in hydrodynamic characteristics of the water column (Escalera et al., 2010; Velo-Suárez et al., 2010). In our TMD-complex simulator, applying a generic encounter-based grazing rate to Dinophysis (Figure 4) at different nutrient loads did not affect the maximum biomass of Synechococcus, Teleaulax, or Mesodinium and only prolonged the bloom phase of Mesodinium for a few days under eutrophic conditions. The timing and duration of blooms of Dinophysis did not change but the maximum cell numbers were lower. While grazing has clear potential to constrain Dinophysis numbers, ultimately, the absence of the Mesodinium required to furnish plastids for acquired phototrophy prevents Dinophysis growth. Our study suggests that the mixoplanktonic nutritional mode is a key to explain the seasonal lag of these blooms reported in the literature (Álvarez-Salgado et al., 2011).

Irradiance

Our simulations indicate that light is not a major factor affecting bloom duration, but it did affect the peak height (Figure 3) and thence, the overall dynamics. The model by Moeller et al. (2016), describing the interaction between the cryptophyte Geminigera cryophila (as the plastid donor) and M. rubrum, found that the need for acquired phototrophy affected the coexistence of the organism with its G. cryophila prey, such that increased irradiance shifted a stable coexistence to boom and bust cycles. Our simulations, portraying the entire TMD-complex, showed that such an interaction extended further to Dinophysis. However, while the maximum biomasses of Teleaulax, Mesodinium, and Dinophysis increased significantly with elevated irradiance (Figure 3), the boom and bust dynamic did not repeat in cycles (Figure 6). For Teleaulax, increased irradiance shortened the bloom duration (Figure 3) by accelerating their growth. Without a recurring reseeding (in nature, e.g., by advection), Teleaulax and Mesodinium were grazed to extinction, and Dinophysis, thus, could not take advantage of the increased nutrients. As noted above, the appearance of other predators grazing these organisms would also prevent a cycle from developing.

Higher irradiance has been found to decrease the time that Dinophysis can retain functional kleptochloroplasts in the absence of reacquisitions (Hansen et al., 2013). Irradiance may impact pSNCM growth by triggering photoacclimation in the acquired photosynthetic apparatus while decreasing their longevity due to photodamage. A higher level of photodamage may raise the demand for their specific prey to resupply plastids. Further developments of the TMD-complex simulator could include the implementation of an irradiance dependence of the parameter crit_IR that controls the longevity of acquired phototrophy (Supplementary Data Sheet 1 Table S3), and photoacclimation could be simulated by changing the chlorophyll to carbon ratio (Chl:C). Photoacclimation is also known to be important to pSNCM growth in variable light (Moeller et al., 2016). Doing justice to this process and the decline of acquired phototrophy linked to both light and nutrient status requires models that better describe food quality that impacts on mixoplankton activity (e.g., Mitra and Flynn, 2010; Mitra et al., 2014; Lundgren et al., 2016; Lin et al., 2018). This requires the use of full variable stoichiometric (C:N:P:Chl) coarse-grain system biology-based ecosystem simulators (e.g., Flynn and Mitra, 2009; Leles et al., 2021).

Temperature

Temperature has a significant effect on the frequency of and extent of Dinophysis blooms (Escalera et al., 2006; Díaz et al., 2016; Pitcher et al., 2017). Temperature impacts are both abiotic and biotic. Warmer waters lead to a more stratified and stable water column that can favor Dinophysis. Such conditions are expected to become more frequent with climate change (Sallée et al., 2021). Temperature also increases metabolic rates, enacted in our TMD-complex food-web model by reference to a Q₁₀ factor (see the “Materials and Methods” section). Temperature affected the time frame of the simulated blooms of Teleaulax, Mesodinium, and Dinophysis and the timing consequentially impacted the likelihood and size of Dinophysis blooms (Figure 7). However, the elevated temperature did not necessarily result in the largest blooms of Dinophysis as it accelerated the growth phase of Mesodinium such that Dinophysis had less exposure time to its specific prey (Figures 3, 7). The optimal matching of Dinophysis appearance and growth with that of Mesodinium appeared to be a critical factor for Dinophysis. While hydrodynamic events of stratification, advection, and upwelling are crucial factors for Dinophysis bloom events (Escalera et al., 2010; Raine et al., 2010; Ajani et al., 2016; Díaz et al., 2016), the components of the TMD-complex need to be present in the appropriate magnitude at the right time for Dinophysis blooms to develop and produce high cell abundances (González-Gil et al., 2010; Sjöqvist and Lindholm, 2011). Mesodinium species and many toxic species of Dinophysis occur globally from polar to tropical regions (Figure 1; Gustafson et al., 2000; Reguera et al., 2012) and are, thus, exposed and adapted to different climates. There are some differences in the distribution of species of Dinophysis. The species D. acuminata is ubiquitous; D. norvegica is more common in temperate to boreal regions while D. caudata bloom in tropical and subtropical regions (Figure 1; Reguera et al., 2012). It is unknown what effect climate change will have on the species distribution of Dinophysis. Furthermore, climate change will have different effects in these regions related to currents, upwelling and downwelling, stratification, and ocean acidification. Integrating the TMD-complex system dynamics model presented in this study, within a climate model could help to address these questions and lead to interesting new results in the research of these organisms.

Mixed Layer Depth and Advection

The MLD affects the way that light, nutrients, organism abundance, and temperature impact plankton ecology. The model simulated average biomass levels, while in reality, both Mesodinium and Dinophysis form aggregations that will affect encounter rates. Dinophysis often aggregate in very high abundances in thin layers (Campbell et al., 2010; González-Gil et al., 2010; Sjöqvist and Lindholm, 2011; Hattenrath-Lehmann et al., 2013; Díaz et al., 2016; Berdalet et al., 2017; Baldrich et al., 2021) that may even resemble monocultures (Farrell et al., 2012). This strategy increases encounter rates with their much faster prey enabling the Dinophysis to launch its ambush-style attacks on Mesodinium (Jiang et al., 2018). Turbulence damages ciliates (Smayda, 2002), and M. rubrum is documented to actively avoid turbulence (Cloern et al., 1994). In fact, the vertical stability of the water column is essential to enable M. rubrum to reach the highest cell abundances (Sanders, 1995; Herfort et al., 2011). Thus, stable water conditions may coincidentally favor the growth of their Dinophysis predator.

Simulations of conditions in a narrow band at depth (Figure 8) show how the formation of hydrodynamic conditions conducive to the formation of Dinophysis layers has the potential to mimic the consequences of bulk-water eutrophication in raising in situ predator-prey encounters. Mixoplanktonic activity provides a route to compensate for a lack of photosynthate in especially eutrophic conditions where the high biomass simultaneously shades light and provides high abundances of potential prey. Both the pSNCM species (Mesodinium and Dinophysis) encountered much more prey in the simulated 1-m thickness layer. For Mesodinium and especially Dinophysis, this led to a much higher biomass yield.

Diarrhetic shellfish poisoning events are often associated with the physical transport of Dinophysis populations by wind-driven advection and coastal counter-currents (Escalera et al., 2010; Raine et al., 2010; Reguera et al., 2012); therefore, the forecasting power of physical models for HAB events are linked to hydrodynamics (Raine et al., 2010; Velo-Suárez et al., 2010; Savvidis et al., 2011; Cusack et al., 2016; Moita et al., 2016; Pinto et al., 2016). In addition to horizontal currents, blooms have also been related to upwellings potentially bringing populations of Dinophysis from shelf bottom waters to the surface water and into (continued or enhanced) contact with their M. rubrum prey (Velo-Suárez et al., 2014). These flows also seed the surface water with different organisms in different proportions. These are conditions that influence factors, which our model suggests, can promote the development of Dinophysis blooms. Impacts of each member of the TMD-complex (and of course additional interactions on those members, such as predation and competition) all play parts in the story. The increase in stratification and MLD in the summertime due to climate change (Sallée et al., 2021) will likely also affect the bloom potential of Dinophysis. A changing ocean may challenge the reliability of existing forecasting models that are based on time-series data of hydrodynamic events. A solution could be mechanistic models that combine hydrodynamics and the biology of the TMD-complex components and allied species.

Conclusion

While our study shows the importance of prey availability for the growth of Dinophysis, it also shows that the timing of prey availability has a significant effect on the spatiotemporal dynamics of the TMD-complex. Simulations of the mixoplanktonic activity of these mixoplankton, configured in this study against data for D. acuminata and its prey M. rubrum, both with their special needs for acquired phototrophy, play an important role in the dynamics of bloom formation. Both high M. rubrum to D. acuminata ratios and the growth phase in which D. acuminata encounters M. rubrum are expected to affect cell abundances of D. acuminata in subsequent blooms. Even though they are not simulated in this study, it is very likely that these effects extend to other species from the genera Mesodinium and Dinophysis that share a similar relationship within cryptophyte-MD-complex. Monitoring programs for Dinophysis could, therefore, benefit from including Mesodinium in their routine observations; this has the advantage that Mesodinium appears in much larger numbers making it easier to detect than the cryptic Dinophysis. The disparity between monitoring and research efforts expended on just the terminal end of the TMD-complex (i.e., on Dinophysis) is notable. The other critical components, namely, Mesodinium and Teleaulax, receive far less attention (Figure 1) and yet their success and timing set the scene for blooms of the toxigenic species.

The dynamics shown by the simulations will perhaps act to prompt more attention in monitoring these other species. For aquaculture, it is not only important to know when the blooms will occur but also their duration. While the main reasons for the decline in Dinophysis blooms are often linked to dispersion (Velo-Suárez et al., 2010), grazing (Kozlowsky-Suzuki et al., 2006), and infections with parasites like Amoebophrya sp. (Velo-Suárez et al., 2014), a collapse in Mesodinium will inevitably lead to a failure in the Dinophysis a few weeks later just like an appearance of a Mesodinium bloom could indicate the potential for development of a harmful bloom of Dinophysis.

The in silico food-web simulator we used in this study may help to test sensitive parameters of different species and theoretically explore their behavior and inspire future experiments. The TMD-complex provides an interesting and closely coupled example of successional interactions, and our model, thus, provides a tool to explore reasons for differences in bloom dynamics between species. As a nitrogen-based model, this simulator is also suitable for implementation in larger nitrogen-based ecosystem models and, therefore, exploration of TMD-complex dynamics and impacts in a larger context. The model presented in this study cannot replace a hydrodynamic model and on its own is not sufficient as a predicting tool for Dinophysis blooms occurring in the field. It is, however, the first biological model featuring Teleaulax, Mesodinium, and Dinophysis. This model could deliver useful results in combination with a hydrodynamic model that can account for stratification and Dinophysis input via currents from the outside (from offshore or up the coastline).

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

AM and KJF conceived the original idea. A-AA ran the simulations under the supervision of KJF and AM. All authors contributed to the construction of the models, analysis of the results, and writing of the manuscript.

Funding

This research was undertaken within the MixITiN project. This project has received funding from the European Union’s Horizon 2020 Research and Innovation Program under the Marie Skłodowska-Curie Grant Agreement No. 766327. AM was supported by MixITiN and Sêr Cymru II WEFO ERDF Programme MixoHUB 82372.

Conflict of Interest

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.799358/full#supplementary-material

Footnotes

1. ^ www.Powersim.com

References