Modeling the Sensitivity of Aquatic Macroinvertebrates to Chemicals Using Traits.

In this study, a trait-based macroinvertebrate sensitivity modeling tool is presented that provides two main outcomes: (1) it constructs a macroinvertebrate sensitivity ranking and, subsequently, a predictive trait model for each one of a diverse set of predefined Modes of Action (MOAs) and (2) it reveals data gaps and restrictions, helping with the direction of future research. Besides revealing taxonomic patterns of species sensitivity, we find that there was not one genus, family, or class which was most sensitive to all MOAs and that common test taxa were often not the most sensitive at all. Traits like life cycle duration and feeding mode were identified as important in explaining species sensitivity. For 71% of the species, no or incomplete trait data were available, making the lack of trait data the main obstacle in model construction. Research focus should therefore be on completing trait databases and enhancing them with finer morphological traits, focusing on the toxicodynamics of the chemical (e.g., target site distribution). Further improved sensitivity models can help with the creation of ecological scenarios by predicting the sensitivity of untested species. Through this development, our approach can help reduce animal testing and contribute toward a new predictive ecotoxicology framework.

Introduction

In the environmental risk assessment (ERA) of chemicals, it is essential to determine the environmental threshold concentration below which ecosystem structure and functioning experience no adverse impacts. In order to set this threshold, a key challenge in ERA remains the extrapolation of effects of toxicants found for a limited number of standard test species to many additional species. Ecosystems are generally populated by hundreds to thousands of species, and each species has the potential to show a different sensitivity toward one of the hundreds or even thousands of different chemical compounds that can be present in our ecosystems.[1,2] Experimental testing of this innumerable amount of species–chemical combinations, plus any possible environmentally realistic mixture of those chemicals, is impossible. We therefore need to improve current modeling approaches and make them flexible for application to any geographic region and any set of abiotic conditions.

Traditional approaches trying to incorporate species diversity into risk assessment include the application of uncertainty factors[3] and the fitting of species sensitivity distributions (SSDs)[4] to available toxicity data. While both methods are extensively used (and frequently combined), they are aimed at being protective rather than predictive and, as such, still maintain large uncertainty due to both a limited knowledge of the mechanisms underlying species sensitivity and a lack of taxonomic diversity. This lack of taxonomic diversity is especially true for uncertainty factors but also holds for SSDs. Regulatory frameworks require SSDs to contain 10 to 15 species in total[5] but divided over the different organism groups (e.g., fish, crustaceans, algae); this results in only 1 or 2 organisms from each organism group, which are often comprised out of the same set of standard test species. For the extrapolation across chemicals, Quantitative Structure–Activity Relationships (QSARs) are commonly used.[6] QSARs use chemical characteristics to predict the toxicity of many chemicals for a certain species. Due to their large demand for experimental toxicity data, however, QSAR models are often built only for specific standard test species like Daphnia magna and Oncorhynchus mykiss and, therefore, fail to account for the large species diversity of real ecosystems. In the past decade, trait-based approaches have been introduced to overcome this lack of realism.[7,8] These approaches incorporate more ecological realism into ERA by considering traits to provide a clear mechanistic link between exposure and effects, making it possible to extrapolate species sensitivities over chemicals acting by the same Mode of Action (MOA).

After the introduction of trait-based approaches as a potential tool in ERA around a decade ago,[8−11] they have been rapidly evolving.[7,12−18] Baird and Van den Brink[8] were among the first to use biological traits to predict species sensitivity. They performed a Principal Components Analysis (PCA)[19] on a species-by-substance matrix (12 macroinvertebrate species, 15 chemicals covering several MOAs), where they introduced a species-by-traits matrix as a set of nominal, passive explanatory variables. They found that up to 71% of the variability in the sensitivity could be explained with only four species traits. In a later study, Rubach et al.[12] developed the approach further, now using single and multiple linear regression instead of PCA and dividing the chemicals into groups according to their MOA. MOA has proven to be a strong determinant of species sensitivity and is, therefore, seen as a promising alternative to chemical class-based predictive toxicity modeling.[20,21] Rubach et al.[12] defined the Mode Specific Sensitivity (MSS) value of each species as the average relative sensitivity of the species to a group of chemicals with the same MOA. Single and multiple linear regressions between the MSS values and species trait data at the family level explained up to 70% of the variation in invertebrate sensitivity to three groups of insecticides. Recently, the approach by Rubach et al.[12] has been extended by incorporating relationships between species traits and sensitivity into predictive models that could potentially be used in risk assessment.[14]

