Hidden drivers of low-dose pharmaceutical pollutant mixtures revealed by the novel GSA-QHTS screening method

See allHide authors and affiliations

Science Advances  07 Sep 2016:
Vol. 2, no. 9, e1601272
DOI: 10.1126/sciadv.1601272


The ecological impacts of emerging pollutants such as pharmaceuticals are not well understood. The lack of experimental approaches for the identification of pollutant effects in realistic settings (that is, low doses, complex mixtures, and variable environmental conditions) supports the widespread perception that these effects are often unpredictable. To address this, we developed a novel screening method (GSA-QHTS) that couples the computational power of global sensitivity analysis (GSA) with the experimental efficiency of quantitative high-throughput screening (QHTS). We present a case study where GSA-QHTS allowed for the identification of the main pharmaceutical pollutants (and their interactions), driving biological effects of low-dose complex mixtures at the microbial population level. The QHTS experiments involved the integrated analysis of nearly 2700 observations from an array of 180 unique low-dose mixtures, representing the most complex and data-rich experimental mixture effect assessment of main pharmaceutical pollutants to date. An ecological scaling-up experiment confirmed that this subset of pollutants also affects typical freshwater microbial community assemblages. Contrary to our expectations and challenging established scientific opinion, the bioactivity of the mixtures was not predicted by the null mixture models, and the main drivers that were identified by GSA-QHTS were overlooked by the current effect assessment scheme. Our results suggest that current chemical effect assessment methods overlook a substantial number of ecologically dangerous chemical pollutants and introduce a new operational framework for their systematic identification.

  • Pollution
  • pharmaceuticals
  • freshwaters
  • microorganisms
  • complexity
  • low-dose
  • mixtures
  • GSA
  • QHTS


How can the effects of long-term exposure to low concentrations of pharmaceutical and personal care product (PPCP) mixtures be assessed on nontarget organisms? This question was recently ranked as number 1 among the 22 scientific priorities regarding PPCPs in the environment as part of a “big question exercise” hosted by the Society of Environmental Toxicology and Chemistry (SETAC), involving more than 500 environmental scientists from 57 countries (1). PPCPs are widely released into the aquatic environment, where they mix at low concentrations over extended periods of time via diverse pathways, resulting in a very complex fate (13). Despite their ubiquity, the effects of PPCPs in the environment are not well understood (1, 2, 4). Part of the difficulty in the elucidation of the risks associated with PPCP pollution arises from the different nature of their effects with respect to classical, acutely toxic pollutants. PPCPs consistently produce sublethal effects even at low concentrations (2, 47) but usually do not present clear evidence of lethality even at unrealistically high doses (2, 4, 6). Sublethal effects are important because they eventually affect natural systems via direct or cascading effects (69). The complexity of the fate of PPCPs and their sublethal effects pose important challenges to our present understanding of pollutant mixture effects (7, 10, 11). Null additive models such as concentration addition (CA) are usually considered to be a reliable hypothesis to accurate mixture effect predictions (1214). CA and its related dose-response theory do not consider low-dose nonlinear/nonadditive sublethal effects (15, 16). The current assumption that effects cannot be understood below a certain threshold, usually the 10 to 20% relative effect termed the “gray zone” (13), has limited the study of the sublethal region of the dose-response curves. Preliminary evidence indicates that significant sublethal effects may occur at doses well below the gray zone, for both individual chemicals and mixtures (5, 10, 1619), although these findings are criticized because of large uncertainties in existing experimental methods (13, 16, 20).

To address the current limitations in the study of sublethal effects of low-dose PPCP mixtures, we propose a new tool consisting of global sensitivity analysis coupled with quantitative high-throughput screening (GSA-QHTS). GSA-QHTS provides a novel perspective for chemical effect assessments by taking advantage of the family of computational GSA techniques to guide experimental design and data analysis of QHTS laboratory experiments (see Fig. 1). In general, GSA apportions the observed variability of the system response (output) onto the system’s drivers (inputs), in terms of both direct (first-order) and interaction (higher-order) effects. This quantifies the relative importance (sensitivity) of the system’s inputs without any a priori assumption on the nature of the system’s internal processes [such as additivity, linearity, or threshold effects (21, 22)]. Among GSA techniques, the elementary effects (EE) method (23) relies on a small number of samples and provides easy-to-interpret statistics that describe the direct effects of input factors (μ and μ*) and their interactions (σ) (21, 22). Here, GSA is used to produce cost-effective (parsimonious) experimental templates that reasonably match present QHTS capacities used in ecotoxicology. As a representative freshwater biological system, we used a high-throughput configuration of a bioluminescent whole-cell biosensor that detects metabolic toxicity based on a freshwater cyanobacterium (Anabaena sp. PCC7120 CPB4337; hereinafter A. CPB4337) (24). Cyanobacteria are focal primary producers in aquatic systems that take part in carbon and nitrogen cycling (25). The EE method was applied to generate a set of realistic low-dose mixtures of commonly found PPCPs in freshwater and characterize their relative importance and interactions (if any). Light intensity was also included in the analysis to illustrate the ability of GSA-QHTS to account for not only chemical factors but also combined effects with other environmental stressors. Finally, novel results of factor importance and interactions obtained with the GSA-QHTS framework were critically compared with those from classical additive effect mixture models and validated on upscaled experimental freshwater benthic microbial communities.

