Eco-evolutionary interaction in competing phytoplankton: genotype sorting likely explains dominance shift and species responses to CO2

How ecological and evolutionary processes interact and together determine species and community responses to climate change is poorly understood. We studied long-term dynamics (over approximately 200 asexual generations) in two phytoplankton species, a coccolithophore (Emiliania huxleyi) and a diatom (Chaetoceros affinis), to increased CO2 growing alone or competing with one another in co-occurrence. To allow for rapid evolutionary responses, the experiment started with a standing genetic variation of nine genotypes in each of the species. Under co-occurrence of both species, we observed a dominance shift from C. affinis to E. huxleyi after about 120 generations in both CO2 treatments, but more pronounced under high CO2. Associated with this shift, we only found weak adaptation to high CO2 in the diatom and none in the coccolithophore in terms of species’ growth rates. In addition, no adaptation to interspecific competition could be observed by comparing the single to the two-species treatments in reciprocal assays, regardless of the CO2 treatment. Nevertheless, highly reproducible genotype sorting left only one genotype remaining for each of the species among all treatments. This strong evolutionary selection coincided with the dominance shift from C. affinis to E. huxleyi. Since all other conditions were kept constant over time, the most parsimonious explanation for the dominance shift is that the strong evolutionary selection potentially altered competitive ability of the two species. Thus, here observed changes in the simplest possible two-species phytoplankton “community” demonstrated that eco-evolutionary interactions can be critical for predicting community responses to climate change in rapidly dividing organisms such as phytoplankton.


