Episodic nutrient enrichments stabilise protist coexistence in planktonic oligotrophic conditions

Seasonal compositional changes in plankton communities are usually considered as species replacements. Given the enormous number of individuals integrating the communities and our limited capacity to count and determine most of them, we likely observe only alternative population peaks of some of the coexisting species. The contemporary coexistence theory addresses coexistence in communities of competing species, considering relative fitness inequalities and stabilising niche differences as components of average long‐term growth rates. Here, we experimentally show that the response patterns predicted by the theory occur when varying nutrient pulses fertilise the planktonic community. We used gently self‐filling 100‐L enclosures to minimise the disturbance of the initial community and different pulse P and N additions to manipulate the apparently species‐poor epilimnetic community of an ultraoligotrophic P‐limited lake. We measured and compared the protist species growth response to gradients of P enrichment and N stoichiometric imbalance. The P and N levels used in most treatments were within the oligotrophic seasonal and inter‐annual variations of the lake and were higher in a few extreme treatments that provided mesotrophic conditions of the remote regions affected by N atmospheric contamination. We alternatively replicated all treatments using ammonium or nitrate as the N source. Most protist species, recorded in this lake across seasons in previous studies, were recovered, indicating a persistent assemblage of species that is seasonally hidden from observation. Recovery included some rare species observed only in the slush layers of the seasonal snow and ice cover. The coexistence‐stabilising mechanisms were indicated by treatment response features, such as frequency‐dependent growth, inverse relationship between fitness inequality and niche differentiation proxies, high‐rank taxonomic levels clustering across the limiting‐nutrient gradient but segregation at the species level according to the type of N supply and resting stage development depending on nutrient conditions. The response similarities between autotrophic and heterotrophic organisms indicate a network of interactions that may reinforce coexistence. Synthesis. The results indicate that many planktonic protist species in oligotrophic waters can show stable long‐term non‐equilibrium coexistence by alternately recovering from very low densities when episodic nutrient enrichments, of varying P and N amounts and composition, occur.


