Phytoplankton mortality in a changing thermal seascape

Predicting spatiotemporal distributions of phytoplankton biomass and community composition heavily relies on experimental studies that document how environmental conditions influence population growth rates. In unicellular phytoplankton, the net population growth rate is the difference between the cell division rate and the death rate. Along with predation and disease, phytoplankton mortality arises from abiotic stress. Although the effect of temperature on the net population growth rate is well understood, studies examining thermally induced death rates in phytoplankton are scarce. We investigated how cell division and death rates of the diatom Phaeodactylum tricornutum varied within its thermal tolerance limits (thermal niche), and at temperatures just above its upper thermal tolerance limit. We show that death rates were largely independent of temperature when P. tricornutum was grown within its thermal niche, but increased significantly at temperatures that approached or exceeded its upper thermal tolerance limit. Furthermore, the sensitivity of mortality increased with the duration of exposure to heat stress and was affected by the pre‐acclimation temperature. Heat waves can be expected to significantly affect phytoplankton mortality episodically. The increasing frequency of heat waves accompanying global warming can be expected to drive changes in phytoplankton community structure due to interspecific variability of thermal niches with potential implications for food web dynamics and biogeochemical cycles.


| INTRODUC TI ON
Phytoplankton growth generates half of the atmosphere's oxygen through photosynthesis, thus providing the energy that fuels pelagic food chains and playing key roles in biogeochemical cycles of oxygen, carbon and nutrient elements. Predicting spatiotemporal distributions of phytoplankton biomass and community composition heavily relies on experimental studies that document how environmental conditions influence net population growth rates, the difference between cell division rate and the death rate. Although the environmental controls on phytoplankton growth are well studied, the factors governing cell death rates remain unclear.
Death rates in natural phytoplankton communities range from <0.01 to >1 day −1 (Marbá et al., 2007). Studies that have calculated death rates from measurements of cell lysis (Agustí & Duarte, 2000;Agustí et al., 1998) include contributions to death due to biotic factors such as viral infection as well as abiotic stress-induced physiological death. The death rate due to physiological stress is typically low under optimal growth conditions, but can increase significantly in response to physiological stresses associated with nutrient starvation and prolonged darkness, although the tolerance to different stresses varies markedly among species (Berges & Falkowski, 1998;Brussaard et al., 1997;Timmermans et al., 2007).
Across the plant and animal kingdoms, ambient temperature, due to its influence on metabolic rate and body size, has been shown to govern both net growth rates and death rates, whereby organisms with higher metabolism have higher birth rates but shorter lifespans (McCoy & Gillooly, 2008). Although these relationships are useful for describing current patterns in natural systems, to forecast how species distributions will respond to a changing climate we need to know whether death rates are influenced by abiotic stress.
Due to the increasing frequency of heat waves accompanying global warming in both terrestrial and aquatic systems (Frölicher et al., 2018;Ruthrof et al., 2018), thermally induced mortality plays a greater role in shaping community structure (Remy et al., 2017). Consequently, a large number of studies have focused on how thermal conditions (i.e. acute or chronic heat stress) induce organism and/or cell death across a wide range of taxa.
Inopportunely, many studies such as these do not report both cell division and death rates, and hence the impact of physiologically induced death on population responses to environmental stressors remains unclear. Experiments that report the effect of thermal stress on mortality as the proportional change in dead or live cells (e.g., Bouchard & Yamasaki, 2008;Samuels et al., 2021;Zuppini et al., 2007) cannot be used to calculate death rates. Hence, studies that directly quantify both cell division and death rates from changes in the abundances of live and dead cells in response to heat stress are needed if we are to understand how warming will alter populations and community composition.
In response to warming (or cooling), an individual organism is able to make transient cellular adjustments to physiology and biochemistry via process known as acclimation (Angilletta Jr. et al., 2003), and in doing so, can alter their thermal niche (Luhring & DeLong, 2017;Padfield et al., 2016;Schulte et al., 2011). Indeed, thermal history has been shown to influence phytoplankton survivorship in the context of heat waves (Samuels et al., 2021) and understanding whether acclimation plays a role in governing death rates is warranted.

