Stoichiometric mismatch causes a warming-induced regime shift in experimental plankton communities.

In many ecosystems, consumers respond to warming differently than their resources, sometimes leading to temporal mismatches between seasonal maxima in consumer demand and resource availability. A potentially equally pervasive, but less acknowledged threat to the temporal coherence of consumer-resource interactions is mismatch in food quality. Many plant and algal communities respond to warming with shifts toward more carbon-rich species and growth forms, thereby diluting essential elements in their biomass and intensifying the stoichiometric mismatch with herbivore nutrient requirements. Here we report on a mesocosm experiment on the spring succession of an assembled plankton community in which we manipulated temperature (ambient vs. +3.6°C) and presence versus absence of two types of grazers (ciliates and Daphnia), and where warming caused a dramatic regime shift that coincided with extreme stoichiometric mismatch. At ambient temperatures, a typical spring succession developed, where a moderate bloom of nutritionally adequate phytoplankton was grazed down to a clear-water phase by a developing Daphnia population. While warming accelerated initial Daphnia population growth, it speeded up algal growth rates even more, triggering a massive phytoplankton bloom of poor food quality. Consistent with the predictions of a stoichiometric producer-grazer model, accelerated phytoplankton growth promoted the emergence of an alternative system attractor, where the extremely low phosphorus content of the abundant algal food eventually drove Daphnia to extinction. Where present, ciliates slowed down the phytoplankton bloom and the deterioration of its nutritional value, but this only delayed the regime shift. Eventually, phytoplankton also grew out of grazer control in the presence of ciliates, and the Daphnia population crashed. To our knowledge, the experiment is the first empirical demonstration of the "paradox of energy enrichment" (grazer starvation in an abundance of energy-rich but nutritionally imbalanced food) in a multispecies phytoplankton community. More generally, our results support the notion that warming can exacerbate the stoichiometric mismatch at the plant-herbivore interface and limit energy transfer to higher trophic levels.

INTRODUCTION

The earliest and most pervasive evidence of biological responses to contemporary climate change comes from advances in the timing of seasonal events such as budburst and flowering of plants, the onset of breeding activities of various animals, and plankton spring blooms (Parmesan, 2006; Thackeray et al., 2010). In many ecosystems, consumers respond to warming differently than their resources (Kharouba et al., 2018), sometimes leading to temporal mismatches between seasonal maxima in consumer demand and resource availability (Durant et al., 2007; Edwards & Richardson, 2004; Kerby & Post, 2013; Thomas et al., 2001). A potentially equally pervasive, but less acknowledged threat to the temporal coherence of consumer‐resource interactions is mismatch in food quality. Most herbivores have higher body contents of essential nutrients such as phosphorus and nitrogen than their foods, leading to stoichiometric imbalances between herbivore needs and plant nutrient content (Elser et al., 2000; Sterner & Hessen, 1994). Many plant and algal communities respond to warming with increased carbon fixation (Toseland et al., 2013; Zvereva & Koslov, 2006), accelerated senescence (Doiron et al., 2014), or shifts toward more carbon‐rich, less palatable species and growth forms (Craine et al., 2010; Klein et al., 2007; Woods et al., 2003; Yvon‐Durocher et al., 2015, 2017). These responses dilute the concentration of essential elements in the forage for herbivores, potentially intensifying the stoichiometric mismatch with herbivore nutrient requirements (Craine et al., 2010; Doiron et al., 2015; Urabe & Sterner, 1996).

One of the most conspicuous seasonal events on Earth is the phytoplankton spring bloom in lakes and oceans of temperate to polar regions (Siegel et al., 2002; Sommer et al., 2012). Its timing is linked to the life cycles of commercially important fish species (Cushing, 1990), and its magnitude is a main determinant of the annual export of photosynthetically fixed carbon to deep waters and sediments (Turner, 2015). Phytoplankton blooms are demographic phenomena resulting from the balance of births and deaths and are therefore intimately linked to the population dynamics of planktonic grazers. It has been proposed that phytoplankton spring blooms occur predominantly in cold waters because, at low temperatures, the development rates of metazoan grazers are slow relative to phytoplankton growth rates, creating a “loophole” where phytoplankton can temporarily outgrow their grazers when light levels seasonally increase (Irigoien et al., 2005; Rose & Caron, 2007). Because the population growth rates of heterotrophs increase more steeply with temperature than the growth rates of autotrophs (Lopez‐Urrutia et al., 2006), grazers frequently catch up with their phytoplankton prey as water temperatures rise. A phytoplankton spring bloom is therefore often followed by an increase in planktonic grazer abundance and the transition to a clear‐water phase (Sommer, 2009; Sommer et al., 2012).

In temperate freshwater lakes, the termination of the phytoplankton spring bloom has traditionally been attributed to grazing by crustaceans of the genus Daphnia (Sterner, 1989). In short, because Daphnia growth and development is insignificant in winter but increases steeply with temperature above 8°C (Vijverberg, 1980), both the spring peak in Daphnia abundance and the subsequent clear‐water phase are tightly synchronized with seasonal warming across Northern Hemisphere lakes (Straile, 2002; Straile et al., 2012). More recently, it has been recognized that microzooplankton, in particular ciliates, can also be major grazers in freshwater systems (Sommer et al., 2012). Because many ciliates are less sensitive to low temperatures than crustaceans, they can cause substantial grazing losses in early phases of a phytoplankton spring bloom (Weisse, 2006). Yet ciliates typically cannot fully suppress a developing bloom but rather track it demographically or alter phytoplankton community size structure (Berger et al., 2010, 2014; Peeters et al., 2007; Sherr & Sherr, 2009). This raises the question of the extent to which ciliates and warming interact in shaping the timing and magnitude of seasonal events during the spring succession of plankton and how these factors play out in terms of mismatches with the energetic and nutritional needs of crustacean zooplankton such as Daphnia.

