Dose Addition in the Induction of Craniofacial Malformations in Zebrafish Embryos Exposed to a Complex Mixture of Food-Relevant Chemicals with Dissimilar Modes of Action

Abstract

Background: Humans are exposed to combinations of chemicals. In cumulative risk assessment (CRA), regulatory bodies such as the European Food Safety Authority consider dose addition as a default and sufficiently conservative approach. The principle of dose addition was confirmed previously for inducing craniofacial malformations in zebrafish embryos in binary mixtures of chemicals with either similar or dissimilar modes of action (MOAs). Objectives: In this study, we explored a workflow to select and experimentally test multiple compounds as a complex mixture with each of the compounds at or below its no observed adverse effect level (NOAEL), in the same zebrafish embryo model. Methods: Selection of candidate compounds that potentially induce craniofacial malformations was done using in silico methods-structural similarity, molecular docking, and quantitative structure-activity relationships-applied to a database of chemicals relevant for oral exposure in humans via food (EuroMix inventory, n = 1,598). A final subselection was made manually to represent different regulatory fields (e.g., food additives, industrial chemicals, plant protection products), different chemical families, and different MOAs. Results: A final selection of eight compounds was examined in the zebrafish embryo model, and craniofacial malformations were observed in embryos exposed to each of the compounds, thus confirming the developmental toxicity as predicted by the in silico methods. When exposed to a mixture of the eight compounds, each at its NOAEL, substantial craniofacial malformations were observed; according to a dose-response analysis, even embryos exposed to a 7-fold dilution of this mixture still exhibited a slight abnormal phenotype. The cumulative effect of the compounds in the mixture was in accordance with dose addition (added doses of the individual compounds after adjustment for relative potencies), despite different MOAs of the compounds involved. Discussion: This case study of a complex mixture inducing craniofacial malformations in zebrafish embryos shows that dose addition can adequately predicted the cumulative effect of a mixture of multiple substances at low doses, irrespective of the (expected) MOA. The applied workflow may be useful as an approach for CRA in general. https://doi.org/10.1289/EHP9888

Headnote

Background: Humans are exposed to combinations of chemicals. In cumulative risk assessment (CRA), regulatory bodies such as the European Food Safety Authority consider dose addition as a default and sufficiently conservative approach. The principle of dose addition was confirmed previously for inducing craniofacial malformations in zebrafish embryos in binary mixtures of chemicals with either similar or dissimilar modes of action (MOAs).

Objectives: In this study, we explored a workflow to select and experimentally test multiple compounds as a complex mixture with each of the compounds at or below its no observed adverse effect level (NOAEL), in the same zebrafish embryo model.

Methods: Selection of candidate compounds that potentially induce craniofacial malformations was done using in silico methods-structural similarity, molecular docking, and quantitative structure-activity relationships-applied to a database of chemicals relevant for oral exposure in humans via food (EuroMix inventory, n = 1,598). A final subselection was made manually to represent different regulatory fields (e.g., food additives, industrial chemicals, plant protection products), different chemical families, and different MOAs.

Results: A final selection of eight compounds was examined in the zebrafish embryo model, and craniofacial malformations were observed in embryos exposed to each of the compounds, thus confirming the developmental toxicity as predicted by the in silico methods. When exposed to a mixture of the eight compounds, each at its NOAEL, substantial craniofacial malformations were observed; according to a dose-response analysis, even embryos exposed to a 7-fold dilution of this mixture still exhibited a slight abnormal phenotype. The cumulative effect of the compounds in the mixture was in accordance with dose addition (added doses of the individual compounds after adjustment for relative potencies), despite different MOAs of the compounds involved.

Discussion: This case study of a complex mixture inducing craniofacial malformations in zebrafish embryos shows that dose addition can adequately predicted the cumulative effect of a mixture of multiple substances at low doses, irrespective of the (expected) MOA. The applied workflow may be useful as an approach for CRA in general. https://doi.org/10.1289/EHP9888

Introduction

Current human risk assessment of toxicological effects of exposure to chemical compounds is mainly focused on single compounds. However, both humans and other organisms are actually exposed to a multitude of different compounds from environmental sources at the same time, and often even to chemicals inducing a same effect. The importance of combined exposures and ensuing cumulative risk assessment (CRA) has long been recognized, leading to calls for action from regulatory bodies. As a response, protocols were drafted to address CRA. The European Food Safety Authority (EFSA) also developed such an initiative, and proposed to apply dose addition, which is a basic assumption to model combined effects as a default to achieve CRA in humans, regardless of (dis)similarity in mode of action (MOA).4'5'6'7 This was justified because effects greater than those predicted by dose addition through, for example, synergistic action, have only been reported at low rates (5-0.8%) in models for human toxicology.8'9 In the opposite case, that is, actual combination effects lower than those predicted by dose addition, dose addition would provide a conservative approach.

In the approach proposed for CRA, the EFSA recommendation is to only combine doses of compounds that induce the same toxicological effect, and such compounds are hence clustered into cumulative assessment groups (CAGs).10 CAGs are hierarchically organized, considering the target organ or target tissue, the specific phenotype of the effect, the MOA, and the mechanism of action. Clustering of compounds into CAGs reduces the number of compounds to be considered for CRA of a given effect, making experimental analysis more manageable. It also enables researchers to distinguish combinations of compounds with similar and dissimilar MOA, which is a key issue in CRA.11 In its Opinion on dissimilarity of MOA, the EFSA Panel on Plant Protection Products and their Residues recognized the lack of knowledge on mode- and mechanism of action of most compounds, and at the same time concluded that there are no criteria on how to define dissimilarity. While acknowledging such applicability drawbacks, we here define dissimilar MOA as the probability that one or more major steps in the toxicological pathway leading to the induction of the specific pathological phenotype under study, including initial molecular interactions of a compound with biological targets, differ among compounds.

To evaluate and substantiate the EFSA approach for CRA, and, in particular, assess the accuracy of the dose addition assumption, several case studies were conducted in the context of the EU-Horizon2020 project EuroMix (https://www.euromixproject.eu/), mainly focusing on pesticides in food. The EuroMix project aimed to contribute to a CRA approach through developing pragmatic tools using in silico and in vitro models, combined in an openly accessible web-based toolbox.7' One of the example CAGs considered in EuroMix was Developmental toxicity, specifically focusing on craniofacial malformations.13 A test model for this purpose is the zebrafish embryo (ZFE), in which a specific measure in the head skeleton, that is, the angle formed by the Meckel's and the palatoquadrate cartilages (M-PQ angle; Figure 1), is used as a proxy for craniofacial malformations in mammals.14 The ZFE was included in a battery of complementary models informing the adverse outcome pathway (AOP) that describes the disruption of retinoic acid metabolism leading to developmental craniofacial defects, including cleft palate, in mammals.

ZFE have been used for >60 y in toxicology. A particular advantage of the model is the external embryonal development, which, combined with their transparency, makes zebrafish embryos accessible for morphological study of teratological effects. For the purpose of this study, visualization of the skeleton, and, in particular, accurate quantification of the M-PQ angle, was enhanced through Alcian blue staining (Figure 1). 7' Zebrafish are relatively easy to house and have a large offspring, making them suitable for high-throughput experiments. Zebrafish share a high degree of genetic similarity to humans given that >70% of human genes have one or more orthologous zebrafish genes, and the similarity in developmental mechanisms underlying skeletal morphogenesis compared with humans renders the ZFE a particularly relevant model for this study. Another practical benefit of the model is that 0-5 d post fertilization ZFE are not considered as laboratory animals. ' In the past, several binary mixtures of reference compounds inducing craniofacial malformations, including cleft palate in mammals, were tested in this model. Combinations of either similar or dissimilar MOA were shown to produce their effects according to the principle of dose addition, regardless of MOA. Similarly, dose addition was observed with binary mixtures of the same compounds in an additional model, namely, differentiating embryonal stem cells. 5

Although dose addition may thus explain the combined effect in binary mixtures of compounds at effect doses, a next challenge is to confirm that dose addition also holds for complex mixtures

at low doses, that is, combinations of (much) more than two compounds, at doses that do not induce an observable effect in toxico-logical studies, because such mixtures provide a more realistic representation of human exposure. An accepted statistical measure for a such a low dose is the highest dose without an observable effect in a dose range in a given toxicological model, that is, the NOAEL. Operationalization of this challenge into a research question central to the study described here would then be as follows: Does dose addition also hold for combined exposure to substances at or below their individual NOAELs? Answering this question in the context of real-life human exposure requires a well-defined approach for selection of relevant compounds, and application of such an approach was therefore a secondary aim of the study. Given that real-life mixtures will usually consist of unrelated compounds, with different MO As, this condition should be taken in account in the study design.

Selection Procedure

Our proposed approach for CRA consists of the steps shown in Table 1. In the present study, we performed the first five steps, as a proof of principle for the approach; the final step, the actual risk assessment was not conducted but is explained in the "Discussion" section. First, we decided to use the above-explained craniofacial malformations as a subcluster in the CAG Developmental Toxicity, using ZFE as the experimental model. At the same time, we used the EuroMix Chemical Inventory (Excel Table SI) as the source for compounds relevant to human exposure. This inventory consists of pesticides identified by the EFSA (e.g., Nielsen et al. ) and was supplemented with nonpesticide substances with food relevance. The resulting list contains 1,598 nonnatural substances, representing a variety of chemical classes and a wide range of applications, such as food additives, plant protection products (PPPs), PPP-metabolites as residue on food, food contact materials, industrial chemicals, environmental pollutants, and pharmaceuticals. All these substances are currently allowed to be present in food at (individual) acceptable levels by EU regulations. From this inventory, we selected potential test compounds relevant to the testcase CAG by applying an in silico-based protocol with the aim of predicting potential membership of a substance to the CAG for Developmental Toxicity, subcluster Craniofacial Malformations. This protocol was based on a) calculations of chemical structure similarity to reference compounds representing different MOAs leading to craniofacial malformations in zebrafish embryos,27 b) modeling of binding to selected receptors (molecular docking), and c) quantitative structure-activity relationship (QSAR) analysis. From this selection, an experimentally manageable subset of eight substances was chosen, representing different MOAs and chemical classes. The ability of the selected compounds to induce craniofacial malformations was then confirmed through individual testing in the ZFE. Furthermore, the prediction of MOAs among compounds was underpinned using a previously defined specific gene expression marker set detecting MOAs related to craniofacial malformations.27 Finally, these eight compounds were tested as a complex mixture of their highest concentrations without an observable effect in the model, that is, their NOAELs, as derived from the single compound experiments. In addition, various dilutions of that NOAEL mixture were examined in the mixture experiment. The present study may help in developing an approach for CRA, by combining the selection of CAG compounds using in silico methods and estimating the response from human mixture exposures by dose addition based on estimated RPFs using the ZFE test.

