Printer Friendly



Coastal blooms of microalgae are a regular worldwide phenomenon. Some species occasionally produce toxic compounds that can be harmful to other organisms. The accumulation of algal toxins in bivalve molluscs constitutes a serious problem to consumers and producers, provoking the closure of shellfish culture areas for long periods and causing important economic losses. Numerous studies have tried to solve this problem or at least to decrease their effect. Most investigations have focused on improving methods to identify and quantify the toxins, and their main objective was to establish monitoring systems to prevent the possible poisoning of shellfish consumers. In fact, almost nothing is known about the metabolism of biotoxins in bivalve molluscs, so it is of great importance to understand the molecular basis of detoxification to design strategies that can accelerate detoxification processes. This is also essential to select populations that might accumulate fewer toxins or eliminate them faster.

The study of gene expression patterns may provide insight into the genes involved in detoxification metabolic pathways. Real-time quantitative polymerase chain reaction (RT-qPCR), a widely used tool to quantify the expression of target genes, is already being used to examine different pathological and toxicological situations in molluscs (Dondero et al. 2005, Araya et al. 2008, Siah et al. 2008, Morga et al. 2010, Wan et al. 2011, Miao et al. 2014, Chen et al. 2015) and other marine animals (Olsvik et al. 2005, 2008, Spinsanti et al. 2006, Sellars et al. 2007, Infante et al. 2008, Mitter et al. 2009, Zheng & Sun 2011).

Real-time quantitative polymerase chain reaction is an extremely sensitive, accurate, reproducible, and rapid technique to measure mRNA levels (Bustin 2000, 2002, Huggett et al. 2005, Kubista et al. 2006). The data obtained by RT-qPCR are typically normalized with an internal control or housekeeping gene. The use of an internal control allows the correction of some systematic errors that may occur during the process such as initial size of the sample and pipetting errors, quantity and quality of RNA obtained, changes in enzymatic efficiencies, or presence of inhibitors (Bustin 2002, Vandesompele et al. 2002).

Normalization with a nonstable reference gene could produce incorrect results (Bustin 2000). Therefore, some mathematical methods have been developed to evaluate the expression stability of candidate genes. The most suitable reference genes can be identified using geNorm (Vandesompele et al. 2002), NormFinder (Andersen et al. 2004), and BestKeeper (Pfaffl et al. 2004) software programs. These selected genes can be then used to normalize the real-time PCR data obtained under the specific experimental conditions. The use of various reference genes to accurately measure the expression of genes under study has been recommended by several authors (Tricarico et al. 2002, Vandesompele et al. 2002, Pfaffl et al. 2004, Huggett et al. 2005).

Highly expressed genes such as [beta]-actin (act), glyceraldehyde-3-phosphate dehydrogenase (gapdh), or IHS ribosomal RNA (18S) have been commonly used as reference genes; however, several authors have demonstrated that their expression under different experimental conditions is not stable, leading to a misinterpretation of the results (Thellin et al. 1999, Hansen et al. 2001, Selvey et al. 2001, Deindl et al. 2002, Dheda et al. 2004, Barber et al. 2005, Volland et al. 2017). This fact also became evident in marine organisms, including molluscs, where their expression varies under different conditions (Aoki et al. 2000, He et al. 2005, Olsvik et al. 2005, 2008, Pan et al. 2005, Ingerslev et al. 2006, Spinsanti et al. 2006, Zhang & Hu 2007, Araya et al. 2008, Siah et al. 2008, Morga et al. 2010, Wan et al. 2011, Zheng & Sun 2011, Volland et al. 2017). In summary, there is not one universal reference gene and it should always be prevalidated for each study conditions (Vandesompele et al. 2002, Dheda et al. 2004, Wong & Medrano 2005). Several investigations of gene expression in marine invertebrates have used traditional reference genes without any evaluation (Dondero et al. 2005, 2006, Franzelliti & Fabri 2005, Cellura et al. 2006, Puinean & Rotchell 2006, Banni et al. 2007, Zorita et al. 2007, Brooks et al. 2009, Chen et al. 2017). Real-time quantitative polymerase chain reaction has been widely performed using a single reference gene, frequently not validated or using genes previously validated under different conditions. Less than 40% of these studies attempted to validate the expression stability of the genes used (see Volland et al. 2017 for review).

The aim of the present study is to identify candidate genes for toxic tides studies with the presence of diarrheic shellfish toxin (okadaic acid) in the mussel Mytilus galloprovincialis (Lamarck 1819), which could be used as reference genes in RT-qPCR to normalize the expression of the studied genes. Eight commonly used reference genes were selected: nicotinamide adenine dinucleotide dehydrogenase, gapdh, cytochrome c oxidase subunit 1 (cox1), 40S rihosomal protein S27 (rps27), translation initiation factor 5 A (tif5a), 40S rihosomal protein S4 (rps4), act, and 18S. The best reference genes in RT-qPCR Mytilus toxicology studies will be selected according to the specifications proposed by Thellin et al. (1999) and Dheda et al. (2005). These genes will be also assessed by performing the relative expression profiles for the M. galloprovincialis mrp1 and mrp2 (Mg mrp1 and Mg mrp2), genes involved in the detoxification of okadaic acid in M. galloprovincialis.


Sample Collection

Wild mussels (Mytilus galloprovincialis) were collected from the culture rafts in Galicia (Spain). Eighteen toxin-free mussels were collected 2 mo after the disappearance of a toxic plankton bloom of Dinophysis microalgae that generate okadaic acid and were used as a control. Fifteen samples (five mussels from different depths: 1, 5, and 10 m) were collected during the early development of the bloom, and 2, 10, 16, and 79 days after the first sampling.

Mussels were dissected to separate the digestive gland (DG), gills (GI), and mantle (MT). Tissue samples were treated with RN Alater (Ambion, Applied Biosystems) according to the manufacturer's instructions and stored at -20[degrees]C. Okadaic acid concentration was measured in all the samples. The total concentration of okadaic acid (ng/g) in the DG of these mussels, after methanol extraction and alkaline hydrolysis of the extracts, was determined by HPLC-MS/MS with online SPE chromatography under alkaline conditions and detection by mass spectrometry triple quadrupole.

