Identification of Suitable Reference Genes for Real Time Quantitative Polymerase Chain Reaction Assays on Pectoralis major Muscle in Chicken (Gallus gallus )

Thirteen reference genes were investigated to determine their stability to be used as a housekeeping in gene expression studies in skeletal muscle of chickens. Five different algorithms were used for ranking of reference genes and results suggested that individual rankings of the genes differed among them. The stability of the expression of reference genes were validated using samples obtained from the Pectoralis major muscle in chicken. Samples were obtained from chickens in different development periods post hatch and under different nutritional diets. For gene expression calculation the ΔΔCt approach was applied to compare relative expression of pairs of genes within each of 52 samples when normalized to mitochondrially encoded cytochrome c oxidase II (MT-CO2) target gene. Our findings showed that hydroxymethylbilane synthase (HMBS) and hypoxanthine phosphoribosyl transferase 1 (HPRT1) are the most stable reference genes while transferrin receptor (TFRC) and beta-2-microglobulin (B2M) ranked as the least stable genes in the Pectoralis major muscle of chickens. Moreover, our results revealed that HMBS and HPRT1 gene expression did not change due to dietary variations and thus it is recommended for accurate normalization of RT-qPCR data in chicken Pectoralis major muscle.


Introduction
Quantitative real-time PCR (RT-qPCR) is a gold-standard biotechnique for gene expression analyses but requires a robust reference genes to correct for technical variation within the experiment [1]. Expression data from RT-qPCR analysis requires reliable normalization to avoid false positive results that may lead to misinterpretations and imprecise conclusions. One of the major difficulties in obtaining accurate expression patterns is the removal of experimentally induced non-biological variation from the true biological variation. A widely used method for normalization involves measuring the expression of a reference gene [2]. A gene can be used as a reference when it is highly and stably expressed in all samples under investigation and is not co-regulated with a target gene [3]. Thus, normalization of gene expression based on a varying reference gene is likely to produce misleading results [4]. According to the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines, reference genes are now selected based on their specificity in interactions between a species or tissue subjected to different experimental treatments. Several studies have reported evaluations of reference gene in animals such as cattle [5], pig [6], sheep [7], goat [8], horse [9], fish [10], and birds [11]. However, a limited number of studies have investigated identification of reference genes used in chicken (Gallus gallus). Moreover, the systematic evaluations of reference genes for chicken specie other than the immune [12] and fibroblast [13] cells are scarce.
The expression stability of reference genes was investigated using five different algorithms. This has enabled us to assess the suitability of 13 common reference genes in Pectoralis major muscle. Samples were obtained from chickens in different development periods post hatch and under different nutritional diets Our results revealed that HMBS and HPRT1 gene are most stable and thus recommended for accurate normalization of RT-qPCR data while TFRC and B2M varied across the experimental conditions and thus are not recommended to normalize gene expression data.

Ethical approval
All experiments were approved by the Universidade Federal de Sergipe Animal Ethics Committee and Brazilian national standards for animal research.

Broiler, experimental condition and muscle sampling
A total of one hundred sixty-eight 1-day-old Cobb male broilers were obtained from a commercial hatchery after being vaccinated against Marek's disease and Newcastle disease. Broilers were allocated in floor pens until day 7. The floor pen had wood shavings and was equipped with tube feeder and nipple drinkers. Temperature was maintained at 32°C at placement, and was gradually reduced to ensure comfort by using a thermostatically controlled heater and fans. To evaluate the stability of expression of the reference genes tested, broilers were submitted to different diets, and slaughtered at different days. At 8 d post hatch, animals were randomly assigned into four groups with seven replicate cages of six birds each. Broilers were fed ad libitum four different growing diets and subsequently four finishing diets for 14 d each. Experimental diets were based on corn and soybean meal and were formulated using the food chemical composition values and nutritional requirements for the average performance of male broilers [14]. Experimental diets were maize-soyabean meal-based supplemented with four levels of L-Lys (12.71, 11.82, 10.99 or 10.22 g/kg) during the growing phase, and four levels of L-Lys (11.25, 10.46, 9.72 or 9.030 g/kg g/kg) during the finishing phase. A total of 52 broilers (4 male in each group) were randomly sampled at 8, 21, 35 and 42 days post hatch and slaughtered. The control group was composed by broilers supplemented with 0 g/kg of L-Lys and slaughtered at 8 d post hatch. Samples of Pectoralis major muscle were isolated and placed in sterile tubes containing RNAlater (Ambion, RNA Carlsbad, CA, USA). Samples were placed at 4°C for 12 h and then stored at 70°C prior to total RNA isolation.