Here we report on a combined experimental and theoretical study testing the following two pairs of contrasting hypotheses. First, because the growth rates of all planktonic organisms respond positively to warming at low late‐winter temperatures, warming accelerates all phases of the spring succession of the plankton, that is, warming causes earlier abundance peaks of first phytoplankton and then grazers, accompanied by an earlier onset of the clear‐water phase (Berger et al., 2010; Straile, 2002; Straile et al., 2012) (Hypothesis A1). Alternatively, warming promotes a build‐up of carbon‐rich phytoplankton biomass of low food quality (Yvon‐Durocher et al., 2017), potentially slowing down the growth of grazers and delaying the clear‐water phase (Hypothesis A2). Second, since ciliates become active at lower temperatures than crustaceans (Weisse, 2006), grazing by ciliates during early bloom phases should slow down phytoplankton growth and, thus, the deterioration of its nutritional value, promoting the subsequent growth of Daphnia and advancing the clear‐water phase (Jäger et al., 2008) (Hypothesis B1). Alternatively, selective ciliate grazing may shift the phytoplankton community toward more grazing‐resistant forms (Tirok & Gaedke, 2006), reducing Daphnia population growth and delaying the initiation of a clear‐water phase (Hypothesis B2). We addressed these hypotheses with a mesocosm experiment on the spring succession of the plankton in which we manipulated temperature (ambient vs. warmed), presence versus absence of ciliates, and presence versus absence of Daphnia in a fully factorial design. The experiment supported Hypotheses A2 and B1, but the food quality response to warming was much stronger than expected and caused a regime shift. To gain insight into the underlying mechanisms of the regime shift, we therefore supplemented the experiment with numerical simulations of a stoichiometrically explicit plankton community model that was tailored to the experimental conditions.

METHODS

Experimental design and setup

We performed the experiment in a fenced outdoor mesocosm facility at the Planegg‐Martinsried campus of Ludwig Maximilian University of Munich from late winter to late spring 2008. We exposed four different plankton communities (a phytoplankton community in the presence vs. absence of ciliates and Daphnia in all possible combinations) to two temperature treatments (ambient vs. warmed by 3.6°C) in a fully factorial design. The eight treatments were subsequently labeled “aN” (ambient, no grazers), “aC” (ambient, ciliates), “aD” (ambient, Daphnia), “aCD” (ambient, ciliates and Daphnia), “wN” (warmed, no grazers), “wC” (warmed, ciliates), “wD” (warmed, Daphnia), “wCD” (warmed, ciliates and Daphnia). Three replicates of each treatment yielded a total of 24 mesocosms.

The mesocosms consisted of cylindrical bags (diameter 0.95 m) made from transparent Tricoron (RKW AG, Wasserburg, Germany), were open to the atmosphere and heat‐sealed in the bottom at a depth of 1.5 m, enclosing a water volume of ~1 m3. Mesocosms were made optically deep by surrounding them with black silage film. The resulting background attenuation coefficient (K bg) for photosynthetically active radiation (PAR) was approximately 1.0 m−1. It was determined from vertical PAR profiles with a spherical sensor (LI‐139SA, Licor, Lincoln, Nebraska, USA) in the absence of phytoplankton. The underwater light environment in the mesocosms thus corresponded to that in a 6‐m‐deep water column of a typical clear‐water lake of the region with K bg = 0.25 m−1. We covered each mesocosm with a transparent lid mounted approximately 2 cm above the mesocosm frames to allow gas exchange with the atmosphere.

The mesocosms were suspended in two adjacent, unshaded concrete pools that acted as water baths for the two temperature treatments. Both pools had an area of 50 m2 and were brushed clean and filled with well water prior to the experiment. The water temperature in the ambient pool was not manipulated, whereas the warmed pool contained an industrial heating element (ISA‐Heinrich‐Industrietechnik, Falkensee, Germany). An electronic circuit tightly regulated the temperature in the warmed pool at +3.6°C above the ambient pool. Water in each pool was continuously mixed by pumping compressed air to the bottom of the pool at six locations and circulating the water horizontally with a submersible water pump.

The experiment started on 22 February (day 1 of the experiment) and ran until 27 May (day 96). On day 1, all the mesocosms were filled with 5 μm filtered well water containing high concentrations of macronutrients except for phosphorus. To experimentally simulate the nutrient pulse that typically occurs during spring overturn, we enriched all the mesocosms with KH2PO4 to 28 mg P m−3 at the beginning of the experiment.

Because microzooplankton cannot be selectively removed from natural plankton communities without also removing most phytoplankton, we created the experimental communities by repeatedly stocking all the mesocosms with small amounts of precultured algae, Daphnia, or ciliates. The stocking regime was intended to mimic the transient conditions at the onset of spring succession, when population densities can change rapidly from initially low late‐winter densities. On days 1, 3, 11, and 21, we inoculated the chlorophytes Monoraphidium minutum and Scenedesmus obliquus, the tebouxiophyte Choricystis minor, the baccilariophytes Nitzschia palea and Cyclotella sp., and the cryptophyte Cryptomonas phaseolus at a cumulative density of 1.3 mg chl a m−3 (Appendix S1: Table S1), which is within the range of late‐winter chl a values in lakes of the region. These algal species cover major taxonomic groups and are edible to zooplankton. On days 3, 11, and 21, we set up the microzooplankton treatments by inoculating the aC, aCD, wC and wCD mesocosms with the ciliate Urotricha furcata from laboratory cultures at densities of 4.2 × 106, 3.4 × 106, and 1.2 × 106 individuals m−3, respectively. This species is very widespread, and the cumulative stocking density was within the range of typical late‐winter densities of this genus. The aN, aD, wN, and wD mesocosms received no ciliates. On day 4, all the aD, aCD, wD, and wCD mesocosms received Daphnia hyalina from a stock culture at a density of 2.1 × 103 individuals m−3, simulating the overwintering Daphnia population in Lake Brunnsee, from which the stocked population had been originally isolated. We deliberately stocked Daphnia at a relatively high density because, in our experience, mortality in the first week after stocking exceeds 50%. Prior to stocking, all organisms were acclimated from room temperature to ambient outdoor water temperature in a stepwise manner over a period of 1 week.

Sampling and laboratory analyses

Water temperature was logged every 30 min at four locations in each pool and in two mesocosms of each treatment (Voltcraft K 204 Datalogger, Conrad Elektronik, Germany). We measured phytoplankton biomass twice per week as chl a concentration and determined its taxonomic composition and biovolume every other week. The abundances of micro‐ and mesozooplankton and the concentrations of soluble reactive phosphorus (SRP), particulate phosphorus (PP), and particulate organic carbon (POC) were measured weekly.