Materials and Methods

In Silico Compound Selection

The EuroMix chemical inventory (Excel Table SI), containing 1,598 substances from several regulatory silos (see the "Introduction" section), was used as starting point for the screening for compounds with a high probability to be part of the CAG Developmental Toxicity, subcluster Craniofacial Malformations. Selection of compounds predicted to be developmental toxicants, and distinguishing these from compounds not likely to be CAG members, was done using the in silico tools available for developmental toxicity. The selection started with calculation of chemical similarity based on Tanimoto coefficients using the publicly available Molecular Access System (MACCS) fingerprints, as compared with reference toxicants, representing different MOAs as previously defined. ' ' ' For this study, the reference compounds comprised triazoles, retinoids, dioxins, aryl hydrocarbon receptor (AhR) activators, histone deacetylase inhibitors, and dithiocarbamates, as well as 2,4-dinitrophenol (2,4-DNP)-like, and boric acid-like and fenpropimorph-like substances (Table 2 shows the represented MOAs; Figure 2, the chemical structures). In view of high similarities among AhR-activating substances [such as poly-chlorinated dioxins and polychlorinated biphenyls (PCBs)], the substance selection for this MOA was refined by using computer-generated receptor binding information (molecular docking calculated receptor binding energies) for human AhR.27 Molecular docking calculations were also considered for human retinoic acid related receptors (RARs and RXRs), and computationally determined binding energies to the RARy were specifically used as a supportive asset for activity in the retinoid pathway. AhR and RARy structures were retrieved from the Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB) [PDB entry: 5NJ8]29 and [PDB entry: 1FD0].30 In silico molecular docking was carried out with the Molecular Operating Environment Dock Program (version 2019). "Triangle Matcher" was selected as the ligand placement methodology. Results of the chemical similarity and molecular docking methods were complemented with a QSAR analysis, applying four relevant QSAR models, that is, Computer Assisted Evaluation of industrial chemical Substances According to Regulations-Developmental Toxicity model (version 2.1.7) (CAESAR) and Proctor & Gamble Developmental/Reproductive Toxicity library (version 1.1.0) (PG) [both models are implemented in the stand-alone software VEGA QSAR (version 1.1.5 48) ' ], Derek Nexus, and Toxicity Estimation Software Tool (T.E.S.T.; version 5.1.1).36 These models indicate the potential of a substance to be developmentally toxic in general. The chemical structure representation of all 1,598 entries in the EuroMix inventory is given in Excel Table S1 in the form of a Simplified Molecular

Input Line Entry System (SMILES) code. This SMILES code can be used directly as input to the respective QS AR software (see html links above). All results data (QSAR predictions, Molecular Docking calculated receptor binding energies and chemical similarity coefficients) are given in Excel Table S1. Candidate substances with a chemical similarity score of at least 50% similarity to one of the reference substances or molecular docking energies of < - 7Kcal/mol to the human AhR were assessed for the presence of at least two of four positive QSAR predictions. This criterion of at least two of four QSARs predicting developmental toxicity was based on the minimum number of positive QSARs observed for the eight MOA reference substances (Table 2). The resulting selection of compounds was the basis for a further manual subselection of eight compounds, starting with the substances most similar to the MOA reference substances, taking into account practical considerations (e.g., commercial availability, solubility) and aiming at a balanced representation of MOAs, as well chemical classes and uses (e.g., not selecting only PPPs in our final selection), for validating the applicability of dose addition in the ZFE experiments. The in silico predictions were validated by referencing QSAR predictions against existing data, further indicated as "available historical data." Thus, from the compounds for which all 4 QSARS predicted no activity, zebrafish embryo toxicity tests (ZFETs; see the section "Experimental Analysis") of three substances, namely, penicillin-G (unpublished result), d-mannitol,37 and thiacloprid,38 were found and could be compared with the QSAR data. Furthermore, unpublished ZFET data of 7 PCB congeners [unpublished results from the EU-FP6-project Assessing the Toxicity and Hazard of Non-dioxin-like PCBs present in food (ATHON); https://cordis.europa.eu/ project/id/22923; accessed 31 January 2022] could be compared with QSAR results for these compounds.

Experimental Analysis

The final selection of eight compounds was tested in zebrafish embryos in subsequent experiments (Table 3) to a) examine general toxicity in the ZFET; b) confirm induction of specific effects on craniofacial formation through M-PQ angle measurement after Alcian blue staining and derive the (relative) potency of each compound; c) examine dose addition, and analyze combined effects in a mixture at NOAEL of each compound and lower; and d) analyze the expression of specific marker genes at the critical effect dose at the 20% effect level (CED2o). All experimental procedures are explained below. All compounds were purchased from Sigma Aldrich except for sodium 5-nitroguaiacolate, which was obtained from Santa Cruz Biotechnology, Inc.

Maintenance ofZebrafish

Zebrafish (Danio Rerio) adults of the AB-strain (obtained from the European Zebrafish Resource Center) were maintained in the Dutch National Institute of Public Health and Environment (RIVM) zebrafish facility. The fish where kept in 8-L Tecniplast ZebTec tanks with gravel images at the bottom as cage enrichment, under a 14:10 h light:dark cycle. The automatic ZebTec flow-through system kept the temperature at 27°C ± 1°C, the conductivity at 500 ± 100 us, and the pH at 7.5 ± 0.5. The fish were fed three times a day, twice with SDS400 (Special Diet Services; Tecnilab-BMI BV) and once with fresh Artemia (Special Diet Services; Tecnilab-BMI BV).

For breeding, three males and three females were placed into one 1.7-L sloped breeding tank (Tecniplast) in the afternoon before spawning, which was triggered by dimming of the light the next morning. Collected eggs were washed with experimental medium [Dutch standard water (DSW): demineralized water with 200 mg/L calcium chloride dihydrate, 20 mg/L potassium bicarbonate, 180 mg/L magnesium sulfate heptahydrate, and 100 mg/L sodium bicarbonate], and eggs at 4 to 64 cell stages were then selected for quality (no coagulation, asymmetry, or formation of vesicles or damaged membranes) under a stereo light microscope.

Embryo Exposure and Morphological Analysis

