Normal operating range (NOR) in Enchytraeus albidus – Transcriptional responses to control conditions

Applied Soil Ecology 85 (2014) 1–10

Contents lists available at ScienceDirect

Applied Soil Ecology journal homepage: www.elsevier.com/locate/apsoil

Normal operating range (NOR) in Enchytraeus albidus – Transcriptional responses to control conditions Sara C. Novais *, Mónica J.B. Amorim Department of Biology and CESAM, University of Aveiro, 3810-193 Aveiro,Portugal

A R T I C L E I N F O

A B S T R A C T

Article history: Received 23 May 2014 Received in revised form 1 August 2014 Accepted 12 August 2014 Available online xxx

Gene expression studies provided a major step forward in ecotoxicology and ecological functional genomics, among other areas. How to classify changes in gene expression to be significant compared to its normal operating range (NOR) requires a number of biological and statistical criteria. Doubts have been casted about when a change in gene expression has ecologically relevant meaning and what is the normal biological variation or NOR. Because genes are highly sensitive and self-regulating pieces, how can we design our experiments to control normal variation? The main aim of the present paper is to discuss such questions, while providing suggestions for improvement and interpretation. We have analysed several gene control data sets of Enchytraeus albidus comparing standard control test conditions, i.e. (1) from experiments performed at different occasions in time, (2) to culture media, (3) to acetone control, (4) from exposure time series (2, 4 and 8 days) and (5) between replicates. We discriminated between some normal variation sources. For instance, the transference of organisms from cultures to standard testing condition suggests a lowered general metabolism. Chemical carriers like acetone seem to cause changes related with soil structure and hypoxia conditions. Organisms have a “time-like” transcription signature when sampled along a time series. The coefficient of variance (CV) between biological replicates was 0.1 < CV < 1.3 for 17% of the genes. Variability will partly be inherent to the species, but most important is to understand the ground for interpretations as discussed. Here, the first attempt for the establishment of a NOR for E. albidus gene expression was made which represents a naturally varying ecological window – this contributes for the discrimination between stress responses and natural ecological variations in soil studies. ã 2014 Elsevier B.V. All rights reserved.

Keywords: Soil Microarray Transcriptional variation Controls Acetone

1. Introduction Transcription studies have the advantage of being fast, sensitive, more chemical specific than population tests (e.g. survival, reproduction), and highly informative of the organisms’ mechanisms of response (van Straalen and Roelofs, 2008). Ecotoxicogenomics helps to identify fitness mechanisms, particularly under stress, or to unravel chemical modes of action and predict effects. In the terrestrial environment microarrays have been successfully applied for a maximum of six species of soil invertebrates – Caenorhabditis elegans, Eisenia fetida, Lumbricus rubellus, Folsomia candida (Brulle et al., 2010), Enchytraeus albidus (Amorim et al., 2011; Novais et al., 2012a) and Enchytraeus crypticus (Castro-Ferreira et al., 2014). Enchytraeids (Oligochaeta) are soil invertebrates, distributed worldwide, that have a key role in the soil organic matter

* Corresponding author. Tel.: +351 234 370790; fax: +351 234 372587. E-mail address: [email protected] (S.C. Novais). http://dx.doi.org/10.1016/j.apsoil.2014.08.005 0929-1393/ ã 2014 Elsevier B.V. All rights reserved.

degradation. Advances in the toxicogenomics field using the enchytraeid E. albidus, a standard ecotoxicology species (ISO, 2004; OECD, 2004), include various chemicals namely metals and pesticides (Novais et al., 2012b,c; Gomes et al., 2014), nanoparticles (Gomes et al., 2012b, 2013), as well as non-chemical factors such as soil properties or exposure time (Gomes et al., 2011; Novais et al., 2012d). Nevertheless, no attention has been given to the range of gene expression variation within control or reference conditions, often absent in microarray results interpretations. De Boer et al. (2011) developed one of the first attempts to obtain a gene normal operating range (NOR) for F. candida under environmental varying factors, i.e. soil types and land-use (dairy farming, agriculture, natural grassland and forest). There are also several other environmental/experimental factors, e.g. changing of culture media or duration of the experiment, which may bias gene expression. Further, even in controlled laboratory experiments, most species will have differences among individuals that are possibly reflected in gene expression. Biological variation, which is in fact one of the pillars of evolution, should of course be considered (quantified) when interpreting transcription results.

2

S.C. Novais, M.J.B. Amorim / Applied Soil Ecology 85 (2014) 1–10

photoperiod 16:8 h light:dark. The organisms are fed twice a week with finely ground and autoclaved rolled oats.

This is even more crucial when testing conditions with low effect on the organisms, e.g. an EC10 (10% effect on reproduction) since the differences in expression can be expected minimal when compared to control and can either be under- or over-estimated in regard to the differences among replicates (biological variation). Hence, a series of unresolved questions persist among scientists, like how ecologically relevant can a change in gene expression be, or how to identify which changes are beyond normal biological variation i.e. will cause an expression in function. To further understand such aspects it is important to identify how much can gene expression vary within organisms (biological replicates). Some authors have established a database pinpointing the most variable genes of (their) species of interest, e.g. for mice up to 3.3% of the gene expression was found varying between organs (Pritchard et al., 2001) and for the bacteria Bacillus subtilis 11% of the genes showed high variation between replicates (Zhou et al., 2011). Moreover, in e.g. the fish Fundulus heteroclitus (Oleksiak et al., 2002) 18% of the genes were found to significantly vary in expression among individuals. On the other hand, other species seem to have very low variation like F. candida (de Boer et al., 2011) with a coefficient of variation in expression per gene between 0.05 and 0.09. Another interesting aspect that is usually underestimated in microarray experiments refers to the effect of carriers (solvents) used for non-water soluble chemicals spiking. Although effects of non-water soluble pesticides have been assessed in E. albidus (Novais et al., 2012b) and compared to the acetone control, the effect of this solvent alone was not assessed in gene expression and should be addressed given the possible interaction effects. Moreover, the question of comparability of results between similar exposures performed at different occasions remains. The main aim of the present paper is to promote critical discussion on these questions and concerns, while providing suggestions which may improve the interpretation of results, helping to discriminate between responses with true positive effects on organisms’ from basal responses to their naturally dynamic environment, this being a false positive (or a true ecological measure). Secondly the aim is to specifically investigate the NOR in E. albidus gene expression to provide an insight basis for the previously described aspect (true vs false positives). To support the discussion we have analysed several gene control data sets of E. albidus comparing standard control test conditions, i.e. (1) from same experiments performed at different occasions in time (e.g. January and October); (2) to culture media; (3) to acetone control; (4) from exposure time series (2, 4 and 8 days) and (5) between replicates.

