Survival data analyses in ecotoxicology: critical effect concentrations, methods and models. What should we use?

hb or not hb - when and why accounting for background mortality in toxicological survival models matters?

Decisions in Environmental Risk Assessment (ERA) about impacts of chemical compounds on different species are based on critical effect indicators such as the 50% lethal concentration (LC₅₀). Regulatory documents recommend concentration-response (or concentration-effect) model fitting on standard toxicity test data to get LC₅₀ values. However, toxicokinetic-toxicodynamic (TKTD) models proved their efficiency to better exploit toxicity test data, at Tier-2 but also at Tier-1, delivering time-independent indicators. In particular, LC₅₀ values can be obtained from the reduced General Unified Threshold model of Survival (GUTS-RED) with both variants, Stochastic Death and Individual Tolerance, that include parameter h_b, the background mortality. Estimating h_b during the fitting process or not depends on studies and fitting habits, while it may strongly influence the other GUTS-RED parameters, and consequently the LC₅₀ estimate. We hypothesized that estimating h_b from all data in all replicates over time should provide more precise LC₅₀ estimates. We then explored how estimating h_b impacted: (i) GUTS-RED model parameters; (ii) goodness-of-fit criteria (fitting plot, posterior predictive check, parameter correlations); (iii) LC₅₀ accuracy and precision. We finally show that estimating h_b does not impact the LC₅₀ precision while providing more accurate and precise GUTS parameter estimates. Hence, estimating h_b would lead to a more protective ERA.

Sensitivity of phytoplankton, zooplankton and macroinvertebrates to hydrogen peroxide treatments of cyanobacterial blooms

Addition of hydrogen peroxide (H₂O₂) is a promising method to acutely suppress cyanobacterial blooms in lakes. However, a reliable H₂O₂ risk assessment to identify potential effects on non-target species is currently hampered by a lack of appropriate ecotoxicity data. The aim of the present study was therefore to quantify the responses of a wide diversity of freshwater phytoplankton, zooplankton and macroinvertebrates to H₂O₂ treatments of cyanobacterial blooms. To this end, we applied a multifaceted approach. First, we investigated the 24-h toxicity of H₂O₂ to three cyanobacteria (Planktothrix agardhii, Microcystis aeruginosa, Anabaena sp.) and 23 non-target species (six green algae, eight zooplankton and nine macroinvertebrate taxa), using EC₅₀ values based on photosynthetic yield for phytoplankton and LC₅₀ values based on mortality for the other organisms. The most sensitive species included all three cyanobacterial taxa, but also the rotifer Brachionus calyciflores and the cladocerans Ceriodaphnia dubia and Daphnia pulex. Next, the EC₅₀ and LC₅₀ values obtained from the laboratory toxicity tests were used to construct a species sensitivity distribution (SSD) for H₂O₂. Finally, the species predicted to be at risk by the SSD were compared with the responses of phytoplankton, zooplankton and macroinvertebrates to two whole-lake treatments with H₂O₂. The predictions of the laboratory-based SSD matched well with the responses of the different taxa to H₂O₂ in the lake. The first lake treatment, with a relatively low H₂O₂ concentration and short residence time, successfully suppressed cyanobacteria without major effects on non-target species. The second lake treatment had a higher H₂O₂ concentration with a longer residence time, which resulted in partial suppression of cyanobacteria, but also in a major collapse of rotifers and decreased abundance of small cladocerans. Our results thus revealed a trade-off between the successful suppression of cyanobacteria at the expense of adverse effects on part of the zooplankton community. This delicate balance strongly depends on the applied H₂O₂ dosage and may affect the decision whether to treat a lake or not.

A critical review of effect modeling for ecological risk assessment of plant protection products