The toxic potency of the test compounds was first established in a ZFET19 in a dose range with at least eight concentrations per compound (Table 4),N=10 embryos per concentration, 1 embryo per well, with 1.5 mL of test medium in a 24-well plate. The plates were placed in an incubator at 27°C under a 14:10 h light:dark cycle. All compounds were dissolved with 0.1% dimethyl sulfoxide (DMSO; Sigma Aldrich) as a solvent. Exposure started immediately after quality assessment of the eggs [before 2 h post fertilization (hpf)], and was continued until 72 hpf. At that time, developmental effects were morphologically scored in a standardized way, with tail detachment, somite formation, eye development, movement of the embryo, heartbeat, blood circulation, embryo pigmentation, pectoral fin, protruding mouth and hatching as morphological end points for embryo development. A balanced assessment was made with lethality and morphological effects to determine the highest applicable concentration in further testing, although, generally, lethality was used as the leading factor, with information from morphological scoring as supporting information for developmental delays. Scoring of teratological effects, which is part of the original scoring system, was not required to derive exposure ranges for further experiments in this study, and is therefore not reported. ZFET in available historical data (penicillin-G, d-mannitol, PCBs) was performed exactly as described above, with the exception of using zebrafish, which were bred from commercially obtained wild-type zebrafish (Ruinemans B V). 7 ZFET with thiacloprid was performed in a parallel study within the EuroMix project, and therefore applying the same method as above. These historical data were included for comparison with in silico predictions of developmental toxicity, and therefore reporting both morphological scores (developmental delays) and teratological scores was considered relevant.

The dose range for the M-PQ analysis was based on the knowledge about concentrations inducing lethality and morphological developmental effects from the first series of experiments (ZFET results), targeting the estimated highest sublethal concentration during the ZFET with five to six subsequent half-logarithmic dilutions and occasionally an intermediate concentration (Table 4); a blank solvent control was added, as well as one higher concentration to account for potential interexperimental variation. Ambient conditions for M-PQ analysis were as for the ZFET. Each dose group contained TV = 10 embryos (target, lower numbers occurred as a result of lethality during exposure or due to loss during processing; see Excel Table S2 for exact numbers), exposure started immediately after quality assessment of the eggs (usually within 2 hpf), and exposure duration was continued until 120 hpf. These dose-range analyses yielded no observed adverse effect levels (NOAELs) per compound, which were used as the basis to compose the mixture, that is, the eight compounds were combined at the ratio of their NOAELs for skeletal malformations (Table 4). The mixture itself was designed and tested as 100% NOAEL, together with additional mixtures in the range of 3-1,000% of that mixture, all at the same ratio. The 100% NOAEL concentration of each single compound conditions was replicated together with the mixture experiment to confirm the absence of observable effect at that concentration (and to exclude such input into the mixture due to, for example, slight dilution differences between experiments 2 and 3).

For optimal assessment of head skeleton malformation, visualization of cartilage structures was enhanced by Alcian blue staining at 120 hpf. In short, this protocol included manual dechoriona-tion of unhatched embryos, euthanasia through rapid cooling on ice, fixation in freshly prepared 4% paraformaldehyde in phosphate-buffered saline (PBS) overnight at room temperature, removal of the fixative, and three consecutive washes in PBS. Pigmentation was removed by a bleaching solution (1 part 2 M potassium hydroxide with 2 parts 3% hydrogen peroxide) until visual transparency of the eye was achieved, usually after ~ 1 h. After three subsequent washes with PBS, the actual staining was through overnight incubation at room temperature in an Alcian blue solution [a filtered solution of 100 mg of Alcian blue 8GX in 70 mL of 100% ethanol, 25 mL of Milli-Q water, and 5 mL of 37% hydrochloric acid (HC1), mixed for 3 h at 50°C]. This was followed by removal of the staining solution and three 20-min washes with destaining solution (210 mL of 100% ethanol, 75 mL of Milli-Q water, and 25 mL of 37% HC1). Properly destained embryos were stored in 100% glycerol. Stained embryos were visualized with an Olympus C5050 ZOOM digital camera mounted on a Leica Labovert FS microscope camera, in a fixed position in a glass capillary that was submerged in glycerol to prevent light distortion. The M-PQ angles (Figure 1) were measured with the angle tool in Photoshop CC 2018 (Adobe Systems Inc.).

Gene Expression Analysis, Real-Time Quantitative Polymerase Chain Reaction

Gene expression analysis was performed as described previously. 7 In short, replicate pools of ZFE (N = 6-10) per compound, each pool consisting of 12-15 embryos, were exposed at a concentration close to the compounds CED20 for M-PQ angle (see the "Results" section and Table 4). Exposure was from ±2to72hpf, after which ZFE were euthanized, snap frozen, and stored at -80°C. Further processing of ZFE, RNA isolation, RNA quality check, cDNA synthesis from RNA isolates, and subsequent realtime quantitative polymerase chain reaction (RT-PCR) were performed according to standard procedures,27 using Taqman gene expression assays (Thermo Fisher Scientific). A limited gene expression marker set was derived as most informative for the action of the concerning compound families from a much wider array in a previous study. 7 Table 5 lists the analyzed genes, including the targeted marker genes and the genes used for reference, namely, glyceraldehyde 3-phosphate dehydrogenase (gapdh), hy-poxanthine phosphoribosyltransferase 1 (hprtl), and actin beta 1 (actbl). The RT-qPCR results are presented as relative expression levels, calculated as 2 - AACT47

Statistical Analysis

