- Research article
- Open Access
'Systems toxicology' approach identifies coordinated metabolic responses to copper in a terrestrial non-model invertebrate, the earthworm Lumbricus rubellus
BMC Biology volume 6, Article number: 25 (2008)
New methods are needed for research into non-model organisms, to monitor the effects of toxic disruption at both the molecular and functional organism level. We exposed earthworms (Lumbricus rubellus Hoffmeister) to sub-lethal levels of copper (10–480 mg/kg soil) for 70 days as a real-world situation, and monitored both molecular (cDNA transcript microarrays and nuclear magnetic resonance-based metabolic profiling: metabolomics) and ecological/functional endpoints (reproduction rate and weight change, which have direct relevance to population-level impacts).
Both of the molecular endpoints, metabolomics and transcriptomics, were highly sensitive, with clear copper-induced differences even at levels below those that caused a reduction in reproductive parameters. The microarray and metabolomic data provided evidence that the copper exposure led to a disruption of energy metabolism: transcripts of enzymes from oxidative phosphorylation were significantly over-represented, and increases in transcripts of carbohydrate metabolising enzymes (maltase-glucoamylase, mannosidase) had corresponding decreases in small-molecule metabolites (glucose, mannose). Treating both enzymes and metabolites as functional cohorts led to clear inferences about changes in energetic metabolism (carbohydrate use and oxidative phosphorylation), which would not have been possible by taking a 'biomarker' approach to data analysis.
Multiple post-genomic techniques can be combined to provide mechanistic information about the toxic effects of chemical contaminants, even for non-model organisms with few additional mechanistic toxicological data. With 70-day no-observed-effect and lowest-observed-effect concentrations (NOEC and LOEC) of 10 and 40 mg kg-1 for metabolomic and microarray profiles, copper is shown to interfere with energy metabolism in an important soil organism at an ecologically and functionally relevant level.
Understanding biological responses to individual toxic chemicals and chemical classes is clearly of key importance for pollution assessment, both for monitoring exposure to existing environmental contamination and for informing the risk assessment of off-target effects. However, ecotoxicological research frequently focuses only on easily measurable endpoints, typically mortality, although more sensitive tests on effect endpoints such as reproduction and growth are also used widely. Thus, a major challenge for ecotoxicology is understanding toxic mechanisms at a molecular level, and how these molecular changes relate to functional changes at the organism and population level . The 'ecotoxicogenomic' post-genomic approach has clear benefits, and is currently generating interest from end users such as regulatory authorities as well as from research scientists [2, 3]. In order for this potential to be realised, a solid bedrock of research is needed to characterise the fundamental responses of important test organisms to a range of model toxins covering a wide chemical space. It will be important to determine just how specific omic fingerprints of toxicity are, and whether they can be used successfully to distinguish between different modes of toxic action, and hence yield novel information on mechanistic toxicology. This 'systems toxicology' approach has been applied in widely used model organisms such as the laboratory rat and other vertebrates [4–7]. However, these animal models have the benefit of many more existing data [8, 9]. In addition, it is often easier to perform manipulative experiments, and there is a much greater scope for complementary mechanistic cell-based work, such as histopathology. In contrast, the situation with non-model, ecologically relevant species is quite different.
The term 'ecologically relevant' is not precisely defined: clearly the most relevant level for studying the effects of chemicals is the community and/or ecosystem, and there are approaches which aim to understand, or at least quantify, responses to pollution at this level (see, for example, [10–14]). Here, however, we refer to controlled studies on single species that may already be widely studied but are not classic model organisms; for example, animals used in regulatory ecotoxicity tests fall into this category, such as the earthworm Eisenia fetida, the enchytraeid Enchytraeus albidus, and collembolans Folsomia candida and Orchesella cincta for terrestrial, and Daphnia magna, Gammarus pulex, chironomid larvae and Mytilus species for aquatic testing. Working with these animals presents some common challenges: none has a fully sequenced genome; it is not generally possible to obtain antibodies against specific molecular targets; they are often so small as to preclude ready dissection of internal organs or tissues; it is impossible or extremely difficult to modulate gene activity, for example by creating knockout strains; and there is in general much less knowledge about fundamental biological systems, such as signalling pathways or gene regulation, in these organisms. Modern omic approaches offer a potential opportunity to circumvent some of these drawbacks [15–22]. In particular, metabolomics and metabonomics have one great advantage for work with non-model organisms: because metabolites are detected directly, and primary metabolites at least are identical across different species, samples can trivially be analysed with no need for prior knowledge of the gene and protein sequences . Metabolomics also reports on the final integrated phenotype of an organism, as metabolism is the final downstream product of gene and enzyme regulation [24–27]. As a consequence, we decided to carry out an integrative study of the metabolic response of Lumbricus rubellus to copper, using both nuclear magnetic resonance (NMR)-based metabolic profiling and cDNA microarrays for transcript profiling.
The earthworm L. rubellus is a common species with a worldwide distribution [28–31]. It is found even in shallow and contaminated soils, so it is appropriate for both laboratory and field studies . It has also been the subject of an expressed sequence tag (EST) sequencing project, permitting the construction of cDNA microarrays . Copper is an essential element that is also highly toxic to soil invertebrates in high concentrations. Hence, as well as inducing general toxic-response pathways, there will also be specific biological mechanisms for copper handling that may be expected to be perturbed, and thus copper is an excellent model toxin for demonstrating an integrative ecotoxicogenomics approach. We exposed worms to sub-lethal levels of copper in a semi-field situation using buried mesocosms, and monitored the dose response using transcriptomics (using a cDNA microarray fabricated with 8,129 EST reporters representative of all consensus gene objects generated from 17,225 high-quality ESTs) and metabolomics (using proton NMR spectroscopy). We also measured changes in reproduction and general condition as measured by body weight and by a well-characterised cellular bioassay; monitoring these functional endpoints is important for phenotypic anchoring of the omic data [15, 34]. Our aim was to look for links between the different levels of information (metabolism, transcription and functional) to obtain stronger inferences about the mechanistic effects of copper than could be gathered from the individual datasets alone. A secondary aim was to test the hypothesis that copper exposure up-regulates histidine metabolism in L. rubellus [35, 36].
Results and discussion
Ecological and functional endpoints (survival, weight change, reproduction rate and neutral red retention by coelomocytes) have been reported in a previous study . In brief, there was no effect on mortality, confirming that the exposure was appropriate for probing sub-lethal molecular responses, but the other functional assays all showed a response at medium to high levels of copper (40 to 160 mg/kg and above).
We were able to reliably profile 42 different small-molecule metabolites across all samples using 1H NMR spectroscopy and software for assisted manual fitting of chemical standards (Table 1). The manual fitting approach has been shown to give high-quality data on single compound concentrations . The metabolomic analysis rests critically on the quality of the data obtained from this step; given the high degree of spectral overlap in one-dimensional proton spectra, and consequent difficulty in fitting the data, care is needed for reliable assignment. As well as assignment, the quality and reproducibility of the spectral fitting step is also a valid concern. We acquired data for one sample for five instrumental replicates; hierarchical cluster analysis (HCA) of all compound concentrations shows that these replicates are more similar than any other two sample spectra (Figure 1). This confirms not only, as expected, the high instrument precision of NMR , but also that our fitting of compounds was very reproducible. Compound assignments were made on the basis of the chemical shift and multiplicity of standards in the Chenomx software library and from online databases . In addition, a two-dimensional correlated spectroscopy (COSY) spectrum was acquired for a representative sample, and used to help confirm assignments. We have assigned 2-hexyl-5-ethyl-furan-3-sulfonic acid (HEFS) previously, and lombricine is assigned on the basis of its structure and expected high concentration in earthworms [41–44]. We also ran spectra of authentic compounds and spiked them into the sample when assignments were doubtful, for example, for purine nucleotides, which have singlet resonances in the aromatic region that are close in frequency, and are thus usually unreliable to assign from database values alone. In addition, a compound which we had previously tentatively assigned as N-α-methylhistidine  was found to be N,N-α-dimethylhistidine (DMH) after comparison with authentic standards. However, we could not obtain a sufficiently pure sample of DMH to use for quantitation, so we averaged concentrations for 1-methylhistidine and dimethylglycine, fitted to the imidazole and N-methyl protons of DMH, respectively. We were able to fit almost all of the peaks in the spectra (Additional file 1). The few remaining unfitted resonances include some at around 8.00 and 6.00 ppm (probable pyrimidine nucleotides), 6.60 and 1.16 ppm (believed to be an uncharacterised metabolite/breakdown product of HEFS) and 5.50 ppm (probable phosphosugar on the basis of its multiplicity and chemical shift, neither glucose-1-phosphate nor glucose-1,6-bisphosphate by comparison with authentic standards). However, these unfitted peaks represent only a tiny fraction of the overall intensity. Finally, one sample in the 480 mg/kg dose group had a very high xanthine concentration (more than 20 times higher than the average for all other worms), which had high leverage in multivariate analyses, and so this one data point (not the whole sample) was treated as a missing value.
Factor analysis showed a clear relationship with copper on axes 1 and 2; in addition, it was obvious that at the highest-concentration dose (480 mg/kg soil) the worms were very different to the others along the second orthogonal axis (Figure 2). This was also true following varimax axis rotation, although, in this case, the copper-related differences dropped down to axes 2 and 3. HCA showed that the high-concentration samples (480 and 160 mg/kg) formed a separate cluster to the low-concentration samples (0 and 10 mg/kg), with the intermediate 40 mg/kg worms falling into both of these clusters (Figure 1). Given that not many samples were available for cross-validation, and that the metabolic responses to copper were clearly non-linear, we chose not to use additional supervised multivariate methods.
There is currently no 'metabolite ontology' database that allows one to systematically annotate individual compounds into different groups; approaches based on network topology for the definition of elementary modes are exciting, but are still under development, and cannot be applied at the whole-organism level for metazoans . Hence, we performed a simple annotation based on existing biochemical knowledge. We note that our annotations were not performed systematically; it is also important to realise that many metabolites should be annotated to multiple categories, which we did not attempt to do. Nevertheless, this simple approach may be a useful first step in interpreting metabolomics data.
Several metabolite groups exhibited a coordinated response; this was especially clear for the lipophilic amino acids (Figure 2). Several other possible coordinated responses were identified, including sugars, nucleotides, organic acids (including several synthesised through the citric acid cycle) and, possibly most interestingly, 'cell membrane compounds' (Figure 2). It should be noted that attempting to assign 'metabolite ontologies', at however limited a level, offers an option that simply does not exist to any great extent for genes: that of chemical similarity. (Some kind of chemical classification of genes might also be made, for example, on the basis of percentage occurrence of DNA bases, but this would be of very limited applicability.) Of the groups we have labelled in Figure 2, the 'cell membrane' group does not include any compounds that directly form cell membranes, but, instead, those that could be involved in either the anabolism or catabolism of lipid compounds. It should be noted that this is the most structurally diverse group that we labelled, containing an amino acid (glycine, which is a direct precursor to serine and, hence, sphingolipid metabolism, and is a breakdown product of choline and betaine), as well as metabolites that form polar lipid head groups (inositol compounds and phosphoethanolamine). In addition, we assigned the earthworm metabolite HEFS to this group: the very high tissue concentration of HEFS implies it has some kind of structural or stabilising role. As amphiphilic compounds are known to stabilise cell membranes during desiccation stress , and desiccation is a common and severe hazard faced by earthworms, we argue that the high HEFS concentrations probably reflect a role in membrane stabilisation.
Having identified these potential groups through multivariate pattern-recognition analysis, we examined the actual data for coordinated functional group responses in more detail; selected groups are shown in Figure 3, which represents the percentage difference compared with controls (that is, directly equivalent to fold change). The coordinated response of the lipophilic amino acids is, again, very clear, and shows a non-linear response to copper: all of these metabolites are reduced at intermediate concentrations, and increased again at the highest level of 480 mg/kg (Figure 3A). In contrast, the membrane compound group shows a definite increase in response to copper (Figure 3B). The HEFS response also fits well within this group, supporting its classification as a 'membrane compound'. Even within these functional groups, fine-scale inter-metabolite relationships can be observed: myo-inositol and scyllo-inositol are the two most highly correlated metabolites within this group, with an apparently sigmoidal response to copper, presumably indicating metabolism through a common enzyme from lipid breakdown/turnover. The organic acids also show a sigmoidal-type relationship, with the compounds directly produced by the citric acid cycle the most closely related; lactate should probably not be considered part of the same functional group (Figure 3C). The remaining plots do not show grouping to the same extent: the sugars glucose and mannose both decrease in response to copper, but glucose-6-phosphate is not evidently correlated with these in any way (Figure 3D). Similarly, there is no obvious correlated response for the nucleosides as a group (Figure 3E); although nicotinate and uridine are both very highly correlated, the fold changes for these two compounds are very small, less than 10%. Larger changes are seen for inosine and xanthine, which decreased in response to copper; both of these are breakdown products, of purine and pyrimidine nucleotides respectively. Finally, the amino acids (excluding glycine and the lipophilic amino acids) showed little overall grouping. 3-Methylhistidine decreases in response to copper, while glutamine shows large variations (Figure 3F). In general, the fold changes are very small, with the large majority being less than twofold up or down. This contrasts with the transcript data, where larger fold changes are observed.
Non-targeted metabolite profiling has often been explicitly categorised as a biomarker discovery tool [48, 49], with the implication that the aim of a metabolomic study is to discover a small number of biomarkers in amongst the biological noise. This approach may be advantageous when the desired outcome is a screening tool (for instance, for disease, or for environmental pollution), when a single (or few) robust biomarker(s) will enable the use of simpler and more robust predictive models (for example, classical linear discriminant analysis ) and targeted detection. An additional advantage not often made explicit is that this greatly enhances the transferability between different laboratories and, hence, the overall potential use and scientific value of the biomarkers. However, this approach is limited when seeking to integrate large datasets and to relate these to biological functions. In the current study, considering the coordinated responses of functionally related metabolites showed clear group responses, which would not have been obvious if only the most strongly changing metabolites had been investigated.
Histidine compound metabolism
Gibb et al  reported that histidine increased in earthworm tissue in response to copper contamination at levels comparable to the present study (up to 160 mg/kg soil), and speculated that this might represent a metabolite-level cellular response for direct detoxification of copper, as observed for the increase in free histidine seen in plant tissues in response to nickel . In addition, histidine and DMH in L. rubellus were affected by metal contamination in worms sampled from the field . We thus had a particular interest in histidine and related compounds, and examined their responses to copper in detail in addition to the general analysis above. We looked at both tissue and normalised concentrations (relative to alanine); it should be noted that the calculation of tissue concentrations assumed 100% extraction of metabolites and, hence, these values are likely to slightly underestimate true concentrations. There was a weak negative correlation with copper for histidine and 3-methylhistidine, and no apparent relationship with DMH (Figure 4). The earlier study used an aqueous tissue extraction , which would not prevent metabolic activity, especially as earthworm proteases are particularly active , and therefore their results may have included contributions from protein-bound histidine. It is certainly likely that a histidine-rich protein could be responsible for detoxification by binding copper . Thus, we also extracted earthworm tissue in water, which was left at room temperature for 24 hours to allow any enzymatic activity to go to completion, and re-analysed the samples using NMR spectroscopy. Aqueous extraction causes extensive proteolysis, inferred from a large increase in free amino acid concentrations, which is effectively complete after one or two hours (data not shown). In general, the observed concentrations of amino acids were approximately two orders of magnitude higher than in chloroform/methanol extracts, and the separation of the different copper dose groups by factor analysis had completely disappeared (data not shown). Considering the histidine compounds only, there were essentially no dose-related differences compared with the chloroform/methanol extracts, indicating that there was no increase in protein-bound histidine that might reflect an increase in metal-binding proteins. We currently have no sure explanation of why our results are different to those previously observed. However, a note of caution has been sounded about the genetic variability of Lumbricus species used in laboratory tests , and we have previously observed high variability in histidine levels in L. rubellus that may have resulted from genetic differences between populations . The population used by Gibb et al  may have had a different genetic background.
In addition to the polar fraction, we also analysed the total lipid fraction by 1H NMR. Lipid spectra are dominated by resonances from fatty acids; because these are all chemically similar, it is not usually possible to identify individual compounds by NMR, and thus it should be thought of as a way of obtaining information on the relative distribution of different moieties, which still constitutes useful biochemical information . Given this, it was not possible to fit individual compound concentrations as for the polar compound data, and so we used a different approach for data analysis. A set of 57 integral regions was selected manually (Additional file 2); note that this does not equate to 57 separate compounds. A heat map of fold change relative to the control group shows quite clearly that the worms at the highest dose level, 480 mg/kg, are very different from the other groups (Figure 5). This is also visible on a single-sample (that is, not averaged) basis (Additional file 3). Factor analysis showed that one sample (one of the 480 mg/kg dose group) was an extreme outlier, which was also apparent on examining the original spectra, and so we excluded this sample and re-analysed the data. This showed clearly, using two different forms of scaling, that the 480 mg/kg worms were separated from the control and low-dose group worms along axis 2, while the 160 mg/kg worms were intermediate (Figure 5). Inspection of the original spectra confirmed that there were real and visible differences in spectral features (Additional file 4). Regions A and C-I represent signals that were decreased by copper exposure; A represents vinylic protons from unsaturated fatty acids, E represents protons allylic to two double bonds from polyunsaturated lipids and G probably represents protons from allylic methylenes that are not adjacent to another double bond. C represents glycerol peaks from triacylglycerols and D represents glycerol peaks from glycerophospholipids. F possibly represents plasmalogens, and H and I both represent signals from terminal methyls. (All assignments based on Sze and Jardetzky .) Region B contains a peak from an unassigned metabolite that is only present in the copper-treated worms, and is absent in the controls. Even though this peak is of very low intensity compared with the largest signals in the spectrum (for example, 0.0075% of the intensity of the main unsaturated methylene resonance following local baseline correction), it is a very specific response and so could still prove a useful future biomarker, possibly using some more sensitive detection method following compound identification.
There were 8,029 cDNA reporters on the microarray. Of these, 7,107 consistently yielded high-quality data and 1,705 exhibited a change in the level of transcription of more than twofold in one treatment group (Additional file 5). One-way analysis of variance (ANOVA), with correction of the error rate for multiple tests following the approach of Benjamini and Hochberg , revealed 329 of these transcripts as significantly altered (p < 0.05; see Additional file 6). Direct functional information relating to earthworm transcripts is negligible. However, approximately one-third display a significant homology (BlastX significance greater than 10-10) to human orthologues, and so we concentrated on these transcripts for further analyses. There were clear differences of expression in the different dose groups (Figures 6A and 6B). Pattern recognition analysis using principal components analysis (PCA) indicated that the major effect was at the two highest dose groups, which were clearly separated from the other samples (Figure 6C).
We used a K-means clustering approach to show that the transcripts fell into four response groups (clusters), which displayed similar dose-dependent profiles with respect to copper exposure. One group comprised of genes showing a general dose-dependent increase with respect to copper (Figure 6D), and two groups displayed the converse relationship being separated by the degree of down-regulation observed in the two lower copper doses (Figure 6E). The most intriguing cluster, however, were those transcripts that were up-regulated when exposed to a low dose of copper with a subsequent down-regulation at higher doses (Figure 6F). This profile indicates that the organisms' hormetic response to copper is also manifest at the level of transcription .
To complement the purely statistical approaches that provide insight into the coordination of transcript response, we also analysed the bias in gene function within transcripts (selected by unsophisticated filtering) to identify those that displayed more than twofold changes in transcript level in response to any of the copper exposures. We chose this group to be 'inclusive' because, although we recognise that it includes a small quantity of noise due to the quantity of transcripts analysed, it provided a more complete representation of functional group effects. This was especially desirable given that, because functional assignment was based on the top human orthologues, only about one-third of the reporters could be assigned confidently. We used a background created by those human orthologues assigned to the complete reporter set to calculate the representation bias for gene ontology terms for those orthologues assigned to transcripts exhibiting copper-induced twofold change in expression . The ontology groups that were significantly overrepresented (p < 0.1) are summarised in Additional file 7. The clear disruption in genes associated with mitochondrial electron transport strongly implies a copper-induced mitochondrial dysfunction, consistent with previous observations in other organisms [59–61]. Even a slight reduction in the capacity for adenosine triphosphate (ATP) production by aerobic respiration would result in a redistribution of energy production through anaerobic processes, and this is indeed evident through the changes in genes involved with sugar mobilisation for energy production (Additional file 7). We also observed transcript changes that we interpret as representing functional interactions between copper and other essential metal ions, in particular changes in expression of calcium- and iron-binding proteins. Finally, there were clear alterations in lipid metabolism at the transcript level, complementing the observations made at the metabolite level.
The power of calculating functional representational bias is inherent in the cumulative nature of the data obtained from a profiling approach. However, it is also important to extract specific transcript profiles of genes with an established and expected functional link to copper (Figure 7, and Additional file 8). For example, metallothioneins are widely responsive to metal exposure, including copper for particular isoforms, in a large number of organisms. Hence, we were interested in determining any such specific responses to copper in L. rubellus. The observed up-regulation of multiple reporters (Figure 7B) indicates that earthworm metallothionein, surprisingly [62, 63], is induced by copper at a far lower exposure concentration level than by cadmium. Copper exposure also increased general toxic stress in the worms, as shown by the induction of HSP70 and HSP40 (Figure 7C); our observations here from the transcript profiles are consonant with previous targeted immunochemical analyses of HSP70 levels in earthworms exposed to copper and metalliferous soils [64, 65]. Copper toxicity also causes an increase in generation of reactive oxygen species, as a result of the mitochondrial dysfunction discussed above. This could explain the increase we observed in specific glutathione-S-transferases (Figure 7D), also previously seen in response to copper toxicity . These reactive oxygen species also damage DNA [67, 68]. Similar DNA damage is also highly likely to have been a factor in the current study, as reflected by the alteration in levels of DNA repair enzymes and of enzymes implicated in cell cycle control, with both increases and decreases in transcript levels observed in response to copper (Figure 7E). Ultimately, copper-induced damage leads to apoptosis [69, 70], which again is consistent with the transcriptomic data, with down-regulation of apoptotic regulators (Figure 7F). Summarising, the responses discussed here all represent plausible modes of cellular disruption, and all of them have been demonstrated previously to be induced by copper exposure. In sharp contrast, when we consider the response of a number of established control genes  and compare them with the other groups we have selected here, we see that they are not copper responsive (Figure 7A). This strongly supports the technical and biological validity of the transcript data. Further external validation is given by comparison with the metabolomic data, where the observations of metabolites are highly consistent with the metabolic changes suggested by the microarray data.
Samples for both microarray and metabolomic analysis were chosen by selecting worms with the best condition from each independent mesocosm, thus avoiding obviously diseased and/or parasitised worms. This has two major advantages for the molecular data. First, it avoids major confounding factors, for example, from chance-affected worms. Secondly, even in the case that there would be an interaction with copper treatment, for example, through copper toxicity decreasing the resistance to external factors, one would expect the molecular mechanisms induced to be general and stress-related. Thus, even though we may have selected individuals that did not provide the most accurate representation of ecological fitness under metal stress, the molecular endpoints will be more closely related to specific copper-induced mechanisms.
Integration of the different datasets
There are multiple possibilities for integrating omic datasets, but these can be considered essentially to fall into two categories: statistical, that is, relying wholly on data-driven associations between variables (for example, [4, 72]); and knowledge-based, that is, using prior knowledge about biological organisation and pathways/networks (for example, [73, 74]). Clearly, an optimal solution would use information from both of these approaches; but then this, in the limit, effectively becomes a full model of metabolism, and this is not currently achievable, even for highly controllable unicellular model organisms. Here we have chosen a separate analysis of datasets and subsequent combination of observations.
Specific genes that are up- or down-regulated were identified from the transcriptomic data (Table 2). Relatively few could be unambiguously annotated to metabolic enzymes; two of these are from carbohydrate metabolism, mannosidase and maltase-glucoamylase. Fortuitously, the metabolic substrate and product of these enzymes, respectively, could be identified in the NMR spectra. Figure 8 shows that in both cases there is a negative correlation between transcript and metabolite levels. This can be readily rationalised for the mannose/mannosidase relationship, that is, increase in the catabolic enzyme results in a decrease in substrate concentrations, but the opposite argument should hold for glucose/maltase-glucoamylase. This is a useful reminder that it may be dangerous to over-interpret metabolite pair relationships such as these, as the most important factor, metabolic fluxes, cannot be inferred directly from gene transcript levels or metabolite concentrations alone. More interestingly, these molecular data could be very clearly related to ecologically important endpoints. It is clear that a high-sugar/low-transcript level is, in both cases, associated with the low copper doses, and there is a distinct change in both metabolite and transcript concentration for the high-dose samples (Figure 8). Moreover, the worms have a positive energy balance for the first condition (low-dose/high-sugar/low-catabolic mRNA level), that is, they gain weight over the course of the experiment, but a negative energy balance, losing weight, for the second condition (high-dose/low-sugar/high-catabolic mRNA level).
In addition, a more powerful approach than looking for alterations in single genes is to search for functional gene groups that are co-ordinately regulated. This is achievable for the transcript data, given that Gene Ontology (GO) classifications are available. Figure 9 shows that a number of transcripts represented with overrepresented ontology groups mapped onto specific Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and it illustrated the impact of copper on genes involved in the electron transport pathway and converse impact on glycolysis, thus indicating the severe remodelling in energetic metabolism upon exposure to copper. The majority of transcripts for electron transport were up-regulated (Figure 9C), and this picture was even clearer for the glycolytic transcripts, which were all up-regulated (Figure 9D). Conversely, phosphoenolpyruvate carboxykinase (PEPCK) transcript levels were decreased by copper (PEPCK was a member of the GO class on which the analysis was based, even though it is not shown on the pathway fragment in Figure 9). Although this represents only one enzyme, it indicates that the gluconeogenic enzymes, as would be expected, were regulated oppositely to the glycolytic enzymes.
These observations are directly reflected by the metabolomic data, where a set of organic acids that are citric acid cycle intermediates were altered as a functional group, and there was a decrease in free sugar concentrations (Figure 3). In addition, when these observations are considered together with the links to specific transcripts and functional endpoints (Figure 7), and the changes in metabolites directly involved in maintaining cellular energy levels (vide infra), it is clear that energy metabolism in L. rubellus was profoundly impacted by copper treatment.
Interpretation of energy reserve data
Some copper-responsive metabolites were those related to maintenance of cellular energy levels, in particular adenosine phosphates. In general, adenosine monophosphate (AMP) and adenosine diphosphate (ADP) levels were both high and correlated to each other, and also showed some positive correlation to copper. This, combined with the lack of any detectable 1H signal for phospholombricine, implies that the worms' cellular energy reserves were exhausted, and thus not truly representative of in vivo metabolism. This is not a simple artefact of the extraction procedure: for comparison, a neutralised 6% perchloric acid extract tells a similar story about the energetic state, that is, AMP and ADP both higher in concentration than ATP, and also completely lacking a signal from phospholombricine (Additional file 9). The worms were sampled under conditions that would have preserved metabolic integrity, by flash-freezing into liquid nitrogen, storage at -80°C, and lyophilised tissues extracted directly into ice-cold solvent/water mixtures. In addition, the high quality of the mRNA extracted for the microarray analyses confirms the preservation of the biological integrity of the samples. Thus, we conclude that the general trend of the adenosine phosphates (AXP) data (ATP decreasing, and AMP and ADP increasing, with increased soil copper) is a reflection of a general disturbance in energy metabolism, as also shown by decreases in free sugars, and up-regulation of transcripts for breakdown of these sugars and catabolism of energy reserves, even if it may not be a true 'snapshot' of in vivo metabolism. Indeed, any extraction procedure is necessarily selective and different protocols will give different windows on the metabolome. Techniques such as high-resolution magic-angle-spinning of tissue biopsies (and even larger samples) are increasingly used for biochemical and biomedical investigations , even though this will lead to far more extensive enzymatic changes than the extractions used here.
There was a clear impact on both metabolomic and microarray profiles at the second-lowest additional copper concentration of 40 mg/kg. This is either as sensitive or more sensitive than the functional parameters measured, such as weight change, reproduction and the neutral red retention time (NRRT) bioassay, all of which had lowest-observed-effect concentrations (LOECs) of 160 mg/kg or higher, and much more sensitive than an endpoint based on mortality for which no significant effect was found at any of the tested exposure concentrations (that is, no-observed-effect concentration (NOEC) greater than 480 mg/kg) . This sensitivity seems a fairly general rule, at least for laboratory experiments, where metabolomics has been shown to be extremely sensitive to individual chemicals compared with traditional endpoints . Such sensitivity provides an interesting insight into the likely value of such analyses for regulatory testing regimes. Here omic approaches may provide insights on low effect level metabolic dysfunctions which may have implication over a life-time exposure, but may not be revealed in the time-limited chronic bioassays that form a cornerstone of modern chemicals policy .
Like some of the observed phenotypic responses (for example, reproduction rate and weight change), the omic data (both transcripts and metabolites) very clearly demonstrate a non-linear response to copper. This was also true for the weight change data, with an initial increase in body weight followed by a clear decrease at higher concentrations. It has been argued that hormesis is a very general response to potential toxins, with low-level exposure either being beneficial or stimulating a compensatory response . While the response seen for copper may merely represent the fact that this metal is both essential and toxic for earthworms, it is currently the case more generally that a mechanistic basis for understanding such 'hormetic' responses is lacking. While initial suggestions have been made to explain hormesis in terms of an overcompensation of homeostasis to disruption that is mediated by the different affinities of stimulatory and inhibitory regulatory pathways , such suggestions have yet to be confirmed or refuted in empirical studies. Our data clearly show that the hormetic response in this case is manifest at a molecular level, indicating a possible role for metabolomics in understanding the mechanistic nature of the response.
Drawing reliable inferences from omic data is often difficult, especially for non-model organisms which may be less well annotated than standard laboratory models. Here we have used a dual metabolomic and transcriptomic strategy, with each element providing complementary data, thus reinforcing the inferences made from the independent datasets. Interpreting the data in terms of functional groups, for both transcript and metabolite data, showed clear molecular group responses to copper. The integrated dataset indicated a clear alteration of energy metabolism as a result of sub-lethal copper exposure, with an increased switch to metabolism of stored carbohydrates, presumably as a consequence of copper interfering with mitochondrial function and reducing the amount of energy available from oxidative phosphorylation. These molecular impacts resulted in higher-level endpoints, including a reduction in body weight at high copper levels, reflecting the changes in energy metabolism, and a decrease in lysosomal integrity, indicating effects on membranes that could parallel the observed impacts on mitochondrial function.
Exposure and sampling
A full description of the earthworm exposure and sampling conditions is given by Spurgeon et al . Briefly, worms were exposed for 70 days under field conditions in mesocosms to controlled levels of soil copper (0, 10, 40, 160 and 480 mg additional copper per kilogram dry weight soil; original copper content of the soil was 16.1 mg/kg), with four replicate mesocosms per dose level. A number of functional endpoints were measured including survival, weight change (average value per mesocosm), NRRT and cocoon production rate. Fifteen adult worms were exposed per mesocosm, and at the end of the experiment, three worms were used to measure tissue copper concentrations; four worms were used to measure NRRT and the three worms maintaining highest condition in each replicate were pooled and used for both transcriptomic and metabolomic analysis. This avoided inclusion of diseased or heavily parasitised worms within the omic experiments that may have otherwise confounded results. The worms were flash-frozen in liquid nitrogen, and ground to powder in a mortar and pestle under liquid nitrogen to obtain a single sample for each mesocosm. Total mRNAs from the samples were then isolated using established protocols [79, 80]. The remaining sample was lyophilised without allowing the sample to thaw, and stored at -80°C until extracted for metabolomic analysis.
The L. rubellus EST project  has established a database of over 17,000 ESTs from unexposed and chemically exposed earthworms (see  for details). All sequences were assembled into clusters and annotated using the PartiGene pipeline . To fabricate the glass slide cDNA microarray, a representative EST (usually the longest) was selected from each of the 8,029 clusters assembled from the ESTs (see  for full details). This sequence was polymerase chain reaction (PCR)-amplified and aliquots (5 μl) of concentrated products mixed in 384 well plates with an equal volume of dimethyl sulphoxide (DMSO). These were then printed onto Ultra-GAP glass slides (Corning) using 48 SMP3 pins (Telecham) mounted in a Spotarray 72 (Perkin-Elmer). Landmarks were introduced at the left-hand corner of each sub-array by the introduction of five replicates of the Lucida Scorecard (Amersham) gene reporters, which show no cross-reactivity to earthworm transcripts. All reporters were cross-linked to the slide by baking at 80°C for 2 hours, and UV cross-linking.
Lucida Scorecard test spike (Amersham Life Sciences) was added to 10 μg of total RNA prior to oligo-dT reverse transcription and coupling to Cy3 using an indirect amino amyl procedure. Clean-up of labelled targets, yield and integrity were all measured according to Owen et al . A reference design was used in which approximately 30 pmol of Cy3 labelled target RNA was hybridised against the common oligonucleotide reference (representing 30 pmol of Cy5). The reference used was a 65–70 mer oligonucleotide designed against the vector sequence between the amplification primer binding site and cDNA insert. Use of this universal reference design allowed a comparison of slides to be made both within and between experiments. After hybridisation (18 hours), slides were washed and imaged according to Owen et al . Array images were subjectively quality controlled for artefacts that would compromise quantification such as background effects and spot morphology prior to image analysis with Imagene 5.0 (Biodiscovery). Subsequently, the calibration standards from the Lucida Scorecard were analysed to objectively assess the sensitivity range and to define both saturation and background readings (Additional file 10). The microarray data can be accessed through the NEBC file store  (EnvBase accession number [EGCAT:4024]).
Statistical analysis of microarray data
Numeric data were imported into LimmaGUI [85, 86] which allowed the subtraction of background measurements, generation of additional quality control plots and subsequent normalisation using a within-array Tiplowess manipulation. Abnormally distributed samples were excluded from further analysis. Normalised data were subsequently imported into GeneSpring 7.3 (Agilent Technologies, Stockport, UK) and represented relative to the median expression within the control group. Final, quantity assessment was performed by visualising a box plot of the normalised data together with the generation of MA plots for the average data from each dose (Additional file 11). Abnormally distributed samples were excluded from the analysis.
Using GeneSpring, the three control slides and three slides for the exposure concentration that most closely matched the cocoon production EC10 for each chemical were used to identify chemically responsive genes. The data were first filtered to include only spots flagged as present in 3 out of the 15 slides. This filtered dataset was then used to generate three gene lists for each dataset. These were: (1) genes with a more than twofold difference in mean expression between control and exposed samples; (2) genes with significantly different (p < 0.05) expression between control and exposed samples after correction for multiple sample testing .
Tissue samples (20–30 mg) were homogenised into 3 ml of extraction solvent (Heidolph SilentCrusher S) using a modified Bligh and Dyer protocol [87, 88], in which a monophasic chloroform/methanol/water extraction is followed by the addition of water and chloroform, splitting the sample into two phases. The proportion of water was adjusted for the fact that lyophilised tissue was used, assuming a 90% water content in normal tissue. Both fractions were dried at 40°C using a rotary vacuum concentrator. The polar fraction was then resuspended in 0.65 ml of NMR buffer (100 mM phosphate buffer pH 7.0 in 2H2O, also containing 0.98 mM sodium trimethylsilyl-2,2,3,3-2H4-propionate (TSP); the 2H2O provided a field frequency lock for the spectrometer and reduced the signal from water protons) and centrifuged (5 minutes, 16,000 g). The supernatant (0.6 ml) was then transferred to 5 mm NMR tubes. The aqueous samples were analysed at 300 K on a 14.1 T DRX600 Avance NMR spectrometer (Bruker BioSpin, Rheinstetten, Germany) with a 600 MHz proton resonance frequency, equipped with a 5 mm broadband inverse probe. The samples were run using an automatic tube changer, over a period of about 12 hours, during which time they were kept at room temperature; the samples were loaded onto the tube changer in randomised blocks. A one-dimensional NOESY sequence with a mixing time of 100 ms was used for water suppression of the residual HOD in the NMR buffer, using a 30 Hz presaturation pulse. The spectra were acquired for 128 transients, with four dummy scans, into 32 K data points over a 12 kHz spectral width. A 3.5-second longitudinal relaxation recovery delay was added for each transient, giving a recycle time of 5 seconds. The chloroform fraction was resuspended in 0.65 ml of CDCl3 containing 0.03% tetramethylsilane (TMS) and transferred to 5 mm NMR tubes. The lipid samples were then analysed using a Carr-Purcell-Meiboom-Gill (CPMG) sequence, with a 8.65-second longitudinal relaxation delay giving an approximately 10-second recycle time (the CPMG sequence gave improved baselines and slight reduction of broad resonances compared to a simple pulse-acquire experiment). The raw data (free induction decays (FIDs)) for all spectra, and the fitted compound data for the polar extracts (vide infra), are available through the NEBC file store  (EnvBase accession number [EGCAT:4024]). The fitted compound data are also provided as Additional file 12.
NMR processing and data analysis
The spectra were initially processed using iNMR 2.2.7 (Nucleomatica, Molfetta, Italy). The summed FIDs were multiplied by an exponential window function equivalent to 0.5 or 1 Hz line broadening (aqueous and lipid spectra, respectively). They were then zero-filled by a factor of 1.5, and Fourier transformed. Phasing was carried out using the automatic 'metabolomic phase correction' option, and adjusted manually where necessary; baseline correction was performed using an automatic first-order polynomial fit. All spectra were referenced to TSP/TMS at 0 ppm. The spectra from polar (aqueous) samples were further analysed using Chenomx NMR Suite 4.6 (Chenomx, Edmonton, AB, Canada). This software fits idealised spectra made up of combinations of Lorentzian line shapes, based on spectra of authentic standards, and estimates compound concentrations using TSP as an internal quantitation standard. Compounds were fitted using the proprietary Chenomx 600 MHz library, to which we added standards for glucose-6-phosphate, HEFS and lombricine. HEFS and lombricine are not commercially available; we purified HEFS from earthworm tissues using solid-phase extraction with a mixed-mode C18/anion exchange phase. Only a single peak for lombricine was fitted, based on a spectrum of an existing tissue extract acquired under fully relaxed conditions, and concentration estimated by comparing this with the integral of a known concentration of TSP. Thus, it is likely that the absolute accuracy of the lombricine concentrations will be lower than for the other compounds reported in this study, although the relative levels (precision) will be comparable. The fitted compound concentrations are available in Additional file 12.
For pattern recognition, data were normalised following the method of Dieterle et al  in which each profile is compared with a representative reference (in our case, a median of all sample profiles). The relative fold change for each variable in turn is calculated, and all values for a spectrum are then divided by the median fold change for that spectrum. Data were log-transformed by log(n i + 0.018) for concentrations (expressed as mM); the value of 0.018 was chosen because it effectively removed the correlation between intensity and standard deviation for a series of five technical replicates, that is, increasing homoscedasticity (the principle is discussed elsewhere [90, 91]). Factor analysis and hierarchical cluster analysis on non-centred data was carried out using Aabel 2.2 (Gigawiz, Tulsa, OH, USA). (Note that the factor analysis here is exactly equivalent to PCA carried out using the correlation matrix, but we have retained the term 'factor analysis' for consistency with the terms used in the software package.) The lipid spectra were integrated within selected regions using iNMR (integral boundaries given in Table 2).
Moore MN: Biocomplexity: the post-genome challenge in ecotoxicology. Aquat Toxicol. 2002, 59: 1-15.
Ankley GT, Daston GP, Degitz SJ, Denslow ND, Hoke RA, Kennedy SW, Miracle AL, Perkins EJ, Snape JR, Tillitt DE, Tyler CR, Versteeg D: Toxicogenomics in regulatory ecotoxicology. Environ Sci Technol. 2006, 40: 4055-4065.
Miracle AL, Ankley GT: Ecotoxicogenomics: linkages between exposure and effects in assessing risks of aquatic contaminants to fish. Reprod Toxicol. 2005, 19: 321-326.
Craig A, Sidaway J, Holmes E, Orton T, Jackson D, Rowlinson R, Nickson J, Tonge R, Wilson I, Nicholson J: Systems toxicology: integrated genomic, proteomic and metabonomic analysis of methapyrilene induced hepatotoxicity in the rat. J Proteome Res. 2006, 5: 1586-1601.
Waters MD, Fostel JM: Toxicogenomics and systems toxicology: aims and prospects. Nat Rev Genet. 2004, 5: 936-948.
Waters M, Yauk C: Consensus recommendations to promote and advance predictive systems toxicology and toxicogenomics. Environ Mol Mutagen. 2007, 48: 400-403.
Heijne WH, Kienhuis AS, van Ommen B, Stierum RH, Groten JP: Systems toxicology: applications of toxicogenomics, transcriptomics, proteomics and metabolomics in toxicology. Expert Rev Proteomics. 2005, 2: 767-780.
Robertson DG: Metabonomics in toxicology: a review. Toxicol Sci. 2005, 85: 809-822.
Lindon JC, Keun HC, Ebbels TMD, Pearce JMT, Holmes E, Nicholson JK: The Consortium for Metabonomic Toxicology (COMET): aims, activities and achievements. Pharmacogenomics. 2005, 6: 691-699.
Berg van den M, Tamis WLM, van Straalen NM: The food web approach in ecotoxicological risk assessment. Hum Ecol Risk Assessment. 1998, 4: 49-55.
Parmelee RW, Wentsel RS, Phillips CT, Simini M, Checkai RT: Soil microcosm for testing the effects of chemical pollutants on soil fauna communities and trophic structure. Environ Toxicol Chem. 1993, 12: 1477-1486.
Forbes VE, Calow P: Species sensitivity distributions revisited: a critical appraisal. Hum Ecol Risk Assessment. 2002, 8: 473-492.
Newman MC, Ownby DR, Mezin LCA, Powell DC, Christensen TRL, Lerberg SB, Anderson BA: Applying species-sensitivity distributions in ecological risk assessment: assumptions of distribution type and sufficient numbers of species. Environ Toxicol Chem. 2000, 19: 508-515.
Blanck H, Wangberg SA: Validity of an ecotoxicological test system – short-term and long-term effects of arsenate on marine periphyton communities in laboratory systems. Can J Fish Aquat Sci. 1988, 45: 1807-1815.
Hines A, Oladiran GS, Bignell JP, Stentiford GD, Viant MR: Direct sampling of organisms from the field and knowledge of their phenotype: key recommendations for environmental metabolomics. Environ Sci Technol. 2007, 41: 3375-3381.
Amelina H, Apraiz I, Sun W, Cristobal S: Proteomics-based method for the assessment of marine pollution using liquid chromatography coupled with two-dimensional electrophoresis. J Proteome Res. 2007, 6: 2094-2104.
Travers SE, Smith MD, Bai JF, Hulbert SH, Leach JE, Schnable PS, Knapp AK, Milliken GA, Fay PA, Saleh A, Garrett KA: Ecological genomics: making the leap from model systems in the lab to native populations in the field. Front Ecol Environ. 2007, 5: 19-24.
Dowling VA, Sheehan D: Proteomics as a route to identification of toxicity targets in environmental toxicology. Proteomics. 2006, 6: 5597-5604.
Soetaert A, Moens LN, Ven Van der K, Van Leemput K, Naudts B, Blust R, De Coen WM: Molecular impact of propiconazole on Daphnia magna using a reproduction-related cDNA array. Comp Biochem Physiol C. 2006, 142: 66-76.
Lee SE, Yoo DH, Son J, Cho K: Proteomic evaluation of cadmium toxicity on the midge Chironomus riparius Meigen larvae. Proteomics. 2006, 6: 945-957.
Knigge T, Monsinjon T, Andersen OK: Surface-enhanced laser desorption/ionization-time of flight-mass spectrometry approach to biomarker discovery in blue mussels (Mytilus edulis) exposed to polyaromatic hydrocarbons and heavy metals under field conditions. Proteomics. 2004, 4: 2722-2727.
Kuperman RG, Checkai RT, Ruth LM, Henry T, Simini M, Kimmel DG, Phillips CT, Bradley BP: A proteome-based assessment of the earthworm Eisenia fetida : response to chemical warfare agents in a sandy loam soil. Pedobiologia. 2003, 4: 617-621.
Bundy JG, Spurgeon DJ, Svendsen C, Hankard PK, Osborn D, Lindon JC, Nicholson JK: Earthworm species of the genus Eisenia can be phenotypically differentiated by metabolic profiling. FEBS Lett. 2002, 521: 115-120.
Raamsdonk LM, Teusink B, Broadhurst D, Zhang NS, Hayes A, Walsh MC, Berden JA, Brindle KM, Kell DB, Rowland JJ, Westerhoff HV, van Dam K, Oliver SG: A functional genomics strategy that uses metabolome data to reveal the phenotype of silent mutations. Nat Biotechnol. 2001, 19: 45-50.
Goodacre R, Vaidyanathan S, Dunn WB, Harrigan GG, Kell DB: Metabolomics by numbers: acquiring and understanding global metabolite data. Trends Biotechnol. 2004, 22: 245-252.
Griffin JL: Metabolic profiles to define the genome: can we hear the phenotypes?. Philos Trans R Soc Lond B Biol Sci. 2004, 359: 857-871.
Fiehn O, Weckwerth W: Deciphering metabolic networks. Eur J Biochem. 2003, 270: 579-588.
Tiunov A, Hale C, Holdsworth A, Vsevolodova-Perel T: Invasion patterns of Lumbricidae into the previously earthworm-free areas of northeastern Europe and the western great lakes region of North America. Biol Invasions. 2006, 8: 1223-1234.
Butt KR, Lowe CN: Anthropic influences on earthworm distribution, Isle of Rum National Nature Reserve, Scotland. Eur J Soil Biol. 2004, 40: 63-72.
Deibert EJ, Utter RA: Earthworm (Lumbricidae) survey of North Dakota fields placed in the US Conservation Reserve Program. J Soil Water Conserv. 2003, 58: 39-45.
Baker GH, Thumlert TA, Meisel LS, Carter PJ, Kilpin GP: "Earthworms Downunder": A survey of the earthworm fauna of urban and agricultural soils in Australia. Soil Biol Biochem. 1997, 29: 589-597.
Spurgeon DJ, Weeks JM, Van Gestel CAM: A summary of eleven years progress in earthworm ecotoxicology. Pedobiologia. 2003, 47: 588-606.
Sturzenbaum SR, Parkinson J, Blaxter M, Morgan AJ, Kille P, Georgiev O: The earthworm Expressed Sequence Tag project. Pedobiologia. 2003, 47: 447-451.
Paules R: Phenotypic anchoring: linking cause and effect. Environ Health Perspect. 2003, 111: A338-A339.
Gibb JOT, Svendsen C, Weeks JM, Nicholson JK: 1H NMR spectroscopic investigations of tissue metabolite biomarker response to Cu(II) exposure in terrestrial invertebrates: identification of free histidine as a novel biomarker of exposure to copper in earthworms. Biomarkers. 1997, 2: 295-302.
Bundy JG, Spurgeon DJ, Svendsen C, Hankard PK, Weeks JM, Osborn D, Lindon JC, Nicholson JK: Environmental metabonomics: applying combination biomarker analysis in earthworms at a metal contaminated site. Ecotoxicol. 2004, 13: 797-806.
Spurgeon DJ, Svendsen C, Lister LJ, Hankard PK, Kille P: Earthworm responses to Cd and Cu under fluctuating environmental conditions: a comparison with results from laboratory exposures. Environ Pollut. 2005, 136: 443-452.
Weljie AM, Newton J, Mercier P, Carlson E, Slupsky CM: Targeted profiling: quantitative analysis of 1H NMR metabolomics data. Anal Chem. 2006, 78: 4430-4442.
Keun HC, Ebbels TMD, Antti H, Bollard ME, Beckonert O, Schlotterbeck G, Senn H, Niederhauser U, Holmes E, Lindon JC, Nicholson JK: Analytical reproducibility in 1H NMR-based metabonomic urinalysis. Chem Res Toxicol. 2002, 15: 1380-1386.
Biological Magnetic Resonance Data Bank – Metabolomics. [http://www.bmrb.wisc.edu/metabolomics]
Bundy JG, Lenz EM, Bailey NJ, Gavaghan CL, Svendsen C, Spurgeon D, Hankard PK, Osborn D, Weeks JA, Trauger SA, Spier P, Sanders I, Lindon JC, Nicholson JK, Tang H: Metabonomic assessment of toxicity of 4-fluoroaniline, 3,5-difluoroaniline and 2-fluoro-4-methylaniline to the earthworm Eisenia veneta (Rosa): identification of new endogenous biomarkers. Environ Toxicol Chem. 2002, 21: 1966-1972.
Euerby MR, Partridge LZ, Gibbons WA: High-performance liquid chromatographic determination of lombricine and N-phosphoryl lombricine in the earthworm by pre-column fluorescence derivatisation with o-phthaldialdehyde-ethanethiol. J Chromatogr. 1988, 445: 433-440.
Ennor AH, Rosenberg H: The isolation of N-phosphoryl-lombricine from earthworms. Biochem J. 1962, 83: 14-17.
Pant R: Isolation of lombricine and its enzymic phosphorylation. Biochem J. 1959, 73: 30-33.
Bundy JG, Keun HC, Sidhu JK, Spurgeon DJ, Svendsen C, Kille P, Morgan AJ: Metabolic profile biomarkers of metal contamination in a sentinel terrestrial species are applicable across multiple sites. Environ Sci Technol. 2007, 41: 4458-4464.
Schwartz J-M, Gaugain C, Nacher J, de Daruvar A, Kanehisa M: Observing metabolic functions at the genome scale. Genome Biol. 2007, 8: R123-
Golovina EA, Hoekstra FA, Hemminga MA: Drying increases intracellular partitioning of amphiphilic substances into the lipid phase – impact on membrane permeability and significance for desiccation tolerance. Plant Physiol. 1998, 118: 975-986.
Kenny LC, Dunn WB, Ellis DI, Myers JE, Robinson AE, Kell DB, Baker PN: Novel biomarkers for preeclampsia detected using metabolomics and machine learning. J Soc Gynecol Investig. 2005, 12: 355-
Holmes E, Tsang TM, Huang JTJ, Leweke FM, Koethe D, Gerth CW, Nolden BM, Gross S, Schreiber D, Nicholson JK, Bahn S: Metabolic profiling of CSF: evidence that early intervention may impact on disease progression and outcome in schizophrenia. PLoS Med. 2006, 3: e327-
Fisher RA: The use of multiple measures in taxonomic problems. Ann Eugen. 1936, 7: 179-188.
Kramer U, Cotter-Howells JD, Charnock JM, Baker AJM, Smith JAC: Free histidine as a metal chelator in plants that accumulate nickel. Nature. 1996, 379: 635-638.
Whiston RA, Seal KJ: Rapid production of axenic specimens of the earthworm Eisenia foetida using microcrystalline cellulose as a carrier medium for antibiotics. Soil Biol Biochem. 1988, 20: 407-408.
Devoid SJ, Etter R, Sugumaran M, Wallace GT, Robinson WE: Histidine-rich glycoprotein from the hemolymph of the marine mussel Mytilus edulis L. binds class A, class B, and borderline metals. Environ Toxicol Chem. 2007, 26: 872-877.
Lowe CN, Butt KR: Culture of commercially obtained Lumbricus terrestris L.: implications for sub-lethal ecotoxicological testing. Soil Biol Biochem. 2007, 39: 1674-1679.
Sze DY, Jardetzky O: High-resolution proton NMR studies of lymphocyte extracts. Immunomethods. 1994, 4: 113-126.
Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995, 57: 289-300.
Calabrese EJ: Paradigm lost, paradigm found: The re-emergence of hormesis as a fundamental dose response model in the toxicological sciences. Environ Pollut. 2005, 138: 379-
Dennis G, Sherman BT, Hosack DA, Yang J, Gao W, Lane HC, Lempicki RA: DAVID: Database for annotation, visualization, and integrated discovery. Genome Biol. 2003, 4 (9):
Heron P, Cousins K, Boyd C, Daya S: Paradoxical effects of copper and manganese on brain mitochondrial function. Life Sci. 2001, 68: 1575-1583.
Levenson CW, Song Y, Narayanan VS, Fitch CA, Yeiser EC: Regulation of mitochondrial cytochrome b mRNA by copper in cultured human hepatoma cells and rat liver. Biol Trace Elem Res. 1999, 70: 149-164.
Totaro EA, Pisanti FA, Glees P, Continillo A: The effect of copper pollution on mitochondrial degeneration. Mar Environ Res. 1986, 18: 245-253.
Morgan AJ, Sturzenbaum SR, Winters C, Grime GW, Abd Aziz NA, Kille P: Differential metallothionein expression in earthworm (Lumbricus rubellus) tissues. Ecotoxicol Environ Saf. 2004, 57: 11-19.
Spurgeon DJ, Sturzenbaum SR, Svendsen C, Hankard PK, Morgan AJ, Weeks JM, Kille P: Toxicological, cellular and gene expression responses in earthworms exposed to copper and cadmium. Comp Biochem Physiol C. 2004, 138: 11-21.
Nadeau D, Corneau S, Plante I, Morrow G, Tanguay RM: Evaluation for Hsp70 as a biomarker of effect of pollutants on the earthworm Lumbricus terrestris. Cell Stress Chaperones. 2001, 6: 153-163.
Marino F, Winters C, Morgan AJ: Heat shock protein (hsp60, hsp70, hsp90) expression in earthworms exposed to metal stressors in the field and laboratory. Pedobiologia. 1999, 43: 615-624.
Letelier ME, Martinez M, Gonzalez-Lira V, Faundez M, Aracena-Parks P: Inhibition of cytosolic glutathione S-transferase activity from rat liver by copper. Chem Biol Interact. 2006, 164: 39-48.
Cervantes-Cervantes MP, Calderon-Salinas JV, Albores A, Munoz-Sanchez JL: Copper increases the damage to DNA and proteins caused by reactive oxygen species. Biol Trace Elem Res. 2005, 103: 229-248.
Sagripanti JL: DNA damage mediated by metal ions with special reference to copper and iron. Metal Ions in Biological Systems. Edited by: Sigel A, Sigel H. 1999, New York: Marcel Dekker, 36: 179-209.
Zhai QW, Ji HB, Zheng ZC, Yu X, Sun LY, Liu XY: Copper induces apoptosis in BA/F3 beta cells: Bax, reactive oxygen species, and NF kappa B are involved. J Cell Physiol. 2000, 184: 161-170.
Raes H, Braekman BP, Criel GRJ, Rzeznik U, Vanfleteren JR: Copper induces apoptosis in Aedes C6/36 cells. J Exp Zool. 2000, 286: 1-12.
Sturzenbaum SR, Kille P: Control genes in quantitative molecular biological techniques: the variability of invariance. Comp Biochem Physiol B. 2001, 130: 281-289.
Urbanczyk-Wochniak E, Luedemann A, Kopka J, Selbig J, Roessner-Tunali U, Willmitzer L, Fernie AR: Parallel analysis of transcript and metabolic profiles: a new approach in systems biology. EMBO Rep. 2003, 4: 989-993.
Hirai MY, Yano M, Goodenowe DB, Kanaya S, Kimura T, Awazuhara M, Arita M, Fujiwara T, Saito K: Integration of transcriptomics and metabolomics for understanding of global responses to nutritional stresses in Arabidopsis thaliana. Proc Natl Acad Sci USA. 2004, 101: 10205-10210.
Ippolito JE, Xu J, Jain S, Moulder K, Mennerick S, Crowley JR, Townsend RR, Gordon JI: An integrated functional genomics and metabolomics approach for defining poor prognosis in human neuroendocrine cancers. Proc Natl Acad Sci USA. 2005, 102: 9901-9906.
Griffin JL: Metabonomics: NMR spectroscopy and pattern recognition analysis of body fluids and tissues for characterisation of xenobiotic toxicity and disease diagnosis. Curr Opin Chem Biol. 2003, 7: 648-654.
Viant MR, Pincetich CA, Tjeerdema RS: Metabolic effects of dinoseb, diazinon and esfenvalerate in eyed eggs and alevins of Chinook salmon (Oncorhynchus tshawytscha) determined by 1H NMR metabolomics. Aquat Toxicol. 2006, 77: 359-
Ankley GT, Miracle AL, Perkins EJ, Daston GP: Genomics in Regulatory Ecotoxicology: Applications and Challenges. 2007, Boca Raton, FL: CRC Press
Calabrese EJ, Baldwin LA: Hormesis: U-shaped dose responses and their centrality in toxicology. Trends Pharmacol Sci. 2001, 22: 285-291.
Chomczynski P, Sacchi N: Single step method of RNA isolation by acid guanidinium thiocyanate phenol chloroform extraction. Anal Biochem. 1987, 162: 156-159.
Galay-Burgos M, Spurgeon DJ, Weeks JM, Stürzenbaum SR, Morgan AJ, Kille P: Developing a new method for soil pollution monitoring using molecular genetic biomarkers. Biomarkers. 2003, 8: 229-239.
Environmental genomics of earthworms. [http://www.earthworms.org]
Owen J, Hedley BA, Svendsen C, Wren J, Jonker MJ, Hankard PK, Lister LJ, Stürzenbaum SR, Morgan AJ, Spurgeon DJ, Blaxter ML, Kille P: Transcriptome profiling of developmental and xenobiotic responses in a keystone soil animal, the oligochaete annelid Lumbricus rubellus. BMC Genomics. 2008, 9: 226-
Parkinson J, Anthony A, Wasmuth J, Schmid R, Hedley A, Blaxter M: PartiGene – constructing partial genomes. Bioinformatics. 2004, 20: 1398-1404.
NEBC file store. [http://nebc.nox.ac.uk/]
Smyth GK: Limma: linear models for microarray data. Bioinformatics and Computational Biology Solutions using R and Bioconductor. Edited by: Gentleman R, Carey V, Dudoit S, Irizarry R, Huber W. 2005, New York: Springer, 397-420.
Wettenhall JM, Smyth GK: limmaGUI: A graphical user interface for linear modeling of microarray data. Bioinformatics. 2004, 20: 3705-3706.
Bligh EG, Dyer WJ: A rapid method of total lipid extraction and purification. Can J Biochem Physiol. 1959, 37: 911-917.
Le Belle JE, Harris NG, Williams SR, Bhakoo KK: A comparison of cell and tissue extraction techniques using high-resolution 1H-NMR spectroscopy. NMR Biomed. 2002, 15: 37-44.
Dieterle F, Ross A, Schlotterbeck G, Senn H: Probabilistic quotient normalization as robust method to account for dilution of complex biological mixtures. Application in 1H NMR metabonomics. Anal Chem. 2006, 78: 4281-4290.
Purohit PV, Rocke DM, Viant MR, Woodruff DL: Discrimination models using variance-stabilizing transformation of metabolomic NMR data. OMICS J Integr Biol. 2004, 8: 118-130.
Parsons HM, Ludwig C, Gunther UL, Viant MR: Improved classification accuracy in 1- and 2-dimensional NMR metabolomics data using the variance stabilising generalised logarithm transformation. BMC Bioinformatics. 2007, 8: 234-
Kanehisa M, Goto S, Hattori M, Aoki-Kinoshita KF, Itoh M, Kawashima S, Katayama T, Araki M, Hirakawa M: From genomics to chemical genomics: new developments in KEGG. Nucleic Acids Res. 2006, 34: D354-D357.
We thank Chenomx, Inc., for access to and assistance with Chenomx NMR Suite 4.6 software. Funding was provided by NERC (grant references NER/T/S/2002/0021 and NE/D007755/1). Jennifer Owen is thanked for the fabrication of the L. rubellus cDNA microarrays. Hector Keun and Tim Ebbels are thanked for valuable discussions.
JGB carried out NMR spectroscopy, metabolomic data analysis, and drafted the manuscript. JKS developed experimental methods. FR carried out spectroscopic analysis of lipid samples. DJS and CS were responsible for earthworm exposure, sampling and functional assays. JW carried out hybridisations and microarray analysis. PK was responsible for analysis of microarray data and integrative pathway-level analysis. DJS, CS, SRS, AJM and PK participated in conception and design of the study. All authors read and approved the final manuscript.
Electronic supplementary material
Additional file 1: A 600 MHz 1H NMR spectrum of typical earthworm extract, polar fraction. (A) and (B) have an expanded vertical scale compared with (C) and (D). Resonance from HEFS (compound 19) at 6.19 ppm is not represented at its full height. Metabolite labels correspond to numbers given in Table 1. * represents an unknown compound that is a probable breakdown product of HEFS. (PDF 277 KB)
Additional file 2: Chemical shift regions for integrals of lipid extract spectral data ('missing' int01 was for internal standard TMS). (DOC 36 KB)
Additional file 4: 600 MHz 1H spectra of lipid extracts. Samples only shown from control (blue) and highest (red) dose groups. One spectrum from red group was excluded as an outlier and is not shown here. (A) Vinylic protons from unsaturated fatty acids; (B) unassigned; (C) glycerol protons from triacylglycerols; (D) glycerol peaks from glycerophospholipids; (E) protons allylic to two double bonds; (G) protons allylic to one double bond; (H) and (I) terminal methyls. (PDF 612 KB)
Additional file 5: Table of transcripts showing more than twofold change in expression as a consequence of copper exposure. (XLS 664 KB)
Additional file 6: Table of transcripts significantly changed in earthworms exposed to copper. ANOVA, p < 0.05, Benjamini and Hochberg  FDR correction. (XLS 186 KB)
Additional file 7: Table of GO terms overrepresented in transcripts whose expression is altered by copper exposure by more than twofold. (XLS 21 KB)
Additional file 8: Table of the targeted functional transcript changes observed during copper exposure of adult earthworms. (XLS 21 KB)
Additional file 9: Comparison of 6% perchloric acid (red) and chloroform/methanol (blue) extraction methods for earthworm tissue: (A) aromatic region; (B) aliphatic region. (PDF 421 KB)
Additional file 10: Assessment of micro-array sensitivity and signal linearity. Representative analysis of the fluorescent signal generated by 10 RNAs introduced at known concentrations prior to labelling and detected by complementary reporter (10 replicates of each reporter spotted on the array). Data were generated from representative arrays selected from transcript analyses performed on RNA extracted from control and copper-dosed samples. (A)-(E) represent data from copper exposures of 0, 10, 40, 160 and 480 mg/kg, respectively. The average signal is indicated by closed circles with technical error bars representing the standard error of the measurements. A fitted regression line is shown for the linear portion of the response together with the R 2 value for the fitted line. (PDF 430 KB)
Additional file 11: Graphical representations of relative gene expression against fluorescence intensity from control and copper-exposed samples. Array data were normalised and filtered (as described in Methods) and the log2 of the average fold change (M) plotted against the log2 of the average mean signal intensity (A). (A)-(E) represent data from copper exposure of 0, 10, 40, 160 and 480 mg/kg, respectively. (PDF 5 MB)
Additional file 12: Table of metabolite concentrations as determined by 600 MHz 1H spectroscopy (μmol/mg tissue dry weight). (XLS 36 KB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
Bundy, J.G., Sidhu, J.K., Rana, F. et al. 'Systems toxicology' approach identifies coordinated metabolic responses to copper in a terrestrial non-model invertebrate, the earthworm Lumbricus rubellus . BMC Biol 6, 25 (2008). https://doi.org/10.1186/1741-7007-6-25
- Nuclear Magnetic Resonance
- Nuclear Magnetic Resonance Spectroscopy
- Citric Acid Cycle
- Metabolomic Analysis
- Copper Exposure