Fig. 1 GSA-QHTS experimental framework.

The novel approach couples GSA with QHTS (GSA-QHTS) to understand the main effects and interactions of combinations of diverse input factors, such as chemicals, biotic or abiotic factors, etc. Here, GSA-QHTS was applied to study the effects on A. CPB4337 of mixtures of 16 PPCPs (C1 to C16) at environmentally realistic low doses (D1 to D3) and the influence of light intensity, an abiotic factor. The steps of the framework are highlighted: parsimonious GSA sampling (to generate an experimental design template), QHTS (with randomized biological replication), and GSA screening of the importance and interactions of input factors controlling biological response.


Generating a low-dose PPCP mixture experimental design based on the EE method

For a total of 17 input factors (16 PPCPs and light intensity), the EE sampling produced an array of 180 unique mixtures. The selected set of PPCPs included most families of pharmaceutical pollutants occurring broadly in freshwaters, such as antibiotics, lipid regulators, stimulants, analgesics, hypertension regulators, and psychiatric drugs (26, 27). Individual PPCP doses in the mixtures ranged from 25.5 to 40,780.0 ng liter−1, differing by nearly three orders of magnitude along their three selected environmentally realistic discrete dose levels (median of means, mean of maxima, and maximum of maxima; Fig. 2A). The frequency of input factor levels in the mixtures was near uniform (Fig. 2B), honoring the discrete uniform (DU) sampling probability distribution used in GSA sampling (see Materials and Methods). In addition, mean dose distribution of the 16 PPCPs along the 180 mixtures occurred in the nanogram per liter to microgram per liter range, defining a robust exposure-oriented array of environmentally realistic dose combinations (Fig. 2C). The total PPCP doses of the mixtures presented a bimodal frequency distribution (Fig. 2D) with a minimum total dose of 25,490 ng liter−1 and a maximum dose of 123,500 ng liter−1, indicating a slight bias toward the maximum values of environmental occurrence (26, 27). This is a drawback of the interlevel uniformity of dose distance required to optimize GSA screening computation. However, it was considered appropriate based on the short-term exposure duration (24 hours) of the QHTS experiment involving A. CPB4337.

Fig. 2 Low-dose PPCP mixture experimental design template.

(A) Discrete dose levels (in nanograms per liter) of each of the 16 PPCPs (C1 to C16) included in the low-dose mixtures. The three discrete levels were the median of means, the mean of maxima, and the maximum of maxima of each PPCP in freshwater. (B) Sampling frequency of the three discrete levels for the 16 PPCPs (C1 to C16) across the 180 mixtures. (C) Mean dose of each PPCP along the 180 low-dose mixtures. (D) Frequency distribution of sum of PPCP doses (in nanograms per liter) within each of the 180 low-dose mixtures.

Exposure to low doses of PPCPs produced significant sublethal effects

The QHTS results showed statistically significant alterations of bioluminescence [surrogate end point of metabolic toxicity (24, 28)] of A. CPB4337 [Fig. 3A; one-way analysis of variance (ANOVA), P < 0.001]. The response of the organism to both PPCPs alone and in mixtures was significantly different [one-way ANOVA with Tukey’s honest significant difference (HSD), P < 0.001] (Fig. 3B). The median response of A. CPB4337 after exposure to individual PPCPs shifted to an increased bioluminescence [hormesis (15)], whereas exposure to PPCP mixtures produced a median response toward a reduction of the bioluminescence signal (metabolic toxicity). The median effect of a substantial number of mixtures (67 of 180) was significantly different from the control median response as shown by the non-overlapping of the bootstrapped 95% confidence intervals (CIs) (Fig. 3C) (29). Overall, the observed sublethal mixture exposure to PPCPs was in the gray zone of the dose-response curves (median and maximum bioluminescence inhibition of 6.25 and 21.5%, respectively) (13). The magnitude of the observed mixture sublethal effects (lower than 30% effect; Fig. 3B) is in agreement with a growing body of literature reporting significant individual or mixture sublethal effects of PPCPs at low realistic doses (nanograms per liter to micrograms per liter) in a variety of organisms (5, 6, 10, 18, 19, 3033). The nonchemical factor light intensity did not affect individual exposure to PPCPs (one-way ANOVA, P = 0.132) (Fig. 3C) but did so in mixture exposures (P = 0.000198) (Fig. 3D), indicating that light intensity may be a factor that plays a significant role in mixture effects of PPCPs to A. CPB4337. Light intensity is a critical factor not only for the metabolic activity of primary producers but also for the environmental fate, transformation, and bioactivity of many chemicals including PPCPs (3436). However, the interaction of light, similar to that of other environmental factors, on the ecological effects of PPCPs and other chemical pollutants remains an unresolved issue in terms of global significance (11, 37).

Fig. 3 Exposure to low doses of PPCPs produced significant sublethal effects.

