Identification of candidate reference genes for quantitative expression analysis by real-time PCR for hypoxic stress in Indian catfish, Clarias batrachus (Linnaeus, 1758)

The selection of most appropriate reference gene(s) is a crucial step in studies quantifying gene expression by quantitative real-time PCR (qPCR). As a need reference gene(s) should be unaffected at transcription level by experimental conditions and/or tissue types. In this study, 11 candidate reference genes were tested in Clarias batrachus for their expression stability using the geNorm, NormFinder and BestKeeper statistical algorithms, and were compared in different tissues (brain, heart, liver, muscle, spleen and head kidney) and treatments (before and after hypoxia exposure). The results indicated that in brain, heart, liver, muscle, spleen and head kidney, respectively, 28S rRNA/TUB, 28S rRNA/TUB, RPL30/28S rRNA, RPL30/TUB, ELF-1A/28S rRNA and ELF-1A/TUB gene pairs were highly stable and were suitable as reference genes to study oxidative stress, while ACTB and B2M were the least stable genes in examined tissues under normoxic and/or hypoxic conditions. The observations suggested the consideration of tissue types and use of at least two reference genes, instead of one in accurate normalization of qPCR data.


Progressive hypoxia Cq Quantification cycle
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.

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 N 2, 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 PrimerQuest SM 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.
Quantitative real-time PCR and data analysis All qPCR reactions were performed in 20 ll total reaction volume (18 ll master mix and 2 ll undiluted cDNA made from 1 lg RNA). The master mix contained 7.2 ll H 2 O, 0.8 ll of each primer (0.4 lM final concentration) and 10.0 ll 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, a 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).

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).

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 ( (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.

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 ( 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 BHYP before hypoxia treatment, AHYP after hypoxia treatment 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 cutoff 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).

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 Fig. 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 \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 Norm-Finder 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).  The numbers in brackets indicate stability values BHYP before hypoxia treatment, AHYP after hypoxia treatment

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 tissuedependent 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 Fig. 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 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).

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. Ranking order is deduced from interpretation of ranking of reference genes in normoxic and hypoxic conditions and best gene pair was ranked accordingly when at least two algorithms were showing the similar results