Linking toxicant physiological mode of action with induced gene expression changes in Caenorhabditis elegans

Background Physiologically based modelling using DEBtox (dynamic energy budget in toxicology) and transcriptional profiling were used in Caenorhabditis elegans to identify how physiological modes of action, as indicated by effects on system level resource allocation were associated with changes in gene expression following exposure to three toxic chemicals: cadmium, fluoranthene (FA) and atrazine (AZ). Results For Cd, the physiological mode of action as indicated by DEBtox model fitting was an effect on energy assimilation from food, suggesting that the transcriptional response to exposure should be dominated by changes in the expression of transcripts associated with energy metabolism and the mitochondria. While evidence for effect on genes associated with energy production were seen, an ontological analysis also indicated an effect of Cd exposure on DNA integrity and transcriptional activity. DEBtox modelling showed an effect of FA on costs for growth and reproduction (i.e. for production of new and differentiated biomass). The microarray analysis supported this effect, showing an effect of FA on protein integrity and turnover that would be expected to have consequences for rates of somatic growth. For AZ, the physiological mode of action predicted by DEBtox was increased cost for maintenance. The transcriptional analysis demonstrated that this increase resulted from effects on DNA integrity as indicated by changes in the expression of genes chromosomal repair. Conclusions Our results have established that outputs from process based models and transcriptomics analyses can help to link mechanisms of action of toxic chemicals with resulting demographic effects. Such complimentary analyses can assist in the categorisation of chemicals for risk assessment purposes.


Background
Organisms display physiologies and life-histories that are determined by a set of interacting biochemical networks [1][2][3][4]. These in turn are controlled through regulation of the expression of proteins translated from RNA that itself is ultimately transcribed from the genome. Understanding how gene/protein networks control complex traits and their responses to external cues is one of the great challenges in biology [5][6][7][8]. Within this discipline of systems biology, emphasis is currently on coupling bottom-up approaches such as qualitative and quantitative detection of DNA/RNA, proteins and metabolites with data mining and visualisation approaches, thereby allowing elucidation of the networks that control complex biological processes [9][10][11][12]. The hypotheses that are generated through this approach can then be validated by focussed studies using loss or gain of function mutants and/or RNAi.
Whilst bottom-up techniques can aid understanding of how gene and protein networks control complex traits, such as those associated with life-history, we are far from fully elucidating the mechanisms through which life-history phenotypes change in response to environmental stimuli. In the case of environmental stress (e.g. toxic chemical exposure), an organism will typically respond by upregulation of a network of protective (e.g. detoxification) and compensatory (e.g. stress proteins) mechanisms [13]. These adaptive changes can alter energy allocation between biological processes, which in conjunction with direct toxic effects, may result in changes in the rate and timing of life-history events.
In physiological ecology, it is well known that life-history traits are not independent but instead are tightly linked as described through rules of metabolic organisation. As an example, a change in feeding rate will have very specific consequences for growth, development and reproduction because organisms are subject to the law of conservation of mass and energy. To understand the basis of such effects, theoretical ecologists have sought to develop top-down models that use phenotypic data to provide an understanding of the physiological mode of action of stress effects.
One of the most prominent ecophysiological systems for modelling the inter-relational links between life-history traits is based on dynamic energy budget (DEB) theory [14,15]. DEB theory explains how organisms acquire energy and allocate resources to key physiological processes such as maintenance, growth and reproduction based on a set of rules for metabolic organisation. The physiological processes that are modelled within the DEB framework all have a clear physiological basis. For example, animals obtain energy from food via the oxidation and reduction reactions of glycolysis, the Krebs cycle, and protein and fatty acid metabolism. Somatic maintenance, i.e. tissue and cell regeneration and recovery from injury, involves stress responses such as protection against oxidative stress, protein chaperoning and DNA repair. Somatic growth occurs through cell division, while reproduction generally requires gamete and egg production. In this way, ecophysiological models such as DEB can provide a topdown approach to predict toxicant impacts on the physiology of individuals that when combined with the bottom-up molecular approach can enhance our knowledge of systemic responses to environmental stress [16][17][18].
Since its inception, DEB theory has been used to describe how toxicants can alter energy usage in an attempt to bridge the gap between individual and ecosystem-based approaches to toxicology [16]. By examining responses of individuals to chemical exposure in the context of their energy balance, DEB theory can be used to explain how chemical toxicokinetics affects the allocation of energy resources; the extent of associated cellular damage and ultimately the performance of a range of life-history traits [17]. DEBtox can be used to derive a number of parameters including the no effect concentration (NEC), parameters describing chemical kinetics and a killing rate that is linked to internal concentrations [17,19]. All of which provide useful information about mechanisms of effect.
Here we use DEB theory and explicitly a recently developed version of the DEBtox model [17,19], to investigate the physiological mode of action of three different types of chemical stressors (the metal cadmium, the non-polar organic fluoranthene, and the herbicide atrazine) during a whole life-cycle exposure of the nematode Caenorhabditis elegans. The three chemicals were selected as representatives of major pollutant subgroups: cadmium is teratogenic, mutagenic and potentially carcinogenic metals; fluoranthene a non-polar polycyclic aromatic hydrocarbon (PAH) which can be expected to show mainly narcotic toxicity and atrazine a herbicide that is known to inhibit electron transport. Transcriptional profiling is then used to identify the changes in gene expression associated with toxicity that underlie observed changes in life-history traits (experiments designed according to the scheme shown in Fig.  1). The aim was to associate gene expression change for each chemical with the physiological mode of action as indicated by DEBtox modelling. Thus, a comprehensive picture of responses to exposure from transcriptional remodelling through to life-history change in worms exposed to different types of chemical stress was sought.

