The shape of a defense-growth trade-off governs seasonal trait dynamics in natural phytoplankton.

Theory predicts that trade-offs, quantifying costs of functional trait adjustments, crucially affect community trait adaptation to altered environmental conditions, but empirical verification is scarce. We evaluated trait dynamics (antipredator defense, maximum growth rate, and phosphate affinity) of a lake phytoplankton community in a seasonally changing environment, using literature trait data and 21 years of species-resolved high-frequency biomass measurements. The trait data indicated a concave defense-growth trade-off, promoting fast-growing species with intermediate defense. With seasonally increasing grazing pressure, the community shifted toward higher defense levels at the cost of lower growth rates along the trade-off curve, while phosphate affinity explained some deviations from it. We discuss how low fitness differences of species, inferred from model simulations, in concert with stabilizing mechanisms, e.g., arising from further trait dimensions, may lead to the observed phytoplankton diversity. In conclusion, quantifying trade-offs is key for predictions of community trait adaptation and biodiversity under environmental change.

with decreased grazing pressure, the summer phytoplankton bloom starts. (v) Summer is marked by severe nutrient depletion leading to strong competition within the phytoplankton community, and the relevance of different zooplankton groups (ciliates, rotifers, cladocerans and copepods) with different feeding strategies and grazing on different groups of phytoplankton. (vi) An increase of the mixing depth as autumn begins leads to a minor reduction of algal biomass and replenishing of nutrients from deeper water. The increase in nutrients may give rise to an autumn phytoplankton and crustacean maximum, paralleled by shifts in algal species composition. (vii) Early winter starts in mid of November and is characterized by an increasing intensity of deep mixing and low irradiance.