(A) Notched box plot for relative bioluminescence of A. CPB4337 control observations (n = 880), individual PPCP exposure (Ind) (n = 768), and exposure to mixtures (Mix) (n = 1080). (B) Median relative bioluminescence level of A. CPB4337 exposed to the 180 mixtures with respect to control levels (relative bioluminescence, 1). Mixtures are ranked in decreasing order based on median relative bioluminescence values. Vertical lines are bootstrapped (n = 999) 95% CIs for the median. Horizontal red dashed lines are bootstrapped 95% CIs for control median relative bioluminescence level (n = 999). The most potent mixture, Mix 16, is highlighted. (C and D) Notched box plots sorted by light intensity (L1 and L2) for relative bioluminescence of A. CPB4337 exposed to individual PPCPs (n = 768) and mixtures of PPCPs (n = 1080), respectively. The notches extend to ±1.58 IQR (interquartile range)/n1/2, where no overlapping of notches among boxes offers evidence of statistically significant differences among their medians (71). Statistically significant differences were tested by one-way ANOVA: ***P < 0.001.

The classical additive mixture approach does not predict the sublethal effects of low-dose PPCP mixtures

According to Loewe additivity (38), the foundation of the mixture CA null model (13), each mixture component contributes to the total mixture effect with fractions of doses that can be summed up if the shapes of each mixture component’s individual dose-response curve are known (39). A key theoretical assumption of CA is that infinitely low doses result in infinitely decreasing effects. Therefore, mixture components below a statistical threshold of effect (no observed effect concentration) can still contribute to mixture effects but in a linear and additive way (13). An important practical limitation of CA is that only chemicals individually displaying dose-response curves can be considered in the model (40, 41). A first-tier, high-dose (10 mg liter−1) benchmarking of the toxicity of the 16 PPCPs was performed to identify candidate PPCPs for dose-response curve derivation (Fig. 4A). From the 16 PPCPs, only two antibiotics, C10 (erythromycin) and C11 (ofloxacin), resulted in strong (near 90%) bioluminescence inhibition of A. CPB4337. C14 (venlafaxine) and C8 (ketoprofen) resulted in weak (less than 20%) inhibition and hormesis, respectively (one-way ANOVA with Dunnett’s test, P < 0.05). Dose-response curves were derived for PPCPs exhibiting inhibition (C10, C11, and C14; Fig. 4B), and the rest of the PPCPs were ignored for further mixture effect modeling based on CA (14). On the basis of the unique C10/C11/C14 ratio in each of the 180 low-dose mixtures, predicted mixture dose-response curves were modeled (Fig. 4C), and additive mixture effects were predicted and compared to the observed experimental mixture effects (Fig. 4D). CA predicted nearly no effect (from 0 to 0.026% inhibition) for all the 180 mixtures in contrast with the observed range of experimental effects (inset in Fig. 4D). The departure of data from the 1:1 reference line (Fig. 4D) confirms the inability of the additive null models to predict the observed effects (NSE, −1.766), where an NSE of 1 represents perfect agreement (42).

Fig. 4 The null additive mixture models do not predict the sublethal effects of low-dose PPCP mixtures.

(A) Individual effects of PPCPs (C1 to C16) on A. CPB4337 relative bioluminescence at 10 mg liter−1. Horizontal lines represent the control’s median, Q25 and Q75 relative bioluminescence. Significance values: *P < 0.05 and ***P < 0.001. (B) Dose-response curves for factors C10 (erythromycin), C11 (ofloxacin), and C14 (venlafaxine). Lines and symbols represent fitted nonlinear models (five-parameter log-logistic models) and experimental data, respectively. Shaded areas represent 95% CIs of model predictions. (C) Modeled mixture dose-response curves for 21 unique C10/C11/C14 ratios present in the sampled 180 PPCP mixtures. The box plot shows the 180 low-dose mixtures’ dose range. (D) Observed versus predicted 1:1 plot (42) for the 180 low-dose mixture effects as predicted by the CA additivity model. Note that only the contributions of C10 (erythromycin), C11 (ofloxacin), and C14 (venlafaxine) are considered by the models. NSE, Nash-Sutcliffe efficiency coefficient with 95% CI in brackets. The inset shows the distribution of observed and predicted luminescence values for the 180 mixtures.

GSA-QHTS identified hidden drivers of low-dose mixture sublethal effects

The graphical representation of GSA statistics (μ, σ) or (μ*, σ) in a Cartesian plane provides a ranking of input factor importance or direct effects (that is, separation from origin along the μ axis) and interaction effects (that is, separation from origin along the σ axis). According to the general distribution pattern of the input factors on the μ*-σ and μ-σ Cartesian planes, systems can be generally classified as linear/additive, mixed, and nonlinear/nonadditive (Fig. 5A) (23). The transition to a nonlinear/nonadditive system is denoted by all the input factors systematically scoring above the 45° line in the μ*-σ plot (43) and in between the μ ± 2σ/Embedded Image reference lines in the μ-σ plot (right panels in Fig. 5A) (23). Thus, it is clear from Fig. 5 (B and C) that the response to low-dose mixtures of PPCPs was nonlinear/nonadditive in our experiments. The comparison between direct EEs calculated as average of absolute values (μ*; Fig. 5B) and those with their individual signs (μ; Fig. 5C) indicates that the effects of some of the PPCPs were also nonmonotonic (that is, could increase or decrease bioluminescence consecutively at increasing doses). For example, C6 in Fig. 5B is shown among the group of the most important factors (rightmost in x axis) based on absolute value effects (μ*), whereas in Fig. 5C, it is shown close to 0 (no direct effects) based on the arithmetic mean of effects (μ), that is, ± effects cancel out in the mean. According to their ranked μ* direct effects (Fig. 5D), only half of the chemicals were most important in controlling low-dose sublethal effects of the PPCP mixtures, and one of them (C13) was negligible. Half of the chemicals exhibited an overall effect of bioluminescence inhibition (μ < 0), and the other half exhibited an overall effect of bioluminescence induction or hormesis (μ > 0).