Total RNA Extraction and First Strand cDNA Synthesis

Total RNA was isolated from [approximately equal to]20 mg of tissues with Nucleo-Spin RNA II kit (Macherey-Nagel) following the manufacturer's instructions, with an additional precipitation with LiCl and finally diluted in "RNA storage solution" (Ambion). Total RNA was treated with TURBO DNA-free (Ambion) to remove any DNA contamination and to avoid amplification of genomic DNA. RNA samples were stored at -80[degrees]C.

The concentration and purity of all RNA samples were examined by measuring the absorbance ratio at 260/280 nm and at 230/260 nm in a NanoDrop 1000 spectrophotometer (Thermo Scientific). Only pure, protein-free, and organic pollutant--free samples were used. RNA integrity was evaluated by denatured 1% (w/v) agarose gel electrophoresis; a sharp and intense 18S band without smears was observed. The 28S ribosomal RNA band is not usually observed in mussels because, as it has been described in other protostome species (Mauriz et al. 2012), there is dissociation into two equally sized subunits under denaturing conditions which reach the same level as the 18S band.

Total RNA (0.6 [micro]g) from each tissue was reverse transcribed using the iScript cDNA Synthesis Kit (Bio-Rad) in a final volume of 20 [micro]L: 0.6 [micro]g total RNA, 4 [micro]L 5X iScript Reaction Mix, 1 [micro]L iScript Reverse Transcriptase, and DEPC water to 20 [micro]L with a temperature cycle at 25[degrees]C for 5 min, 42[degrees]C for 30 min, and 85[degrees]C for 5 min. All cDNA samples were diluted 1:5 with DNase and RNase-free water and stored at -20[degrees]C.

Candidate Reference Gene Selection and Primer Design