| MODELLING
Pelagic ecosystem models typically focus on the ecological interactions (disease and predation) that affect phytoplankton mortality, neglecting mortality that arises spontaneously or in response to adverse abiotic conditions (Follows et al., 2007). As far as we are aware, three alternative models have been used previously to model the temperature dependencies of cell division rate (μ) and death rate (d) in microalgae, with the specific growth rate of a population (r) calculated as the difference between the cell division rate and the death rate: where r, μ and d are calculated as described in Section 3.3.
In the first approach, death rates are assumed to be independent of temperature, whereas the cell division rate is assumed to follow a unimodal dependence on temperature, increasing with temperature below the thermal optimum (T opt ) and declining with temperature above T opt (Figure 1a). This is an approach used in many pelagic ecosystem models (Follows et al., 2007).
F I G U R E 1 Different approaches that currently exist to account for the temperature dependence of net growth rate (solid black line). (a) Death rate (dashed line) is assumed to be temperature independent, hence is 'fixed', whereas the cell division rate (gross growth rate; grey sold line) is assumed to follow a unimodal dependence on temperature, whereby growth rate increases with temperature below the thermal optima (T opt ) and declines with temperature above T opt . In an alternative approach depicted in (b), both cell division rate and death rate are assumed to increase with temperature, but with the death rate increasing more rapidly than cell division rate. Both approaches result in a net growth rate (black solid line) that exhibits a unimodal, asymmetrical curve known as the thermal performance curve (TPC) Alternatively, in the second approach, Serra-Maia et al. (2016) assumed that both the cell division rate and death rate obey the Arrhenius equation: where T is the absolute temperature (K), A 1 and A 2 are pre-exponential factors (units of day −1 ), E 1 and E 2 are the activation energies for the cell division and cell death terms and R is the universal gas constant.
Although there is a mechanistic basis for this equation for chemical reactions, it should be considered empirical when applied to the complex processes of cell growth and cell death.
In the third approach, Thomas et al. (2017) proposed a double exponential model (Figure 1b), in which the cell division rate is an exponential function of temperature, and mortality includes a baseline death rate (d 0 ) and an exponential increasing rate, but with the death rate increasing more rapidly than cell division rate: where μ 1 is the cell division rate at 0℃, and μ 2 is the exponential change in this rate with increasing temperature, d 0 is a temperatureindependent death rate, while d 1 and d 2 jointly describe the exponential increase in death rate with temperature.
Although all of these approaches yield similar TPCs for net growth, zero net growth rates are achieved with differing assumptions for the temperature dependence of mortality and consequently depict strikingly different physiological responses to heat stress.
Although both the Serra-Maia et al. (2016) and Thomas et al. (2017) models assume high death rates within the thermal niche, datasets that quantify rates of phytoplankton mortality are scarce and would be valuable for those ecosystem models that incorporate mortality.
Furthermore, whether the decline in net population growth rate (r) due to high temperature stress at temperatures between T opt and maximum temperature limit for growth (T max ) is primarily due to a decline of cell division or an increase in morality remains to be quantified. In this study, we specifically address the following key ques-
Cultures were acclimated for more than 10 generations before an inoculum from an exponentially growing batch culture was transferred into fresh medium to establish biological replicates (n ≥ 3).
Similarly, during exponential growth, an inoculum from each parent culture was then diluted 1:100-1000 (depending on the growth rate) into sterile, acid-washed borosilicate culture tubes (Fisherbrand™, catalogue number: 11537423, Fisher Scientific UK Ltd) to initiate a number of daughter cultures (5 ml working volume). These daughter cultures were then distributed along a temperature gradient that was established using a tube-based thermal gradient block, where hot and cold water generated using separate immersion heater chillers (LT Ecocool™ 100; Grant Instruments) was circulated through milled channels at opposite ends of the block. Once housed in the thermal gradient block, individual tubes were monitored daily by tracking the abundances of live and dead cells as described below.
We found that the inoculation process induced some cell death, and in some tubes could contribute up to 20% of the initial 'dead' population. Consequently, we carried out the inoculation process and incubated the daughter tubes overnight at the parent culture conditions using the thermo-gradient block. This allowed cultures to recover from shock induced by transfer before they were subsequently heat stressed and resulted in starting 'live' populations greater than 90%.