Fig. 5 GSA-QHTS is able to characterize global drivers of low-dose mixture sublethal effects.

(A) GSA schematic representation of input factor distributions along μ*-σ and μ-σ Cartesian planes for linear/additive, mixed, and nonlinear/nonadditive systems. (B and C) GSA results in the μ*-σ and μ-σ Cartesian plots, respectively, for the 17 input factors [16 PPCPs (C1 to C16) and light intensity (L)] analyzed in the present study. Red lines indicate the limits for linear additive systems. For clarity, only important and nonimportant input factors are identified in the figure. Hollow symbols indicate important factors (larger separation from the μ*-σ or μ-σ plane origin). (D) Ranked input factors by importance (μ*) showing the proposed limits (red dashed lines) for important, moderately important, and nonimportant factors.

As a complementary line of evidence, two-way ANOVA was performed for each group of input factors (most important, C16, C3, C10, C5, C4, C14, C1, and C6; less important, C11, C7, L, C12, C9, C8, C15, and C13). In the most important factors group, seven of the eight factors were found to participate in first- or second-order significant terms (ANOVA, P < 0.05), whereas no significant term was found to participate in any first- or second-order term for the less important factors group (section 2.2.4 in S1). Although ANOVA only explored first- and second-order terms, the results corroborate the ranking of factor importance and interactions identified by EEs, where the latter also simultaneously consider higher-order interactions (22).

Both the EE method and two-way ANOVA confirm that there was structured information in the sublethal effects observed with the 180 low-dose PPCP mixtures and not just stochastic noise as assumed by the gray zone hypothesis (13). In addition, the EE method revealed (and ANOVA confirmed) that the nature of sublethal effects from low-dose PPCP mixtures was nonlinear/nonadditive, which explains why additive null models for chemical mixtures failed in their predictions (Fig. 4D). According to Fig. 4A, only 2 of the 16 PPCPs [erythromycin (C10) and ofloxacin (C11)] and, to a lesser extent, venlafaxine (C14) were potential candidates to drive any low-dose mixture effect. However, despite being important factors (especially erythromycin, which ranked third in global importance according to GSA-QHTS; Fig. 5C), seven other PPCPs were found to be as (if not more) important in low-dose mixtures (Fig. 5D). These mixture drivers remained hidden to the high-dose benchmarking (Fig. 4A) and therefore hinder the predictions of the null additive models (Fig. 4D). These results indicate that a simple linear/additive chemical risk assessment approach may overlook a substantial number of chemicals that may be critical under real low-dose environmental conditions. GSA-QHTS advances preliminary observations, suggesting nonlinear/nonadditive effects resulting from low-dose mixtures of PPCPs (5, 10, 33). The new framework adds substantial value by allowing a robust, global, and comparative understanding of chemical drivers and their effects that remained defective until now (1, 5, 10, 44). Finally, the EE method was effective when discriminating the importance of chemical factors and other nonchemical stressors. For example, despite the apparent significant contribution of light to low-dose PPCP mixture effects (Fig. 3D), its global importance was negligible (Fig. 5D). This provides a direct comparison of diverse factors in a common scale metric, addressing an unresolved, yet urgent, demand in ecotoxicology (11, 45, 46).

Important drivers of mixture bioactivity identified by GSA-QHTS are also relevant at higher ecological complexity scales

We performed an ecological scaling-up mesocosm experiment to assess the generalizability of the factors identified as important by the GSA-QHTS (Fig. 6). The most potent mixture found in the low-dose mixture experiment (mixture number 16, hereafter Mix 16; Fig. 3B) was selected with two variations (Mix 16-4 and Mix 16/10; see Materials and Methods). Mix 16-4 lacked the four most important PPCPs identified by the EE method (C16, C3, C10, and C5; see Fig. 5D), whereas Mix 16/10 was a 10-fold dilution of Mix 16. The effect of these three low-dose mixtures was evaluated in a series of community-level metabolic end points (see Materials and Methods) (47). Multivariate generalized linear models (GLMMV) (48) revealed a statistically significant interaction between the factors treatment (control, Mix 16, Mix 16-4, and Mix 16/10) and time (model 1 in Table 1; ANODEV, Wald P = 0.03) but no effects from treatment or time separately. Changes in model communities were expressed mainly on Yeff and β-Glu activity (section 4 in S1), indicating relevant effects on both the autotrophic and heterotrophic components of the microbial community. Mixture effects were unraveled by fitting individual GLMMV (models 2 to 5 in Table 1) with time as an explanatory factor. We found that time was significant in control and Mix 16-4 treatments (ANODEV, Wald P = 0.01 and 0.04, respectively) but not in Mix 16 and Mix 16/10 (P = 0.54 and P = 0.61, respectively), indicating that the four most important PPCPs in the full Mix 16 blocked the temporal evolution of the metabolic processes of the freshwater benthic communities. The insensitivity to time reflects reduced dynamic behavior and increased temporal autocorrelation, common symptoms of microbial community stress, lethality, and ecological shifts (4951). The removal of the four important PPCPs was more effective in preventing effects than a 10-fold dilution (Mix 16/10), although they accounted for less than 3% of the total PPCP dose of the mixture (Fig. 6). These results confirm that sublethal effects at the population level can drive disruptive effects at the community level (8, 10).