A wide diversity of plant protection products (PPP) is used for crop protection leading to the contamination of soil, water, and air, which can have ecotoxicological impacts on living organisms. It is inconceivable to study the effects of each compound on each species from each compartment, experimental studies being time consuming and cost prohibitive, and animal testing having to be avoided. Therefore, numerous models are developed to assess PPP ecotoxicological effects. Our objective was to provide an overview of the modeling approaches enabling the assessment of PPP effects (including biopesticides) on the biota. Six categories of models were inventoried: (Q)SAR, DR and TKTD, population, multi-species, landscape, and mixture models. They were developed for various species (terrestrial and aquatic vertebrates and invertebrates, primary producers, micro-organisms) belonging to diverse environmental compartments, to address different goals (e.g., species sensitivity or PPP bioaccumulation assessment, ecosystem services protection). Among them, mechanistic models are increasingly recognized by EFSA for PPP regulatory risk assessment but, to date, remain not considered in notified guidance documents. The strengths and limits of the reviewed models are discussed together with improvement avenues (multigenerational effects, multiple biotic and abiotic stressors). This review also underlines a lack of model testing by means of field data and of sensitivity and uncertainty analyses. Accurate and robust modeling of PPP effects and other stressors on living organisms, from their application in the field to their functional consequences on the ecosystems at different scales of time and space, would help going toward a more sustainable management of the environment. Graphical Abstract Combination of the keyword lists composing the first bibliographic query. Columns were joined together with the logical operator AND. All keyword lists are available in Supplementary Information at https://doi.org/10.5281/zenodo.5775038 (Larras et al. 2021).

Taking full advantage of modelling to better assess environmental risk due to xenobiotics-the all-in-one facility MOSAIC

In the European Union, more than 100,000 man-made chemical substances are awaiting an environmental risk assessment (ERA). Simultaneously, ERA of these chemicals has now entered a new era requiring determination of risks for physiologically diverse species exposed to several chemicals, often in mixtures. Additionally, recent recommendations from regulatory bodies underline a crucial need for the use of mechanistic effect models, allowing assessments that are not only ecologically relevant, but also more integrative, consistent and efficient. At the individual level, toxicokinetic-toxicodynamic (TKTD) models are particularly encouraged for the regulatory assessment of pesticide-related risks on aquatic organisms. In this paper, we first briefly present a classical dose-response model to showcase the on-line MOSAIC tool, which offers all necessary services in a turnkey web platform, whatever the type of data analyzed. Secondly, we focus on the necessity to account for the time-dimension of the exposure by illustrating how MOSAIC can support a robust calculation of bioaccumulation metrics. Finally, we show how MOSAIC can be of valuable help to fully complete the EFSA workflow regarding the use of TKTD models, especially with GUTS models, providing a user-friendly interface for calibrating, validating and predicting survival over time under any time-variable exposure scenario of interest. Our conclusion proposes a few lines of thought for an easier use of modelling in ERA.

In Silico Methods for Environmental Risk Assessment: Principles, Tiered Approaches, Applications, and Future Perspectives

This chapter aims to introduce the reader to the basic principles of environmental risk assessment of chemicals and highlights the usefulness of tiered approaches within weight of evidence approaches in relation to problem formulation i.e., data availability, time and resource availability. In silico models are then introduced and include quantitative structure-activity relationship (QSAR) models, which support filling data gaps when no chemical property or ecotoxicological data are available. In addition, biologically-based models can be applied in more data rich situations and these include generic or species-specific models such as toxicokinetic-toxicodynamic models, dynamic energy budget models, physiologically based models, and models for ecosystem hazard assessment i.e. species sensitivity distributions and ultimately for landscape assessment i.e. landscape-based modeling approaches. Throughout this chapter, particular attention is given to provide practical examples supporting the application of such in silico models in real-world settings. Future perspectives are discussed to address environmental risk assessment in a more holistic manner particularly for relevant complex questions, such as the risk assessment of multiple stressors and the development of harmonized approaches to ultimately quantify the relative contribution and impact of single chemicals, multiple chemicals and multiple stressors on living organisms.

Comparing statistical analyses to estimate thresholds in ecotoxicology