2.2. Experimental procedure The experimental procedure consists of five parts as related to the objectives, divided in two experiments (Experiments 1 and 2), including to evaluate: (1) the difference between experiments performed separately in time (Experiments 1 and 2); (2) the effect of transferring animals from culture conditions onto standard control test conditions: control (CT) vs culture (CULT) – Experiment 1; (3) the effect of acetone (frequently used as carrier) compared to water: acetone (ACT) vs water control (CT) – Experiment 1; (4) the effect of different exposure times in CT conditions: 2, 4 and 8d – Experiment 2; (5) the variation within CT biological replicates (normal operating range – NOR) – Experiment 1. All comparisons were performed using three replicates. Tests followed the common procedures in the standard guidelines (ISO, 2004; OECD, 2004) with due adaptations. The soil used for all testing was the natural standard soil LUFA 2.2 (Speyer, Germany) which has the following properties: pH (CaCl2) = 5.5; organic matter = 4.4%; texture = 6% clay, 17% silt, 77% sand. To perform the acetone exposure, soil was spiked with the same amount of acetone used as carrier with non-water soluble pesticides in previous experiments (Novais et al., 2012b) and left to evaporate overnight. Before test start, soils were moistened to 50% of the WHCm. Each replicate consisted of a test vessel with 25 g of moist soil and ten adult worms with well developed clitellum (reproductively active) randomly selected from the cultures. Three replicates per treatment were used. The exposure ran at 20  C and at 16:8 h photoperiod. Further, organisms were sampled directly from the cultures (CULT – 10 adult worms with well developed clitellum for each of the three biological replicates) at day 0 to assess the effects of standard control test conditions – CT. The various control exposure conditions are summarised in Table 1. At test end, organisms of each replicate were carefully removed from soil, rinsed in deionised water and snap frozen in liquid nitrogen in RNA later (Ambion, USA) containing criotubes. All samples were stored at 80  C till further analysis. Results from a previous experiment (Novais et al., 2012c) were used for the 2, 4 and 8 days control comparisons – Experiment 2. 2.3. RNA extraction, labelling and hybridizations The following procedures were the same as described in Novais et al. (2012b). In short, total RNA was isolated with the Trizol extraction method followed by a DNAse treatment and further phenol/chloroform extractions for purification purposes. A singlecolour design was used starting from 500 ng of total RNA as input for amplification and labelling with the Low Input Quick Amp Labelling Kit (Agilent Technologies), following the manufacturer protocols. Positive controls were added with the one-colour RNA Spike-In Kit (Agilent Technologies). Purification of the amplified and cyanine 3 labelled cRNA was performed with the RNeasy columns (Qiagen). The 8  15 k custom 60-mer oligonucleotide

2. Material and methods 2.1. Test species E. albidus (Henle, 1837) used in the present study were maintained in laboratory cultures under controlled conditions in a mixture of LUFA 2.2 natural soil and OECD artificial soil (3:1), moistened to 40–60% of the soil maximum water holding capacity (WHCm). Cultures are kept at a temperature of 18  C with a

Table 1 Summary of the different control samples used in Experiments 1 and 2, using Enchytraeus albidus, with the respective sampling times and procedure followed. Samples from Experiment 2 are underlined. Controls (codes) Culture media (CULT) LUFA-Ct (CT) LUFA-Ct solvent (ACT)

Exposure time (days) 0 – –

– 2/2 2

– 4 –

– 8 –

Procedure

Experiment

Direct sampling ISO/OECD ISO/OECD

1 1 and 2 1

S.C. Novais, M.J.B. Amorim / Applied Soil Ecology 85 (2014) 1–10