Fig. 6 Ecological scaling-up experiment.

(A) River benthic microbial community inocula were obtained from a nearby unpolluted stream in the Llémena River (Girona, Spain). Schematic benthic microbial community modified from the study of Egan et al. (72). (B) A set of cobbles was scraped for their microbial communities, which were used to colonize rough glass substrata under laboratory conditions. (C) Experimental microbial communities were exposed to three low-dose PPCP mixtures (Mix 16, Mix 16-4, and Mix 16/10; see Materials and Methods). (D) Selected community-level end points (F0, F, Ymax, Yeff, β-Glu, and Phos; see Materials and Methods) covered both autotrophic and heterotrophic global fitness indicators and were monitored as a function of time.

Table 1 ANODEV table summary of GLMMV models.

GLMMV models were fitted to the response of experimental freshwater benthic communities exposed to selected PPCP low-dose mixtures and sequential analysis of deviance (ANODEV) was performed. MV data used to build each GLMMV include the six community-level metabolic end points measured: F0, the dark-adapted basal fluorescence; F, the light-adapted steady state fluorescence; Ymax, the maximum photosynthetic efficiency of photosystem II (PSII); Yeff, the effective quantum yield of PSII; β-Glu, β-glucosidase; and Phos, alkaline phosphatase. The experiment included two factors: treatment (four levels) and time of exposure (three levels). Treatment levels were as follows: control (n = 5); Mix 16 (n = 3), PPCP mixture 16; Mix 16-4 (n = 3), mixture 16 without the four most important PPCPs from GSA results; and Mix 16/10 (n = 3), mixture 16 diluted 10 times. Time levels were as follows: 24, 36, and 120 hours of exposure for each treatment (n = 4). Therefore, the total number of samples was n = 42. The null hypothesis (H0) for the Wald test is that the reduction in model residual deviance is 0. Res. df, residual degrees of freedom; df diff., degrees of freedom difference added by the sequential inclusion of each term of the model.

View this table:


Main findings

Contrary to current scientific practice (1214), typical high-dose benchmarking failed to identify the complete set of low-dose mixture effect drivers, and the bioactivity of the mixtures was not predicted by the null additive mixture models. The novel GSA-QHTS framework was effective at screening the importance of PPCP pollutants at environmentally realistic low-dose mixtures and at identifying the main drivers of PPCP sublethal effects. GSA-QHTS revealed that the observed mixture effects were not random but driven by a specific subset of PPCPs via nonlinear and/or nonadditive effects. This confirms previous evidence of unexpected low-dose effects of PPCPs (5, 10, 18, 19, 52) and offers an alternative explanation to the source of the mixture effects typically unaccounted for by additivity assumptions in environmental samples (53). The ecological scaling-up experiment illustrated how these hidden PPCPs (identified by GSA-QHTS at the population level) were also responsible for disruptive effects at the community level. Because real-world exposure to chemical pollution of natural systems is often dominated by low-dose complex mixtures combined with other biotic and abiotic stressors, we may underestimate the actual risk of chemical pollution with current methods.

Limitations and further research

The results of this work are limited to the specific biological system used as a model and provide a first-tier phenomenological screening for the effect of low-dose PPCPs in freshwater systems. Further research should include other biological systems to get a wider picture of the generalizability of our findings. An identification of the mechanisms behind the effects of the focal pharmaceuticals requires extension of the GSA-QHTS to receptor-mediated biological tests designed to unravel mechanistic pathways of stressor effects (54). GSA results are sensitive to a meaningful a priori and realistic selection of both input factors and ranges, or GSA may result in misleading conclusions (21, 22). A limitation of the present application of the GSA-QHTS is that the EE method, because of its reduced sampling intensity, provides a qualitative ranking of the relative importance of the drivers studied. As life sciences’ high-throughput capacity improves, alternative quantitative GSA variance decomposition methods (5557) that require a larger number of samples can be implemented. This will allow quantitative apportioning of the total variance of the effects onto the individual or combined input factors studied (21, 22). Armed with the identification of main drivers provided by the GSA-QHTS approach and mechanistic understanding based on receptor-mediated analysis, further research should include the development of new conceptual and modeling approaches to predict and integrate nonlinear/nonadditive effects of chemical pollution within the “big picture” of large-scale ecological systems (11, 46).

Broader impacts

As a reliable and robust approach, GSA-QHTS could be of interest to scientists designing experiments that test the effects of combined stressors, not only chemical but also physical and biological, therefore crossing the border of different disciplines. The present work has foundational implications not only for ecotoxicology but also for biological and medical sciences because it provides a new experimental framework to approach complexity under realistic conditions.