Life-history responses
Comparable values for mean lifespan of control worms (9.0, 8.8 and 9.7 days for Cd, FA and AZ) and control brood size (203, 204 and 211 eggs/worm for Cd, FA and AZ) were found between the three experiments. For Cd, DEBtox fits indicated that exposure affected survival patterns in time ( Fig. 2A). Worms exposed to low Cd concentrations showed an increase in lifespan (significant at 1.69 mg/L, Tukey's post-hoc, p < 0.001), although a non-significant reduction was found at high concentrations. DEBtox indicated an effect of Cd on temporal patterns and rates of egg production (Fig. 2, Panel B) and growth (Fig. 2, Panel C). Total broodsize at 9 mg/L was reduced by 93% compared to controls (Fig 2, Panel B) and the total broodsize EC 50 for Cd was 3.93 (95%: CI 3.15-4.94) mg L -1 . A reduction in growth rate and final body length was apparent in worms exposed to >2.5 mg/L of Cd (Fig. 2, Panel C).
DEBtox fits showed no clear dose associated effect of FA on survival patterns over the exposure period (Fig. 2, Panel D). FA did, however, have a negative affect on egg production (Fig. 2, Panel E), although an EC 50 for broodsize could not be calculated because the maximum reduction found was only 42% at 4 mg/L (n.b. 37% and 39% at 1 and 2 mg/L respectively). FA had an effect on growth patterns over time, with model fits describing slower growth at high FA concentrations, although ultimately worms at all concentrations reached the same body size. (Fig. 2, Panel F).
DEBtox fits indicated that AZ affected survival over time, with steepest mortality curves found at the highest exposure concentrations (Fig. 2, Panel G). DEBtox also highlighted that AZ affected temporal patterns of egg production (Fig. 2, Panel H). Greatest reduction of 56% was at 300 mg/L and the EC 50 for broodsize was 224 (175 -273) mg L -1 . Temporal patterns of weight change indicated slow growth and at higher concentrations and a reduction in final body size (Fig. 2, Panel I).
All growth, reproduction and lifespan data for all exposure concentrations and for all measured time points for DEBtox parameterisation for each compound were used (Fig. 2). Parameter estimates and the best-fitting physiological mode of action are given in Table 1. Log-likelihood ratio values for each of the fitted models are shown in Additional file 1. These values indicate the fit of the different tested models to the data, with the model showing the lowest log-liklihood ratio value having the best fit.
From the DEBtox fits, a set of single parameters for each physiological (e.g. growth rate, maximum reproduction rate), aging related (e.g. damage killing rate, damage amplification) and toxicological (elimination rate, killing rate) parameters were estimated for each of the tested chemicals (Table 1). Of these perhaps the most toxicologically meaningful at the system level is the likely physiological (energetic) mode of action of the compound. This parameter is determined by separate fitting of the DEBtox model for each possible physiological mode of action to identify the effect that gives the best fit of the model to the measured data (see Additional file 1). For Cd, the physiological mode of action that best fitted the data was an effect on energy assimilation. For FA, the best fitting mode of action was costs for growth and reproduction; while for AZ cost for maintenance of somatic tissue gave the best fit (Additional file 1). The DEBtox modelling, thus, suggests that the three tested chemicals which cause some similar life-history changes (e.g. concentration dependent reduction in egg production) and also some distinct changes (e.g. strongly reduced final body size for Cd, reduced growth rate to similar final body size for FA), do so through different physiological mechanisms. The mechanistic aetiology of the physiological changes in resource allocation associated with exposure to each chemical, thus, became the focal point of the transcriptional analysis.
Transcriptional profiling by whole genome microarray Successful and repeatable hybridization was indicated by a high correlation (r 2 > 0.9) and linear relationship  . Initially the effect of exposure to each chemical on the transcriptome was analyzed by PCA of each unfiltered normalized data-set. Distribution and clustering of samples within PCA models indicated that there was no specific impact of exposure on between sample variation for any chemical (Fig. 3). For Cd, the scores plots for PC1 and PC2 of a PCA model with an R 2 = 0.4 indicated a separation of the control and 1.25 m g/L samples from the three highest treatments along an axes across the first two PCs (Fig. 3a). The scores plots for PC1 and PC2 of a PCA model for FA with an R 2 = 0.39 indicated a partial separation of samples along the first PC, with the group of samples exposed to 0.5 mg/L identified as outliers (Fig. 3b).
Re-analysis of the data-set excluding the 0.5 mg/L samples produced a weaker PCA model (R 2 = 0.31), this model did, however, support the presence of a concentration related separation of transcript profiles (see Additional file 3). Because, however, all samples exposed to 0.5 mg/L showed a high similarity and there was no evidence through normalization of a systematic bias in the preparation or hybridization of these arrays, there is no clear and justified reason to exclude these transcript profiles from further analysis. The PC1 and PC2 score plots of a PCA model for AZ with an R 2 = 0.31 showed that the control, 5 and 25 mg/L exposed worms were separated from the 75 and 150 mg/L atrazine exposed worms along PC1 (Fig. 3c).
To identify gene expression changes that were associated with the toxicosis of the three chemicals, detailed analysis was undertaken to identify pathways linked to chemical effects at the test concentration that approximated most closely to the EC 50 for total broodsize for each chemical (2.25 mg/L Cd, 1 mg/L for FA, 150 mg/L for AZ). The predominant approach used was an analysis based on annotation enrichment within a set of differentially regulated genes. This was designed to allow analysis of the changes in expression that were associated with the putative physiological mechanisms of action identified by DEBtox modelling.
For Cd, a PLS-DA analysis of normalised individual gene expression levels for controls and 2.25 mg/L Cd exposed worms indicated a distinct separation of the samples along PC1 of a highly robust model (Q 2 = 0.753) (see Additional file 4a). PLS-DA analysis of the relative percentage of changing genes by biological process GO term indicated a separation of control and four successfully hybridised 2.25 mg/L exposed samples within a highly robust model (Q 2 = 0.99) (Fig. 4a). GO terms highly correlated with PC1 included terms associated with cellular organisation (GO:0016043), DNA metabolism (GO:0006259) and also further processes downstream of these effects associated with cell cycle (e.g. GO:0007049) and reproduction (GO:0000003). An indication of effect on energy metabolism was evident by changes in the expression of genes associated with lipid metabolism (GO:0006629)( Table 2). A second approach to annotation enrichment analysis used the Database for Annotation, Visualization and Integrated Discovery (DAVID) tool. This web based resource provides a set of functional annotation tools for identification of biological meaning within  differentially regulated gene lists. DAVID input for Cd was a list of 2326 genes identified as significantly up or down regulated between the 0 and 2.25 mg/L Cd exposure without multiple sample correction (against an expectation of 937 by chance) (n.b. with false discovery correction [20] applied, only one gene (predicted gene ID Y59A8B.12) of unknown function was significant). The list without multiple sample correction was used as DAVID input to limit type II errors in significant gene identification. The fact that the DAVID pathway enrichment analysis system itself includes an algorithm for multiple sample correction effectively limits type 1 errors during identification of enriched pathways. The DAVID analysis for Cd identified three KEGG pathways that were significantly over represented in the 0 vs 2.25 mg/L gene-list. These related to glutathione metabolism (5 of 12 present genes), taurine and hypotaurine metabolism (2 of 2 present genes) and protein export (4 of 10 present genes). DAVID analysis also identified a number of significantly enriched terms contained within the annotation clusters (Table 3). A number of the most highly enriched of these clusters, such as membrane association, cation binding, glycoprotein signalling, transcription and protein kinase activity were also significantly enriched in both the FA and AZ exposed worms. Cd also shared two annotation clusters containing significant terms with FA. These related to EGF domains and isomerase activity. Four annotation clusters were shared between Cd and AZ exposed worms, of which protein chaperoning and mitochondria associated genes were the most strongly enriched. Three of the significant annotation clusters were unique to Cd exposed worms. These related to transcription factor activity, transmembrane transport and thiol proteases.
Expression patterns of individual genes that were associated with the significant annotation terms were investigated across the full exposure range, This was done to identify a delineative sub-set of annotated changing genes that showed a significant linear relationship with Cd exposure concentration as opposed to only possibly idiosyncratic expression change at the approximate EC 50 concentration for which the gene list for pathway analysis was generated (Fig. 5). Such gene level scrutiny, therefore, provided insight into the mechanistic effects of Cd exposure that allowed a deeper assessment of the transcriptional changes associated with Cd exposure.
One category of genes within which a number of transcripts showed a significant linear relationship of expression with Cd exposure was the transcription Figure 4 Partial least squares discriminant of percentage changing genes per GO term at concentrations close to reproduction EC 50 . Scores plot for PC1 and PC2 from a PLS-DA analysis conducted using the percentage of genes associated with each GO term changing in expression by 1.5 fold at a concentration approximating to the brood size EC 50 when compared to controls for (A) Cd, (B) FA, and (C) AZ. Only terms represented by 10 or more genes were used for the analysis. Table 2 GO term classifications correlated with PC1 of a PLS-DA analysis. assessed by GO term classification counts of the 20 Gene Ontology Biological process terms most strongly correlated with PC 1 within a PLS-DA analysis that separates control worms from worms exposed to an approximate EC50 for each of the three tested chemicals. factors. Genes within this category showed both significant up and down regulation. These included Zn finger type factors such as athp-1 and also nuclear hormone receptors such as nhr-162, nhr-268 and nhr-274 ( Fig 5). The effect of Cd exposure on these multiple transcriptional factors was indicative of a wider remodelling of the transcriptome as shown by the changes in expression profiles across treatments illustrated within the PCA (Fig. 3a). Among the metal binding genes responding to Cd were mtl-1 and mtl-2, both of which code putative metallothioneins. These two genes were upregulated in all Cd treatments (Fig. 5). Expression of a number of structural proteins was also related to Cd exposure. These include the cadhedrin cdh-5 which was down-regulated and the tektin cytoskeleton protein ecoded by sequence R02E12.4 which was strongly upregulated. Expression of multiple meprin metalloproteases (e.g. nas-18) and E3 ubiquitin ligases that are involved in protein tagging and degradation were also significantly affected by Cd ( Fig 5). Another category of genes highlighted as being responsive to Cd exposure were protein kinases.
Here of nine significant genes, eight such wnk-1 were upregulated. Expression of genes with specific involvement in mitochondrial biogenesis and functioning, such as dap-3 and ADP/ATP carrier protein was also significantly related to Cd concentration indicating an effect of exposure on mitochondrial integrity. For FA exposed worms, PLS-DA analysis of individual normalised gene expression data of the controls and 1 mg/L exposed samples indicated a clear separation along PC1 of a highly robust model (Q 2 = 0.993) (Additional file 4b). The PLS-DA analysis of relative percentage of changing genes per GO biological process term, also showed a clear separation of controls from the 1 mg/L exposed samples (Q 2 = 0.99) (Fig. 4b). GO biological processes represented among the terms most negatively correlated with PC1 for FA exposed worms included ontologies associated with catabolism (GO:0009056) and cell organisation and biogenesis (GO:0016043): terms which were also identified by the same analysis for Cd and AZ ( Table 2). Terms that were particular associated with FA exposure included protein trafficking (GO: 0015031) and carbohydrate metabolism (GO:0005975).
Functional annotation analysis using DAVID was conducted with a list of 1807 genes significantly differentially expressed (p < 0.05) without multiple sample correction between the control and 1 mg/L FA treatments. False discovery correction according to Benjamini and Hochberg [20] produced a list of 38 significant transcripts. Based on the uncorrected list, DAVID identified six KEGG pathways that were significantly over represented in differentially expressed genes. These related to glutathione metabolism (5 of 12 genes present), metabolism of xenobiotics by P450 (16 of 22 genes present), lycine, serine and threonine metabolism (6 of 20 genes present), nitrogen metabolism (4 of 10 genes present), alkaloid biosynthesis (2 of 2 genes present) and the ribosome (20 of 74 genes present).  Figure 5 Life-history and expression of key genes in Cd exposed nematodes. Transcriptional responses for a series of selected genes associated with ontology clusters significantly enriched for differentially expressed genes between adult C. elegans exposed to a control and 2.25 mg Cd/L agar. Array data were normalised per chip and per gene median polishing and are expressed relative to the control samples. Error bars indicate the SD of the mean.
DAVID analysis also identified a number of annotation clusters that contained a significant overrepresentation of genes differentially expressed following FA exposure (Table 3). In addition to clusters shared with both Cd and AZ exposed worms and also with Cd exposed worms as outlined above, annotation terms were identified that were shared only by FA and AZ exposed worms. These included annotations associated with DNA structure and also haem/monooxygenases: which suggest a role of cytochrome P450s in the metabolism of these two organic chemicals. Six annotation clusters with significant terms were also identified that were unique to the FA exposure. Among these oxidoreductases were most enriched, while clusters associated with both mRNA transcript processing and the biosynthesis of amino acid and proteins that are symptomatic of an effect of FA on macromolecule turnover were also significantly enriched.
For annotation terms identified as significantly enriched for genes differentially regulated by FA, patterns of expression of individual genes associated with each term were investigated across the exposure range (Fig. 6). A number of the genes showing dose related expression change were involved in metal ion binding (a common term with Cd and AZ exposed worms). Genes in this category responsive to FA included transcripts associated with macromolecule turnover, such as the fatty acid desaturase fat-6 and the meperin metalloprotease nas-6, immunoglobulins such as him-4 and substrate transport, such as the P-glycoprotein ABCtransporter pgp-9 (Fig 6). This latter gene is known to be involved in xenobiotic transport within the cell. A further category of genes responsive to FA, as well as Cd and AZ, exposure were the protein kinases. Three genes in the category were down regulated,(e.g. cAMPdependent protein kinase) and two upregulated (e.g. Serine/threonine protein kinase) at 1 mg/L FA. Genes involved in oxidoreductase activity differentially expressed following FA exposure included transcripts coding flavin containing monoxygenase that are involved in oxidative stress protection, as well as downregulated genes involved in dehydrogenation such as D-aspartate oxidase. Genes involved in mRNA processing that were differentially expressed included multiple small nuclear ribonucleoprotein (e.g. snr-4 and snr-5), all of which were upregulated at 1 mg/L FA (Fig. 6). For AZ exposed worms, PLS-DA analysis of normalised gene expression data for the controls and 150 mg/L AZ exposed worm indicated a separation along PC1 of a highly robust model (Q 2 = 0.998) (Additional file 4c). PLS-DA analysis of the relative percentage of changing genes per biological process GO term also indicated a clear separation of control and 150 mg/L exposed samples along PC1 (Q 2 = 0.985) (Fig. 4c). Analysis of PC1 loading to identify the most negatively correlated terms indicated a strong representation of terms associated with nucleic acid metabolism (GO:0006139) and organelle organisation and biogenesis (GO:0006996), a response that was shared with Cd exposed worms ( Table 2). Unique terms that were associated with AZ exposure included protein biosynthesis (GO:0006412) and metabolism (GO:0019538), as well the metabolism of individual amino acids (GO:0006519). These responses are all indicative of a potential effect of AZ on protein turnover.
Functional annotation analysis using DAVID was conducted using a list of 1333 genes significantly differently expressed (p < 0.05, without false discovery correction) between the 0 and 150 mg/L AZ treatments. Multiple sample correction according to Benjamini and Hochberg [20] reduced this list to 47 significant genes. From the longer gene list, DAVID analysis identified three KEGG pathways that were significantly over represented. These were glycolysis (7 of 19 present genes), pyruvate metabolism (8 of 14 present genes) and the ribosome (18 of 74 present genes).
A number of significantly enriched annotation clusters were also identified by DAVID analysis (Table 3). Clusters shared with both Cd and FA, for which AZ worms also showed a high level of enrichment of differentially regulated genes were for cation binding, protein synthesis and intracellular signalling. Enriched terms shared with Cd exposed worms were associated with mitochondrial function, body patterning and stress and heat shock, while those shared with FA were associated with oxidoreductase activity and DNA integrity. Seven enriched annotation clusters were uniquely associated with AZ exposure. Notable among these was larval development which was the most strongly enriched cluster among the significant terms. Glycolysis, and also terms associated with protein tagging and degradation were also significant.
Specific genes that were associated with annotation terms highlighted for AZ that showed a significant linear relationship of expression with exposure concentration included genes associated with metal ion binding such as numerous zinc finger transcription factors, as well as transcripts encoding proteins involved in protein tagging and degradation such as Meperin metalloproteases and ubiquitin ligases (Fig. 7). Protein kinases responding to AZ exposure, included both down regulated genes, such as the cdk activating kinase encoded by sequence H01G02.2 and upregulated genes such as the serine/ threonine protein kinase nekl-1. A large number of changing genes were associated with the larval development term. These included histones, transcription factors such as taf-1 which encodes a factor that possesses histone acetyl transferase (HAT) activity as well as Figure 6 Life-history and expression of key genes in FA exposed nematodes. Transcriptional responses for a series of selected genes associated with ontology clusters significantly enriched for differentially expressed genes between adult C. elegans exposed to a control and 1 mg FA/L agar. Array data were normalised per chip and per gene median polishing and are expressed relative to the control samples. Error bars indicate the SD of the mean.  Life-history and expression of key genes in AZ exposed nematodes. Transcriptional responses for a series of selected genes associated with ontology clusters significantly enriched for differentially expressed genes between adult C. elegans exposed to a control and 150 mg AZ/L agar. Array data were normalised per chip and per gene median polishing and are expressed relative to the control samples. Error bars indicate the SD of the mean. genes involved in mRNA splicing, the control of apoptosis and macromolecule (e.g. lipid) metabolism (Fig. 7). A number of changing genes were also associated with nuclear activity, including numerous transcription factors, such as lin-11 and nuclear hormone receptors, such as nhr-23. Differentially expressed monooxygenases included peptidylglycine alpha-amidating monooxygenase (Fig 7), which showed a significant increase with exposure, indicating an effect of AZ on the nematode p450 system.

Discussion
The genetic control of life-history and the responses of traits to stress is a fundamental aspect of biology that is of physiological and ecological significance [2,21]. The approach taken by molecular geneticists to understand life-history events like growth, reproductive timing, reproductive output and life-span has mostly been through reverse genetics. In particular, thanks to deletion strains and the use of knock down methods, it has been possible to identify a range of genes and small metabolite determinants that effect the life-history of laboratory model species such as Saccharomyces cerevisiae, Drosophila melanogaster, Arabidopsis thaliana and C. elegans, [3,[22][23][24]. As an example, loss of function mutants for genes in the insulin signaling pathway have been shown to extend lifespan in C. elegans [25,26] and also confer reduced sensitivity to environmental stress [27][28][29][30]. This suggests the presence of a functional link between stress responses, such as the processing and detoxification of toxic chemicals, and the timing of lifecycle events including death [31].
For ecologists, the notion of an association between (chemical) stress exposure and the control and timing of life-history events is consistent with existing and well developed concepts. From an ecological perspective, lifehistory can be understood through the use of physiologically based models that describe how organisms acquire resources from food and subsequently allocate these to key biological processes, such as growth, maintenance, and reproduction [14]. Further, by linking rules for metabolic organization within dynamic energy budget theory with pharmacokinetic models that describe chemical uptake and elimination, it is possible to link the hazards associated with chemical exposure to resulting toxic effects [32]. This link is described within the DEBtox framework which provides a physiologically based context for understanding the mechanistic and energetic basis of toxicant effects.
Here, data from three continuously monitored lifecycle toxicity tests with C. elegans were used to parameterize DEBtox models from which a number of biologically meaningful parameters describing toxicity, including toxicokinetic and life-history associated parameters and physiological mode of action, were estimated. The life-cycle toxicity tests themselves and subsequent DEB modeling indicated that the three tested chemicals each had effects on a range of life-history traits. Further, even though time dependent patterns of survival, growth and reproduction were affected by each chemical, the DEBtox analysis indicated that the underlying physiological mechanisms causing these effects were different in each case.
To provide initial insight into the energetic basis of life-history effects, DEBtox model fits were used to allow calculation of a number of physiological, ageing related and toxicological parameters ( Table 1). Of these perhaps the most toxicologically meaningful is the physiological (energetic) mode of action of the compound which is determined by serial fitting of separate models. For Cd the physiological mode of action that best fit the data was an effect on energy assimilation. This supports the results of previous DEBtox analysis of Cd effect in C. elegans [33,34], as well as the results of similarly designed studies of Cd toxicity for the springtail Folsomia candida [18], and Daphnia magna [35]. An effect of Cd on energy acquisition in C. elegans is also consistent with the observation of reduced feeding following Cd exposure [36,37].
For FA, the best fitting physiological mode of action within DEBtox was costs for growth and reproduction. This mode of action suggests that the production of somatic and reproductive tissue becomes most costly due to FA exposure, meaning that growth (and egg production) is slowed, but that also ultimately the same final body size of adult worms is reached. These model predictions for the an effect of FA exposure on costs for growth and reproduction are consistent with the mode of action observed for pentachlorobenzene in C. elegans [37], suggesting this may be a generic physiological mode of action for narcotic organic chemicals in this species.
For AZ, the best fitting physiological mode of action within DEBtox was costs for somatic and/or maturity maintenance. These results, to our knowledge, represent the first DEBtox based assessment of the physiological mode of action of any triazine herbicide. This mode of action suggests that AZ exposed nematodes are required to divert an increased proportion of available energy to protect or repair somatic and reproductive tissues.
To address the specific mechanistic aetiology of the physiological changes in resource allocation associated with the toxicity of each chemical, transcriptional profiling using a whole genome oligonucleotide microarray was used. This was not the first study that has used microarray analysis to study the effects of chemical exposure on the C. elegans transcriptome [38][39][40][41]. Uniquely, however, in this study we were able to link observed transcriptional changes with modelled energy budget responses as detailed through DEBtox fits. This combined analysis provides better phenotypic anchoring for the interpretation of the transcriptional changes that are associated with chemical exposure.
Of course while the DEBtox analysis did indicate particular chemical effects, some common responses were seen as reflected by the identification of a number of common pathways and annotation terms that were enriched for genes differentially regulated following exposure to all chemicals. These common enriched annotations included membrane proteins, cation binding, protein kinase activity/signal transduction, ribosomal proteins and protein degradation. Many of these pathways, such as changes in cation trafficking [42], mitogen-activated protein kinase signalling [39,43] and protein chaperoning [40,44], have been previously identified as responsive to chemical exposure in C. elegans. Moreover in yeast, studies have identified a group of "common environmental response" genes, the expression of which were ubiquitously altered by a range of environmental stressors. When viewed at the pathway rather than individual gene level [45], many of the common responses found in this study are consistent with pathways such as heat shock, organic compound metabolism, oxidative stress, protein turnover, metal transport, thioredoxin/glutathione regulation and energy generation that were highlighted as ubiquitously stress responsive in yeast [46,47].
In addition to transcriptional responses shared between all chemicals, further differentially expressed gene enriched annotation terms were also identified by DAVID analysis that were common between particular chemical pairs. These included the KEGG pathway for glutathione metabolism, a peptide that is involved in protection of cells from chemical injury resulting from oxidative stress [13] for both Cd and FA; and haem/ monoxygenase metabolism for the FA and AZ. The later indicates a likely role of cytochrome P450 in the detoxification of these two organic chemicals. For FA this potential role for P450 was supported by the differential regulation of many (16 of 22) genes in C. elegans that are associated with the KEGG pathway for metabolism of xenobiotics by P450. Both the glutathione and monooxygenase pathways have been previously implicated as core to the responses of many species to a range of chemical stressors [13].
As well as common responses, unique transcriptional changes were also seen that when culminating in overrepresentation at the pathway level potentially contribute to the specific aetiology of the life-history effects of each chemical. For Cd, DEBtox indicated that the most plausible mode of action was a reduction in energy assimilation. At the mechanistic level, the most obvious manifestation of this effect would be down-regulation of biological processes involved in sugar and lipid metabolism, electron transport, ATP turnover and mitochondrial function. Although the mitochondria were highlighted as a target of Cd (and also FA), there were other annotation clusters that actually showed a higher degree of enrichment for differentially expressed genes in Cd exposed worms. One such category of genes were associated with cation binding. These included transcription factors and structural proteins, as well as the two putative metallothionein mtl-1 and mtl-2 which have been functionally linked to Cd responses in C. elegans [48,49]. Other gene expression response associated with Cd exposure included membrane associated genes and also transcriptional activity, with nucleic acid metabolism also highlighted. Previous work has implicated Cd in DNA damage in earthworms [50], nematodes [43] and mammalian cell lines [51]. Further, a global down regulation of RNA production has also been identified in C. elegans exposed to Cd [52]. These findings, coupled with the indications from our transcriptional analyses, raise the possibility that Cd may affect DNA integrity and transcription, thereby causing a reduction in transcription and protein synthesis that ultimately leads to a suppression of energy demand. Such effects would be manifested as an overall change in energy assimilation, but with the effect coming from suppressed demand rather than reduced energy production.
For FA, the best fitting DEBtox physiological mode of action was costs for growth and reproduction. This suggests that worms exposed to FA may have to expend a greater amount of their available energy to the production of new biomass. The ontology term and pathways analysis for FA exposed worms provided convincing support for an effect of FA on such production costs. Thus DAVID analysis indicated a specific effect of FA on the synthesis of amino acids and proteins and a strong effect on the ribosome, the site of protein production; while the PLS-DA analysis highlighted effects on catabolism, including genes associated with the KEGG pathway for protein transport. These effects all suggest a possible effect of FA on protein synthesis and turnover.
A previous metabolomic study of the effects of pyrene (another PAH) on earthworms indicated increased protein breakdown as a result of decreased glucose and fatty acid metabolism [53] and similar responses have been seen in rats during toxicological studies [54]. This suggests that energy metabolism may be switched from carbohydrates to protein possibly as a result of an increase in the amount of free amino acid present in cells as a result on an increase in the degradation of protein compromised by chemical exposure. This will result in the need for frequent re-synthesis of key structural proteins during growth, development and egg production at a considerable energetic cost to the organism.
For AZ, the best fitting physiological mode of action within DEBtox was costs for maintenance of somatic and reproductive tissue. Without further insight it may be hypothesized that these costs are likely to be linked with an increase in the protection of the range of C. elegans proteins involved in detoxification, antioxidant protection and protein chaperoning [55]. While AZ worms did show significantly enriched differential expression of monooxygenases, including a number of genes involved in antioxidant defense, it was noticeable that many of the primary components of the antioxidant protection systems, such as SOD, catalase, glutatione peroxidase, peroxiredoxins, were not differentially expressed as a result of AZ exposure: at least at the transcriptional level at the sampled time-point.
Based on the annotation analysis of genes differentially expressed by AZ exposure, the most strongly enriched categories were transcripts associated with larval development. Within this category, many of the differentially expressed genes were linked with DNA integrity, such as the expression of histones, and genes involved in transcriptional regulation and mDNA splicing. Further within both the PLS-DA and DAVID analysis, an enrichment of differentially expressed genes associated with DNA integrity and repair was indicated. This is all consistent with previous work on atrazine which has shown a link between atrazine exposure and DNA damage [50,56,57]. This suggests that increased demands to maintain somatic and reproductive tissue may be associated not only with demands for organic compound metabolism, but also with the repair of DNA damage. Such repair, which is necessary to maintain DNA integrity and reduce rates of cellular apoptosis, requires the use of energy resources that ultimately need to be diverted from other processes.
To fully understand the effects of chemical stress on life-history responses, such as growth, reproduction and longevity, there is a need to link the effects of chemicals at the detailed mechanistic level with higher organisation changes in resource allocation and trait performance. Process based models such as DEBtox can provide a framework to understand the physiological basis of life-history in terms of energy allocation. Physiological mode of action predictions from DEBtox can provide a useful indication of the basis of the toxic effect. However, these modes of action still represent changes in very broad processes. Transcriptional profiling on the other hand can provide detailed understanding of the target site and resulting effects of a chemical, but not how these expression changes combine to cause effects on the biological response of major life-cycle traits for the whole organism. Combined, though, as has been done here, the two approaches can provide a unique view of both the range of gene expression response associated with chemical exposure and the culminating effects on resource allocation and life-history.
The work that is presented here illustrates a stepped approach that utilises both physiological based modelling and gene expression profiling to probe the energetic basis and aetiology of chemical effects on demographic traits. The analysis as conducted essentially represents a two stage process in which DEBtox analysis is first used to generate hypotheses concerning the physiological mode of action of the tested chemicals. Gene expression profiling is then used to investigate the validity of these hypotheses in terms of biological processes and individual genes.
Clearly the approach to data modelling and interpretation used represents only one of a number of possible ways that the comprehensive demographic and transcriptional data-sets can be analysed. One potential alternative approach could have been to include all expression data and life-history measurements within a single analysis to look for correlation between trait measurements and gene expression change. While this method could have been used to look for associations between gene expression and individual traits, such an approach is not yet possible within a DEBtox framework since all trait data is included within the model to generate only a single set of DEBtox parameters per chemical. To undertake an analysis that allows correlations between DEBtox parameters and gene expression, similar data-sets for further chemicals would be required to allow robust correlations to be made. To allow the wider toxicogenomic community the opportunity to apply different approach to the analysis of the results described here, all data has been made accessible in a MIAME compliant format using the advanced path of data dissemination, namely via MIAME/Env (for details see: nebc.nox.ac.uk/miame/miame_env.html). Storage in the MIAME/Env format allows unrestricted (open access) in support of further deep mining of the data.

Conclusions
Linking molecular mechanisms of toxicant to their physiological mode of action is not just of academic interest. Within applied toxicology understanding the cascade through which chemicals exert their effects can be important for toxicant categorisation, development of structure activity relationships, and identification of the most appropriate models for describing the joint effects of chemical mixtures [58]. Here we have demonstrated the value of adopting a combination of top down process based modelling and a bottom up mechanistic approaches to identify physiological modes of action of toxicants and associated gene expression changes. Through process based modelling it was identified that three xenobiotics from different chemical classes, produce some similar effects, but also caused specific lifehistory responses that could be linked through process modelling to effects on resource allocation. Transcriptional profiling allows the further identification of these effects in terms of induced gene expression change and their relationship with exposure concentration. Such detailed assessment of life-history effects and their causes can provide a platform for a more detailed categorisation of chemicals for appropriate risk assessment and a sound mechanistic basis for the identification of effect linked markers of chemical injury.

Experimental design
Three separated experiments were conducted, one with each chemical, with each experiment designed to measure both the effects of exposure on a range of life-history traits and also to obtain a sample of exposed nematode for transcriptional profiling (see Fig. 1 for schematic of experimental design). All experiments were conducted with C. elegans strain cib-1(e2300)I., a temperature sensitive mutant that is wild-type below 25°C, but displays an embryonic lethal phenotype at 25°C and above.
For the exposure worms were reared during pre-exposure phases at 15°C, a temperature with no effect on embryonic lethality. The main exposure was conducted using eggs synchronised by bleaching from the preexposed worms with the temperature switched to 25°C, at temperature at which embryonic lethality is induced within the first few cell divisions. when the resultant offspring reached L1 stage. This design meant that a cohort of worms could be reared throughout the main exposure for their entire life-span, thereby provide suitable data for DEBtox modelling, without the population used for transcriptional profiling becoming infiltrated by hatched and growing juveniles as would have been the case if a wildtype strain was used. The advantage of this was that the transcriptional profile gained from the analysis of exposed replicated cohorts of worms reared to adult could be linked to measurements of life history traits (growth, maturation, reproduction and mortality) recorded from similarly kept individuals without any influence of the presence of different worm stages in the sample on measured expression pattern. Thus, all samples generated differed only in relation to chemical exposure rather than stage specific effects within the population.

Life-cycle toxicity tests
Dosing of the NGM agar and OP50 were conducted using analytical grade cadmium chloride (Cd) (Sigma Chemicals, Poole, UK), fluoranthene (FA) (Fisher Chemicals, Loughborough, UK) and atrazine (AZ) (Sigma, UK). Stock solutions were prepared in bi-distilled water, ethanol and dimethyl sulfoxide (DMSO) respectively. Stock solutions were added to NGM cooled to 55°C, to give a total of 12 test concentrations per experiment. For each experiment, initially a stage synchronised population was obtained by bleaching culture plates of gravid adults. Aliquots of extracted eggs were then transferred to multiple NGM plates dosed to each test concentrations and previously inoculated with OP50. These eggs were allowed to hatch and the worms grown to adult at 15°C. Each exposed population was then synchronised again by bleaching. These resultant eggs (derived from life-time exposed parents) were then transferred to each of five replicate freshly dosed 90 mm diameter NGM plates for each of the 12 test concentrations per experiment (seven for FA controls). This time the eggs were hatched and at L1 stage switched to 25°C and grown to L2/L3 stage. For each test concentration, 48 individual worms were then taken at randomly from the five replicate plates per treatment and each transferred to a separate well of a 12-well plate dosed to the same test concentration. These plates were then kept also at 25°C and subsequently worms transferred daily to a fresh similarly dosed well. At each transfer, worm survival was checked and the number of eggs counted. Every 6 hours throughout the exposure, the growth of 10 random individuals per treatment was recorded by photographing worms using a Nikon Coolpix 4500 digital camera. Size was measured using Image Pro Express software (Media Cybernetics, Marlow, UK) and volumetric length (cubic root of body volume) calculated [34].

Microarray analysis
In parallel with monitoring of the single worms, the remaining approximately 10,000 worms that were on the five initial replicate large plates per treatment were kept also at 25°C and periodically checked to ensure they developed simultaneously with the individual worms at the same concentration. At 12 hours after onset of egg production, to ensure stage synchronisation across all samples, all replicate plates from five concentrations that spanned the exposure range used in the life-history study (Cd: 0, 1.125, 2.25, 4.5, 9 mg/L; FA: 0, 0.1, 0.25, 0.5, 1 mg/L; AZ 0, 5, 25, 75, 150 mg/L), were washed in M9 and worms collected by gravitational isolation. This provided a sample of exposed worms for each replicate. The worm samples were immediately frozen in liquid nitrogen and stored at -80°C.
Global RNA transcript abundance was analysed using oligo microarrays printed using C. elegans oligonucleotide set version 1.1 (Operon™). Probes were designed from the WormPep DNA Release 75 database, maintained and developed by the Sanger Institute (Cambridge, UK). This set comprises 19,873 70 mer oligonucleotides designed to the 3' region of all predicted genes. The oligonucleotides were printed onto UltraGAPS™ (Corning, Barry, UK) slides at a concentration of 15 μM. Appropriate landmarks of 8 replicates of the Lucida Scorecard (Amersham) were across slides sub-array. The Lucida Scorecard is a selection of heterologous gene reporters which show no cross reactivity to C. elegans transcripts. Analysis of calibrators exploited 10 "alien" RNA spikes (components of the Amersham Lucidea Scorecard) introduced at known concentrations, between 1 pmol and 30 nmol, prior to labelling and each hybridized to replicate reporter spots on the array. Image analysis of the signals generated by these reporters was performed for each array to determine the sensitivity and relationship between RNA concentration and fluorescent signal. All oligonucleotides were immobilised by baking at 80°C for 2 hours and UV cross linked. Slides were blocked using a blocking solution of 1% Bovine Serum Albumin, 5 × SSC and 0.1% SDS (all Sigma, UK) at 42°C for 30 minutes, washed five times in sterile water and dried.
Sample RNA was isolated (Tri-reagent, Sigma, UK), cleaned (RNeasy Mini kit, Qiagen, Crawley, UK) and quantified spectroscopically. All experiments were conducted following a reference design with the reference sample compiled from a mixture of RNA extracted from control and Cd, FA and AZ exposed worms. Use of this reference was intended to provide optimal coverage of the spotted genes. For hybridisation, sample and reference mRNA was purified from 15 μg total RNA using Qiagen columns. RNAs were indirectly labelled with fluorescent dyes, reference Cy3 and samples Cy5 respectively and reverse transcribed using Moloney Murine Leukaemia Virus reverse transcriptase (ABgene, Surrey, UK). Cy3 and Cy5 labelled cDNA was purified using the CyScribe™ GFX Purification Kit (Amersham Biosciences, Buckinghamshire, UK). cDNA labelling quality was assessed by scanning 1 μl of sample run on a slide gel under a green (560 BP filter) and red (675 BP filter) laser on a GENETAC™ LS IV scanner (Genomic Solutions, Huntingdon, UK). Hybridisations used 20 pmol of Cy dye molecules for each channel at 42°C for 24 hours. After removal and washing, slides were scanned at 633 nm (Cy5) and 543 nm (Cy3) on a ScanArray™ Express HT microarray scanner (Perkin Elmer, Beaconsfield, UK).

Data modelling and statistical analyses
A modified version of DEBtox [32], a model based on DEB theory [14] was used to analyse effects on life history traits and to discern the physiological mode of action of each toxicant, as well as a range of other physiologically relevant parameters associated with the toxicokinetic and toxicodynamics of the chemical. Within the DEBtox framework five models have been defined, which represent different energetic modes of action through which chemicals may exert their toxic effect. These are assimilation (less energy acquired from food); maintenance (increased energy demand to maintain existing somatic tissue); costs for growth and reproduction (increased energy demand for production of new somatic tissue and reproductive structures); costs for reproduction (increased energy demand for production of progeny); and hazard to embryo (direct effect on reproduction such as a direct toxic effect on embryo survival).
To fit these variants, the basic DEBtox model described by Kooijman and Bedaux [17] was extended to deal with life-cycle experiments; thereby allowing simultaneous modelling of multiple traits (Jager et al., 2004). Modifications were needed to deal with the specific details of the C. elegans life cycle [19,33]. All DEBtox models were fitted using maximum likelihood estimation in Matlab® version 7.0 (Release 14). In all fits, reproduction data only from the period t < 6.5 days were included, as egg production typically stops after this time even in exposed worms and so little information is gained by inclusion of data from further time points. For fluoranthene all models were fitted without including data from the highest exposure concentration (4 mg/L) due to concerns that this concentration exceeded the water solubility limit for the chemical within the test system (crystals observed in system). Results from the fits of the five single DEBtox models were compared, with the lowest log-likelihood values taken to represent the most plausible physiological mode of action for each chemical.
The whole genome microarray data was imported into GeneSpring 7.3 (Agilent Technologies, Stockport, UK). Background measurements were removed and metadata associated with effects on individual worms at each exposure concentration assigned to each array. Normalisation was then performed using Edwards background subtraction and a tip lowess/scaled normalisation. Values were then divided against mean control expression. Quantity assessment was performed by visualising a box plot of the normalised data. This allowed identification of abnormally distributed samples that may be representative of compromised samples or hybridisation. From the exposure concentrations selected for detailed analysis, only one such abnormally distributed sample of 2.25 mg/L Cd was identified. This single sample was not included in any subsequent data analyses.
Initial analyses of the normalised data were conducted using pattern recognition analyses SIMCA-P, version 10 (Umetrics, Umeå, Sweden). Principal component analysis (PCA) was used to identify biologically meaningful patterns of transcriptome remodelling. All PCA models generated were cross validated using an iterative procedure in which the model is rebuilt using only 6/7 of the data as a training set and this is then used to predict the class of the remaining 1/7. From this, R 2 X (the fraction of the sum of squares of the X-matrix explained by the model) can be calculated. Typically, a robust PCA model has an R 2 X > 0.5.
Identification of the major pathways involved in toxicant response was restricted to a comparative analysis at an exposure concentration as close as possible to the EC 50 for total brood size for each chemical. Comparison at close to the brood size EC 50 , rather than at higher or lower effect levels was selected since at this exposure level, substantive chemical induced remodeling of the transcriptome of consequence for life-history traits can be expected without apoptotic and necrotic process dominating the signal. Samples selected for the analysis were all successfully hybridized replicates from the Cd 2.25 mg/L (57% reduction in brood size), FA 1 mg/L (37% reduction in brood size) and AZ 150 mg/L (56% reduction in brood size) treatments, as well as all control samples for each test chemical.
To identify the pathways involved in biological response at approximate EC 50 levels, two analytical approaches were used, each utilisng gene ontology as a tool for interpretation of expression change. Firstly, the percentage of genes associated with each GO term changing in expression by 1.5 fold or more in comparison to mean control expression were initially identified for each slide. Genes in this list were then assigned to relevant GO terms and for those terms represented by 10 or more genes, the percentage of 1.5 fold changing genes per term was then compared to the average percentage change. Values of relative change per GO term against mean were then subject to partial least squared discriminant analyses (PLS-DA) using SIMCA-P. All PLS-DA models were cross validated allowing the Q 2 (the fraction of the variation in both the X-and Y-matrices explained by the model), as well as R 2 X and R 2 Y (the fraction of the sum of squares of the Y-matrix explained by the model) values to be calculated (a robust PLS-DA model typically has a Q 2 score > 0.4, while a Q 2 > 0.7 indicates a highly robust model). After confirming separation of control and exposed samples by scores plot analysis, loading plots were reviewed to identify the 20 GO terms most strongly associated with sample separation along the most relevant PC. These GO terms were then assigned to major parental terms for biological process using the GO term Classification Counter http://www.animalgenome.org/bioinfo/tools/ countgo.
For a second analysis of the expression data, control and EC 50 level exposed samples were compared by t-tests to identify significantly changing genes. False discovery correction was not applied in the selection of the set of differential expressed genes because such approaches, while reducing type I errors for null associations can also increases the type II errors for those associations that are not null resulting in the exclusion of truly chemical responsive genes from the gene list [59,60]. Further, the potential for type I error in identifying overrepresented pathways is also controlled for by the fact that the Database for Annotation, Visualization and Integrated Discovery (DAVID) tool that was used for enrichment analysis itself includes an algorithm for multiple sample correction.
The gene list obtained from the t-test comparison was used as input into the DAVID bioinformatic resource available at http://david.abcc.ncifcrf.gov [61]. Analysis was conducted using the functional annotation clustering tool. This tool identifies annotation categories including GO terms, protein-protein interactions, protein functional domains, disease associations, bio-pathways, sequence general features, homologies, gene functional summaries, gene tissue expressions and literatures that are significantly enriched within the gene list [n.b. following multiple sample correction according to 20]. An algorithm is then used to measure the relationships between terms based on their co-associated genes, and from there to group similar, redundant and heterogeneous annotations into discrete clusters. The obtained annotation clusters can then be ranked according to the statistical significance of cluster enrichment. While pathway analysis can provide initial indications of possible effects, scrutiny of individual genes is also a useful approach that can compliment pathway analysis. Thus the expression pattern of individual annotated genes that were ascribed to significant clusters were analysed to highlight transcripts showing patterns of expression that were significantly related to exposure concentrations within a linear regression analysis. To allow open access for the community, all microarray collected was fully MIAMI compliant and has been submitted with appropriate meta-data to Gene Expression Omnibus in a MIAME complient format (as accessions GSE21008, GSE21010, GSE210011 respectively). Submission in this format allows unrestricted (open access) for future deep mining of the data.