So far, these predictive models have only been built for insecticides, while it remains open whether the same traits are important in explaining invertebrate sensitivity to other groups of chemicals. Therefore, in this study, we develop predictive models of macroinvertebrate sensitivity to a larger and more diverse set of MOAs. We first optimize the sensitivity prediction method of Rubach et al.[12] and Rico and Van den Brink[14] using both literature research and comparative analysis. The methods, results, and discussion of this technical optimization are given in the Supporting Information (S1). Next, and as the main focus of this article, we study the influence of MOA on species sensitivity rankings and, subsequently, on the resulting predictive trait models. By studying which traits are included in the trait-sensitivity models, we aim to get a better mechanistic understanding of how species sensitivity is determined by MOA. In a final evaluation, we test the power of species sensitivity predictions and assess if model performance depends on data availability. As a result of our research, we deliver a trait-based macroinvertebrate sensitivity modeling tool that builds predictive models for potential use in risk assessment for a diverse set of predefined MOAs.

Methods

Overall Methodology

The prediction tool was developed in the R environment and automates sensitivity rankings of aquatic macroinvertebrates and the construction of a corresponding set of trait-based models from a MOA-based chemical grouping (Figure ). A hands-on explanation of the R tool and all R scripts and databases required to run the tool are available in the Supporting Information (S1 and S2, respectively). The tool consists of four different parts: database collection, preliminary input data processing, data processing, and production of output. Each part in turn consists of multiple steps, the essentials of which will be presented briefly.

Figure 1

Structure of the developed R tool, divided into the four different parts: database collection (blue), preliminary input data processing (light gray), data processing (dark gray), and production of output (red). MSS refers to Mode Specific Sensitivity; LOOCV refers to Leave-One-Out-Cross-Validation.

Database Collection and Preliminary Input Data Processing

Several public databases were utilized to obtain data on MOA, toxicity, chemical properties, traits and taxonomy. We accept the limitations and errors of the databases exploited in this study, as it is not within the scope of this study to perform a quality assessment on the databases used.

The extensive (1213 chemicals) MOA database developed by Barron et al.[20] was used to extract a list of chemicals with their CAS numbers and their division into six groups of broad MOA (narcosis, acetylcholinesterase (AChE) inhibition, ion/osmoregulatory/circulatory (IOC) impairment, neurotoxicity, reactivity, and electron transport inhibition) and their subdivision into 31 groups of specific MOA (see Table 1 in ref (20)). The database was cleaned (e.g., chemicals for which no specific MOA was defined were removed) and prepared for further processing (e.g., specific MOAs with long names were replaced by a letter to ensure readability). For more details on preprocessing, see the corresponding R script in the SI (S2).

The U.S. Environmental Protection Agency (EPA) ECOTOX database[22] was selected as the source of toxicity data. See instructions in S1 on how to download the ECOTOX database. Only the ECOTOX tables tests, results, and species are incorporated into the R tool.

A chemical properties database was acquired by batch-running all CAS numbers available in the MOA database in the EPI (Estimation Programs Interface) Suite programs KOWWIN, MPBPVP, and WSKOWWIN, respectively, obtaining data on logKow (octanol–water partitioning coefficient); melting point, boiling point, and vapor pressure; and water solubility.[23] Modeled values were only used when experimental data were lacking. Eventually, only data on water solubility are used in the tool (as a check for realistic concentration values). All data are kept in, however, to enable future flexibility of the tool. Data on molecular weight were obtained and added to the chemical properties database by extracting SMILES (Simplified Molecular-Input Line-Entry System) of all MOA CAS numbers from the SMILECAS database (also available through EPI Suite) and calculating the molecular weight based on these SMILES using the rcdk package in R (version 3.4.5).[24]