Agilent microarrays developed for this species (Novais et al., 2012a) were used for the hybridization of cRNA samples. Hybridization was performed using the Gene Expression Hybridization Kit (Agilent Technologies) and each biological replicate was individually hybridized on one array. After the 17 h hybridization at 65  C with a rotation of 10 rpm, microarrays were washed using Gene Expression Wash Buffer Kit (Agilent Technologies) and scanned with the Agilent DNA microarray scanner G2505B. The custom microarray contains 1035 unique probes, randomly built on-slide in 14 replicates, which were designed based on the EST sequencing project for E. albidus presented in the database EnchyBASE (Novais et al., 2012a) (http://bioinformatics. ua.pt/enchybase). 2.4. Analysis of microarray data Fluorescence intensity data was measured with Feature Extraction software (version 10.5.1.1, Agilent Technologies). Reports on the Agilent Spike-in control probes and box plots of each array were used for quality control. Processing of the data and statistical analysis were performed using BRB Array Tools version 4.2.1 Stable Release (http://linus.nci.nih.gov/BRB-ArrayTools. html). Using background subtracted signal, the replicated spots within each array were averaged and the intensities were log2 transformed. Data was then normalized using quantile normalization method. Raw and processed data are available from Gene Expression Omnibus (GEO) at the NCBI website [platform: GPL14928; series: GSE45790]. Statistical class comparisons were performed using two-sample t-test with 95% confidence level for the assessment of differentially expressed genes, between the array groups: (1) CT2d vs CULT and (2) ACT2d vs CT2d – Experiment 1; (3) CT4d vs CT2d and (4) CT8d vs CT2d – Experiment 2. The log2 ratios of expression between those groups were calculated and used for further analysis. Annotation of the differentially expressed genes (p < 0.05) was performed based on their similarity to sequences in the National Centre for Biotechnology Information (NCBI) database as determined by the Basic Local Alignment Search Tool (BLAST) (Altschul et al., 1990), as well as using InterProScan based on InterPro database (Zdobnov and Apweiler, 2001). The sequences were submitted to Blast2GO (Conesa et al., 2005) being compared with peptide sequence databases using BLASTX analysis (with e-value < 10 3) and InterProScan, both of which were taken into consideration for sequence annotations in Blast2GO. GO term enrichment analysis (Alexa et al., 2006) was performed for the differentially expressed genes using the same Blast2GO software (Fisher’s exact test, p < 0.05). Heat map and clustering analyses were performed using MultiExperiment Viewer (MeV version 4.7.3, TIGR). 2.5. Quantitative real-time PCR validations Total RNA (1 mg) from CT and ACT samples was converted into cDNA through a reverse transcription reaction using the SuperScript First-Strand Synthesis System for qPCR (Invitrogen). Platinum SYBR Green qPCR SuperMix-UDG (Invitrogen) was used for amplification on the 7500 Real-Time PCR System (Applied Biosystems). Oligo ExplorerTM (version 1.1.0) software was used to design primer sets for four target genes (Proteasome subunit alpha type 6 – EAC00263, Beta tubulin – EAC00282, MADS flclike protein 2 – EAC00896, Aquaporin 1 – EAC00867; see Table S3, Supplementary material) and one endogenous control gene (myosin alkali light chain 1 – EAC00302). The standard and melting curves were used to determine each primer sets efficiency and specificity, respectively. Primer sequences can be found in EnchyBASE – http://bioinformatics.ua.pt/enchybase/ (Novais et al., 2012a).

3

Quantitative real-time PCR (qPCR) was performed with three biological replicates of each condition (the same used for the microarray Experiment 1), applied in triplicate on a 96-well optical plate (GeneAmp1, Applied Biosystems). Reaction conditions consisted of one initial cycle at 50  C for 2 min, followed by a denaturation step at 95  C for 2 min, 40 cycles at 95  C for 32 s and 1 cycle at 60  C for 1 min. Finally, a dissociation step was made consisting of 15 s at 95  C, 1 min at 60  C, and 15 s at 95  C. A mean normalized expression value was calculated from the obtained Ct values of the test genes with Relative Expression Software Tool (REST-MSC). 3. Results and discussion Data analysis followed a stepwise approach as follows. 3.1. Comparisons between experiments performed at different occasions in time In approach 1, all microarray data from the 2 experiments were normalized together (Experiment 1: CT, ACT, CULT; Experiment 2: CT2d, CT4d, CT8d). After a hierarchical cluster of samples based on each treatment averaged intensities, the two individual experiments were grouped separately (see Fig. S1 – Supplementary material). The normalized intensity values of arrays between experiments were considerably different, revealing differences between experiments rather than differences between treatments. This was not fully surprising given that the experiments were carried out approximately one year apart, using different batches of organisms, chemicals and set up. Of course the advantages of a single-colour design enable all comparisons between treatments upon analysis, but we realise that comparisons should be made only using exposures from the same experiment in such cases. If we think about it, this is not so different from survival/ reproduction toxicity results, where performance in controls varies and results cannot be combined, unless similar absolute values occur between controls of the two experiments. Variability between tests depends on various factors and also on species, e.g. regarding reproduction tests, F. candida offspring numbers in controls can range between 200 and 2000 juveniles, and E. albidus between 50 and 500 juveniles, this under the exact same control test procedures/conditions. In terms of gene expression, to exemplify this variability, we have compared by means of linear regression, controls from Experiment 1 and 2, which were kept in the exact same standard test conditions and sampled after 2d (Fig. 1). With this regression analysis (R2 = 0.95) we can also observe that the log2 fold change for certain genes was up to 6 between the two controls, showing the possible variation in responses between experiments. Such variability of gene expression in E. albidus was higher than in the springtail F. candida (de Boer et al., 2011). de Boer et al. (2011) have addressed the variation of gene expression using LUFA controls from different experiments performed separately in time and found reproducible gene expression signatures, with a maximum log2 fold change of 0.42. There are a number of reasons that could be added to the explanation of differences in variability between these two species, namely the reproduction (parthenogenetic vs sexual), age (synchronized vs random clitellate adults), and number of pooled animals (30 vs 10). In fact, such a natural transcriptional variation is known to represent a very high advantage and value as it can contribute to differences among individuals in populations which enables populations to better adapt to changing environmental conditions (Roelofs et al., 2009; Whitehead and Crawford, 2006). If we consider the results obtained by de Boer et al. (2011) with F. candida in 26 soils, both reproduction and gene expression

4

S.C. Novais, M.J.B. Amorim / Applied Soil Ecology 85 (2014) 1–10

Fig. 1. Comparison of results from Enchytraeus albidus gene expression in control treatments between Experiments 1 and 2. Linear regression between log2 intensity averaged values of control treatments, sampled after 2 days, from Experiments 1 to 2 (CT2d – Experiment 1; CT2d – Experiment 2).

indicated low variation. Such pattern had been verified previously when testing 19 distinct soils (Amorim et al., 2005a) in F. candida reproduction whereas Enchytraeids (E. albidus and E. luxuriosus) showed much higher sensitivity to the same 19 soils, including no reproduction at all and mortality of adults (Amorim et al., 2005b). Moreover, E. albidus are commonly more sensitive to toxicants compared to F. candida and in line with present observations. Such is mostly related with the ecological differences among species. 3.2. Effects of standard control test conditions (in comparison to cultures) and acetone In approach 2, arrays were normalized within each experiment, i.e. to compare treatments that were done per experiment, excluding comparison between experiments and variation from that source. In Experiment 1, statistical class comparisons (two sample ttests, p < 0.05) between treatments (control standard test conditions = CT; cultures = CULT; acetone = ACT) were performed and yield the number of differentially expressed genes (DEG) shown in Fig. 2 (Experiment 1).

Fig. 2. Number of differentially expressed genes in Enchytraeus albidus in treatments from Experiments 1 and 2. Number of significantly up- and downregulated genes (two-sample t-tests, p < 0.05) in Enchytraeus albidus when comparing organisms of: (1) CT–CULT – control experimental conditions vs culture; (2) ACT–CT – acetone vs control conditions; (3) CT4d – 4d vs 2d control experimental conditions; (4) CT8d – 8d vs 2d control experimental conditions.

Expression profiles of the total 282 DEG are represented in the heat map of Fig. 3 for every replicate of the three conditions tested. The hierarchical clustering of samples was performed based on the log2 intensity values for each gene (Fig. 3B). Results show a larger difference between ACT and the other conditions (CULT and CT), a pattern maintained when using averaged values between replicates (Fig. 3C). Based on cluster analysis it can be seen that expression profiles between CT and CULT are similar as supported by the fact that replicates of each treatment are not always clustered together, e.g. CT_2 and CT_3 are more similar to CULT_1 and CULT_2 than to CT_1 (Fig. 3B). Log2 ratios of gene expression between the comparison groups of interest were calculated. To simplify, the treatment that compares CT with CULT was designated CT–CULT (effect of control test media conditions) and comparison between ACT with CT was designated ACT–CT (effect of acetone). The complete list of affected genes with and without significant blast homologies is given in Table S1 (Supplementary material). The total number of transcripts shared between CT–CULT and ACT– CT treatments is presented in the Venn diagram (Fig. 4A). From genes affected by CT–CULT and ACT–CT treatments there were few (24) common (Fig. 4). Further, only 5 out of the 24 common genes had the same direction (up- or downregulation) of response (Table S1, Supplementary material). In both treatments, more down than up-regulated transcripts were observed. From the 24 genes commonly affected, three had homologies with known proteins in public databases and coded for a beta tubulin (EAC00282), a transcription factor from the MADS-box family (MADS flc-like protein 2 – EAC00896) and the aquaporin 1 (EAC00867), all regulated in opposite directions when comparing CT–CULT and ACT–CT. Aquaporins are water channel proteins that provide pathways for water movement across cell membranes (Gomes et al., 2011) and the over expression of these genes have been associated with hypoxia (Kaur et al., 2006). The process of spiking the soil with acetone, including the solvent evaporation and soil drying, has been shown to change the soil structure and reduce the total soil porosity (Thompson et al., 1985). This change in the test media structure which possibly decreased water bioavailability from soil particles could be the origin for the over expression of water channel protein coding genes (confirmed by qPCR – Table S3, Supplementary material). Similarly, the differences in soil porosity between culture and standard test media would explain the opposite change with down-regulation of these genes in CT–CULT, since the organisms are transferred to a soil media with supposedly increased porosity. Gomes et al. (2011) have also observed an over expression of this gene in E. albidus when exposed to soils with a higher clay content. These findings suggest that expression of aquaporin gene can possibly be a marker of changes related with soil physical structure or changed porosity/ aeration. Although CT–CULT induced more genes than ACT–CT, the range of log2 fold inductions were lower ( 2 to +2) compared to the solvent treatment ( 3 to +6). Moreover, in ACT–CT 20% of genes were regulated more than three times higher or lower in comparison to control, whereas in CT–CULT only 2% of the genes reached that level (Table S1, Supplementary material). The results clearly indicate that the organisms’ stress/change level is higher in the soil spiked with acetone than due to a change from culture to the test condition/LUFA soil. The majority of these high-foldchange expression transcripts cannot be discussed given their unknown function. Avoidance tests with E. albidus (unpublished data) showed that 80% of the organisms avoid the soil spiked with the same amount of acetone (10 ml per 100 g of soil) whereas at the reproduction level (Novais et al., 2010) a 30% increased reproduction rate was

S.C. Novais, M.J.B. Amorim / Applied Soil Ecology 85 (2014) 1–10

5

Fig. 3. Heat map and clustering of genes and samples in Enchytraeus albidus from Experiment 1. A – heat map and hierarchical clustering of samples and genes (Pearson’s uncentered) based on the differentially expressed genes (two sample t-test, p < 0.05) in Enchytraeus albidus when comparing organisms: (1) of control experimental conditions (CT) vs culture (CULT); (2) exposed to acetone (ACT) vs control conditions (CT). The log2 intensity values represent the averaged background subtracted fluorescence expression values of genes in each sample. B – detail from the hierarchical cluster of samples shown in (A). C – hierarchical cluster of samples using averaged values between replicates.

Fig. 4. Venn diagrams of differentially expressed genes between treatments in Enchytraeus albidus from Experiments 1 and 2. Venn diagram representing the total number of common and uniquely differentially expressed genes in Enchytraeus albidus between the two statistical class comparisons (two sample t-tests, p < 0.05) performed in Experiment 1 (A): CT–CULT – control experimental conditions vs culture; ACT–CT – acetone vs control conditions; and Experiment 2 (B): CT4d – 4d vs 2d control experimental conditions; CT8d – 8d vs 2d control experimental conditions.

6

S.C. Novais, M.J.B. Amorim / Applied Soil Ecology 85 (2014) 1–10

observed, possibly a hormesis like effect. Hence, such level of gene response towards ACT treatment has a physiological consequence and also ecological meaning. These results highlight the importance of comparing with a solvent control on microarray experiments when carriers are required for spiking procedures (e.g. Novais et al., 2012b), and also the importance of studying the effect of the carrier itself by comparing vs a water control. A gene set enrichment analysis of gene ontology terms (Alexa et al., 2006) was independently applied to the lists of up- and down-regulated transcripts in CT–CULT and ACT–CT treatments to identify the biological processes significantly more abundant in the gene lists. The most specific terms of biological processes significantly affected by either CT–CULT or ACT–CT are listed in Table 2. All significant differential transcripts within each GO term are given in Table S2 (Supplementary material). Based on this enrichment analysis, the transference of the organisms to control conditions of a standard test (CT_CULT) affected several biological pathways such as the repression of processes related with the electron transport system – cellular respiration (GO:0045333), proton transport (GO:0015992) and ATP biosynthetic process (ATP biosynthetic process), or related with cellular division – mitosis (GO:0007067) and spindle organization (GO:0007051). The repression of these processes along with inhibited mechanisms of macromolecules biogenesis (cellular component biogenesis – GO:0044085) gives an indication that the organisms might be somehow lowering their metabolism in general – this could be due to reduced locomotion in the test vessel (reduced size compared to culture box), saving or reallocating energy. This lowered metabolism is in accordance with measurements of the electron transport system activity (Novais and Amorim, 2013) where 4d after transfer from culture to control standard test conditions the organisms are consuming less energy per hour (1.62  0.37 mJ/mg Org) than the organisms in the culture (2.64  1.20 mJ/mg Org). On the other hand, CT–CULT treatment induced biological processes related to cellular defence

responses (GO:0006968), cuticle development (GO:0042335) and superoxide anion generation (GO:0042554), through the upregulation of dual oxidase 2 (EAC00863) and uromodulin precursor (EAC00835) (Table S2, Supplementary material). Dual oxidase 1 belongs to a family of NADPH-oxidases (NOX/DUOX family) that generally regulates production of reactive oxygen species (ROS) (Gomes et al., 2011). In Caenorhabditis elegans, this gene was found to be associated with the production of superoxide to be used in tyrosine cross-linking during the formation of cuticle (Donko et al., 2005). Uromodulin gene has been indicated as a regulator of the water permeability of tissues (Mori et al., 2009; Sedor, 2010). The up-regulation of these two genes related with cuticle development and water permeability can possibly be linked to the improvement in oxygen levels and general gas exchange as a response to changes in soil media. As for the ACT–CT treatment, repressed biological processes resulted from the down-regulation of several b-actin, b-tubulin (confirmed by qPCR – Table S3, Supplementary material) and tropomyosin genes that are involved in processes such as protein polymerization (GO:0051258), striated muscle tissue development (GO:0014706) or microtubule-based movement (GO:0007018). Also the protein degradation process seems to be compromised by the down-regulation of proteasome subunit genes (confirmed by qPCR – Table S3, Supplementary material), among the ones responsible for repressed regulation of cellular protein metabolic processes (GO:0032270) and apoptotic processes (GO:0006915). The role of proteasomes in apoptosis is rather complex, with examples being reported of proteasome inhibitors inducing apoptosis in some cells, whereas in others preventing it (Wójcik, 1999, 2002). Such complexity does not allow inferring about which association exists between acetone and apoptotic processes. The induced biological processes by acetone refer to the up-regulation of the gene coding for aquaporin 1 (Table S2, Supplementary material) involved in water transport (GO:0006833), as previously discussed, and the expression was confirmed by qPCR (Table S3,

Table 2 Enrichment analysis of gene ontology (GO) terms. Significant enriched GO terms (p < 0.05) in Enchytraeus albidus within the significantly down- and up-regulated genes in: (1) CT–CULT – organisms from control experimental conditions in comparison with organisms from the culture; (2) ACT–CT – organisms exposed to acetone in comparison with organisms from control conditions. Only the most specific biological process terms are given. GO ID

GO term definition

Genes differentially expressed by CT–CULT Repressed GO:0015992 Proton transport GO:0035193 Larval central nervous system remodelling GO:0006754 ATP biosynthetic process GO:0007067 Mitosis GO:0007051 Spindle organization GO:0045333 Cellular respiration GO:0044085 Cellular component biogenesis – Induced GO:0006968 Cellular defence response GO:0042335 Cuticle development GO:0042554 Superoxide anion generation – Genes differentially expressed by ACT–CT Repressed GO:0051258 Protein polymerization GO:0014706 Striated muscle tissue development GO:0006915 Apoptotic process GO:0007018 Microtubule-based movement GO:0007517 Muscle organ development GO:0032270 Positive regulation of cellular protein metabolic process – Induced GO:0035626 Juvenile hormone mediated signalling pathway GO:0006833 Water transport a b

Library annotateda

Significantb

p-Value

7 2 5 5 5 13 30

4 2 3 3 3 5 8

0.009 0.021 0.023 0.023 0.023 0.027 0.049

1 1 1

1 1 1

0.024 0.024 0.024

4 4 11 5 5 7

2 2 3 2 2 2

0.014 0.014 0.015 0.023 0.023 0.045

1 1

1 1

0.017 0.017

The number of annotated ESTs present in the library within the identified GO terms is indicated for comparison. The number of annotated differentially expressed ESTs within each GO term.

S.C. Novais, M.J.B. Amorim / Applied Soil Ecology 85 (2014) 1–10

Supplementary material). The hypoxia hypothesis is also supported by other evidences that resemble mechanisms known to occur under oxygen deprivation (Hochachka et al., 1996) such as the inhibition of protein synthesis seen by the down regulation of several ribosomal proteins or suppression of the proteasomal protein degradation (confirmed by qPCR – Table S3, Supplementary material) (Nota et al., 2008) also verified in the acetone exposure (Table S1, Supplementary material).

7

for significant alterations in gene expression, with log2 fold changes up to 7. Changes in the organisms’ physiological levels in control along time, had already been verified in studies where oxidative stress and neuronal biochemical parameters were measured (Gomes et al., 2012a; Novais et al., 2011, 2014). Similarly, this also occurs in other species e.g. in oysters Crassostrea virginica (Conners and Ringwood, 2000), earthworm E. fetida (Saint-Denis et al., 1999) or mussels Mytilus galloprovincialis (Vlahogianni and Valavanidis, 2007).

3.3. Quantitative real-time PCR 3.5. Biological variation in gene expression To validate the microarray results, a selection of four genes differentially expressed after the acetone treatment were quantified with qPCR. The selected genes coded for proteasome subunit alpha type 6 (EAC00263), beta tubulin (EAC00282), MADS flc-like protein 2 (EAC00896) and aquaporin 1 (EAC00867), using myosin alkali light chain 1 as housekeeping (EAC00302). The housekeeping gene was selected from a set of 3 tested genes based on the least variation in expression levels between water and acetone controls in the microarray experiment. A significant correlation of 0.903 (Pearson’s correlation, p = 0.0966) was obtained between microarray and qPCR platforms. All response patterns were confirmed with only small variations in fold change magnitudes, with the direction of expression response (up- or down-regulation) being always the same (Table S3, Supplementary material). 3.4. Effects of exposure time in control test conditions In Experiment 2, microarray results were analysed from organisms exposed in control standard test conditions exposed for 2, 4 and 8 days (Novais et al., 2012c), comparing day 2 with the other. Statistical class comparison (two sample t-tests, p < 0.05) between CT4d–CT2d and CT8d–CT2d resulted in the differentially expressed genes observed in Fig. 2 (Experiment 2). The complete list of significant altered genes is given in Table S4 (Supplementary material). The total number of significant transcripts shared between CT4d and CT8d treatments, when compared to CT2d, is presented in the Venn diagram of Fig. 4B. From the total DEGs only 18 had altered expressions both in CT4d and CT8d when compared to CT2d (Fig. 4B). In these commonly affected genes the direction of response (up- or down-regulation) was always the same independently of the sampling time (Table S4 – Supplementary material). From the 5 commonly affected genes with known homology, 3 were down-regulated and coded for a protease, a heat shock protein and a calcium binding protein, while 2 genes were upregulated and coded for a transcription factor and a translation elongation factor (Table S4 – Supplementary material). Maintaining the organisms for 4d in standard control testing conditions (CT4d) repressed the expression of genes mainly related with ubiquitin-dependent protein catabolic processes, ATP synthesis, DNA repair and response to stimulus, while enhancing expression of genes mainly involved in carbohydrate metabolism (Table S4 – Supplementary material). These results suggest that at 4d the organisms are more adjusted to the new test conditions than at day 2, as they seem to be lowering their response mechanisms. At 8d fewer differences in gene response were observed compared to 2d showing an increased response again after 4d which indicates that such variations may also be part of the normal organisms’ gene changes status. The repressed genes were mainly related with ATP synthesis and translation. As for the induced genes, these were mostly involved in the regulation of transcription (Table S4 – Supplementary material). To summarise, the different sampling periods between 2 and 8 days in standard control test conditions were responsible alone

One of the major concerns when interpreting microarray experiments is to discriminate differentially expressed genes due to a treatment from normal biological variation. According to Zhou et al. (2011) the coefficient of variance (CV) between genes from one dataset to another is useful to assess variation, even when there are large differences in their means. Based on the biological variation in Bacillus subtilis gene expression, Zhou et al. (2011) established a CV cut-off of 0.1, defining genes with CV > 0.1 as highly variable. Based on this we used the CT samples from Experiment 1 to calculate the CV between replicates. The CV values of all genes between the three biological replicates were calculated as the ratio of the standard deviation to the mean. Overall, the CV ranged from 0 to 1.3 with 83% of the genes showing a CV < 0.1 (Fig. 5). If we take CV > 0.1 as an arbitrary value, from the 1035 genes present in the array, 179 (17%) had high variance levels (CV > 0.1). From these genes, 19 (1.8%) are genes varying more than 50% and 3 have CV higher than 1 (Fig. 5). This relatively high variance levels in gene expression between individuals are not very surprising taking into account that these represent the normal biological variation when analysing whole organism RNA. Compared to tissue-specific RNA, whole organism variation is more expressive. Such gene expression variations can thus be considered within the normal operating range (NOR) for E. albidus. The list of these genes with higher variability (and known homology) is presented in Table 3. Among these genes, the most representative are related with cell motility processes e.g. actins and myosins, as well as with transcription with the presence of several transcription factors. There are also several varying genes involved in proteolysis e.g. carboxypeptidases, fibrinolytic enzyme, cathepsin l-like, proteasome subunit alpha type-4. Some of these proteolytic genes are also involved in anti-inflammatory or defence mechanisms which,

Fig. 5. Average variance in Enchytraeus albidus gene expression, per gene, calculated from the standard control replicates. The genes are ordered from the ones with low variance to high variance. The dashed line marks the 10% of variation limit.

8

S.C. Novais, M.J.B. Amorim / Applied Soil Ecology 85 (2014) 1–10 Table 3 List of the most varying genes in control conditions. List of genes of Enchytraeus albidus microarray in which expression greatly varies between replicates of control organisms [coefficient of variance (CV) > 0.1]. Only genes with known homology are shown. Cluster IDa

Blast homology

CV

EAC01014 EAC00432 EAC01069 EAC00254 EAC00358 EAC00910 EAC00112 EAC00484 EAC00265 EAC01056 EAC01009 EAC00323 EAC00492 EAC00267 EAC00311 EAC00657 EAC00671 EAC00829 EAC01024 EAC00481 EAC00097 EAC00710 EAC00180 EAC00397 EAC00516 EAC00290 EAC00933 EAC00962 EAC01003 EAC01049 EAC01061 EAC01083 EAC01116 EAC00463 EAC00433 EAC00656 EAC00732 EAC00154 EAC01030 EAC01093 EAC01105 EAC01119 EAC00178 EAC00351 EAC00987 EAC00179 EAC00313 EAC00927 EAC00214 EAC00809 EAC00880 EAC00901 EAC00508 EAC00549 EAC01057 EAC00373 EAC00681 EAC00404 EAC01016 EAC00847 EAC01050 EAC00513 EAC00325 EAC00284 EAC00286

40 s ribosomal protein Actin partial Actin partial Actin-partial Actin-partial Adenosine kinase Alpha skeletal muscle-like isoform 4 Aspartate cytoplasmic Beta-glucanase Beta-actin Carboxypeptidase Cathepsin l-like Chorion b-zip transcription factor Chymotrypsinogen 2-like Cysteine-rich protein 1 Cytochrome b-c1 complex subunit 10-like Dpy-30 homolog Elongation factor 1 alpha Fibrinolytic enzyme Heat shock protein 60 Heat shock protein gp96 Hypothetical protein TGME49_026570 [Toxoplasma gondii ME49] Immunlectin a Keratinocyte-associated protein 2 Kif21a protein Lysozyme MADS flc-like protein 2 MADS flc-like protein 2 MADS flc-like protein 2 MADS flc-like protein 2 MADS flc-like protein 2 MADS flc-like protein 2 MADS flc-like protein 2 Mam and ldl-receptor class a domain-containing protein c10orf112-like Methyltransferase mb3374 Mitochondrial ribosomal protein l48 Monosaccharide-transporting atpase Myosin heavy chain Myosin heavy chain Myosin heavy chain Myosin heavy chain Myosin heavy chain Neuronal acetylcholine receptor subunit beta-2 precursor Neuronal calcium sensor 2 Non-structural polyprotein Plastin-1 Procollagen type vi alpha 4 Proteasome subunit alpha type-4 Protein Protein Protein yippee-like 5 Puroindoline b protein Ran binding protein 1 Replication factor a protein 3 domain-containing protein Ribonuclease uk114 Selenium binding protein 1 Short branched chain specific acyl- mitochondrial-like Short-chain collagen c4- partial Solute carrier family 22 (organic cation transporter) member partial Splicing factor 3a subunit 3 Tektin 2 Waprin-phi1-like isoform 2 Xk-related protein 6-like Zinc carboxypeptidase a 1-like Zinc carboxypeptidase a 1-like

0.573617 0.17314 0.642924 0.112318 0.159172 0.370498 0.248055 0.343675 0.49539 0.110012 0.144632 0.269702 0.177713 0.619052 0.192786 0.193423 0.523849 0.153314 0.152171 0.792422 0.100765 0.430356 0.554655 0.311622 0.371754 0.104175 0.114475 0.116795 0.155812 0.170853 0.226904 0.238756 0.275994 0.350037 0.356216 0.318379 0.473379 0.118249 0.18092 0.19698 0.200431 0.221426 0.13305 0.359536 0.128043 0.221899 0.356671 0.101221 0.266175 0.349127 0.105324 0.382635 0.105967 0.440606 0.115418 0.207511 0.529729 0.196915 0.272171 0.112163 0.162089 0.101961 0.149249 0.12448 0.126144

a

Cluster ID as can be found in EnchyBASE (http://bioinformatics.ua.pt/enchybase/)

together with the presence of other genes like lysozyme, immunlectin or chaperones, point to different immune status of the organisms at the time of sampling. This same observation that differences in immune status may be one of the largest sources of gene expression variation was made by Pritchard et al. (2001) when

studying the normal variance in mouse gene expression under control situations. The high variation in expression level of these genes under normal conditions increases their chance of being scanned as false positives/negatives when testing a stressor. The present results

S.C. Novais, M.J.B. Amorim / Applied Soil Ecology 85 (2014) 1–10

require further consolidation but represent a first indicative basis. The fact that 3 replicates were used for this particular step could be considered smaller than optimal, although it actually means that data comes from 30 animals (each replicate has a pool of 10), often higher number of individuals than the one used for similar studies with e.g. rats or cell lines. Hence, for example, if reviewing the gene expression from ACT– CT, it is possible to observe that 34 of the 130 DEGs were genes having a CV > 0.1 in CT (Table S1, Supplementary material; Table 3). From that list of 34 genes, no gene with known homology was selected after gene enrichment analysis and therefore was not considered in the discussion related with acetone effects. The same for CT–CULT, where 6 of the 176 DEGs could be false positives based on the CV but again these were not included after enrichment analysis and therefore were not considered in the discussion. The housekeeping gene used for qPCR validations in Experiment 1 had a CV of 0.005 in CT. In relation to Experiment 2, if considering the list of genes with CV > 0.1 (Table 3), these represent 19% of the total 161 DEGs due to exposure time in CT (Table S4 – Supplementary material), which seems to indicate the presence of more than biological variation alone. The organisms seem to have “time” in the transcription signatures even if all other factors are kept constant within 0–8 days. 4. Conclusions Overall, the present gene analysis was targeted to measure variations in controls including changes in exposure time and test media, i.e., ecological (non-stress) factors, having provided for the first time an insight on levels under such factors. There are many studies and species with larger annotation coverage available from which this information could be extracted and analysed, providing an important comparison line or NOR guidance. It was clear that exposure time and the use of carriers should be regarded if included in the test design. Results from experiments performed at different occasions in laboratory can have differences larger than the treatments in each experiment. Overall, organisms showed a gene expression CV < 1.4 which represents the first guidance NOR in E. albidus transcription. This information should be considered in future microarray experiments to support/improve the interpretation of levels of significant variation and hence improve the ecological relevancy. Further, other techniques within next generation sequencing, e.g. RNA-Seq can highly advance the screening and refining of the NOR for nonmodel organisms given the quantitative real time of multiple samples. Retrieving information from the control samples about the most highly variable (and stable) genes can help filtering off genes that are true from false positives within normal biological variation. The choice of the corresponding solvent control is a key aspect for a direct comparison and the effect of the carrier alone should be studied. We recommend that to answer the question of how ecologically relevant can a change in gene expression be or how to quantify at which level changes are beyond normal biological variation will benefit from (a) establishing a more solid NOR (species specific) to use as reference, and (b) anchoring gene expression level studies with responses at other levels, e.g. effects on oxidative stress markers, reproduction, survival or energy measurements. Acknowledgements This study was sponsored by the funding FEDER-Fundo Europeu de Desenvolvimento Regional through COMPETE-Programa Operacional Factores de Competitividade, and by National funding through FCT-Fundação para a Ciência e Tecnologia, within the

9

research project NANOkA FCOMP-01-0124-FEDER-008944 (Refa. FCT PTDC/BIA-BEC/103716/2008) and through a postdoctoral grant to Sara Novais (SFRH/BPD/80337/2011). Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.apsoil.2014.08.005. References Alexa, A., Rahnenfuhrer, J., Lengauer, T., 2006. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics 22, 1600–1607. Altschul, S., Gish, W., Miller, W., Meyers, E., Lipman, D., 1990. Basic local alignment search tool. J. Mol. Biol. 215, 403–410. Amorim, M.J.B., Nogueira, A.J.A., Soares, A.M.V.M., Römbke, J., Scheffczyk, A., 2005a. Effects of different soil types on the Collembolans Folsomia candida and Hypogastrura assimilis using the herbicide Phenmedipham. Arch. Environ. Contam. Toxicol. 49, 343–352. Amorim, M.J.B., Rombke, J., Scheffczyk, A., Soares, A.M.V.M., 2005b. Effect of different soil types on the enchytraeids Enchytraeus albidus and Enchytraeus luxuriosus using the herbicide Phenmedipham. Chemosphere 61, 1102–1114. Amorim, M.J.B., Novais, S.C., Van der Ven, K., Vandenbrouck, T., Soares, A.M.V.M., De Coen, W., 2011. Development of a microarray for Enchytraeus albidus (Oligochaeta): preliminary tool with diverse applications. Environ. Toxicol. Chem. 30, 1395–1402. Brulle, F., Morgan, A.J., Cocquerelle, C., Vandenbulcke, F., 2010. Transcriptomic underpinning of toxicant-mediated physiological function alterations in three terrestrial invertebrate taxa: a review. Environ. Pollut. 158, 2793–2808. Castro-Ferreira, M.P., de Boer, T.E., Colbourne, J.K., Vooijs, R., van Gestel, C.A.M., van Straalen, N.M., Amorim, M.J.B., Roelofs, D., 2014. Transcriptome assembly and microarray construction for Enchytraeus crypticus, a model oligochaete to assess stress response mechanisms derived from soil conditions. BMC Genomics 15, 302. Conesa, A., Gotz, S., Garcia-Gomez, J.M., Terol, J., Talon, M., Robles, M., 2005. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 21, 3674–3676. Conners, D.E., Ringwood, A.H., 2000. Effects of glutathione depletion on copper cytotoxicity in oysters (Crassostrea virginica). Aquat. Toxicol. 50, 341–349. de Boer, T.E., Birlutiu, A., Bochdanovits, Z., Timmermans, M.J.T.N., Dijkstra, T.M.H., Van Straalen, N.M., Ylstra, B., Roelofs, D., 2011. Transcriptional plasticity of a soil arthropod across different ecological conditions. Mol. Ecol. 20, 1144–1154. Donko, A., Peterfi, Z., Sum, A., Leto, T., Geiszt, M., 2005. Dual oxidases. Philos. Trans. Royal Soc. B Biol. Sci. 360, 2301–2308. Gomes, S.I.L., Novais, S.C., Soares, A.M.V.M., Amorim, M.J.B., 2011. Effects of soil properties and time of exposure on gene expression of Enchytraeus albidus (Oligochaeta). Soil Biol. Biochem. 43, 2078–2084. Gomes, S.I.L., Novais, S.C., Gravato, C., Guilhermino, L., Scott-Fordsmand, J.J., Soares, A.M.V.M., Amorim, M.J.B., 2012a. Effect of Cu-nanoparticles versus one Cu-salt: analysis of stress biomarkers response in Enchytraeus albidus (Oligochaeta). Nanotoxicology 6, 134–143. Gomes, S.I.L., Novais, S.C., Scott-Fordsmand, J.J., de Coen, W., Soares, A.M.V.M., Amorim, M.J.B., 2012b. Effect of Cu-nanoparticles versus Cu-salt in Enchytraeus albidus (Oligochaeta): differential gene expression through microarray analysis. Comp. Biochem. Physiol. C Toxicol. Pharmacol. 155, 219–227. Gomes, S.I.L., Soares, A.M.V.M., Scott-Fordsmand, J.J., Amorim, M.J.B., 2013. Mechanisms of response to silver nanoparticles on Enchytraeus albidus (Oligochaeta): survival, reproduction and gene expression profile. J. Hazard. Mater. 254–255, 336–344. Gomes, S.I.L., Scott-Fordsmand, J., Amorim, M.J.B., 2014. Profiling transcriptomic response of Enchytraeus albidus to Cu and Ni: comparison with Cd and Zn. Environ. Pollut. 186, 75–82. Henle, F.G.J., 1837. Ueber Enchytraeus, eine neue Anneliden-Gattung. Arch. Anat. Physiol. Med. 74–90. Hochachka, P.W., Buck, L.T., Doll, C.J., Land, S.C., 1996. Unifying theory of hypoxia tolerance: molecular metabolic defense and rescue mechanisms for surviving oxygen lack. Proc. Natl. Acad. Sci. U. S. A. 93, 9493–9498. ISO, 2004. Soil Quality – Effects of Pollutants on Enchytraeidae (Enchytraeus sp.). Determination of Effects on Reproduction and Survival, Guideline no. 16387. ISO (International Organization for Standardization), Geneve, Switzerland. Kaur, C., Sivakumar, V., Zhang, Y., Ling, E.A., 2006. Hypoxia-induced astrocytic reaction and increased vascular permeability in the rat cerebellum. Glia 54, 826–839. Mori, T., Kawachi, H., Imai, C., Sugiyama, M., Kurata, Y., Kishida, O., Nishimura, K., 2009. Identification of a novel uromodulin-like gene related to predatorinduced bulgy morph in Anuran tadpoles by functional microarray analysis. PLoS One 4. Nota, B., Timmermans, M.J.T.N., Franken, C., Montagne-Wajer, K., Marien, J., De Boer, M.E., De Boer, T.E., Ylstra, B., Van Straalen, N.M., Roelofs, D., 2008. Gene expression analysis of Collembola in cadmium containing soil. Environ. Sci. Technol. 42, 8152–8157.

10

S.C. Novais, M.J.B. Amorim / Applied Soil Ecology 85 (2014) 1–10

Novais, S.C., Amorim, M.J.B., 2013. Changes in cellular energy allocation in Enchytraeus albidus when exposed to dimethoate, atrazine and carbendazim. Environ. Toxicol. Chem. 32, 2800–2807. Novais, S.C., Soares, A.M.V.M., Amorim, M.J.B., 2010. Can avoidance in Enchytraeus albidus be used as a screening parameter for pesticides testing? Chemosphere 79, 233–237. Novais, S.C., Gomes, S.I.L., Gravato, C., Guilhermino, L., de Coen, W., Soares, A.M.V.M., Amorim, M.J.B., 2011. Reproduction and biochemical responses in Enchytraeus albidus (Oligochaeta) to zinc or cadmium exposures. Environ. Pollut. 159, 1836– 1843. Novais, S.C., Arrais, J., Lopes, P., Vandenbrouck, T., De Coen, W., Roelofs, D., Soares, A. M.V.M., Amorim, M.J.B., 2012a. Enchytraeus albidus microarray: enrichment, design, annotation and database (EnchyBASE). PLoS One 7, e34266. Novais, S.C., De Coen, W., Amorim, M.J.B., 2012b. Gene expression responses linked to reproduction effect concentrations (EC 10, 20, 50, 90) of dimethoate, atrazine and carbendazim, in Enchytraeus albidus. PLoS One 7, e36068. Novais, S.C., De Coen, W., Amorim, M.J.B., 2012c. Transcriptional responses in Enchytraeus albidus: comparison between cadmium and zinc exposure and linkage to reproduction effects. Environ. Toxicol. Chem. 31, 2289–2299. Novais, S.C., Howcroft, C.F., Carreto, L., Pereira, P.M., Santos, M.A.S., De Coen, W., Soares, A.M.V.M., Amorim, M.J.B., 2012d. Differential gene expression analysis in Enchytraeus albidus exposed to natural and chemical stressors at different exposure periods. Ecotoxicology 21, 213–224. Novais, S.C., Gomes, N.C., Soares, A.M.V.M., Amorim, M.J.B., 2014. Antioxidant and neurotoxicity markers in the model organism Enchytraeus albidus (Oligochaeta): mechanisms of response to atrazine, dimethoate and carbendazim. Ecotoxicology 23, 1220–1233. OECD, 2004. Guidelines for the Testing of Chemicals. Enchytraeid Reproduction Test, No. 220. OECD (Organization for Economic Cooperation and Development), Paris, France.

Oleksiak, M.F., Churchill, G.A., Crawford, D.L., 2002. Variation in gene expression within and among natural populations. Nat. Genet. 32, 261. Pritchard, C.C., Hsu, L., Delrow, J., Nelson, P.S., 2001. Project normal: defining normal variance in mouse gene expression. Proc. Natl. Acad. Sci. 98, 13266–13271. Roelofs, D., Janssens, T.K.S., Timmermans, M.J.T.N., Nota, B., Marien, J., Bochdanovits, Z., Ylstra, B., Van Straalen, N.M., 2009. Adaptive differences in gene expression associated with heavy metal tolerance in the soil arthropod Orchesella cincta. Mol. Ecol. 18, 3227–3239. Saint-Denis, M., Narbonne, J.F., Arnaud, C., Thybaud, E., Ribera, D., 1999. Biochemical responses of the earthworm Eisenia fetida andrei exposed to contaminated artificial soil: effects of benzo(a) pyrene. Soil Biol. Biochem. 31, 1837–1846. Sedor, J.R., 2010. Uromodulin and translational medicine: will the SNPs bring zip to clinical practice? J. Am. Soc. Nephrol. 21, 204–206. Thompson, M.L., McBride, J.F., Horton, R., 1985. Effects of drying treatments on porosity of soil materials. Soil Sci. Soc. Am. J. 49, 1360–1364. van Straalen, N., Roelofs, D., 2008. Genomics technology for assessing soil pollution. J. Biol. 7, 19. Vlahogianni, T.H., Valavanidis, A., 2007. Heavy-metal effects on lipid peroxidation and antioxidant defence enzymes in mussels Mytilus galloprovincialis. Chem. Ecol. 23, 361–371. Whitehead, A., Crawford, D.L., 2006. Variation within and among species in gene expression: raw material for evolution. Mol. Ecol. 15, 1197–1211. Wójcik, C., 1999. Proteasomes in apoptosis: villains or guardians? Cell. Mol. Life Sci. 56, 908–917. Wójcik, C., 2002. Regulation of apoptosis by the ubiquitin and proteasome pathway. J. Cell. Mol. Med. 6, 25–48. Zdobnov, E.M., Apweiler, R., 2001. InterProScan – an integration platform for the signature-recognition methods in InterPro. Bioinformatics 17, 847–848. Zhou, Y., Yu, W.-B., Ye, B.-C., 2011. Variation of gene expression in Bacillus subtilis samples of fermentation replicates. Bioprocess Biosyst. Eng. 34, 569–579.