Standardized time
The duration of the seasonal phases varied among years. To account for this meteorological year-to-year variation, we aligned the sampling dates to a standardized time axis. First, each sampling date (e.g., day 25 in phase 2 in 1986) was scaled relative to the duration of the respective phase in that year (e.g. phase 2 lasts 50 days in 1986) resulting in the relative sampling day (e.g. 25/50). Multiplying the relative sampling day with the inter-annual mean duration of the respective phase (e.g. 46 days) yields the standardized day number of that sampling date (e.g. 25/50*46=23). Based on this method each sampling date can be assigned to a certain week in a standardized year (data see https://doi.org/10.6084/m9.figshare.11830464.v1). To display the seasonal biomass dynamics (shown in Fig. 2), we took the inter-annual median and quartiles of the biomass data for every standardized week and then smoothed the data by averaging the medians/quartiles of two adjacent weeks (moving average).

Aggregation of species into morphotypes
We used an intermediate level of taxonomic resolution distinguishing 36 morphotypes to achieve taxonomic consistency across the long sampling period. Each of these morphotypes contributed at least 5% to the biovolume of total phytoplankton at an individual sampling date during 1979-1982, i.e. the information about rare species got lost for the years 1979-1982. Considering the years 1979-1999, these 36 morphotypes comprise about 92% of total phytoplankton biomass. We omitted counts of morphotypes which could either not be identified, were very rare or encountered only during individual sampling events or short periods. Given its improved reliability for long-term studies we used this dataset in previous studies as well (1,(4)(5)(6). Details can also be found in LakeBase (https://fred.igb-berlin.de/Lakebase).

Mean relative biomasses
To evaluate the relative importance of a phytoplankton morphotype over the 21 years of sampling, we derived its mean annual relative biomass as follows: First, we calculated the relative biomass of each morphotype for every sampling date. Second, we averaged these relative biomasses among all dates within each year which yields the corresponding annual relative biomass of each morphotype. Finally, we derived the mean of these annual relative biomasses across years. This procedure reduces the influence of outliers at single dates, and gives equal weight to all sampling dates per year and all years, which partly differed in their total biomass and sampling resolution. The relative importance of each morphotype during distinct seasonal phases (e.g., early spring) were computed accordingly by considering only the relative biomasses of the dates within that phase (data see https://doi.org/10.6084/m9.figshare.11830464.v1). The calculated mean relative biomasses allowed to infer the respective biomass-trait distributions since each morphotype represented a specific trait combination.

C:P ratio
We used the cellular carbon to phosphate mass ratio of phytoplankton (C:P) as an indicator for nutrient depletion which was measured at the standard sampling site in 1995 (data see https://doi.org/10.6084/m9.figshare.11830464.v1) (7). The cellular C:P is more informative than the ambient phosphorous concentration in the water as phytoplankton can store phosphorous. Furthermore, most phytoplankton species can take up substantial amounts of phosphorous even at concentrations below the detection limit as they prevail in Lake Constance throughout summer.

Vertical mixing intensity
Given the large depth of the lake, most biological activity and thus sampling effort was concentrated on the upper 0-20 m depth. Thus, we report here the mean biomasses averaged across this water layer. As the surface phytoplankton concentration is usually much higher than in deep strata, deep vertical mixing (i.e. down to 60 or 100m depth) implies a net export from the surface to larger depths. To quantify this phytoplankton export, the vertical mixing intensity was inferred from a one-dimensional hydrodynamic K -ε turbulent exchange model (8,9) and expressed as net exchange rate from the uppermost layer  Appendix 2: Trait data Tab. S1: Morphotype number and name, its assigned trait values of defense δ, maximum growth rate r ( −1 ), phosphate affinity ( −1 µ −1 ) and cell volume (µ 3 ) according to Bruggeman(10)  The concave trade-off between δ and r can be also found when including all trait data from Bruggeman(10)

Appendix 3: Supporting results
We describe supporting results here, (1) correlations between phosphate affinity P and defence δ resp. maximum growth rate r, (2) the relationship of all three traits to cell volume, (3) the seasonal biomass-trait distribution considering all seven phases, the influence of re-oligotrophication on (4) the seasonal dynamics of r and δ, and (5) on the seasonal dynamics of P.

Correlations of defense and maximum growth rate with phosphate affinity
A concave trade-off curve was most obvious for the trade-off between δ and r (Fig. 3). In the trait space of δ and P, a similar pattern might be seen albeit with more scatter and one exception, Ceratium hirundinella being very defended and highly phosphate affine (Fig. S2a, ρ = -0.28, p = 0.09, Spearman rank correlation coefficient never biomass-weighted). The pattern for P and r is even more scattered (Fig. S2b, ρ = -0.12, p = 0.5).

Trait correlations with cell volume
The correlation between each trait and cell volume as a master trait was tested (Fig. S3). Maximum growth rate r was negatively correlated to cell volume (ρ = -0.59, p = 10⁻ 4 ), defense was positively correlated to cell volume (ρ = 0.49, p = 10⁻ 3 ), while we found no correlation for phosphate affinity (ρ = 0.22, p = 0.19). Thus, the trade-off between δ and r may partially arise from the weak correlations of both traits with cell size.
However, the large scatter in the relationship between δ and cell size (Fig. S3b) shows that other defense mechanism are important as well, e.g. cell shape, formation of filaments. Hence we used trait data instead of the approximation cell size.
Furthermore, cell size as a functional trait is harder to link directly to a certain environmental driver being sensitive to multiple factors as a 'master trait', compared to, e.g. defense being selected by high grazing pressure, or phosphate affinity being selected by phosphate depletion. Therefore, we do not include cell size in our main consideration of how the community trait composition responds to seasonal environmental changes, but we use cell size more as a 'master trait' providing to some extent a potential mechanism for the observed trade-off between defense and maximum growth rate.

Biomass-trait distributions for all seasonal phases
The biomass distribution in the δ-r trait space responded to the seasonally changing environment in a remarkably gradual and consistent way (Fig. 2, Fig. S4). In late winter and early spring, vertical mixing and the resulting high export of phytoplankton from the euphotic zone to deep water layers was a dominant driver of the phytoplankton community in deep Lake Constance while grazing pressure and nutrient depletion were very low. Morphotypes with high r being able to compensate for high losses and to exploit the high nutrient concentrations dominated, whereas morphotypes with low r and high were almost absent (Fig. S4a, b). This is reflected in the community average trait values (late winter: ̅ = 0.51, ̅ = 1.56 −1 ; early spring: ̅ = 0.52, ̅ = 1.57 −1 ). Morphotypes with low or high phosphate affinities had high biomasses, indicating the absence of a selection pressure on this trait. During late spring, grazing pressure increased mostly by ciliates ( Fig. 2) but did not initiate a shift of the overall biomass distribution towards higher ( ̅ = 0.48, ̅ = 1.55 −1 ) (Fig. S4c). During the clear-water phase (CWP), the grazing pressure was at its annual maximum (Fig. 2).
The community average maximum growth rate decreased slightly ( ̅ = 1.35 −1 ) while the mean defense level did not change ( ̅ = 0.48) despite the high grazing pressure (Fig. S4d), probably due to a delayed numerical response of highly defended but slowly growing morphotypes. In summer, nutrient depletion and grazing pressure were the dominant drivers of phytoplankton (Fig. 2).

Influence of re-oligotrophication on seasonal dynamics of defense and maximum growth rate
Concurrent to the re-oligotrophication of the lake during the study period, the pattern of the seasonal biomasstrait distribution only marginally changed (Fig. S5). In early spring in 1979 -1988, the morphotypes were on average slightly less defended than on average in 1989-1999 due to Rhodomonas spp. being more abundant, while intermediately defended diatoms were less common. These changes did not alter the average maximum growth rate. In summer in 1979-1988, the average maximum growth rate was higher than in 1989-1999, e.g. due to a higher share of the relative fast-growing Fragilaria crotonensis, whereas in 1989-1999 slow-growing Ceratium hirundinella was relatively more abundant. These changes did not affect the average level of defense. To conclude, the changes in the species composition observed during the reoligotrophication did not change the overall seasonal pattern with a dominance of fast growing, undefended morphotypes in early spring and of slow growing, highly defended morphotypes in summer.
Our explanation for the lack of a clear responsiveness of the defence-growth rate trait distribution is that the grazing pressure likely changed little during most of the investigation period lasting from 1979-1999. We know that from 1987-1998 neither total ciliate biomass nor species composition changed significantly (11). The crustaceans, i.e. the other important group of herbivores, had lower abundances in 1997-1998 than in 1979-1996 during July to September. These fairly constant biomasses of herbivores fit with the measurements of 14C primary production from 1980-1996. During this time primary production declined (only) during summer by only 25% (12). We presume that potential effects on the trait distribution are too small to be clearly visible.
On the other hand, the persistence of the trait distribution and its seasonal dynamics suggest that trophic interactions played a major role in this lake during the whole study period, as supported by numerous other studies(13). Influence of re-oligotrophication on seasonal dynamics of phosphate affinity In summer, SRP dropped always below 2 µg P/l in the surface layer and we measured a strong increase in the cellular phytoplankton C:P ratio in summer 1995. Nevertheless we found no distinct seasonal pattern in community average phosphate affinity when considering all years together or the years 1979-1988 (Fig.   S6a,b). For the years 1989-1999 we see a slight increase in phosphate affinity during the season, which results in somewhat higher values in summer and autumn compared to 1979-1988. The overall changes of the community average trait values are small relative to the entire trait range (3 -1600 −1 µ −1 ).

Appendix 4: Model description
The model The following model equations describe the biomass dynamics of phytoplankton species and one zooplankton group : The phytoplankton growth is limited by nutrients, described by a Monod term. By assuming a fixed pool of nutrients, we can write the available (dissolved) nutrient concentration as = − ∑ =1 − 1 , that is, the total amount of nutrients subtracted by the nutrients fixed in biomass of phytoplankton and zooplankton (14). The nutrient concentration is written in units of phytoplankton biomass, i.e. represents the maximum phytoplankton biomass obtainable from the nutrient pool in the absence of mortality. The phytoplankton species differ in their maximum growth rates , but share the same half-saturation constant for nutrient uptake (in units of phytoplankton biomass, see above) and natural mortality . Hence, the species with the highest , performing well at high resource availability, is also the superior competitor under strong resource depletion (i.e. it has the lowest * ) in the model. The grazing of zooplankton on phytoplankton is described by a Holling type II function with the maximum grazing rate and the half-saturation constant . Phytoplankton species have different values of defense against zooplankton. We assume that defended phytoplankton cells also demand handling time of the predator equal to that of undefended phytoplankton but without energy gain because unselective feeders, which dominate in Lake Constance, are probably not able to discriminate between them and attack both (15). Accordingly, gives the probability of not being consumed (i.e., not ingested or digested) and surviving when attacked with values ranging between 0 (undefended) and 1 (completely defended). The probability of being consumed is then given by 1 − , which scales the maximum grazing rate (see Eq. S1) and corresponds to the 'edibility', typically used in a limnetic context (10).
The conversion efficiency of consumed phytoplankton into zooplankton biomass is assumed to be equal among the phytoplankton species. represents the zooplankton mortality.

Trade-off curve
The trade-off curve between defense and maximum growth rate is given by the function where denotes the shape parameter, the slope parameter and the maximum growth rate of the most defended species ( = 0.9). If < 1, the trade-off curve is concave. > 1 gives a convex trade-off curve and = 1 a linear one. We assume a concave trade-off curve with = 0.2, = 1.6 −1 and = 0.5 −1 approximately reflecting the trade-off curve found in the trait data (Fig. 3a). For comparison, we consider also a convex trade-off curve with = 2, = 1.92 −1 and = 0.5 −1 , that crosses the concave trade-off curve at minimal and maximal defense levels (i.e. shares the same endpoints).

Parametrization and initialization
We considered different phytoplankton species with trait values spanning the whole feasible trait space. We Based on measurements conducted at Lake Constance (13,16,17), we parametrized the model as follows:

Numerical integration
The numerical integrations of the model were done with the ode45 solver of the deSolve package in R(18).
We run the simulations for 10,000 days and calculated the mean biomasses of the last 1000 days to detect the phytoplankton species dominating in the long term. Furthermore, we checked which species survive in the short term, that is, within the first 100 days. The extinction threshold was set to 10 −4 −3 . We