Physiological and Molecular Effect Assessment versus Physico

Jul 25, 2011 - Edward J. Perkins , Philipp Antczak , Lyle Burgoon , Francesco Falciani , Natàlia Garcia-Reyero , Steve Gutsell , Geoff Hodges , Aude ...
0 downloads 0 Views 2MB Size
Download PDF
ARTICLE pubs.acs.org/est

Physiological and Molecular Effect Assessment versus Physico-Chemistry Based Mode of Action Schemes: Daphnia magna Exposed to Narcotics and Polar Narcotics Nathalie Dom,*,† Lucia Vergauwen,† Tine Vandenbrouck,† Mieke Jansen,‡ Ronny Blust,† and Dries Knapen†,§ †

Laboratory for Ecophysiology, Biochemistry and Toxicology, Department of Biology, University of Antwerp, Groenenborgerlaan 171, 2020 Antwerp, Belgium ‡ Laboratory of Aquatic Ecology and Evolutionary Biology, Katholieke Universiteit Leuven, Ch. Deberiotstraat 32, 3000 Leuven, Belgium § Physiology and Biochemistry of Domestic Animals, Department of Veterinary Sciences, University of Antwerp, Universiteitsplein 1, 2610 Wilrijk, Belgium

bS Supporting Information ABSTRACT: Structural analogues are assumed to elicit toxicity via similar predominant modes of action (MOAs). Currently, MOA categorization of chemicals in environmental risk assessment is mainly based on the physicochemical properties of potential toxicants. It is often not known whether such classification schemes are also supported by mechanistic biological data. In this study, the toxic effects of two groups of structural analogues (alcohols and anilines) with predefined MOA (narcotics and polar narcotics) were investigated at different levels of biological organization (gene transcription, energy reserves, and growth). Chemical similarity was not indicative of a comparable degree of toxicity and a similar biological response. Categorization of the test chemicals based on the different biological responses (growth, energy use, and gene transcription) did not result in a classification of the predefined narcotics versus the predefined polar narcotics. Moreover, gene transcription based clustering profiles were indicative of the observed effects at higher level of biological organization. Furthermore, a small set of classifier genes could be identified that was discriminative for the clustering pattern. These classifier genes covaried with the organismal and physiological responses. Compared to the physico-chemistry based MOA classification, integrated biological multilevel effect assessment can provide the necessary MOA information that is crucial in high-quality environmental risk assessment. Our findings support the view that transcriptomics tools hold considerable promise to be used in biological response based mechanistic profiling of potential (eco)toxicants.

1. INTRODUCTION Environmental risk assessment (ERA) is mainly based on standardized (acute) toxicity tests on individual species from different trophic levels (e.g., algae, daphnids, fish) in which effect concentrations are determined. In case of data gaps (when no or little effect data are available), computational prediction and modeling techniques such as quantitative structure activity relationships (QSARs), read across, and category formation are applied to predict acute effect concentrations of potential toxicants.1 These empirical data sets and/or the results of the predictive models can only offer limited toxicity information on defined end points without providing real mechanistic information. However, mechanistic information is considered to be of great importance in ecotoxicological studies and in ERA. This information is not only useful to fundamentally understand the impact of pollutants on living organisms, it can also be used to deduce general principles for the categorization and assessment of potential toxicants. In other words, chemicals can also be r 2011 American Chemical Society

grouped according to their biological mode of action (MOA). Combining chemicals into groups with defined properties will improve read across possibilities and priority setting. Furthermore, additional testing of chemicals may become redundant or may be performed in a more targeted way.2,3 Structural analogues are assumed to act via (a) comparable MOA(s). Accordingly, different methods have been developed that categorize chemicals into MOA categories based on their intrinsic physicochemical characteristics.4 6 In general, two main MOA categories are distinguished: baseline toxicants (narcotic compounds) and excess toxicants (reactive and specifically acting Special Issue: Ecogenomics: Environmental Received: April 8, 2011 Accepted: July 25, 2011 Revised: June 22, 2011 Published: July 25, 2011 10

dx.doi.org/10.1021/es201095r | Environ. Sci. Technol. 2012, 46, 10–18

Environmental Science & Technology

ARTICLE

compounds). Baseline toxicants are chemicals that do not interact with a specific receptor in the organism. They can be divided into two subgroups: non polar narcotics and polar narcotic compounds. Toxicity is assumed to be caused by perturbation of the structure and functioning of biological membranes. The toxic potency of these chemicals correlates strongly with their hydrophobicity (log Kow, the octanol/water partitioning coefficient), as a measure for membrane damage. Polar narcotics are assumed to be slightly more toxic than non polar narcotic compounds due to the hydrogen bond capacity of the functional groups. Alternatively, excess toxicants are considered to evoke toxicity via additional, often more specific, mechanisms such as the presence or absence of certain specific chemical functionalities.5,7 Although narcosis is the least specific and nonreactive MOA, it is still particularly important in ecotoxicology since approximately 70% of all industrial organic chemicals are considered to be narcotic substances. Thus far, MOA determination and categorization of these compounds is predominantly based on the structure and resulting physical and chemical properties of potential toxicants. The exact toxic mechanism is still an area of active research.1 The physicochemical properties (e.g., certain structural alerts, high Log Kow, etc.) of a substance are key elements in triggering toxicity, but molecular and physiological events are just as important. It often remains an open question whether this chemically oriented MOA categorization is also supported by biological and mechanistic effect assessment at different levels of biological organization. Hence, there is a need for comprehensive MOA studies at different levels of biological organization in order to fully comprehend how organisms respond to chemical stress.2,8 Information on the molecular mechanisms and pathways leading to a toxic effect can be obtained via transcriptomics analyses. Studying gene transcription alterations assists in the understanding of the biological mechanistic basis of an adverse effect. Simultaneously, transcript fingerprinting offers an efficient means to group chemicals based on discriminative molecular signatures, without prior knowledge of gene annotations and functions.9,10 Chemical effect assessment at the molecular level is considered key to gain insight into the MOA. However, there are many successive and interrelated steps between a molecular initiating event and an adverse outcome at a higher level of biological organization relevant to risk assessment. Therefore, linking effects at the molecular level with higher effect levels is considered to be crucial in the elucidation of a particular toxicant’s MOA.11,12 This study was designed to contribute to the investigation of the underlying molecular and biochemical mechanisms of baseline toxicants in the model organism Daphnia magna. Two sets of structural analogues, defined as narcotics (alcohols: ethanol, methanol, and isopropanol) and polar narcotics (chlorinated anilines: aniline, 4-chloroaniline, 3,5-dichloroaniline, and 2,3,4-trichloroaniline) were selected for short-term toxicity experiments. Effects were measured at different levels of biological organization (gene transcription, energy budget, growth, and immobilization). The transcriptomics data was used both for transcript fingerprinting and for investigating the MOA(s). The MOA classification which is mainly based on chemical structure was biologically challenged: the discriminative potential of transcript fingerprints to distinguish between the predefined narcotics (alcohols) and polar narcotics (anilines) was assessed. In addition, the transcript fingerprints were used