PPCPs and dose levels were selected on the basis of a meta-analysis of PPCP occurrence in Spanish freshwaters (details are in S2). The 16 PPCPs included in the present study encompassed (i) antibiotics: erythromycin (ERY, C10; >95%; Fluka) and ofloxacin (OFLO, C11; 99.9%; Sigma Aldrich); (ii) β-blockers: atenolol (ATE, C1; 91%; Sigma Aldrich), bezafibrate (BEZA, C2; ≥98%; Sigma Aldrich), and gemfibrozil (GEM, C4; >99%; Sigma Aldrich); (iii) stimulants: nicotine (NICO, C12; 99%; Sigma Aldrich) and caffeine (CAFE, C13; >99%; Sigma Aldrich); (iv) analgesics: ketoprofen (KETO, C8; >98%; TCI), ibuprofen (IBU, C7; ≥98%; Sigma Aldrich), paracetamol (PARA, C9; 99%; Sigma Aldrich), and diclofenac [DICLO, C6; ≥98% (sodium salt); Sigma Aldrich]; (v) diuretics: furosemide (FURO, C3; ≥98%; Sigma Aldrich) and hydrochlorothiazide (HYDRO, C5; ≥99%; Sigma Aldrich); (vi) an antidepressant: venlafaxine (VENLA, C14; >98%; TCI); and (vii) anticonvulsants: carbamazepine (CARBA, C16; ≥98%; Sigma Aldrich) and primidone (PRIMI, C15; ≥98%; Sigma Aldrich). The Chemical Abstracts Service number and pharmacological family are summarized in S3. Additional information on preparation, handling of stock solutions, and chemical stability can be found in S2.

Global sensitivity analysis coupled with quantitative high-throughput screening

The GSA technique known as the EE method, or the “Morris method” (2123), was used as the GSA framework (Fig. 1) to (i) guide the experimental plan in establishing a parsimonious set of k = 17 input factors (mixtures and experimental light condition) to perform QHTS (GSA sampling) and (ii) compute EE sensitivity screening measures μ (or μ*) and σ (GSA screening). The EE method (23) scores the relative importance and interactions of k input factors according to the marginal changes (EEs) that they produce in the output variable (relative bioluminescence) when they are changed one at a time at discrete levels with all the other factors present (21, 22). The EE method is qualitative in nature and can be applied to assess the relative importance of each of the k input factors in the context of all other factors varying at the same time. Additional information on the EE method is summarized in previous studies (21, 22, 58).

GSA sampling

To compute EEs, Morris (23) proposed to sample t random trajectories across the k-dimensional space of the input factors, varying one factor at a time at r discrete levels within the input factor probability distributions. Typically, t = 10 produces satisfactory results (59). The number of experimental units (N) to perform the EE analysis is given byEmbedded Image(1)

In our case, GSA sampling was performed for a total of 17 input factors (16 pharmaceutical pollutants and 1 abiotic factor, light intensity), resulting in N = 10 (17 + 1) = 180 experimental units. As discrete levels for the input factors, three doses (“D” in Fig. 1) of each PPCP were selected on the basis of statistical descriptors (median of means, mean of maxima, and maximum of maxima) of their environmental concentrations. Light intensity (L) was fixed at two physiological levels (see QHTS below). DU probability distribution functions were used to represent input factor level probabilities. An enhanced sampling for uniformity [eSU (58)] method with an oversampling size of 1000 was used to optimize the sampling quality (uniformity, spread, and time). The experimental design template of the 180 experimental units used in the QHTS experiment can be found in S4.

QHTS: Bioactivity determinations based on A. CPB4337

Once the experimental template is defined by GSA sampling, a biological receptor is tested for each of the 180 mixtures with the corresponding light intensity (Fig. 1). A. CPB4337 was used as a biological receptor (24, 60). A. CPB4337 is a genetically engineered whole-cell cyanobacterial biosensor that constitutively produces bioluminescence (24). Changes in the bioluminescence signal were used as a surrogate end point of metabolic toxicity (24, 28). Standard growth and maintenance conditions were followed (24, 60). Exposure experiments were performed in continuously shaken, 24-well, transparent plastic microtiter plates (60). The treatments assigned to each plate and well were totally randomized for each independent experiment, ensuring that each plate contained 6 control wells and 18 treatment wells. A total of six independent experiments were performed, with M = 11 plates per experiment. This resulted in treatments with interexperimental and intraexperimental replications of 6 and 1, respectively, and controls with interexperimental and intraexperimental replications of 6 and at least 6M, respectively. Plates were incubated for 24 hours in physiological light conditions: L1 (~35 μmol [photons] m−2 s−1) and L2 (~65 μmol [photons] m−2 s−1). After 24 hours, 100 μl of each biological sample was transferred to 96-well, white microtiter plates for bioluminescence measurements (60). Bioluminescence was normalized to the mean control response of each individual plate (n = 6) and expressed as relative bioluminescence.

GSA screening

For each factor (xi, i = 1…k) and level (j = 1…r), QHTS generated an experimental output value Y. A factor’s EE (uj,i) is computed asEmbedded Image(2)

The mean and SD of the EEs ui over the r levels for each factor produce three sensitivity measurements: μ, the mean of the EEs; μ*, the mean of the absolute value of the EEs; and σ, the SD of the EEs. μ* and μ estimate the overall direct (first-order) effect of a factor. Both statistics are calculated because μ accounts for the directionality of the effects (positive or negative) but can suffer from compensation of opposing sign effects, leading to apparent low importance, whereas the absolute values in μ* compensate for this artifact but do not provide directionality. σ estimates the higher-order characteristics of the parameter. Although EEs (ui) are local sensitivity measurements, their moments μ, μ*, and σ are considered global measurements because they are computed across the variation space of all input factors, defined by their probability distributions (21, 22).

