Introduction

Many environmental and pathogenic bacteria exhibit colonial forms of life, either as surface biofilms [1, 2] or non-attached aggregates [3–5]. In most situations, these colonies are subjected to hydrodynamic stresses, which may lead to flow-induced fragmentation and aggregation [6]. These physical processes, combined with biological growth, determine the fate of aggregates and the microbial lifestyle. Nonetheless, despite their abundance and fundamental importance, a generic framework to study the effects of fluid flow on bacterial colonies, as well as their fragmentation and aggregation across scales, is still in its infancy [7–10]. Here, we focus on colonies of bloom-forming cyanobacteria and ask the question: how flow-induced aggregation and fragmentation regulate dynamic changes in colony size?

Toxic cyanobacterial blooms are a major nuisance for freshwater ecosystems worldwide [11]. These blooms often lead to increased turbidity, cause harm to other aquatic organisms, and deteriorate water quality. Recent evidence indicates that the severity and frequency of cyanobacterial blooms are increasing due to eutrophication and climate change [12–15]. Consequently, water management authorities show growing interest in techniques to prevent, monitor, and control these toxic blooms [16]. A particular effort is being taken to combat blooms of Microcystis, a non-motile colony-forming freshwater cyanobacterial genus that dominates many aquatic ecosystems worldwide [17, 18]. Fluid dynamics is at the heart of Microcystis bloom formation, but also on some of the methods to control them. A large colony size and positive buoyancy give Microcystis a higher average light dose in comparison to non-buoyant algae [19, 20]. Control methods such as deep artificial mixing can disrupt this competitive advantage by increasing the turbulence level in lakes [21]. Yet, the mechanisms of Microcystis colony formation, and in particular, the impact of fluid flow on the formation of colonies are poorly understood.

Microcystis cells are held together by a layer of extracellular polymeric substances (EPS), facilitating the formation of large colonies of up to a millimeter in diameter with diverse morphological arrangements [22–24]. Two main mechanisms for colony formation have been proposed [17]: (i) cell division, i.e., clonal expansion of cells which remain attached by their encapsulating EPS layer after division, and (ii) aggregation, i.e., initially independent cells/colonies collide and adhere to each other. On the one hand, sequencing of the genetic diversity within Microcystis colonies supports the hypothesis that colony formation under natural conditions is primarily driven by cell division [25]. On the other hand, cell aggregation can occur on a shorter time scale and may offer improved protection against high grazing pressure [26]. Aggregation becomes significant when the collision frequency is high enough for aggregation rates to overcome cell division rate [6], and the stickiness of colonies is sufficiently large (controlled by the EPS composition, pH, and ionic strength of water or medium [27–29]). In addition to aggregation, fluid flow also influences colony formation through fragmentation, potentially serving as a mechanism for colony reproduction [30] to overcome the growth limitations of large colonies [31]. Beyond modulating colony size through aggregation and fragmentation, fluid flow has also been proposed as a driving factor in changes in colony morphology and growth rate [32–34].

Fig. 1

Flow-induced alterations in colony size and morphology impact the vertical migration and grazing resistance of colonies, thereby influencing the success of the Microcystis population. This bears significance for forecasting models which have often assumed a fixed colony size [35–37]. Information on colony size is also important for bloom control via artificial mixing systems, which rely on estimates of colony flotation velocity, determined by both colony shape and size [21, 38]. Moreover, it remains uncertain whether mixing systems, such as bubble plumes, can enhance their efficiency by fragmenting colonies or whether artificial mixing may inadvertently promote aggregation. Colony fragmentation has been observed in grid-stirred [39] and propeller mixed tanks [32] at energy dissipation rates comparable to bubble plumes [40]. However, these studies reported averaged dissipation rates, while local values near the fast moving surfaces can be orders of magnitude higher [41]. Furthermore, the size distribution of the colony fragments has not been characterized.

In the present study, we investigate to what extent Microcystis colony fragmentation and/or aggregation occurs at levels of turbulence observed for naturally or artificially mixed lakes. To this end, colonies of a laboratory culture of Microcystis were subjected to various intensities of shear flow (Fig. 1). We identified the mechanism of colony fragmentation using direct measurements of dynamic changes in colony size distribution under controlled laboratory conditions and analyze the characteristics of colonies formed by aggregation and cell division. Furthermore, the fragmentation behavior of the laboratory culture was compared with field samples of Microcystis spp. We also developed a numerical model based on a two-category discrete population to simulate changes in the size of Microcystis colonies via aggregation and fragmentation. Finally, we provide a phase map based on flow characteristics and cyanobacterial abundance to highlight different regimes in which fluid flow plays a significant role in cyanobacterial colony formation.

Fragmentation of Large Microcystis Colonies Occurs through Erosion

A culture of Microcystis strain V163 was filtered to collect large colonies, and colonies on the filter were subsequently re-dispersed in fresh BG-11 medium up to the desired total biovolume fraction ϕ (Materials and Methods). The flow was generated by a cone-and-plate shear operating in an inertial regime, and it was characterized by the average energy dissipation rate , hereafter called dissipation rate (Materials and Methods). The size distribution of the colonies was measured over time by microscopic imaging of the suspension through the transparent lower surface (Fig. 1). The colony size is described here by its relative diameter l = d/d₁, where d is the Feret’s diameter (maximum distance between two boundary points of the colony) and d₁ is the single-cell diameter. The number of cells in a colony, N_c , was calculated from its relative diameter via the relation N_c = l^f, where f = 2.09 is the average fractal dimension measured for colonies of Microcystis strain V163 (Materials and Methods). The size distribution was computed as the normalized biovolume distribution n(l), i.e., the number of cells in colonies in size bin {l; l + Δl}, divided by the total number of cells in the suspension.