The observed ZFET scores and M-PQ angles were analyzed by dose-response analysis using the PROAST software (version 70.3; Slob 2002;62 https://www .rivm.nl/en/proast) and appeared to be best described with the four-parameter exponential model of a set of nested models (in PROAST: E5-CED): y = aX (c(1_exp(~-I/fc) )), with parameters a, b, c, and d describing the response at dose 0 (background value), the potency of the chemical, maximum fold change in response compared with background response (upper plateau), and steepness of the curve (on a log-dose scale), respectively. For M-PQ analysis, we used the three-parameter exponential model, by fixing the fourth parameter, c (the maximum fold change), to the maximum value possible: that is, a change by a factor 5 because five times the background M-PQ angle of ~ 35° is ~ 180°. The reason for omitting parameter c from the model (as a free parameter) was that the maximum M-PQ angle is not well defined for some compounds (Figure 1), making the assumption that it would be the same among compounds doubtful. Thus, the remaining parameters to be estimated are a (mean response at dose 0; background), b) (CED, potency), and d (a steepness parameter). PROAST allows for analyzing the combined dose-response data from multiple compounds in a single analysis, that is, the model is fitted to all compounds in one fit, by including compound as a covariate in the dose-response model. In this analysis, parameter d is assumed to be equal for all compounds. After fitting the model, this assumption was visually checked and appeared reasonable, in line with earlier experience with toxicological data in general (e.g., Slob and Setzer). In fitting the model, parameter a and CED or RPF were made dependent on the compound as a covariate.

In all single compound dose-response experiments, the significance of a craniofacial effect was evaluated by a trend test, in particular, by statistically testing the fitted model against the no-response model (y = a) using the log-likelihood ratio test or the Akaike information criterion (AIC) in PROAST. CEDs, together with their 90% confidence intervals (expressed as the lower and upper bounds: CEDL, CEDU), were calculated at a predefined critical effect size (CES) of a 20% increase in MPQ angle, a value which was also used in previous work with this model. 7 The 20% value is assumed to represent an effect well beyond the background variation while accommodating the maximally achieved effect with all compounds at the same time. Given that the CEDs for the various compounds represent equipotent doses, the ratio of the CED of a given compound compared with the reference compound is the relative potency factor (RPF). Note that each compound can serve as a reference compound because RPFs will be proportional with either choice. Also note that the RPF does not depend on the value of CES, given that the dose-response curves are parallel on a log-dose scale (Figure 3). The RPF can be used to express the dose of a given compound as the equivalent dose of the reference compound. PROAST includes the option for RPF estimation based on dose-response modeling and can provide confidence intervals for the RPFs (model 46).

Whereas the CED20 is a more precise estimate of the potency of a compound than the NOAEL, the latter was used as an accepted measure for the purpose of designing a mixture at low doses. NOAELs were calculated based on the data obtained from the M-PQ data of each single compound (experiment 2), using Graphpad Prism (version 9.1.0). First, potential differences between the means of the dose groups were analyzed with a Kruskal-Wallis test. Next, a post hoc Dunn's multiple comparisons test was performed to detect significant differences of the dose groups as compared with the nonexposed control (0.1% DMSO). The highest concentration with no statistical significant difference compared with the blank solvent control (i.e., 0.1% DMSO) of each compound was considered as the NOAEL for skeletal malformation. A 100-mgp < 0.05 was considered statistically significant.

To evaluate whether the cumulative effect of the mixture can be explained through dose addition even at low doses (i.e., at or below the NOAELs of the individual compounds), a range of dilutions of the mixture consisting of NOAELs of the compounds was applied to the embryos. In case of dose addition, the mixture dose response can be predicted by the sum of the RPF-adjusted doses of the single compounds. Evaluation of dose addition is through visual comparison of the mixture responses with the fitted curve, which predicts dose addition. Dose addition is likely if the pattern of mixture response and the confidence intervals per mixture dose are in agreement with the fitted curve. Systematic deviation from the fitted curve may, apart from experimental errors, indicate a higher (left-shift) or lower (right-shift) of the mixture, which may result from interaction of the compounds through, for example, synergy or antagonism, respectively. In such cases, this visual inspection will directly provide understanding of the magnitude of the deviation and, thereby, its impact on risk assessment as compared with the situation based on dose addition.

For each compound, the estimated dose (CED20) associated with a 20% change in M-PQ angle was applied to the zebrafish embryos in separate experiments (experiment 4) to measure gene expression at that dose. The statistical significance of differential gene expression between compound-exposed and control samples was tested with a f-test [p < 0.05; R, version 3.6.0 (R Development Core Team)] on log-transformed expression values per gene.

Results

In Silico Compound Selection

Application of the specific selection criteria for chemical similarity, complemented with molecular docking modeling binding energies of the substances to AhR and RAR/RXR, and criteria for available developmental toxicity QSAR model predictions to determine likely developmental toxicants (Excel Table S3) predicted 79 candidates from the complete inventory of 1,598 substances to be very likely developmentally toxic, giving craniofacial malformations (Excel Table S4). Furthermore, of the 1,389 substances for which the QSAR models could generate a prediction, only 33 substances were predicted by all four QSAR models to be developmental toxicants (Excel Table S6), whereas 229 were considered to be very unlikely developmentally toxic, given that none of the four QSAR models gave a positive result (Excel Table S5). The remaining substances for which QSAR predictions could be generated (1,127) had partially positive results from the battery of four QSARs (1-3 models positive). From the 305 substances with 50% or more structural similarity to the reference compounds (Table 2), or AhR binding energy of < - 7 kcal/mol, a total of 79 substances were also predicted by two or more of the four QSAR models to be likely developmental toxicants. Two of these 79 mixture candidate compounds (tebuconazole and triadimenol, which is the metabolite of the reference substance triadimefon) overlapped with the EFS A lists of developmental toxicants potentially related to craniofacial malformations (see Tables 25.26 and 25.25 in

Nielsen et al.13). These were excluded for the final subselection because they had already been experimentally confirmed as developmental toxicants inducing craniofacial malformations.

For the purpose of validating the applicability of dose addition in a complex mixture consisting of doses at their NOAELs, an experimentally manageable selection of 8 compounds was made from the list of 79 candidate toxicants potentially leading to craniofacial malformations in zebrafish embryos (Table 6, and in Excel Table S4 indicated in green), aiming at a balanced representation of chemical families and MOAs. The candidate compounds with the highest similarity to the reference compounds (given in order of decreasing similarity in Excel Table S4) were initially selected. One first choice compound, methenammonium chloride (MOA: AhR binding), was not commercially available and was, therefore, replaced by two other AhR binders; benzalkonium chloride, with a comparably high calculated AhR binding while accepting its lower score of positive QSARs, and propyl gallate with much lower predicted AhR binding but more positive predictions from the QSARs. Two categories were left out: boric acid-like compounds, because no compounds with sufficient similarity were retrieved from the inventory, and RA-like compounds, because the selected compound with highest similarity to RA and a high predicted RARy binding energy, beta-apo-8'-carotenal (food additive E160e), could not be sufficiently dissolved. Other retinoid candidates, such as capsanthin and capsarubin (E160c) and Annatto (E160b), were also expected to suffer from extremely low water solubility. Other candidates with sufficient similarity to retinoic acid were calculated to have low or nonbinding receptor docking energies (Excel Table S1).

To further evaluate the predictivity of the in silico selection procedure, the availability of existing ZFETs in our historical

ZFET collection (available historical data) was checked for any of the 229 substances that according to the in silico criteria were considered very unlikely to induce developmental effects (Excel Table S5). Three substances, penicillin-G, D-mannitol, and thiacloprid could thus be retrieved, all three confirming the absence of developmental effects up to very high exposure doses (Table 7; Excel Figure SI). In addition, ZFET developmental toxicity screening results were available for seven PCB-congeners, each at a single fixed concentration of 1 uM (Table 8; Excel Figure SI). Of these, 1 |±M PCB180 induced no observable effect in the ZFET and did not get positive predictions from any of the four QSARs. Further, 1 |±M PCB126, which was used as the reference substance for AhR-activation related to developmental toxicity (cf. Table 2), showed the highest score in this PCB collection for teratological effects in the ZFET while also resulting in the highest number of QSARs predicting developmental toxicity. The other five PCB congeners showed intermediate teratological scores in the ZFET and an intermediate number of positive QSARs. Altogether, ZFET and QS AR results were well related among these PCB s.

Experimental Analysis of Single Compounds: Experiments 1 and 2

According to the experimental scheme (Table 3), analysis started with the determination of effects of each of the eight compounds in a ZFET, with the purpose of establishing a dose range in the ensuing experiments. These ZFET results indicated large potency differences among the selected compounds in inducing lethality and morphological developmental effects, whereas even no lethality was observed with paclobutrazol (PBZ) at maximal solubility (detailed results are shown in Excel Figure S2; see also the concluded highest applicable concentrations in Table 9, second column).

In the second series of experiments (Table 3, experiment 2), dose-response data were generated for head skeleton malformations with each compound tested individually (Table 9). An example of increased M-PQ angle as a measure for head skeleton malformations is shown in Figure 1, after exposure to 2-ethylhexa-noic acid (EHA). Figure 1 (right) also illustrates that high concentrations of some compounds, in this case again EHA, produced severe agenesis of the head skeleton, leading to failure to measure the M-PQ angle. For all compounds, there was a significant trend between dose and increasing M-PQ angle (combined plot in Figure 3, per compound plots in Excel Figure S3). The first purpose of this second series of single compound experiments was to find the NOAELs to be used in the next experiment (Table 3, experiment 3), that is, as the starting point for composing the NOAEL mixtures. The NOAELs thus identified in these experiments are listed in Table 9 (for full analysis, see Excel Figure S4), and from these data, the compounds sodium-dimethyldithiocarbamate (SDC) and sodium-5-nitroguaiacolate (S5N) were estimated to be the most and least potent compound, respectively. The second purpose of these experiments was to find equipotent doses for the marker gene experiment. The CED20S, which were calculated for that purpose, showed a high correlation with the NOAEL (Table 9). The confidence intervals of these CD20S were relatively small (Table 9; Excel Figure S3), with CEDU/CEDL ratios < 1.3 (Table 9).

Next, the dose-response data of all eight compounds were combined and analyzed in a single dose-response analysis, with compound as a covariate, and 2-EH as the reference compound for expressing RPFs of the other compounds. The resulting fit of the model (Figure 3A; individual observations per embryo are shown in the separate curves for each compound in Excel Figure S3) revealed three clusters of compounds, with SDC as the most potent compound (RPF ~ 1,550), benzalkonium chloride (BAC), fenpro-pidin (FEN), PBZ with intermediate potencies (RPFs ~ 15), and 2-ethyl-1-hexanol (2EH), propyl gallate (PG), EHA, and S5N as the least potent compounds (RPFs close to 1; note that the RPF of the reference 2-EH is by definition equal to 1). The resulting confidence intervals for the RPFs (see the legend to Figure 3, Table 9) were narrow (less than a factor of 1.25). Figure 3B shows the same model fit, but now against dose in terms of the reference compound 2-EH, that is, by adjusting each compounds doses by the relevant RPF of that compound.

Mixture Exposure

In the third experiment, a range of mixtures of the eight compounds was composed based on their NOAELs (shown in Tables 4 and 9). However, the mixtures consisting of 300% and 1,000% of the NOAELs were lethal, leaving 100% as the highest usable dose for analysis.

Although the effects of the eight compounds at the NOAEL remained undetected in the individual compound analyses, embryos exposed to the combination of compounds at their NOAELs had substantially greater M-PQ angles than those in the control and single compound groups (Figure 4). Figure 5 shows the observed M-PQ angles against the dilution series of the NOAEL mixture. Dose-response analysis resulted in observable effects down to 30% NOAEL dilution (= 1.5 on a log scale), and a CED05 calculation produced a value of ~ 14% NOAEL mixture, that is, a 7-fold dilution of the 100% NOAEL mixture was estimated to result in a 5% increase of the M-PQ angle.

Next, the combined data set of the single compound experiments (as analyzed in Figure 3) was combined with the data from the mixtures, and reanalyzed in the same way as in Figure 3B, but now with the response from the mixtures included as well. This analysis (Figure 6) shows that the mixtures overlapped with the individual compound responses. Note that in this model fit, the mixtures contributed to the estimation of the RPFs; therefore, the estimated RPFs were slightly different from those in Figure 3.

Gene Expression Analysis

In this additional series of experiments, we aimed to confirm differential toxicological action of the eight compounds on the level of gene expression. Application of the defined gene expression marker set (Table 5) showed that the eight selected compounds had partly overlapping and partly nonoverlapping profiles (Figure 7, Table 10; see Excel Table S7 for underlying data). Particularly distinctive profiles could be observed for 2EH (relative strong down-regulation of cyp26al and igfbpl), EHA (relative strong down-regulation of cyp26al), SDC (relative strong down-regulation of cyp26bl, dlx5a, and loxl3b), and PBZ (relative strong up-regulation of cyp26al, cyp26bl), whereas BAC, FEN, S5N, and, to a lesser extent, PG showed much overlap (relative strong up-regulation in most genes, although with differential down-regulation of cyp26al and ntla) and thereby separated from the other compounds.

Discussion

At this time, there is no internationally harmonized approach for the CRA of chemical substances relevant in view of human exposure, such as food-contaminating pesticide residues or environmental pollutants.7 To fill this gap, the EFS A has proposed to apply dose addition as a default for CRA, mainly focusing on plant protection products present as food contaminants. In the EFSA strategy, CAGs are defined by toxicological phenotypes, where developmental toxicity is one such a CAG, with craniofacial malformations, including cleft palate, as a subcluster. For this subcluster, dose addition was previously tested and confirmed as a valid prediction of combination effects of two reference compounds with both similar and dissimilar MO As, in zebrafish embryos, ' and in the embryonal stem cell test. The ensuing exercise was to apply the principle of dose addition to real-world conditions, where human exposure is characterized by combined exposure to multiple chemicals, all (or most) below a NOAEL identified in toxicological models. Our results clearly confirmed the hypothesis that dose addition holds for low doses just as well, including doses below the NOAEL.

In practice, it will not be possible, nor ethically justified, to apply in vivo testing to "all" compounds with potential exposure in humans for the purpose of establishing potential effects regarding the end point considered. As an alternative, we explored the application of in silico methods to select candidate compounds that are likely to have a significant contribution to health risk for the CAG of craniofacial malformations. Although in silico methods are useful as a first screening and selection of chemicals with a potential effect, they are generally not suitable to establish (relative) potencies. The latter needs to be examined experimentally, in this case M-PQ measurement in zebrafish embryos. The resulting RPFs can be used in a CRA, by adding realistic human exposures according to the dose addition principle, that is, by adjusting each single compound exposure by its RPF, resulting in an exposure expressed in terms of the reference compound.

Selection of Chemicals

For the purpose of compound selection for this mixture experiment, the applied combined in silico tools identified both likely developmental toxicants and substances that are not likely to give observable developmental toxicity effects. The subselection of 8 compounds taken from the initial selection of 79 positive compounds, as well as the additionally selected AhR-activating substance propyl gallate, were all confirmed as developmental toxicants in the ZFET, whereas some of the predicted negative compounds were in compliance with available historical ZFET data (Table 7; Excel Figure SI). This indicates that the applied in silico methods may be an appropriate tool for the purpose of compound selection into CAGs, at least for the case of the CAG of developmental toxicity. More specifically, M-PQ analysis experimentally confirmed the potential of all 8 selected compounds to induce craniofacial malformations, as predicted through in silico analysis.

Regarding the QSAR models, it should be noticed that predictions of developmental toxicity (like other specific

Environmental Health Perspectives

047003-11

130(4) April 2022 toxicological end points) are based on associations of both the chemical functionalities and the physicochemical properties and topological structure indicators associated with in vivo effects in mammals, notably in rats. Vice versa, predictions of the absence of developmental toxicity are based on the absence of chemical functionality and similarity to the reference data set of substances that are not developmentally toxic in mammalian toxicity experiments. The confirmative results of the additional comparison of negative predictions by in silico methods with effects of such substances in the ZFET, and of comparative toxicity ranking of a set of PCB congeners by in silico analysis and ZFET, are therefore in support of the applicability of in silico methods in ZFE. However, formal validation of the method would require an analysis of a much larger data set, and it should be kept in mind that the in silico procedure is not a tool for predicting the toxicity potency of potential CAG members.

The selection criteria both for the chemical similarity calculation (50% Tanimoto score) and for the number of positive QS ARs (2 of 4) were set arbitrarily, aiming at a high probability to select positive candidate substances for the mixture experiment. Structural similarity to previously identified reference compounds (Table 2) was introduced in our selection procedure to allow for a mixture selection that explicitly represents different toxicity pathways leading to craniofacial malformations. However, MO A (toxicity pathway) is only auxiliary when identifying CAG membership, and therefore the cutoff criterion of 50% structure similarity is not an absolute value that can be objectively evaluated. Such a cutoff value is highly dependent on the method of calculation, that is, which structural fingerprint is used, and on the specific similarity coefficient applied (in this case Tanimoto). Furthermore, it should be kept in mind that when structural similarity is used for other purposes, particularly for the risk assessment of real-life mixtures, it should be confirmed that the reference database contains all relevant chemical structures.

Similarly, the criterion to select substances with at least two of the four available QSAR models predicting toxicity for a specific phenotype appeared to be useful for the exercise in this study. However, the additional inclusion of propyl gallate, which was predicted to be a strong AhR binder, but scored only one positive QSAR, indicates that a minimum of two positive QSAR models is possibly too restrictive for CAG-membership identification, and a single positive QSAR may be considered sufficient for a more inclusive approach to CAG membership definition. For that purpose, a weight-of-evidence evaluation of available QSAR models, considering their specificity, is recommended.

At the time of this study, AhR, RAR, and RXR were the only available molecular docking receptor models representing a complementary in silico method for selecting molecules relevant to the CAG. Meanwhile, this collection has been further expanded to models for other relevant nuclear receptors and enzymes, such as cyp26, involved in retinoic acid metabolism.27 These could be used in future extensions of this selection strategy.

Using chemical similarity and read-across to identify CAG membership is limited to identification of substances similar to the known references. References used herein were derived from the previous testing of a wide variety of compounds in ZFE. ' which in turn were selected on the basis of an inventory of inducers of skeletal malformations and cleft palate in rodents,13 supplemented through literature search that included human data. Although major chemical families and associated toxicological actions 7 were covered in this way, the collection is probably not complete, with subsequent exclusion of potentially relevant substances. As an example, diethylstilbestrol (DES) was predicted to be developmentally toxic by all four QSAR models but not identified by similarity analysis because of a missing appropriate reference structure.

Dose Addition

Dose addition is generally accepted for mixtures with similar MOAs, but disputed for dissimilar MO As. We have previously shown, both in ZFE and in EST, that dose addition predicts observed cumulative effects in binary mixtures, irrespective of MOA.14'24'25 The implication of this result is that MOA is not an important criterion in CRA, which thereby can be considerably simplified. It also supports mixture modeling of the many compounds that have multiple or even completely unknown MOAs.

It has been argued that mixtures consisting of individual doses below the NO AEL will not result in effects, based on the argument that the sum of zero effects is again zero (e.g., Cassee et al., and reviewed by Kortenkamp et al.2 ). In this line, a later report of three joint EU Scientific Committees stated that it was unlikely that mixtures of chemicals with different MOAs, when individually present at or below their zero-effect levels, would be of health or environmental concern.49 This paradigm was challenged soon thereafter, with evidence that the "something from nothing" principle is valid in the context of endocrine disrupters. Similarly, a more recent study showed cumulative effects in male rats of mixtures of environmental contaminants with relevance for male sexual development at dilutions of their NOAELs as low as a factor of 15, on marks related to this end point. In this context, it should be kept in mind that although an effect at the NOAEL is commonly assumed to be zero, there may very well be an effect at the NOAEL that was, however, too small to be detected owing to the lack of observational or statistical power for that small effect size (e.g., Crepet at al. ). Here, we provided unambiguous experimental evidence for the fact that a NOAEL cannot not be interpreted as a zero-effect dose (Figures 4 and 5). Furthermore, we showed that the effect of a mixture of doses at or below the NOAELs can be predicted by dose addition, using information on the RPFs, which are estimated in suitable dose-response studies. In this way, the overlapping dose-response curves of mixtures and individual compounds (Figure 6) supported prediction of dose addition, that is, the mixtures showed no shift to the left (indicating synergy), nor to the right (indicating antagonism).

Marker Genes and Toxicological Action

Composition of the mixture in this study was as diverse as possible, in terms of chemical families, based on molecular docking results and a previous gene expression study.27 In that expression study, a limited set of marker genes was generated from a comprehensive array of ~ 70 marker genes, derived from literature and microarray studies in ZFE that were known to be associated with pathways relevant to craniofacial malformations in humans, rodents, and zebra-fish. Selection of key genes was achieved through mapping to toxicological pathways relevant to impaired craniofacial development and subsequent testing in ZFE exposed to a panel of reference compounds known to activate major relevant pathways. The remaining 9 genes, used in the present study, were informative and in combination discriminative among the targeted pathways. Thus, when compared with regulation by reference compounds (Table 2), it appears that triazole-like activity was induced by mainly PBZ and to a lesser extent by FEN, BAC, and S5N; RAR-like activity by 2EH; DTC-like activity, although modestly, by 2EH, EHA, BAC, FEN, S5N, and SDC; AhR activity by SDC; and valproic acid (VPA)-like activity by 2EH, EHA, BAC, FEN, S5N, and SDC. Clearly, the limitations of this analysis are that a) pathways interact, and some genes are therefore regulated in one way or another by activation of a single pathway; b) this limited set of genes is insufficient to conclude unambiguously on precise pathway activation; and c) many compounds have multiple MOAs and may thus activate multiple pathways. 7 On the other hand, the number of pathways relevant to craniofacial malformations is limited,27 and although supportive indications for activated pathways can be obtained, the purpose of testing the set of 9 marker genes is not to exactly identify activated pathways but, rather, more to show differential toxicological action among the tested compounds. Indeed, the observed differences of expression profiles showed that no two compounds induced a fully identical expression profile, and the analysis was thus indicative for the presence of dissimilar acting compounds in the mixture. This differential gene expression in our array of compounds in the mixture therefore further supports that similar MOAs are not required for dose addition.

