Algal blooms as a reactive dynamic response to seasonal perturbation in an experimental system

Algal blooms are typical of many aquatic freshwater ecosystems in seasonal environments. Such blooms could derive from transient reactive dynamics of algae and limiting nutrients following seasonal perturbation events. Linking parameter estimates derived from previously published lab experiments with empirical estimates of algal density dependence, we modeled dynamic interactions between nutrients and the green algal species Chlorella vulgaris and tested model predictions in a dozen 140 L mesocosms supplied with bi-weekly inputs of liquid fertilizer. Consistent with the reactive nutrient-driven model, Chlorella populations exhibited an initial surge in abundance over the first month followed by collapse as they rapidly converged on stable equilibria. The reactive model suggests that the magnitude of transient blooms is positively related to the augmentation of nutrients and depression of algae over the winter period. The magnitude of both algal peaks and equilibrium abundance was positively related to fertilizer loading, as predicted by the reactive model. Our results suggest that transient reactive responses to climate-driven perturbation events can be an important contributor to seasonal algal blooms observed in many temperate freshwater ecosystems. Controlled experimental studies such as ours may be helpful in understanding and potentially mediating the impact of fertilizer run-off on freshwater systems in temperate agricultural landscapes.


Introduction
Massive seasonal fluctuations in the abundance of phytoplankton populations are a common feature of many temperate aquatic ecosystems. While many factors are known to affect algal growth, a well-established body of theory suggests that algal growth is strongly shaped by temperature, light availability, and dynamic interactions with their limiting nutrients (O'Brien 1974;Grover 1989Grover , 1991aDeAngelis 1992;Huisman and Weissing 1995;Huppert et al. 2002Huppert et al. , 2005. Here, we evaluate the hypothesis that seasonal algal blooms represent a transient "reactive" response to perturbation of nutrient levels and algal abundance occurring during the non-growing season that has effects on dynamics during the subsequent growing season (Neubert and Caswell 1997;Neubert et al. 2004;Caswell and Neubert 2005). The defining characteristic of reactive models is that perturbation from stable equilibria can cause short-term amplification of the disturbance, resulting in massive fluctuation in state variables before the system settles back into an equilibrium (Neubert and Caswell 1997).
While 12-month periodicity in phytoplankton abundance presumably due to annual climatic variation has been observed in 48% of 125-time series datasets analyzed by Winder and Cloern (2010), cycle lengths of 6 months (20%) or irregular patterns of fluctuation (30%) were also common. Following on previous models of resource pulse dynamics (Holt 2008, Steiner et al. 2009, Klausmeier 2010, Holland et al. 2015, we accordingly assume bi-phasic (square-wave) alternation between a non-growing season with inadequate environmental conditions to support phytoplankton growth and a growing season in which conditions are sufficient for phytoplankton to flourish. The alternation between growing and non-growing phases of growth, caused by seasonal variation in solar radiation and temperature, has been clearly demonstrated as a diagnostic feature of algal fluctuations in one of the most intensively studied lake systems, Lake Washington (Arhonditsis et al. 2004). In keeping with common ecological usage (Bender et al. 1984), we treat the changes in phytoplankton and nutrient abundance that occur during the non-growing season as a pulse perturbation away from equilibria conditions that apply during the growing season.
Such perturbation should be routinely expected in temperate freshwater ecosystems, for example, due to algal dieoff during winter (Arhonditsis et al. 2004) accompanied by accumulation of nutrients due to fertilizer run-off during precipitation events (Novak et al. 2003;Kim et al. 2006). The logical consequence of such perturbations is that algal populations at the beginning of the growing season would be substantially depressed at the same time as nutrient abundance is substantially augmented. This combination of perturbations could be particularly potent in generating reactive algal response as algal growth begins anew at the beginning of the next growing season. There could be important practical benefits in better understanding environmental drivers of such reactive nutrient-driven algal blooms since the frequency and magnitude of harmful algal blooms are thought to have increased in recent years (Michalak et al. 2013;Paerl and Otten 2013). This increasing trend is pronounced in global datasets for large lakes (Ho et al. 2019), although data for smaller lakes in North America show little trend over a 10-year period spanning the beginning of the twenty-first century (Wilkinson et al. 2021).
Here, we use time series trials of 190 days length in 140 L mesocosms to experimentally test for reactive transient dynamics using a nutrient-driven model of population dynamics of Chlorella vulgaris, a globally common species of green algae often used for biofuel projects. Detailed experimental data were already available on the effect of nutrient concentration on algal growth rates, which we linked with density-dependent feedbacks parameterized for our experimental system. We use this parameterized lab model to estimate the impact of nutrient and algal perturbations on the reactivity, resilience, and stability of the nutrient-Chlorella system. We have also experimentally tested the ability of the nutrient-driven algal model to predict the impact of augmented fertilizer input on the amplitude of blooms and equilibrium abundance of Chlorella populations.

Model
Our nutrient-driven model was based on a substantial body of previous theoretical and empirical work on nutrientphytoplankton dynamics (e.g. O'Brien 1974;Grover 1989;DeAngelis 1992;Huppert et al. 2002Huppert et al. , 2005. Our model links a closed population of Chlorella vulgaris (P) whose growth rate is shaped by both density-dependent processes, estimated in our experiments, and potentially limiting levels of nutrients (N), estimated in earlier benchtop experiments (Pannikov and Pirt 1978).
where r max is the maximum per capita phytoplankton growth rate (cells per day), a is the daily rate of nutrient influx (μg/ mL), b max is the maximum rate of nutrient absorption per phytoplankton cell (μg/cell-day), c is the nutrient concentration (μg/mL) at which nutrient absorption is equal to b max /2, d is the nutrient concentration (μg/mL) at which phytoplankton growth rate is equal to r max /2, φ is a scaling constant linking maximum rates of nutrient absorption and phytoplankton growth to local environmental conditions, and γ is the the log-linear (Gompertz logistic) effect of algal density on the exponential rate of algal population growth. Parameters b max , c, r max , and d were obtained from Panikov and Pirt's (1978) benchtop experiments spanning a wide range of concentrations of nitrogen or phosphorus in a flow cytometer with high flushing rates to maintain extremely low population density of Chlorella vulgaris cells.
Although 6 parameters are required by our model of nutrient-phytoplankton dynamics, one parameter was experimentally controlled (a) and three other parameters have been estimated directly from previously published experimental trials for Chlorella vulgaris (b, c, and d). That left just two parameters to be estimated from our phytoplankton time series experiments: φ (specifying the maximum rate of phytoplankton growth and nutrient absorption under our experimental conditions) and γ (specifying the densitydependent effect of phytoplankton on per capita mortality rates). For simplicity, we assumed that phytoplankton growth and absorption are scaled in the same way by local environmental conditions, so we applied the same coefficient. We systematically sampled across a wide range of parameter values and used least-squares minimization to identify the most parsimonious set of parameters to explain algal population density over time, when grown in isolation. The likelihood surfaces for these 2 parameters were smooth and unimodal, so the fitting procedure was quite straightforward. Algal density and nutrient levels were initialized in the model with the same starting values as we used in the mesocosm experiments (P 0 = 50,000 cells/mL and N 0 = 2 μg/ mL). R code for parameter estimation and model simulation is available on figshare.
A useful first step towards understanding consumerresource models such as that represented by Eqs. (1) and (2) is by evaluating a linearized version of the model in the vicinity of the system equilibrium, whose rates of change are specified by the Jacobian matrix A, whose elements are calculated in the following manner: The stability and reactivity of the system are then evaluated by applying the Jacobian matrix to the linearized system around the equilibrium.
where n = N − N eq and p = P − P eq represent deviations of nutrient concentration (N) and phytoplankton density (P) from system equilibria (N eq and P eq ). If the dominant eigenvalue (λ 1 ) of A < 0, then, the system is locally stable and in the long term, the system will recover equilibrium levels of nutrients and algae following perturbation. The speed of long-term recovery in such a stable system is scaled with the magnitude of λ 1 (A), often referred to as system resilience (Pimm and Lawton 1977;Neubert and Caswell 1997;Ives et al. 2003).
The magnitude of λ 1 (A) tells us little, however, about the immediate response to perturbation, and this is relevant to seasonal algal blooms because many dynamic systems can experience further deviation in the short term before they eventually converge on the stable equilibrium over the longer term (Neubert and Caswell 1997). To evaluate that possibility, we calculated the dominant eigenvalue λ 1 (H) of the Hermitian part of the Jacobian matrix: H = (A + A T )/2, where A T refers to the transpose of matrix A, often referred to as the reactivity of the system (Neubert and Caswell 1997). A system will be reactive provided that λ 1 (H) > 0, implying that at least some perturbations will result in amplified deviation from the stable equilibrium in the short term, rather than reduction. The degree of amplification scales with the magnitude of positive values of λ 1 (H). We used experimental estimates of the matrices A and H to evaluate the long-term stability and short-term reactivity of nutrientlimited Chlorella populations following perturbation from equilibrium.

Materials and methods
The Aqualab facility at the University of Guelph provided a temperature-and light-controlled environment for the experimental array of 12 fiberglass tanks of 140 L volume used in this study. The mesocosms were cylindrical, with a funnel-shaped bottom, and had a diameter of 60 cm and height of 70 cm. An independent LED lighting system suspended within each tank provided sufficient photosynthetically active radiation for rapid phytoplankton growth. We used filtered well water as our medium. This water was first passed through three Pentek-pleated filters (25, 1, and 0.2 µm) to remove both large and small sediment particles and through a UV-light unit to remove any remaining bacteria and cysts. After filling with well water, each tank was allowed to stabilize for 5-7 days to a temperature of 20 °C before algae were added. To mimic natural patterns of nutrient flushing in freshwater systems, tanks were provided with either 0.25 or 0.5 mL of liquid fertilizer (Scotts Miracle-Gro LiquaFeed © all-purpose 12-4-8 fertilizer) twice a week, which works out to 0.06 μg/mLday of nitrogen and 0.02 μg/mL-day of phosphorus. Fertilizer was added only after sampling was completed on any given day. Chlorella vulgaris cultures were obtained from the Canadian Phycological Culture Centre (CPCC) at the University of Waterloo, Canada, (strain #90) and kept in COMBO at 20 °C with 16 h:8 h dark:light cycle. Twice a week (Mondays and Thursdays), we sampled algal density in each tank. Algal density was calculated as the average of three readings (3 mL of water for each reading) obtained with a fluorometer (Trilogy Laboratory Fluorometer, Turner Designs, San Jose, CA, USA) cued to the concentration of chlorophyll in vivo (excitation 441/82 nm; emission 660-710 nm). We transformed fluorometer readings into algal density (cell/mL) using a calibration curve obtained by counting number of cells under the microscope (n = 44; calibration curve: cells/mL = 0.57618 × average fluorometer reading + 5.57760, R 2 = 0.64). We inoculated each tank with 0.050 × 10 6 cells of algae and monitored subsequent algal population changes in algae over the next 190 days with bi-weekly counts.

Results
Exponential growth rates of green algae declined with density (Fig. 1a) and a log-linear model ( Fig. 1c; y = 1.344 -0.104log(x), F 1,316 = 83.28, p < 0.001, R 2 = 0.206) provided substantially better fit to our data on the exponential growth rates of green algae (ΔAIC = 45.3) than a linear density-dependent model ( Fig. 1b; y = 0.0742 -0.000000136x, F 1,316 = 30.27, p < 0.001, R 2 = 0.085). Water quality samples taken at the beginning (mean (N/P) = 143.7, sd(N/P) = 96.2) vs. end of our experiments (mean(N/P) = 5.9, sd(N/P) = 3.1) exhibited 24-fold decrease in the ratio of nitrogen to phosphorus, suggesting that nitrogen was being taken up by algae at a much faster rate than phosphorus, implying nitrogen as a limiting resource. Given this trend, we assumed nitrogen limitation was most appropriate for our system and incorporated a log-linear density-dependent effect on algal growth. We note that alternate models based on phosphorus limitation yielded qualitatively similar outcomes.
Green algal abundance rapidly increased following initiation of the trials, typically reaching a peak during weeks 2-4 of the trial (Fig. 1a). Maximum densities across 6 replicates averaged 1.47 × 10 6 cells/mL. After 65 days, algal densities dropped dramatically, leveling off at an average density of 0.31 × 10 6 cell/mL (Fig. 1a). After fitting the 2 free parameters (φ and γ) and initiating the model with the same algal density (0.050 × 10 6 cells) with nitrogen concentration (2 μg/L) observed at the beginning of lab trials, our nutrient-driven model (Eqs. 1 and 2) echoed key features observed in the algal time series, predicting a sharp peak of 1.49 × 10 6 cell/mL in the 3rd week, followed by rapid decline and stabilization by day 60 at a density of ~ 0.47 × 10 6 cells/ mL (Fig. 1a).
For the set of parameters used to model the nutrient-Chlorella system (a = 0.06 μg nitrogen/mL per day, b max = 7.44 10 −6 μg nitrogen/mL per day, c = 0.094 μg nitrogen/mL, d = 0.082 μg nitrogen/mL, r max = 6.00 per day, φ = 0.05, and γ = 0.0086 per day), the elements of the Jacobian matrix A were calculated as follows: a 11 = − 0.800, a 12 = − 1.28 × 10 −7 , a 21 = 3.86 × 10 6 , and a 22 = − 0.0086. The dominant eigenvalue of the resultant transition matrix (A) was less than 0 (λ = − 0.40), so our model suggests that the nutrient-Chlorella system should be locally stable. On the other hand, the dominant eigenvalue of the Hermetian part of the transition matrix (H) was much larger than 0 (ρ = 1,931,265), so our model suggests that the system should be highly reactive and therefore prone to large excursions in the short term before eventually converging on the system equilibrium. Both predictions were consistent with our experimental results.
We explored the starting conditions most likely to stimulate a massive reactive response by perturbing nutrient Chlorella densities are decreased to 1/2, 1/4, 1/8, or 1/16 the equilibrium level; c, d nutrient levels are augmented to 2 × , 4 × , 8 × , or 16 × the equilibrium level; or e, f both nutrients and Chlorella are perturbed simultaneously levels, Chlorella levels, and both nutrient and Chlorella levels over a range of initial values and then modeling subsequent transient dynamics (Appendix 1). Reactivity was evident in time series realizations of the model following perturbations in which algal density was reduced below equilibrium levels (Fig. 2a), nutrient level was augmented above the equilibrium (Fig. 2c), or algal density was reduced while nutrient levels were augmented (Fig. 2e). Such perturbations are typical of seasonal temperate systems and was mimicked by our experimental protocol.
The parameterized Chlorella growth model predicts that minor perturbations of either algal density (e.g., 1/2 or 1/4 the equilibrium level of algal density) or nutrient levels (e.g., 2 × the equilibrium concentration of nutrients) should not induce algal blooms as the system equilibrates, whereas perturbations of greater magnitude are capable of generating blooms (Fig. 2a, c). For sufficiently large perturbation of algal density, nutrients, or both, the amplitude of the resulting algal bloom positively scales with the magnitude of perturbation (Fig. 3). The sensitivity of the system response differs, however, between state variables (Fig. 3). Our model predicts that the magnitude of bloom response to a given magnitude of nutrient perturbation should be roughly double that of algal perturbation of identical magnitude (i.e., comparison of reduction of algal density to 1/16 the equilibrium value vs. 16-fold increase in nutrient level). In other words, nutrient supply is a substantially stronger driver of bloom magnitude in this system than is initial algal density (Fig. 3). Pairing algal reduction with nutrient augmentation should result in blooms of greater magnitude than either algal reduction or nutrient augmentation in isolation (Fig. 3) but have little effect on the time to equilibration (Fig. 2f) compared to algal perturbation by itself (Fig. 2b). In other words, nutrient perturbation increases the degree of system reactivity beyond that of algal reduction alone (Fig. 2f), contributing synergistically to increased bloom amplitude.
Nutrient perturbation by itself (Fig. 2d) should result in blooms of much shorter duration than equivalent perturbations in algal density (Fig. 2b). Interestingly, nutrient perturbation by itself does not generate a reactive system response, since the deviation from the equilibrium is predicted to immediately decay following perturbation rather than growing for the first while (2d). Both of these characteristics stem from differences in the growth constraints of each state variable. When algae are uncommon, nutrients increase in concentration at a constant rate. Algal growth, on the other hand, is constrained by population size even when nutrients are super-abundant. At the beginning of the growing season, following perturbation from equilibrium levels that were in play at the end of the previous growing season, nutrients increase at a maximum rate, whereas algal numbers are incremented much more slowly simply because the population is small. This fundamental difference in initial growth rates results in a reactive system, whereby deviations initially grow before stabilizing and eventually decaying as the system equilibrates.
Our model predicts that doubling the rate of nitrogen influx to 0.12 μg nitrogen per mL-day should lead to a substantial increase in peak (2.10 × 10 6 vs. 1.49 × 10 6 cells/mL) and equilibrium algal abundance (0.89 × 10 6 vs. 0.47 × 10 6 cells/mL) relative to outcomes with nutrient influx = 0.06 μg nitrogen per mL-day (Fig. 3a). Our experimental trials with double the daily influx of fertilizer experienced fourfold increase in peak algal density (5.82 × 10 6 cells/mL), significantly higher (t = 23.80, df = 8.04, p < 0.001) than maximum algal densities in the low fertilizer trials. Experimental trials with double the daily influx of fertilizer experienced a 2.8-fold increase in algal density over the last 120 days of the experiment (0.536 × 10 6 cells/mL), significantly higher (t = 5.66, df = 7.55, p < 0.001) than mean equilibrium density in the latter part of the low fertilizer trials.

Discussion
Results from our study suggest that transient interactions with their nutrients can play an important role in shaping the magnitude of seasonal algal blooms in freshwater ecosystems. As predicted by the nutrient-driven model, Chlorella populations invariably went through a massive bloom when placed in well water and supplied with a steady influx of fertilizer. The magnitude of our experimental algal blooms was strongly responsive to increased levels of nutrient loading. Nutrient levels in the water column are well-known to be associated with mean levels of both plankton productivity and biomass (Tilman et al. 1982;DeAngelis 1992;Elser et al. 2000;Carpenter et al. 2001;Carpenter 2016). There is also compelling evidence of algal variation over time being linked to seasonal variation in nutrient levels in some lakes (DeAngelis 1992;Vasseur et al. 2005) and rivers (Huppert et al. 2002).
Our phytoplankton model parameterized for Chlorellanutrient interactions in our experimental microcosms builds on a substantial body of previous theoretical and experimental work showing that water turbidity, incident light intensity, water depth, cell size and nutrient concentration within cells, nutrient influx, and nutrient loss rates shape phytoplankton dynamics (O'Brien 1974;Grover 1991a, b;DeAngelis 1992;Huisman and Weissing 1995;Huisman 1999;Passarge et al. 2006). Algal-nutrient models suggest that algal populations presented with a surge in nutrient availability at the beginning of the growing season combined with depressed algal abundance due to over-winter mortality should be expected to exhibit pronounced blooms of initial abundance followed by a subsequent collapse to greatly reduced densities as the nutrient pool equilibrates with phytoplankton demand (O'Brien 1974;DeAngelis 1992;Huppert et al. 2002). Similar patterns are often observed in both laboratory and natural ecosystems (Vasseur et al. 2005;Huppert et al. 2002). Even at the elevated growth rates displayed by Chlorella vulgaris, such transient blooms can last several weeks of the growing season in temperate environments, so could potentially impact food web dynamics.
Our model assumes density-dependent feedback on per capita growth rates by Chlorella vulgaris. A well-developed body of theory (Huisman and Weissing 1995) and experimental work with Chlorella vulgaris both as a monoculture (Huisman 1999) and in competition with other green algal species (Passarge et al. 2006) suggests that competition for light was responsible for the density-dependent effect on algal growth rates that we observed. This interpretation is supported by additional benchtop batch experiments (to be analyzed in depth in a separate publication) with Chlorella raised in vials, in which one set of replicates was positioned behind the second set of replicates, such that 1/2 the vials were shaded from full exposure to light. This experiment resulted in depressed abundance in the shaded vials, suggesting that light competition was a strong contributor to density-dependent effects. Such so-called "closure" due to density-dependent feedbacks in consumer-resource models of this type can be of critical dynamical importance (Steele and Henderson 1992), in part because it can potentially suppress reactivity (Neubert et al. 2004). For small perturbations of nutrients or algae, density-dependence might explain why blooms are unlikely. For larger perturbations, however, the strength of density-dependent regulation was apparently of insufficient magnitude to prevent extreme transient excursions in algal abundance. Such transient algal population dynamics are particularly likely under the kind of stochastic perturbations likely to occur in temperate seasonal ecosystems due to variation in temperature during the non-growing season, length of the non-growing period, frequency of precipitation events, and variation in nutrient load arriving from the surrounding environment.
It is likely that more precise model predictions could be obtained by building a more complex model controlling for the impact of nutrient uptake in relation to cell size, nutrient concentration within algal cells and in the water column, light incidence, water turbidity, stoichiometry, and spatial distribution of both light and plankton at depth (Grover 1991a;Huisman 1999;Elser et al. 2000;Heggerud et al. 2020;Wang et al. 2022). Both the Monod and Droop 1 3 formulations are commonly used to model nutrient-phytoplankton interactions (Burmaster 1979;DeAngelis 1992;Wang et al. 2022). Both Droop and Monod formulations generate similar equilibrium outcomes (Burmaster 1979) but differ in transient responses if phytoplankton cells are capable of storing excess quantities of nutrients within the cell, thus allowing cells to continue to divide after nutrient levels have declined (Wang et al. 2022). For analytical tractability in this instance, however, we opted for the more stream-lined Monod formulation to allow us to estimate density-dependence effects from time series data as well as simplify the mathematical analysis of the impact of initial nutrient and plankton abundance and nutrient influx rates on system reactivity.
Empirical evidence from other perturbed ecosystems suggests that transient blooms induced by reactive dynamics between algae and their limiting nutrients could have important dynamic ramifications across the rest of the food web (Holt 2008;Hastings 2010;Yang et al. 2010). Like many other food web models, a wide range of outcomes are possible. Herbivorous consumers should benefit from algal blooms, leading to short-term bursts in consumer abundance. While mild surges in consumer abundance would be expected to eventually flatten out as the system converges again on the stable equilibrium, particularly strong population responses by consumers could lead to subsequent collapse by consumer populations due to overcompensatory density-dependent responses (Holt 2008). Consumer-resource models based on a Lotka-Volterra framework predict that the magnitude of the pulsed resource response should positively influence the strength of consumer response (Holt 2008). This prediction is consistent with a meta-analysis of 161 field and experimental studies of resource pulse effects on consumer-resource interactions (Yang et al. 2010). The strength of consumer response is typically less than that of the resource pulse (Yang et al. 2010), suggesting that massive algal blooms might be expected to support a stronger herbivore response than minor blooms.
Resource pulses of long duration would be expected to have a stronger impact on herbivore numerical responses than short-lived pulses (Holt 2008), and this prediction was also supported by Yang et al.'s (2010) meta-analysis. Our models suggest that green algal blooms induced by reductions in algal density and increases in nutrient load during the non-growing season should equilibrate on the order of a month or so, whereas surges in nutrient load should equilibrate much more quickly. We conjecture that reductions in algal density combined with nutrient augmentation should have a stronger effect on herbivore populations than perturbations of nutrients by themselves, due to both increased bloom magnitude and the increased duration of the perturbed resource response.
Carnivorous zooplankton and fish populations at even higher trophic levels would be expected to respond to seasonally forced algal perturbations through their own demographic responses to temporary surges in herbivore abundance (Steele and Henderson 1992;Carpenter et al. 2001;Yang et al. 2008). Interestingly, Yang et al.'s (2010) meta-analysis suggests that population responses to resource perturbation by species that are 2 trophic levels above the perturbed resource are often stronger than that of species 1 link above the resource. This suggests that carnivorous zooplankton and fish populations might be expected to respond even more strongly than herbivorous zooplankton to seasonally-forced algal blooms, all other things being equal. Depending on the degree of closure of top predators, such interactions might themselves be highly reactive, leading to cascading episodes of transient population surges across the food web (Yang et al. 2010) and potentially even more complex community dynamics.
Given that multiple plankton species compete for light and the same set of limiting nutrients, blooms of variable magnitude and duration due to stochastic variation in environmental conditions could influence the outcome of lottery competition and thereby alter species coexistence (Chesson 1994). Particularly strong blooms also have the potential to push systems across critical thresholds, leading to alternative stable states (Scheffer et al. 2001). Since there are plausible theoretical reasons to expect multiple equilibria in freshwater ecosystems (Scheffer et al. 1993, Scheffer and, such episodic stochastic events could have far-reaching effects on community structure and long-term dynamics. It is increasingly obvious that coherent conservation and management planning will depend in large part on our ability to understand both short-and long-term responses to the inevitable perturbations that characterize real ecosystems (Hastings 2010). Data Availability Source data and R code for Fig. 1 and 2

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.