RNA isolation and cDNA synthesis
Total RNA was isolated using the SV Total RNA Isolation System Kit (Promega, Madison, WI, USA) according to the manufacturer's protocol. The RNA concentration was determined by analyzing 2.5 μl of extract using AstraGene UV/Vis Spectrophotometer (AstraNet Inc., Bath, UK). Integrity of RNA was evaluated using 1% agarose gel. As quality control assay, the absence of contaminant genomic DNA in RNA preparations was verified using RNA as a template in real-time PCR assays (minus RT control, i.e. RNA not reverse-transcribed to cDNA). One microgram of total RNA from each group was subjected to oligo (dT) reverse transcription using the GoScript Reverse Transcription System (Promega, Madison, WI, USA) following the manufacturer's instruction. The resulting cDNA samples were stored at -20°C until RT-qPCR analysis was performed.

Selection of reference genes and RT-qPCR optimization
Thirteen reference genes were selected for investigation to identify the most stably expressed reference gene(s) to be used in RT-qPCR studies. Reference genes were chosen based on literature and nucleotide sequences were recovered from the Ensembl database (www.ensembl.org/). These sequences were used to design primers by the PrimerQuest software provided by Integrated DNA Technologies, Inc (IDT, Coralville, IA, USA). Additionally, BLAST analyses were performed to verify their specificity. The primers sequence information used on this study is shown in Table 1 and the primer specificity in the Fig 1. Prior to the RT-qPCR amplification, part of cDNA was pooled into a new tube to make one composed sample which was then diluted to construct the standard curves for PCR conditions optimization and calculation of PCR efficiency. For this purpose, four cDNA quantity (1, 5, 15, 45 ng) and four primers concentrations (200, 400, 800 and 1000 nM) were tested. The following experimental RT-qPCR conditions were used: 1 cycle of 95°C for 10 min, 40 cycles of 95°C for 10 sec and 60 sec at 60-62°C. Additional steps with a gradual increase in temperature from 60-62 to 95°C were used to obtain the dissociation curve. After the analysis of efficiency, the most adequate annealing temperature and primer concentration was used to perform PCR reaction. A reaction mix without template was also used to detected possible reagent contamination. The PCR amplification reaction was performed at different wells and in duplicates. Each primer pair was checked for size specificity of the amplicon by 1.5% agarose gel electrophoresis and ethidium bromide staining. The PCR amplification efficiencies were calculated for each target and reference gene assay using the formula E = (10-1/slope-1) × 100 [15].

Quantitative real-time PCR and preprocessing
Quantitative real-time PCR reactions were performed using SYBR Green detection with GoTaq qPCR Master Mix (Promega, Madison, WI, USA), using gene specific primers. All reactions were performed on CFX96 Real-Time PCR Detection System. A blank (No Template Control) was also added in each assay. The thermal cycling conditions for the RT-qPCR were as follows: 1 cycle at 95°C for 10 min, 40 cycles of amplification at 95°C for 15 s and annealing at 60-62°C for 1 min. The average of the quantification cycle (Cq), defined as the fractional cycle number at which the fluorescence passes the fixed threshold, was determined using manual quantification settings. Values of Cq for control wells were excluded from further analysis, as the values were greater than 35 or not detected.

Determination of expression stability of reference genes
To prepare data set input, Cq values were transformed to relative quantities using the formula (1+E) -ΔCq [16]. The sample with Cq maximum expression was used as the calibrator with a set value of 1 ((1+E) -Cqoriginal-Cqmax ). Thus, the lowest value was zero. To calculate the stability of the candidate genes we have used the BestKeeper [17], NormFinder [18], geNorm Excel [19], geNorm SAS [20] and ΔCt method [21] were used. Relative quantities was used as the input to perform the analysis on NormFinder, geNorm and ΔCt method. No transformed Cq values are required for BestKeeper analysis. In addition, geNorm was used to determine the optimal number and the reference genes required for reliable normalization of RT-qPCR data. All data were analyzed according to the instructions for each software tested. Relative expression of MT-CO2 gene Relative expression of MT-CO2 target gene was normalized by the reference genes through Pfaffl method. The relative expression of target gene was calculated from the geometric average of two best and worst reference genes. The geometric average of the Cq values with a standard deviation of less than 0.5 cycles was used in further calculations. Significance values were obtained by using analysis of variance (ANOVA) for all the tested reference genes between different contrasting groups. Relative expression at 8 days was used as control. A "t" test was applied to detect statistical difference between mean of groups and differences were considered at α = 0.05.