The results show that the initial size distribution had a bimodal shape, with a subpopulation of small colonies peaking around l ~ 1 – 2 (single cells, dimers, and trimers) and a subpopulation of large colonies, peaking around l ~ 20 – 40 (Fig. 2A). The cut-off size between the two subpopulations is defined here as the minimum of the bimodal distribution. Some small colonies (l ~ 1 – 2), up to 20% of the total biovolume, were always present in the initial distribution even after filtration, due to pre-stirring required to keep the suspension homogeneous. However, medium-sized colonies (l = 4 – 8) accounted for a negligible fraction of the distribution.

When subjected to an intense dissipation rate (), the large colonies of Microcystis strain V163 displayed considerable fragmentation, with their median diameter (50 % percentile of the biovolume distribution of the large colonies) decaying to about half of the initial value after a few hours (Fig. 2B). In contrast, the median diameter of small colonies showed little change in size and remained below l ≤ 2 during the entire experiment. The size distribution after 1 hour indicated that most colonies larger than l ≥ 30 were fragmented, leading to a rise in the subpopulation of small colonies (Fig. 2C). Notably, the bimodal shape of the size distribution was preserved, with very few medium-sized colonies. This suggests that aggregation of single cells into larger colonies is negligible under this intense shear flow condition, and the main effect of the shear flow on the size distribution is the decrease in size of the large colonies concomitant with an increase in the fraction of small colonies. After a few hours of shear, about 80 % of the total biovolume was composed of small colonies (Fig. 2D). This time scale was much shorter than the cell division time, such that changes in the colony size distribution by division could be neglected within the experimental time [20].

Direct measurement of the fragment distribution function g(l,l’) (i.e., the frequency of fragments with size l arising from the fragmentation of a colony of size l’) is an experimental challenge [42–44]. Individual fragmentation events could not be visualized as the residence time of colonies inside the field of view is too short. Instead, a suitable model can be identified from the rate of change in the biovolume distribution (Fig. 2E). A negative rate was observed for large colonies in the range l ≥ 25, indicating a loss in that size range by fragmentation. In contrast, the newly created fragments appeared as positive peaks, with a distinct peak close to the single-cell size. This result is compatible with an erosion mechanism, in which small fragments with one or a few cells are detached from the outer layer of the colony. For an ideal binary erosion, a fragmentation event of a colony with i cells results in one cell being eroded from the surface of the colony with a complementary fragment of i – 1 cells remaining. Small deviations from an ideal binary erosion are observed in Figure 2E, with most eroded fragments lying in the range 1 ≤ l ≤ 2. Nevertheless, ideal erosion provides a good approximation for the results. Accurate identification of complementary fragments from the rate of change in the biovolume distribution was more difficult, because the broad initial distribution of our suspension leaded to a broad distribution of complementary fragments. In this case, the rate of change in the biovolume distribution shifted from a net loss of eroded colonies at larger sizes (l ≥ 25) to a net gain of complementary fragments at intermediate sizes (5 ≤ l ≤ 20).

Assuming an ideal binary erosion mechanism, the fragmentation frequency κ (l) was computed from the temporal dependence of the biovolume distribution (Materials and Methods). Here, κ(l) is defined as the average number of fragmentation events (loss of a single cell) that one colony of size l suffers per unit of time. The fragmentation frequency of Microcystis strain V163 was calculated as a function of colony size for three values of dissipation rate (Fig. 2F). The results show that the fragmentation frequency increases with colony size. Furthermore, an increase in dissipation rate sharply intensifies the erosion of colonies, with a two order of magnitude increase in κ(l) when varies from 1.1 to 9.6 m²/s³.

Aggregation-Formed Colonies are Less Resistant than Division-Formed Colonies

In addition to inducing fragmentation, a turbulent flow may also increase the collision frequency of colonies. If the collisions are adhesive and the suspension is concentrated enough, then new colonies will be formed by aggregation. We subjected a single-cell suspension to the cone-and-plate shear flow to investigate whether aggregation can be a dominant colony-forming mechanism compared to cell division. Cultures of the Microcystis strain V163 were filtered to obtain a single-cell suspension, and the total biovolume fraction was adjusted with fresh BG-11 medium (Materials and Methods). A typical initial bio-volume fraction distribution is presented in Figure 3A. In addition to single cells, some small colonies up to l = 3 were present in the filtered suspension due to imperfect filtering and/or aggregation during the pre-stirring. The size distribution was unimodal (single peak), in contrast with the bimodal distribution (double peak) in the suspension of division-formed colonies (Fig. 2A). After about 1 hour at a moderate dissipation rate () and a total biovolume fraction of ϕ = 10⁻⁴, the single-cell suspension showed aggregation of cells into small colonies, with a steady-state median diameter of l = 1.8 (Figs. 3B and 3C). Note that cell division can be neglected as the measurement time was much shorter than the division time. An increase in the collision rates due to a single-cell suspension led to an increase in the median diameter, although the increase was minor. At a high total biovolume fraction (ϕ = 5 • 10⁻⁴), a steady-state median diameter of l = 2.4 was reached, with a negligible fraction of colonies above l ≥ 5 (SI Appendix, Fig. S4C). Therefore, the aggregation of single cells is likely limited by weak bonds between cells. At the steady-state regime of median diameter, the newly-formed aggregates with a large diameter would be quickly fragmented back into smaller parts.

Fig. 2

Fig. 3

A suspension of division-formed colonies under the same moderate dissipation rate () did not suffer noticeable fragmentation (Fig. 3D). The culture was grown in quiescent medium, without shaking, thus the colonies are assumed to have formed solely by cell division. The erosion of division-formed colonies was negligible, as the fraction of small colonies remained below 20% (SI Appendix, Fig. S2B). Furthermore, the nearly constant median diameter of the large colonies (observed also for other total biovolume fraction values, SI Appendix, Fig. S2) shows that the division-formed colonies did not aggregate at this moderate dissipation rate.