Cumulative Risk Assessment

This study was designed to evaluate whether dose addition applies in a complex mixture of compounds, with all doses at or below their NOAELs. It shows that such low doses in combination can lead to a substantial effect of the mixture, and should therefore not be ignored in CRA. In a CRA, the information on RPFs of the selected compounds in a given C AG can be combined with specific human exposure levels, no matter how low they might be, using the dose addition concept. Dietary mixture exposures can be calculated from data from food and drinking water surveys, as was done in a case study in the EuroMix project, for pesticides with relevance to the CAG Steatosis, eventually refined through modeling of aggregate exposure from nondietary sources in addition to dietary sources. ' Another approach to estimate exposure is human bio-monitoring, through analysis of CAG-specific contaminants and their metabolites in human biomaterial, particularly urine and blood.56'57'58 Such information, combined with analysis of potential sources of the measured contaminants, can help to build and refine physiologically based pharmacokinetic models for estimating internal human exposure.59 As explained above, exposure data of multiple contaminants can be converted into a mixture dose using RPFs, expressed as a total dose of a reference compound, where all (also low) doses of individual compounds should be included in the calculation. In this way, the CRA can be further operationalized in line with usual risk assessment of a single compound, but now with the resulting risk estimate relating to the mixture exposure. Preferably, the uncertainties in the RPFs as well as in the estimated exposure levels of each compound are taken into account, in a probabilistic risk assessment (e.g., Bosgra at al. ). In addition, further assessment of human relevance of the nonhuman, reductionistic, experimental models should be implemented in the application for human risk assessment (see Veltman et al. ).