Samples for phytoplankton, water, and seston chemistry were taken after mixing each mesocosm with a Secchi disk. Chlorophyll a was measured on 250 μm filtered water by in vivo fluorescence (TD 700, Turner Designs, Sunnyvale, California, USA) within a few hours after sampling. The concentrations of SRP and PP were measured using standard methods (Wetzel & Likens, 1991). POC was filtered on pre‐combusted glass fiber filters (GF6, Schleicher & Schüll, Germany) and determined with a C/N analyzer (EA 1110, C3‐Analysetechnik GmbH, Munich, Germany). The taxonomic composition of the phytoplankton was assessed by microscopic counts of Lugol‐preserved samples in an inverted microscope. A minimum of 600 individual units of the most abundant taxa were counted. We measured cell or colony size for the determination of phytoplankton biovolume on digital images with the software ImageJ 1.40 g (National Institutes of Health, Bethesda, Maryland, USA). Phytoplankton biovolume estimates were based on approximations of cell shapes by simple geometrical bodies (Hillebrand et al., 1999).

Mesozooplankton (almost exclusively D. hyalina) were sampled with a 1.5 m vertical haul of a 105‐μm mesh net (20‐cm opening diameter), preserved in 4% formalin–25% sugar solution, and subsequently counted in a dissecting scope. Ciliates were sampled with a 1.5‐m vertical haul of a 10‐μm mesh net (10‐cm opening diameter), preserved in Lugol's solution, and subsequently counted in settling chambers in an inverted microscope.

Treatment effectiveness

The warmed mesocosms had a 3.6°C higher water temperature than the ambient mesocosms throughout the experiment (Figure 1a). The ciliate treatments were also maintained throughout the experiment, that is, no ciliates were found in mesocosms in which we had not stocked U. furcata. The inoculation with Daphnia failed in one replicate of the ambient, Daphnia treatment, which was therefore excluded from all analyses. Four phytoplankton taxa spontaneously colonized the mesocosms and became abundant in some treatments (Figure 2, Appendix S1: Table S1). All of these taxa made their first appearances around the same time in all treatments.

FIGURE 1

Time course of (a) soluble reactive phosphorus (SRP), (b) phytoplankton as chl a, (c) total phosphorus (TP), (d) the seston carbon:phosphorus ratio (C:P), (e) Daphnia abundance, and (f) ciliate abundance. Shown are mean values of the eight experimental treatments as indicated in the legend. N = 3 for all treatments except for ambient, Daphnia, where N = 2 (see Appendix S1: Figure S1 for separate plots of all replicates). Dotted lines in panel (a) show water temperatures in the ambient and warmed treatments. The results of 3‐way repeated measures analyses of variance (ANOVAs) of the time‐dependent effects of temperature, ciliate, and Daphnia treatments on the variables in panels (a–f) are summarized in Appendix S1: Tables S2–S7. The horizontal broken line in panel (b) indicates the threshold for a clear water phase

FIGURE 2

Taxonomic composition of the phytoplankton community over the course of experiment. Panels (a–h) show the mean percentage contribution of 12 different taxa to total phytoplankton biovolume (N = 3 replicate mesocosms, except for aD, where N = 2). Six of the taxa were initially stocked (Choricystis minor, Monoraphidium minutum, Scenedesmus obliquus, Cryptomonas phaseolus, Cyclotella sp., Nitzschia palea), and the remaining taxa colonized spontaneously (Closteriopsis acicularis, Achnanthes coarctata, Synedra acus, Pseudanabaena sp., an unidentified flagellate, and an unidentified coccoid green alga). Horizontal bars above each panel mark the bloom period (when phytoplankton biomass was >2 mg chl a m−3). As an indication of within‐treatment variance, white numbers at the top of each panel show the percentage community overlap among the replicates of each treatment (calculated with Equation 1), where values >80% indicate low variance

Statistical analyses

All statistical analyses were performed in SPSS version 25 and are reported in detail in Appendix S1: Tables S2–S11. Data on the temporal trajectories of state variables were log‐transformed and analyzed using 3‐way repeated measures analysis of variance (ANOVA), using temperature (ambient vs. warmed) and the presence/absence of ciliates and Daphnia, respectively, as treatment factors and sampling date as the within‐subject effect. We used 3‐way ANOVA to test for treatment effects on the timing of the clear‐water phase and on the timing and magnitude of the peaks in phytoplankton biomass and the abundances of Daphnia and ciliates. We defined peak magnitude as the maximum value of chl a, Daphnia, and ciliate abundance observed in each replicate and peak timing as the day on which a maximum was observed. Timing of the clear‐water phase was defined as the first day following the spring bloom on which chl a concentration decreased below 2 mg m−3.

Treatment effects on phytoplankton community composition were investigated in two ways. First, we used 3‐way repeated measures ANOVA to test for treatment effects on the proportional contribution of taxa >30 μm to total phytoplankton biovolume. Taxa with a longest axial dimension that exceeds 30 μm are difficult to handle by most planktonic grazers and are often categorized as grazing resistant (Bell, 2002; Sprules & Knoechel, 1984). Second, we assessed the similarity in phytoplankton community composition within and between treatments. We defined similarity as the percentage overlap in (biovolume‐based) community composition between pairs of mesocosms, calculated as where p and p are the proportional contributions of taxon k to total phytoplankton biovolume in mesocosms i and j, respectively. For each sampling date, we calculated percentage overlap in community composition between all pairs of mesocosms according to the preceding formula. We then calculated the mean percentage overlap within a given treatment by averaging the pairwise overlap values between the three replicates of that treatment. We also calculated the mean percentage overlap between any two different treatments by averaging the pairwise overlap values between the replicates from one treatment with each of the replicates from the other treatment. These calculations indicated that ambient and warmed communities diverged quickly over time, whereas communities within the same temperature treatment remained similar. We tested this further by calculating the mean percentage overlap within a given temperature treatment by averaging over all pairwise overlap values between mesocosms sharing the same temperature. We compared this to the mean percentage overlap between temperature treatments by averaging over all pairwise overlap values between the 12 warmed mesocosms with each of the 11 ambient mesocosms.

Model structure