The results above indicate that bonds formed during flow-induced aggregation (whether between two single cells or two division-formed colonies) are weak and cannot withstand intense dissipation rates. In contrast, the bonds between cells within a single division-formed colony are significantly more resistant. This difference in mechanical strength of the cellular binding in aggregation-formed versus division-formed colonies can be attributed to the structure of the EPS layer that surrounds the colonies. Aggregated colonies would have interstices between cells that are not filled with EPS. In contrast, cell division allows enough time for the gaps to be filled with secreted EPS, effectively increasing the bond strength between cells (Figs. 3B and 3D). However, we could not confirm this hypothesis visually, as the EPS layer in Microcystis strain V163 is too thin (< 1 μm) for a clear identification using optical microscopy.

Dynamical Changes in Colony Size Modeled by a Two-Category Distribution

Our experimental results for colonies of Microcystis strain V163 (Fig. 3) suggest the presence of two categories of colonies: aggregated colonies (denoted as C₁) and division-formed colonies (C₂). The total biovolume distribution n(l) is given by the sum of both distributions n₁(l) + n₂(l), normalized such that . Within times shorter than the cell division time, the biovolume distribution of each colony category is described by a population model [27]:

where A_i(l, t) and F_i(l, t) describe the aggregation and fragmentation rates as a function of colony size for each of the two categories, and T_i (l, t) is a transfer rate of biovolume between the two categories. The turbulent shear-induced collisions between colonies shift the distribution towards larger aggregates, and the aggregation rate is given by

where l_o = (l^f – l’^f)^1/f. The collision frequency is described by the kernel β(l, l’) and for a turbulent shear flow, we have , where ν is kinematic viscosity of the fluid [45]. The ratio of aggregative collisions over total collisions is given by the stickiness parameter α_i. For the category C₁ (aggregated colonies), the stickiness is positive, α₁ > 0. In contrast, we neglect the aggregation of C₂ colonies (collisions between C₂−C₂ and C₁−C₂), since the bonds formed by aggregation would be weaker and quickly rupture again, such that α₂ = 0. This assumption is corroborated by the absence of aggregation between large colonies in our experiments (Fig. 3D).

The colony formation will be limited by fragmentation if the colony exceeds a critical size , where S_i and q_i are constitutive parameters that depend on the aggregate structure and bond strength, with a lower bound at [46, 47]. We assume that aggregation-formed colonies have a small critical size, as the overlap between the EPS layer of each cell is narrow, leading to weak bonds. Meanwhile, division-formed colonies have an EPS layer that fills the pores between cells, thus increasing the resistance to hydrodynamic stress. The fragmentation rate is given by a first-order kinetic equation [46]:

where κ_i (l) is the fragmentation frequency, which is a function of the probability distribution function of the hydrodynamic stress. Assuming a Gaussian distribution, the fragmentation frequency is [46]:

where . Due to the scarcity of experimental data in the literature, simplified models are usually assumed for the fragment distribution function g_i(l,l’), such as erosion, equal fragments, or a uniform distribution [48]. Assuming a binary fragmentation (i.e., two fragments per fragmentation event), the fragment size distribution can be written in the form g_i(l,l’) = δ(l – l_a) + δ(l – l_b), where l_a and l_b are the sizes of fragments. Based on the results from division-formed colonies under intense dissipation rate (Fig. 2E), an ideal erosion mechanism is adopted for category C₂ colonies, such that . In contrast, the experimental results for the small colonies (Fig. 3C) show a wider size distribution of eroded fragments which can not be captured by an ideal erosion mechanism, so an equal fragments mechanism is adopted for category C₁ colonies, i.e., . The size distributions of C₁ and C₂ colonies are coupled through a transfer rate of biovolume T_i (l, t). We propose that the cells eroded from C₂ colonies are transferred to category C₁, such that:

while the complementary fragments of C₂ colonies remain in that category.

We solved the balance equation (1) using a zero-order discrete population balance, or sectional method, in which the size distribution is described by uniformly sized bins where the size distribution is evaluated at a representative size centered on the bin [48, 49]. The simulations were initialized with the size distributions at t = 0 obtained in our experiments. The two pairs of parameters for critical size ({q₁,S₁} and {q₂,S₂}) and the stickiness parameter (a₁) were used to fit the model to the experimental data (Materials and Methods).

For an initial suspension composed mainly of large colonies at an intense dissipation rate (), ), the model prediction of a decrease in the median size of C₂ colonies up to a steady-state diameter is in good agreement with the experimental findings (Fig. 2B). Moreover, the erosion mechanism accurately predicts the biovolume transfer from C₂ to C₁ colonies (Fig. 2D). The use of the erosion mechanism for the fragment distribution function is further corroborated by previous studies simulating bacterial aggregates [50]. For C₁ colonies at this intense dissipation rate, the model predicts that their median diameter remains near l* = 1 (i.e., the lower bound for the critical size). However, the experimental results show a slightly wider size distribution of both the small and large colonies than predicted by an ideal erosion mechanism (Fig. 2C). Nevertheless, the global behavior of the size distribution is well described by the two-category model. Furthermore, a very good agreement is also observed between the measured fragmentation frequency and the mo del predictions given by Eq.(4) (Fig. 2F).

To complement the mathematical model and to further investigate the mechanical origins of the resistance of division-formed colonies to intense shear, we performed rheology tests on the cyanobacterial colonies (Materials and Methods). Simple and oscillatory shear measurements revealed elasto-viscoplastic mechanics [51, 52] of the concentrated aggregate with a yield stress of τ_y = 4.3 ± 0.3 Pa. The low volume fraction of cells in aggregate (~ 12 %) indicates that the mechanical properties are mainly dictated by the EPS layer. Equating the average hydrodynamic shear stress to the yield stress results in a critical colony size of (Eq. 10 in Materials and Methods). The eroded fragments observed in the fragmentation experiments (Fig. 2E) are mostly smaller than this critical size. This suggests that large colonies will suffer erosion when the local hydrodynamic stress exceeds the yield stress of the EPS layer. Although the fracture behavior of colonies in suspension is expected to be different from a concentrated aggregate under shear rheology, the latter provides a reasonable and convenient estimate of the material properties of the colonies. Related to the critical size of the colony and their mechanical properties is the exponent q₂ (see equation 4). The best fit of the results in Fig. 2 provided q₂ = 4.1. This value lies in the upper limit of previous results for colloidal aggregates [42, 46]. Note that for colloids, q is related to the stress concentration around cracks in the fractal structure of the aggregate [47], and low values of q are associated with brittle deformation, in which the stress concentration grows quickly with aggregate size. In contrast, cyanobacterial colonies have a soft plastic behavior that is less sensitive to large cracks, thus justifying the large value of q.