to investigate the perturbed molecular pathways and processes. The (sub)organismal (growth and energy reserves) responses were measured to obtain an integrated overview of the physiological impact of exposure to two groups of structural analogues.

2. MATERIALS AND METHODS 2.1. Chemicals. Details on the purity and suppliers of the test chemicals can be found in the experimental section of the Supporting Information. 2.2. Daphnia magna. Cultures of single clone Daphnia magna were kept in 1-L glass beakers filled with aerated and biofiltered tap water, each containing 20 25 specimens (as described in ref 8). During all exposures, fed (same algae mixture as during culturing) Daphnia neonates were exposed for 48 h and 96 h per exposure condition at a density of 10 organism/50 mL in OECD standard test water of medium hardness (OECD guideline 203, annex 2). Test media were renewed every 24 h for the alcohols (to account for their high(er) volatility) and every 48 h for the anilines. 2.3. Preliminary Experiments. Preliminary immobilization experiments (detailed protocol in experimental section of Supporting Information), in which ten neonates per condition were exposed for 96 h, were conducted to determine every chemical’s specific EC1 96 h and EC10 96 h (US EPA Probit software (http://www.epa.gov/eerd/stat2.htm) (results of these calculations are reported in Table S1). 2.4. Effect Assessment at Different Levels of Biological Organization. Subsequently, daphnia neonates (between 0 and 24 h old) were exposed in glass recipients to each chemical’s respective EC1 96 h and EC10 96 h values (see Table S1) to allow multilevel effect assessment (growth, energy reserves, and gene transcription were assessed) of the seven chemicals at the same level of acute toxicity (exposure at equitoxic acute concentration). The experiments were performed in triplicate at a density of 1 daphnid/5 mL of reconstituted fresh water. The necessary numbers of organisms for the respective end points (growth and energy use) were collected after 48 h and 96 h of exposure and for gene transcription analyses after 96 h of exposures: per replica, 10 daphnids were used for growth, 50 for energy measurements (a pool 25 daphnids was used for lipid measurements, another pool of 25 daphnids was used for glycogen and protein measurements), and 30 for gene transcription analyses. Growth was monitored over a 4-d period. Carapax length measurements (the length from the top of the daphnid’s head to the base of the spine (n = 10)) were performed by means of a microprojector (Projectina, Switzerland). The available energy reserves (lipids, carbohydrates, proteins) were measured on pooled body homogenates (n = 25) collected per condition after 48 h and 96 h of exposure.13,14 Energy fractions were determined spectrophotometrically and transformed into energetic equivalents using their respective energy combustion (39.5 kJ/g lipid, 24 kJ/g protein, 17 kJ/g carbohydrates).13 Lipids were extracted based on the chloroform methanol extraction protocol according to the Bligh and Dyer 15 and absorption was measured at 375 nm. A tripalmitine (Acros Organics, Belgium) standard curve in chloroform was constructed and used for total lipid content calculations. 11

dx.doi.org/10.1021/es201095r |Environ. Sci. Technol. 2012, 46, 10–18

Environmental Science & Technology

ARTICLE

Protein content was measured using the Bradford method.16 A calibration curve of bovine serum albumin (Sigma, Belgium) was used. Concentrations were measured at 600 nm. On the same body homogenates, the total carbohydrate content was determined using Anthron reagents17 and measured at 630 nm. To calculate the carbohydrate content, a glycogen (Sigma, Belgium) calibration curve was made. After testing for normality, one way analyses of variance (ANOVA) were performed, followed by a Bonferroni post hoc test when significant differences (p < 0.05) were observed (SigmaPlot v11.0, Systat Software Inc.). Gene expression analysis was performed on RNA of pooled samples (30 daphnids per replicate (n = 3)) that was extracted using the Trizol RNA extraction procedure (Invitrogen, Belgium), followed by DNase treatment (Fermentas, Germany) and subsequent phenol/chloroform extraction. Aminoallyl dUTP labeled cDNA was produced starting from 5 μg of total RNA from each replicate. A total amount of 50 pmol labeled cDNA with a frequency of incorporation between 20 and 50 was vacuum-dried to use as probe (detailed protocol described in ref 14). The Daphnia magna microarray platform used in this study is an extension of a previously developed custom array that consisted of 2445 life stage, molting, and energy metabolism related gene fragments selected through Suppression Subtractive Hybridization (SSH) procedures.18,19 Sequencing enabled the selection of the unique gene fragments (n = 1189), which were supplemented with another 734 gene fragments related to other stress factors (fish, Pasteuria ramosa, and carbaryl stress) selected through SSH experiments.20 This new custom Daphnia magna cDNA microarray platform is composed of 1734 gene fragments (NCBI GEO platform GPL13463). The cDNA fragments were spotted in triplicate on aminosilane coated slides (Asper Biotech, Estonia) using the Q-Array Mini (Genetix, UK). Hybridization, washing, and scanning of the microarray slides were performed according to the protocol described in ref 14. The hybridization design was a universal reference design (a mixture of aliquots from control and exposed samples), which is recommended when class discovery is the main purpose of the experiment.21 One of the three biological replicates of each exposure condition was labeled with one dye, the remaining two samples were labeled with the second dye.14 Statistical analyses were performed using the R package limma (linear models for microarray data).22 Spots for which red and green median foreground intensity < median background intensity +2 SD on all arrays were deleted before analysis. Median intensity data were background corrected using a normal-exponential convolution model.23 Subsequently, normalization between-arrays was performed using Variance Stabilization Normalization (vsn).24 For each probe, a linear model was fitted to intensity ratios, after which probes were ranked in order of evidence of differential expression using an empirical Bayes method. Both high and low concentrations of each chemical stressor were contrasted against their respective controls. Genes with p < 0.05 and log2FC < 0.75 or log2FC > 0.75 (log2 fold change) were considered as differentially transcribed gene fragments. MIAME-compliant25 raw microarray data have been deposited with NCBI’s Gene Expression Omnibus26 and are accessible through GEO series accession number GSE29993. Clustering analyses were performed using MultiExperiment Viewer (MeV, v4.5.1, TIGR).27 Genes that were differentially expressed in one of the 14 conditions (EC1 and EC10 treatment of every tested chemical) were all included in the analyses. The

corresponding log2-ratios of the other conditions were included in the analyses as well (n = 392). To investigate the degree of similarity between the multivariate gene transcription data sets of the different chemical treatments, hierarchical clustering (HCL) analysis (Euclidean distance, average linkage) was performed with statistical support (bootstrapping; n = 1000)) for the nodes of the trees. Significance Analysis of Microarray (SAM) was used to identify key genes that allow distinguishing among compounds with a different (molecular) MOA. Furthermore, a multivariate analyses, partial least-squares (PLS) was performed using the Simca P11.5 software package (Umetrics) to evaluate the correspondence between the transcriptomics data (x-variables) and the higher level effects (lengths and energy reserves; y-variables). R2 and Q2 values are reported: R2 is the fraction of variation explained by the model. Q2 indicates the fraction of total variation of the x-variables that can be predicted by a component, as estimated by cross validation of R2 values closer to 1 indicate a better model. The genes that contributed the most to the model can be selected based on a Variable ImPortance (VIP) score. Genes (x-variables) with a VIP > 1 are assumed to be the most relevant for explaining the y-variables. In addition, the sequence information of the key genes (identified by the SAM methods and the PLS) was submitted to Blast2GO28 for biological gene function assessment of the differentially transcribed genes.

3. RESULTS AND DISCUSSION In this study, the effects of two groups of compounds with predefined MOA (by means of physico-chemistry based MOA classification) were investigated at different levels of biological organization. Alcohols are defined as non polar narcotic compounds, whereas chlorinated anilines are defined as polar narcotic compounds.7 Because it is generally assumed that structural analogues elicit toxicity via a similar MOA, one could expect that the multilevel effects would be indicative for each predefined MOA group. However, the results of this study did not support this assumption. 3.1. Growth and Energy Reserves. Under normal conditions, individuals use energy for their metabolism, growth, and reproduction. When facing additional costs such as preventing or restoring damage from pollutant stress, organisms need to relocate some of the acquired energy to repair and/or maintain their physiological integrity, leaving less energy for growth and eventually for reproduction.13 In addition, since chemical exposure can induce a decreased food intake, the organism may also be less efficient in acquiring energy. Effects on the energy stores differed among the different chemical exposures (Figure 1). Exposures of 48 h and 96 h to AN, CA, DCA, and isopropanol (further named group A compounds in this manuscript) did not elicit significant effects on the total caloric content or on the growth of the exposed animals compared to their respective controls (p > 0.05). Thus, normal homeostasis of the animals under toxicant stress did not seem to be disturbed. In contrast, the compensatory changes in energy supply after exposure to TCA, ethanol, and methanol (further named group B compounds; Figure 1A, C, E) were clearly at the expense of growth since significant differences in length were already observed after 48 h of exposure (Figure 1B, D, F). Exposure to TCA EC1 significantly decreased the total lipid content (p = 0.035) and the total protein content (p = 0.004) but the amount 12

dx.doi.org/10.1021/es201095r |Environ. Sci. Technol. 2012, 46, 10–18

Environmental Science & Technology

ARTICLE

Figure 1. Effects on total caloric content (expressed as a percentage of the respective control) after 96 h and on growth of Daphnia magna during a 96-h exposure period to TCA (A, B), ethanol (C, D) and methanol (E, F). Significant differences are indicated by (*) if p < 0.05 and (**) if p < 0.001.

of carbohydrates was not significantly affected (p = 0.168) compared to the control. Ethanol EC10 exposure caused a significant decrease in the total caloric content. Although the total lipid content was significantly decreased (p = 0.016) after ethanol EC1 exposure, this had no significant impact on the total caloric content (p > 0.05). 3.2. Differential Gene Transcription. 3.2.1. Chemical Similarity versus Gene Transcript Fingerprinting. A total of 392 gene fragments was identified as differentially transcribed in at least one of the chemical treatments. On this set of 392 gene fragments,

HCL was applied to assess whether it was possible to discriminate between nonpolar and polar narcotic compounds based on their molecular fingerprints (Figure 2). For every compound except for ethanol, the low (EC1) and the higher (EC10) exposures clustered together per compound. Not the exposure concentration but the chemical compound was the main discriminator in our experiments. The high bootstrap values (>70% in every case) indicate the high reliability of the results. The dendrogram of the HCL analysis clearly divided the chemicals into two separate clusters. With the exception of the low exposure concentration of ethanol (EC1), 13

dx.doi.org/10.1021/es201095r |Environ. Sci. Technol. 2012, 46, 10–18

Environmental Science & Technology

ARTICLE

Figure 2. Hierarchical clustering (Euclidean distance, average linkage) based on log2-ratios of the differentially transcribed genes in Daphnia magna after 96 h exposure to the different test chemicals. When a gene was differentially expressed in one of the 14 conditions the corresponding ratios of the other conditions were included in the analyses (n = 392). Bootstrap values (in %) are indicated on the branches.

these two clusters were identical to the two groups (group A and group B) of compounds that were identified based on the results of growth and energy assessment. These group B compounds (TCA, ethanol, and methanol) triggered a higher amount of differentially transcribed genes (174, 145, and 179 genes respectively, combination of EC1 and EC10) in comparison to aniline, CA, DCA, and isopropanol exposure (11, 50, 30, and 14 differentially transcribed genes respectively). According to Verhaar et al.,29 a narcotic MOA can be assumed for simple organics such as aliphatic alcohols or aromatic hydrocarbons. A generalization of this approach (e.g., by the implementation in the Toxtree software and OECD QSAR Toolbox) has led to a simple structural rule-based compound MOA classification. Chemicals that elicit a narcotic level of toxicity were identified as compounds containing any combination of C, H, O, and halogen, excluding R,β-unsaturated carbonyl groups such as electrophilic functionalities. Chemicals having additional hydrogen bond donor acidity are considered to be polar narcotic compounds (e.g., phenols or anilines). However, this chemical structure based approach was not consistent with the results obtained in this study. No clear discrimination between the predefined narcotics and polar narcotics could be made based on gene fingerprints or based on higher effect levels (growth and energy use). However, the Verhaar classification scheme (which is now generally used in aquatic ecotoxicology) was initially based on fish data obtained under acute exposure regimes.29 In this context, we already illustrated in a previous study that, while (chlorinated) anilines act as polar narcotics in fish (and also algae and bacteria), the toxicity responses in D. magna were clearly different.8 This deviating acute toxicity pattern might be the underlying reason why a clear clustering of two groups of structural analogues was not observed. Furthermore, in a study of Russom et al.,4 which not only addressed the chemical class but also the direct toxicity of a large set of chemicals, TCA was categorized as a nonpolar narcotic. TCA was found to be additive to both a reference nonpolar narcotic chemical (octanol) and a reference polar narcotic chemical (phenol). Chemicals with a relatively high log Kow (> 2.7) that would normally be assigned to the polar narcotic class solely based on their chemical structure, tended to behave more like narcotic compounds.30 Thus for compounds such as TCA with a log Kow > 2.7, the polar narcotic effect of hydrogen bonding was moderate for these compounds. Consequently, it is not that surprising the predefined narcotics and polar narcotics group differently in this mechanistic study. Obviously, these divergent MOA categorizations (i.e., MOA based on chemical similarity versus MOA based on biological

response data) can have certain implications on ERA procedures. As already stated in the introduction, QSARs are regarded as important and recommended tools in the reduction of animal use in ERA. Baseline toxicity QSARs are considered to predict the aquatic toxicity of either narcotic or polar narcotic compounds with reasonable accuracy.31 However, a wrong MOA classification will ultimately lead to low quality toxicity predictions. Therefore, it has already been suggested that the inclusion of biological data will be essential to improve compound MOA categorization schemes.32,33 3.2.2. Discriminative Gene Set. In a next step, we tried to identify the crucial genes responsible for that particular grouping (group A versus group B). SAM identified 45 unique indicator or marker genes with a FDR < 0.05 (Figure S1 and Table S2 of Supporting Information). A general tendency of down regulation (95% of the indicator genes) was observed after exposure to group B chemicals (chemicals that induced growth impairment and a significant decrease in caloric content). In contrast, a tendency of up regulation of these marker genes was observed in group A chemicals (chemicals that did not induce growth impairment nor a decrease in energy stocks). The discriminative power of these 45 genes was illustrated by the HCL results: the dendogram of the HCL based on these marker genes (Figure S1 and Table S2) discriminated the chemicals in a manner comparable to the HCL based on all the differentially expressed gene transcripts (n = 392) (Figure 2). These results illustrate the potential of using transcriptomics tools for future compound categorization purposes. Such indicator or marker genes can in fact be used in a classifier to biologically screen and discriminate chemicals based on similarities in their response pattern: the specific gene transcription patterns of those indicator genes per specific molecular MOA can form the basis of a classifier tool. At this point, the 45 identified marker genes are only able to distinguish between the two tested compound groups that elicit a different biological response, but such gene sets can easily be extended and/or specified for other groups. For actual mechanistic and predictive applications in ERA, a broader reference framework comprising additional baseline (and excess) toxicity compounds will be essential. After selection of specific indicator genes of a classifier, the next step would be the evaluation of these genes in a validation concept. Simultaneously, it will be important to select the most appropriate gene set (as highly correlated genes do not necessarily increase the prediction outcome), leading to a minimum amount of descriptive marker genes that will allow for high throughput screening, without compromising the classification power. As 14

dx.doi.org/10.1021/es201095r |Environ. Sci. Technol. 2012, 46, 10–18

Environmental Science & Technology

ARTICLE

certain physicochemical descriptors are used to determine the activity of various chemicals, molecular biological descriptors (e.g., up or down regulation of a particular set of genes) could also be used in the future. Such gene activity does not necessarily have to be determined using expensive and time-consuming microarrays. For differential gene expression analysis of a relatively small set of key genes, other techniques, such as quantitative real time PCR (qRT-PCR), are faster and more cost-effective alternatives. 3.2.3. Discriminative Indicator Genes and Their Biological Function. The function of the indicator genes and whether or not they are (causally) linked to the toxic effects at higher level of biological organization is rather subordinate for compound classification purposes based on gene fingerprinting. Hence, gene function is less important if the main purpose of the classifier genes is to be used in chemical grouping and categorization. The fact that even unannotated cDNA libraries and microarrays can be used for this purpose is a significant advantage.34 However, interpretation of the biological functions of these genes can be very valuable because biological interpretation can assist in gaining biological insight into the specific differences in the (toxic) MOAs among different compounds. Therefore, the sequence information of the 45 marker genes was submitted to Blast2Go.28 Only for 28 sequences of these 45 key genes, GO information could be retrieved: 10 sequences could not be annotated and 7 gene transcripts were described as “hypothetical proteins” of which no further GO annotations were given (Table S2). Complementary to the (sub)organismal effects, it is noteworthy that most of the annotated marker genes were involved in processes essential for growth and energy supply (e.g., cuticle formation, transcription/translation, metabolic processes such as glycolysis, gluconeogensis, lipolysis, and the citric acid cycle, transport and binding of essential biomolecules such as ATP, lipids, and proteins). First, indicator genes encoding for cuticula and ribosomal proteins were identified. Because of their central role in the synthesis of a new exoskeleton, cuticula proteins are essential for daphnia growth. The differences in their gene transcription levels can either be a direct consequence of the exposure or can be related to the actual molting behavior of the daphnids.35 Ribosomes also play a major role in growth of the organisms, as they are key actors in protein synthesis. An induction of genes coding for these proteins may represent a general attempt to overcome toxicant stress and/or to restore the damage elicited by the chemical exposure.36 Another important group of identified marker genes was involved in basic energy metabolism related pathways (such as carbohydrate, lipid, and protein metabolism). Carbohydrates are crucial in any organism’s energy metabolism. They fuel the energy-demanding maintenance and synthesis processes but also provide the necessary building blocks for various anabolic pathways.13 One of the indicator genes was a gene coding for the bifunctional enzyme, 6-phophofructo-2-kinase/fructose-2,6biphosphatase that plays a unique role in the control of glucose homeostasis by allowing to switch from glycolysis to gluconeogenesis.37 Proteins and lipids are mobilized as well to supply sufficient energy in order to overcome toxicant exposure and to maintain homeostasis during toxicant exposure.13 The gene encoding for cathepsine L, a lysosomal cysteine proteinase that plays a major role in intracellular protein catabolism, was up regulated after exposure to group A compounds, possibly in an attempt to maintain the basal metabolism during exposure.

Furthermore, genes encoding for fatty acid binding proteins (see Table S2 for Genbank accession numbers) were identified as indicator genes. They are involved in transporting free fatty acids to the mitochondria for degradation into acetyl-CoA. Acetyl CoA functions as input to the citric acid. Moreover, it is also a precursor for ketone bodies that are produced under stress conditions such as chemical exposure. Additionally, gene transcripts involved in the electron transport system and oxidative phosphorylation were identified as marker genes, e.g., a gene coding for cytochrome c, an electroncarrying mitochondrial protein. Its main function in the cellular respiration is to transport electrons from cytochrome c reductase (Complex III) to cytochrome oxidase (Complex IV). Another important enzyme in the cellular respiration is the mitochondrial transporter, ATP/ADP translocase. The gene encoding for this enzyme was also up regulated after group A exposure. ATP/ADP translocase assists in moving ATP out of the mitochondrial inner membrane, while assisting ADP to enter the inner mitochondrial membrane as a precursor to synthesize ATP. Finally, four vitellogenin (Vtg) fragments (Vtg 1 and 2, Vtg fused with superoxide dismutase, and Vtg family member vit-2) were identified as marker genes. They were up regulated after exposure to group A compounds and down regulated after exposure to group B compounds. In previous studies, down regulation of Vtg related genes has been suggested to be an early marker of reproduction impairment (through suppressed energy supply to the developing embryo), whereas the up regulation of these genes was suggested to be indicative of enhanced energy demands to cope with toxicant stress and consequent mobilization of lipid reserves.36 A decrease in the organism’s growth due to energy depletion will ultimately impact survival and reproduction at the population level. Especially in case of crustacean reproduction, fecundity is directly related to the daphnid’s size since a large size allows rapid maturation and consequently shortens the time to first reproduction.13,38 In general, transcripts involved in the global energy metabolism were up regulated after exposure to group A chemicals (AN, CA, DCA, and isopropanol), and exposure to these chemicals did not cause the organisms’ energy reserves to be depleted. On the other hand, energy metabolism related transcripts were down regulated after exposure to group B chemicals (TCA, ethanol, and methanol), which had a clear negative effect on the energy budget and the growth of the organisms. This may sound contradictory at first, as high transcript levels of these genes may be thought to result in increased protein levels and increased activity of the involved pathways, thereby more efficiently oxidizing fuel and thus more quickly depleting the organism’s energy reserves. However, there is not always a direct link between transcript and protein levels. Gene transcription does not only directly influence protein levels but in turn can also be influenced by the presence of certain proteins.39 For example, transcript levels may be increased simply to maintain protein concentrations when protein degradation is increased or when protein synthesis is impaired. Such maintenance efforts may or may not be successful, so even decreased protein levels may be observed in conjunction with increased transcript levels. More importantly, even when the link between transcript and protein levels is known, the regulation of the energy metabolism involves many more complex processes and the overall energy status is, of course, dependent on energy intake as well. It is possible, for example, that the up regulated transcripts in group A allow compensatory processes to proceed more efficiently, resulting in 15

dx.doi.org/10.1021/es201095r |Environ. Sci. Technol. 2012, 46, 10–18

Environmental Science & Technology

ARTICLE

Figure 3. Partial least-squares (PLS) analysis of the different treatments using differentially expressed genes as predictors and (sub)organismal effects as response data. R2 and Q2 for this two-compartment model are, respectively, 0.892 and 0.711. The first component explained 30.0% of the variability, the second component explained 19.8%. The genes with a Variable Importance (VIP) score above 1.5 for both components in the model are indicated in blue and in red. The genes in blue are numbered according to their appearance in Figure 3.

less toxic effects and lower strain levels on the energy metabolism in general. Therefore, transcriptional effects should be interpreted with great care. They can be very valuable to point out the major metabolic and biochemical routes that are affected by chemical exposure, and certainly allow formulation of hypotheses to be tested in new experiments, but a simple and linear causeand-effect logic is very hard to deduce in most cases. 3.3. Use of Integrated Multilevel Responses in ERA. It is therefore clear that the linkage of molecular data with phenotypic responses of ecotoxicological importance remains one of the major difficulties in the interpretation and evaluation of microarray data. In the previous paragraphs, the effects of exposure were in fact discussed per biological level: effects were described separately on the organismal (growth), biochemical/physiological (energy reserves), and molecular (gene transcription) levels. In an effort to integrate these multilevel effects a PLS analysis was performed. Genes with a VIP score higher than 1 contribute the most to the model and are assumed to be the most relevant for explaining the y-variables. However, using this cutoff (of VIP > 1) resulted in a rather high amount of genes (n = 159). All 45 previously identified indicator genes (by the SAM analysis) were an integral part of this list. To fine-tune the selection of genes and to only focus on the most important ones, the cutoff was set to VIP > 1.5. This resulted in 47 genes. Most of these predictor genes were plotted in close proximity of the organismal and physiological response data, indicating covariance of the genes with the higher effect levels (see Figure 3). Two of these genes had an inverse relationship with the higher effect levels. From the 47 genes, 26 genes were previously identified as key genes in the SAM analysis (tagged in Figure 3 according to the numbering in Figure S1 and Table S2). The remaining genes were identified as unknowns, hypothetical proteins, and Vtg(-like) genes.

This analysis illustrated that the majority of the discriminative indicator genes responsible for the clustering of group A versus group B compounds were correlated, and therefore possibly linked to the higher level effects (growth and energy reserves). Although a causal relationship still remains to be demonstrated by testing specific hypotheses that result from this study, it is clear that the identified indicator genes can have potential to be used as molecular biomarkers. In agreement with previous studies,40,41 our findings support the view that transcriptomics tools hold considerable promise to be used in biological response based mechanistic profiling of potential (eco)toxicants. In fact, transcriptomics data can provide the missing biological information that is needed to make a reliable biology-based MOA class categorization.

’ ASSOCIATED CONTENT

bS

Supporting Information. Additional information on the preliminary experiments together with an overview of every chemical’s respective EC1 96 h and EC10 96 h (Table S1) that allowed effect assessment at the different levels of biological organization at equitoxic concentration; Figure S1, the results of the clustering based on the 45 key genes (identified by the SAM analyses) and a heatmap in which the functions of these key genes are represented; Table S2 summarizes the GO annotation results of these respective 45 key genes.This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*Tel: +32 (0)3 265 35 32; fax: +32 (0)3 265 34 97; e-mail: [email protected]. 16

dx.doi.org/10.1021/es201095r |Environ. Sci. Technol. 2012, 46, 10–18

Environmental Science & Technology

ARTICLE

’ ACKNOWLEDGMENT We thank Karin Van den Bergh for technical assistance. This research was funded by a PhD grant of the Agency for Innovation by Science and Technology in Belgium (IWT, project 73275) and by the OSIRIS - EU project (GOCE-CT-2007-037017).

(17) Roe, J. H.; Dailey, R. E. Determination of glycogen with Anthron reagent. Anal. Biochem. 1966, 15 (2), 245–&. (18) Soetaert, A.; van der Ven, K.; Moens, L. N.; Vandenbrouck, T.; van Remortel, P.; De Coen, W. M. Daphnia magna and ecotoxicogenomics: Gene expression profiles of the anti-ecdysteroidal fungicide fenarimol using energy-, molting- and life stage-related cDNA libraries. Chemosphere 2007, 67 (1), 60–71. (19) Soetaert, A.; Moens, L. N.; Van der Ven, K.; Van Leemput, K.; Naudts, B.; Blust, R.; De Coen, W. M. Molecular impact of propiconazole on Daphnia magna using a reproduction-related cDNA array. Comp. Biochem. Physiol., C: Toxicol. Pharmacol. 2006, 142 (1 2), 66–76. (20) Orsini, L.; Jansen, M.; Souche, E.; Geldof, S.; De Meester, L. EST characterization and SNPs discovery in the waterflea Daphnia magna. BMC Genomics 2011, 12, 309. (21) Knapen, D.; Vergauwen, L.; Laukens, K.; Blust, R. Best practices for hybridization design in two-colour microarray analysis. Trends Biotechnol. 2009, 27 (7), 406–414. (22) Wettenhall, J. M.; Smyth, G. K. limmaGUI: A graphical user interface for linear modeling of microarray data. Bioinformatics 2004, 20 (18), 3705–3706. (23) Ritchie, M. E.; Silver, J.; Oshlack, A.; Holmes, M.; Diyagama, D.; Holloway, A.; Smyth, G. K. A comparison of background correction methods for two-colour microarrays. Bioinformatics 2007, 23 (20), 2700–2707. (24) Huber, W.; von Heydenbreck, A.; Sultmann, H.; Poustka, A.; Vingron, M. Variance stabilization applied to microarray data calibration and to the quantification of differential expression. Bioinformatics 2002, 18, S96–104. (25) Brazma, A. Minimum Information About a Microarray Experiment (MIAME) - Successes, Failures, Challenges. Thescientificworld 2009, 9, 420–423. (26) Barrett, T.; Edgar, R. Gene expression omnibus: Microarray data storage, submission, retrieval, and analysis. DNA Microarrays, Part B: Databases Stat. 2006, 411, 352–369. (27) Saeed, A. I.; Sharov, V.; White, J.; Li, J.; Liang, W.; Bhagabati, N.; Braisted, J.; Klapa, M.; Currier, T.; Thiagarajan, M.; Sturn, A.; Snuffin, M.; Rezantsev, A.; Popov, D.; Ryltsov, A.; Kostukovich, E.; Borisovsky, I.; Liu, Z.; Vinsavich, A.; Trush, V.; Quackenbush, J. TM4: A free, open-source system for microarray data management and analysis. Biotechniques 2003, 34 (2), 374. (28) Conesa, A.; Gotz, S.; Garcia-Gomez, J. M.; Terol, J.; Talon, M.; Robles, M. Blast2GO: A universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 2005, 21 (18), 3674–3676. (29) Verhaar, H. J. M.; Vanleeuwen, C. J.; Hermens, J. L. M. Classifying environmental pollutants. 1. Structure-Activity Relationships for prediction of aquatic toxicity. Chemosphere 1992, 25 (4), 471–491. (30) Veith, G. D.; Broderius, S. J. Rules for Distinguishing Toxicants That Cause Type-I and Type-II Narcosis Syndromes. Environ. Health Perspect. 1990, 87, 207–211. (31) von der Ohe, P. C.; Kuhne, R.; Ebert, R. U.; Altenburger, R.; Liess, M.; Schuurmann, G. Structural alerts - A new classification model to discriminate excess toxicity from narcotic effect levels of organic compounds in the acute daphnid assay. Chem. Res. Toxicol. 2005, 18 (3), 536–555. (32) Cronin, M. T. D. The current status and future applicability of quantitative structure-activity relationships (QSARs) in predicting toxicity. ATLA, Altern. Lab. Anim. 2002, 30, 81–84. (33) Dom, N.; Nobels, I.; Knapen, D.; Blust, R. Bacterial gene profiling assay applied as an alternative method for mode of action classification: Pilot study using chlorinated anilines. Environ. Toxicol. Chem. 2011, 30 (5), 1059–1068. (34) Tyler, C. T.; Filby, A. L.; Taisen, I.; Kramer, V. J.; Larson, J. D.; van Aggelen, G.; van Leeuwen, K.; Viant, M. R.; Tillitt, D. E. Application of genomics to tiered testing. In Genomics in Regulatory Ecotoxicology; Ankley, G. T., Miracle, A. L., Perkins, E. J., Daston, G. P., Eds.; CRC Press: London, 2008. (35) De Schamphelaere, K. A. C.; Vandenbrouck, T.; Muyssen, B. T. A.; Soetaert, A.; Blust, R.; De Coen, W.; Janssen, C. R. Integration

’ REFERENCES (1) Bradbury, S. P.; Russom, C. L.; Ankley, G. T.; Schultz, T. W.; Walker, J. D. Overview of data and conceptual approaches for derivation of quantitative structure-activity relationships for ecotoxicological effects of organic chemicals. Environ. Toxicol. Chem. 2003, 22 (8), 1789–1798. (2) Escher, B. I.; Hermens, J. L. M. Modes of action in ecotoxicology: Their role in body burdens, species sensitivity, QSARs, and mixture effects. Environ. Sci. Technol. 2002, 36 (20), 4201–4217. (3) Vonk, J. A.; Benigni, R.; Hewitt, M.; Nendza, M.; Segner, H.; van de Meent, D.; Cronin, M. T. D. The Use of Mechanisms and Modes of Toxic Action in Integrated Testing Strategies: The Report and Recommendations of a Workshop held as part of the European Union OSIRIS Integrated Project. ATLA, Altern. Lab. Anim. 2009, 37 (5), 557–571. (4) Russom, C. L.; Bradbury, S. P.; Broderius, S. J.; Hammermeister, D. E.; Drummond, R. A. Predicting modes of toxic action from chemical structure: Acute toxicity in the fathead minnow (Pimephales promelas). Environ. Toxicol. Chem. 1997, 16 (5), 948–967. (5) Verhaar, H. J. M.; Vanleeuwen, C. J.; Hermens, J. L. M. Classifying environmental pollutants. 1. Structure-Activity Relationships for predidiction of aquatic toxicity. Chemosphere 1992, 25 (4), 471–491. (6) Netzeva, T.; Pavan, M.; Worth, A. Review of Data Sources, QSARs and Integrated Testing Strategies for Aquatic Toxicity; EUR 22943; European Commission, Joint Research Centre: Ispra, Italy, 2007. (7) Verhaar, H. J. M.; Ramos, E. U.; Hermens, J. L. M. Classifying environmental pollutants 0.2. Separation of class 1 (baseline toxicity) and class 2 ('polar narcosis') type compounds based on chemical descriptors. J. Chemom. 1996, 10 (2), 149–162. (8) Dom, N.; Knapen, D.; Benoot, D.; Nobels, I.; Blust, R. Aquatic multi-species acute toxicity of (chlorinated) anilines: Experimental versus predicted data. Chemosphere 2010, 81 (2), 177–186. (9) Gatzidou, E. T.; Zira, A. N.; Theocharis, S. E. Toxicogenomics: A pivotal piece in the puzzle of toxicological research. J. Appl. Toxicol. 2007, 27 (4), 302–309. (10) Ankley, G. T.; Daston, G. P.; Degitz, S. J.; Denslow, N. D.; Hoke, R. A.; Kennedy, S. W.; Miracle, A. L.; Perkins, E. J.; Snape, J.; Tillitt, D. E.; Tyler, C. R.; Versteeg, D. Toxicogenomics in regulatory ecotoxicology. Environ. Sci. Technol. 2006, 40 (13), 4055–4065. (11) Schirmer, K.; Fischer, B. B.; Madureira, D. J.; Pillai, S. Transcriptomics in ecotoxicology. Anal. Bioanal. Chem. 2010, 397 (3), 917–923. (12) Ankley, G. T.; Bennett, R. S.; Erickson, R. J.; Hoff, D. J.; Hornung, M. W.; Johnson, R. D.; Mount, D. R.; Nichols, J. W.; Russom, C. L.; Schmieder, P. K.; Serrrano, J. A.; Tietge, J. E.; Villeneuve, D. L. Adverse outcome pathways: A conceptual framework to support ecotoxicology research and risk assessment. Environ. Toxicol. Chem. 2010, 29 (3), 730–741. (13) De Coen, W. M.; Janssen, C. R. The missing biomarker link: Relationships between effects on the cellular energy allocation biomarker of toxicant-stressed Daphnia magna and corresponding population characteristics. Environ. Toxicol. Chem. 2003, 22 (7), 1632–1641. (14) Vandenbrouck, T.; Soetaert, A.; van der Ven, K.; Blust, R.; De Coen, W. Nickel and binary metal mixture responses in Daphnia magna: Molecular fingerprints and (sub)organismal effects. Aquat. Toxicol. 2009, 92 (1), 18–29. (15) Bligh, E. G.; Dyer, W. J. A rapid method of total lipid extraction and purification. Can. J. Biochem. Physiol. 1959, 37 (8), 911–917. (16) Bradford, M. M. Rapid and sensitive method for quantitation of microgram quantities of protein utilizing principle of protein-dye binding. Anal. Biochem. 1976, 72 (1 2), 248–254. 17

dx.doi.org/10.1021/es201095r |Environ. Sci. Technol. 2012, 46, 10–18

Environmental Science & Technology

ARTICLE

of molecular with higher-level effects of dietary zinc exposure in Daphnia magna. Comp. Biochem. Physiol., D: Genomics Proteomics 2008, 3 (4), 307–314. (36) Pereira, J. L.; Hill, C. J.; Sibly, R. M.; Bolshakov, V. N.; Goncalves, F.; Heckmann, L. H.; Callaghan, A. Gene transcription in Daphnia magna: Effects of acute exposure to a carbamate insecticide and an acetanilide herbicide. Aquat. Toxicol. 2010, 97 (3), 268–276. (37) Rider, M. H.; Bertrand, L.; Vertommen, D.; Michels, P. A.; Rousseau, G. G.; Hue, L. 6-Phosphofructo-2-kinase/fructose-2,6-bisphosphatase: Head-to-head with a bifunctional enzyme that controls glycolysis. Biochem. J. 2004, 381, 561–579. (38) Koivisto, S. Is Daphnia magna an ecologically representative zooplankton species in toxicity tests? Environ. Pollut. 1995, 90 (2), 263–267. (39) Brazhnik, P.; de la Fuente, A.; Mendes, P. Gene networks: How to put the function in genomics. Trends Biotechnol. 2002, 20 (11), 467–472. (40) Moens, L. N.; van der Ven, K.; Van Remortel, P.; Del-Favero, J.; De Coen, W. M. Expression profiling of endocrine-disrupting compounds using a customized Cyprinus carpio cDNA microarray. Toxicol. Sci. 2006, 93 (2), 298–310. (41) Nota, B.; Verweij, R. A.; Molenaar, D.; Ylstra, B.; van Straalen, N. M.; Roelofs, D. Gene Expression Analysis Reveals a Gene Set Discriminatory to Different Metals in Soil. Toxicol. Sci. 2010, 115 (1), 34–40.

18

dx.doi.org/10.1021/es201095r |Environ. Sci. Technol. 2012, 46, 10–18