| Cell staining and flow cytometric measurements
To assess mortality and estimate the growth, subsamples (0.12 ml) from each tube were removed and incubated in the presence of a nucleic acid stain (SYTOX green, Molecular Probes) to quantify live and dead cells following a protocol based on Peperzak and Brussaard (2011). Specifically, a working solution of SYTOX green (50 μM) was prepared in dimethylsulfoxide (DMSO) and stored at −20℃ until just before use. This working solution (2 µl) was then added to the 0.12 ml volumes that had been removed from individual tubes (0.5 μM final concentration) and incubated in the dark at room temperature for 20 min. Incubated samples were then analysed by a flow cytometer (BD Accuri C6 Plus, Becton Dickinson) equipped with a 20 mW blue (488 nm) laser.
Using instrument-specific software (BD Accuri™ C6 Plus Software, Version 1.0.27.1), populations of P. tricornutum were first discriminated on bi-plots of forward-scatter (FSC; a proxy for cell size) versus chlorophyll-a (>670 nm) fluorescence. Gated diatoms were then discriminated as viable 'live' and non-viable 'dead' using green fluorescence (533/30 nm) histograms, where the number of cells permeable and impermeable to SYTOX-green were counted and their respective proportions of the total population calculated. Representative gating is shown in Figure S1. Heat-killed cells of P. tricornutum (60℃ for 60 min) were used as positive controls for SYTOX staining; non-viable cells were more than two decades brighter than the unstained, viable cells. The gating of the diatom populations and subsequent gating of 'live' and 'dead' populations were post-processed as both the FSC and chlorophyll fluorescence shifted across experimental treatments. All data were acquired with the same instrument settings.

| Calculation of net, gross growth and death rates
The population dynamics of live and dead cells is assumed to be given by: where μ is the cell division rate (units of day −1 ), d is the death rate (units of day −1 ), N L is the abundance of live cells and N D is the abundance of dead cells.
The specific growth rate (r), also known as the instantaneous rate of increase (Wood et al., 2005), is the difference between cell division and death rates: The specific growth rate, r, was calculated from: The change in the number of dead cells per ml was calculated as the abundance of dead cells at time t, subtracted from the abundance of dead cells at time t + 1: The average abundance of live cells, N L in a time period was calculated from Equation (5), where r is the specific growth rate that results in new live cells: Death rate (d) is commonly defined as the proportion or percentage of a population that dies within a given time interval. For some organisms, the time interval is usually a year, but for microalgae a more appropriate time interval is a day. Within the thermal niche (T < T max ), death rate (day −1 ), d, was calculated as: This approach gives a minimum estimate of death rate because a proportion of dead cells may be underestimated due to their change in cell properties. Phytoplankton death results in degradation of photopigments followed by disintegration into unrecognizable debris. Hence, dead cells may lie outside the gating strategy outlined in Figure S1. As such, death rates may be underestimated, particularly at the most extreme temperatures as this process can occur within a day (Veldhuis et al., 2001).
At temperatures where net population growth was negative (T > T max ), death rate was calculated from the decline in the abundance of live cells.

| Modelling temperature dependence of growth
We fit 12 equations from the literature available in the R pack- The parameter values in Equation (8) were found by maximum likelihood estimation and nonlinear least-squares fits (using the Levenberg-Marquardt method), respectively, and confidence intervals were found by bootstrapping using the R package 'modelr' (Wickham, 2019) to generate 10,000 replicates of the data, each of which was randomly sampled with replacement. The T opt , μ max , T min (lower temperature limit for growth) and T max were then calculated for each bootstrap replicate. The 95% CI were then calculated as the range between the 2.5th and 97.5th quantiles from the distributions of the derived parameters.

| Modelling temperature dependence of mortality: Post-hoc analysis
We visually evaluated our dataset post-collection and found that the temperature dependence of death rates did not appear to support any of the three approaches outlined in the introduction above.
Instead, death rates were essentially zero for the majority of temperatures examined below a threshold; following this threshold, death rates increased significantly. Consequently, we conducted a post-hoc exploratory analysis in the form of piecewise regression.
Piecewise regression techniques have been used as a statistical method to model ecological thresholds (Toms & Lesperance, 2003).
These models are also known as broken-stick models, the simplest of which is where two straight lines are joined sharply at a breakpoint (to represent a threshold) and can be described as: where d i is the mortality for the ith temperature, T i is the corresponding temperature, T mort is the thermal threshold for mortality (i.e. the breakpoint) and e i are residual errors. The slopes of the lines are β 1 and β 1 + β 2 , and hence β 2 can be interpreted as the increase in the sensitivity to thermal stress.
Here, we use piecewise regression to identify the temperature (℃), that is, breakpoint(s) (T mort ) at which death rates change abruptly (from approximately zero). Both linear and piecewise regression models were fitted to the relationship between temperature (℃) and death rate. Segmented regression analysis was conducted using the 'segmented' package (Muggeo, 2015) of the R software (Version 1.1.442). The significance of this breakpoint was assessed using sequential hypothesis testing via the Score test.