Eight commonly used reference genes were selected in this study (Table 1): NADH dehydrogenase subunit 4 (nd4), gapdh, cox1, rps27, tif5a, rps4, act, and 18S. The Mg mrp1 and Mg mrp2 genes were used to assess the quality of these reference genes. Oligonucleotide primers were designed with OligoAnalyzer 3.1 software ( paying special attention to primer length, annealing temperature, base composition, and 3'-end stability. Primers were synthesized by Thermo Scientific (Germany). To ensure optimal reaction efficiency, amplicons length ranged between 84 and 171 bp. The specificity of amplicons was confirmed by the presence of a single peak in the melting curve and by the presence of a single band of the expected size on 2% agarose gels. The PCR products were subsequently cloned and their identity confirmed by sequencing. The primer sequences, amplicon lengths, and efficiencies are listed in Table 2.

Real- Time Quantitative Polymerase Chain Reaction

SsoFast EvaGreen Supermix (Bio-Rad) real-time PCR was run in duplicate in 96-well reaction plates with the iCycler iQ machine (Bio-Rad). The PCR final volume was 20 [micro]L, containing 10 [micro]L SsoFast EvaGreen Supermix (Bio-Rad), 400 nM forward and reverse primers, and 4 [micro]L of 1:5 diluted cDNA. The cycling conditions were as follows: initial template de-naturation at 95[degrees]C for 30 sec followed by 40 cycles of de-naturation at 95[degrees]C for 5 sec, annealing/elongation at 60[degrees]C for 10 sec, and at 75[degrees]C for 10 sec for fluorescence measurement. These cycles were followed by a melting-curve analysis, ranging from 60[degrees]C to 100[degrees]C, with temperature increases in steps of 0.5[degrees]C every 10 sec. In all cases, nontemplate controls were included to rule out DNA contamination. Melting curve, gel analysis, and sequences of each candidate gene were analyzed to verify that a single PCR product was amplified for each set of primers.

Baseline values were calculated for all plates using the iCycler iQ Software. To ensure comparability between data obtained from different experimental plates, the threshold value was manually set at 200 RFU to automatically calculate Cq values.

To generate the data bases for the determination of PCR amplification efficiency (E) of each transcript, it is recommended to use various dilutions in triplets of a pool of all available cDNAs (Pfaffl et al. 2004). Therefore, for each primer pair, a standard curve was obtained based on known quantities of cDNA (five-fold serial dilutions corresponding to cDNA transcribed from 50 to 0.08 ng of total RNA). Each dilution was assayed in triplicate (Table 2). Polymerase chain reaction amplification efficiency can be either defined as percentage (from 0 to 1, or from 0 to 100) or as folds of PCR product increase per cycle (from 1 to 2). Polymerase chain reaction efficiency (defined as percentage) was calculated with the Bio-Rad iQ software V3.1 from the slope of the standard curve:

E = [10.sup.(-1/slope)] - 1 E(%) = ([10.sup.(-1/slope)] - 1) * 100.

Analysis of Gene Expression Stability

The data obtained for each tissue were analyzed using the three most commonly Excel-based software applications: geNorm, NormFinder, and BestKeeper. For geNorm and NormFinder analysis, Cq values were transformed to relative quantities (Hellemans et al. 2007):

[mathematical expression not reproducible]

geNorm ( (Vandesompele et al. 2002) is a mathematical method based on the principle that the expression ratio of two proper reference genes should be constant across samples. For each reference gene, the pairwise variation with all other reference genes is calculated as the SD of the logarithmic transformed expression ratios, followed by the calculation of a reference gene stability value (M value) as the average pairwise variation of a particular reference gene with all other tested candidate reference genes (Vandesompele et al. 2002). Genes with the lowest M values have the most stable expression. The gene with the highest M is then excluded from the analysis and the calculation is repeated in a stepwise fashion until the best two genes are found.

NormFinder ( (Andersen et al. 2004) is another program which orders the candidate genes according to their mRNA expression stability value (genes with the lowest stability value have the most stable expression). It calculates gene expression stability using a model-based approach to rank best genes in accordance with their minimal combined inter-and intragroup expression variation. This software estimates all the variations of the candidate reference genes and the variation that can occur between sample subgroups of the whole data set. The stability value represents a measure of the systematic error that will be introduced in the final results when this particular gene is used as a reference gene.

BestKeeper ( (Pfaffl et al. 2004) calculates descriptive statistics of the Cq values for each gene and also calculates Pearson correlation coefficients for each candidate reference gene pair. The identification of stably expressed reference genes is based on the fact that their expression levels have to be highly correlated. From the considered stably expressed genes, the BestKeeper Index is calculated as the geometric mean of its candidate reference genes' Cq values (Pfaffl et al. 2004). Then, correlation between each candidate reference gene and the index is calculated. BestKeeper determines the most stably expressed genes based on SD of Cq values and Pearson correlation coefficient (r) of each gene to BestKeeper Index. Stable genes display low SD and high r values.

Gene Expression and Statistical Analysis

The Cq values were transformed to quantities (non-normalized expression) by using the equation:

Q = [(1 +E).sup.-Cq].

The normalized gene expression was calculated as the ratio between Q and the normalization factor (NF). The NF used was the geometric mean of the quantities (Q) of the three most stable reference genes. To facilitate comparisons, these normalized expression levels were divided, for each gene and tissues, by their minimum value (this means that the lowest normalized expression value is equal to 1).

Statistical analyses were conducted with the IBM SPSS Statistics 19.0 package. Data were tested for normality (Shapiro--Wilk test) and for homogeneity of variances. The Student's t-test was used to compare the gene expression levels in the three mussel tissues under two different conditions, with okadaic acid and okadaic acid-free conditions.


Amplification Specificity and PCR Efficiency

Efficiency of RT-qPCR (E), slope values and correlation coefficients (r) were determined using serial 1:5 dilutions of template cDNA [18S starts with an initial dilution of 1:10,000 with respect to the other genes) on a iCycler iQ machine (Bio-Rad). RNA, nontemplate controls, and non-RT controls were included. Polymerase chain reaction amplification efficiency can be calculated with the Bio-Rad iQ software V3.1 from the slope of the standard curve using the following equation: E = [10.sup.(-1/slope)] - 1 (Kubista et al. 2006), expressed as percentage E(%) = [[10.sup.(-1/slope)] - 1] x 100 (values between 0 and 100) or as time of PCR product increase per cycle (from 1 to 2). Hundred percent must be considered the maximum efficiency. Efficiency results were ranged between 111.3% for act in GI and 85.9% for gapdh in the DG, and correlation coefficients varied from 0.988 to 0.999 for tif5a in the DG and 18S in GI, respectively (Table 2).

Transcription Levels of Candidate Reference Genes

The mRNA expression of the eight candidate reference genes (nd4, gapdh, cox1, rps27, tif5a, rps4, act, and 18S) was studied in three different tissues; DG, GI, and MT of Mytilus galloprovincialis. The Cq value indicates the number of cycles at which the fluorescence intensity reaches a specific threshold level of detection. Cq values are inversely associated with the abundance of a specific transcript in the sample. The individual transcription levels of the tested genes showed reduced variations across the studied tissue samples and treatments (Fig. 1). Cq values of candidate genes in M. galloprovincialis tissues ranged between 16.1 ([+ or -]0.6) for the most transcribed gene cox1 and 26.5 ([+ or -]1.4) for the least transcribed gene act. Most genes showed a small difference between the maximum and minimum Cq values (five cycles or less) with the exception of the DG nd4 (Fig. 1).

Expression Stability of Candidate Reference Genes

The aim of this research was to identify reference genes with minimal variability in our samples. To determine the most stable reference genes, the use of raw Cq values alone is not enough to calculate gene expression stability. Therefore, the data obtained were further analyzed and the stabilities of the eight candidate genes were calculated using the three most commonly Excel-based tools: geNorm, NormFinder, and Best-Keeper. Efficiencies, which are obtained from the standard curves for each reaction, were used to transform the Cq values into concentrations and then used as data file for these software programs.

The geNorm suggests that for the three given tissues (DG, GI, and MT), an accurate NF of RT-qPCR data can be calculated by using the two most stably expressed genes, rps4 and gapdh (Table 3). The pairwise variation [V.sub.n/n+1] performed (Fig. 2) shows that the addition of supplementary reference genes will benefit the consistency of the determined NF. According to the geNorm stability rank, the use of a third gene has a significant effect and should preferably be included for final calculation (Vandesompele et al. 2002). This extra gene could be cox1 for the DG and GI, and rps27 for the MT. Curiously, the expression of the genes nd4 and 18S, used as reference genes in many studies, appears to be less stable than other genes.

The arrangement obtained using NormFinder program (Table 3) appears to be similar to the one determined using geNorm. Some differences were observed in the positions but the most stable genes were practically the same using both programs. NADH dehydrogenase subunit 4 and 18S were defined as the least consistent controls.

The results obtained after BestKeeper statistical and regression analyses are shown on Table 3. The best ranked genes selected by the three programs BestKeeper, geNorm, and NormFinder were quite similar.

Expression Profiles of Mg mrp1 and Mg mrp2

The expression profiles of the Mg mrp1 and Mg mrp2 genes in the DG of mussels under a toxic plankton bloom (that generates an okadaic acid concentration of 6,383.3 ng [g.sup.-1]) are shown in Figure 3. Expression profiles varied clearly when expression is normalized using a set of three selected reference genes (gapdh, rps4, and cox1) or the widely used and less-stable act and 18S genes. The use of act and 18S as reference genes in the presence of okadaic acid leaves to a misinterpretation of the results.


Real-time polymerase chain reaction is widely used for the quantification of gene expression because of its high sensitivity and robustness. Accurate normalization of gene expression levels is a sine qua non condition for reliable results. Many studies have also shown that the use of a single reference gene is not appropriate for normalizing the quantification of gene expression and that the use of a set of several reference genes is more suitable for this purpose (Vandesompele et al. 2002). Unfortunately, there is no stable universal reference gene expressed in different tissues and under different biological and environmental conditions. More than 60% of the reviewed studies by Volland et al. (2017) did not experimentally validate the used reference genes. This warrants the search for stably expressed genes in each experimental system, and for the development of an accurate normalization strategy. In this study, candidate genes for RT-qPCR assays were evaluated in mussels under a toxic plankton bloom. It is difficult to recognize correct genes to use for normalization. The selection of the eight candidate genes was based on known recognized reference genes from other organisms. Ideally, reference genes should show a stable and regular pattern irrespective of how the organism is being treated. Normalizing genes should not be involved in pathways which could be activated under the experimental situation because they will probably be unstable genes. Some of the most frequently used algorithms for the evaluation of reference genes are highly dependent on the proposition that their expression is not coregulated under experimental conditions. That is why the eight genes that have been selected by us represent different functional classes and are involved in different biological processes (Table 1), so the effect of such coregulation should theoretically be minor.

Because the reference genes must be stable under all experimental conditions in which the expression of target genes is studied, they must be carefully sampled to evaluate the reference genes. Therefore, in this study, toxin-free mussels and okadaic acid-poisoned mussels at four different times within the intoxication-detoxification curve have been sampled. Ninety-three samples of each tissue were analyzed. The Cp values obtained (Fig. 1) show that the candidate genes are quite stable in all experimental conditions with the exception of nd4 in the DG, in both toxin-free mussels and in mussels with toxin. In addition, all Cp values oscillate in a fairly narrow range that would allow the use of all these genes in a single experiment.

The geNorm, NormFinder, and BestKeeper software programs are commonly used to determine the most stable housekeeping genes from a group of candidates. Each of these approaches has different advantages and weaknesses. In geNorm, the stability of a candidate gene is determined by pair-wise comparison of variation of expression ratios. Therefore, the quality of a reference gene is influenced by both the number and the set of candidate genes that are included in the analysis. Instead of the most stably expressed gene, geNorm tends to select the gene with the highest degree of similarity to the expression pattern of candidate genes in the whole data set (Hibbeler et al. 2008). BestKeeper algorithm can only be used under limited conditions. It uses BestKeeper index, which is the geometric mean of the candidate genes, that are expressed with an SD lower than 1. This condition may be too strict for many candidate genes under different experimental situations. Furthermore, the quality of each reference candidate gene is determined only by the SD of its expression in different samples. NormFinder software determines the stability of the candidate genes based on an estimate of the inter- and intragroup variation combining the advantages of the two other approaches. The best solution is probably to use different algorithms to assess the stability of the candidate genes and select those well ranked in all the approaches. In this study, both geNorm and NormFinder identified gapdh, rps4, and cox1 as the most stable combination of reference genes to use as normalizing genes in the DG and GI tissues, whereas for the MT, cox1 has to be changed for rps27.

The gapdh gene is frequently used as a reference gene because of its key role in the glycolytic pathway; however, gapdh is associated with several functions during evolution such as exocytotic membrane fusion, nuclear RNA export, DNA replication and repair, apoptosis, or viral pathogenesis (Bustin 2002, Huggett et al. 2005). Glyceraldehyde-3-phosphate dehydrogenase mRNA levels might be regulated under common situations, and in some cases, gapdh could be inappropriate as a reference gene for some experimental conditions. Our expression data confirmed gapdh as a good internal standard for the study of gene expression in mussels under toxic tides. Similar results were obtained with gapdh in the striped dolphin Stenella coeruleoalba (Gray 1866) (Spinsanti et al. 2006), Japanese flounder Paralichthys olivaceus (Temminck & Schlegel 1846) (Aoki et al. 2000), and more recently was validated as a suitable reference gene in the flat oyster Ostrea edulis (Linnaeus 1758) (Morga et al. 2010).

The rps4 and rps27 genes encode proteins of small ribosomal subunit. The rps4 and rps27 genes are structural components of this ribosomal subunit. Although genes encoding ribosomal proteins were previously recommended for use in only less-sensitive detection methods such as Northern blot (Thellin et al. 1999), many recent validation studies have reported that ribosomal protein genes showed remarkable expression stabilities in different cell lines and tissues of mammals (Brinkhof et al. 2006, Spinsanti et al. 2006, Janovick-Guretzky et al. 2007), fishes (Infante et al. 2008), shellfish (Siah et al. 2008), and plants (Barsalobres-Cavallari et al. 2009). The results obtained in this study also suggest that genes encoding ribosomal proteins may be good candidates to replace traditional reference genes as internal controls in RT-qPCR assays.

Cytochrome c oxidase subunit 1 is a mitochondrial gene encoding a component of the mitochondrial respiratory chain complex IV. This gene is one of the most popular markers for population genetics and phylogeographic studies across the animal kingdom. This gene is frequently used as a reference gene because its expression is considered stable. Our results in the DG and GI of Mytilus confirm this observation.

Results also suggest that 18S, nicotinamide adenine dinucleotide dehydrogenase, and act, commonly used as normalizing genes, should be used with caution as a reference gene in mussel RT-qPCR studies because they are usually ranked as the least stable genes in the three tissues analyzed. The 18S gene was usually used as an ideal internal control gene for RT-qPCR and many times for mussel expression works (Dondero et al. 2005, 2006, Franzelliti & Fabri 2005, Puinean & Rotchell 2006), but the abundance in comparison with mRNA makes it unsuitable for normalization (Vandesompele et al. 2002). Many studies evaluating rRNA-level variations (Hansen et al. 2001, Tricarico et al. 2002, Olsvik et al. 2005, Sellars et al. 2007, Araya et al. 2008, Siah et al. 2008) had confirmed that 18S was quite unstable. Our results using Mg mrp1 and Mg mrp2 genes confirmed that normalization with 18S and act is questionable for gene expression data of the Mytilus DG exposed to okadaic acid (Fig. 3). Mytilus galloprovincialis mrp1 and Mg mrp2 encode two multidrug resistance-associated proteins that are involved in the detoxification of okadaic acid in M. galloprovincialis (Lozano et al. 2015).

The act gene encodes for the most abundant intracellular protein in eukaryotic cells and is an important component implicated in the cytoskeleton structure that supports the cell and determines its shape, playing an essential role in the phagocytosis and encapsulation. Frequently, act genes are used as reference genes in numerous RT-qPCR expression studies in molluscs (Cellura et al. 2006, Dondero et al. 2006, Banni et al. 2007), but many researches had demonstrated that act is an inappropriate internal control for RT-qPCR (Selvey et al. 2001, Deindl et al. 2002, Huggett et al. 2005, Zhang & Hu 2007, Lozano et al. 2014, 2015). Our results also suggest that the expression of act is affected under the presence of okadaic acid and is therefore an unstable gene.

It is difficult to compare our results in molluscs with others studies because there is a lack of studies about the selection of reference genes in those organisms (Wan et al. 2011). To date, there are few reports about reference genes in bivalves (Araya et al. 2008, Siah et al. 2008, Morga et al. 2010, Dheilly et al. 2011, Cubero-Leon et al. 2012, Mauriz et al. 2012, Lozano et al. 2014, 2015) and most of them refer to hemocytes (Araya et al. 2008, Siah et al. 2008, Morga et al. 2010, Hurtado-Oliva et al. 2015, Volland et al. 2017). In the Pacific oyster Crassostrea gigas (Thunberg 1793, Dheilly et al. 2011), gapdh and adp-rihosylation factor 1 were the most stable genes. Glyceraldehyde-3-phosphate dehydrogenase is also well ranked in the flat oyster Ostrea edulis. Several authors (Araya et al. 2008, Siah et al. 2008, Morga et al. 2010) had evaluated different control genes in two situations for soft-shell clams [Mya arenaria (Linnaeus 1758)] and found that elongation factor 1 and ribosomal protein (S-18) were the most stable genes. In the disk abalone Haliotis discus discus (Reeve 1846), ribosomal protein L5, elongation factor 1, and succinate dehydrogenase were considered the most stable genes (Wan et al. 2011). These observations suggest the possibility that molluscs of different species or under different experimental conditions require different internal controls for real-time polymerase chain reaction normalization.

The optimal number of reference genes used for normalization is a trade-off between practical considerations and accuracy. Therefore, three reference genes have been suggested to be sufficient in most cases (Vandesompele et al. 2002) when only a few target genes are going to be studied. geNorm software (Vandesompele et al. 2002) can be used to determine the optimal number of reference genes for normalization by means of pairwise variation ([V.sub.n/n+1]), calculated between the normalization factors [NF.sub.n] and [NF.sub.n+1] (using n or n + 1 reference genes, respectively). Vandesompele et al. (2002) have suggested an arbitrary cutoff value of 0.15 below which the inclusion of more reference genes is not required. Ling and Salvaterra (2011) propose that the minimal [V.sub.n/n+1] represents the most stable NF achievable, thus corresponding to the optimal number of reference genes. According to these authors, five reference genes would be required in the DG and MT, and six in the GI because [V.sub.5/6] achieves the lowest values in the DG and MT, and [V.sub.6/7] shows the lowest value in GI (Fig. 2).

This study is the first, to our knowledge, that evaluates a set of candidate genes as reference genes for mRNA expression analysis in mussels under a toxic plankton bloom. It has been shown here that genes traditionally used for reference purpose may lead to variable expression profiles under different experimental conditions and demonstrate, in agreement with previous researches (Wong & Medrano 2005, Volland et al. 2017), the importance of an accurate reference gene validation for each species or tissue and for each new experimental situation. From the panel of the eight genes tested, gapdh, rps4, and cox1 was the most suitable combination of reference genes for the normalization in the DG and gill tissues, whereas gapdh, rps4, and rps27 was the best for the MT.


We are grateful to the Instituto Tecnoloxico do Mar and the Centro de Investigations Marinas (Xunta de Galicia) for providing the samples and to Mary Vazquez for helpful comments on the English version of the manuscript. This research was funded by the Conselleria de Innovation e Industria, Xunta de Galicia (Spain), through the cooperation agreement EPITOX (2008-2011). R. Martinez-Escauriaza and V. Lozano were financially supported by the EPITOX project.


Andersen, C. L., J. L. Jensen & T. F. Orntoft. 2004. Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 64:5245-5250.

Aoki, T.. H. Naka, T. Katagiri & 1. Hirono. 2000. Cloning and characterization of glyceraldehyde-3-phosphate dehydrogenase cDNA of Japanese flounder Paralichthys olivaceus. Fish. Sci. 66:737-742.

Araya, M. T., A. Siah, D. Mateo, F. Markham, P. McKenna, G. Johnson & F. C. Berthe. 2008. Selection and evaluation of housekeeping genes for haemocytes of soft-shell clams Mya arenaria. J. Invertebr. Pathol. 99:326-331.

Banni, M., F. Dondero, J. Jebali, Ft. Guerbej, H. Boussetta & A. Viarengo. 2007. Assessment of heavy metal contamination using real-time PCR analysis of mussel metallothionein mt10 and mt20 expression: a validation along the Tunisian coast. Biomarkers 12:369-383.

Barber, R. D., D. W. Harmer, R. A. Coleman & B. J. Clark. 2005. GAPDH as a housekeeping gene: analysis of GAPDH mRNA expression in a panel of 72 human tissues. Physiol. Genomics 21:389-395.

Barsalobres-Cavallari, C. F., F. E. Severino, M. P. Maluf & I. G. Maia. 2009. Identification of suitable internal control genes for expression studies in Coffea arabica under different experimental conditions. BMC Mol. Biol. 10:1.

Brinkhof, B., B. Spee, J. Rothuizen & L. C. Penning. 2006. Development and evaluation of canine reference genes for accurate quantification of gene expression. Anal. Biochem. 356:36-43.

Brooks, S., B. Lyons, F. Goodsir, J. Bignell & J. Thain. 2009. Biomarker responses in mussels, an integrated approach to biological effects measurements. J. Toxicol. Environ. Health A 72:196-208.

Bustin, S. A. 2000. Absolute quantification of mRNA using real-time reverse transcription polymerase chain reaction assays. J. Mol. Endocrinol. 25:169-193.

Bustin, S. A. 2002. Quantification of mRNA using real-time reverse transcription PCR (RT-PCR): trends and problems. J. Mol. Endocrinol. 29:23-39.

Cellura, C, M. Toubiana, N. Parrinello & P. Roch. 2006. HSP70 gene expression in Mytilus galloprovincialis hemocytes is triggered by moderate heat shock and Vibrio anguillarum, but not by V. splendidus or Micrococcus lysodeikticus. Dev. Comp. Immunol. 30:984--997.

Chen, H., J. Zha, L. Yuan & Z. Wang. 2015. Effects of fluoxetine on behaviour, antioxidant enzyme systems, and multixenobiotic resistance in the Asian clam Corbicula fluminea. Chemosphere 119:856-862.

Chen, K., E. Li. T. Li, C. Xu, Z. Xu, J. G. Qin & L. Chen. 2017. The expression of the [DELTA]6 fatty acyl desaturase-like gene from Pacific white shrimp (Litopenaeus vannamei) under different salinities and dietary lipid compositions. J. Shellfish Res. 36:501-509.

Cubero-Leon, E., C. M. Ciocan, C. Miner & J. M. Rotchell. 2012. Reference gene selection for qPCR in mussel, Mytilus edulis, during gametogenesis and exogenous estrogen exposure. Environ. Sci. Pollut. Res. Int. 19:2728-2733.

Deindl. E, K. Boengler, N. van Royen & W. Schaper. 2002. Differential expression of GAPDH and beta3-actin in growing collateral arteries. Mol. Cell. Biochem. 236:139-146.

Dheda, K., J. F. Huggett, S. A. Bustin, M. A. Johnson, G. Rook & A. Zumla. 2004. Validation of housekeeping genes for normalizing RNA expression in real-time PCR. Biotcchniques 37:112-119.

Dheda, K., J. F. Huggett, J. S. Chang, L. U. Kim, S. A. Bustin & M. A. Johnson. 2005. The implications of using an inappropriate reference gene for real-time reverse transcription PCR data normalization. Anal. Biochem. 344:141-143.

Dheilly. N. M.,C. Lelong, A. Huvet & P. Favrel. 2011. Development of a Pacific oyster Crassostrea gigas 31, 918-feature microarray: identification of reference genes and tissue-enriched expression patterns. BMC Genomics 12:468.

Dondero, F., A. Dagnino, H. Jonsson, F. Capri, L. Gastaldi & A. Viarengo. 2006. Assessing the occurrence of a stress syndrome in mussels Mytilus edulis using a combined biomarker/gene expression approach. Aquat. Toxicol. 78S:13-24.

Dondero, F., L. Piacentinia, M. Banni, M. Rebelo, B. Burlando & A. Viarengo. 2005. Quantitative PCR analysis of two molluscan metallothionein genes unveils differential expression and regulation. Gene 345:259-270.

Franzelliti, S. & E. Fabri. 2005. Differential HSP70 gene expression in the Mediterranean mussel exposed to various stressors. Biochem. Biophys. Res. Commun. 336:1157-1163.

Hansen, M. C, A. K. Nielsen, S. Molin, K. Hammer & M. Kilstrup. 2001. Changes in rRNA levels during stress invalidate results from mRNA blotting: fluorescence in situ rRNA hybridization permits renormalization for estimation of cellular mRNA levels. J. Bacterial. 183:4747-4751.

He, N., Q. Qin & X. Xu. 2005. Differential profile of genes expressed in hemocytes of white spot syndrome virus-resistant shrimp (Penaeus japonicus) by combining suppression subtractive hybridization and differential hybridization. Antiviral Res. 66:39-45.

Hellemans, J., G. Mortier, A. Paepe, F. Speleman & J. Vandesompele. 2007. qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol. 8:R19.

Hibbeler, S., J. P. Scharsack & S. Becker. 2008. Housekeeping genes for quantitative expression studies in the three-spined stickleback Gasterosteus aculeatus. BMC Mol. Biol. 9:18.

Huggett, J., K. Dheda, S. Bustin & A. Zumla. 2005. Real-time RT-PCR normalisation; strategies and considerations. Genes Immun. 6:279-284.

Hurtado-Oliva, M. A., S. J. Gomez-Hernandez, J. N. Gutierrez-Rivera, N. Estrada, P. Pina-Valdez, M. Nieves-Soto & M. A. Medina-Jasso. 2015. Gender differences and short-term exposure to mechanical, thermic, and mechanical--thermic stress conditions on hemocyte functional characteristics and hsp70 gene expression in oyster Crassostrea corteziensis (Hertlein, 1951). J. Shellfish Res. 34:849-859.

Infante, C, M. P. Matsuoka, E. Asensio, J. P. Canavate, M. Reith & M. Manchado. 2008. Selection of housekeeping genes for gene expression studies in larvae from flatfish using real-time PCR. BMC Mol. Biol. 9:28.

Ingerslev, H. C, E. F. Pettersen, R. A. Jakobsen, C. B. Petersen & H. I. Wergeland. 2006. Expression profiling and validation of reference gene candidates in immune relevant tissues and cells from Atlantic salmon (Salmo salar L.). Mol. Immunol. 43:1194-1201.

Janovick Guretzky, N. A., H. M. Dann. D. B. Carlson, M. R. Murphy, J. J. Loor & J. K. Drackley. 2007. Housekeeping gene expression in bovine liver is affected by physiological state, feed intake, and dietary treatment. J. Dairy Sci. 90:2246-2252.

Kubista, M., J. M. Andrade, M. Bengtsson, A. Forootan, J. Jonak, K. Lind, R. Sindelka, R. Sjoback, B. Sjogreen & L. Strombom. 2006. The real-time polymerase chain reaction. Mol. Aspects Med. 27:95-125.

Ling, D. & P. M. Salvaterra. 2011. Robust RT-qPCR data normalization: validation and selection of internal reference genes during post-experimental data analysis. PLoS One 6:e17762.

Lozano, V., R. Martfnez-Escauriaza, C. Bernardo-Castineira, C. Mesias-Gansbiller, A. J. Pazos, J. L. Sanchez & M. L. Perez-Paralle. 2014. A novel class of Pecten maximus POU gene, Pma-POU-IV: characterization and expression in adult tissues. J. Exp. Mar. Biol. Ecol. 453:154-161.

Lozano, V., R. Martinez-Escauriaza, M. L. Perez-Paralle, A. J. Pazos & J. L. Sanchez. 2015. Two novel multidrug resistance associated proteins (MRP/ABCC) from the Mediterranean mussel (Mytilus galloprovincialis): characterization and expression patterns in detoxifying tissues. Can. J. Zool. 93:567-578.

Mauriz, O.. V. Maneiro, M. L. Perez-Paralle, J. L. Sanchez & A. J. Pazos. 2012. Selection of reference genes for quantitative RT-PCR studies on the gonad of the bivalve mollusc Pecten maximus L. Aquaculture 370-371:158-165.

Miao, J., L. Pan, W. Zhang, D. Liu, Y. Cai & Z. Li. 2014. Identification of differentially expressed genes in the digestive gland of manila clam Ruditapes philippinurum exposed to BDE-47. Comp. Biochcm. Physiol. C Toxicol. Pharmacol. 161:15-20.

Mitter, K., G. Kotoulas, A. Magoulas, V. Mulero, P. Sepulcre, A. Figueras, B. Novoa & E. Sarropoulou. 2009. Evaluation of candidate reference genes for QPCR during ontogenesis and of immune-relevant tissues of European seabass (Dicentrarchus labrax). Comp. Biochem. Physiol. B Biochcm. Mol. Biol. 4:340-347.

Morga, B., I. Arzul, N. Faury & T. Renault. 2010. Identification of genes from flat oyster Ostrea edulis as suitable housekeeping genes for quantitative real time PCR. Fish Shellfish Immunol. 29:937-945.

Olsvik, P. A., K. K. Lie, A. E. O. Jordal, T. O. Nilsen & I. Hordvik. 2005. Evaluation of potential reference genes in real-time RT-PCR studies of Atlantic salmon. BMC Mol. Biol. 6:21.

Olsvik, P. A., L. Softeland & K. Lie. 2008. Selection of reference genes for qRT-PCR examination of wild Atlantic cod Gadus morhua. BMC Res. Notes 1:47.

Pan, D., N. He, Z. Yang, H. Liu & X. Xu. 2005. Differential gene expression profile in hepatopancreas of WSSV-resistant shrimp (Penaeus japonicus) by suppression subtractive hybridization. Dev. Comp. Immunol. 29:103-112.

Pfaffl, M. W., A. Tichopad, C. Prgomet & T. Neuvians. 2004. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper--Excel-based tool using pair-wise correlations. Biotechnol. Lett. 26:509-515.

Puinean, A. M. & J. M. Rotchell. 2006. Vitellogenin gene expression as a biomarker of endocrine disruption in the invertebrate, Mytilus edulis. Mar. Environ. Res. 62:S211-S214.

Sellars, M. J., T. Vuocolo, A. L. Leeton, J. G. Comana, M. B. Degnan & N. P. Preston. 2007. Real-time RT-PCR quantification of Kuruma shrimp transcripts: a comparison of relative and absolute quantification procedures. J. Biotechnol. 129:391-399.

Selvey, S., E. W. Thompson, K. Matthaei, R. A. Lea, M. G. Irving & L. R. Griffiths. 2001. Beta-actin--an unsuitable internal control for RTPCR. Mol. Cell. Probes 15:307-311.

Siah, A., C. Dohoo, P. McKenna, M. Delaporte & F. C. Berthe. 2008. Selecting a set of housekeeping genes for quantitative real-time PCR in normal and tetraploid haemocytes of soft-shell clams, Mya arenaria. Fish Shellfish Immunol. 25:202-207.

Spinsanti, G., C. Panti, E. Lazzeri, L. Marsili, S. Casini, F. Frati & C. M. Fossi. 2006. Selection of reference genes for quantitative RT-PCR studies in striped dolphin (Stcnella coeruleoalba) skin biopsies. BMC Mol. Biol. 7:32.

Thellin, O., W. Zorzi, B. Lakaye, B. De Borman, B. Coumans, G. Hennen, T. Grisar, A. Igout & E. Heinen. 1999. Housekeeping genes as internal standards: use and limits. J. Biotechnol. 75:291-295.

Tricarico, C, P. Pinzani, S. Bianchi, M. Paglierani, V. Distante, M. Pazzagli, S. A. Bustin & C. Orlando. 2002. Quantitative realtime reverse transcription polymerase chain reaction: normalization to rRNA or single housekeeping genes is inappropriate for human tissue biopsies. Anal. Biochem. 309:293-300.

Vandesompele, J., K. De Preter, F. Pattyn, B. Poppe, N. Van Roy, A. De Paepe & F. Speleman. 2002. Accurate normalization of realtime quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 3:0034.1-0034.11.

Volland, M., J. Blasco & M. Hampel. 2017. Validation of reference genes for RT-qPCR in marive bivalve ecotoxicology: systematic review and case study using copper treated primary Ruditapes philippinarum hemocytes. Aquat. Toxicol. 185:86-94.

Wan, Q., I. Whang, C. Y. Cho, J. S. Lee & J. Lee. 2011. Validation of housekeeping genes as internal controls for studying biomarkers of endocrine-disrupting chemicals in disk abalone by real-time PCR. Comp. Biochem. Physiol. C Toxicol Pharmacol 153:259-268.

Wong, M. & J. F. Medrano. 2005. Real-time PCR for mRNA quantitation. Biotechniques 39:75-85.

Zhang, Z. & J. Hu. 2007. Development and validation of endogenous reference genes for expression profiling of medaka (Oryzias latipes) exposed to endocrine disrupting chemicals by quantitative real-time RT-PCR. Toxicol. Sci. 95:3563-3568.

Zheng, W. & L. Sun. 2011. Evaluation of housekeeping genes as references for quantitative real time PCR analysis of gene expression in Japanese flounder (Paralichthys olivaceus). Fish Shellfish Immunol. 30:638-645.

Zorita, E., E. Bilbao, A. Schad, I. Cancio, M. Soto & M. P. Cajaraville. 2007. Tissue- and cell-specific expression of metallothionein genes in cadmium- and copper-exposed mussels analysed by in situ hybridization and RT-PCR. Toxicol. Appl. Pharmacol. 220:186-196.


Laboratorio de Biologia Molecular y del Desarrollo, Departamento de Bioquimica y Biologia Molecular, Instituto de Acuicultura, Universidade de Santiago de Compostela, Constantino Candeira s/n. Santiago de Compostela 15782, Spain

(*) Corresponding author. E-mail:

DOI: 10.2983/035.037.0108
Name, symbol, and function of reference genes.

Symbol    Gene name                   Function

nd4       NADH dehydrogenase          Respiratory chain (NADH
          subunit 4                   dehydrogenase)
gapdh     Glyceraldehyde-3-           Glycolytic enzyme
cox1      Cytochrome c oxidase        Respiratory chain
          subunit 1                   (cytochrome c oxidase)
rps27     40S ribosomal protein S27   Structural component of 40S
                                      ribosomal subunit
tif5a     Translation initiation      Translation initiation factor
          factor 5A                   activity
rps4      40S ribosomal protein S4    Structural component of 40S
                                      ribosomal subunit
act       [beta]-actin                Cytoskeleton structural
18S       18S ribosomal RNA           Structural constituent of
Mg mrp2   Multidrug resistance--      Detoxification processes
          associated protein

Symbol, primer sequences, amplicon size, and accession number of
candidate reference genes.

Symbol  Primer forward 5'-3'    Primer reverse 5'-3'     Accession


                         Efficiency (%)       Correlation
Symbol   Amplicon (bp)   DG    GI     MT      DG      GI     MT

nd4      114            105.8  100.5  100.9   0.996   1      0.997
gapdh    114             99.3  107.3   96.1   0.992   0.99   0.993
cox1     151             85.9   97.6   91.6   0.996   0.99   0.996
rps27    114             92.0  101.2   97.3   0.996   0.99   0.99
tif5a    171             96.9   99.3  104.6   0.988   0.99   0.991
rps4     121             91.1   93.5   95.8   0.994   1      0.998
act      138             91.0  111.3  100.3   0.993   0.99   0.996
18S       84             99.2   87.7   94.0   0.990   1      0.996
mrp2     105            100.1   --     --     0.990   --     --

Results from geNorm, NormFinder, and BestKeeper analyses for reference
genes in RT-qPCR in the DG, GI, and MT.

Rank  geNorm       Mean M value  NormFinder  Stability value

 1    gapdh jrps4  0.56          gapdh       0.088
 2    gapdh/rps4   0.56          rp.s4       0.105
 3    cox1         0.9           cox1        0.150
 4    act          1.07          tif5a       0.180
 5    tif5a        1.17          rps27/act   0.188
 6    rps27        1.26          rps27/act   0.188
 7    18S          1.41          18S         0.200
 8    nd4          1.69          nd4         0.318
 1    gapdh/rps4   0.63          cox1        0.088
 2    gapdh/rps4   0.63          rps4        0.092
 3    cox1         0.73          gapdh       0.093
 4    act          0.82          act         0.106
 5    rps27        0.88          tif5a       0.136
 6    tif5a        1             rps27       0.137
 7    nd4          1.08          nd4         0.139
 8    18S          1.16          18S         0.188
 1    gapdh/rps4   0.8           rps4        0.094
 2    gapdh/rps4   0.8           rps27       0.114
 3    rps27        1.06          gapdh       0.117
 4    act          1.2           18S         0.118
 5    tif5a        1.34          cox1        0.126
 6    18S          1.43          act         0.149
 7    cox1         1.55          tif5a       0.172
 8    nd4          1.89          nd4         0.339

Rank  geNorm       BestKeeper  r     BestKeeper  SD

 1    gapdh jrps4  gapdh       0.78  gapdh       0.76
 2    gapdh/rps4   rps4        0.77  cox1        0.84
 3    cox1         tif5a       0.71  rps4        0.88
 4    act          act         0.67  18S         1.06
 5    tif5a        rps27/nd4   0.58  rps27       1.11
 6    rps27        rps27/nd4   0.58  act         1.12
 7    18S          cox1        0.52  tif5a       1.19
 8    nd4          18S         0.38  nd4         2.15
 1    gapdh/rps4   cox1/rps4   0.87  cox1        0.56
 2    gapdh/rps4   cox1/rps4   0.87  rps4        0.67
 3    cox1         rps27       0.86  gapdh/18S   0.78
 4    act          gapdh       0.84  gapdh/18S   0.78
 5    rps27        act         0.81  rps27       0.87
 6    tif5a        tif5a       0.77  act         0.90
 7    nd4          nd4         0.65  nd4         0.92
 8    18S          18S         0.33  tif5a       1.08
 1    gapdh/rps4   rps4        0.76  18S         0.89
 2    gapdh/rps4   tif5a       0.74  rps4        0.97
 3    rps27        rps27       0.70  gapdh       1.15
 4    act          gapdh       0.69  rps27       1.18
 5    tif5a        act/18S     0.68  tif5a       1.25
 6    18S          act/18S     0.68  cox1        1.32
 7    cox1         cox1        0.65  act         1.42
 8    nd4          nd4         0.49  nd4         2.37
COPYRIGHT 2018 National Shellfisheries Association, Inc.
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2018 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Martinez-Escauriaza, Roi; Lozano, Vanessa; Perez-Paralle, Maria Luz; Pazos, Antonio J.; Sanchez, Jos
Publication:Journal of Shellfish Research
Article Type:Report
Date:Apr 1, 2018

Terms of use | Privacy policy | Copyright © 2021 Farlex, Inc. | Feedback | For webmasters |