The Tachet database provides a coding of 22 biological and ecological traits describing 472 species known to live in French freshwaters.[25,26] The database is based on “a very large and scattered published expert knowledge and diverse literature sources··· We also included unpublished observations of ourselves and colleagues.”[26] From this, we deduced that the database represents the average trait state of European species and is, therefore, suitable for our study. Since our interest lies in predicting species sensitivity and we want to avoid overfitting, only the following traits for which we could hypothesize a mechanistic relation with sensitivity were extracted:[12] maximum potential size (i.e., ≤ 0.25 cm, > 0.25–0.5 cm), life cycle duration (i.e., ≤ 1 year, > 1 year), potential number of cycles per year (i.e., < 1, > 1), dispersal mode (i.e., aquatic passive, aerial active), respiration mode (i.e., tegument, gill), feeding mode (i.e., shredder, scraper), current velocity preferendum (i.e., slow, fast), salinity preferendum (i.e., freshwater, brackish water), temperature preferendum (i.e., < 15 °C, > 15 °C), and pH preferendum (i.e., ≤ 4, > 5.5–6; see Table S1 for an overview of all the traits and corresponding trait categories).

To facilitate the cross-linking of information among the different databases, the Taxonomy database of the NCBI (National Centre for Biotechnology Information)[27,28] was used to extract the scientific names, along with the taxonomic rank and unique id, of all the species present in both the ECOTOX and the Tachet database. For this, we used the taxize package in R (version 0.9.0).[29] A copy of the ECOTOX species table and the Tachet database, both with updated taxonomy, are provided in the SI (S2).

Data Processing

Data processing has been executed once for each broad and once for each specific MOA (36 times in total). The data processing part of the tool consists of eight consecutive steps (Figure ): (i) first Mode Specific Sensitivity (MSS) calculation, (ii) removal of species without trait information, (iii) second MSS calculation, (iv) species–traits matching, (v) preparation of trait data, (vi) optimization, (vii) removal of redundant traits, and (viii) Leave-One-Out-Cross-Validation (LOOCV).

The MSS values are calculated twice: once including all species for which sufficient toxicity data are available (step i), and once only including species for which we have ascertained that also trait data are available (step iii). In between the two MSS calculations, species for which no trait data are available are removed (step ii). Performing the MSS calculation twice is necessary, because the MSS values depend on relative sensitivities and therefore are influenced by the in- or exclusion of species. The MSS calculations are implemented as described by Rubach et al.[12] In short, only toxicity tests lasting between 1 and 4 days and studying the effect of chemical stress on mortality are included in the analysis. Within each chemical, first the log transformed LC50 or EC50 values are normalized using the mean and standard deviation of all sensitivity values found for that chemical, resulting in a relative sensitivity of each species toward that chemical. Subsequently, these relative sensitivities are averaged over all chemicals belonging to the same MOA, resulting in an MSS value. Importantly, there were three fundamental differences in our method compared to the method of Rubach et al.:[12] the unit of the toxicity data (we use mol/L instead of μg/L), an extra check for realism of concentration values (should not exceed solubility), and the selection of tests with the longest exposure duration to enlarge the chance of reaching equilibrium concentration (see discussion for more extensive explanation).

After the second MSS calculation, species–traits matching is executed at the lowest taxonomic level possible (step iv). Since species level trait data are scarce, and traits measured at the genus and species levels are strongly correlated,[30] species–traits matching is done at the genus level. Prior to species–traits matching, MSS values are converted to genus level by averaging the MSS values of all species belonging to the same genus. Traits data are also converted to genus level by taking the median of the original fuzzy codes of all taxa belonging to that genus. This fuzzy coding scheme is specifically developed for describing species traits data and is ideal for differentiating species by their affinities to different trait modalities (i.e., categories) belonging to several traits (for further explanation of fuzzy codes, see ref (31)).

Next, the trait data are prepared for linear regression by expressing continuous traits (e.g., size) as weighted averages of the different trait modalities[12,14,25] and factorial traits (e.g., mode of respiration) as fixed within species (step v).[12,14] The latter means that for each of the factorial traits, the modality for which the taxon has the highest affinity is selected as the modality this taxon carries and uses throughout its entire lifespan. In the case of equal affinity to more than one trait modality (e.g., 40% gills, 40% skin respiration), a missing value is inserted. These missing values are problematic for (multiple) linear regression, because regression requires closed data sets with no missing values. To solve this, an optimization step is performed (step vi), which removes all gaps from the data set by deleting any species with missing trait values.