Individual effect of PPCPs and mixture effect prediction based on the additive CA model

The bioactivity of each of the 16 PPCPs at the two different light conditions was evaluated individually at each of the three dose levels defined by their statistical descriptors (median of means, mean of maxima, and maximum of maxima) (see QHTS: Bioactivity determinations based on A. CPB4337 section). In addition, their bioactivity was also evaluated at 10 mg liter−1 at the two different light conditions as a first-tier toxicity screening. Complete dose-response curves were assessed for erythromycin (C10), ofloxacin (C11), and venlafaxine (C14). Six independent experiments were performed with randomized experimental designs, interexperimental and intraexperimental replication schemes, experimental conditions, and bioluminescence measurements identical to those explained for QHTS before. Mixture effect predictions were computed on the basis of the well-known predictive arrangement of the CA equation as described by Faust et al. (39) (see details in section 3.3 of S1).

Effects of selected mixtures on experimental benthic microbial communities

An overview of the ecological scaling-up experiment involving experimental benthic microbial communities can be found in Fig. 6. River benthic microbial community inocula were obtained from a nearby unpolluted stream (Llémena River, Girona, Spain). A set of 10 to 12 cobbles was collected from a riffle area and scraped for their biofilm microbial communities. The inocula were used to colonize rough glass substrata under laboratory conditions for 15 days. Nutrients, light, and temperature were controlled (see S2 for details). Microbial communities were exposed for 7 days to three low-dose PPCP mixtures (Mix 16, Mix 16-4, and Mix 16/10). Mix 16 was the most potent of the 180 mixtures (see section 2.2.1 of S1). Mix 16-4 was identical to Mix 16 in components and doses except without the PPCPs C16, C3, C10 and C5, whereas Mix 16/10 was a 10-fold dilution of Mix 16. The effect of the mixtures was evaluated on a series of community-level metabolic end points (the photosynthetic parameters F0, the dark-adapted basal fluorescence; F, the light-adapted steady-state fluorescence; Ymax, the maximum photosynthetic efficiency of PSII; and Yeff, the effective quantum yield of PSII; as well as the extracellular enzymatic activities β-Glu and Phos) that covered both autotrophic and heterotrophic global fitness indicators suited to study the effects of chemical pollution on freshwater benthic microbial communities (47, 61). See S2 for further details on experimental methods.

Experimental design and statistical analysis

Basic quality control of data and data analysis were performed in R version 3.1.2 and RStudio version 0.98.1091 (see section 1.2.1 of S1 for details). For the QHTS experiments, sample sizes for treatments were fixed to n = 6. ANOVA was used to study difference of means. A two-sided test was used with a minimum testing level for statistical significance (α) set to 0.05, with exact P values reported in the text and figures. When multiple testing was performed, Tukey’s HSD method with a correction for unbalance experimental design (62) was used. When multiple testing was limited to the comparison of treatments with a single reference, the correction of Dunnett’s test was used (63). Power (1 − β) in all ANOVAs was >0.8 for relevant effect sizes (see S1 for details). Bootstrapped 95% CIs were calculated by the bias-corrected and accelerated interval method (see section 2.2.1 of S1 for details) (64). EE method computation was performed using the eSU approach (58) with a publicly available MATLAB toolbox (see section 2.2.3 of S1 for details) (65). Dose-response curves and derived parameters were computed using the drc package for R (66). The mixture effect predictions according to the CA model were computed using a freely available R script developed by the authors (see section 3.4 of S1 for details) (67). Analysis of CA model fitness with observed data was performed using the publicly available FITEVAL software implemented in MATLAB (42, 68). For the model benthic microbial community experiments, sample sizes for treatments were fixed to n = 3. Time-dependent MV responses of the freshwater benthic experimental communities were analyzed using GLMMV (48) with the “mvabund” package for R, version 3.10.4 (69). Statistical significance of individual or nested GLMMV model fits was evaluated by ANODEV (70), and MV test statistics were constructed using a Wald statistic (69), where P values were approximated by resampling rows using residual permutation as described by Szöcs et al. (48) (see section 4 of S1 for details). Power was, in general, >0.76 for relevant effect sizes. Extended information on the analysis performed as well as a nonexhaustive list of R packages used in the study can be found in S1.


Supplementary material for this article is available at

S1. Description of data sets and guide for data analysis

S2. Supplementary methods

S3. Table summary of pharmaceutically active pollutants

S4. Composition of the 180 mixtures

table S1.1. Summary of observations by Case and Exp. Group for the GSA.csv data set.

table S1.2. Summary of the variables included in the GSA.csv data set.

table S1.3. Summary of observations by Case and Exp. Group for GSA_MIX_MEDIAN data set.

table S1.4. Summary of the variables included in the GSA_MIX_MEDIAN.csv data set.

table S1.5. Summary of observations by Exp. Group for ANTIBIOTIC_EC50.csv data set.

table S1.6. ANTIBIOTIC_EC50.csv data set.

table S1.7. Summary of observations by Exp. Group for biofilm.csv data set.

table S1.8. biofilm.csv data set.