The plankton community model describes two grazer populations with fixed nutrient:carbon stoichiometry feeding on a phytoplankton population with flexible nutrient:carbon stoichiometry in a well‐mixed water column of uniform depth. Phytoplankton and grazer production can be limited by energy (PAR and food carbon, respectively) or the phosphorus content of algal biomass, where the growth rates of both algae and grazers are increasing functions of the flexible phosphorus: carbon quota of phytoplankton biomass. The model describes the dynamics of seven dynamical state variables: the vertical gradient of PAR in the water column, dissolved mineral phosphorus concentration, phytoplankton biomass, the phosphorus:carbon quota of phytoplankton biomass, the stock of sedimented particulate phosphorus, and the biomasses of ciliates and Daphnia. Like the experimental mesocosms, the model system is closed for phosphorus. It harbors a fixed total amount of phosphorus (the sum of phosphorus in dissolved and particulate forms), and mineral phosphorus supply to the water column originates exclusively from an initial pulse of dissolved phosphorus (mimicking spring turnover), excretion by living organisms, and remineralization of dead biomass.

The model extends an existing stoichiometrically explicit phytoplankton–Daphnia model (Diehl, 2007) in three ways. First, microzooplankton was included as a second grazer of phytoplankton and a second food of Daphnia. Second, based on recent experimental evidence (Uszko et al., 2015, 2017), the functional response of Daphnia changed from type 2 to type 3. Third, several rate parameters were assumed to be temperature dependent, including all metabolic and loss rates, the maximum production rates of phytoplankton, and the attack rates and handling times of grazers. All model assumptions, dynamical equations, and choices of parameter values are described and motivated in detail in Appendix S2. The definitions and units of all variables, functions, and parameters are listed in Appendix S2: Tables S1 and S2.

Model analyses

We tailored our numerical analyses to the mesocosm experiment by setting all environmental parameters (incident radiation, water column depth, background light extinction, total nutrients) and the initial conditions of all state variables to values approximating the experiment. Furthermore, three phytoplankton parameters were adjusted to experimental observations. The minimum algal nutrient quota was set to the observed minimum seston P:C ratio, and the maximum algal production and sinking rates (p max and v) were set to values representative of algae in ambient and warmed mesocosms, respectively. In accordance with observations, we assumed that p max < p max and v  > v , where the suffixes a and w stand for the ambient and warmed treatments, respectively. All other phytoplankton traits were assumed to be equal between “ambient” and “warmed” communities (see Appendix S2 for details). We subsequently investigated the transient spring dynamics of ambient and warmed communities using the following three approaches:

To study the transient dynamics of ambient and warmed plankton communities during late winter to late spring, we performed numerical simulations of the full dynamical system described by Appendix S2: Equations S1–S7 (Appendix S2: Table S1). Simulations using ambient versus warmed phytoplankton traits were driven by the actual, measured temperatures in the ambient and warmed mesocosms (Figure 1a), respectively, and were started from the approximate initial densities of phytoplankton at the time of grazer inoculation (Appendix S2: Table S2). Simulations were performed for all combinations of presence versus absence of ciliates and Daphnia, respectively, to explore the full range of experimental treatments. All simulations were run in MATLAB version R2018b. Annotated MATLAB code and the driving temperatures are provided in Diehl et al. (2021b) at https://doi.org/10.5281/zenodo.4715500.

To gain analytical insights into possible attractors of the system, we performed a zero net growth isocline (ZNGI) analysis of a simplified version of the model. We simplified the model in two ways. First, we excluded micrograzers and focused exclusively on phytoplankton–Daphnia dynamics. Second, we assumed that the mineralization of sedimented particles and nutrient uptake by phytoplankton were instantaneous processes. Under these circumstances, the pool of sedimented nutrients becomes zero, and the concentration of dissolved mineral nutrients is measurably above zero only when algae are nutrient replete (i.e., the nutrient quota is at its maximum). The system then simplifies to a two‐dimensional set of differential equations describing the dynamics of phytoplankton and Daphnia, the remaining state variables being determined by algebraic equations, as described in Appendix S2. The transient and long‐term dynamics of this system can be understood from an analysis of the ZNGIs in the phytoplankton–Daphnia space. ZNGIs are combinations of phytoplankton and Daphnia densities for which the net rate of change of either population is zero. Intersections of these isoclines with the axes and with each other mark system states that are either stable or unstable equilibria. Isocline analyses were performed under different scenarios of constant environmental temperatures ranging from 6°C to 25°C.

We performed bifurcation analyses to search for thresholds in parameter values that engender qualitatively different system behaviors. To do so, we simulated the full dynamical system to equilibrium along gradients of environmental conditions (temperature, nutrient enrichment, water column depth) and algal traits (maximum growth rate, minimum nutrient quota, sinking velocity), varying one parameter at a time and starting from different initial conditions.

EXPERIMENTAL RESULTS

At ambient temperatures, a typical spring succession developed. After a short lag phase, phytoplankton grew exponentially at an average rate of 0.08/day from day 8 onward, producing a moderate spring bloom (Figure 1b, blue lines). The phytoplankton community went through a succession of diverse taxa (Figure 2a–d), most of which suffered from substantial sinking losses, as indicated by the fast drawdown of the limiting nutrient; total phosphorus, which is the sum of dissolved and particular phosphorus, declined rapidly and continuously throughout the experiment (Figure 1c, blue lines; Appendix S1: Table S4, effect of time and time × temperature interaction p < 0.001). In treatments with Daphnia, phytoplankton net growth came to a halt and biomass peaked around day 40–50 (Figures 1b and 3a, dark blue lines and bars; Appendix S1: Table S8a, temperature × Daphnia interaction p < 0.001), when water temperature exceeded 8°C for the first time (Figure 1a, dotted blue line). From approximately that day on, the population of D. hyalina, which had been declining during the earliest phase of the experiment, grew exponentially at an average rate of 0.11/day and reached a peak on days 76–83 (Figures 1e and 3c, dark blue lines and bars). During the exponential growth phase of Daphnia, the phytoplankton bloom crashed and reached a clear‐water state (<2 mg chl a m−3) between days 65 and 72, which lasted until the end of the experiment (Figures 1b and 3d, dark blue lines and bars). In the absence of Daphnia, phytoplankton biomass peaked later (day 60–65) and at higher biomass; it subsequently declined somewhat, leveled off at a moderate biomass in most replicates and reached a clear‐water state in only two out of six mesocosms without Daphnia (Figures 1b and 3a,d,e, light blue lines and bars; Appendix S1: Table S8a–c, temperature × Daphnia interactions p < 0.001).

FIGURE 3

