Expression of target and reference genes in Daphnia magna exposed to ibuprofen

Background Transcriptomic techniques are now being applied in ecotoxicology and toxicology to measure the impact of stressors and develop understanding of mechanisms of toxicity. Microarray technology in particular offers the potential to measure thousands of gene responses simultaneously. However, it is important that microarrays responses should be validated, at least initially, using real-time quantitative polymerase chain reaction (QPCR). The accurate measurement of target gene expression requires normalisation to an invariant internal control e.g., total RNA or reference genes. Reference genes are preferable, as they control for variation inherent in the cDNA synthesis and PCR. However, reference gene expression can vary between tissues and experimental conditions, which makes it crucial to validate them prior to application. Results We evaluated 10 candidate reference genes for QPCR in Daphnia magna following a 24 h exposure to the non-steroidal anti-inflammatory drug (NSAID) ibuprofen (IB) at 0, 20, 40 and 80 mg IB l-1. Six of the 10 candidates appeared suitable for use as reference genes. As a robust approach, we used a combination normalisation factor (NF), calculated using the geNorm application, based on the geometric mean of three selected reference genes: glyceraldehyde-3-phosphate dehydrogenase, ubiquitin conjugating enzyme and actin. The effects of normalisation are illustrated using as target gene leukotriene B4 12-hydroxydehydrogenase (Ltb4dh), which was up-regulated following 24 h exposure to 63–81 mg IB l-1. Conclusions As anticipated, use of the NF clarified the response of Ltb4dh in daphnids exposed to sublethal levels of ibuprofen. Our findings emphasise the importance in toxicogenomics of finding and applying invariant internal QPCR control(s) relevant to the study conditions.


Background
Toxicogenomics is an evolving discipline investigating stressor impact at the genome level [1]. Microarrays are state-of-the-art tools for global gene expression profiling of the response of an organism in a particular biological context. However, expression levels of key genes responding on the microarray need to be validated [2]. The preferred technique of validation is by real-time quantitative polymerase chain reaction (QPCR) [3,4]. Although QPCR is a fairly reliable technique, amplification can vary depending on factors such as RNA integrity, reverse transcriptase (RT) efficiencies, sample-to-sample variations in amplification efficiency, and variation in cDNA sample loading. Using the same sample sizes, assessing RNA integrity and equalising RNA concentrations prior to RT are some of the basic normalisation steps in QPCR [5]. However, normalisation to some internal control is essential for accurate QPCR to balance sample-to-sample variations within the RT and PCR reactions. Internal control is generally achieved using reference genes, also known as housekeeping genes.
Some of the most commonly used reference genes, such as β-actin and glyceraldehyde-3-phosphate dehydrogenase (GAPDH), are historical carryovers from northern blotting, a "predecessor" of QPCR [5]. It is good practice to rigorously validate the suitability of these reference genes under the specific experimental condition to which they are applied [5][6][7]. Total RNA and rRNA have also been suggested as internal controls, although again differences between tissue, individuals and experimental conditions apply [8]. Moreover, using RNA as an internal control has the disadvantage of not controlling for variation inherent in the RT and PCR reactions [5]. None of these genes or RNA products are completely invariable, and thus each needs validation before use [5][6][7]. A further problem of using only one reference gene for normalisation, is that it theoretically should have been normalised itself prior to target gene normalisation [9]. Various methods have been proposed to overcome these problems. One involves the use of a normalisation factor (NF) based on the geometric mean of multiple carefully selected reference genes. This NF can be calculated by geNorm (see Vandesompele et al., 2002 [10] and the geNorm website [11]) or other freely available Excel (Microsoft) based software applications, e.g. BestKeeper [12] (see Pfaffl et al., 2004 [9]). Vandesompele et al., 2002 [10] recommend that NFs should be estimated from a minimum of three reference genes. The advantage of this approach is that it allows for the cumulative error of the entire process from RNA extraction to the QPCR.
Here we evaluate 10 candidate reference genes for QPCR in Daphnia magna following a 24 h exposure to the nonsteroidal anti-inflammatory drug (NSAID) ibuprofen (IB); 20, 40 and 80 mg IB l -1 . This concentration range represents a sublethal exposure based on our earlier work showing that the effect concentration needed to immobilise 50% (EC 50 ) of the exposed individuals was 107.7 mg IB l -1 following 48 h [13]. The tested candidate reference genes were chosen to cover a range of gene ontologies (Table 1, Methods), but mostly representing what may be considered as classical reference genes such as actin and GAPDH [6].
We use the geNorm algorithm [10] to estimate the variability of the reference genes, and to discover an optimal normalisation factor (NF), based on the geometric mean of three of them. geNorm estimates reference gene variability by calculating an expression level ratio for any two candidate reference genes. The variation between the expression ratios is then used as an inverse measure to estimate the variability of the analysed reference genes (please see Methods or Vandesompele et al., (2002) [10] for further information). The approach is illustrated using as target gene the D. magna ortholog of leukotriene B4 12hydroxydehydrogenase (Ltb4dh), which was identified as being up-regulated in a suppressive subtractive hybridisation following 24 h exposure to 63-81 mg IB l -1 [13]. Ltb4dh is of special interest because NSAIDs are known to inhibit the mammalian biosynthesis of various eicosanoids (e.g. prostaglandins, epoxyeicosatrienoic acids and leukotriene) that play important regulatory and signalling functions, for instance regulation of ion flux [14,15]. As anticipated, use of the NF clarified the response of Ltb4dh to ibuprofen.