table S1.9. P values of Wald test for significance of time variable of the ANODEV test for each univariate GLMs fitted to each end point (Phos, Gluc, F0, Ymax, F, Yeff) for each treatment with time as explanatory variable.

table S2.1. Number of observations for calculating each occurrence of statistical descriptors.

table S2.2. Stability of PPCPs under experimental conditions.

fig. S1.1. Q-Q plot for normal distribution for control observations (n = 880).

fig. S1.2. Notched box plot for Lum of control observations by Experiment.

fig. S1.3. hovPlots for bioluminescence as a function of the Experiment for Control observations.

fig. S1.4. hovPlots for bioluminescence as a function of the Experiment for Treatment observations.

fig. S1.5. Notched box plot for Lum of controls (n = 880), individual exposure to PPCPs (n = 768), and exposure to mixtures of PPCPs (n = 1080).

fig. S1.6. 95% family-wise confidence level for the difference of means according to anova.2_glht model.

fig. S1.7. Notched box plot for Lum of control observations (n = 880), divided into two aleatory groups C1 (n = 440) and C2 (n = 1080).

fig. S1.8. Notched box plot for Lum for treatment observations shorted by L level (1, 2).

fig. S1.9. Boxplots of Lum as a function of PPCPs (C1 to C16).

fig. S1.10. Histogram (left panel), normal Q-Q plot (central panel) and jack after jackknife plot (right panel) for bootstrapped medians of Lum values (R = 999) for “treatment” = 1.

fig. S1.11. Bootstrapped (R = 999) medians and 95% CIs for the 180 treatments.

fig. S1.12. Tolerance level in the estimation of median values for lum for the 180 mixture effects.

fig. S1.13. EE μ*-σ plot for the 17 studied input factors (16 PPCPs and light intensity).

fig. S1.14. Ranked input factors by importance (μ*).

fig. S1.15. EE μ*-σ plot for the 17 studied input factors (16 PPCPs and light intensity).

fig. S1.16. Diagnostic plots of Anov.1 model.

fig. S1.17. Diagnostic plots of Anov.4 model.

fig. S1.18. Dose ranges of PPCPs (in nanograms per liter).

fig. S1.19. Histogram of the frequency distribution of the total sum of the 16 PPCPs in the 180 mixtures.

fig. S1.20. Rose plot presenting the relative abundance (counts) of each PPCP (variable) and level (value) of each PPCP in the 180 mixtures.

fig. S1.21. Rose plot presenting the mean concentration (value) of each PPCP (variable) in the 180 mixtures.

fig. S1.22. Dose-response data and fitted LL5 drm models (chem.1) to the experimental data.

fig. S1.23. Schematic representation of the calculation sequence required to predict chemical mixture effects according to CA model.

fig. S1.24. The 21 dose-response models (LL.5 models) fitted to the in silico–predicted dose-response patterns of the 21 unique combinations of C10, C11, and C14.

fig. S1.25. Observed versus predicted bioluminescence values of A. CPB4337 to the 180 low-dose mixtures of PPCPs.

fig. S1.26. FITEVAL report for the predicted versus experimental low-dose mixture effects of PPCPs.

fig. S1.27. Response of six community-level end points measured along time as a function of treatment on model freshwater benthic microbial communities.

fig. S1.28. Residual versus fits plot to check the quadratic mean-variance assumption of negative binomial regression (with different metabolic end points coded in different colors).

References (73109)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank Y. Khare for the software used for eSU that was applied in the study and A. L. Barran-Berdon for assistance with Fig. 6B. Funding: This research was supported by grants CTM2013-45775-C2-1-R, CTM2013-45775-C2-2-R, and CGL2010-15675 from MINECO; by the Dirección General de Universidades e Investigación de la Comunidad de Madrid; by Research Networks grants S2013/MAE-2716 and 0505/AMB-0395; and by the NEREUS COST Action ES1403 of the EU Framework Programme Horizon 2020. I.R.-P. is a Ramón Areces Foundation Postdoctoral Fellow. M.G.-P. holds an FPI grant from MINECO. R.M.-C. acknowledges support as a Fellow of the University of Florida Water Institute. I.R.-P. acknowledges support from the Ramón Areces Foundation (Becas Fundación Ramón Areces para Estudios Postdoctorales. XXVII Convocatoria de Becas para Ampliación de Estudios en el Extranjero en Ciencias de la Vida y de la Materia). Author contributions: I.R.-P., R.M.-C., R.R., and F.F.-P. identified the research question. I.R.-P., R.M.-C., F.F.-P., M.G.-P., S.G., and F.L. formulated the hypothesis. R.M.-C., I.R.-P., M.G.-P., S.G., and S.S. developed the experimental design. S.G., I.R.-P., M.G.-P., S.S., and M.C. planned the experiments. S.G., M.G.-P., I.R.-P., and M.C. performed the experiments. S.G., M.G.-P., M.C., and I.R.-P. processed the data and compiled the data sets. I.R.-P., M.G.-P., and R.M.-C. analyzed the data. I.R.-P. and R.M.-C. prepared graphical material. I.R.-P., M.G.-P., R.M.-C., F.F.-P., S.S., and R.R. discussed the results. I.R.-P., F.F.-P., and R.M.-C. wrote the paper. F.F.-P., R.M.-C., S.S., R.R., and F.L. revised the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: Data sets, R scripts, and additional material are freely available at Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article