Timing and magnitude of cardinal events during spring succession. The upper row of panels shows the timing of (a) the phytoplankton bloom peak, (b) the peak in ciliate abundance, (c) the peak in Daphnia abundance, and (d) the onset of the clear‐water phase. A clear‐water phase was never reached in any of the warmed treatments. In the ambient treatments without Daphnia, a clear‐water phase was reached in only one replicate each. These treatments are “flagged” with a bleaker, hatched coloring, to signal that the bars are based on a single replicate. The lower row of panels shows the peak biomass or peak abundance of (e) phytoplankton, (f) ciliates, and (g) Daphnia. Bars show mean values ±1 SE of the experimental treatments as indicated in the legend. N = 3 for all treatments except for ambient, Daphnia, where N = 2. (h) Peak abundance of Daphnia plotted against peak chl a concentration. Each data point represents one mesocosm (N = 11)

In contrast, the warmed mesocosms experienced a massive, long‐lasting phytoplankton bloom and the complete absence of a clear‐water phase, irrespective of the presence/absence of Daphnia (Figures 1b and 3d,e, red and orange lines and bars; Appendix S1: Tables S3 and S8c, main temperature effects and temperature × Daphnia interactions p < 0.001). From day 8 to 27, phytoplankton biomass increased on average 2.5 times faster in the warmed mesocosms (average growth rate 0.21/day) than in the ambient mesocosms (Figure 1b; Appendix S1: Table S3, time × temperature interaction p < 0.001). Within a few weeks, the bloom became dominated by a single taxon—the small, single‐celled chlorophyte M. minutum–which contributed 75%–95% to total phytoplankton biomass between days 27 and 69 and >30% thereafter, when it was partly replaced by an unidentified coccoid green alga (Figure 2e–h). The Monoraphidium‐dominated bloom experienced only minor sinking losses, as indicated by the constancy of total phosphorus from day 30 onwards (Figure 1c, red and orange lines; Appendix S1: Table S4, time × temperature interaction p < 0.001). Phytoplankton biomass continued to increase exponentially at an average rate of 0.016/day from day 36 until the end of the experiment, when all warmed treatments exceeded 140 mg chl a m−3 (Figure 1b, red and orange lines; Appendix S1: Table S3, time × temperature interaction p < 0.001).

As expected, Daphnia initially benefited from early warming to temperatures above 8°C, reaching average densities 17 times higher in the warmed than in the ambient mesocosms on day 27 (Figure 1e; Appendix S1: Table S5, time × temperature interaction p < 0.001). Yet Daphnia population growth ceased in the warmed mesocosms already between days 27–48 at peak densities, which were on average 3–10 times lower than the Daphnia peak densities subsequently reached at ambient temperatures (Figures 1e and 3c,g; Appendix S1: Table S9a,b, main temperature effects p < 0.001). Remarkably, Daphnia in the warmed mesocosms declined thereafter to densities near or below detection level even though its food was abundant and of edible size (Figure 1b,e, red and orange lines). Poor Daphnia performance under warming conditions coincided with a continuous decline in the nutritional quality of the phytoplankton. After depletion of mineral phosphorus (Figure 1a), the seston C:P ratio in the warmed mesocosms increased almost linearly throughout the experiment, suggesting that a large fraction of the phytoplankton growth was fueled by cell‐internal phosphorus stores (Figure 1d, red and orange lines; Appendix S1: Table S6, main temperature effect and time × temperature interaction p < 0.001). This produced an increasing stoichiometric mismatch between Daphnia's nutritional needs and the phosphorus content of phytoplankton. Daphnia population growth in the warmed mesocosms ceased at seston C:P ratios of 800–900 and turned negative above seston C:P ratios of 1000–1100 (Figure 1d,e, red and orange lines). Seston C:P ratios in the warmed mesocosms exceeded extreme values of 1300 on day 69 and reached up to 2200 toward the end of the experiment (Figure 1d, red and orange lines). In contrast, at ambient temperatures, seston C:P ratios stayed below 800 throughout the experiment, and Daphnia population growth came to a halt first during the clear‐water phase, when food quantity became limiting, as indicated by the chl a concentration, which dropped below 2 mg m−3 (Figures 1b,d,e and 3c,d, blue lines and bars; Appendix S1: Table S9a, main temperature effect p < 0.001).

In agreement with expectations, microzooplankton responded faster than Daphnia to early phytoplankton bloom development at low temperatures. Average growth rates during exponential growth phases were 0.2/day for the ciliate U. furcata and 0.1/day for D. hyalina (Figure 1e,f). Ciliates reduced the magnitude of the phytoplankton bloom (Figure 3e; Appendix S1: Table S8b, main ciliate effect p = 0.031) but did not suppress it (Figure 1b). As a consequence, grazing by ciliates delayed the deterioration of algal food quality (Figure 1d; Appendix S1: Table S6, main ciliate effect p = 0.002, time × ciliate interaction p = 0.027) and transiently benefited Daphnia, which reached higher peak densities in the presence than in the absence of ciliates (Figures 1e and 3g; Appendix S1: Table S9b, main ciliate effect p = 0.009). Conversely, ciliates were negatively affected by Daphnia, as indicated by the lower ciliate peak abundances in the presence of Daphnia (Figure 3f; Appendix S1: Table S10b, main Daphnia effect p < 0.001).

Phytoplankton taxa whose longest axial dimension exceeds 30 μm (categorized as grazing resistant) were initially absent, but three such taxa spontaneously colonized the mesocosms during the experiment: Pseudanabaena sp., Closteriopsis acicularis, and Synedra acus, Pseudanabaena and Closteriopsis occurred in all mesocosms from day 12 and 26, respectively, onward. Closteriopsis became abundant in the ambient mesocosms on days 43 to 68 and Pseudanabaena became abundant in the ambient mesocosms with Daphnia on days 54 to 68; both taxa remained rare in the warmed mesocosms throughout the experiment (Figure 2). Synedra was first recorded in a warmed mesocosm on day 68 but became abundant only in the ambient mesocosms with Daphnia from day 82 onward (Figure 2c,d). The proportional contribution of taxa exceeding 30 μm to total phytoplankton biomass was negligible in the warmed mesocosms but exceeded 50% in most of the ambient mesocosms from day 43 onward (Figure 4a; Appendix S1: Table S11, main temperature effect and time × temperature interaction p < 0.001).

FIGURE 4