We further tested the model with the aggregation experiments. A suspension of C₁ colonies with an initial distribution close to the single-cell size (Fig. 3A) was simulated. The model agrees well with the experimental results for the median colony size, accurately predicting a steady-state value when there is a balance between collisions and fragmentation (Fig. 3B). The best fit of the results provided a small value for the stickiness parameter (α₁ = 0.023). The aggregation rate is linearly dependent on the stickiness and the total biovolume fraction (Eq. 2), thus an increase in these parameters (α₁ and ϕ) allows the suspension to reach the steady-state size faster. However, the steady-state size has low sensitivity to these parameters (α₁ and ϕ), as the fragmentation frequency has a higher-order dependence with colony size (Fig. 2F). The equal fragments hypothesis leads to a single-peaked size distribution, whose steady-state median increases with the total biovolume fraction. If an ideal erosion mechanism was assumed for the C₁ colonies, the peak would be fixed around the single-cell size (SI Appendix, Fig. S6). The experimental results for the C₁ colonies are in better agreement with the equal fragments hypothesis (Fig. 3C).

Note that the ideal erosion and equal fragments mechanisms are idealized scenarios, and the real distribution is expected to be more dispersed. Furthermore, a fractal description of small colonies is limited, as a continuum function between biovolume and diameter of colonies may not hold. Moreover, the turbulent model used for fragmentation (Eq. 4) poses two limitations: (i) the theory was developed for homogeneous Gaussian turbulence [46], which differs considerably from the boundary-layer turbulence generated in a cone-and-plate shear [53], and (ii) the colonies were larger than the Kolmogorov scale for the intense dissipation rate used in our fragmentation measurements (η = 20 μm for a dissipation rate of , SI Appendix), meaning that the colony size is comparable to the smallest eddies in the flow. Despite these limitations, the proposed two-category model provides a good approximation of the fragmentation process observed in our experiments. Simulations that take into account spatially inhomogeneous and non-Gaussian turbulence are demanding and lack an explicit analytic expression for κ(l) that can be applied to a population model [54].

Fig. 4

Field Samples of Microcystis are more Resistant to Fragmentation due to Thicker EPS Layer

Samples of Microcystis spp. were collected from Lake Gaasperplas - the Netherlands (Materials and Methods). A major difference in fragmentation behavior is observed between the laboratory cultures of Microcystis strain V163 and the field samples of Microcystis spp. While the laboratory cultures suffered intense fragmentation at a dissipation rate of , the colonies collected from the field showed negligible fragmentation at this dissipation rate (Fig. 4A). The size distribution of field colonies also had a bimodal shape both initially and after 1 hour of flow (Figures 4C and E). The small colonies displayed an increase in size due to aggregation up to 2.5x times the initial median size (Fig. 4D), although the fraction of small colonies remained below 11% even after a few hours of shear (Fig. 4F). The large colonies displayed only a small increase in size due to aggregation and quickly reached a steady-state size.

The higher resistance of field colonies to shear can be attributed to the thick EPS layer that surrounds the colony, which could be seen under brightfield microscopy using a negative staining (Fig. 4B, Materials and Methods). In contrast, this layer was barely visible in colonies of laboratory cultures of Microcystis strain V163. Nevertheless, the limited amount of aggregation observed for large field colonies indicates that aggregation-formed colonies have a smaller resistance to shear than division-formed colonies, similar to the results observed for laboratory colonies (Fig. 3). Differences in the morphology of colonies from field samples and laboratory cultures are frequently observed in the literature, with the latter often presenting only single cells or small colonies [55]. Colony formation may be influenced by various factors [17], where especially high calcium and low nitrogen concentrations are identified as factors promoting high EPS production [56–58]. This is consistent with the calcium and nitrogen concentrations measured in Lake Gaasperplas (71 mg/L Ca²⁺ and 0.04 mg/L N) versus BG-11 medium (9.8 mg/L Ca²⁺ and 246 mg/L N) (Materials and Methods).

Aggregation, Fragmentation, or Cell Division? A Guideline

Under which conditions are aggregation, fragmentation or cell division the dominant processes in colony formation? In a quiescent fluid, cyanobacterial colonies are expected to grow solely due to cell division with a specific growth rate , which is a function of various environmental factors such as light intensity and nutrient availability [17, 20]. Let be the average colony size of a suspension. Assuming that the cells remain attached after division, the rate of increase in the average colony size will be

The increase in colony size will progress until a maximum size l_m is reached. The upper limit in colony size is still a subject of intense discussion, with species variability, reduced division rate at large sizes, and dissolution of the EPS layer being important factors [17, 31]. If the suspension is concentrated enough, collisions between colonies will also increase the average colony size due to aggregation. Let us consider initially a monodisperse suspension with a single size . At short times, the rate of increase in the average colony size due to aggregation will be (SI Appendix)

Assuming f = 2 for simplicity (which is similar to the measured value for Microcystis strain V163), we can use Eqs. (6) and (7) to define a dimensionless number Ag, quantifying the ratio between the rates of increase in the average colony size due to aggregation and cell division, given by

