Background

Quantitative real-time PCR (qPCR) is presently the most popular method for the mRNA gene expression analysis, because of its high specificity, sensitivity and reproducibility (Wong and Medrano 2005; Bustin et al. 2009). However, the qPCR results require accurate normalization using internal reference genes (RGs) to avoid the quantification errors caused by experimental variations, as several factors, including variations in amount of starting material, enzymatic efficiency and presence of inhibitors, may cause such varying results (Mehta et al. 2010). As a requirement, the RGs should be consistently expressed in all the tissues/cells, in which the target gene is examined and must be unaffected at transcription level by the experimental conditions (Xing et al. 2009). Traditionally, highly expressed genes such as beta actin (ACTB), beta-2 microglobulin (B2M), elongation factor-1 alpha (ELF-1A) and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) were among some of the most commonly used RGs (Thellin et al. 1999; Catalan et al. 2007). However, various studies had shown that expression of these common RGs varies with species, tissue type, physiological state and/or type of treatment (Thellin et al. 1999; Zheng and Sun 2011). Thus, the suitability of particular RGs depends on the system being investigated and the inherent experimental conditions (Suzuki et al. 2000). Apart from these factors, the evaluation of potential RGs stability also depends on the method used for analysis (Cao et al. 2012). The most commonly used normalization algorithms to derive the more appropriate set of RGs are NormFinder (Andersen et al. 2004), geNorm (Vandesompele et al. 2002) and BestKeeper (Pfaffl et al. 2004). These algorithms employ different approaches that allow evaluation of the suitability of gene or set of genes as a normalization factor (Piehler et al. 2010). Hence, the selection of the most stable gene or set of genes under the particular experimental conditions as well as analytical approach is crucial to accurate profiling of gene expression.

The Indian catfish, Clarias batrachus (present name Clarias magur (Hamilton-Buchanan 1822), Ng and Kottelat (2008)), is an air-breathing teleost, endemic to Indian subcontinent (Chonder 1999). Usually the fish inhabits various habitats (fresh and brackish water bodies, muddy marshes, swamps and wet lands), where it frequently encounters low dissolved oxygen (hypoxia) and tolerates temporary water deprivation during summer periods (Saha and Ratha 2007). These features make C. batrachus a useful model for studying gene modulation responses, which may provide insights to better understand its adaptation mechanisms to survive under these adverse conditions. As hypoxia is a strong stress factor, genes from different functional categories are expected to be transcriptionally regulated in the course of hypoxic insult (Rytkonen et al. 2010) and the same is also true for potential RGs. Based on the literature, to shortlist recommendable RGs candidates for hypoxia studies are challenging. To date, there are no reports of established suitable RG or set of genes for expression studies in C. batrachus for any physiological conditions.

Thus, the present study was aimed to assess the appropriateness of the commonly used RGs and to develop a set of genes that can be used for normalization of qPCR gene expression data in the profiling of gene expression in C. batrachus under control and experimental (hypoxic) conditions. The study had tested 11 potential RGs, in 6 tissues, which are involved in the maintenance of basic cellular functions and thus are assumed to be constitutively expressed at constant levels.

Methods

Animals and hypoxic conditions

Live fish of size 16–20 cm (30–80 g) were collected from commercial catches and brought to laboratory for acclimatization for at least for a month, at 5.00 ± 0.1 mg/l dissolved oxygen (DO), with access to air at 22 ± 3 °C. The details of experimental treatments to induce oxidative stress are given in Tripathi et al. (2013). Briefly, fish were divided into 12 batches (six control and six experimental) of three fish each. All the experimental groups (without access to air) were brought to 0.98 ± 0.1 mg/l DO, by fish own respiration after which one group was sacrificed immediately (referred as progressive hypoxia; PH) and the others were incubated for 1, 2, 3, 6 and 12 h before sacrifice. The corresponding control batches were kept under normoxia (5.00 ± 0.1 mg/l DO) for same time periods as that of experimental once in similar manner. Before each sacrifice, fish were weighed and euthanized in 300 mg/l concentration of tricaine methanesulfonate (MS-222, Sigma) to ameliorate suffering. Tissue samples of brain, muscle, liver, spleen, heart and head kidney were collected and snap frozen in liquid N2, until further analysis. The protocols followed were approved by Institutional Animal Ethics Committee (IAEC).