Water chemistry and quantification of ibuprofen
The chemical parameters pH and conductivity remained stable throughout the exposure, ranging between 7.8-7.9 and 422-452 µS cm -1 , respectively. Water temperatures ranged from 20.8 to 21.0°C. Ibuprofen-sodium increased measured pH by approximately 9% and conductivity by approximately 7% in the highest treatment (80 mg IB l -1 ) compared with the control (data not shown).
Quantification of IB revealed that nominal concentrations were within ± 10% of the measured concentration; except for 20 mg IB l -1 replicate three that was 10.6% lower than the nominal (data not shown). No mortality was observed in D. magna in any of the treatments following exposure.

Relative expression and variability of candidate reference genes
Relative gene expression was estimated using DART-PCR (see Methods). There was a significant (p < 0.05) downregulation of Atb and Cyp at 80 mg IB l -1 compared with the control ( Table 2). The eight candidate reference genes unaffected by the IB treatment were analysed by geNorm (see Methods) ranking the least variable genes as; UBC = GAPDH <Act <WARS <SDH <TBP <18S <28S (Table 2). Subsequently, geNorm was applied to estimate seven normalisation factors,NF 2 to NF 8 , based on the geometric mean of the relative expression of the included reference genes. The geNorm algorithm allocates the two least vari- a Several of the candidate reference genes are classified within multiple biological processes (P), molecular functions (F) or associated with one or more cellular components (C). A biological function, accompanied by a GO number, has been listed here for each gene (for more information visit http://www.geneontology.org/). b Target gene also known as 15-oxoprostaglandin 13-reductase or more generally as NADP(H)-dependent alkenal/one oxidoreductase able genes, here GAPDH and UBC, to NF 2 . Then the sequential NF 3 was based on the genes from NF 2 and the third least variable gene Act, and so on.

Variability of normalisation factors
Pairwise comparisons of sequential normalisation factors computed by geNorm revealed a similar level of variability between NF 2 vs. NF 3 , NF 3 vs. NF 4 , NF 4 vs. NF 5 and NF 5 vs. NF 6 ( Fig. 1). This suggests that the use of either of these NFs as internal control would be valid. But pairwise comparisons of NF 6 vs. NF 7 and NF 7 vs. NF 8 , (28S (NF 7 ) and 18S (NF 8 )), more than doubled the variability (Fig. 1). Vandesompele et al., 2002 [10] recommend using as few reference genes as feasible, but a minimum of three. Thus, although it gave a similar variability to other combina-tions, we chose NF 3 , (GAPDH, UBC and Act) since it gave an optimal NF with a low level of variation using the lowest possible number of reference genes.

Relative expression of reference genes following normalisation
There is a circular problem with normalising a target gene with reference genes that should be normalised themselves prior to target gene normalisation [9]. Thus, in addition to normalising the target gene, we also normalised all 10 reference genes by NF 3 . Atb had the same significant (p < 0.05) differences as prior to normalisation, and should therefore be regarded as a target gene. None of the other reference genes, including Cyp, were significantly (p > 0.05) affected by IB following normalisation (data not shown).

Relative expression of the target gene Ltb4dh following normalisation
Without normalisation, the target gene Ltb4dh was significantly (p < 0.05) up-regulated at 80 mg IB l -1 compared with any of the other treatment groups having a 3.24-fold higher expression than the control (Fig. 2). The relative expression of Ltb4dh followed a significant (p < 0.05) dose-dependent relationship (data not shown). Normalisation with NF 3 produced a similar expression (Fig. 2). The relative expression of Ltb4dh was normalised individually by each of the 10 candidate reference genes and by NF 3 to assess how their variability influenced target gene expression. Overall, the mean fold difference in relative expression between 80 mg IB l -1 and the control was 3.34 ± 0.85 SD ranging from 1.66 (18S) to 5.15 (Atb) (Fig. 3). Target gene expression normalised to any of the reference genes was compared to normalisation with NF 3 revealing that there was a significant (p < 0.05) difference within the Pairwise variation of sequential normalisation factors (NF n vs. NF n+1 ) estimated by geNorm Figure 1 Pairwise variation of sequential normalisation factors (NF n vs. NF n+1 ) estimated by geNorm. NF 2 is based on the geometric mean of GAPDH and UBC; NF 3 is NF 2 and Act; NF 4 is NF 3 and WARS; NF 5 is NF 4 and SDH; NF 6 is NF 5 and TBP; NF 7 is NF 6 and 28S; NF 8 is NF 7 and 18S (gene name abbreviations as Table 1). 80 mg IB l -1 treatment between normalisation of Ltb4dh to Atb and Cyp respectively, and Ltb4dh normalised to NF 3 (Fig. 3). Although not significant (p > 0.05), normalisation of Ltb4dh to either 18S or 28S caused a large increase in the variation of target gene expression in any of the IB treatments (Fig. 3). Both rRNA genes reduced significant differences observed in Ltb4dh, with or without normalisation to NF 3 (data not shown).

Discussion
We investigated the variability of 10 candidate reference genes in D. magna following a 24 h exposure to IB to discover the least variable internal control(s) for QPCR normalisation. A comparison of the reference genes, using the geNorm software, ranked the least variable genes as; UBC = GAPDH <Act <WARS <SDH <TBP <18S <28S. Atb and Cyp were not included in the geNorm analysis as they were down-regulated by IB. Furthermore, geNorm identified the optimal normalisation factor as NF 3 based on the geometric mean expression of UBC, GAPDH and Act. This NF was based on the lowest recommended number of reference genes with the lowest level of variation [10]. The response to IB of an example target gene, Ltb4dh, was little changed by application of NF 3 (Fig. 2), although at the highest concentration of IB, there was a slight increase in the response, and a diminution of the variation between replicates. Such reduction in the variation between replicates is exactly what one hopes to achieve through normalisation.
Relative gene expression of target gene Ltb4dh in Daphnia magna following a 24 h exposure to ibuprofen  Table 1). Expression levels are displayed relative to the mean control level. An asterisk (*) indicates a significant difference (p < 0.05) between the highest treatment and any other treatment.
Relative gene expression of target gene Ltb4dh in Daphnia magna following a 24 h exposure to ibuprofen  Table 1). An asterisk (*) indicates a significant difference (p < 0.05) from Ltb4dh normalised by NF 3 .
QPCR confirmed that Ltb4dh was up-regulated following exposure to IB [13]. NSAIDs are known to inhibit the biosynthesis of various eicosanoids that play important regulatory and signalling functions, e.g. regulation of ion flux. In mammals, Ltb4dh is involved in the metabolism of leukotriene B 4 an eicosanoid that is formed in the lipoxygenase pathway [15]. There is evidence of eicosanoid pathways in invertebrates being similar to the mammalian pathways [16]. Further experimentation, involving a global expression profile based on cDNA microarrays, is underway to support this hypothesis and reveal the overall molecular stress response of D. magna to IB. To our knowledge, this is one of the first genomic studies to look at the impact of a NSAID on an aquatic invertebrate. The only non-mammalian studies are on yeast Saccharomyces cerevisiae [17], a Drosophila cell line [18] and Bacillus megaterium [19].
Following normalisation to NF 3 it transpired that Atb, assayed as a candidate reference gene, was affected by IB, and so should be regarded as an IB target gene. In a review on the molecular mechanism of NSAIDs in mammals, Tarnawski and Jones, 2003 [20] reveal that it is likely that NSAIDs exert many inhibitory effects, one of which is the inhibition of cell division. The observed suppression of Atb following exposure to IB could be linked to this, since alpha-tubulin helps form microtubules that are involved in cell division, cell structure, and transportation of vesicles and organelles.
Of the remaining candidate reference genes, Act, GAPDH, UBC, WARS, SDH and TBP had a constant expression following exposure to IB. But although, GAPDH and UBC were the most stable reference genes they both had a weakly significant (p < 0.10) linear trend of down-regulation. Thus, relying solely on either of these two reference genes for normalisation could obscure the output.
Using reference genes as an internal control in QPCR is currently the recommended approach [5]. However, Bustin and Nolan, 2004 [21] advocate the introduction of more standard analysis and reporting procedures in QPCR similar to the MIAME guidelines [22] established for microarray technology. In union with other authors [5,10,21,23] we would like to recommend some basic applications for QPCR: (i) apply a dissociation curve in every QPCR, and omit any samples diverging from the curve; (ii) use a minimum of two non-template controls (NTCs) in every QPCR, and omit any samples having a cycle threshold (C t ) value close to or higher than the mean NTC C t -value; (iii) avoid inter-assay comparison if possible, as it is more appropriate to duplicate the samples in one QPCR than using triplicates of each sample, but having to run two QPCRs (S.N. Peirson, personal communication); (iv) use a NF based on a minimum of three validated reference genes for target gene normalisation [10]; (v) integrating the DART-PCR approach [23], which verifies sample-to-sample variations in amplification efficiency with a normalisation factor software package, e.g. geNorm [10], creates a powerful tool for analysing relative gene expression of target genes.