The dimensionless ratio A_g reveals the dominant mechanism for colony formation (cell division is dominant for A_g ≪ 1, while aggregation is dominant for A_g ≫ 1) and depends on the total biovolume fraction ϕ and dissipation rate , as shown by the phase diagram in Figure 5. For a lake with a total biovolume fraction corresponding to an early bloom (ϕ = 3 · 10⁻⁷, according to WHO alert level 1 [59]), a low dissipation rate typical of wind-induced turbulence () [60], a specific growth rate of μ ~ 0.008 h⁻¹ [20], a kinematic viscosity of ν = 10⁻⁶m²s⁻¹ and a stickiness of α ~ 0.02, we estimate a ratio A_g ~ 10⁻³. Under these conditions, the colony size increases due to cell division while the role of aggregation is negligible (region I in Fig. 5). In contrast, during a dense bloom event, the scum layer can reach a very high biovolume fraction of ϕ ~ 10⁻² [61, 62], which leads to A_g ~ 10₂ (keeping all remaining parameters constant). Therefore, colony formation in dense scum layers is dominated by colony aggregation (region II in Fig. 5). The aggregation process can form larger colonies which will have an enhanced buoyancy, thus assisting in the persistence of the scum layer.

Fig. 5

Fragmentation limits the maximum colony size to a critical value l*. Our experimental results (Fig. 5) indicate that aggregated colonies have a lower critical size than division-formed colonies for the same turbulence intensity. Conversely, a higher critical dissipation rate is needed for division-formed colonies () than for aggregated colonies () to achieve the same critical size l* = l_m. This creates an interesting regime at high biovolume fractions and dissipation rates in the range , in which only division-formed colonies can survive the turbulent flow, even though aggregation is faster than division (A_g ≫ 1) (region III in Fig. 5). Typical values of critical dissipation rate for our experiments are and , assuming a maximum colony size of l_m = 50. Under natural turbulent conditions (e.g. surface wind-mixing), the dissipation rate is normally below [60], which is insufficient to induce colony fragmentation. Artificial mixing systems can achieve levels of turbulence that are many orders of magnitude higher. In the core of bubble plumes, dissipation rates of are reached [40]. This value was sufficient to prevent aggregation of large colonies in our experiments but was not high enough to induce fragmentation of division-formed colonies. Moreover, colonies remain only a short time (~ 1 min) in the bubble plume core, after which they are dispersed to weaker turbulent regions (dissipation rates of just a few meters away from the bubble plume core). The intense hydrodynamic stress required to fragment division-formed colonies (region IV in Fig. 5) can only be achieved by mixers with fast-moving surfaces, such as the cone-and-plate setup used in our experiments () or propellers used in stirred tanks () [41, 63]. These mixing devices can provide intense dissipation rates in small fluid volumes, but their implementation across entire lakes would be excessively energy-consuming.

Conclusions

Despite the widespread presence of bacterial colonies across species and environmental habitats, the influence of fluid flows on colony formation remains a fertile ground for continued exploration. A significant challenge is the lack of a generic framework in which the aggregation and fragmentation of colonies can be studied under various hydrodynamic stresses. We take a major step towards overcoming this challenge by presenting a combination of experiments and theory to elucidate the underlying mechanisms that determine the size distribution of cyanobacterial colonies over a wide range of flow intensities. Our experiments show that division-formed Microcystis colonies can be fragmented under intense dissipation rates , following an erosion mechanism. Previous studies [32, 39] have observed significant fragmentation of Microcystis colonies in field samples at dissipation rates as low as . However, such studies primarily report volume averages of the dissipation rate, overlooking the likelihood of fragmentation occurring near fast-moving solid surfaces (grid stirrer [39] and propeller [32]), where the local hydrodynamic stresses are orders of magnitude higher. Instead, our results indicate that fragmentation of division-formed colonies can not be achieved with wind mixing or aeration systems.

We have demonstrated that flow-induced aggregated colonies of Microcystis form weaker structures than division-formed colonies, likely due to the weak EPS bonds between cells. Furthermore, we showed that colony aggregation is unlikely in lakes under early bloom conditions, as the collision frequency of colonies is too low. Our results, based on the mechanical properties of colonies, are corroborated by previous studies based on the sequencing of individual Microcystis colonies [25, 64]. These studies have identified that the genetic diversity within Microcystis colonies is much smaller than the diversity between Microcystis colonies, thus supporting the hypothesis that colony formation is dominated by cell division rather than by aggregation of colonies. Nonetheless, our results demonstrate that during scum formation in very dense blooms or phytoplankton patches, flow-induced aggregation likely plays a role and may quickly increase colony size. A summary of these outcomes is presented in the phase map illustrated in Figure 5.

From a broader perspective, our findings have relevance beyond predicting and mitigating cyanobacterial blooms. Since fluid flow affects the colonies and aggregates of many species [6–8, 65–67], the experimental methods and theoretical framework presented here, along with the identified mechanisms, could serve as a foundation to better understand the relationship between mechanics and biological functions and traits of colonies under hydrodynamic stresses and their morphodynamics and lifestyle. Other examples include non-natural settings, such as algae or cyanobacteria cultured in bioreactors, where dense cell populations experience intense mixing [68–70]. The results presented here could serve as design guidelines for such applications.

Materials and Methods

A. Sample preparation

Microcystis strain V163 cultured under laboratory conditions, as well as samples collected at Lake Gaasperplas - The Netherlands containing Microcystis spp. colonies were used for the experiments (Fig. 1). Microcystis strain V163 was originally isolated from Lake Volkerak - The Netherlands [71], grows in colonies and was obtained from the algal collection of the Institute for Biodiversity and Ecosystem Dynamics, University of Amsterdam. The Microcystis strain V163 was batch-cultured in Erlenmeyer flasks with standard BG-11 medium. The cultures were kept in an incubation chamber at a temperature of 20 °C. The light intensity was kept at 12 μmol photon m⁻²s⁻¹ with a light:dark cycle of 14h:10h. The cultures were diluted with new BG-11 medium every week and kept at a maximum total biovolume fraction of ϕ = 5 · 10^–4. The cultures were not shaken except after the addition of medium and before filtration. The Microcystis strain V163 has a single-cell diameter of d₁ = 6.8 ± 1.4 μm (SD, N = 189). We categorize partially divided cells also as single cells, as these also behave as single units in a fragmentation/aggregation process. When considering only spherical cells not undergoing division, a diameter of 5.1 ± 0.8 μm (SD, N = 328) was obtained. Colonies of up to 500 μm were present in the culture under the incubation conditions as described before.