Different methods are used in ecotoxicology to estimate thresholds in survival data. This paper uses Monte Carlo simulations to evaluate the accuracy of three methods (maximum likelihood (MLE) and Markov Chain Monte Carlo estimates (Bayesian) of the no-effect concentration (NEC) model and Piecewise regression) in estimating true and apparent thresholds in survival experiments with datasets having different slopes, background mortalities, and experimental designs. Datasets were generated with models that include a threshold parameter (NEC) or not (log-logistic). Accuracy was estimated using root-mean square errors (RMSEs), and RMSE ratios were used to estimate the relative improvement in accuracy by each design and method. All methods had poor performances in shallow and intermediate curves, and accuracy increased with the slope of the curve. The EC5 was generally the most accurate method to estimate true and apparent thresholds, except for steep curves with a true threshold. In that case, the EC5 underestimated the threshold, and MLE and Bayesian estimates were more accurate. In most cases, information criteria weights did not provide strong evidence in support of the true model, suggesting that identifying the true model is a difficult task. Piecewise regression was the only method where the information criteria weights had high support for the threshold model; however, the rate of spurious threshold model selection was also high. Even though thresholds are an attractive concept from a regulatory and practical point of view, threshold estimates, under the experimental conditions evaluated in this work, should be carefully used in survival analysis or when there are any biological reasons to support the existence of a threshold.

Comparison of the General Threshold Model of Survival and Dose-Response Models in Simulating the Acute Toxicity of Metals to Danio rerio

We exposed zebrafish (Danio rerio) to different concentrations of lead and cadmium, and monitored them for survival at 24, 48, 72, and 96 h. Metal toxicity was predicted and compared using the dose-response and general threshold survival models in terms of required data sets, fit performance, and applicability. Environ Toxicol Chem 2019;38:2169-2177. © 2019 SETAC.

DRomics: A Turnkey Tool to Support the Use of the Dose-Response Framework for Omics Data in Ecological Risk Assessment

Omics approaches (e.g., transcriptomics, metabolomics) are promising for ecological risk assessment (ERA) since they provide mechanistic information and early warning signals. A crucial step in the analysis of omics data is the modeling of concentration-dependency which may have different trends including monotonic (e.g., linear, exponential) or biphasic (e.g., U shape, bell shape) forms. The diversity of responses raises challenges concerning detection and modeling of significant responses and effect concentration (EC) derivation. Furthermore, handling high-throughput data sets is time-consuming and requires effective and automated processing routines. Thus, we developed an open source tool (DRomics, available as an R-package and as a web-based service) which, after elimination of molecular responses (e.g., gene expressions from microarrays) with no concentration-dependency and/or high variability, identifies the best model for concentration-response curve description. Subsequently, an EC (e.g., a benchmark dose) is estimated from each curve, and curves are classified based on their model parameters. This tool is especially dedicated to manage data obtained from an experimental design favoring a great number of tested doses rather than a great number of replicates and also to handle properly monotonic and biphasic trends. The tool finally provides restitution for a table of results that can be directly used to perform ERA approaches.

Effects of the herbicide Roundup® on the metabolic activity of Gammarus fossarum Koch, 1836 (Crustacea; Amphipoda)

Pesticides can easily reach surface waters via runoff and their potential to have detrimental impacts on freshwater organisms is high. Not much is known about how macroinvertebrates react to glyphosate contamination. In this study we investigated lethal and sublethal effects of the exposure of Gammarus fossarum to Roundup®, a glyphosate-based herbicide. The LC10 and LC50 values after 96 h were determined to be 0.65 ml/L Roundup® (230 mg/L glyphosate) and 0.96 ml/L Roundup® (340 mg/L glyphosate), respectively. As a sublethal measure of toxicity we conducted eight experiments with the feeding activity and the respiratory electron transport system (ETS) activity as endpoints. All experiments lasted seven days. Although the LC10 concentration of Roundup® was used for the feeding activity tests, 49% of the gammarids died before the end of the experiments, which is inconsistent with the calculated LC10-values. The feeding activity was significantly higher in Roundup®-enriched water (mean = 0.18 mg/mg x d) in comparison to pure spring water (mean = 0.079 mg/mg x d). No significant difference was observed between the ETS activity, which was determined after 24, 48 or 96 h after the start of the experiment, of the gammarids in Roundup® solution and in the control. The LC-values determined here are rather high, and exceed background glyphosate concentrations in most anthropogenically influenced surface waters. The increased feeding activity when exposed to Roundup® in combination with an unchanged ETS activity suggests effects on the metabolic efficiency of G. fossarum. We argue that Roundup® enhances the anabolic activity (feeding activity) in order to maintain the catabolic activity (ETS activity).