Conclusions
Investigators engaging in toxicogenomic research using invertebrates, or vertebrate cell lines, exposed to NSAIDs may be able to apply our findings directly under similar experimental conditions, e.g. acute exposure. However, we strongly recommend validating reference genes prior to commencing target gene expression profiling under experimental conditions different than those described herein. Our findings stress the importance of performing fundamental research, by validating the most invariant internal QPCR control(s) of a particular experimental condition, prior to investigating the expression of target genes examined by QPCR only, or as QPCR validation of microarray data.

Test species
Daphnia magna Straus (Clone Type 5 -IRCHA), originally obtained from the Water Research Centre (WRc), Medmenham, UK, were investigated. For further information on culturing conditions of D. magna please see Hooper et al., 2006 [24].

Experimental design and analytical chemistry
The experimental design consisted of three replicates of a control and three treatments with ibuprofen-sodium (Sigma-Aldrich, CAS no. 31121-93-4) containing 20, 40 and 80 mg IB l -1 , respectively. Each replicate consisted of 100 third brood neonates (< 24 h old) placed in a 1000 ml glass beaker containing one litre of reconstituted freshwater with or without the addition of IB. Neonates were exposed to IB for 24 h under the same conditions as stock cultures (see Hooper et al., 2006 [24]), but without feeding. To quantify IB, 1.5 ml was sampled from each replicate at 0 h and 24 h and stored at -20°C. Subsequently, the free base concentration of IB was measured at 217 nm by UV-spectrophotometry using an Ultrospec3000 (Biochrom) after the method of Pascoe et al., 2003 [25]. Conductivity, pH and water temperature were checked at the beginning and at the end of the exposure to verify stable water chemistry.