Primer efficiency and specificity
The annealing temperature of 60-62°C was effective for all primers. The RT-qPCR assays had amplification efficiencies between 91% and 104%. The determination coefficients (R 2 ) values were higher than 0.9909 except for UBC (0.9705; Table 1). Primer specificity was evaluated with a single peak in all melting curves and no primer-dimer was detectable indicating optimal performance of the primers (Fig 1). The no-template controls (NTC) that was included in all RT-qPCR assays did not show any positive Cq value throughout the study. Comparing reference genes by descriptive statistics A RT-qPCR assay based on SYBR Green detection was designed for the transcriptional profiling of 13 genes commonly used as reference in experiments of RT-qPCR analysis ( Table 2). Evaluation of quantification cycle (Cq) values by descriptive statistical analysis showed a variability among the different reference genes. The 13 genes were clearly distributed into two different categories of expression level. Six genes (ACTA1, EEF1, GAPDH, LDHA, RPL5, UBC) coded for highly expressed mRNAs with the majority of Cq values between 13 and 20 cycles; seven other genes (ACTB, HMBS, HPRT1, MRPS27, MRPS30, TFRC, B2M) were moderately expressed with Cq values between 21 and 28 cycles ( Table 2). As shown in the Table 2, MRPS27 (CV = 2.27%) and HPRT1 (CV = 2.34%) had the narrowest variance (lowest Cq dispersion), while TFRC (CV = 5.18%) and B2M (CV = 4.57%) had the widest variance (highest dispersion).

Quantitative analysis and determination of expression stability of reference genes
Before stability analysis, the Shapiro-Wilk test was used to determine if the residuals from the ANOVA followed a normal distribution. Data sets that failed the normality test (P < 0.05) were transformed using the logarithm transformation. The resulting transformed variables were consistent with a normal distribution (Table 3). To evaluate the reference genes expression stability, five different algorithms were used: BestKeeper, NormFinder, geNorm Excel, geNorm SAS and ΔCt method. For each method and gene a ranking of stability values was obtained and the most stable gene was indicated. The results are summarized on Table 4. BestKeeper. The best reference genes are those that have the lowest coefficient of variance and standard deviation. Reference genes with SD values >1 are considered not stable and should be avoided [17]. In the present study only UBC had a SD > 1 and thus was removed from the analysis. As a limit of ten genes can be analyzed at once, TFRC and B2M were also not included in the analysis although these genes had SDs below one. In addition, B2M had the highest SD of the remaining genes so was unlikely to be highly ranked as previous describe on the Table 2. The BestKeeper also combines highly correlated genes into an index and the software then compares the correlation between the BestKeeper index and the candidate reference gene to calculate the correlation coefficient (r) value. Genes that have r values close to 1.0 represent the most stable genes. The Pearson correlation coefficient between BestKeeper index and each reference gene gave the highest stability expression to HMBS (r = 0.96) and MRPS30 (r = 0.94); the lowest value was observed to HPRT1 (r = 0.87) followed by ACTA1 (r = 0.88) and LDHA (r = 0.88) ( Table 2).
NormFinder. NormFinder is a Microsoft Excel-based Visual Basic application that employs a model-based approach. Genes showing high variance between the analyzed samples should be avoided. A low value, indicating a low combined intra-and inter-group variation, indicates high expression stability [18]. By NormFinder analysis, the most stable reference genes HMBS (= 0.17) followed by ACTA1 (= 0.21) both ranked among the two most stable reference genes; TFRC (= 1.87) and B2M (= 1.25) were defined as the least stable (Table 4 and Fig 2).
geNorm. geNorm ranks genes based on their average expression stability (M), and the candidate gene possessing the lowest M value is the most stably expressed gene in that set. To identify reference genes with stable expression, geNorm indicates genes with M values below  Table 4. Results of ranking of thirteen reference genes obtained using five different algorithms and overall rank of candidate reference genes for chicken muscle. Three overall best references genes are indicated in bold. the threshold of 1.5. We set the threshold of 1.0 to ensure the selection of the most stable genes. In this work, two version of geNorm algorithms were used: geNorm Excel [19] and geNorm SAS [20].  Table 4). ΔCq method. The comparative ΔCq method assesses the most stable reference genes by comparing relative expression of "pairs of genes" within each tissue sample reach treatment [21]. If the ΔCq value between pairs of genes remains constant for all samples tested, those reference genes are either stably expressed. The results were very similar to those obtained by BestKeeper and NormFinder, which indicated HMBS (0.80) gene as most stable; but slightly different from the results observed on the geNorm SAS and highly different from geNorm Excel. On the other hand, TFRC was the most variable and had the highest ranking in all algorithms tested ( Table 4).