(a) Time course of the proportional contribution of taxa with longest axial dimension >30 μm to total phytoplankton biomass. Shown are mean values ±1 SE of the eight experimental treatments as indicated in the legend. N = 3 per treatment, except for ambient, Daphnia, where N = 2. (b) Similarity in phytoplankton community composition within and between temperature treatments. Shown is the mean percentage overlap (±95% confidence interval) in (biovolume‐based) community composition in all possible pairwise comparisons among the 11 ambient mesocosms (blue, N = 55), among the 12 warmed mesocosms (red, N = 66), and between each of the 11 ambient with each of the 12 warmed mesocosms (purple, N = 132). Appendix S1: Table S12 gives additional statistics on overlap in community composition within and between treatments

Overall, temperature was the main determinant of community composition. The ambient and warmed mesocosms diverged quickly (average community overlap 9%–26% from day 26 onward), whereas all the ambient mesocosms (average overlap >40%) and in particular all the warmed mesocosms (average overlap >81%) remained much more similar to each other throughout the experiment (Figure 4b, Appendix S1: Table S12). In contrast, grazers had almost no effect on the taxonomic composition of the phytoplankton community. For each temperature treatment and throughout the experiment, community overlap within grazer treatments did not differ from community overlap across grazer treatments, with one exception (Appendix S1: Figure S2). Ambient treatments with and without Daphnia started to diverge on day 54 and differed on days 82–95, when first Pseudanabaena and then Synedra became abundant in the presence of Daphnia (Figure 2a–d, Appendix S1: Figure S2A). In opposition to Hypothesis B2, ciliates did not promote grazing‐resistant taxa >30 μm (Figure 4a). Given that Daphnia performed better in the ambient treatments, where taxa >30 μm contributed more than 50% of the phytoplankton biomass from day 43 onward, it seems likely that Daphnia did ingest and grow on these taxa.

MODEL ANALYSES

When model systems with ambient and warmed phytoplankton traits are forced with the experimental temperature regimes (dotted lines in Figure 1a), their plankton communities enter contrasting seasonal trajectories very similar to those observed in the experiment (Figures 1 and 5). In ambient simulations, system dynamics depend on whether Daphnia are present or not. In the presence of Daphnia, a typical spring succession occurs with a moderate phytoplankton bloom followed by a ciliate peak, then a Daphnia peak, and finally a clear‐water phase (Figure 5b,e,f, dark blue lines); in the absence of Daphnia, phytoplankton and ciliates decline only weakly after their peaks, phytoplankton biomass settles to a moderately high C:P ratio, and a clear‐water phase does not occur (Figure 5b,d–f, light blue lines). In contrast, in warm simulations, a massive phytoplankton bloom develops rapidly and persists at extremely high C:P ratios throughout the entire spring period, regardless of the presence/absence of Daphnia (Figure 5b,d, red and orange lines). When Daphnia are present, they show a rapid but short, transient increase to a peak level considerably lower than in ambient simulations and subsequently decline to near extinction (Figure 5e, red lines). The presence of ciliates also reproduces most of the experimentally observed patterns: Ciliates show temperature dependent bloom dynamics (reaching an earlier peak under warming, Figure 5f), reduce the peak biomass of phytoplankton and their C:P ratio (Figure 5b,d), and weakly promote Daphnia population growth (Figure 5e).

FIGURE 5

Numerical simulations mimicking the eight experimental treatments as indicated in the legend. Time course of (a) soluble reactive phosphorus (SRP), (b) phytoplankton biomass, (c) total phosphorus (TP), (d) seston carbon:phosphorus ratio (C:P), (e) Daphnia biomass, and (f) ciliate biomass

ZNGI analysis of the simplified model without ciliates reveals the mechanisms causing these contrasting dynamics. A system with ambient phytoplankton traits has a single attractor at all temperatures, which is a globally stable equilibrium where Daphnia controls phytoplankton at low density and in a state of high nutritional value (Figure 6a–c). Starting from low late‐winter densities, this equilibrium is reached via a damped oscillation similar to a phytoplankton spring bloom followed by a clear‐water phase (Figure 6a–c). In contrast, a system with warmed phytoplankton traits has such a single attractor only at temperatures ≤10.6°C (Figure 6d). At higher temperatures, a second stable coexistence attractor emerges where phytoplankton biomass is 50–100 times higher but of very low nutritional value (Figure 6e). At this second attractor, increasing temperature causes phytoplankton biomass to increase and its nutritional value to decrease even further. Consequently, Daphnia biomass declines with warming until, at temperatures ≥15.1°C, the second attractor turns into a phytoplankton‐only system and Daphnia goes extinct (Figure 6f). Starting from low late‐winter densities, only that second attractor can be approached (Figure 6e,f). Thus, sufficiently rapid early season warming can force Daphnia‐phytoplankton dynamics onto a trajectory toward this alternative, high phytoplankton–Daphnia extinction attractor.

FIGURE 6

Attractors of the simplified phytoplankton–Daphnia model system in absence of microzooplankton at 10°C (panels a and d), 14°C (panels b and e), and 20°C (panels c and f). Shown are zero net growth isoclines (ZNGI) of Daphnia (blue lines) and of phytoplankton populations (green lines) with production and sinking traits representative of ambient (panels a–c) and warmed (panels d–f) communities. Solid black lines are phytoplankton–Daphnia trajectories starting from low initial densities (green triangles, phytoplankton 50 mg C m−3, Daphnia 10 mg C m−3). In panels (a–d), the system goes through one or several transient cycles before it converges to a stable, low phytoplankton–Daphnia equilibrium (red circles). In panel (e), the system reaches a stable, high phytoplankton–Daphnia equilibrium and in panel (f) a stable, high phytoplankton–Daphnia extinction equilibrium (black circles). The warmed 14°C and 20°C systems have a stable, low phytoplankton–Daphnia equilibrium as an alternative attractor (panels e and f, red circles) that would be reached from initial conditions above the dotted line (which separates the basins of attraction of the alternative stable equilibria). The open circles in panels (e) and (f) denote an unstable equilibrium (a saddle point). In each panel, the nutrient quota of phytoplankton is at its maximum (molar C:P = 65) in the dark gray area and decreases continuously toward the broken line, which marks system states where the phytoplankton nutrient quota is at its minimum (C:P = 2600). Boundaries between gray scales indicate, from left to right, C:P ratios of 65, 250, and 850, respectively. System states beyond the broken line are impossible because the total amount of nutrient in the system cannot support a higher carbon biomass. Daphnia population growth is positive below the grazer isocline and negative outside of it, where phytoplankton biomass is either too low (to the left) or too nutrient deplete (above and to the right) for positive Daphnia growth. Phytoplankton growth is positive below its isocline and negative above it. The phytoplankton isocline has an inflection point because of Daphnia's type 3 functional response. Note the log scale of the x‐axes