MOSAIC: a web-interface for statistical analyses in ecotoxicology

In ecotoxicology, bioassays are standardly conducted in order to measure acute or chronic effects of potentially toxic substances on reproduction, growth, and/or survival of living animals. MOSAIC, standing for MOdeling and StAtistical tools for ecotoxICology, is a user-friendly web interface dedicated to the mathematical and statistical modelling of such standard bioassay data. Its simple use makes MOSAIC a turnkey decision-making tool for ecotoxicologists and regulators. Without wasting time on extensive mathematical and statistical technicalities, users are provided with advanced and innovative methods for a valuable quantitative environmental risk assessment. MOSAIC is available at http://pbil.univ-lyon1.fr/software/mosaic/ .

New Insights to Compare and Choose TKTD Models for Survival Based on an Interlaboratory Study for Lymnaea stagnalis Exposed to Cd

Toxicokinetic-toxicodynamic (TKTD) models, as the General Unified Threshold model of Survival (GUTS), provide a consistent process-based framework compared to classical dose-response models to analyze both time and concentration-dependent data sets. However, the extent to which GUTS models (Stochastic Death (SD) and Individual Tolerance (IT)) lead to a better fitting than classical dose-response model at a given target time (TT) has poorly been investigated. Our paper highlights that GUTS estimates are generally more conservative and have a reduced uncertainty through smaller credible intervals for the studied data sets than classical TT approaches. Also, GUTS models enable estimating any x% lethal concentration at any time (LC_x,t), and provide biological information on the internal processes occurring during the experiments. While both GUTS-SD and GUTS-IT models outcompete classical TT approaches, choosing one preferentially to the other is still challenging. Indeed, the estimates of survival rate over time and LC_x,t are very close between both models, but our study also points out that the joint posterior distributions of SD model parameters are sometimes bimodal, while two parameters of the IT model seems strongly correlated. Therefore, the selection between these two models has to be supported by the experimental design and the biological objectives, and this paper provides some insights to drive this choice.

Robust Fit of Toxicokinetic-Toxicodynamic Models Using Prior Knowledge Contained in the Design of Survival Toxicity Tests

Toxicokinetics-toxicodynamic (TKTD) models have emerged as a powerful means to describe survival as a function of time and concentration in ecotoxicology. They are especially powerful to extrapolate survival observed under constant exposure conditions to survival predicted under realistic fluctuating exposure conditions. But despite their obvious benefits, these models have not yet been adopted as a standard to analyze data of survival toxicity tests. Instead simple dose-response models are still often used although they only exploit data observed at the end of the experiment. We believe a reason precluding a wider adoption of TKTD models is that available software still requires strong expertise in model fitting. In this work, we propose a fully automated fitting procedure that extracts prior knowledge on parameters of the model from the design of the toxicity test (tested concentrations and observation times). We evaluated our procedure on three experimental and 300 simulated data sets and showed that it provides robust fits of the model, both in the frequentist and the Bayesian framework, with a better robustness of the Bayesian approach for the sparsest data sets.

Combined effect of temperature and ammonia on molecular response and survival of the freshwater crustacean Gammarus pulex

Freshwater ecosystems are experiencing mounting pressures from agriculture, urbanization, and climate change, which could drastically impair aquatic biodiversity. As nutrient inputs increase and temperatures rise, ammonia (NH₃) concentration is likely to be associated with stressful temperatures. To investigate the interaction between NH₃ and temperature on aquatic invertebrate survival, we performed a factorial experiment on the survival and molecular response of Gammarus pulex, with temperature (10, 15, 20, and 25°C) and NH₃ (0, 0.5, 1, 2, 3, and 4mg NH₃/L) treatments. We observed an unexpected antagonistic interaction between temperature and NH₃ concentration, meaning survival in the 4mg NH₃/L treatment was higher at 25°C than at the control temperature of 10°C. A toxicokinetic-toxicodynamic (TK-TD) model was built to describe this antagonistic interaction. While the No Effect Concentration showed no significant variation across temperatures, the 50% lethal concentration at the end of the experiment increased from 2.7 (2.1-3.6) at 10°C to 5.5 (3.5- 23.4) mg NH₃/L at 25°C. Based on qPCR data, we associated these survival patterns to variations in the expression of the hsp70 gene, a generic biomarker of stress. However, though there was a 14-fold increase in hsp70 mRNA expression for gammarids exposed to 25°C compared to controls, NH₃ concentration had no effect on hsp70 mRNA synthesis across temperatures. Our results demonstrate that the effects of combined environmental stressors, like temperature and NH₃, may strongly differ from simple additive effects, and that stress response to temperature can actually increase resilience to nutrient pollution in some circumstances.