| RE SULTS
Segmented regression analysis identified a relationship between mortality and temperature that was dependent on the thermal range examined. Death rates were largely independent of temperature when P. tricornutum was exposed to temperatures within the thermal niche, and only become significant when cultures were exposed to temperatures that approached or exceeded the T max . Preacclimation temperatures and the duration of exposure can alter the thermal threshold for mortality (T mort ) that, in turn, has a relation to T max for growth.
When P. tricornutum was exposed to temperatures within its known thermal niche ~10-30℃ (Siegel et al., 2020), we found no evidence of a relationship between mortality and temperature ( Figure 2). In this thermal range, death rates were negligible even at the longest exposure periods and irrespective of pre-acclimation conditions (β 1 ≤ 0.001 day −1 · ℃ for 10, 20 and 26.5℃ at 72 h;  Figure S2). We observed no mortality until temperatures equalled or exceeded 27.1℃ (95% CI: 26.7, 27.6)-a thermal threshold which triggered a breakpoint, that is, a change in slope. At these breakpoints, mortality increased significantly with temperature (β 2 ; p < 0.001) with death rates in excess of 0.35 day −1 · ℃ and the slope of the increase (β 2 ) was highest for cultures acclimated to supraoptimal conditions (0.50 day −1 · ℃ for 26.5℃ at 72 h; Table 1).
Notably, we observed that T mort (breakpoint values) was dependent on thermal history (Table S2). Cells grown at optimal conditions of 20℃ had a significantly colder breakpoint of 27.1℃ (95% CI: 26.7, 27.6) in comparison to those grown at sub-or supra-optimal conditions with T mort of 29.0℃ (95% CI: 28.8, 29.1) and 30.4℃ (95% CI: 29.9, 30.9), respectively. However, unlike cells grown at 10℃ whereby T mort was observed after 24 h, death rates were entirely independent from temperature cultures grown at 26.5℃ as a T mort was not observed until 72 h, and hence T mort were dependent on both thermal history and exposure times. Exposure times also altered death rates when exposed to temperatures beyond the breakpoint,

F I G U R E 2
Temperature-dependent death and net population growth rates (closed and open symbols, respectively) of Phaeodactylum tricornutum acclimated to (a) sub-optimal (10℃), (b) optimal (20℃) and (c) supra-optimal (26.5℃) temperatures. Death rates were calculated following 72 h of exposure. Symbols depict the mean rate of n = 3 biological replicates and are plotted along with their associated standard deviation. Solid lines represent nonlinear least-squares fits (using the Levenberg-Marquardt method) modelled to growth rate data using Equation (8) and death rate data using Equation (9). Broken lines correspond to the 95% confidence intervals of Equations (8) and (9), estimated by parametric bootstrapping. Death and growth rates after 24 and 48 h of exposure are depicted alongside 72 h in Figure S2. Note that bootstrapping was not carried out on growth rate data obtained from supra-optimal cultures due to low number of sampling points (n = 12) [Colour figure can be viewed at wileyonlinelibrary.com] whereby temperature-dependent mortality rates increased with increasing duration. For example, in cultures grown at 10 and 20℃, β 2 was 4-and 10-fold higher, respectively, after 72 h in comparison to 24 h (Table 1). Interestingly and as mentioned above, β 2 was only significant after 72-h meaning acclimation to supra-optimal conditions delayed the onset of death at supra-optimal temperatures.
Finally, the T max for growth appeared to correlate relatively well with T mort , the thermal threshold for death. For cultures preacclimated to optimal conditions, the T mort was well predicted from the T max (i.e. close to the 1:1 line; Figure 3) and was not significantly affected by the duration of exposure (i.e. confidence intervals overlap for 24-and 72-h exposure). Similarly, a strong relationship between T mort and T max for sub-optimal cultures was evident only following 72 h of exposure (95% CI includes 1:1 line; Figure 3). In contrast to sub-optimal and optimal cultures, following 24 h of exposure, we observed no T mort in cultures pre-acclimated to supraoptimal temperatures within the range that we examined, meaning that during brief exposure to heat stress T mort was greater than 32.3℃ (represented by open diamond symbol in Figure 3). With longer exposure to heat stress (72 h), a T mort was detected, but of all treatments was the furthest from the 1:1 line and hence showed the greatest disparity between T mort and T max . Indeed unlike cultures acclimated to sub-optimal and optimal temperatures, T max in supraoptimal cultures underestimated T mort as growth ceased at 27.6℃ but cells remained viable until temperatures exceeded 30.4℃ (95% CI: 29.9, 30.9).