Bifurcation analyses reveal how system dynamics depend on environmental conditions and algal traits. A Daphnia extinction equilibrium (either as the only possible outcome or as one of two alternative states) is facilitated by lower nutrient enrichment, a shallower water column, and lower background light attenuation (Figure 7a,b). Under given environmental conditions, Daphnia extinction is promoted by a high maximum growth rate and a low minimum nutrient quota (high C:P) of the algal food (Figure 7c–f). Algal sinking velocity has a weaker impact (Figure 7g,h) because both low and high sinking losses can contribute to low algal food quality. Low sinking losses sustain high algal carbon biomass, while high sinking losses reduce nutrient availability in the water. Both mechanisms reduce algal nutritional value. Finally, algal carrying capacity—and thus algal C:P—increases with temperature (see intersections of algal isoclines with the x‐axes in Figure 6) because the algal growth rate increases faster with warming than algal losses (Appendix S2: Table S2). In contrast, the net growth of Daphnia at such saturating algal food densities is hardly affected by temperature because food processing (handling time) and loss rate have similar temperature dependence (Appendix S2: Table S2). Consequently, warming increases food abundance but decreases food quality, promoting the existence of an alternative Daphnia extinction equilibrium at higher temperatures (Figure 7c–h).

FIGURE 7

Attractors of the phytoplankton–ciliate–Daphnia model system as functions of environmental conditions (water column depth, total nutrient concentration, temperature) and algal traits (maximum growth rate, minimum nutrient quota [expressed as maximum C:P ratio], sinking velocity). Colors indicate parameter combinations for which a low phytoplankton–low ciliate–Daphnia equilibrium is (i) the only attractor (yellow), (ii) coexists with a high phytoplankton–high ciliate–Daphnia equilibrium (blue), or (iii) coexists with a high phytoplankton–high ciliate–Daphnia extinction equilibrium (green); (iv) in black areas, Daphnia cannot persist at all. With the exception of two parameters and the parameters that are varied along the plot axes, all parameter values are set to the environmental conditions of the experiment (water column depth = 1.5 m; initial total nutrient concentration = 28 mg P m−3; background light attenuation coefficient = 1.0 m−1) and to the default values listed in Appendix S2: Table S2, assuming trait values of the warmed phytoplankton communities (maximum algal growth rate at 20°C = 1.9 day−1; maximum algal molar C:P ratio = 2600; sinking velocity = 0.05 m day−1). Panel (b) differs in that it mimics a clear water (background light attenuation coefficient = 0.25 m−1). Panels (c, e, and g) differ in that they mimic a nutrient‐poor system (initial total nutrient concentration = 14 mg P m−3). The black diamond in panel (a) and the broken vertical lines in panels (d, f, and g) indicate conditions representative of the warmed experimental treatments

DISCUSSION

Our experimental data suggest that processes involved in three of the four hypotheses described in the introduction acted simultaneously on the plankton community. In accordance with Hypothesis A1, all planktonic organisms (phytoplankton, ciliates, Daphnia) responded positively to warming at the initially low late‐winter temperatures. Yet, in contrast to Hypothesis A1, warming did not accelerate the phases of a typical spring succession of the plankton. Instead, and in agreement with Hypothesis A2, warming promoted the build‐up of carbon‐rich phytoplankton biomass of low nutritional quality that subsequently reversed the direct positive effect of warming on Daphnia population growth. The consequences for the seasonal trajectory of the plankton community were, however, much more dramatic than anticipated by Hypothesis A2. The population growth of Daphnia and a resulting clear‐water phase were not just delayed; rather, warming caused a regime shift: Daphnia went near‐extinct and a clear‐water phase did not occur at all.

As expected, ciliates grew faster than Daphnia at low late‐winter temperatures. In accordance with Hypothesis B1, ciliates slowed down the population growth of phytoplankton and, thus, the deterioration of its nutritional value. At ambient temperature, this had all of the consequences predicted by Hypothesis B1: phytoplankton biomass reached an earlier and lower peak, Daphnia abundance reached an earlier and higher peak, and the clear‐water phase was slightly advanced. In the warmed mesocosms, ciliates also slowed down phytoplankton growth and temporarily increased Daphnia abundance but could only delay the regime shift. Eventually, phytoplankton grew out of grazer control also in the presence of ciliates, and the Daphnia population crashed. In contrast to Hypothesis B2, ciliates did not shift the phytoplankton community toward more grazing‐resistant forms (defined as taxa >30 μm) and had an overall positive effect on Daphnia. It is possible that direct consumption of ciliates contributed to this positive effect on Daphnia (Stoecker & Capuzzo, 1990; Wickham, 1998). Ciliates have low C:P ratios and can thus “trophically upgrade” low‐quality algal food for omnivorous mesozooplankton (Golz et al., 2015; Martin‐Creuzburg et al., 2005).

The model analyses strongly suggest two key mechanisms that, together, triggered the observed regime shift under the warming scenario. First, because ambient and warmed phytoplankton were assumed to be equally susceptible to ingestion by grazers but reached much higher biomass in warmed than ambient model runs (Figure 5b), the Daphnia decline in warmed simulations was entirely caused by the developing stoichiometric mismatch between algal nutrient content and grazer nutritional needs (Figure 5d). The resulting dynamics are fundamentally different from those expected under a mismatch with food quantity, where the timing of peak resource abundance and peak energetic needs of consumers is mismatched (Durant et al., 2007). In contrast, food of edible size was abundant during the period of Daphnia decline in warmed treatments, in both the simulations and the experiment (Figures 1b and 5b, red and orange lines). The concomitant food quality mismatch explains why, in the experiment, peak abundance of Daphnia was strongly negatively correlated with peak chl a concentration (Figure 3h, Pearson correlation coefficient = −0.91, p < 0.001). In laboratory experiments, Sommer (1992) observed a positive population growth of D. hyalina when the C:P ratio of its food (Scenedesmus acutus) was ≤885, whereas the stocked Daphnia population declined to extinction when the C:P ratio was ≥980. This is very similar to our observation that Daphnia population growth turned negative at seston C:P ratios of 1000–1100. To our knowledge, our experiment is the first empirical demonstration of the “paradox of energy enrichment” (i.e., grazer starvation in an abundance of energy‐rich but nutritionally imbalanced food) (Andersen et al., 2004; Diehl, 2007; Loladze et al., 2000) in a multispecies phytoplankton community. The phenomenon was previously observed in simpler laboratory settings where monocultures of S. acutus at steady state were inoculated with Daphnia (Sommer, 1992; Urabe et al., 2002).