Basagran® induces developmental malformations and changes the bacterial community of zebrafish embryos

This study aimed to assess the effects of Basagran^® on zebrafish (Danio rerio) embryos. The embryos were exposed to Basagran^® at concentrations ranging from 120.0 to 480.6 mg/L, and the effects on embryo development (up to 96 h) and bacterial communities of 96 h-larvae were assessed. The embryo development response was time-dependent and concentration-dependent (106.35 < EC₅₀ < 421.58 mg/L). The sensitivity of embryo-related endpoints decreased as follows: blood clotting in the head and/or around the yolk sac > delay or anomaly in yolk sac absorption > change in swimming equilibrium > development of pericardial and/or yolk sac oedema > scoliosis. A PCR-DGGE analysis was used to evaluate changes in the structure, richness, evenness and diversity of bacterial communities after herbicide exposure. A herbicide-induced structural adjustment of bacterial community was observed. In this study, it was successfully demonstrated that Basagran^® affected zebrafish embryos and associated bacterial communities, showing time-dependent and concentration-dependent embryos' developmental response and structural changes in bacterial community. Thus, this work provides for the first time a complementary approach, which is useful to derive robust toxicity thresholds considering the embryo-microbiota system as a whole. The aquatic hazard assessment will be strengthened by combining current ecotoxicological tests with molecular microbiology tools.

Modelling algae-duckweed interaction under chemical pressure within a laboratory microcosm

Contaminant effects on species are generally assessed with single-species bioassays. As a consequence, interactions between species that occur in ecosystems are not taken into account. To investigate the effects of contaminants on interacting species dynamics, our study describes the functioning of a 2-L laboratory microcosm with two species, the duckweed Lemna minor and the microalgae Pseudokirchneriella subcapitata, exposed to cadmium contamination. We modelled the dynamics of both species and their interactions using a mechanistic model based on coupled ordinary differential equations. The main processes occurring in this two-species microcosm were thus formalised, including growth and settling of algae, growth of duckweeds, interspecific competition between the two species and cadmium effects. We estimated model parameters by Bayesian inference, using simultaneously all the data issued from multiple laboratory experiments specifically conducted for this study. Cadmium concentrations ranged between 0 and 50 μg·L(-1). For all parameters of our model, we obtained biologically realistic values and reasonable uncertainties. Only duckweed dynamics was affected by interspecific competition, while algal dynamics was not impaired. Growth rate of both species decreased with cadmium concentration, as well as competition intensity showing that the interspecific competition pressure on duckweed decreased with cadmium concentration. This innovative combination of mechanistic modelling and model-guided experiments was successful to understand the algae-duckweed microcosm functioning without and with contaminant. This approach appears promising to include interactions between species when studying contaminant effects on ecosystem functioning.

Constructing Time-Resolved Species Sensitivity Distributions Using a Hierarchical Toxico-Dynamic Model