Next, redundant traits are removed to avoid overfitting (step vii). Traits are considered redundant (i) when they are clearly aliased with other traits (collinearity) and (ii) when they do not show enough variation in their different trait modalities. The tool tackles these two issues through (i) a collinearity maximum of 0.7 (see ref (32) for a review of different methods to deal with collinearity), removing all traits that exceed this maximum, starting with those that correlate with the largest number of other traits or, in the case of multiple traits, match this criterion with the largest exceedance of the collinearity maximum, and through (ii) a minimum on the trait modality diversity index (derived from the Shannon diversity index[33]), removing all traits with a trait modality diversity below this minimum.

Finally, a Leave-One-Out-Cross-Validation (LOOCV) is executed to quantify the predictive power of the model (step viii). LOOCV is done by successively leaving out one species from the training data set and building the MSS model based on the remaining species (for an explanation of model building, see section ). LOOCV was preferred above k-fold cross-validation, because LOOCV is better at giving a reliable picture of the true R2 than k-fold cross-validation.[34] With the produced model, the MSS value of the left-out species is predicted and afterward compared to the known value by calculating the squared error. From the LOOCV, the mean squared prediction error (MSPE) and a prediction coefficient (P2) were calculated as follows:where ŷ is the estimated MSS value, y is the calculated MSS value, and s2 is the variance of all MSS values in the data set. When all predictions perfectly match observations, P2 equals 1. A negative P2 value indicates that prediction errors exceed the total variance of the sensitivity values, and therefore the model has poor accuracy.

Production of Output

Each MOA output consists of the MSS ranking and the best MSS model. The MSS ranking is obtained by sorting the MSS values found in the first MSS calculation from low (most sensitive) to high (least sensitive). To visualize the taxonomic patterns in species sensitivity, we made heat maps of the MSS rankings, dividing the MSS values into four bins, ranging from sensitive to tolerant (MSS ≤ −1; −1 to 0; 0 to 1; ≥ 1). The best MSS model is selected in two steps. First, all possible linear models are constructed using the regsubset function of the R package leaps (version 3.0).[35] Next, the best model is selected from all constructed models based on the small sample unbiased Akaike’s Information Criterion (AICc), which takes both model fit and model complexity into consideration and which is additionally extended with a bias correction term for small sample size.[36]

Results and Discussion

Figures of all MSS rankings resulting from the first MSS calculation and tables comparing within and among broad MOA relative sensitivities can be found in S1 and S3, respectively. For the interpretation of the results, it is important to realize that the MSS values represent the relative sensitivity of a species to a group of chemicals with the same MOA. The results of the MSS rankings (section 3.1) and the MSS models (section 3.2) are not described together, because the MSS models result from a further reduced and transformed MSS data set (Figure ).

MSS Rankings

We found large differences in the MSS rankings of species depending on MOA, both within the broad MOAs (between the specific MOAs belonging to one broad MOA) and between the different broad MOAs. Differences between broad MOAs were, however, larger than differences within broad MOAs. This can be seen when comparing the standard deviation (SD) of the MSS values across the broad MOAs (0.68), with the SD of the MSS values within the broad MOAs (0.26, 0.38, 0.4, 0.52, 0.59, and 0.63 for within, respectively, the broad MOAs AChE inhibition, narcosis, IOC impairment, electron transport inhibition, neurotoxicity, and reactivity; S3). Having a higher SD between the broad MOAs compared to within the broad MOAs confirms that grouping into specific MOAs is helpful for modeling species sensitivity, as we expected.[20,21,37]

A heat map of the MSS rankings shows the general taxonomic pattern of species sensitivity and tolerance to the different MOAs (Figure ). For AChE inhibition, for instance, arthropods are most sensitive and mollusks are most tolerant. Whether crustaceans or insects are the most sensitive arthropods to AChE inhibition remains unclear, since both groups contain comparable numbers of sensitive and tolerant genera. We see similar results in two closely related studies[12,14] as well as in a review of semifield experiments.[38] The sensitivity pattern for narcotic chemicals is closely related to the sensitivity pattern of AChE inhibition, with crustacean and insect genera containing the largest number of sensitive genera. We compared this result with two outdoor microcosm studies performed with the fungicide azoxystrobin, classified as a narcotic chemical.[20] In the study of Zafar et al.,[39] significant effects were found for only one insect species (Chaoborus obscuripes) and for none of the crustacean species tested. In the study from Cole et al.[40] (summarized in ref (41)), the mollusk Sphaeriidae was the only species showing negative effects on occurrence after a constant exposure of 10 μg/L; followed by negative effects on Gammaridae, Oligochaeta, and Planorbidae at 30 μg/L; and on Asellidae at 100 μg/L. No negative effects on any of the insect groups tested (Hemiptera, Chaoboridae, Chironomidae) were found. The discrepancy between our results and the two microcosm studies might be due to the wrong assignment of narcosis as the MOA of azoxystrobin. Indeed, although two studies[20,42] classified azoxystrobin as a narcotic chemical, two other classification schemes (the QSAR Toolbox developed by the OECD, and Toxtree) assigned a non-narcotic MOA to azoxystrobin.[21] In the last paragraph of section , we discuss the causes and effects of MOA misclassification in more detail.