| INTRODUC TI ON
Planktonic protist communities undergo species replacement at relatively short time-scales (e.g. seasonally in lakes). However, this view might be conditioned by our limited empirical capacity to penetrate the likely high diversity of rare species coexisting in these protist communities at any time (Smetacek, 2018). The number of protist individuals is extremely high; in oligotrophic lakes and an open ocean, there are approximately 10 3 -10 4 protist cells/ ml (Reynolds, 2006), which is more than a billion per cubic metre, and approximately 10 16 in the photic zone of a lake. Consequently, scientific observations consider only a small fraction of the whole community. Even with current molecular techniques, the sequencing depth is still limited to account for the total diversity of local (a few litres) planktonic microbial communities (Pedros-Alio, 2012).
Given the large number of individuals, an alternative view to the seasonal replacement of species, which implies colonisation from other parts and local extinction, is the stable non-equilibrium coexistence of many species in a fluctuating environment (Li & Chesson, 2016). Fundamentally, the necessary condition for coexistence is that all species should show positive long-term average growth rates. Any negative difference in competitive ability (i.e. relative fitness [RF] inequality) for a common limiting resource between two species has to be compensated by stabilising growth components related to niche differentiation (ND; Chesson, 2018).
In planktonic oligotrophic conditions, we may hypothesise that a significant part of ND is related to the amount and composition of nutrients. Lower affinity for a limiting nutrient, resulting in lower competitive ability, could be compensated for by different growth performance characteristics related to fluctuating non-limiting nutrients. Accordingly, we may conjecture that if a sufficiently large water volume was experimentally confined to prevent immigration and submitted to different nutrient pulse manipulations, within the range naturally occurring in the water mass, a large part of the currently unseen species of the coexisting community should be recovered, particularly species blooming during other periods of the year. Similarly, treatments simulating more infrequent environmental conditions should likely recover rarer species.
In the Episodic Nutrient Enrichment Experiment (ENEX), we used self-filling columnar enclosures that collected the species-poor epilimnetic community of an ultraoligotrophic P-limited lake and manipulated the similar initial community of each enclosure with a variety of P and N nutrient additions to simulate fluctuations in the environment. According to the coexistence hypothesis, we expected to recover, with different treatments, most of the species that had been observed in the lake in previous seasonal studies . However, species recovery would not be enough to conclude stable coexistence. The patterns of species growing in a set of experimental treatments should agree with the requirements and predictions of the contemporary coexistence theory (HilleRisLambers et al., 2012). Several phenomenological indicators would support the long-term non-equilibrium coexistence of protist species based on the findings of our experiment: 1. Under theoretical models, RF and ND have precise definitions (Li & Chesson, 2016), which have been estimated experimentally by manipulating species densities, for instance, in plots of annual species  and planktonic cyanobacteria cultures (Gallego et al., 2019). However, because of the large number of species involved in our study, we considered diffuse competition that involved comparable interaction strengths for all pairs of species (Chesson, 2000). The RF and ND values of individual species are assessed with respect to the average fitness and ND of a community. RF refers to the scaled difference in the per capita growth rates of species in the absence of stabilisation (Adler et al., 2007). The scaling terms reflect differences in how species growth rates respond to shared limiting factors (Chesson, 2000). Consequently, the RF of a species was approximated assuming that it was proportional to its maximum growth rate in any of the treatments, which indicates an intrinsic competitive ability at some P concentration.
The ND was approximated as the negative average correlation between the growth of the target species and each of the remaining species in the complete set of treatments. Both RF and ND contribute to positive long-term average growth rates of coexisting species. Species with RF values above the community average need not show large ND values. Conversely, a higher ND should correspond to a lower RF. A small advantage in RF can be easily compensated for by a small ND.
2. Negative frequency-dependent growth rates are the hallmark of stabilising niche differences (HilleRisLambers et al., 2012), meaning that the per capita population growth declines with the relative abundance of the species in a community (Adler et al., 2007).
In our study, each species had only one initial relative density; however, we could consider that species sharing a similar nutrient niche, that is, responding preferably to the same treatments, should show, as a group, a negative frequency-dependent growth rate relationship, if they share similar stabilising mechanisms. The recovering from very low densities when episodic nutrient enrichments, of varying P and N amounts and composition, occur. initial rare species should achieve higher per capita growth rates than the most common species in nutrient manipulations. The degree of stabilisation could differ between species of different nutrient niches.
3. The stabilising ND can result from several mechanisms, including resource partitioning, density-dependent predation and storage effects (Chesson, 2000). The storage effect combines speciesspecific responses to the environment and population buffering by the presence of a persistent stage during periods of low-competitive ability. The reproductive effort is reduced, and energy investment goes to survivorship. This persistent stage can take the form of resting stages that sustain the population over poorrecruitment periods and fuel population growth during favourable periods (Li & Chesson, 2016). Many planktonic protists produce resting stages, but chrysophytes are particularly remarkable in that sense (Ellegaard & Ribeiro, 2018). Our experiment took place in a lake where chrysophytes dominated. Therefore, the patterns of cyst formation and occurrence among treatments were expected to be informative about stabilising mechanisms. 4. Competition may drive either phylogenetic clustering or overdispersion (Mayfield & Levine, 2010). Considering autotrophic organisms in the framework of the coexistence theory, high-rank taxa should cluster across the limiting resource (i.e. P) axis because P affinity is a fairly conserved trait among the main phytoplankton groups (Reynolds, 2006). In contrast, the overdispersion of lower rank taxa (i.e. species) concerning N:P stoichiometry or N sources would indicate stabilising niche differences.