1
Introduction 38 Recent studies have repeatedly shown that ecological and evolutionary processes happen on similar 39 time scales (Carroll et al., 2007;Reznick, 2013). Understanding how both processes interact has 40 sparked the emergence of the new field of eco-evolutionary dynamics (Fussmann et al., 2007) aiming 41 to understand how rapid evolutionary adaptation influences ecological processes, and vice versa 42 (Hendry, 2016). Strong eco-evolutionary coupling is particularly expected in phytoplankton species, 43 photoautotrophic aquatic microbes that can form large populations. With one cell division per day, 44 many species lend themselves for experimental studies allowing for appropriate replication and 45 hundreds of generations of evolutionary change (Reusch and Boyd, 2013). Moreover, their ecological 46 importance cannot be overemphasized, as marine phytoplankton species are responsible for half of all 47 global photosynthesis (Falkowski et al., 2008). 48 Here we studied eco-evolutionary dynamics in response to increased seawater CO2 concentration, one 49 prominent aspect of ocean global change (Doney et al., 2009). As primary producers, all phytoplankton 50 species depend on the availability of inorganic carbon. On time scales too short for evolutionary 51 adaptation, non-calcifying species mostly respond by enhanced growth (Li et al., 2017;Schaum et al., 52 2012) leading to higher abundances (Sommer et al., 2015) ("ecological winners", e.g. diatoms), while 53 most calcifying species react negatively ("ecological losers", e.g. coccolithophores) to CO2 enrichment 54 (Doney et al., 2009;Riebesell, 2004). Nested within such broad functional categories of sensitivity is on which selection can operate (Lohbeck et al., 2012). 57 Within phytoplankton, competition for abiotic resources is ubiquitous (Tilman, 1977 (Lawrence et al., 2012), competition was shown to constitute one major biotic 60 driver of adaptation. At the same time experimental evolution studies using phytoplankton and 61 addressing the presence and absence of competition as a selection factor are largely absent (Scheinin 62 et al., 2015). 63 Here, we set out to study how the evolutionary response to increased CO2 of two bloom forming and 64 geographically co-occurring phytoplankton species, Emiliania huxleyi (a coccolithophore) and 65 Chaetoceros affinis (a diatom), was altered by competition and assessed putative effects on ecological 66 dynamics. Studying such eco-evolutionary processes requires that at least two species can be kept over 67 long-term (i.e. more than 10s of generations such that evolutionary change can happen) in co-occurring 68 experimental settings, delaying or preventing Gause's principle of competitive exclusion (Hardin 69 Garret, 1960;Tilman, 1977). Hence, we established long-term coexistence in semi-continuous batch 70 cycles by taking advantage of the species´ different nutrient uptake related strategies (Doney et al.,71 2009; Riebesell, 2004;Sommer, 1984). Diatoms have high nutrient uptake rates and are consequently 72 characterized by high maximum growth rates (Litchman et al., 2007). At the same time, diatoms have 73 high half saturation constants and low affinity for nutrient uptake (Litchman et al., 2007). As such they 74 represent the 'velocity-adapted' species (Sommer, 1984) and thrive in nutrient replete and fluctuating 75 conditions, being able to rapidly monopolize nutrients. Coccolithophores in contrast, have lower 76 maximum nutrient uptake rates. Hence, after a nutrient pulse they initially lose out against diatoms. 77 Their lower half-saturation constants and higher affinity, however, make them the better competitors 78 under nutrient poor conditions thriving at later successional stages (Sommer, 1984). In this system, we 79 ran the experiment for ca. 200 asexual generations under fully factorial selection conditions of 80 adaptation to ambient and increased CO2, in combination with and without competition, i.e. the co-81 occurrence of the respective other species (Fig. 1 a). To allow for rapid evolutionary adaptation via genotype sorting (Lohbeck et al., 2012), populations of both species were assembled using nine 83 genotypes in equal frequencies (Fig. 1 a) that could be traced via microsatellite genotyping. Reciprocal  84  adaptation assays were conducted to test for adaptation of both species to enhanced CO2 and the  85  presence of the respective other species in the first and second half of the experiment (after ca. 50 and  86 200 generations) (Fig. 1 b). 87 We hypothesize that on the short term, species contribution to the "community" would be diverging 88 between CO2 treatments due to varying CO2 responses of the target species. Specifically, the relative 89 contribution of E. huxleyi should be lower in the high CO2 compared to the ambient treatment due to 90 potential negative or neutral effects of CO2. In contrast, the relative contribution of C. affinis should 91 be higher in high CO2 treatment due to a potential fertilizing effect of increased CO2 in combination 92 with the effect of CO2 on E. huxleyi. In the long-term, we hypothesize that evolutionary dynamics (i.e. 93 genotype sorting and eventually adaptation) diverge between both target species. Specifically, the 94 potentially negatively affected coccolithophore should be selected towards more stress tolerance, while 95 the possible fertilizing effect on the diatom could lead to selection for the diatom's potential to 96 efficiently utilize increased levels of CO2 and in turn achieve higher growth rates. Experimental set up and culturing conditions 120 The selection treatments were set-up in a fully orthogonal way with the single factors CO2 (two levels, 121 ambient and high) and co-occurrence or absence of a 2 nd species (Mono and Mix) resulting in 4 122 treatment combinations per species (Fig. 1 a) µmol*L -1 nitrate, 0.99±0.02 µmol*L -1 phosphate and 4.34±0.14 µmol*L -1 silicate. In order to simulate 133 a bloom situation with resource competition under depleted nutrients (Paulton, 1991), we ran the study 134 as semi-continuous batch cultures with both species reaching the stationary phase (supplementary 135 Figure S1 growth graphs of the different species). This resulted in the almost full depletion of all 136 macronutrients (nitrate, phosphate and silicate) at the end of the batch cycles (supplementary Fig S2). 137 Each experimental unit was inoculated with an initial total biomass of 8280 µm 3 *mL −1 . When the two 138 species were growing together, they were added at a ratio of 1:1 of the biomass to start the experiment. 139 The nutrient concentrations and starting conditions were determined in prior pilot studies (Listmann 140 and Hattich unpublished data) to ensure coexistence for at least 50-100 generations. The experiment 141 was carried out under constant rotation (0.75 min −1 ) at ~22 °C and a 17:7 day:night cycle reaching a 142 maximum light intensity of 350 μmol*m -2 *s -1 three hours after dusk and dawn. At the end of each 143 batch cycle cell numbers and cell volume were determined using microscopy (Zeiss Axiovert 144 Observer) to calculate the contribution of each species to the total biomass and the amount needed to 145 transfer to the next batch cycle; again 8280 µm 3 *mL -1 . 146

Frequency assessments of genotypes 147
In order to follow genotype frequencies, we re-isolated cells of both species after 8, 32, 64, 160 and 148 288 days. Genotypes of both species could be unambiguously identified via microsatellite genotyping 149 (see (Hattich et al., 2017), for detailed information on E. huxleyi and for C. affinis genotype 150 identification see supplementary material). For the quantification of E. huxleyi and C. affinis genotypes 151 a maximum of 20 cells per culture were re-isolated by dilution in 48 well plates. This provided a 152 theoretical detection limit of 5% difference between the contributions of the genotypes to both species´ 153 populations. A lower detection limit could in theory have been achieved via isolation and identification 154 of more cells but was not the goal of this study. Details on the reisolation and quantification via 155 microsatellite analysis are given in supplementary material. For the genotype composition of C. affinis 156 we only had enough data for analysis two time points in the mixed cultures. 157

Reciprocal adaptation assays 158
We carried out reciprocal adaptation assays to test for adaptation in both species to enhanced CO2, the 159 competition of the respective other species, and both factors in combination ( Fig. 1 b). The assays on 160 day 64 and 288 assessed how adaptation played out over time. Those adaptation assays compared 161 evolved populations in control and new environments rather than evolved and ancestral populations 162 because methods such as cryopreservation are not readily available for our target study species 163 (Collins, 2011b). Specifically, every evolved culture was tested in all four treatment combinations in a 164 full factorial way. For example, E. huxleyi that were long-term treated with high CO2 conditions 165 (treatment combination Mono High), were in the assay exposed to both ambient and high CO2 166 concentrations and then also to the co-occurring diatom (already long-term exposed to ambient and 167 high CO2 and competition to avoid confounding responses) in ambient and high CO2 concentrations. 168 Thus, each treatment combination from the selection leads to 4 assay treatment combinations. This 169 resulted in a total of 16 assay treatment combinations per species (Fig. 1 b; for the details of the 170 reciprocal assay set up see supplementary material). Using a full factorial adaptation assay allows us 171 to test on the one hand for the adaptation of the single factor CO2 (i.e. in the statistical analysis, this 172 could be identified as significant interactions of the respective selection and assay treatment factors 173 (e.g. selection CO2 x assay CO2)). On the other hand, we can additionally test for the adaption to 174 increased CO2 in combination with competition but then also how competition itself affected adaptive 175 responses. This allows disentangling the single and combined treatment factors at the same time. Here 176 we focused on growth rates because this response is directly related to Darwinian fitness in an asexual 177 population (Elena and Lenski, 2003) and is independent of nutrient availability due to the presence of 178 another species (Tilman, 1977). 179

Statistical analyses 180
In all treatments with co-occurring species, the relative species contribution of E. huxleyi to total 181 biomass (Mix) over time was analyzed using a generalized least squares (GLS) model (m0<-gls 182 (relative Biomass ~ Selection CO2 * Time)). Differences in variance structure were adjusted for the 183 factor "Selection CO2", and accounted for autocorrelation over time. For the statistical analysis it was 184 sufficient to only look at the relative contribution of one of the two species, as their respective 185 contribution to the community was complementary. For the analysis of the absolute biomass of each 186 species in all cultures we used the same GLS model. The analysis was done for each species separately 187 and started with the following full model: m0<-gls (Biomass ~ Selection CO2 * Selection Culture * 188 Time). After reduction and analysis of the model structure we had to account for autocorrelation over 189 time and differences in variance structure (varIDENT) for the "Selection Culture" in E. huxleyi and 190 "Selection CO2" and "Selection Culture" in C. affinis. There was a strong change around 160 days of 191 experiment and we found significant interactions with "Time" for all factors in the GLS model. In 192 order to investigate the effects of the different treatments on both species before and after this time 193 point, we divided the time series data into two parts -BC1-BC20 and BC21-BC36 on which we 194 repeated the described GLS model (BC1-20 and BC21-36, respectively<-gls (Biomass ~ Selection CO2 195 * Selection Culture)). 196 A permutational multivariate analysis of variance (permanova, with 999 permutations (using the 197 package "vegan")) was used to test for the differences in genotype compositional change between the 198 treatments. 199 To analyse the reciprocal assay data we used repeated measure analysis of variance (rmANOVA). 200 Since during the experiment, 4 and 3 cultures were lost, we had to omit the data of the "lost" replicates 201 from all statistical analyses of the assays. Before starting the analyses we tested for the homogeneity 202 of variance using a Fligner Killeen Test (Fligner and Killeen, 1976) and for the normality of the 203 residuals using a Shapiro-Wilk Normality Test (Shapiro and Wilk, 1965). The assumption of sphericity 204 was not violated because we only had one repeated measure (Field et al.). All assumptions were met 205 such that we could continue with the analysis without data transformation. The growth rate was first 206 assessed for the overall effects of "Time" (between the two assays), "CO2" and "Culture" in both 207 selection and assay environment and thus we started with a complete data set analysis 208 (aov((muExp)~(Selection CO2 * Selection Culture * Assay CO2 * Assay Culture * Time) + 209 Error(Replicate/Time))) followed by a separate analysis of each assay. Here, interactions of selection 210 treatment * assay treatment would indicate evolutionary adaptation over the course of the experiment. 211 All modelling and statistical analyses were done using the software R (Coreteam, 2016 The outcome of interspecific competition in ambient and high CO2 was mirrored in relative biomass 217 shifts: In the first half of the experiment the diatom dominated the two-species community with a 218 relative biomass (mean over first 20 batch cycles) of ca. 95% in high compared to ambient CO2 219 condition with 80% (Fig. 2

Temporal dynamics of absolute species biomass 227
Similar to the dynamics of relative species contributions, the time course of absolute biomasses of both 228 species was characterized by two distinct phases that changed around 160 days of the experiment (Fig.  229 3 a-d, increase across all treatments (Fig 3a and  cultures whereas in mix-culture increased CO2 even had a slightly positive effect (Fig. 3a, b and  In contrast to E. huxleyi the biomass of C. affinis varied markedly between batch cycles with only a 242 slight decrease of biomass over time (Fig 3c and d The genotype composition changed uniformely over time in both species (Fig. 4,  and B57, respectively). In E. huxleyi there was no effect of CO2 on the genotype sorting (Fig. 4, Table  258 S6, "Selection CO2" R 2 =0, p=0.401, permutations=999) and only a small difference between the single 259 and co-occurring cultures (Fig. 4, Table S6 "Selection culture" R 2 =0, p=0.085, permutations=999). The 260 small difference likely came about due to slightly slower sorting to the dominant genotype C41 in the 261 mix-culture (Fig. 4 a vs b). For C. affinis we could not analyze the effect of either treatment on the 262 genotype composition. 263

Reciprocal adaptation experiments 264
In E. huxleyi we found no significant interaction of selection with the assay treatment in neither of the 265 two assays which would indicate evolutionary adaptation (Table S7, E. huxleyi). Specifically, long-266 term selection to high CO2 did not result in increased growth in the high CO2 assay condition compared 267 to ambient selected populations (Fig 5a and b). However, there were non-adaptive effects of CO2 and 268 Culture selection and assay treatments that varied between the two assay experiments: After 64 days, 269 selection under high CO2 had a negative effect on the growth rates of E. huxleyi when measured in the 270 presence of the diatom ("mix" assay treatment), while "mix" assay conditions led to a decline in growth 271 rates in general ( Fig. 5a and  responses to CO2 with and without competition that even resulted in an seemingly adaptive effect to 282 increased CO2 in the second assay after 288 days, evident as a significant interaction of selection and 283 assay treatment (Fig 5 d "u-shape" pattern in pooled data, Our results show that long-term selection in communities to an abiotic stressor can drastically diverge 291 from the predicted outcomes of short term experiments with the same abiotic driver (Hattich et al.,292 2017; Riebesell, 2004). This divergence potentially arises from the effect of evolutionary on ecological 293 processes, as previously shown in other systems (Hairston et al., 2005;Schoener, 2011). During the 294 first half of the experiment our initial expectations (for biomass) with respect to species being an 295 ecological "loser" and "winner" in response to enhanced CO2 were met for E. huxleyi but not C. affinis. 296 This led to the expected dominance of C. affinis, which was more pronounced under increased CO2 297 conditions. However, over time, the observed effects of increased CO2 on biomass changed for both 298 species so as to revert our predicted ecological "loser" and "winner" outcome. Consequently, the 299 relative contributions in the communities flipped and E. huxleyi became dominant under both CO2 300 conditions. We note that the "dominance-flip" coincided with strong genotype selection in both species, 301 which was however not driven by CO2 or competition. Likewise, neither species showed adaptation to 302 CO2 or competition in the reciprocal assays. 303 The absent selective and adaptive response to CO2 in both species, is for coccolithophores in contrast 304 to previous studies (Lohbeck et al., 2012) and for diatoms both in contrast (Scheinin et al., 2015) and 305 agreement (Li et al., 2017 interspecific competition has the potential to affect evolution in phytoplankton, as already intraspecific 324 competition was shown to alter adaptive responses to CO2 (Collins, 2011a). Here, the dominant C. 325 affinis strongly reduced the biomass and growth rate of E. huxleyi, demonstrating that the effect of 326 competition was strong, especially in the first half of the experiment. In addition, we found a further 327 reduction of growth rates to increased CO2 with the co-occurrence of the diatom indicating that the 328 ecological context did play a role albeit not affecting the evolutionary change. That no effect of 329 competition on evolution, or adaptation to competition could be observed potentially depends on how 330 coexistence was allowed in our experiment. Owing to the species' different nutrient uptake strategies 331 (Sommer, 1984) niche partitioning likely appeared temporary over the course of a batch cycle. Whereas 332 the diatom had an advantage during the first batch cycle days with replete nutrients, the coccolithophore 333 was favored and could grow longer towards the end of each batch cycle when nutrients became limiting 334 (Sommer, 1984) (Fig. S1-1 and S1-2). As such selection on standing genotype composition likely took 335 place within these two niche spaces of the both coexisting species, unaffected by an interspecific 336 competitor. 337 Despite the absence of an adaptive response to increased CO2 or competition, we still observed strong 338 evolutionary change as very reproducible, directional sorting of genotypes across all treatments. 339 Interestingly, these evolutionary dynamics selected for one single genotype of each species that became 340 dominant already after approximately 64-160 days of the experiment. While we know that evolution 341 was independent of the intended CO2 and competition treatments, non-intended laboratory selection 342 was the most likely driver in our study. Laboratory selection is inevitable and often strong and has been 343 previously demonstrated over very long time scales by comparing ancestral bacterial populations to 344 laboratory evolved populations (Lenski and Travisano, 1994). Here, selection became apparent in an 345 overall increased E. huxleyi and decreased C. affinis biomass and changing growth rates over time (Fig.  346  S1-3). Interestingly, an increase of E. huxleyi growth rate was also observed in another long-term 347 laboratory selection study (Schlüter et al., 2016). In these experiments however, laboratory adaptation 348 was only a "background" signal, whereas in our experiments the chosen treatments were too weak to 349 impose a selection force to overcome laboratory selection. Schlüter  selection also occurred in an earlier attempt of the experiment, suggests that the eco-evolutionary 366 dynamics observed are reproducible (Fig. S1-4 and S1-5). We propose that selection by nutrient 367 limitation on standing variability in genotype's competitive abilities was the most likely driver 368 resulting in the eco-evolutionary interaction. However, necessary experimental characterization of the 369 genotypes´ nutrient uptake associated traits is (currently) missing. Other studies reported nutrient 370 uptake related traits in marine protists to vary between genotypes to the same extent as growth rates 371 (Brandenburg et al., 2018). Since we have shown that the variability in growth rates among the 372 genotypes used in this study is even larger than that of the response to increased CO2 (Hattich et al.,373 2017), we postulate that their nutrient uptake related traits vary as well. Another indication that 374 nutrients were driving the change in our system is the observed shift away from C. affinis towards E. 375 huxleyi, the favored species by nutrient limitation (Tyrrell and Merico, 2004). Consequently, future 376 studies should focus more on the consequences on how nutrients select to predict phytoplankton change 377 (Thomas et al., 2017). This will become particularly relevant as nutrient uptake related traits explain 378 competitive ability under different nutrient conditions (Litchman et al., 2007;Sommer, 1984) and there 379 will be more nutrient limitation on phytoplankton in the future ocean (Boyd et  not CO2 environment turned an ecological "loser" (with respect to biomass contribution to the 385 community) into a "winner" and vice versa. As such, this flip demonstrates for the first time that eco-386 evolutionary interactions play out in competing phytoplankton communities. Such interactions can 387 drastically alter the effect of environmental drivers and lead to diverging predictions of future changes 388 compared to such resulting of short-term studies. Our results call for an inclusion of more realistic 389 experimental evolution conditions in future studies, not only using realistic nutrient regimes, but more 390 importantly also including multi-species settings and their underlying mechanisms allowing for stable 391 coexistence to simultaneously investigate ecological and evolutionary processes in phytoplankton. 392

Acknowledgements 393
We thank Thomas