Figure 2

Heat map displaying taxonomic distribution of species sensitivity toward the different MOAs including (a) all species with an MSS value and (b) only species which have been tested on five or six MOAs. MSS values are divided into four bins, ranging from sensitive to tolerant. White indicates an absence of data. Upper cladogram shows similarities between MOAs, whilst the left cladogram shows the taxonomic tree.

Besides general patterns in species sensitivity, the heat map also shows that there is not one genus, family, or class which is sensitive to all MOAs (Figure a). When looking over all MOAs, the arthropod phylum contained the largest fraction of species classified as being more sensitive than average, followed by nematodes, mollusks, and annelids (Figure a). Genera belonging to the phyla Bryozoa, Cnidaria, Platyhelminthes, and Rotifera were never more sensitive than average. At the class level, Hexanauplia contained the largest fraction of genera more sensitive than average. This fits with expectations, because of their relatively small size and, herewith, large surface to volume ratio. Another result is that common test taxa are often not the most sensitive at all (Figure b). We find that Daphnia never belongs to the most sensitive group, and other commonly tested taxa (e.g., Asellus, Procambarus, and Chironomus) are for the majority of the MOAs found to be of medium tolerance. Only the commonly tested taxon Ceriodaphnia shows high sensitivity to two MOAs: AChE inhibition and reactivity.

The information presented in Figure  can also be used to identify accurately for which MOA-taxon combinations data are lacking. Ciliophora have, for instance, been studied more frequently for reactive chemicals, while Nematodes have been studied more frequently for the MOAs IOC impairment and electron transport inhibition (Figure ). At a more general level, the MOAs IOC impairment, electron transport inhibition, and reactivity have been studied to a lesser extent than the others. Especially data on insect species are missing for these MOAs, which might explain the counterintuitive result of seemingly reduced sensitivity of insects to these MOAs. Indeed, a field study performed in the 1970s shows that two mayfly (Ephemeroptera) species significantly reduced in occurrence after application of the electron transport inhibitor Antimycin A, while all other invertebrate taxonomic groups present at the study site remained unaffected.[43] This indicates that certain insect species are indeed sensitive to electron transport inhibition, but these species are not included in the ECOTOX database. That the taxonomic focus of some MOAs differs can lead to, besides a misinterpretation of results as has just been demonstrated, a bias in the sensitivity ranking and, subsequently, a bias in the sensitivity models. That the taxonomic composition of the species assemblage used to construct models is important is well-known. Maltby et al.[44] found, for example, that hazardous concentrations (HC5) derived from arthropod SSDs were significantly lower than those derived from nonarthropod invertebrates. In order to avoid any bias in predictive modeling, we should, therefore, ensure not only high taxonomic coverage but also evenness across the different taxonomic groups.

Sensitivity Models

The final modeling effort resulted in 12 significant (p < 0.05) models explaining 31 to 90% of the variation in MSS values, covering five broad MOAs and seven specific MOAs (Table ). Some of the MOAs show overlap in the traits that explained sensitivity best. For example, life cycle duration (life) and feeding mode (feeding) are included in the models for three of the six broad MOAs. This is a good indicator that these traits are in general important in explaining species sensitivity. Results from several studies[12,14,17] confirm this by including the same or closely related traits in their models. Additionally, an extensive review of potential sensitivity related traits classify feeding mode and life cycle duration as traits known to have an established link with sensitivity for several taxa.[7] Other traits selected by our modeling effect (e.g., temperature and salinity preferendum) are only included in explaining one broad MOA, which might indicate that these traits are less important for determining species sensitivity. For temperature preferendum, this is indeed confirmed by the relatively low prediction coefficient of the models including this trait (Table ). One of the models containing salinity preference (alicyclic GABA antagonism), however, has one of the highest prediction coefficients and should, therefore, be considered important. Indeed, a relationship between salinity preference and the toxicity of a GABA antagonist makes sense, because GABA is one of the most common neurotransmitters[45] and is therefore potentially influenced by the presence or absence of strategies to deal with salt stress (see ref (46) for strategies to deal with salt stress).