Tissue preparation, RNA extraction and reverse transcription
Following exposure, the neonates were immediately placed in 0.2 ml of RNAlater ® (Ambion) and stored at -80°C. Subsequently, total RNA was extracted using the RNeasy Mini kit with on-column DNase treatment (Qia-gen) to remove any traces of genomic DNA following the manufacturer's instructions. RNA concentrations were determined by spectrophotometry using GeneQuant Pro (Biochrom). The integrity of the RNA was verified by 1% agarose gel electrophoresis. Four micrograms of DNasetreated total RNA was reverse transcribed with Oligo(dT) 12-18 primers (Invitrogen) using the Omniscript Reverse Transcriptase kit (Qiagen) following the manufacturer's instructions. Once synthesised, cDNA was diluted 10-fold resulting in total RNA concentrations of 10 ng/µl and stored at -20°C.

Real-time quantitative PCR
Ten expressed sequence tags (ESTs) were obtained from the NCBI website [26]

Data analysis and statistics
Following QPCR, the raw data was exported to an Excel workbook (Microsoft), entitled Data Analysis for Real-Time PCR (DART-PCR) [23]. This programmed workbook enables calculation of threshold cycles and amplification efficiencies for every sample. Differences in amplification efficiency were assessed using one-way analysis of variance (ANOVA), based upon the null hypotheses: (i) that amplification efficiency is comparable within sample groups (outlier detection) and (ii) that amplification efficiency is comparable between sample groups (amplification equivalence) [23]. Outliers identified by DART-PCR and samples diverging from the dissociation curve were omitted prior to further analysis.
Differences in relative expression were assessed on logtransformed data using ANOVA and Tukey's honestly significant difference (HSD) for post hoc comparisons. Equality of variance was tested using Levene's test. The relative expression of the reference genes calculated by DART-PCR was exported as an Excel workbook (Micro-soft) to geNorm version 3.4 [10]. The "geNorm" algorithm is based on the proposition that, in any set of assays, the expression levels of suitable reference genes should be perfectly correlated. However, because of technical difficulties in assessing correlation in circumstances when there is little variation, a modified procedure is adopted as follows [10]. First, the ratio of expression levels is calculated for any two candidate reference genes. Assuming that the expression ratios should not vary between assays, the variation between the expression ratios is used as an inverse measure, here referred to as gene variability, of the adequacy of the reference genes to be used for normalisation. Thus, the lower the variation, the better the reference genes. Variation is here measured as standard deviation. Given a set of assays that covers the treatments of interest, this procedure can be used to identify the best two reference genes from all possible pairs of reference genes available to the experimenter. A normalisation factor (NF 2 ) is then produced for these two genes, calculated as the geometric mean of their expression levels. The procedure is then extended to assess whether stepwise inclusion of additional reference genes -producing NF 3 to NF n -will reduce the average pairwise variation between them. The optimal number of reference genes to be included in a NF is estimated by comparing the pairwise variation between sequential normalisation factors [10], e.g. NF 3 vs. NF 4 . All statistical tests were performed in SPSS 12.0.1 for Windows (SPSS, Chicago, IL, USA) applying a significant level of 5%, except reference gene variability and NF analysis, which was computed in geNorm version 3.4 [10].