Second, the regime shift is promoted by a high phytoplankton growth rate, low minimum algal nutrient quota, low algal sinking losses, high temperature, and low nutrient and high light supply. The former enables phytoplankton to rapidly grow out of grazer control at the onset of spring warming, and the remaining factors together ensure that a high phytoplankton biomass of low nutritional value can be sustained over long periods. Under these circumstances, grazer populations remain small and are therefore unable to mitigate the unfavorable quality of their food through grazing and nutrient cycling (Sommer, 1992; Urabe et al., 2002). In contrast, the slower growth rates implemented in ambient model runs prevent phytoplankton from growing out of control, and grazers keep phytoplankton at a low but nutrient‐replete biomass of high food quality, similar to experimental observations.

Interestingly, the warming‐induced, very rapid taxonomic shift in the phytoplankton community occurred also in the absence of grazers (Figure 2a,e), suggesting that temperature alone was responsible for the differences in community development during the early phase of the experiment up to day 26. This indicates that Monoraphidium had a higher specific growth rate than its competitors at temperatures >8°C (the starting temperature in the warmed mesocosms), but not at lower temperatures. We cannot address this suggestion with data because temperature‐growth curves including low early‐winter temperatures are not available for the taxa from our experiment. Yet among the three taxa that together accounted for the lion's share of total algal biovolume during the first month of the experiment (>75% and 91%, respectively, in the ambient and warmed mesocosms), M. minutum has a considerably higher maximum growth rate than S. obliquus and Cyclotella meneghiniana (Thomas et al., 2016). Typically, a high maximum growth rate at optimum temperature comes at the expense of a low growth rate at low temperatures (Eppley, 1972; Norberg, 2004; Thomas et al., 2016), giving some credibility to the hypothesis that the differential success of Monoraphidium was caused by a temperature trade‐off in its growth performance.

If the preceding scenario is correct, then why did Monoraphidium not come to dominance in the ambient treatments at a later stage of the experiment? When temperature reached 8°C in the ambient mesocosms on day 39, total nutrients were already strongly depleted and continued to decrease, probably because the then dominant, relatively large and nonmotile taxa Scenedesmus, Closteriopsis, and Cyclotella experienced high sedimentation losses. This shortage of nutrients likely prevented the bloom from reaching the extreme biomass levels observed in the warmed treatments. Still, toward the end of the experiment, Monoraphidium became dominant and prevented a clear‐water phase also in ambient treatments (Figure 1b), but only in the absence of Daphnia (Figure 2a–d). In line with model predictions (Figure 6e,f) and laboratory experiments (Sommer, 1992; Urabe et al., 2002), the latter indicates that once Daphnia has established a sufficient population, it can control its algal food and maintain it in a nutritionally adequate state.

Because it is impossible to separate phytoplankton from microzooplankton in natural lake water, we performed our experiment with assembled communities of precultured taxa. We cannot rule out the possibility that this only moderately diverse phytoplankton community was predisposed to the competitive dominance of a single taxon with extremely low nutritional value. Moreover, regardless of the bulk nutrient‐to‐carbon stoichiometry of phytoplankton, herbivore growth can be positively affected by taxonomic diversity of its food per se (Behl & Stibor, 2015; Gamfeldt et al., 2005; Marzetz et al., 2017; Striebel et al., 2012). For example, Urabe and Waki (2009) found that higher algal species richness mitigated the adverse effects of low algal phosphorus content on the growth of Daphnia, and Behl et al. (2012) observed both higher algal C:P ratios and higher growth rates of Daphnia in algal polycultures than in monocultures.

Our assembled community was thus likely more vulnerable to a warming‐induced regime shift than would be a more diverse, natural community. Yet, the observed increase in carbon to nutrient stoichiometry is qualitatively consistent with similar warming trends described from other mesocosm experiments, laboratory studies of the intraspecific temperature dependence of algal carbon to nutrient stoichiometry, and observations across vast geographical temperature gradients (De Senerpont Domis et al., 2014; Martiny et al., 2013; Rhee & Gotham, 1981; Taucher et al., 2012; Toseland et al., 2013; Yvon‐Durocher et al., 2015). The proposed reduction in sinking losses is also consistent with the general observations that warming leads to a reduction in phytoplankton cell size (Daufresne et al., 2009; Yvon‐Durocher et al., 2011) and that sinking velocity declines with cell size (Miklasz & Denny, 2010; Ptacnik et al., 2003).

Because higher temperatures promote thermal stratification and thereby reduce nutrient inputs to the upper mixed layer of lakes and oceans, Van de Waal et al. (2010) proposed that the resulting decrease in nutrient availability and a concomitant increase in carbon availability (driven by elevated CO2 levels) may cause phytoplankton carbon:nutrient ratios to increase in the future. Increased phytoplankton carbon:phosphorus ratios are also predicted by models that take into account effects of warming on cellular physiology and resource allocation (Shuter, 1979; Toseland et al., 2013). Our study suggests that reduced mineral nutrient content in phytoplankton biomass can trigger positive feedbacks similar to those observed in systems where warming benefits grazing‐resistant phytoplankton taxa, such as cyanobacterial colonies, accentuating a food quality mismatch with grazers and thereby promoting undesired or even harmful algal blooms (O'Neil et al., 2012; Paerl & Huisman, 2008; Wilson et al., 2006). Future climate change may thus exacerbate food quality mismatching at the producer–grazer interface (Anderson et al., 2017; Persson et al., 2011), altering the functioning of many ecosystems and impacting ecosystem services.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS

Sebastian Diehl, Stella A. Berger, and Herwig Stibor designed the experiment, and Stella A. Berger, Sebastian Diehl, and Herwig Stibor led and supervised its execution. Sebastian Diehl and Stella A. Berger analyzed the experimental data, and Wojciech Uszko and Sebastian Diehl performed model analyses. Sebastian Diehl wrote the manuscript, and all authors contributed to revisions.

Appendix S1

Click here for additional data file.

Appendix S2

Click here for additional data file.