Table 1

Model Coefficients of the Best Models (Smallest AICc) That Were Found Significant (p ≤ 0.05) for the Different MOAs Using Exhaustive Linear Regression Analysisa

broad MOA	specific MOA	pH pref.	dispersal mode	respiration mode	life cycle dur.	life cycles year–1	feeding mode	temp. pref.	max. potential size	velocity pref.	salinity pref.	R2	P2
narcosis		–0.32	–0.42	0.51								0.33	0.013
	nonpolar				–0.77							0.42	–0.711
	polar		–0.44	0.53								0.36	0.189
neurotoxicity				–0.35		1.26	0.16				2.44	0.31	–0.29
	alicyclic GABA antagonism										6.13	0.44	0.193
AChE inhibition					–1.11	–0.75	–0.17	–0.92				0.41	–0.027
	organophosphate				–0.66		–0.15	–0.83				0.33	–0.33
	carbamate	1.16			–0.74	–1.48						0.62	–0.125
reactivity			0.46									0.67	0.255
	chromate								–0.91	–1.26		0.9	–1.611
ETIb	uncoupling oxidative phosphorylation						0.25		–0.64			0.41	–0.446
IOC impairment						–1.6				1.04		0.48	0.326

Model fit is shown as the adjusted R2 (R2), and predictive power is shown as the prediction coefficient (P2). See Table S1 for an explanation of the traits and trait modalities used in this analysis.

Electron transport inhibition.

For some traits, the correlation to sensitivity is positive for some MOAs and negative for other MOAs (e.g., dispersal, Table ). Species were more sensitive to narcosis when they preferred more neutral waters (pH preferendum), were capable of dispersing actively (dispersal), and breathed through their tegument (see Table S1 for the included the trait modalities). For reactivity, however, we find that the relationship between dispersal and sensitivity is the exact opposite, and the more passive organisms were found to be more sensitive. This can be explained by the large gap in insect genera for the MOA reactivity (Figure ). When comparing the data included in the multiple linear regression, the data set for reactivity only contains some insensitive Odonata and Diptera genera that are classified to disperse actively, while the data set for narcosis contains both sensitive Ephemeroptera and sensitive Trichoptera genera. Because of this imbalance, it remains unclear whether the direction of the trait-sensitivity relationship of the two MOAs is indeed mechanistic, or it is simply an artifact of data availability. We think that the latter is true, because the trait profile of species sensitive to narcotic chemicals matches perfectly with the trait profiles of insects belonging to the orders Ephemeroptera, Plecoptera, and Trichoptera (EPT), which are generally known as sensitive species.[47] It can be hypothesized that our trait selection does not contain the best traits, and the patterns we see only arise because traits are phylogenetically correlated, i.e., belong to trait syndromes.[47,48] Further analyses including phylogenetic information can probably further elucidate this.

Data Gaps and Potential for Improvement

Low availability of trait data has in earlier studies been described as one of the major threats to the successful implementation of trait-based approaches in ecotoxicology.[7,13] We confirm this by showing that the lack of trait data was the main obstacle in model construction. Over all the MOAs, an average of only 12% (SD ± 5%) of the species for which we had sufficient data to calculate the first MSS value were included in the construction of the MSS models. For an average of 56% (±10%) of the species, no match at all could be found in the trait database. For 15% (±5%) of the species, a match could be found in the trait database, but values on one or multiple traits of interest were lacking, and the species were therefore removed during the optimization process. It is generally known that traits databases hold a high number of gaps, and 15% of the species with missing values matches closely to 18% of the species having an incomplete trait description in a comparable trait database for Europe.[49] Due to the large loss of taxa in the analysis because of missing traits information, also only 20% (±7%) of the toxicity tests used to calculate the first MSS value ended up as an underlying data point in the multiple linear regression analysis. This indicates that a large part of the available toxicity data remain unused, merely because of insufficient trait data. Future research should therefore focus on completing and simultaneously extending the taxonomic extent of existing trait databases. This can be achieved by obtaining new data but also by combining existing trait databases, although the latter could become complicated due to existing differences in methodologies and categories between the different trait databases. Baird et al.[50] review the technological challenges in creating and sharing trait data, and Culp et al.[51] show examples of trait databases currently available for different habitats and taxonomic groups (cf. Table 2 in Culp et al.[51]).