Determination of the optimal number of reference genes
In this work, geNorm Excel was also to predict an optimal number of reference genes. As shown in Fig 3, GAPDH/HMBS pairwise combinations were co-ranked as most stable genes (M = 0.037) and the inclusion of the third reference gene did not contribute significantly to the variation of the normalization factor (V2/3< 0.15). Based on the cut-off value of 0.15, the two most stable reference genes of this dataset would be sufficient for accurate normalization.

Overall rank of reference genes
The stability measurements produced by all algorithms were combined to establish a consensus rank of the stable genes. A comparison of the rankings produced by the five approaches revealed differences as a consequence of the type of algorithm applied. The results are summarized in Table 4. The Spearman correlation coefficients was moderate to high except to pairwise combination between geNorm Excel and geNorm SAS (r = 0.07) indicating a high discrepancy between them. The highest correlation was observed between BestKeeper and the ΔCq method (r = 0.85) followed by NormFinder and ΔCq method (r = 0.84), indicating a high concordance in ranking between them (Fig 4). According to overall ranking, the HMBS (1) was ranked as the most stable gene among 13 reference genes (Table 4). Genes ACTA1 (2) and HPRT1 (3) were nearly always the most stably expressed genes when different algorithm was applied for selection of suitable reference gene. As shown in the Table 4, HMBS, ACTA1 and HPRT1 were top ranked by at least two algorithms. Conversely, TFRC (13), B2M (12) and UBC (11) showed the lower stability in expression level by all five algorithms.

Selection of pairs of reference genes to normalized target gene expression
To determine if there was a co-expression pattern in gene expression between all pairs of the tested genes, Pearson's correlation analysis was performed on gene expression values. The threshold for the correlation among pair of genes was decided as ±0.85. The correlation between all reference genes was moderate to high among the applied algorithms indicating that eight pairwise combinations of genes were co-regulated. The highest correlation was observed between ACTB and EEF1 (r = 0.89) and, the lowest correlation was between the TFRC and MRPS30 (r = 0.00) ( Table 5).
According to results describe above (the overall rank, the Pearson coefficient of correlation and the optimal number of reference genes) HMBS and HPRT1 was the most stable of pair of reference genes and choose to normalize the target gene expression.

Evaluation and validation of selected reference genes using MT-CO2
Two different sets of reference genes were used to normalize gene expression. The best-ranked genes HMBS and HPRT1 formed the first set, and a second set was formed by TFRC and B2M that showed the worst values (Table 4). By using the geometric mean of gene expression of best reference genes (HMBS and HPRT1) were observed a significant difference in expression for MT-CO2 in skeletal muscle of chicken slaughtered at 42 d post hatch compared to control group (Ratio = 0.340, P = 0.025; Fig 5). However, when MT-CO2 was normalized using lower ranked reference genes (TFRC and B2M) a significant decrease of gene expression can be detected (Ratio = 0.201, P = 0.032; Fig 5), indicating that TFRC and B2M are inappropriate to be used as reference genes. In addition, results showed when this worst reference gene set is used MT-CO2 expression was also significantly decreased in skeletal muscle of chicken slaughtered at 35 d post hatch (Fig 5).

