Resilience of Natural Phytoplankton Communities to Pulse Disturbances from Micropollutant Exposure and Vertical Mixing

Freshwaters are increasingly exposed to complex mixtures of pharmaceutical and personal care products (PPCPs) from municipal wastewater, which are known to alter freshwater communities’ structure and functioning. However, their interaction with other disturbances and whether their combined effects can impact ecological resilience (i.e., the ability of a system to tolerate disturbances without altering the system's original structure and processes) remain unexplored. Using in situ mesocosms in 2 lakes with different nutrient levels (mesotrophic and eutrophic), we assessed whether a pulse exposure to sublethal concentrations of 12 PPCPs affects the ecological resilience of natural phytoplankton communities that experienced an abrupt environmental change involving the destabilization of the water column through mixing. Such mixing events are predicted to increase as the effects of climate change unfold, leading to more frequent storms, which disrupt stratification in lakes and force communities to restructure. We assessed their combined effects on community metrics (biomass, species richness, and composition) and their relative resilience using 4 indicators (cross‐scale, within‐scale, aggregation length, and gap length), inferred from phytoplankton communities by discontinuity analysis. The mixing disturbance alone had negligible effects on the community metrics, but when combined with chemical contaminants significant changes were measured: reducing total biomass, species richness, and altered community composition of phytoplankton. Once these changes occurred, they persisted until the end of the experiment (day 20), when the communities’ structures from the 2 highest exposure levels diverged from the controls. The resilience indicators were not affected by PPCPs but differed significantly between lakes, with lower resilience found in the eutrophic lake. Thus, PPCPs can significantly alter community structures and reinforce mechanisms that maintain ecosystems in a “degraded state.” Environ Toxicol Chem 2019;38:2197–2208. © 2019 The Authors. Environmental Toxicology and Chemistry published by Wiley Periodicals, Inc. on behalf of SETAC.


INTRODUCTION
The growing diversity of pharmaceuticals and personal care products (PPCPs) found in freshwaters is of concern because these substances have the potential to negatively impact biota, both individually and in mixtures (Loos et al. 2009). of the water column at the expense of normal stratification patterns in lakes, thereby forcing communities to reassemble, in suboptimal conditions (Cotner et al. 2000). Freshwater ecosystems can effectively recover from these extreme events (Niemi et al. 1990;Wallace 1990), but coexposure to different stressors (including PPCPs) might compromise such coping mechanisms, thus increasing their vulnerability to multiple stressors Richmond et al. 2017). This highlights the need to assess whether the combination of pulse-like disturbances in a multistressor scenario can erode the ecological resilience of aquatic ecosystems and increase their vulnerability to anthropogenic changes (Bundschuh et al. 2017).
Ecological resilience is defined as the ability of an ecosystem to tolerate disturbances while maintaining structures, processes, and functions similar to prestress conditions (Holling 1973). This ability is finite. Once critical disturbance limits are exceeded, regime shifts may occur, causing the ecosystem to reorganize in alternative states with new sets of structures, processes, and functions (Scheffer et al. 2001;Beisner et al. 2003). The outcomes of regime shifts driven by anthropogenic stressors are uncertain and generally considered to be detrimental for biodiversity and ecosystem services (Dent et al. 2002;Folke et al. 2004;Jouffray et al. 2015). Quantifying ecological resilience is not trivial and requires holistic approaches that consider the hierarchical organization of ecosystems across multiple dimensions and non-linear behaviors (regime shifts) of ecosystems (Levin 1998;Gunderson 2000;Angeler et al. 2013aAngeler et al. , 2013b. The textural discontinuity hypothesis (Holling 1992) explicitly accounts for the complexity of ecosystems' structures and provides a frame to quantify ecological resilience (Angeler et al. 2015). This hypothesis acknowledges that scale-specific processes create distinctive structures in organism assemblages, where resources are available to support species and ecological functions (Holling 1992). These scaling patterns can be evaluated by assessing organisms' body size distributions (Nash et al. 2014;Spanbauer et al. 2016). According to the discontinuity hypothesis, species belonging to a similar body mass group operate (see conceptual Figure 1) in specific spatial and temporal scales that are fundamentally different from other groups (Holling 1992). The assembling of ecological communities takes place through the sorting of organisms with different characteristics while they strive to occupy available niches in a system (Soberón 2007). Discontinuous body size distributions are observed in many terrestrial and aquatic ecosystems, including birds (Wardwell et al. 2008), reptiles (Allen et al. 1999), mammals (Lambert 2006), and plankton and fish (Havlicek and Carpenter 2001). A few hypotheses have been formulated to explain why discontinuity occurs, which include biotic interactions (Hutchinson 1959) spanning short timescales (i.e., ecological responses) to longer timescales (including evolutionary processes; Smith et al. 2004). Nevertheless, theoretical and empirical evidence has shown that discontinuities in body size structures are stable features in ecosystems and are preserved along broad environmental gradients (Allen et al. 2005;Spanbauer et al. 2016). Strong disturbance can, however, substantially change these conservative scaling features, cause major reorganization (i.e., regime shift), and alter the discontinuity pattern (Allen et al. 2005;Allen and Holling 2008;Spanbauer et al. 2016). For instance, Spanbauer et al. (2016) showed that body size structures (inferred from discontinuity analysis) of diatoms from paleo record changed considerably in response to long-term climatic shifts that influenced hydrological conditions.
Using discontinuity analysis, we quantified 4 indicators of resilience of natural lake phytoplankton communities to the joint effect of 2 pulse-like disturbances: chemical pollution and environmental disturbance exemplified by an artificially mixing event. The 4 indicators, as illustrated in Figure 1, were 1) the number of scales present (cross-scale structure), 2) the species distribution at each scale (within-scale structure), 3) the length of species aggregates, and 4) the gap length between consecutive aggregates. Quantifying the number of scales present in a system provides an assessment of the cross-scale aspect of resilience (Peterson et al. 1998;Allen et al. 2005). Once dominant scales are identified, the distribution of species within and across scales can be evaluated. Resilience is expected to increase with increasing overlap of species performing similar functions at a single scale, and with the frequency they tend to occur across other scales (Peterson et al. 1998;Allen et al. 2005;Angeler et al. 2015). High redundancies of functions provide a robust control over disturbances by compensating for the loss of species and/or functions from multiple scales (Peterson et al. 1998;Allen et al. 2005;Wardwell et al. 2008). In addition, 2 other indicators were included: aggregation length and gap size. Spanbauer et al. (2016) found that aggregation length also changed in response to climatic regime shifts. Empirical evidence has shown that the position of species along body mass aggregates is © 2019 The Authors wileyonlinelibrary.com/ETC FIGURE 1: Conceptualization of the discontinuity approach. Species (individual circles) from the community are first ranked (log-transformed) from the smallest to the largest body size (or cell volume for phytoplankton). Discontinuities or gap and species clusters (i.e., 3 aggregates in the example) having similar size are subsequently identified. These aggregates reflect scale-specific processes, whereas the gaps represent transitional areas where resources are highly variable (Holling 1992;Stow et al. 2007 nonrandom, especially near the edges (Allen et al. 1999). Species located at the margins of aggregates are more prone to extinction, while at the same time these locations are generally preferred by invasive species (Allen et al. 1999). Combining the 4 indicators (i.e., the number of crossscale, species/functional distributions within and across scale, aggregation length, and the gap size) provides a broader characterization of the joint effects of stressors on the resilience of phytoplankton communities (e.g., Baho et al. 2014). Phytoplankton communities are particularly appealing model systems for this type of study because they respond quickly to environmental changes (McCormick and Cairns 1994) and are known to differ in community dynamics across contrasting ecosystem states, for example, eutrophication and acidification (Scheffer et al. 2001;Scheffer and Carpenter 2003;Stendera and Johnson 2008). In the present study, field experiments in 2 lakes were conducted during early summer stratification, where the communities were exposed to a single pulse of a PPCP mixture at 5 exposure levels following an artificial mixing event that disrupted the stability of the water column, which forced the communities to restructure. We tracked and assessed changes in community structures by measuring several metrics (biomass, species richness, and species composition) and the 4 resilience indicators across treatments and lakes. Resilience was quantified using an approach described by Bundschuh et al. (2017) that has until now not been applied in the context of chemical pollution. Based on current knowledge regarding the strength and intensity of disturbances required to induce a systemic change in ecosystems (Raffaelli et al. 2000;Havlicek and Carpenter 2001;Forys and Allen 2002), we adopted the null hypothesis that the combination of the 2 pulse-like disturbances, sublethal PPCP exposures and the artificial mixing event, would cause nonsignificant effects on resilience indicators even though changes in community structures are expected, reflecting differences in species' tolerance (Blanck 2002).

Site
The experiment was carried out in 2 lakes located close to Oslo, southeastern Norway. The lakes are situated a few kilometers apart but differed in nutrient levels. Lake Årungen has a surface area of 1.2 km 2 with eutrophic conditions (total phosphorus ∼30 µg/L) and receives nutrient runoff from the surrounding arable land (Sharma et al. 2008). Lake Gjersjøen has a surface area of 2.6 km 2 with mesotrophic conditions (total phosphorus typically ∼10 µg/L) and is mainly surrounded by forest (Xiao et al. 2014). The experiments were conducted in wind-sheltered areas, Årungen (59°69′68′′N, 10°73′95′′W) and Gjersjøen (59°79′67′′N, 10°77′36′′W), over 3 wk (between 9 and 29 June 2016) during which the lakes were thermally stratified (Supplemental Data, Figure S1A and B). Permission to perform the experiment was formally obtained from the urban administrations of Ås and Oppegård, respectively.

Experimental setup
Sampling phytoplankton communities and microcosms setup. Phytoplankton communities were sampled using a Limnos sampler both at the depth of maximum chlorophyll-a concentration and along a depth profile of 10 m (with 1 sample every meter) to obtain 2 respective communities: the "local community" and a "mixed community" that mimicked conditions following a stormy event (see Supplemental Data, Figure S2). Large grazers were excluded from the communities by filtration using a 60-µm nylon mesh. Once filtered, the communities (local and mixed) were transferred into 2 larger plastic carboys (50 L). Direct sunlight exposure was avoided when handling phytoplankton. Homogenized subsamples were taken from the 2 respective carboys to create 2 experimental controls (local and mixed controls) and 5 treatment levels (L.I-L.V) where only the mixed communities were exposed to 5 increasing concentrations of a PPCP mixture.
PPCP exposure levels. A mixture comprising 12 PPCP compounds (Table 1) was used, based on a previous experiment performed by Pomati et al. (2017). These compounds are commonly found in European freshwaters (Supplemental Data, Table S1). The composition and relative proportion of the PPCPs were kept constant, whereas their concentrations varied by a factor of 2.7 (roughly Euler's number) between each successive level, covering a concentration over 2 orders of magnitude (Table 1). The concentrations used represent environmentally realistic exposure scenarios (Table 1). Up to level IV, the concentrations of the individual compounds were within the range observed in European lakes and rivers. At the highest treatment level (L.V) 3 compounds-atenolol, clarithromycin, and benzophenone-4-exceeded the environmental range (Table 1). Dimethyl sulfoxide was used as a carrier solvent. Preparation and spiking of microcosms were conducted following Organisation for Economic Co-operation and Development guidelines. Chemicals used to prepare the spiking solutions were obtained from Sigma-Aldrich, ICN Biochemicals, and GlaxoSmithKline. Details regarding the PPCP mixture preparations can be found in Supplemental Data, Table S2.
Microcosms using dialysis bags. Well-homogenized subsamples of 908 mL of the local and mixed communities were taken to establish the controls (hereafter referred to as "control A" and "control B," respectively) and treatments where the mixed community was spiked inside glass beakers with the PPCPs mixture (Supplemental Data, Figure S2). Controls and treatments were set in triplicate. Treatments were spiked with 91 µL of the PPCP mixture using 5 stock solutions of increasing concentrations (L.I-L.V), whereas the controls were spiked with equivalent volumes of the carrier solvent (dimethyl sulfoxide).
Immediately after spiking, subsamples of 158 mL were collected from all microcosms to evaluate starting (day 0) conditions. The rest (750 mL) were transferred into 2.5-m-long cellulose ester dialysis bags with a certified pore size of a molecular weight cutoff ranging between 100 and 500 Da and a flat diameter of 3 cm (Spectra/por; Spectrum Europe). Bags wileyonlinelibrary.com/ETC © 2019 The Authors were closed with universal nylon clips (Spectra/por) at both ends. The pore size of the dialysis bags permits the free exchange of nutrients, gases, and chemicals with molecular weight <500 Da with the surrounding lake, while isolating phytoplankton in a realistic environmental setting (Pomati and Nizzetto 2013;Pomati et al. 2017). The behavior of the 12-PPCP mixture inside the dialysis bags had been previously analyzed and described by Pomati et al. (2017). Briefly, after spiking, the compounds diffuse out from the bags with halflives ranging from 12 to 48 h depending on the substance. After typically 3 to 7 d, the concentrations of the PPCP mixture compounds inside the bags become negligible. The exposure scenario reflected therefore a single pulse event.
Once sealed, the dialysis bags were carefully placed in a custom-built submersible protective acrylic rack (length 2.7 m, width 1.5 m, height 0.15 m, transparent to photosynthetically active radiation) and incubated at 3 and 2.5 m below the surface in Årungen and Gjersjøen, respectively. The incubation depths corresponded to the depth of maximum chlorophyll-a concentration (Gjersjøen 3 µg/L and Årungen 5 µg/L). The rack consisted of different chambers to separate dialysis bags of different treatments and was attached to a floating platform (made of PVC pipes) that was secured to the lake floor (Supplemental Data, Figure S3).
Sampling procedure and taxonomy analysis. Samples for species determination were collected and analyzed for 2 time points: at the beginning (day 0) and the end (day 20) of the experiment. To sample from the dialysis bags, the plexiglass frame was raised toward the surface and covered with a darkcolored tarp to minimize light stress during sampling. Samples were taken by cutting a predetermined length of the dialysis bags, based on the length/volume conversion (3.1 mL/cm) provided by the manufacturer, after gently massaging the dialysis bags for homogenization. Samples of 50 mL were collected from each replicate and fixed with Lugol's solution (0.5 mL) to determine taxonomical community composition.
Phytoplankton was identified using the Utermöhl technique (protocol CEN-EN 15204) and an inverted microscope. Taxa were identified to the lowest possible taxonomic unit (generally species) and biomass (mg/m³) was calculated from geometric conversions following a standard protocol (CEN-EN 16695). Information on the initial condition (day 0) was obtained by analyzing controls A and B only because we did not expect an immediate change in species composition after a few hours preceding the addition of diffuse contaminants.

Statistical analyses
Quantifying resilience. Based on taxonomy data, the relative resilience of the phytoplankton communities was quantified using the Bayesian classification and regression tree (BCART) model to identify within-and cross-scale patterns in biomass (Angeler et al. 2012). Species matrices from each replicate, consisting of log-transformed biomasses in ascending order, were created prior to the analysis. Briefly, the analysis (BCART) identifies homogenous groups of individuals with similar biomass  from sequentially splitting groups (termed "branching tree") and performs a random search considering all possible probabilistic combinations that a given split occurs (Chipman et al. 1998). These likelihoods are then ranked, and the final outcome represents the branching tree with the maximum likelihood, where the terminal nodes represent groups of maximum homogeneity (Chipman et al. 1998). Homogenous groups are assumed from ecological theory to operate in defined and distinct scaling structure (Allen et al. 2005;Angeler et al. 2012). The best models were used to assess 4 indicators of resilience defined as follows (see conceptual Figure 1): 1) the number of groups (to define the number of scales and subsequently quantify the cross-scale attributes of resilience), 2) the number of species present in each scale (used to measure the within-scale structure of resilience), 3) the length of each scale given the difference between the highest and the lowest log-transformed biomass of species belonging to a given group (aggregation length), and 4) the distance between successive aggregates or scales (gap length or discontinuities as described by Allen et al. [2005]). The software was developed by Chipman et al. (1998) and is freely available (http://www.rob-mcculloch.org/code/ CART/index.html).
Statistical comparisons. Statistical analyses were performed using R (Ver. 3.3.1) statistical programing (R Core Development Team 2017).
The effects of contaminants on phytoplankton univariate data, total biomass, and species richness were analyzed by a multiple regression analysis using treatments (i.e., the controls and the 5 exposure levels) and lake as factors. Log transformation was applied to fulfill the assumption of normality. The interaction term between treatment and lake was considered for inference. Predictions on effects of increasing the concentration of PPCPs on biomass and species richness were based on the models' coefficients. Restricted maximum likelihood was implemented in the analyses to account for the unequal number of observations between the replicates of the controls and treatments.
Nonmetric multidimensional scaling (NMDS) based on Bray-Curtis similarity was used on square root-transformed species-matrix data to assess the effects on the phytoplankton community composition. The NMDS was complemented with permutational multivariate analysis of variance (ANOVA) using Bray-Curtis similarity and square root-transformed speciesmatrix with 9999 unrestricted permutations.
Multivariate ANOVA (MANOVA) using the replicates' averages of the 4 indicators (see section Quantifying resilience) was performed using treatment and lake as main factors. Log transformation was performed in some cases to fulfill the requirements of parametric tests. The interaction term between main effects was considered. The MANOVAs were complemented with NMDS to contrast potential effects of contaminants on the combined 4 indicators (i.e., the cross-scale structure, within-scale structure, aggregation length, and gap size) of resilience (using Euclidean distance). The NMDS was supported by permutational MANOVA (using Euclidean distance, with 9999 unrestricted permutations).

Phytoplankton community metrics: Total biomass, species richness, and community composition
Differences across lakes in the absence of PPCPs. Control A and control B were compared at the beginning of the experiment (day 0) to assess the effect of the artificial mixing event. The average biomasses in the absence of mixing (control A) were 1571.3 and 533.7 mg/m 3 , whereas those observed after mixing (control B) were 1335.8 and 605.7 mg/m 3 for eutrophic lake Årungen and mesotrophic lake Gjersjøen, respectively. These differences in phytoplankton biomasses between the controls (in each lake) were not statistically significant. The mixing disturbance had a small effect on the community composition where the controls showed a large degree of overlap in ordination space at day 0 (Figure 2A), whereas interlake differences were more apparent. Eutrophic Lake Årungen had a higher total phytoplankton biomass while having numerically lower species richness (Table 2 and Figure  3) than mesotrophic Lake Gjersjøen. Species belonging to the taxonomic class Cyanophyceae (blue-green algae) were predominant in Årungen, whereas Bacillariophyceae (diatoms) were dominant in Gjersjøen (Figure 4). The controls (A and B) had comparable biovolume, species richness, and community compositions (Figure 3, Figure 4). In addition, the community compositions of the controls from the 2 lakes appeared to follow the same trajectory over time (Figure 2A).

Effects of PPCPs.
The presence of the PPCP mixture had significant effects on the total biomass and species richness of phytoplankton at the end of the experiment (Table 2 and Figure 3). The interaction term (lake × treatment) was not significant for total biomass, whereas it was significant for species richness ( Table 2). The coefficients of the multiple regression analyses showed that increasing the concentration of PPCPs significantly reduced the total biomass (β = -0.13, t value = -4.08, p < 0.01, R 2 = 0.98) and species richness (β = -1.99, t value = -4.04, p < 0.01, R 2 = 0.78) of algae. The decrease in total biomass observed in Lake Årungen was twice (on average 50% lower) than that of Lake Gjersjøen (on average 25% lower) for the highest exposure treatment (L.V).
Regardless of the interlake differences, the structure of the phytoplankton communities ( Figure 2B and C, Figure 4) changed considerably in the presence of PPCPs (permutational ANOVA [PERMANOVA]: Årungen F 6,20 = 1.45, p < 0.01; Gjersjøen F 6,20 = 2.18, p < 0.01). The biomass of most taxa decreased with increasing concentration of contaminants ( Figure 4). However, Cyanophyceae and Chrysophyceae were found to increase in relative abundance at the highest treatment level (L.V) in Gjersjøen, whereas a similar pattern was observed in Årungen for unidentified taxa ("other") and Chrysophyceae. The community compositions of the lower treatment levels (L.I and L.II) were comparable to the controls ( Figure 2B and C). The community compositions showed moderate signs of divergence from the wileyonlinelibrary.com/ETC © 2019 The Authors controls (A and B) until treatment level III ( Figure 2B and C), whereas a pronounced segregation was observed at the highest exposure treatment (L.V). By the end of the experiment, the phytoplankton community compositions of the 2 controls (A, local; B, mixed) were overlapping more in Gjersjøen than in Årungen ( Figure 2B and C).
Discontinuity analysis: Quantifying ecological resilience. The resilience of phytoplankton communities to the combined chemical pollution and mixing disturbance was evaluated using the 4 indicators: cross-scale structures, species distribution across scale, aggregation length, and gap size. The artificially induced mixing disturbance had no effect on the indicators of resilience of phytoplankton communities at the start of the experiment, whereas significant differences were found to occur across lakes (day 0; Supplemental Data, Table S3). These differences in resilience indicators across lakes persisted until the end of the experiment (day 20), when 3 out of 4 indicators were significantly different (MANOVA; Table 3). Lake Årungen had statistically lower number of scales (cross-scales) and number of species per scale with larger gap sizes than Gjersjøen, whereas aggregation length was similar ( Figure 5 and Table 3). In contrast, the treatment effect (presence of PPCPs) and the lake × treatment interactions were not significant for any of the resilience indicators ( Figure 5 and Table 3). Similarly, the results of multivariate analyses performed on the 4 combined resilience indicators of the phytoplankton communities showed negligible treatment effects (PERMANOVA: Årungen F 6,20 = 0.99, Gjersjøen F 6,20 = 098, p > 0.05). The controls and treatment levels from Gjersjøen ( Figure 6B  whereas in Årungen ( Figure 6A) the 2 highest exposure levels were moderately separated from the centroids of the controls and lower exposure levels (L.I-L.III).

DISCUSSION
In this study we assessed whether the resilience of natural phytoplankton communities was impacted by 2 pulse-like disturbances in the form of a sublethal PPCP exposure gradient and an artificial mixing disturbance using a quantitative concept outlined by Bundschuh et al. (2017) that is based on discontinuity analysis. We found that although the communities of both lakes suffered minor effects from the artificial mixing event at the start of the experiment (day 0), their ability to converge to the conditions of "undisturbed" control A by the end of the experiment (Figure 2A) was largely unaffected. These negligible effects can be attributed to 1) the fact that phytoplankton are "used" to the restructuring constraints imposed by mixing events, which are important triggers of seasonal dynamics in phytoplankton (Reynolds 1984;Diehl 2002), and 2) the absence of sharp chlorophyll-a (proxy for phytoplankton biomass) gradients along the water column when phytoplankton communities were sampled and mixed to create control B (Supplemental Data, Figure S1A and B). Nevertheless, exposure to PPCPs had a significant effect on species richness, biomass production, and community composition and prevented convergence, as observed from the controls. We observed that these effects occurred at concentration levels (L.IV and L.V; Figure 2B and C) that were at the higher end of measured environmental concentration ranges (Table 1). It is, however, critical to highlight that despite using a pulse exposure design, where the chemicals have been shown to diffuse out of the dialysis bags within 7 d (Pomati et al. 2017), the effects persisted until the end of the experiment, indicating a strong impact across multiple generations of algae. These findings ( Figures 2B and C, 3, and 4) are in line with other recent studies reporting that PPCPs induced changes at the community level of phytoplankton assemblages, including phenotypic diversity (Pomati and Nizzetto 2013;Pomati et al. 2017), community composition (Wilson et al. 2003;Richmond et al. 2016;Lee et al. 2016), and community metabolism (Wilson et al. 2004). The effects of the lower treatment levels (L.I-L.III) on the community endpoints ( Figures 2B and C exposure levels (L.IV-L.V). The applied exposure levels differed only by a factor of 2.7 between each consecutive level. Despite this small difference, a nongradual response was observed in community structures between levels III and IV. This raises the question of whether these findings underpin the existence of a threshold in the community toxic response. The 4 indicators of resilience, however, do not support this assumption. The similarities in cross-scale, within-scale, aggregation length, and gap size across treatments suggest that the resilience indicators of natural phytoplankton communities were not affected by PPCPs. Previous studies focusing on a range of different disturbances (including the impacts of nutrient enrichment, size-selective manipulations, and non-native species invasions) involving marine benthic communities (Raffaelli et al. 2000), freshwater pelagic communities (Havlicek and Carpenter 2001), and terrestrial ecosystems (Forys and Allen 2002) also found that scaling structures, inferred from body mass, marginally changed despite observing significant changes in community structure. The scaling structures (i.e., the 4 resilience indicators) reflect processes and responses occurring at a higher level of ecosystem organization and are less prone to changes induced by moderate stressors or stressors imposed over limited time (Allen et al. 2006). In contrast, fundamental elements such as community composition might vary due to differences in species' tolerance (Raffaelli et al. 2000;Havlicek and Carpenter 2001;Blanck 2002;Forys and Allen 2002;Allen et al. 2006). Combining traditional methods of characterizing community structures (e.g., species composition, size distribution) with resilience assessments provides a broader understanding of the ecological impacts of stressors (Baho et al. 2014;Angeler and Allen 2016). The lack of finer temporal analysis of the resilience indicators, which would have refined the present results, did not overly influence our inferences. It is unlikely that the resilience indicators changed and reverted in between the unaccounted sampling points (days 1-19) because once changes in scaling pattern occur they tend to persist over time (Scheffer et al. 2001;Scheffer and Carpenter 2003;Spanbauer et al. 2016). To confirm this, we compared the discontinuity analysis with an alternative method to quantify resilience (Supplemental Data, Text S1 and Figure S4), based on the scaling relationship between the size and abundance of individual phytoplankton that was measured regularly (4-d intervals). Results (Supplemental Data, Figure S4  treatments, further confirming our inferences and supporting the formulated hypothesis. The similarities in the resilience indicators across treatments do not imply that exposure to PPCPs will be without possible ecological consequences. The 2 lakes, Årungen (eutrophic) and Gjersjøen (mesotrophic), had significantly different indicators of resilience throughout the entire experimental period (Supplemental Data, Table S3). Eutrophic Lake Årungen had a significantly lower number of scales (cross-scale) and number of species per scale (within-scale), while having larger gap sizes compared to the mesotrophic Lake Gjersjøen. The lakes' contrasting trophic statuses can be considered as a spatial analogue of a regime shift from a clear water to a turbid water state (Ibelings et al. 2007). Theory predicts that significant changes in structure, function, and processes occur following a regime shift (Scheffer et al. 2001;Beisner et al. 2003). The difference in the resilience indicators across the 2 lakes may indicate that eutrophication has potentially decreased the resilience of phytoplankton in Lake Årungen and apparently led to a regime shift. The eutrophic Lake Årungen is a clear example of a lake with a poor ecological state with lower phytoplankton diversity and frequent cyanobacterial blooms (Romarheim and Riise 2009). In the present study, PPCPs at the higher dispensed levels tended to increase the difference in resilience indicators between the lakes. The presence of PPCPs can further aggravate the conditions of the eutrophic lake by reducing biodiversity (species richness; Figure 3C) and preferentially selecting certain taxa (including cyanobacteria and chrysophytes; Figure 4). A meta-analysis based on published results regarding chemical pollutants showed that PPCPs can boost the growth of cyanobacteria by suppressing other competing species of algae (Harris and Smith 2016).
In addition, the present results support the general idea that biodiversity enhances ecological stability against disturbances (Yachi and Loreau 1999;Loreau 2004;Tilman et al. 2006). The lake with the highest number of species (Lake Gjersjøen) appeared to perform better in biomass production when exposed to the contaminants (Figure 3). This stability is claimed to be mediated through overlapping of functions (functional redundancy [Yachi and Loreau 1999;Loreau 2004 loss of species or their function can be compensated by others (Peterson et al. 1998). However, a proper assessment of the distribution of functions within and across scales was not feasible because functional guilds for phytoplankton are still uncertain (Angeler et al. 2015). Such information can be beneficial in providing additional details about which functional groups are more likely to be impaired by PPCPs or by other stressors and can be easily implemented using the results of the discontinuity analysis (within-scale structure; Angeler et al. 2015).

CONCLUSIONS
Our experiment addressed the subtle ecological effects of chemical pollution at environmentally realistic exposure levels and their interaction with environmental disturbances. Our field experiments were limited to 2 lake ecosystems that differed in nutrient concentrations; thus, the interpretation of our results cannot be generalized to all aquatic ecosystems. Nonetheless, our results showed that the effect of a single pulse of PPCP superseded the effects of a single mixing disturbance on community structures, whereas the resilience indicators of phytoplankton communities were found to be only marginally affected by both stressors. Despite the lack of effects of PPCPs on resilience indicators, it was observed that within the time frame of the experiment (20 d following pulse disturbance) the community metrics from the 2 highest treatment levels did not converge toward the undisturbed control. The presence of PPCPs can therefore be potentially problematic for ecosystems that experienced previous disturbances. This implies that PPCPs can worsen the conditions of degraded ecosystem such as eutrophic lakes, thereby making restoration goals more difficult to achieve. Though the present study provides insights on the ecological effects of PPCP pulse exposure, the need for long-term continuous exposure scenarios in combination with multiple stressors persists because PPPCs are being continuously released into the environment (Richmond et al. 2017). This information can be useful to guide future assessment of the role of chemical pollution in ecosystems and might be beneficial for environmental management in a multistressor context.