In our analysis, 16% of the traits we deemed important regarding species sensitivity did not enter the multiple linear regression process, because they were highly correlated to another trait. This could indicate the permanent presence of collinearities between some of the sensitivity-related traits we selected (also defined as trait syndromes[48]). However, there was a clear relationship between the number of traits and the number of species going into the multiple linear regression process, indicating that collinearity was primarily found in small data sets (Figure S4). Therefore, traits were in most cases only correlating with each other due to chance, originating from the small size of the remaining data sets after all data processing steps. Besides the sensitivity-related traits selected and evaluated in this study, it can be helpful to enhance trait databases with finer morphological characteristics which are more directly related to bioaccumulation (e.g., lipid content) or to the internal distribution of the chemicals (e.g., target site distribution). Buchwalter and Luoma[52] found, for example, that although well-studied traits like body size or gill size did not explain macroinvertebrate sensitivity to heavy metals, the relative number of ionoregulatory cells did relate to dissolved metal uptake rates. Rubach et al.[53] also tried to improve trait-models by using finer morphological traits. Additionally, they measured the majority of these traits on individuals captured in the same season (although a different year) and in the same biogeographic region (sometimes even the exact same location) as they were collected for toxicity testing, which reduces the presence of trait variability due to phenotypic or seasonal intraspecific variability. However, their results show only a minor increase in model fit by measuring traits, even with the lowered intraspecific variability, both for the single linear regression and for the multiple linear regression (compare Table 2 from ref (12) with Table 3 and Figure 2 from ref (53)).

There are other aspects that could improve trait-based models, three of which we will discuss in more detail. First, measurement of internal tissue concentrations (i.e., the tissue residue) in addition to external exposure concentrations would improve toxicity predictions, because toxic effects are more closely related to the concentration inside the organism than the concentration in the water.[54,55] Although this was already proven in the 1980s, research on the (partial) incorporation of tissue residues into ERA only started the previous decade (e.g., in the biotic ligand model[56]) and has still not been accepted as a standard in chemical regulation. We hope that the recent development and increasing importance of mechanistic effect models like toxicokinetics–toxicodynamics models (e.g., the General Unified Threshold Model of Survival[57] or the Dynamic Energy Budget[58]) will lead to a more frequent measurement and publication of internal concentrations, especially since these process-based models explain the link between internal concentrations and dynamic processes underlying toxic responses, and additionally enable extrapolation from standard test conditions, e.g., to different exposure scenarios.[59−61] Rubach et al.[53] additionally demonstrated that quantitative links between the parameters of these mechanistic effect models fitted on internal concentration data and traits are substantially stronger than quantitative links between classical sensitivity end points (e.g., EC50, LC50) and traits. Measurements of internal concentrations could therefore improve our models in two ways: one, by using internal instead of water concentrations and, two, by trying to explain TKTD parameters, instead of LC50 values, which could subsequently be used to model EC50 values at any exposure profile of interest.

Second, and especially important in the absence of internal concentrations, the exposure duration of toxicity tests should be long enough to ensure approximate equilibrium between external and internal concentrations. Without going too deep into the discussion of acute versus chronic exposure, biotransformation, and elimination, we do want to emphasize that failure to obtain equilibrium concentration can result in an underestimation of effects, simply because under the same exposure time, a smaller organism reaches a higher internal concentration faster than a larger organism.[55] This means that for an acute toxicity test, different organisms may require a different minimum exposure time. We tried to account for this by selecting the toxicity test with longest exposure duration where possible. However, it can still be that some of the variability we see in Figure may be associated with differences in exposure duration (e.g., in one toxicity test, a species was exposed for 48 h to chemical A, while in another toxicity test, the same species was exposed for 96 h to chemical B).