RNA extraction and cDNA synthesis

Homogenization of frozen samples was carried out under liquid nitrogen. Total RNA was isolated from 30 mg tissue samples using the Nucleospin RNA II kit (Macherey–Nagel, Germany) as per manufacturer’s instructions and was treated with DNase I provided within the kit. The quantity of RNA was measured with PICOPET01 spectrophotometer (Picodrop Ltd., Cambridge, UK), while its quality (for possible degradation) was evaluated with denatured gel electrophoresis. One microgram of total RNA was used for cDNA synthesis with the Superscript III first-strand synthesis supermix for qPCR kit (Invitrogen, Carlsbad, CA, USA) according to manufacturer’s recommendations.

Primer design, amplification and efficiency

The selected genes for testing of RGs are given in Table 1. Oligonucleotide primers specific to selected RGs were designed with the PrimerQuestSM tool (Integrated DNA technologies; http://eu.idtdna.com) and were examined for the potential to serve as qPCR references in the study of gene expression in C. batrachus before and after hypoxia exposure (Table 1). The PCR products were examined by gel electrophoresis, where MHCII, eELF4III and 18S rRNA were not amplified. Primer efficiencies (E) for each RG were calculated from tenfold dilution series that were used to make standard curves for each gene. The E value was calculated according to the formula E = 10−1/slope and the specificities of the PCR amplifications were confirmed by melting curve analysis.

Table 1 Primer sequences and amplicon characteristics of selected candidate reference genes analyzed in the present study
Full size table

Quantitative real-time PCR and data analysis

All qPCR reactions were performed in 20 µl total reaction volume (18 µl master mix and 2 µl undiluted cDNA made from 1 μg RNA). The master mix contained 7.2 µl H2O, 0.8 µl of each primer (0.4 µM final concentration) and 10.0 µl of the SYBR Green Mix (Roche Applied Science, Laval, PQ, Canada). The following cycling conditions were used: (1) denaturation, 5 min at 95 °C; (2) amplification repeated 40 times, 10 s at 95 °C, 10 s at 55 °C, and 15 s at 72 °C with ramp rate of 4.4, 2.2, and 4.4 °C/s, respectively; (3) melting curve analysis, 5 s at 95 °C and 1 min at 65 °C with ramp rate of 4.4 and 2.2 °C/s, respectively, then up to 95 °C at a rate of 0.1 °C/s; (4) cooling, 10 s at 40 °C with ramp rate of 2.2 °C/s, in a Light Cycler 480 (Roche Applied Science) system. Quantification cycle values (Cq; Bustin et al. 2009) were obtained by employing the second-derivative maximum method, while absolute quantification was used to derive the concentration values in Light cycler 480 software (version 1.5).

To determine the expression stability of the selected RGs, the transcript level of each gene was measured in six tissues (brain, heart, liver, muscle, spleen and head kidney) of C. batrachus from hypoxia treated and their respective control groups. For all gene, each tissue was tested with 36 biological replicates (fish), i.e., 18 each for before and after hypoxic treatments. There were a total of six hypoxia treatment time intervals and each had six biological replicates. Thus, a total of 37 reactions were studied per reference gene per tissue, including one negative control (template control without cDNA). The tenfold dilution series were set up for ten measuring points for each gene. The Cq values were directly used to compare RNA transcript levels of different genes, while absolute concentration values were used for intra- and inter-group variability analysis using one-way ANOVA. In all cases, α level of 5 % (p < 0.05) was selected to signify statistically significant differences. All statistical tests were performed with SPSS v 12.01 (SPSS, 2003). To select the most stable RG across different samples and experimental conditions, the transcript level of the candidate RGs were used to analyze with the NormFinder (v0.953) (Andersen et al. 2004) and geNorm (v3.5) (Vandesompele et al. 2002) software packages, while the raw Cq values were used with BestKeeper programme (Pfaffl et al. 2004).

Results

Amplification efficiency of potential reference genes

Out of 11 potential RGs tested, 8, i.e., 28S rRNA, TUB, ACTB, B2M, ELF-1A, GAPDH, RPL41 and RPL30, were successfully amplified by qPCR from the brain, heart, muscle, liver, spleen and head kidney tissues of C. batrachus, in control and experimental samples (after hypoxic treatments). The melting curve analysis showed single peak for all above genes (except RPL41), while products from each of the amplifications appeared as a single band of the expected size on agarose gels. Further, for GAPDH, inconsistent expression among the tissues and hypoxic time periods was observed. The Cq range in the expression of GAPDH was 5–40 (in brain and heart), 5–35 (in liver), 5–39.96 (in muscle), 5–39.08 (in spleen) and 5–37.14 (in head kidney). Hence, GAPDH was not included in analysis. So, six genes (28S rRNA, TUB, ACTB, B2M, ELF-1A and RPL30) were included in further study. The E values of these genes were ranged from 1.93 to 1.99 (Table 1).

Expression levels of the reference genes in different tissues

Before hypoxia treatment

Before hypoxic treatment, the Cq values of the six genes ranged between 9.33 ± 0.98 and 34.13 ± 2.02 (Fig. 1a). 28S rRNA was the most highly expressed gene with Cq values ranging between 9.33 ± 0.40 and 11.01 ± 1.35, followed by RPL30 (16.13 ± 1.13–17.57 ± 0.41), ELF-1A (15.76 ± 0.18–18.84 ± 2.07), TUB (16.89 ± 1.02–22.19 ± 0.65) and B2M (19.96 ± 0.61–21.90 ± 0.93), while the lowest expression (highest Cq >30) was exhibited by ACTB (30.36 ± 1.89–34.13 ± 1.56) in all the examined tissues (Table 2). In comparison, expression of TUB was most affected by tissue type and had a significant (p < 0.05) Cq variation of 5.30 between liver and brain. Similarly, ACTB and ELF-1A expression level also exhibited relatively high tissue specificity, with Cq variation of 3.77 (p < 0.05) between heart and head kidney and 3.08 (p < 0.05) between muscle and head kidney, respectively (Table 2). RPL30, 28S rRNA and B2M were the genes that showed the least Cq variation (1.44, 1.68 and 1.94, respectively); however, these inter-group variations were also significantly (p < 0.05) different among tissues.

Fig. 1
figure 1

Expression level of candidate reference genes. Quantification cycle values in different tissue (brain, heart, liver, muscle, spleen and head kidney) samples of Clarias batrachus. Expression levels of selected genes under (a) normoxic and (b) hypoxic conditions are shown as medians (lines), 5th/95th percentile (box) and range (whiskers) n = 72 reaction samples cumulative in each tissue. RPL30 ribosomal protein L30, ACTB β-actin, ELF1A elongation factor 1 alpha, 28S rRNA 28S ribosomal RNA, TUB tubulin, B2M β2-microglobulin

Full size image
Table 2 Crossing point (Cq) values of candidate reference genes in Clarias batrachus tissues before and after hypoxia
Full size table

After hypoxia treatment

Following hypoxic treatments, the Cq values of all the six genes were of similar expression patterns as in control groups (Fig. 1b); however, tissue-specific change in transcript level was observed in all the genes (Table 2). In head kidney, muscle and brain tissues, least changes (Cq variation range from 0.01 to 0.16, p > 0.05) were observed in RPL30; in liver and heart, 28S rRNA exhibited least Cq variation (0.07 and 0.14, p > 0.05), while ACTB was highly affected (Cq variation of 0.92, p > 0.05). In spleen, TUB was the only gene that showed least Cq variation (0.10, p > 0.05), while all other genes showed higher variation (0.43–0.84, p > 0.05). The observed intra-group Cq variations were significantly (p > 0.05) not different from each other.

Expressional stability of the reference genes in C. batrachus tissues

Before hypoxic treatment

The expressional stability analysis of six candidate RGs revealed that the gene-stability measure (M) (range 0.141–0.395) values of all the genes in the six examined tissues were lower than the threshold of 1.5, accepted by geNorm (Fig. 2). geNorm ranked 28S rRNA/TUB as the most stable genes in brain, heart and spleen (Fig. 2a–c). In liver, head kidney and muscle, respectively, ELF-1A/TUB, ELF-1A/28S rRNA, RPL30/TUB gene pairs were found to be most stable (Fig. 2d–f). The pairwise variation (V n/n+1) analysis showed that except for muscle, V 2/3 values in brain, heart, liver and head kidney were less than the proposed geNorm cut-off of 0.15, whereas in spleen, the required M value resulted only after the addition of fifth gene (V 4/5 value) (Table 3). The ranking of stability values from NormFinder analysis showed that in brain, heart and spleen, TUB and 28S rRNA were the best ranked genes, whereas in liver, muscle and head kidney, respectively, RPL30, 28S rRNA/RPL30, TUB and ELF-1A, RPL30 (Table 4) were among the two best ranked genes. The correlation analysis of Cq values performed by BestKeeper showed that 28S rRNA and TUB were two best RG pairs in brain and heart tissues, while 28S rRNA, B2M/RPL30, ELF-1A/B2M, TUB and ELF-1A, B2M in liver, muscle, spleen and head kidney, respectively, were ranked as best stable RGs (Table 5).

Fig. 2
figure 2

Expression stability of candidate reference genes before hypoxia. Expression stability of candidate reference genes in brain (a), heart (b), spleen (c), liver (d), head kidney (e) and muscle (f) tissues of Clarias batrachus, before hypoxic treatment. The expressions of the reference gene(s) were determined by qPCR and data were analyzed with geNorm to calculate the expression stability (M) of each gene. A lower M value indicates more stable expression

Full size image
Table 3 Pairwise variation values of potential candidate reference genes before and after hypoxic treatments as determined by geNorm
Full size table
Table 4 Stability ranking of potential candidate reference genes, before and after hypoxic treatments as determined by NormFinder
Full size table
Table 5 Stability ranking of potential candidate reference genes before and after hypoxic treatments as determined by BestKeeper
Full size table

In relation to hypoxic treatment

The M values of all the genes in the six examined tissues varied between 0.178 and 0.379 (Fig. 3). Based on the least M value the most stable genes in brain, heart, liver, muscle, spleen and head kidney were RPL30/28S rRNA, 28S rRNA/TUB, 28S rRNA/B2M, ELF-1A/TUB, ELF-1A/28S rRNA and ELF-1A/28S rRNA, respectively (Fig. 3a–f). For expression in brain, heart, spleen and head kidney, the V 2/3 values were <0.15, while in liver and muscle V 2/3 values were more than 0.15 (Table 3). However, in liver the addition of third RG lowered the V 3/4 value to less than the proposed cut-off (Table 3), but in muscle even after addition of up to sixth RG, the V value was not found below the proposed cut-off. The results of NormFinder analysis were consistent with that of geNorm, where in brain, heart and liver 28S rRNA was among the highly stable RGs, while in spleen and head kidney, ELF-1A was ranked among stable genes (Table 4). However, in muscle, results of NormFinder were different from geNorm, where ELF-1A was best ranked. The results of BestKeeper analysis showed that in brain, heart and head kidney, 28S rRNA was ranked as highly stable RGs, while in liver, muscle and spleen, respectively, RPL30, TUB and ELF-1A were the most stable genes (Table 5).

Fig. 3
figure 3

Expression stability of candidate reference genes after hypoxia. Expression stability of candidate reference genes in brain (a), spleen (b), head kidney (c), heart (d), liver (e) and muscle (f) tissues of Clarias batrachus, after hypoxic treatment. The expressions of the reference genes were determined by qPCR and data were analyzed with geNorm to calculate the expression stability (M) of each gene. A lower M value indicates more stable expression

Full size image

Discussion

The major aim of this study was to determine a set of stable candidate RGs for catfish, C. batrachus, which can be used for normalizing the qPCR data. As a prerequisite, RGs used not co-regulated or present on the same pathway and would therefore have stable expression in the experiment under consideration (Bower and Johnston 2009). The levels of expression of the RGs have also been suggested for consideration in the process of RG selection (Zhang et al. 2005). In the present study, the candidate RGs selected were involved in a number of different biological processes; TUB and ACTB are both parts of the cytoskeleton (McCurley and Callard 2008), RPL30, RPL41, ELF-1A, 18S and 28S rRNA are all involved in protein synthesis (Martin-Marcos et al. 2007; Nijhof et al. 2009; Thellin et al. 1999), whereas GAPDH, B2M and MHCII are involved in metabolic and immune-related functions (Thellin et al. 1999; Frawley et al. 2011). The expression pattern of the these RGs, under normoxic and hypoxic conditions in different tissues of C. batrachus, showed different levels of expression, as high (28S rRNA, RPL30, ELF-1A), intermediate (TUB and B2M) and low (ACTB). Hence, 28S rRNA, RPL30 and ELF-1A can be suitable candidate for data normalization of high-abundance, TUB and B2M for medium-abundance and ACTB for low-abundance mRNA in expression studies in different tissues. In addition, the least inter- and intra-group variations were observed for RPL30 and 28S rRNA, suggesting that these genes are least affected by tissue types. The ACTB was the least expressed gene and its higher inter- and intra-group variation makes it unsuitable for normalization studies in C. batrachus tissues for experiments involving hypoxic treatments.

Apart from the expression level considerations, stability analysis based on the M values obtained in geNorm, the stability index from NormFinder and descriptive statistics produced by Bestkeeper suggested that all the six genes used in this study would be suitable for normalization of gene expression data in tissue-dependent manner. This was evident from the slight changes in the order of gene stability obtained from each of the three software packages in respect to tissues and hypoxic treatment. However, for C. batrachus, ACTB and B2M, widely used RG internal controls (Lupberger et al. 2002) were the least stable genes identified by all three analysis methods in examined tissues under normoxic and/or hypoxic conditions. These results were consistent with the earlier reports, which found ACTB to be highly regulated by treatments and have unstable expression in hypoxia-cultured human chondrocytes and lipopolysaccharide-stimulated human THP-1 and K562 cells (Cao et al. 2012; Foldager et al. 2009), whereas B2M was the least stable gene in human reticulocytes (Silver et al. 2006).

The present study on C. batrachus revealed novel candidate RGs (RPL30, TUB, 28S rRNA and ELF-1A) with a more stable expression in experimental contexts, in comparison to frequently used RGs (e.g. ACTB, GAPDH and 18S). In catfish, Ictalurus punctatus, 18S gene was found to be the best for normalization under different physiological conditions (LHRH injection, fasting, and low water stress response) (Small et al. 2008). However, there were concerns for the use of 18S gene as a reference, as its much greater level of expression relative to target mRNA (Tang et al. 2007). Moreover, the outcome of the gene stability evaluation differed between the programmes used, as they employ different algorithms (Nijhof et al. 2009). Therefore, in the present study, according to geNorm, NormFinder and BestKeeper, RPL30, TUB, 28S rRNA and ELF-1A were among the three most stable RG candidates in the control as well as in hypoxic treatment group in the six tissues examined. It was also evident that in brain, heart and head kidney tissues, the identified two best ranked RGs would suffice as normalization factors. However, inclusion of second ranked genes, i.e., ELF-1A in liver and muscle, and RPL30 in spleen would be a valid normalization strategy with more accurate and reliable results for the treatments under consideration. The outcome of the present study also emphasizes the consideration of tissue types as an important factor during the selection of the RGs under experimental conditions (Cinar et al. 2012), as the best ranked genes differed across the tissue types. The ranking results produced by three programmes revealed that in brain, heart, liver, muscle, spleen and head kidney, respectively, 28S rRNA/TUB, 28S rRNA/TUB, RPL30/28S rRNA/ELF-1A, RPL30/TUB/ELF-1A, ELF-1A/28S rRNA/RPL30 and ELF-1A/TUB gene pairs were highly stable and may suffice the requirements of RGs (Table 6).

Table 6 The cumulative ranking of candidate reference genes obtained by geNorm, Normfinder and/or BestKeeper programmes
Full size table

Conclusions

The present study clearly demonstrated the importance of employing multiple approaches to select the reliable and stable RGs in conditions, such as normoxia and hypoxia. Using this approach, multiple RG candidates in combination were identified for normalization factors in the tissue-specific manner, under normoxic and hypoxic conditions.