Samples of Microcystis spp. were collected from Lake Gaasperplas - the Netherlands (SI Appendix for water monitoring data and location of sampling point) on 14th August 2022. Water from the surface layer was filtered over a plankton net with a 100 μm mesh size. Samples were kept refrigerated at 4 °C, and the experiments were conducted over the next three days. The plankton composition was dominated by Microcystis spp. (89 ± 7% , SD, N = 231), with Microcystis aeruginosa being the main morphospecies, with an average single-cell diameter of d₁ = 7.0± 1.4 μm (SD, N = 100). Colonies of Dolichospermum spp. were present in lower numbers. However, no differentiation of colony species was made during image processing, as it was not feasible to isolate them.

Negative staining was used to verify for the presence of an EPS layer surrounding the colonies of Microcystis strain V163 and of Microcystis spp. in field samples. Nigrosin was added to the surrounding medium, causing the EPS layer to appear bright under brightfield microscopy.

Upon appropriate filtering of the culture, we obtained a suspension with a specific size range. Filter paper with a 25 μm pore size was used to obtain a suspension primarily consisting of single cells (colonies with 2 or 3 cells also present in small numbers). For a suspension mainly comprising large colonies, we filtered the culture with 53 μm nylon mesh. The colonies collected in the mesh were then redispersed in a new BG-11 medium. To estimate the total biovolume fraction of the suspension, 1 ml of the sample underwent a heat treatment to disaggregate colonies [72]. At least 1000 cells were counted in a Neubauer hemocytometer assisted by the image processing software ImageJ. The total biovolume fraction of the suspension was calculated as ϕ = cv₁, where c is the cell concentration and v₁ is the average single-cell volume. Prior to each experiment, the desired biovolume fraction was adjusted using new BG-11 medium for the Microcystis strain V163 culture and with filtered (0.8μm) water from Lake Gaasperplas for field samples of Microcystis spp..

B. Cone-and-plate shear

A cone-and-plate shear flow was selected for our experiments (Fig. 1B and SI). This geometry is commonly used for standard rheological measurements because it provides a homogeneous shear rate under a laminar flow regime. Although less common in inertial regimes, it has been recommended by previous studies to investigate the effects of hydrodynamic stress in tissues [53, 73, 74]. Furthermore, the mixing provided by turbulence helps maintain a homogeneous concentration throughout the suspension.

Our setup consists of an upper conical probe made from transparent acrylic, with an outer diameter of 42 mm and a cone complementary angle of θ = 5.44° (Fig. SI). The conical probe is attached to a steel shaft, and it is driven by a rheometer head Anton Paar DSR 502. The lower surface consists of a 1mm thick glass slide. A cylindrical acrylic ring confines the fluid, preventing spillage of the sample due to strong centrifugal forces. The inner diameter of the ring is 46 mm, which ensures that the nominal shear rate between the upper probe’s outer surface and the ring is the same as between the conical and lower surfaces. A clearance gap of 10 μm was set between the tip of the cone and the lower surface to avoid solid-solid friction.

At low rotational speeds, the flow is laminar with a homogeneous shear rate , where ω is the angular velocity of the conical probe. The flow transitions to a turbulent regime at faster rotations, as the inertial forces become dominant. Since the shear rate is not uniform at the non-linear regime, the nominal shear rate is no longer an appropriate measure for shear intensity. Instead, the average intensity is given by the energy dissipation rate , which measures the rate of energy loss by viscous dissipation per unit fluid mass. Under steady-state conditions, this dissipation rate balances the power injected into the fluid, therefore:

where ρ and V are, respectively, the density and volume of the fluid sample, and T is the torque on the conical probe, measured by the rheometer head. The flow regime can be identified from the energy dissipation rate curve against the angular velocity, with a ~ ω² dependence in the laminar regime and a ~ ω^5/2 dependence in the turbulent regime (SI Appendix, Fig. S7B). All shear experiments conducted in this work were in the latter regime. The average hydrodynamic shear stress can be calculated from the dissipation rate by

The experimental setup provides a configuration where the dissipation rate can be controlled over a large range, while stress values can be measured with high precision. This, combined with optical access and simultaneous microscopy, presents a well-suited tool to study the response of biological samples under fluid flow.

C. Fragmentation and aggregation experiments

The fragmentation and aggregation experiments involved subjecting a suspension of Microcystis colonies to the cone-and-plate shear and using brightfield microscopy to image the colonies under shear. The lower surface was positioned on top of an inverted microscope Nikon TS100 with a 4x objective lens. The center of the field of view was located 15 mm from the rotation axis and with the focal plane set 0.5 mm above the lower surface. Illumination was provided by a Schott KL 1500 HAL light source equipped with a cyan filter (low pass with a cut-off wavelength of 555 nm) to minimize the photosynthetic activity of the sample. Images were captured with a Prime BSI Express sCMOS camera with a resolution of 1.63 μm per pixel. Prior to each experiment, a reference background image was created using distilled water, which was then discarded. A volume of 3.5 ml of the suspension of Microcystis colonies was pipetted in between the shear surfaces. The shear protocol was as follows: pre-shear at for 1 minute to mix the suspension; 12 shear cycles of 5 minutes each followed by 12 cycles of 20 minutes each at the target dissipation rate (imaging cycles, at low dissipation rate of 0.019 m²/s³ for 1 min, were intercalated with each shear cycle). During the imaging cycles, 200 frames were captured at a rate of 10 fps, ensuring that each frame was sampled over a different set of colonies. Examples of the captured frames are shown in Figure 1.

D. Image analysis

The frames underwent pre-processing with background subtraction and thresholding (SI Appendix). Image segmentation was performed with the Scikit Image library for Python [75]. Colony size was determined by their Feret diameter d, which represents the maximum distance between two boundary points of the colony. The number of cells in a colony was estimated from the diameter using the formula N_c = (d/d₁)^f, where f is the average fractal dimension of the colonies. Since not all colonies were within the depth-of-focus, a correction function for the colony count was estimated using a suspension with a known size distribution. The image magnification was optimized for intermediate and large colonies Given the wide range of colonies (from single cells ≈ 7μm to large colonies ≈ 500μm). Colonies smaller than 25μm were measured with lower precision, necessitating a second correction function for accurate counting (SI Appendix, Eq. S1).