| Experiment description
The experiment was performed in Lake Redon located in Central Pyrenees, at an altitude of 2,240 m a.s.l. (Figure 1a, Supporting Information Appendix S1). It is an ultraoligotrophic high-mountain deep lake (Ventura et al., 2000), in which the planktonic community is P-limited (Camarero & Catalan, 2012). We performed pulse TDP ratio of the low P and N additions (N_P) was close to the initial conditions in the lake, and the other two levels decreased the ratio (N_P+, N_P++) and increased the N imbalance (N+_P, N++_P) respectively. Each N:P combination was assessed using either ammonium (NH + 4 ) or nitrate (NO − 3 ). N was added as NH 4 Cl in five treatments and as KNO 3 in the other five. P was always added as K 2 HPO 4 . To specify the form of DIN added, the code includes an 'H' for NH + 4 additions (NH++_P, NH+_P, NH_P, NH_P+ and NH_P++) and an 'O' for NO − 3 additions (NO++_P, NO+_P, NO_P, NO_P+ and NO_P++). This information is summarised in Table S1. Photographs: Pau Giménez-Grau nutrient perturbations of the summer epilimnetic protist community using enclosures of 100 L volume, in which we could expect 10 8 -10 9 initial protist individuals. We modified P availability, N imbalance with respect to P and N source (ammonium vs. nitrate) to provide a range of conditions for evaluating competitive ability across the limiting resource (P) and potential ND across the stoichiometric variation and N supply.  Figure S1; Felip & Catalan, 2000), allowing the species to grow with less competition when limiting nutrients were added.
Under natural conditions, episodic N increases correspond to both nitrate and ammonium (Ventura et al., 2000). Consequently, all treatments were alternatively replicated using one of these N forms. The total P, nitrate and ammonium levels selected were mostly within the oligotrophic seasonal and inter-annual variations found in the lake and snow cover (Ventura et al., 2000), and above this range in a few treatments, providing mesotrophic conditions according to total P and N imbalance levels of the regions profoundly affected by N atmospheric contamination (Camarero et al., 2009). Both P and N were added to each treatment to avoid the addition of only one nutrient that could induce the other to become markedly limiting at some point during the experiment. Two should partially exclude large zooplankton (Franks, 2001). As expected, the total number of crustacean individuals in each enclosure was low. Diaptomus cyaneus and Daphnia pulicaria showed an avoidance reaction during filling and were absent or extremely low in number in the enclosures. The densities of Cyclops abyssorum in enclosures did not differ from that in the water column, but the results did not show any significant correlation with their random differences between enclosures; juveniles and adults are mainly predating on Daphnia in this lake (Ventura & Catalan, 2008). Two integrated water column samples of the epilimnion and upper metalimnion (0-20 m depth) were collected on 6 August 2013 to assess the initial protist community in the experiment.
The protist abundance and taxonomy were evaluated microscopically, and the biovolume was estimated by adjusting appropriate geometrical forms to the species shape (Hillebrand et al., 1999).
Details on productivity and seston stoichiometry from the experiment can be found in the study by Giménez-Grau et al. (2020).
Nutrients were not exhausted in any enclosure during the experiment. Treatment changes in chlorophyll and particulate C showed patterns matching the bulk biovolume changes reported here. Even with nutrient enrichment, the biomass was low. Therefore, there was little or no self-shading caused by protist growth in the enclosures.
The light environment was likely deeply modified in the enclosures, particularly because of the UV absorption by the enclosure walls.
The two enclosures with no nutrient treatment (NA) were taking account of this 'experimental bias'. As shown in the Results section, the community composition in these latter enclosures was the closest to the initial water column community. Therefore, the deviations observed in the other enclosures can be attributed with confidence to the nutrient treatments.