Classical species sensitivity distribution (SSD) is used to assess the threat to ecological communities posed by a contaminant and derive a safe concentration. It suffers from several well-documented weaknesses regarding its ecological realism and statistical soundness. Criticism includes that SSD does not take time-dependence of the data into account, that safe concentrations obtained from SSD might not be entirely protective of the target communities, and that there are issues of statistical representativity and of uncertainty propagation from the experimental data. We present a hierarchical toxico-dynamic (TD) model to simultaneously address these weaknesses: TD models incorporate time-dependence and allow improvement of the ecological relevance of safe concentrations, while the hierarchical approach affords appropriate propagation of uncertainty from the original data. We develop this model on a published data set containing the salinity tolerance over 72 h of 217 macroinvertebrate taxa, obtained through rapid toxicity testing (RTT). The shrinkage properties of the hierarchical model prove particularly adequate for modeling inhomogeneous RTT data. Taking into account the large variability in the species response, the model fits the whole data set well. Moreover, the model predicts a time-independent safe concentration below that obtained with classical SSD at 72 h, demonstrating under-protectiveness of the classical approach.

Application of a generalized linear mixed model to analyze mixture toxicity: survival of brown trout affected by copper and zinc

Increased concerns about the toxicity of chemical mixtures have led to greater emphasis on analyzing the interactions among the mixture components based on observed effects. The authors applied a generalized linear mixed model (GLMM) to analyze survival of brown trout (Salmo trutta) acutely exposed to metal mixtures that contained copper and zinc. Compared with dominant conventional approaches based on an assumption of concentration addition and the concentration of a chemical that causes x% effect (ECx), the GLMM approach has 2 major advantages. First, binary response variables such as survival can be modeled without any transformations, and thus sample size can be taken into consideration. Second, the importance of the chemical interaction can be tested in a simple statistical manner. Through this application, the authors investigated whether the estimated concentration of the 2 metals binding to humic acid, which is assumed to be a proxy of nonspecific biotic ligand sites, provided a better prediction of survival effects than dissolved and free-ion concentrations of metals. The results suggest that the estimated concentration of metals binding to humic acid is a better predictor of survival effects, and thus the metal competition at the ligands could be an important mechanism responsible for effects of metal mixtures. Application of the GLMM (and the generalized linear model) presents an alternative or complementary approach to analyzing mixture toxicity.

Separating the effects of water physicochemistry and sediment contamination on Chironomus tepperi (Skuse) survival, growth and development: a boosted regression tree approach

More comprehensive ecological risk assessment procedures are needed as the unprecedented rate of anthropogenic disturbances to aquatic ecosystems continues. Identifying the effects of pollutants on aquatic ecosystems is difficult, requiring the individual and joint effects of a range of natural and anthropogenic factors to be isolated, often via the analysis of large, complicated datasets. Ecotoxicologists have traditionally used multiple regression to analyse such datasets, but there are inherent problems with this approach and a need to consider other potentially more suitable methods. Sediment pollution can cause a range of negative effects on aquatic animals, and these are used as the basis for toxicity bioassays to measure the biological impact of pollution and the success of remediation efforts. However, experimental artefacts can also lead to sediments being incorrectly classed as toxic in such studies. Understanding the influence of potentially confounding factors will help more accurate assessments of sediment pollution. In this study, we analysed standardised sediment bioassays conducted using the chironomid Chironomus tepperi, with the aim of modelling the impact of sediment toxicants and water physico-chemistry on four endpoints (survival, growth, median emergence day, and number of emerging adults). We used boosted regression trees (BRT), a method that has a number of advantages over multiple regression, to model bioassay endpoints as a function of water chemistry, sediment quality and underlying geology. Endpoints were generally influenced most strongly by water quality parameters and nutrients, although some metals negatively influenced emergence endpoints. Sub-lethal endpoints were generally better predicted than lethal endpoints; median emergence day was the most sensitive endpoint examined in this study, while the number of emerging adults was the least sensitive. We tested our modelling results by experimentally manipulating sediment and observing the impact on C. tepperi endpoints. For survival, experimental observations were accurately predicted by models, which highlighted the importance of conductivity and dissolved oxygen for this endpoint. In comparison, experimental median emergence day was poorly modelled, most likely due to the influence of a wider range of predictors identified as being important influences on this endpoint in models. To demonstrate how BRT model results compare to more traditional techniques, we analysed survival data using multiple regression. Both models yielded similar results, but boosted regression trees offer important advantages over multiple regression. Our results illustrate how boosted regression trees can be used to analyse complex ecotoxicological datasets, and reinforces the importance of water chemistry in sediment toxicology.