| DISCUSS ION
To predict changes in phytoplankton populations, we must not only understand the conditions that favour growth but also those that contribute to cell death, as these losses are an important component of closing carbon fluxes in aquatic systems (Agustí et al., 1998). In addition to affecting the carbon cycle, the population dynamics of phytoplankton affects carbon available to higher trophic levels as well as phytoplankton blooms, community structure and succession.
Cell death in response to different stress stimuli including darkness and nutrient starvation (Berges & Falkowski, 1998), heat stress (Bouchard & Yamasaki, 2008;Samuels et al., 2021), as well as many others (Zuppini et al., 2007) is often reported in terms of the proportion of dead cells relative to the total population. By contrast in this study, we used viability staining to quantify rates of cell death in response to temperature and in doing so provide data that support the assumption made in many pelagic ecosystem models that mortality rates are independent of temperature. Only when exposed to lethal temperatures, that is, temperatures outside the thermal niche, did TA B L E 1 Parameters of the linear and broken-stick models (Equation 9) for temperature-dependent death rates for cultures acclimated to sub-optimal (10℃), optimal (20℃) and supra-optimal (26.5℃) temperatures following differing exposure times (24, 48 and 72 h) Where a significant difference in slopes was found, only the second equation of the broken-stick model is shown.
we observe significant cell death. Furthermore, the baseline death rate of exponentially growing populations was practically zero, consistent with previous observations across a wide range of species (Timmermans et al., 2007).
The thermal stress response is embedded in two dimensionsintensity and duration (Rezende et al., 2014) (Gibor, 1956). Similarly in Chlorella saccharophila, mortality was evident following 2-h exposure to 44℃ (Zuppini et al., 2007), a species that does not grow above 20℃ (Vona et al., 2004) (Baker et al., 2018). Likewise in the symbiotic dinoflagellate, Symbiodinium spp., Dunn et al. (2004) also demonstrated that mortality is a consequence of intensity and duration. Hence, as was demonstrated in this study, small differences in exposure time or temperature can have significantly different outcomes for a population's survival.
Across a wide range of marine organisms, a population's survival following heat stress is not only dependent on the thermal intensity and duration but also the pre-acclimation (i.e. thermal history; Nguyen et al., 2020;Siegle et al., 2018;Stillman & Tagmount, 2009;Zhang et al., 2021). Similarly, here we found pre-acclimation temperature influenced the population response, whereby P.
tricornutum acclimated to non-optimal temperatures increased thermal tolerance as cells grown at both sub-optimal and supraoptimal temperatures had higher thermal thresholds for mortality than those grown at the temperature optimum. These findings are in keeping with work from previous authors where phytoplankton acclimated to hyperthermic conditions were characterized by traits associated with an increased tolerance to oxidative stress, for example less susceptibility to photobleaching (Takahashi et al., 2013) or lowered intracellular levels of ROS (Buerger et al., 2020).
Increases in ROS levels have been observed to immediately precede death (Franklin et al., 2004) and in some circumstances ROS has been directly linked to programmed cell death (Zhao et al., 2020). It is feasible that thermal acclimation involves mechanisms that serve to rebalance intracellular levels of ROS either by reducing production and/or by increasing scavenging. In this way, thermally acclimated cells may exhibit greater thresholds for oxidative stress, thus tolerating increased doses of temperature stress and consequently being able to survive for longer.
Studies such as ours which aim to understand the nature of cell death in response to stress stimuli and quantify the significance of cell death as a loss process are important for closing carbon budgets in aquatic systems and determining carbon transfer to higher trophic levels. We were able to provide one of the first datasets that quantify the temperature dependence of death rates. But it is likely that the temperature threshold for mortality will be speciesspecific and linked to their upper thermal limits with thermally tolerant species having higher thresholds than thermally sensitive species. Future efforts should aim to quantify death rates beyond the upper thermal limits across a number of key phytoplankton.
Similarly, studies that examine the capacity for cellular repair and regrowth following heat stress would be beneficial for predicting local extinction events.

ACK N OWLED G EM ENTS
The authors would like to thank Tania Cresswell-Maynard for rearing and maintenance of the University of Essex phytoplankton culture collection. Kirralee G. Baker was supported by NERC grant NE/ P002374/1 awarded to Richard J. Geider.

F I G U R E 3
Relationship between the critical maximum temperature (T max ) of net growth and the thermal threshold of mortality (T mort ) for Phaeodactylum tricornutum acclimated to suboptimal (10℃; circles), optimal (20℃; squares) and supra-optimal (26.5℃; diamonds) temperatures following 24 and 72 h of exposure (open and closed symbols, respectively). The solid line indicates the 1:1 relationship. Error is the 95% confidence intervals of equations (8) and (9), estimated by parametric bootstrapping (Tables S1 and  S2). Note no breakpoints (T mort ) were detected following 48 h of exposure and are not shown