Finally, improvements can be made by integrating Bayesian methods into the current approach. Approximate Bayesian computation can distinguish which mechanisms contribute most to patterns observed in the data[62] and can thereby help optimize model complexity. Bayesian model averaging can additionally incorporate model-selection uncertainty into statistically derived predictive models and can therefore provide a more realistic estimation of model uncertainty compared to the approach used in our study.[63,64]

Interestingly, data availability alone could not explain differences in model performance. An increase in the number of toxicity tests, chemicals, or genera did not improve model performance substantially (Figures S1, S2, and S3). Actually, model fit reduces slightly with an increase in the number of tests, chemicals, or genera. However, cross-validation results show an exact opposite trend, with a slight reduction in error with an increase in data availability. We think that model performance does not increase with data availability due to one (or a combination) of the following three reasons. First, the chemical groups (in our case, the MOAs) may insufficiently differentiate the chemicals according to the effects that they cause in invertebrates. As in any kind of grouping, mistakes or missing information could result in a wrong grouping in MOA for the chemicals (see refs (21) and (37) for studies on MOA classification errors). Additionally, MOA is not a constant property of a compound but may vary between species or life stage, depending, for instance, on the availability of target sites, exposure duration or frequency, or the end point of interest.[65] Photosynthetic inhibitors, for example, are specific toxicants toward primary producers but are often baseline toxicants toward invertebrates.[65] The chance that a chemical expresses multiple MOAs in different species becomes smaller, however, when analysis is restricted to species belonging to the same organism group (e.g., invertebrates, as in this study), because the basic biochemical systems and molecular targets affected by each MOA may be generally conserved across many species.[66] A second aspect that could have prevented increased model performance with increased data availability can be that we missed the right traits to mechanistically explain the relationship between MOA and sensitivity. As mentioned before, obtaining additional and finer morphological traits can help improve trait-based models, especially if these traits explain the toxicodynamics of the chemical, since processes related to toxicodynamics are currently not covered by the traits in our selection. Finally, we argue that the models could not achieve a better model performance because, regardless of data availability, trait data alone may be insufficient in explaining species sensitivity. Several studies show that complementing trait data with data on phylogenetics greatly enhances model performance, and that both traits and phylogenetic indicators explain a distinct part of species sensitivity.[67,68] Other attempts combining phylogenetics and physiochemical properties into predictive models have also proven successful,[1,69] although phylogenetics is still impossible to include in an inclusive approach as we performed in this study. Nevertheless, we suggest to further explore the potential of species traits for sensitivity predictions, as they enable large scale applicability and increased mechanistic understanding.

Future Direction

Our prediction tool is an addition to the fast growing and evolving science behind environmental risk assessment, which can be easily amended in the future when more data or better data processing and analysis procedures are identified. It is very flexible and facilitates (i) testing the effectiveness of different kinds of chemical grouping, e.g., based on the (physical) mechanism of action,[70] (ii) testing the predictive value of “new” traits, hypothesized to have a relationship with species sensitivity, (iii) the addition of new predictors, e.g., based on phylogenetics or physiochemical properties, and (iv) repeating the modeling exercise for different groups of organisms (e.g., fish or algae), as long as trait data are available. Analyses comparable to those reported herein would take days to weeks of preparation time to collect all input data, followed by more time to conduct and compile all analyses, requiring multiple types of commercial software. Using the software and algorithms developed for this work, these intensive efforts can now be compiled in only a few hours, already having all main open-source databases incorporated and using only R, a free software environment.

Both the sensitivity rankings and the sensitivity models produced in this study are an important step forward on the challenging road toward predictive ecotoxicology with reduced animal testing.[71] The sensitivity rankings help reduce animal testing by guiding future taxonomic focus of toxicity tests depending on the estimated MOA of the new chemicals. The sensitivity models can predict the sensitivity of species or even entire communities never tested before, avoiding additional animal tests. This allows us to rank the sensitivity of species occurring in any community composition anywhere, or in other words, to determine the worst-case ecological scenario for any aquatic ecosystem. Developing such ecological scenarios will fill one of the remaining gaps in the construction of environmental scenarios[72−74] deemed necessary for future ecological risk assessment frameworks. Together with exposure scenarios, ecological scenarios will form the basis for developing spatial-temporal explicit population-, community-, and ecosystem-level effect models for use in prospective environmental risk assessment for chemicals.