| Numerical methods
The growth rates for each species (i) and treatment (j) were calcu- at the end of the experiment (j) or the average value at the initial water column conditions (o), and t is the experiment duration.
Biovolume (biomass) was used instead of density because we assumed competition for resources. Biomass facilitates comparison between species with large size differences (Li & Chesson, 2016).
Species that were not detected in the initial water column samples were assumed to show an initial density at least one order of magnitude lower than the minimum density measured (other assumptions did not change the main results). As described previously, a proxy for RF was the maximum growth rate found for a species in any of the treatments. The average correlation of the growth rate of a species with the growth rates of other species, multiplied by −1, was used as a proxy for ND. In both cases, for comparison, the values were standardised to z-scores. These proxies are crude approaches to RF and ND, probably with little precision speciesby-species but likely sufficiently informative about the general relationships between the two components of coexistence in the community. Calculations and statistics were performed using the R software version 3.4.3.
The protist community response to the nutrient treatments was evaluated using multivariate regression trees (MRT; De'ath, 2002) using the Hellinger distance (Cieslak et al., 2012). The aim was to obtain a hierarchical classification of the protist community response to the treatments. The method consists of binary site clustering by repeated splitting, with each split defined by a simple rule based on environmental values, in our case nutrient additions. This rule is chosen to maximise the homogeneity of the data and to minimise the impurity (total sum of squares of the response variable values to the node mean), and each cluster (or node) represents a group of associated habitat defined by the species assemblage and environmental values. Experimental (21) and initial (2) samples were included in the analysis. The best MRT was selected using the minimal cross-validation relative error (MCVRE) criteria (Breiman et al., 1984;Tibshirani & Tibshirani, 2009) and 999-permutation tests. The mvpart package in R version 3.2.3 was used.
We used the IndVal index to determine the characteristic species of each MRT cluster at each split level. The procedure compares species fidelity and specificity to groups against randomised species distributions. We used the r package indicspecies (De Caceres et al., 2012). A 999-permutation test was used with no site combination (duleg = TRUE) and the association function 'IndVal.g'. Only species with a p-value < 0.05 were selected as indicators of the cluster. The indicator species of the same cluster can be considered as sharing a similar P and N nutrient niche.

| Community response to nutrient enrichments
Considering the experimental enclosures and the two integrated samples to assess the initial composition, a total of 151 protist taxa were distinguished: 88 autotroph species (although most of them were mixotrophs), 42 heterotrophs and 21 cysts that could not be related to any particular species (Figure 2a) , and a few to higher taxonomic levels. The number of taxa distinguished in each enclosure was between 60 and 90, similar to or higher than in the initial lake samples (50)(51)(52)(53)(54)(55)(56)(57)(58)(59)(60). Approximately 39% of the taxa were observed in the two initial water column samples. Less than 20% of the taxa were found in only one enclosure, and approximately 50% appeared at least in half of the enclosures at the end of the experiment (Figure 2).
More than 50% species exhibited increased growth at many, if not all, P additions, indicating a general P limitation. This feature was evident when the species were pooled predominantly into autotroph and heterotroph guilds (Appendix S2, Figure S2). Both trophic categories showed increasing biomass along the P enrichment gradient and no response or decline in the N imbalance gradient. Prokaryotes responded similarly ( Figure S2C). Therefore, P addition stimulated growth in the whole microbial interaction network, not only in autotrophs.
Two response patterns were found across the P enrichment gradient concerning the high-rank taxonomic groups (Figure 3; Appendix S2; Figure S3). The most substantial biovolume increase in chrysophytes and dinoflagellates occurred from low to medium P enrichment; conversely, cryptophytes, diatoms, chlorophytes and heterotrophs showed a sustained increase across the P enrichment gradient. Two response patterns of the autotroph taxa were also detected in the N imbalance experiment (Figure 3). The biomass of chrysophytes, dinoflagellates and cryptophytes declined with increasing N imbalance, whereas diatoms, chlorophytes and heterotrophs showed similar biomass values along the gradient. Beyond these general gross patterns, many species showed idiosyncratic responses, as revealed by the MCVRE-MRT.
The species response classification using the MCVRE-MRT resulted in five significant progressive splits that accounted for 55% of the total taxon variation (Figure 4). The first split was defined by the degree of P enrichment, splitting the enclosures of the two highest P additions from the rest. The second split isolated the treatments with no additions and the initial water column samples of reference.
Therefore, at this level, there were three groups with high P availability. Interestingly, the third split was no further related to P levels but to the N source present in the two highest P enrichments (ammonium or nitrate). The fourth split divided the group of the N imbalance treatments into two subgroups, one with higher ammonium addition (NH++_P and NH+_P) and the other with the highest nitrate addition (NO++_P). The fifth split divided the low ammonium and nitrate treatments, with little further variation explained. The first and the fourth splits showed the highest increase in explained variation.
Species responding similar to the treatments could be considered to share a similar nutrient niche. Therefore, the indicator taxa for each cluster (i.e. nutrient niche) at distinct levels of the tree were determined using IndVal and permutation tests. Forty-five taxa (30%) were found significant (Figure 4; Appendix S3). Taxa appearing in one enclosure or one treatment at low numbers did not pass the randomisation test. Nonetheless, some of these taxa were quite remarkable. For instance, rare litostomate ciliates (Lacrymaria and Lagynophrya), which previously had only been found during winter in the slush layers of the snow and ice cover of the lake , were recovered in high ammonium and P treatments ( Figure 2).

| Stable coexistence indicators
Most of the taxa showed positive growth in several enclosures ( Figure 2). Only 13 taxa (8.6%) showed declined growth in all F I G U R E 2 Taxa found across the experiment ordered by their maximum growth in any of the treatments. (a) Number of occurrences in the experimental enclosures. Colour bars highlight the indicator value of the taxa to a treatment (see Figure 4); NA refers to the no-addition enclosures; + and * indicate the taxa found in the seasonal Lake Redon studies of Felip, Bartumeus, et al. (1999) and Felip, Camarero, et al. (1999) respectively. (b) Taxon growth rates in each treatment. The abundance of each taxon in the initial samples is shown as a reference (black line) treatments. The frequency-dependent growth was typical among all taxa. The highest growth rates were achieved by species with initial low biovolume, irrespective of whether they were indicator taxa at some treatment level (Figure 2). This pattern was not spurious because of the use of the 'estimated' initial density values for the rare taxa during growth rate calculations. Excluding these initially 'unseen' taxa, the pattern was maintained. Among the cluster indicator taxa, both abundant and rare taxa were present at the beginning of the experiment. Therefore, the response to varying additions was not conditioned by the initial community structure.
Among autotrophic organisms, both indicator and neutral taxa showed nonlinear maximum growth rate decay with increasing initial population biomass ( Figure 5, note log-scales). When the initial population tended to zero biomass, the doubling times tended to 2-4 days, but when the biovolume was above 1,000 µm 3 /ml,

F I G U R E 3
Autotrophic protist response to P enrichment (a-e) and N imbalance (f-j) grouped by the main taxonomic groups (chrysophytes, dinoflagellates, cryptophytes, diatoms, and chlorophytes, respectively). The mean values of the lake initial conditions and the enclosures without nutrient additions are included as references. Bottom bar diagrams indicate the proportion between the groups in each of the treatments. The red and dark blue colours represent the treatments with ammonium as the N source, the orange and light blue colours represent the treatments with nitrate and green represents the enclosures with no addition (NA). LOESS tendency lines are indicated the doubling times quickly increased to several weeks or months ( Figure 5). Hence, the frequency-dependent growth occurred independently of the species nutrient niche.
The indicator species of P enrichment tended to show low ND values and higher RF values (Figure 6), irrespective of whether their N-source preference was ammonium or nitrate. This result indicates that they are better competitors with respect to the limiting nutrient and, thus, require less contribution of stabilising ND.
The indicator species of N imbalance and no-addition treatments showed higher ND values than those of P enrichment (Figure 6), particularly when RF < 0, indicating that stabilising mechanisms related to the N:P stoichiometry may favour their persistence.
All taxa presented a negative correlation between RF and ND (R = −0.246, p = 0.01).
Concerning cyst formation, a direct indicator of storage mechanisms, Dinobryon cylindricum, was the only abundant chrysophyte species that could be identified as both cyst and vegetative forms.
The degree of encystment increased in consonance with the enrichment level in both P enrichment and N imbalance treatments (Figure 7a,b). The highest cyst/vegetative cell ratios were found in the highest P enrichment and N imbalance treatments, regardless of the N source. Furthermore, there was substantial variability in the number of chrysophyte cysts in the water column of the enclosures, but on average, the density was significantly lower in the ammonium F I G U R E 4 Hierarchical classification of the protist response to nutrient additions using multivariate regression trees (MRT) and significant indicator species at five clustering levels assessed by IndVal. Autotrophic taxa are highlighted in bold characters. Asterisks indicate the splits in which a taxon was a significant indicator of a cluster. Colour intensity indicates the degree of indicator value. The accumulated amount of variance explained at each split is indicated P enrichments (Figure 7c,d). The MRT analysis indicated that the chrysophyte cysts of diameters 5 and 3 µm appeared to be related to NO_P+ and NO_P++ treatments and chrysophyte cysts S399 to NH+_P and NH++_P treatments (Figure 4), indicating nutrient ND in the chrysophyte species involved.

| The hidden community
The number and composition of the total taxa recovered from the apparently species-poor epilimnetic waters during the experiment were close to the total taxa found in a 2-year study of the ice-free period (Felip, Bartumeus, et al., 1999) and another study during the ice-covered period in Lake Redon (Felip, Camarero, et al., 1999).
Most of the species that were not observed in the initial samples were those that were found deeper in the water column during summer stratification (i.e. at the hypolimnetic deep chlorophyll maximum) or later in the season during the 1996-1997 study.
Therefore, we can infer that these species recover periodically when episodic nutrient fluctuations occur naturally throughout the year. In Lake Redon, the episodes of nutrient enrichment are related to the water mixing interaction with the sediments during the spring and autumn seasons (Ventura et al., 2000); nitrification under ice-covered and thawing conditions (Catalan, 1992); inflow of snow melting water from the catchment (Catalan, 1989); and atmospheric events of dust deposition (Camarero & Catalan, 2012). These seasonal nutrient fluctuations are common to many cold dimictic lakes throughout the Northern Hemisphere (Catalan et al., 2002).
In fact, episodic nutrient fluctuations are frequent in all lakes (Ding et al., 2017;Wilhelm & Adrian, 2008) and marine systems (Calil et al., 2011;Romero et al., 2013). Unexpectedly, we also recovered some rare litostomate ciliates that were previously observed during F I G U R E 5 Frequency-buffered growth of autotrophic taxa (n = 84). The higher the relative contribution to the total biomass of the community at the beginning of the experiment, the lower the maximum growth rate achieved in any of the treatments. Note that the values are logarithmically transformed. The linear model fitted (black line) considered intercept random effects of the several indicator clusters, which improved the fitting precision with respect to a simple linear model. However, considering random slope effects did not improve the model; thus, the slopes of the relationship for each treatment cluster (colour lines) did not significantly differ from the general response slope. Linear mixed models were fitted using the lme4 r package F I G U R E 6 Relationship between the proxies of relative fitness (RF) and niche differentiation (ND) for the autotrophic organisms in the experiment. RF was assumed to be proportional to the maximum growth rate in any of the treatments and ND to the negative average correlation of the growth rate of a given species with the rates of the rest of the species. Both variables were standardised for comparison F I G U R E 7 Chrysophyte encystment. Abundance of Dinobryon cylindricum cells (light colour) and proportion showing encystment (dark colour) in the P enrichment (a) and N imbalance (b) treatments. Abundance of total chrysophyte cysts in the P enrichment (c) and N imbalance (d) treatments winter in the slush layers of the snow and ice cover of the lake (Felip, Camarero, et al., 1999) that showed high ammonium and P concentrations, similar to the enclosure conditions where they appeared during the experiment. These organisms were phagotrophs; therefore, their association with specific nutrient conditions was either indirect through food web interactions or was related to sensory niche constrictions (Schmidt et al., 2010). In summary, we can conclude that the actual planktonic protist community is richer in species than what can be observed with current methods at any given time. Long-term species coexistence in this community will depend on the balance between RF and ND components.

| Relative fitness inequality
Nutrient addition has been traditionally used to reveal fitness differences among co-occurring autotrophic species. An increase in the supply rate of a limiting resource should lead to an increase in the abundance of the species that are most limited by the resource (Tilman et al., 1982). In our case, the protist community was clearly limited by P, and therefore, the response along the P addition gradient revealed the RF differences. Most species showed positive growth with nutrient additions (Figure 2). There was phylogenetic clustering of the main groups across the P axis (i.e. chrysophytes and dinoflagellates at lower P levels, and cryptophytes, diatoms and chlorophytes at higher P levels). This clustering pattern indicates niche conservatism related to P availability in most taxa of the same group. This feature is well known in freshwater phytoplankton (Reynolds, 2006) and has been used to develop trophic-state indicators (Phillips et al., 2013).
The contemporary coexistence theory interprets this kind of pattern as reflecting RF differences (competitive ability) or environmental filters (Mayfield & Levine, 2010). The two mechanisms are not easily distinguishable and, in fact, are not mutually exclusive. In any case, the RF differences do not promote coexistence but promote competitive exclusion (Tilman et al., 1982). For coexistence being possible, the RF inequalities should be compensated for by stabilising ND.

| Resource partitioning
In contrast to their higher competitive ability at low P concentrations, the sensitivity of chrysophytes and dinoflagellates to N imbalance indicates a stabilising niche difference that may favour coexistence between major groups in an environment of fluctuating inorganic nutrient resource stoichiometry.
In addition to N imbalance, the N source accompanying P levels can be particularly relevant as a stabilising niche difference. In our experiment, once P was available, the main factor of community differentiation was the form of N source, either nitrate or ammonium. The maximum growth rates achieved by the species and the total biomass growth were similar for the two N sources and dictated by the amount of P, but the species involved differed. Resource partitioning is one of the main mechanisms of coexistence and is independent of environmental fluctuations (Chesson, 2000). In Lake Redon and other similar P-limited lakes, nitrate is usually more abundant than ammonium as a N source. Nonetheless, large P influxes are usually associated with ammonium if the P sources are sediments (Ventura et al., 2000) or atmospheric deposition (Camarero & Catalan, 1996). Currently, the N-source preference for phytoplankton has been recognised as a complex issue that can be species-specific (Glibert et al., 2016). The MRT results showed that species do segregate according to the N source with the same P availability. The MRT clusters of high P enrichment with ammonium as the N source showed a few more indicator species than those with nitrate as the N source. Even so, all main taxonomic groups of autotrophic protists showed indicator species in both N-source clusters. In agreement with this finding, nitrate is one of the main factors explaining the variation in chrysophyte cyst composition across Pyrenean lakes (Pla et al., 2003).

| Frequency-dependent growth
Without a stabilising mechanism, any small difference in RF between species will result in the long-term exclusion of the less competitive species. The maximum per capita population growth rates achieved by autotrophic species in the experiment showed a negative relationship with their abundance in the initial lake conditions. The negative slope of this relationship can be used as a quantification of the stabilisation (Adler et al., 2007), if we assume that the species within the same nutrient guild share similar stabilising mechanisms. Interestingly, the slope was similar for species groups with a distinct nutrient response: either indicators of P enrichment with ammonium, P enrichment with nitrate, N imbalance or without specific response to treatments. Consequently, it can be concluded that there is no difference in using one or another N source for long-term sustenance in the community.
Belonging to different N guilds appears to be a stabilising mechanism for P competition. Species less competitive along the P gradient should show higher ND, which was indicated, in general terms, by the opposite tendency between the RF and ND proxies. This negative relationship need not be linear; only lower RF values need to be compensated by higher ND values. However, the frequency-buffered population growth found within each of the nutrient guilds ( Figure 5) indicated the activity of other stabilising niche differences in addition to P and N use. Such mechanisms may include partitioning relating to other nutrients and environmental conditions, storage effects and prey-specific predation.

| Storage effects
Storage effects can be the dominant coexistence mechanisms in fluctuating environments (Angert et al., 2009), and they may take many forms. The simplest persistent stage is shown by a remaining population with low competition because of the very low relative density compared to other species and the specific demographic potential of the species. Nevertheless, this remaining population must be sufficiently high in absolute numbers of individuals to be far from stochastic extinction (Lande, 1993). In a lake protist community, in which the populations oscillate over many orders of magnitude, extremely low relative populations may still include millions of individuals. Low frequency also constitutes a refuge for prey that share predators with other species (Abrams & Matsuda, 1996), thus causing an increase in survival. In summary, storage effects can be easily achieved by protist populations without specific adaptations.
Nevertheless, some protist groups show specialised resting stages.
In our study, chrysophyte cysts were particularly conspicuous.
Chrysophyte biomass did not show a significantly different response to P addition when ammonium or nitrate was the N source, although the most responsive species were not the same with different N sources. However, there were differences in the encystment rate, which were higher with nitrate. Consequently, the stabilising niche differences include not only resource partitioning but also differentiation in the way persistent stages are induced.
Unfortunately, the connection between chrysophyte vegetative forms and their cysts has not been established, except for a few species. Dinobryon cylindricum is one of them, and its cyst development has been studied experimentally, which led to the determination of a population density threshold that triggers cyst formation (Sandgren, 1983(Sandgren, , 1988. In our experiment, the encystment rate also increased at higher densities in the P enrichment treatments. At first glance, it could be interpreted as a link between encystment and intraspecific competition that buffered population growth and facilitated coexistence. However, the encystment proportion also increased in high N imbalance treatments, without high population densities. This case was a clear indication of encystment by limited environmental tolerance and, therefore, there was no reason to interpret the P encystment response differently. In fact, the predominance of chrysophytes over cryptophytes and diatoms at low P concentrations was an indication of chrysophyte niche conservatism concerning nutrients, being highly competitive at low levels but showing intolerance to high resource availability. Encystment provides a storage effect that enhances survival outside the optimal nutrient conditions as well as frequency dependence of growth.

| Predation
We found that the heterotrophic protist biomass was highly coupled to the autotrophic response to the treatments. The coupling could be mediated partially by direct interaction and partially by prokaryotic growth, in which biomass also changed coherently with the autotrophic response. The parallel autotroph-heterotroph biomass response suggested low differential predation. However, the IndVal tests revealed many heterotrophic taxa indicators of specific MRT nutrient clusters, which indicate specialised predation and, thus, a source of ND. In fact, the biomass response to P additions was similar from prokaryotes to ciliates; the P-enhanced growth of primary producers propagated to the whole microbial food web. Therefore, even when we consider guilds of a similar life-form as competitors (e.g. phytoplankton), the fact is that they are not entirely independent of other guilds (e.g. nanozooplankton), and their interactions can act as stabilising niche differences. Indeed, bacterivores include heterotrophic and autotrophic (mixotrophic) organisms, with the latter exhibiting selective predation (Ballen-Segura et al., 2017). A low overlap in predators acts as a stabilising niche mechanism. Effects of predation on coexistence does not occur at the level of the outcome of resource competition but at the same level as resources (Chesson & Kuang, 2008). Species limited solely by predation, rather than resources, can coexist if they differ sufficiently on how different predators limit them relative to how much they differ in their competitive abilities (Kuang & Chesson, 2010).

| Perspective
The coexistence of planktonic protist species cannot be seen exclusively through pairwise competition mechanisms within closed guilds; instead, it is a complex structure of diverse competitive networks that partially overlap . It includes networks of weak and strong interactions; strong interactions related to the use of limiting nutrients or shared prey might define RF differences, whereas the myriad of weak interactions (e.g. alternative N source, stoichiometric requirements and prey selection) might contribute to stabilising niche differences. The merging of network and coexistence theory seems to be the way forward (Saavedra et al., 2017). Very closely related species may show small RF differences, but they also require small ND to maintain coexistence . For instance, in highly oligotrophic conditions, such as in the mountain lake studied, the coexistence of many chrysophycean species with similar P-and N-source requirements need not be seen as a paradox. Small differences in sensitivity to physical constraints, N imbalance, cyst formation regulation, specific predators and mixotrophy level may easily balance the RF differences concerning P uptake ability. The environment and protist populations in mountain lakes fluctuate enormously, allowing species to coexist, rather than limiting them, despite P or N scarcity. More productive or extensive aquatic systems with alternating or co-limiting nutrients may sustain even richer protist communities because the average per capita growth rate differences between species over prolonged periods would be harder to achieve. Changes in the limiting resource may occur more often and environmental and competitive covariation should be lower in such aquatic systems than in oligotrophic sites, enhancing stabilising niche differences.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data supporting the results are available from Dryad Digital Repository: https://doi.org/10.5061/dryad.tx95x 69vh