To estimate the colony fractal dimension f, a reference sample of colonies had its size distribution characterized in a cylindrical counting chamber, followed by measurement of the total cell count in a hemocytometer. The fractal dimension was fitted such that the estimated cell count was equal to the measured count. For Microcystis strain V163, a value of f = 2.09 was obtained, while for the field samples of Microcystis spp. it was 1.64.

For the direct computation of the fragmentation frequency, we assume a pure fragmentation with an ideal binary erosion, such that the time derivative of the colony size distribution is

where N_p [i](t), κ[i], and l[i] refer to the bin count, fragmentation frequency, and relative diameter at the center of bin i. This formula is only valid for the bins larger than the single cell, as the latter does not suffer fragmentation. The bin count is initialized with the experimental value and evolved in time using an implicit Radau scheme (SI Appendix). The fragmentation frequency of K_b of each bin was obtained by minimizing the squared error between the experimental and numerical bin counts. The script for image analysis is available here. Detected features for all datasets and raw images for the dataset in Fig. 2A-E are available here.

E. Yield stress of cyanobacterial colonies

The mechanical resistance of Microcystis colonies can be estimated by separating the colonies from the medium and measuring the yield stress of the concentrated aggregate. A suspension of Microcystis strain V163 was filtered through a 53 μm mesh, and the large colonies were resuspended in fresh BG-11 medium. The sample was centrifuged for 30 min at 1000 g, and the pellet was separated. A dense gelatinous phase composed of cells and EPS was obtained, with cells accounting for 12% of the total biovolume fraction. The yield stress of the concentrated aggregate was measured with an Anton Paar MCR 301 rheometer. The upper probe was a sand-blasted 50 mm cone with 1 ° angle, while the lower surface was covered with sandpaper to avoid wall slip. Each sample was pre-sheared at 10 s⁻¹ for 2 minutes. Steady-shear and large amplitude oscillatory shear tests were performed (SI Appendix, Fig. S5). From a Herschel-Bulkley fit over the steady shear data, , the dynamical yield stress τ_y = 4.3 ± 0.3 Pa (SD from best fit) was estimated. Due to the limited volume of samples available, only one repetition was possible for steady shear and oscillatory shear tests each.

F. Simulations

The balance equation (1) for the total normalized biovolume distribution was solved numerically using a zeroorder discrete population balance, or sectional method, in which the size distribution is described by uniformly sized bins, and the biovolume distribution is evaluated at a representative size centered on the bin[48]. The integrals for the aggregation (Eq. 2), fragmentation (Eq. 3), and transfer (Eq. 5) rates were discretized with the approach by [49]. The biovolume distribution evolved in time with a Runge-Kutta 4-5 adaptive time step scheme written in Fortran. The simulation results were fitted with the experimental data by minimizing the absolute error of each category’s 25th, 50th, and 75th percentiles of the biovolume distribution, using a Nelder-Mead algorithm [76]. The stickiness α₁ and for the pair of constitutive parameters S₁ and q₁ for the critical size of small colonies (category C₁ ) were obtained by fitting the model to the dataset of a single-cell suspension under and ϕ = 0.01. The pair of constitutive parameters S₂ and q₂ for the critical size of large colonies (category C₂) were obtained by fitting the model to the dataset of a large colony suspension under = 5.8 m²/s³ and ϕ = 0.01. These datasets were chosen as their aggregation/fragmentation times were similar to the measurement time, while for the other datasets, the variations were either too fast or slow to be captured in detail. The script for the simulation is available here.

Supporting Information Text

Field sample collection

Samples of Microcystis spp. were collected from Lake Gaasperplas - the Netherlands - at the sampling point 52°18’31.0”N 4°59’59.2”E on 14th August 2022, at 14:00. This is an artificial lake used for recreational purposes. No visible surface scums were present at the time of sampling. Water quality measurements of Lake Gaasperplas were obtained from the Dutch Water Authority (Rijkswaterstaat). Figure S1 shows micrographic images of the phytoplankton samples.

Image analysis

This section describes the routine used to analyze raw images of the aggregation and fragmentation process of Microcystis colonies under a cone-and-plate shear setup. Each frame has dimensions of 2048x2048 pixels. Prior to each experiment, 3.5 ml of distilled water was pipette in the setup, and 200 frames were captured at and 10 fps. A background image was generated by averaging over the frames and applying a Gaussian blur filter. After that, the distilled water was discarded, and 3.5 ml of cyanobacteria suspension was pipetted into the setup. Each time point in the data sets corresponds to a stack of 200 frames. Two treatments are made, one to detect small colonies(or single cells) and another to detect large colonies:

• Small colonies (l < 5): Image is cropped into a smaller section (500x500 pixels), to speed the processing. The background is subtracted using the reference image. The image intensity is renormalized and then thresholded. Small features (0.5 < l < 5) are segmented, measured, and counted.

• Large colonies (l > 5): A Gaussian blur is applied to filter small colonies. The background is subtracted using a reference image. The image intensity is renormalized and then thresholded. Small features (l > 5) are segmented, measured, and counted.

For each time point, the counted features of the 200 frames were distributed in bins of width Δl = 0.5. Since not all colonies are within the depth-of-focus, a correction function for the featured count was estimated using a suspension with a known size distribution, given by:

Here, N_p [i] and N_c [i] are the number of features before and after correction in a bin centered at l[i]. To obtain the calibration parameters λ_i, a sample of filtered colonies of Microcystis strain V163 was subjected to the cone-and-plate shear setup for 2 min at , and then transferred to a counting chamber. The latter was used as a standard for the size distribution. The calculated parameters were λ₁ = 1, λ₂ = 73.13, λ₃ = 13.71 and λ₄ = 1.24. The small size range (l ≤ 4) has low precision due to the pixel resolution, therefore the parameter λ₁ was adjusted to ensure biovolume conservation. The value of λ = 1 was used for t = 0, and it was fitted in each subsequent time point so that the total cell count is constant through time. Finally, the normalized biovolume distribution for each time point is calculated.