Discussion
The use of different reference genes to evaluate relative expression data has an important impact on the final normalized results. Here, to improve the gene expression analysis by RT-qPCR in skeletal muscle of domestic chicken, we evaluated the suitability of 13 common reference genes in Pectoralis major muscle obtained from chickens in different development periods post hatch and under different nutritional diets Then, their expression was analyzed by using five different algorithms: BestKeeper, NormFinder, GeNorm Excel, GeNorm SAS and ΔCt method.
Our study revealed considerable variation in the correlation among the applied algorithms. We expected that geNorm Excel e geNorm SAS would show the most consistency compared with other algorithms in terms of ranking of reference genes. For example, the ACTA1 gene was within the second or third the most stable genes in four of the five algorithms used in this study, except to geNorm SAS. Although they are the same algorithm, geNorm Excel use the real efficiency calculated by standard curve and accept PCR efficiency calculated above 100%. On the other hand, geNorm SAS use the primer efficiency with maximum value up to 100% (Table 6), which explains the divergent results among both algorithms. Ranking inconsistence was also observed between others algorithms. An interesting fact is that UBC (1) is clearly pointed as more stable genes in geNorm Excel results but only ranked the ninth in geNorm SAS and seven in overall rank. Moreover, despite the relatively high correlation between the BestKeeper and the ΔCt algorithms (r = 0.85), the application of these two algorithms delivered identical ranking in only 2 out of the investigated 13 reference genes while between NormFinder and ΔCt algorithms (r = 0.84) seven genes were ranked in the same position. However, Beekeeper ranked four genes (EEF1, GAPDH, MRPS27 and RPL5) in the same position (4) which may explain this apparent inconsistency between this algorithm and others ones. Therefore, based on our observations we recommend the use of at least Δct, NormFinder and geNorm Excel algorithms for selecting the most stably expressed reference genes.
Although this work was not designed to measure the effect of lysine supplementation or periods post hatch in the reference gene expression, it appears that treatments affects the stability rankings of some specific reference genes. Our data showed that expression stability varied considerably among the 13 reference genes tested in chicken Pectoralis major muscle. This was particularly the case with TFRC and B2M, and this would therefore be completely inappropriate to use as reference genes. From all genes evaluated, those that scored highest in their requirements as reference genes were HMBS, ACTA1 and HPRT1. These three genes fulfilled most criteria as suitable reference genes in that it was strongly stable and displayed minimal fluctuation. For an accurate measure of expression levels normalization by multiple reference genes is suggested. However, the strongest correlation was observed between HMBS and ACTA1 expression. Therefore, we selected HMBS and HPRT1 or ACTA1 and HPRT1 as the most appropriated pair of reference genes. HMBS gene has already shown stability in, salmon muscle [22] and bovine muscle [23]. Interestingly, ACTB was ranked as the tenth reference gene and clearly performed with lower stability compared to other reference genes. This result is in accordance with previous report that observed variations in ACTB expression throughout RT-qPCR experiments [8].
Our results revealed an appropriate reference gene set to normalize the MT-CO2 target gene. The Mitochondrially encoded cytochrome c oxidase II is a protein-coding gene component key of the respiratory chain that catalyzes the reduction of oxygen to water and ATP production. Our results clearly demonstrated that using a commonly uncorrected reference gene without verification would yield to misleading gene expression. The results of relative gene expression analysis demonstrate the importance of normalizing to appropriate reference gene set that is unaffected by experimental condition may differ in the expression stability of the reference genes.
A large number of RT-qPCR studies have been carried out concerning the validation of reference genes in birds [11]. However, it has been difficult to find information on appropriate reference genes for use in chicken. To our knowledge, this is the first study that investigated different candidate reference genes for gene expression analyses and seems to be useful in guiding researchers performing gene expression analyses in this specie.

Conclusion
We present experimentally validated comparisons of reference genes for the normalization of RT-qPCR expression data in studies involving lysine supplementation in chicken at different Suitable Reference Genes for Skeletal Muscle of Chickens days post hatch. Our findings highlighted HMBS and HPRT1 as the most stably expressed genes that can be used for normalization of expression data in Pectoralis major muscle of chickens. The identification of these two reference genes seems to be useful in guiding researchers performing gene expression analyses in chicken Pectoralis major muscle and will benefit future expression studies by providing stability in gene expression analyses.
Supporting Information S1 File. Raw data of Ct values according to experimental treatments. (XLS)