Ideally, RPFs are generated for all nonnatural substances to which humans may be exposed, and for all test models in each CAG. However, when dealing with large numbers of chemicals, limitations in test capacity will hamper derivation of (experimental, in vivo) RPFs for all compounds. Therefore, the in silico prioritization procedure is useful for selecting substances that are most likely to contribute substantially to the CRA, thereby limiting the number of compounds for which potency information will have to be generated. Another driver for selection could be (expected) exposure concentrations; however, this involves the risk that substances will be excluded from CRA based on a (relatively) low exposure while important information on potency of contributing substances is missing. Accordingly, the present study was conducted with only a subselection from 79 candidates relevant to the CAG Developmental Toxicity/Craniofacial Malformations, for proof of principle. In practice, it will be impossible to test all compounds that might cause the CAG effect, but for a comprehensive CRA, the list of eight compounds subjected to RPF testing and exposure modeling should be expanded. A most relevant subselection would then be based on highest in silico prediction scores and highest exposures, instead of variation of MOA and chemical classes, as done in this study.

Conclusions

Overall, it can be concluded that the applied in silico methods evaluating compounds for potential effects on the CAG Developmental Toxicity, subcluster Craniofacial Malformations, resulted in a successful selection of eight compounds for further experimental evaluation. Indeed, all selected compounds did induce craniofacial malformations in the zebrafish embryos, albeit with varying potency. Experimental evaluation of a mixture consisting of the NOAELs of these eight compounds resulted in a substantial effect, and effects were even observed at dilutions of the NOAEL mixture. The responses observed at the various dilutions of the NOAEL mixture could be predicted by dose addition, which supports that dose addition holds at doses at or below the NOAELs just as it does at higher doses. Because the eight compounds had distinct toxicological actions, these findings once more confirm that dose addition does not depend on the MOA of the constituent compounds. Dose addition thus can be used for CRA of compounds that evoke craniofacial malformations. The combined steps of a) CAG definition, b) in silico prediction of CAG membership and substance selection, and c) estimation of RPFs provide a useful basis for a strategy for CRA.

Acknowledgments

We acknowledge F. den Ouden for her support of the experimental work. We are grateful to M. Luijten for her critical and constructive review of the manuscript.