For the direct computation of the fragmentation frequency, a pure fragmentation with an ideal binary erosion was assumed, such that the time derivative of the colony size distribution is

where N_p [i](t), κ[i], and l[i] refer to the bin count, fragmentation frequency, and relative diameter at the center of bin i. This formula is only valid for the bins larger than the single cell, as the latter does not suffer fragmentation. The bin count distribution N_p [i] is initialized with the experimental value, and an initial guess fragmentation frequency κ[i] is made for all bins larger than a single cell. The values of N_p [i] are evolved in time using the Python solver scipy.integrate.solve_ivp, with an implicit Radau method. For each bin, the cumulative squared error between the experimental and numerical bin counts is computed. The error is minimized using the values of κ[i] as the optimization parameter and the Python solver scipy.optimize.minimize. The uncertainty of κ[i] is propagated from (S2) talking into account the bin count uncertainty. The bin size for each dataset was chosen in order to minimize the bin count uncertainty. The script for image analysis is available here.

Fragmentation measurements with Microcystis strain V163

Additional measurements of fragmentation of large colonies of Microcystis strain V163 were performed. Colonies grown by cell division were filtered through a mesh of 53 μm, and the collected colonies were redispersed in BG-11 medium. Under moderate dissipation rate of the fragmentation was negligible (Fig. S2), with the median size and fraction of large colonies remaining stable over 5 hours of flow. The behavior is independent of the total biovolume fraction, showing that aggregation rates are negligible for large colonies. In contrast, the steady median size for small colonies increases with the total biovolume fraction, although still an order of magnitude smaller than division-formed colonies. Significant fragmentation can be seen from a dissipation rate of (Fig. S3). For a large value of , the small colonies become the largest fraction in below 5 min (the smallest time resolution in our setup).

Aggregation measurements with Microcystis strain V163

Additional measurements of aggregation of single cells of Microcystis strain V163 were performed. The suspension was filtered through a 25 μm filter paper to select mostly single cells. Under moderate dissipation rate of , the aggregation increased with the total biovolume fraction ϕ (Fig. S4). Under an intense dissipation rate of , no aggregation was observed, and the median size of colonies remained equal to the initial value (near the single-cell size).

Rate of increase in the average colony size due to pure aggregation

Consider a suspension of colonies subjected only to aggregation, neglecting cell division and fragmentation. In this case, the time derivative of the normalized biovolume distribution n(l, t) is given by [1]:

where l_o = (l^f – l’^f)^1/f, and the collision frequency kernel is β(l,l’). The average colony size is defined as . Applying this definition to Eq. (S3), we obtain that the time derivate of (t) follows:

If we assume that the size distribution remains monodisperse for short times, we can write that n(l, t) = δ(l – (t)). Substituting this hypothesis into Eq. (S4) results that

Applying a variable change, the integral can be rewritten as:

Solving the integral and applying the model for the collision kernel in a turbulent shear flow [2], β(l, l’) = 0.163 (l + l’)³ , results that the rate of increase in the average colony size is:

In the limit of fractal dimension f = 2, the growth rate reduces to;

which is independent of .

Rheological characterization of concentrated colonies of Microcystis strain V163

A steady-shear (Fig. S5A) and a large-amplitude-oscillatory-shear (Fig. S5B) tests were performed in a sample of concentrated colonies of Microcystis strain V163. Due to the limited volume of samples available after separation, only one repetition was possible for each test. Before each test, the sample was subject to a preshear at for 200s, followed by 100s of rest. For the steady-shear test, the minimum transient time was determined in a separate test (40s for to 200s for ). For the oscillatory shear, the minimum transient time was determined by the rheometer software. A Herschel-Bulkley fit was applied to the steady shear data, , from which a dynamical yield stress τ_y = 4.3 ± 0.3 Pa (SD from best fit) was estimated.

Comparison of hypothesis for the fragment distribution function

Assuming a binary fragmentation (i.e., two fragments per fragmentation event), the fragment size distribution can be written in the form g_i(l, l’) = δ(l – l_a) + δ(l – l_b), where l_a and l_b are the size of fragments. Two ideal models for binary fragmentation are erosion and equal fragments, described by the following functions:

Due to a lack of experimental evidence of single colony fragmentation events, both models were tested for the category C₂ colonies. The equal fragments hypothesis was chosen as it better captures the experimental results for the size distribution of Microcystis strain V163 colonies aggregated from single cells (Fig. S6).

Cone-and-plate shear setup

The turbulent shear was generated by a cone-and-plate shear, with the dimensions depicted in Fig. S7A. The flow regime can be identified from the average energy dissipation rate curve against the angular velocity ω, with a ~ ω² dependence in the laminar regime and a ~ ω^5/2 dependence in the inertial region (Fig. S7B) [3, 4]. The Kolmogorov scale n, which is an estimate of the size of the smallest eddies in a turbulent flow, can be calculated from the relation . The Kolmogorov scale in our results range from η ~ 85 μm for the lowest dissipation rate () to η ~ 18 μm for the highest dissipation rate ().

Fig. S1

Fig. S2

Fig. S3

Fig. S4

Fig. S5

Fig. S6

Fig. S7

Acknowledgements

We thank Maria van Herk, Pieter Slot, and Merijn Schuurmans for their support with laboratory cultures. We are grateful to Axel Gunderson for his assistance with fieldwork, and to Daan Giesen for his support with experimental setup design. We especially thank Johan Oosterbaan and Richard Steel for their contributions during the project’s conception. We acknowledge funding from the Hoogheemraadschap van Rijnland. MJ acknowledges support from the ERC grant no. “2023-StG-101117025, FluMAB.