This work was undertaken within the H2020-project EuroMix (https://www.euromixproject.eu/), funded by the European Commission (grant agreement 633172 to J.D.v.K.), and the Dutch Ministry of Health, Welfare and Sports (project 5.1.2: Knowledge base and policy advise on CMRS substances). Zebrafish embryo results with PCBs were generated for the EU-FP6-project Assessing the Toxicity and Hazard of Non-dioxin-like PCBs present in food (ATHON; grant agreement 22923).

References

References

1. Kienzler A, Bopp SK, van der Linden S, Berggren E, Worth A. 2016. Regulatory assessment of chemical mixtures: requirements, current approaches and future perspectives. Regul Toxicol Pharmacol 80:321-334, PMID: 27211294, https://doi.Org/10.1016/j.yrtph.2016.05.020.

2. European Commission. 2012. Communication from the Commission to the Council: the combination effects of chemicals. Chemical mixtures. 31 May 2012. COM(2012) 252 final. Brussels, Belgium: European Commission. https://www. eumonitor.eu/9353000/1/j4nvke1fm2yd1u0j9wik7m1c3gyxp/vkcwedxwhlas/v=s7z/ f=/com(2012j252_en.pdf [accessed 23 November 2021].

3. Bosgra S, van Eijkeren JC, Slob W. 2009. Dose addition and the isobole method as approaches for predicting the cumulative effect of non-interacting chemicals: a critical evaluation. Crit Rev Toxicol 39(5):418-426, PMID: 19514914, https://doi.org/10.1080/10408440902787592.

4. Cassee FR, Groten JP, van Bladeren PJ, Feron VJ. 1998. Toxicological evaluation and risk assessment of chemical mixtures. Crit Rev Toxicol 28(1 ):73 101, PMID: 9493762, https://doi.org/10.1080/10408449891344164.

5. EFSA (European Food Safety Authority). 2013. International frameworks dealing with human risk assessment of combined exposure to multiple chemicals. EFSA J 11(7):3313, https://doi.Org/10.2903/j.efsa.2013.3313.

6. Gao Y, Feng J, Kang L, Xu X, Zhu L 2018. Concentration addition and independent action model: which is better in predicting the toxicity for metal mixtures on zebrafish larvae. Sci Total Environ 610-611:442-450, PMID: 28820978, https://doi.Org/10.1016/j.scitotenv.2017.08.058.

7. Rotter S, Beronius A, Boobis AR, Hanberg A, van Klaveren J, Luijten M, et al. 2018. Overview on legislation and scientific approaches for risk assessment of combined exposure to multiple chemicals: the potential EuroMix contribution. Crit Rev Toxicol 48(9)796-814, PMID: 30632445, https://doi.org/10.1080/10408444.2018.1541964.

8. Cedergreen N. 2014. Quantifying synergy: a systematic review of mixture toxicity studies within environmental toxicology. PLoS One 9(5):e96580, PMID: 24794244, https://doi.org/10.1371/journal.pone.0096580.

9. Martin 0, Scholze M, ErmlerS, McPhie J, Bopp SK, Kienzler A, et al. 2021. Ten years of research on synergisms and antagonisms in chemical mixtures: a systematic review and quantitative reappraisal of mixture studies. Environ Int 146:106206, PMID: 33120228, https://doi.Org/10.1016/j.envint.2020.106206.

10. EFSA. 2013. Scientific opinion on the identification of pesticides to be included in cumulative assessment groups on the basis of their toxicological profile. EFSA J 11(7)3293, https://doi.Org/10.2903/j.efsa.2013.3293.

11. EFSA. 2013. Scientific opinion on the relevance of dissimilar mode of action and its appropriate application for cumulative risk assessment of pesticides residues in food. EFSA J 11(12)3472, https://doi.Org/10.2903/j.efsa.2013.3472.

12. van der Voet H, Kruisselbrink JW, de Boer WJ, van Lenthe MS, van den Heuvel JJBH, Crepet A, et al. 2020. The MCRA toolbox of models and data to support chemical mixture risk assessment. Food Chem Toxicol 138:111185, PMID: 32058012, https://doi.Org/10.1016/j.fct.2020.111185.

13. Nielsen E, Norhede P, Boberg J, Krag Isling L, Kroghsbo S, Hadrup N, et al. 2012. Identification of cumulative assessment groups of pesticides. EFSA support publ 9(4):269E, https://doi.org/10.2903/sp.efsa.2012.EN-269.

14. Staal YCM, Meijer J, van der Kris RJC, de Bruijn AC, Boersma AY, Gremmer ER, et al. 2018. Head skeleton malformations in zebrafish {Danio rerio) to assess adverse effects of mixtures of compounds. Arch Toxicol 92(12)3549-3564, PMID: 30288550, https://doi.org/10.1007/s00204-018-2320-y.

15. Menegola E, Veltman CHJ, Battistoni M, Di Renzo F, Moretto A, Metruccio F, et al. 2021. An adverse outcome pathway on the disruption of retinoic acid metabolism leading to developmental craniofacial defects. Toxicology 458:152843, PMID: 34186166, https://doi.Org/10.1016/j.tox.2021.152843.

16. Gamse JT, Gorelick DA. 2016. Mixtures, metabolites, and mechanisms: understanding toxicology using zebrafish. Zebrafish 13(5)377-378, PMID: 27618129, https://doi.org/10.1089/zeb.2016.1370.

17. Kimmel CB, Miller CT, Moens CB. 2001. Specification and morphogenesis of the zebrafish larval head skeleton. Dev Biol 233(2):239-257, PMID: 11336493, https://doi.org/10.1006/dbio.2001.0201.

18. van Boxtel AL, Pieterse B, Cenijn P, Kamstra JH, Brouwer A, van Wieringen W, et al. 2010. Dithiocarbamates induce craniofacial abnormalities and downregu-late sox9a during zebrafish development. Toxicol Sci 117(1):209-217, PMID: 20530235, https://doi.org/10.1093/toxsci/kfq169.

19. Hermsen SA,van den Brandhof EJ,van derVen LT, Piersma AH. 2011. Relative embryotoxicity of two classes of chemicals in a modified zebrafish embryotox-icity test and comparison with their in vivo potencies. Toxicol In Vitro 25(3):745-753, PMID: 21238576, https://doi.Org/10.1016/j.tiv.2011.01.005.

20. Bradford YM, Toro S, Ramachandran S, Ruzicka L, Howe DG, Eagle A, et al. 2017. Zebrafish models of human disease: gaining insight into human disease atZFIN. ILAR J 58(1):4-16, PMID: 28838067, https://doi.org/10.1093/ilar/ilw040.

21. Javidan Y, Schilling TF. 2004. Development of cartilage and bone. Methods Cell Biol 76:415-436, PMID: 15602885, https://doi.org/10.1016/s0091-679x(04)76018-5.

22. European Commission. 2019. On the protection of animals used for scientific purposes, Directive 2010/63/EU of the European Parliament and of the Council. Brussels, Belgium: European Commission.

23. Sukardi H, Chng HT, Chan EC, Gong Z, Lam SH. 2011. Zebrafish for drug toxicity screening: bridging the in vitro cell-based models and in vivo mammalian models. Expert Opin Drug Metab Toxicol 7(5)579-589, PMID: 21345150, https://doi.org/10. 1517/17425255.2011.562197.

24. Zoupa M, Zwart EP, Gremmer ER, Nugraha A, Compeer S, Slob W, et al. 2020. Dose addition in chemical mixtures inducing craniofacial malformations in zebrafish {Danio rerio) embryos. Food Chem Toxicol 137:111117, PMID: 31927004, https://doi.Org/10.1016/j.fct.2020.111117.

25. van Oostrom CT, Slob W, van der Ven LT. 2020. Defining embryonic developmental effects of chemical mixtures using the embryonic stem cell test. Food Chem Toxicol 140:111284, PMID: 32205227, https://doi.Org/10.1016/j.fct.2020. 111284.

26. Kortenkamp A, Faust M, Scholze M, Backhaus T. 2007. Low-level exposure to multiple chemicals: reason for human health concerns? Environ Health Perspect 115(suppl 1):106-114, PMID: 18174958, https://doi.org/10.1289/ehp.9358.

27. Heusinkveld HJ, Schoonen WG, Hodemaekers HM, Nugraha A, Sirks JJ, Veenma V, et al. 2020. Distinguishing mode of action of compounds inducing craniofacial malformations in zebrafish embryos to support dose-response modeling in combined exposures. Reprod Toxicol 96:114-127, PMID: 32553615, https://doi.Org/10.1016/j.reprotox.2020.06.002.

28. Lhasa Limited. 2020. Derek Nexus v.6.0.1, Nexus v.2.2.2 (Build 282, July 2018) with Derek Knowledgebase, https://www.lhasalimited.org/products/derek-nexus. htm [accessed 23 November 2021].

29. Schulte KW, Green E, Wilz A, Platten M, Daumke 0. 2017. Structural basis for aryl hydrocarbon receptor-mediated gene activation. Structure 25(7):1025-1033.e3, PMID: 28602820, https://doi.Org/10.1016/j.str.2017.05.008.

30. Klaholz B, Moras D. 2002. C-H...0 hydrogen bonds in the nuclear receptor RARgamma-a potential tool for drug selectivity. Structure 10(9):1197-1204, PMID: 12220491, https://doi.org/10.1016/S0969-2126(02)00828-6.

31. Chemical Computing Group. 2019. Molecular Operating Environment: Integrated Computer-Aided Molecular Design Platform, version 2019. https:// www.chemcomp.com/Products.htm [accessed 23 November 2021].

32. Cassano A, Manganaro A, Martin T, Young D, Piclin N, Pintore M, et al. 2010. CAESAR models for developmental toxicity. Chem Cent J 4(suppl 1):S4, PMID: 20678183, https://doi.org/10.1186/1752-153X-4-S1-S4.

33. Wu S, Fisher J, Naciff J, Laufersweiler M, Lester C, Daston G, et al. 2013. Framework for identifying chemicals with structural features associated with the potential to act as developmental or reproductive toxicants. Chem Res Toxicol 26(12):1840-1861, PMID: 24206190, https://doi.org/10.1021/tx400226u.

34. Benfenati E, Manganaro A, Gini G. 2013. VEGA-QSAR: Al inside a platform for predictive toxicology. Proceedings of the 2nd Workshop on Popularize Artificial Intelligence (PAI2013). 5 December 2013. Turin, Italy: Ceur-WS.

35. VEGA Software. 2020. VEGA QSAR Download. Version 1.1.5 48. https://www. vegahub.eu/download/vega-qsar-download/ [accessed 23 November 2021].

36. U.S. Environmental Protection Agency. 2021. Toxicity Estimation Software Tool (TEST): A program to estimate toxicity from molecular structure. https://www. epa.gov/chemical-research/toxicity-estimation-software-tool-test [accessed 23 November 2021].

37. Piersma AH, Bosgra S, van Duursen MBM, Hermsen SAB, Jonker LRA, Kroese ED, et al. 2013. Evaluation of an alternative in vitro test battery for detecting reproductive toxicants. Reprod Toxicol 38:53-64, PMID: 23511061, https://doi.org/ 10.1016/j.reprotox.2013.03.002.

38. Heusinkveld HJ, Zwart EP, de Haan A, Braeuning A, van der Ven LTM, Alarcan J. In press. Dose-addition in chemically-induced steatosis in zebrafish embryos.

39. Petzold AM, Bedell VM, Boczek NJ, Essner JJ, Balciunas D, Clark KJ, et al. 2010. SCORE imaging: specimen in a corrected optical rotational enclosure. Zebrafish 7(2):149-154. Jun, PMID: 20528262, https://doi.org/10.1089/zeb.2010.0660.

40. Hermsen SA, Pronk TE, van den Brandhof EJ, van der Ven LTM, Piersma AH. 2012. Triazole-induced gene expression changes in the zebrafish embryo. Reprod Toxicol 34(2):216-224, PMID: 22664267, https://doi.0rg/IO.IOI6/j.reprotox. 2012.05.093.

41. Abe M, Maeda T, Wakisaka S. 2008. Retinoic acid affects craniofacial patterning by changing Fgf8 expression in the pharyngeal ectoderm. Dev Growth Differ 50(9):717-729. Dec, PMID: 19046160, https://doi.Org/10.1111/j.1440-169X. 2008.01069.x.

42. Shimomura T, Kawakami M, Okuda H, Tatsumi K, Morita S, Nochioka K, et al. 2015. Retinoic acid regulates Lhx8 expression via FGF-8b to the upper jaw development of chick embryo. J Biosci Bioeng 119(3):260-266, PMID: 25239070, https://doi.Org/10.1016/j.jbiosc.2014.08.010.

43. Funato N, Nakamura M. 2017. Identification of shared and unique gene families associated with oral clefts. Int J Oral Sci 9(2):104-109, PMID: 28106045, https://doi.org/10.1038/ijos.2016.56.

44. Hutson MS, Leung MCK, Baker NC, Spencer RM, Knudsen TB. 2017. Computational model of secondary palate fusion and disruption. Chem Res Toxicol 30(4):965-979, PMID: 28045533, https://doi.org/10.1021/acs.chemrestox. 6b00350.

45. Menegola E, Broccia ML, Di Renzo F, Giavini E. 2006. Postulated pathogenic pathway in triazole fungicide induced dysmorphogenic effects. Reprod Toxicol 22(2): 186-195. Aug, PMID: 16781842, https://doi.Org/10.1016/j.reprotox.2006.04.008.

46. van Boxtel AL, Kamstra JH, Fluitsma DM, Legler J. 2010. Dithiocarbamates are teratogenic to developing zebrafish through inhibition of lysyl oxidase activity. Toxicol Appl Pharmacol 244(2):156-161, PMID: 20045017, https://doi.org/10. 1016/j.taap.2009.12.028.

47. Schmittgen TD, Livak KJ. 2008. Analyzing real-time PCR data by the comparative CT method. Nat Protoc 3(6):1101 1108, PMID: 18546601, https://doi.org/10. 1038/nprot.2008.73.

48. SlobW, Setzer RW. 2014. Shape and steepness of toxicological dose-response relationships of continuous endpoints. Crit Rev Toxicol 44(3):270-297, PMID: 24252121, https://doi.org/10.3109/10408444.2013.853726.

49. European Commission. 2012. Opinion on the Toxicity and Assessment of Chemical Mixtures. SCHER, SCCS, and SCENIHR. Brussels, Belgium: European Commission.https://ec.europa.eu/health/scientific_committees/environmental_ risks/docs/scher_o_150.pdf [accessed 23 November 2021].

50. Kortenkamp A. 2014. Low dose mixture effects of endocrine disrupters and their implications for regulatory thresholds in chemical risk assessment. Curr Opin Pharmacol 19:105-111, PMID: 25244397, https://doi.Org/10.1016/j.coph.2014. 08.006.

51. Conley JM, Lambright CS, Evans N, Cardon M, Medlock-Kakaley E, Wilson VS, et al. 2021. A mixture of 15 phthalates and pesticides below individual chemical no observed adverse effect levels (NOAELs) produces reproductive tract malformations in the male rat. Environ Int 156:106615, PMID: 34000504, https://doi.org/ 10.1016/j.envint.2021.106615.

52. Slob W. 2014. Benchmark dose and the three Rs. Part I. Getting more information from the same number of animals. Crit Rev Toxicol 44(7):557-567, PMID: 25000332, https://doi.org/10.3109/10408444.2014.925423.

53. Crepet A, Vanacker M, Sprong C, de Boer W, Blaznik U, Kennedy M, et al. 2019. Selecting mixtures on the basis of dietary exposure and hazard data: application to pesticide exposure in the European population in relation to steatosis. Int J Hyg Environ Health 222(2):291-306. Mar, PMID: 30579770, https://doi.org/10.1016/ j.ijheh.2018.12.002.

54. Kennedy MC, Garthwaite DG, de Boer WJ, Kruisselbrink JW. 2019. Modelling aggregate exposure to pesticides from dietary and crop spray sources in UK residents. Environ Sci Pollut Res Int 26(10):9892-9907, PMID: 30734911, https://doi.org/10.1007/s11356-019-04440-7.

55. Norwegian Institute of Public Health, Gkrillas A, Dirven H, Andreassen M, Papadopoulou E, Hjertholm H, et al. 2020. Risk assessment of phthalates based on aggregated exposure from foods and personal care products and comparison with biomonitoring data. EFSA J 18(suppl 1):e181105, PMID: 33294044, https://doi.Org/10.2903/j.efsa.2020.e181105.

56. HusoyT, Andreassen M, Hjertholm H, Carlsen MH, Norberg N, Sprong C, et al. 2019. The Norwegian biomonitoring study from the EU project EuroMix: levels of phenols and phthalates in 24-hour urine samples and exposure sources from food and personal care products. Environ Int 132:105103, PMID: 31470218, https://doi.Org/10.1016/j.envint.2019.105103.

57. Karrer C, Andreassen M, von Goetz N, Sonnet F, Sakhi AK, Hungerbiihler K, et al. 2020. The EuroMix human biomonitoring study: source-to-dose modeling of cumulative and aggregate exposure for the bisphenols BPA, BPS, and BPF and comparison with measured urinary levels. Environ Int 136:105397, PMID: 31884417, https://doi.Org/10.1016/j.envint.2019.105397.

58. Thepaut E, Dirven HAAM, Haug LS, Lindeman B, Poothong S, Andreassen M, et al. 2021. Per- and polyfluoroalkyl substances in serum and associations with food consumption and use of personal care products in the Norwegian biomonitoring study from the EU project EuroMix. Environ Res 195:110795, PMID: 33524335, https://doi.Org/10.1016/j.envres.2021.110795.

59. HusoyT, Martinez MA, Sharma RP, Kumar V, Andreassen M, Sakhi AK, et al. 2020. Comparison of aggregated exposure to di(2-ethylhexyl) phthalate from diet and personal care products with urinary concentrations of metabolites using a PBPK model-results from the Norwegian biomonitoring study in EuroMix. Food Chem Toxicol 143:111510, PMID: 32615240, https://doi.org/10. 1016/j.f ct.2020.111510.

60. Bosgra S, van der Voet H, Boon PE, Slob W. 2009. An integrated probabilistic framework for cumulative risk assessment of common mechanism chemicals in food: an example with organophosphorus pesticides. Regul Toxicol Pharmacol 54(2): 124-133, PMID: 19303907, https://doi.Org/10.1016/j.yrtph.2009.03.004.

61. Veltman CHJ, van der Ven LTM, Luijten M. In press. A framework to assess human relevance of alternative testing strategies-a case study with triazole exposure leading to craniofacial malformations.

62. Slob W. 2002. Dose-response modeling of continuous endpoints. Toxicol Sci 66:298-312, PMID: 11896297, https://doi.Org/https://doi.org/10.1093/toxsci/66.2.298.