Influence of heat stress, sex and genetic groups on reference genes stability in muscle tissue of chicken

Quantitative RT-PCR is an important technique for assessing gene expression. However, a proper normalization of reference genes prior to expression analyses of target genes is necessary. The best normalizer is that gene which remains stable in all samples from different treatments. The aim of this study was to identify stable reference genes for normalization of target genes in muscle tissue from three genetically divergent chickens groups (Peloco, Cobb 500® and Caneluda) under environmental (heat stress and comfort) and sex influence. Expressions of ten reference genes were tested for stability in breast muscular tissue (Pectoralis major muscle). Samples were obtained from 36 males and females of two backyard breeds (Caneluda and Peloco) and one commercial line (Cobb 500®) under two environments. The heat stress and comfort temperature were 39 and 23°C, respectively. Animals were housed in the Animal Science Department at Universidade Estadual do Sudoeste da Bahia. We analyzed the expression data by four statistical tools (SLqPCR, NormFinder, Bestkeeper and Comparative CT). According to these tools, genes stability varied according to sex, genetic group and environment, however, some genes remained stable in all analyzes. There was no difference between the most stable genes for sex effect, being MRPS27 more stable for both males and females. In general, MRPS27 was the most stable gene. Within the three genetic groups, the most stable genes were RPL5, HMBS and EEF1 to Cobb 500®, Peloco and Caneluda, respectively. Within the environment, the most stable gene under comfort and heat stress conditions was HMBS and MRPS27, respectively. BestKeeper and Comparative Ct were less correlated (28%) and SLqPCR and NormFinder were the most correlated (98%). MRPS27, RPL5 and MRPS30 genes were considered stable according the overall ranking and can be used as normalizer of relative expression of target genes in muscle tissue of chickens under heat stress.


Introduction
Quantitative reverse transcription PCR (RT-qPCR) has been used as a reliable technique for gene expression quantification, in terms that it is more sensitive than standard PCR analysis for transcripts [1][2][3]. By this technique, gene quantification is obtained as the result of a gene expression of interest (target gene, TG) subtracted from the expression of a normalizer gene (reference gene, RG). Such quantitation can be very useful, but the results are highly dependent on the stability of the used RG for normalization. This procedure aims to minimize technical variations inherent to the quality and quantity of RNA introduced during the extraction procedure, reverse transcription and PCR efficiency [4]. Therefore, standardization of TG expression data is crucial to minimize technical variations that may obscure and affect the results of gene expression quantification [5].
For this technique, RG involved in basic cellular coding processes of protein such as cytoskeleton construction (actin), glycolysis (glyceraldehyde 3-phosphate dehydrogenase-GAPDH), protein folding (cyclophilin), synthesis of ribosomal subunits (rRNA), electron transport and ubiquitin degradation (ubiquitin-UBC) are typically used. Also, ribosomal RNA are used, such as 18S and 28S rRNA. However, an ideal RG must have stable level of expression under various experimental conditions at different stages of development, tissue type and stimuli from external environment (e.g. heat and thermal stress, immunological challenge, among others) [6].
The stability of RG can be analyzed by software that use statistical tools, such as those implemented in geNorm [7], NormFinder [8] e Bestkeeper [9] programs. These programs allow analyzing expression data obtained by RT-qPCR technique and aim to assess the stability of these genes indicating the most appropriate. Normalization procedure requires suitable amplification efficiencies of both reference and target genes [2], therefore, an appropriate choice of RG are essential. Furthermore, the use of a single RG has not been recommended, in terms that it may lead to a poor standardization [7]. According to these authors, multiple RG should be used for proper gene expression quantification. Recently, the use of multiple RG has been used and recommended as a suitable approach for TG normalization in chicken [10].
Several studies have indicated suitable RG in cattle [11], pigs [12], ovine [13], goats [14], horse [15] and fish [16,17]. Also, studies have been carried out to identify suitable RG in chickens (Gallus gallus) [10,18]. However, to date, there are no studies reported that have evaluated suitable RG as normalizer under similar conditions as proposed in this study (genetic group, sex and environment).
Considering the above, we aimed to identify stable reference genes for normalization of target genes in muscle tissue from three genetically divergent groups of chicken (Peloco, Caneluda and Cobb 500 1 ) considering the influence of environment (heat stress and comfort) and sex.

Animals
In this study, we used 36 males and females birds, being 12 chicks of each breed (Peloco and Caneluda, backyard breeds, and Cobb 500 1 , commercial line). Commercial animals were purchased a week after the birth of backyard birds in the Universidade Estadual do Sudoeste da Bahia (UESB), Itapetinga Campus, and raised under the same environmental conditions from November 2 to December 2 of 2015 with an average temperature of 26.5˚C. The predominant climate of Itapetinga region is the semi-arid, in which the temperature increases during the day and decrease during the night. Nutritional diets followed the requirements of the Brazilian Tables for Poultry and Swine [19] and the feed was produced in the poultry sector of UESB (Table 1). All birds were raised in semi-open stalls and lined with wood shavings (wood chips).

Heat stress
Heat stress was provided in two stages so that all birds had the same slaughter age (30 days). First, six birds of Peloco breed and six birds of Caneluda breed were subjected to heat stress under an average temperature of 39.5˚C and environmental relative humidity of 60% for one hour. In the second stage, six birds of Cobb 500 1 line were subjected to heat stress with the same conditions of temperature and humidity for 30 minutes. During the heat stress period, animals had ad libitum access to food and water.
During the heat stress, birds were constantly observed for behavioral changes, in order to avoid deaths caused by excessive temperature. The acute heat stress was determined at the moment that most of the birds (±90%) were prostrate (lying with the abdominal faced down), and with increased respiratory rate. All birds were slaughtered by cervical dislocation after heat stress period.
All control birds (six birds of each genetic group) were slaughtered by cervical dislocation at the second stage of the experiment at 4 am (local time) to ensure thermal comfort temperature (23˚C).

Tissue sampling
After slaughter, muscle tissue samples were collected from the breast (Pectoralis major), placed in cryogenic tubes, identified and stored in liquid nitrogen. Samples were transported to the Veterinary Genetics Laboratory of the Universidade Estadual de Santa Cruz (UESC), separated and stored in Ultrafreezer (-80˚C).

Reference genes selection and RT-qPCR optimization
Sequences of ten reference genes were selected from literature [18] and used for analysis of expression stability using RT-qPCR (Table 2). To obtain the standard curve, we used a cDNA pool of all treatments aiming to optimize and calculate the PCR efficiency. We used three cDNA concentrations (15,45 and 135ng/μl) and three primer concentrations (200, 400, 800 mM). RT-qPCR reaction conditions were set with initial denaturation temperature at 95˚C for two minutes, and 40 cycles of denaturation at 95˚C for 15 seconds. The extension temperature was individually standardized for each pair of primer for 60 seconds. At the end of amplification reaction, we included an additional step with gradual temperature increasing from 60 to 95˚C for dissociation curve analysis. Amplification of all genes was performed in duplicate in a 7500 Fast Real Time PCR System (Applied Biosystems, Foster City, CA, USA). Results were obtained by using the Sequence Detection Systems software (V. 2.0.6) (Applied Biosystems Foster City, CA, USA) that generated the cycle threshold (Ct) parameter. The Ct values of duplicates were obtained directly from the above program and used to calculate the average Ct and standard deviation. PCR amplification efficiency was calculated for each reference gene using the following formula: E = (10^( -1/slope)-1 ) x 100 [20]. After efficiency analysis, the most appropriate annealing temperature and primer concentration were used in PCR reactions.

Real time quantitative PCR
The RT-qPCR reaction was performed using SYBR Green detection kit with GoTaq qPCR Master Mix (Promega, Madison, WI, EUA), using specific primers. Gene amplification was performed in duplicate using the Real Time PCR 7500 Fast system (Applied Biosystems, Foster City, CA, EUA) and the results obtained with the Sequence Detection Systems program (V. 2.0.6) (Applied Biosystems, Foster City, CA, EUA) that generated the cycle threshold (Ct) parameters.
The Ct values were exported to Microsoft Excel to calculate the Ct mean, standard deviation and the standard curve for each gene. A negative control (ultra-pure water) also was added in each assay. The qPCR reaction conditions were defined as follow: Initial denaturation at 95˚C during ten minutes and 40 cycles of denaturation at 95˚C for 15 seconds. The extension temperature between 60 and 64˚C during one minute was ideal for all primers. Ct values of control wells were excluded from subsequent analyzes as well as the undetectable values.

Determination of expression stability of reference genes
The average Ct values of reference genes were used to set the input files according to each software. To calculate the stability of reference genes, four different analysis methods were used: SLqPCR [21], NormFinder [8], Comparative Ct (ΔCt) [22], and BestKeeper [9]. All analyzes were performed in the statistical environment R [23], except Bestkeeper, which was analyzed in Microsoft Excel. Analyses were performed according to authors' recommendations of each tool.
The SLqPCR is a statistical analysis package for R software [23] using the method of VAN-DESOMPELE et al. (2002). To perform the stability analysis, it is necessary to transform Ct values in relative values (Q). This method automatically calculates a measure (M) of expression stability for each control gene in a given group of samples. The lower the M value, more stable, and thus recommending a pair of RG as normalizer.
The NormFinder was used as an algorithm for R statistical environment developed by Andersen et al., (2004) and provides information on intra and inter-group variability indicating the best endogenous control [8]. This algorithm indicates the expression of the most stable Table 2. Description of Gallus gallus reference genes, their specific primers used in RT-qPCR analysis and parameters derived from RT-qPCR analysis. All primers were designed by NASCIMENTO et al., (2015). gene by the lower mean values of stability (S), identifying the reference genes according to their groups of inter and intra-expression variation, obtaining a single gene with a stable expression.

ACTA1
BestKeeper analyzes the variability of reference genes expression through the Ct calculation [24]. It uses a comparative Ct that calculates gene expression changes based on relative difference between an experimental sample and normalizer [25]. Data interpretation is based on standard deviation [±CP] that should be less than 1, on standard deviation [± X-Fold] being less than 2 and coefficient of correlation [r], which should be as large as possible [9]. The correlation coefficient (r) and standard deviation (SD) are the main parameters to evaluate the expression of reference genes and must be large and low, respectively [9]. Furthermore, authors suggest excluding genes that have SD above 1.5. They do not indicate the importance of each measurement, leaving to researcher's criteria the decision to choose which parameter will be used. Based on these previous studies, it was decided to use only the SD to select the most stable RG.
The method of Comparative Ct (ΔCt) determines the most stable gene by comparison of relative expression of pairs of genes in each sample of the treatment [22]. The stable reference genes are those that remain constant ΔCt values for all tested samples. For each tool, a ranking of stability was obtained.
RankAggreg package [26] was used in R software by Monte Carlo algorithm to calculate the Spearman distance and obtain the overall ranking. This package sorts from the most stable to the less stable genes, considering the stability values and the frequency that each gene is displayed according to the algorithms of stability analysis tools (SLqPCR, NormFinder, Best-Keeper and Comparative Ct).
Ct data was converted to linear values (2 -ct ) as recommended by LIVAK and SCHMITT-GEN [27]. ANOVA was performed to identify the difference among effects and their interactions (P 0.05).

Efficiency and specificity of primers
To validate the reference genes used in this study, we tested primers efficiency to check their main features prior to RT-qPCR. The annealing temperature of primers ranged from 60 to 64˚C and RG expression efficiency varied between 94 and 109% corresponding to a slope of -3.25 to -3.11. The coefficient of determination values (R 2 ) were greater than 0.99 (Table 3). Primers specificity was assessed by dissociation curve, which showed only one peak indicating that no primer dimer was detected, presenting excellent performance (Fig 1).

Descriptive statistics of reference genes
Ct values < 29 are strong positive reactions indicating abundance of nucleic acid in the sample, Cts of 30-37 are positive reactions indicating moderate amounts of nucleic acid, Cts of 38-40 are weak reactions indicating minimal amounts of nucleic acid [28].
Ten reference genes were analyzed by RT-qPCR. It is clear to note that there is an expression variability among the quantification cycles of the 10 genes. Based on their expression, these genes were grouped into two categories (strong, moderate). Nine genes (ACTA1, LDHA, EEF1, RPL5, UBC, HPRT1, MRPS27, HMBS, MRPS30,) showed strong expression with Ct values ranging between 12 and 29 cycles, and only one gene (TFRC) presented moderate expression, ranging from 25 to 35 cycles. Bestkeeper was used to obtain the range of variation coefficient. The lowest value was observed for MRPS27 gene (CV = 3.2%) and the highest for ACTA1 (CV = 7.8%), as showed on Table 4.

Stability analysis of reference genes expression
Aiming to cover all factors on stability analysis, data were divided into 4 groups of analysis: General group (not considering differences between genetic groups, environment and sex), Genetic group (considering the backyard breeds and commercial line), Environment (acute heat stress and thermal comfort) and Sex (male and female).
SLqPCR. According to this tool, the best reference genes for each treatment were: General,  Table 6.
BestKeeper. TFRC and ACTA1 genes obtained SD greater than 1.5 and were excluded from analysis. The most stable gene recommended by BestKeeper for each group was: General,  Table 7.

Spearman correlation between tools
In order to verify the correlation between tools, Spearman's rank correlation analysis was used by SAS 1 software. Tools were compared by the general order of reference genes. The  (Table 9).

Overall rank of reference genes
The overall ranking of the most stable RG obtained by this package is observed on Table 10. RPL5 gene expression was statistically different only between genetic groups (P = 0.0294). Nevertheless, MRPS30 and MRPS27 were not different among the studied effects. Interactions were not statistically significantly. According to these results, the normalization factor was calculated as the geometric mean of MRPS30, MRPS27 and RPL5 genes. The normalization factor was not different between breed, treatment or sex (Table 11).

Discussion
Analyses of gene expression have optimized the selection of traits of economic interest in animal breeding. RT-qPCR has been the method of choice and recommended for RG analyses in terms that it has been used by many studies in different species, tissues and treatments [11-16, 18, 27]. Recommended values of amplification efficiency were observed for all analyzed primers in this study ranging from 94 to 109%, and R 2 values above 0.99. Annealing temperature between 60 and 64˚C was efficient for amplification of all genes ( Table 2). These efficiency values and temperatures differ from those proposed by NASCIMENTO et al., (2015) which used the same primers in muscle tissues from G. gallus under different treatments. This is due to different conditions in the experimental design and equipment used, highlighting the importance of performing efficiency tests for each experiment.
According to our results, MRPS27 gene was ranked as the most stable on General group, not considering Genetic group, Sex and Environment factors. In contrast to General group,   (7) 1.06 (9) 1.26 (7) 1.46 (9) 0.87 (7)  considering the other conditions (Peloco, Caneluda, Cobb 500 1 , Male, Female, Heat stress, thermal comfort), the most stable genes varied between RPL5 (Peloco), HMBS (Cobb 500 1 and Heat stress) and EEF1 (Caneluda). Besides the General group, MRPS27 gene was ranked as the most stable in terms of thermal comfort, males and females. The evaluated tools presented a good correlation between each other considering the most stable RG, except for Bestkeeper and Comparative Ct, which displayed greater discrepancy (28% of correlation), contrary to results obtained by NASCIMENTO et al., (2015). These tools use the standard deviation as stability parameter, however, Comparative Ct considers values transformed by ΔCt method and Bestkeeper uses the real values of Ct [7,20,22]. The most correlated tools were SLqPCR and NormFinder (95%). A high correlation between these two tools was expected as they are the most used for analyzing the stability of reference genes, and both perform more specific calculations to determine the most stable RG.
Interestingly, for most of the tools (SLqPCR, NormFinder and Bestkeeper) ACTA1 gene was above the sixth position in the overall ranking, however, on Comparative Ct (ΔCt) analysis this same gene was considered as the most stable for all analyzed factors. Comparative Ct method uses VANDESOMPELE et al., (2002) methodology, but without the same accuracy, only by comparing the expression of pairs of genes [22], which may explain the observed variation in RG ranking. Influence of heat stress in muscle tissue of chicken As there are variations among breeds, temperatures and physiology between males and females, it is expected to occur such variations beyond statistical differences of algorithms. MRPS27 and RPL5 genes remained the most stable for all factors, however, due to the variation between experimental conditions, it is recommended to use a third RG (MRPS30) as a control for target genes normalization. MRPS27 and MRPS30 gene encodes the 28S subunit protein and are related to death associated protein 3 (DAP3), while RPL5 encodes a small protein that is a 60S subunit component and is responsible for transporting rRNA 5S to the nucleolus [29].
Results presented in this study point to different RG than those observed by NASCI-MENTO et al. (2015), which indicated HMBS and HPRT1 genes, among the same pairs of primer, as the most stable to be used as normalizers for gene expression studies in chickens. However, in agreement to NASCIMENTO et al. (2015), our results showed that TFRC gene remains the worst among all tested genes for all categories, except for Cobb 500 1 animals and heat stress environment, and thus not recommended as normalizer. These different results are due the difference between environments and genetic groups, reinforcing the importance of stability analysis study before any experiment aiming to analyze the expression of target genes for certain treatment. Influence of heat stress in muscle tissue of chicken The importance in choosing ideal reference genes is reflected in the normalization of the expression of interest genes. Furthermore, the use of a non-stable reference gene may negatively affect the results of Fold-Change on the ΔΔCT calculation [27], leading to a misinterpretation of the results. The three reference genes most stable in this study are suitable to normalize gene expression data from different genetic groups of chicken that exhibited divergent behavior, as induced by heat stress. Therefore, these genes must be calculated as a normalization factor as the geometrical mean of selected RG and may be applicable for data normalization in a wide range of gene expression studies using different chicken breeds and environmental factors.
Many RT-qPCR studies have tried to validate reference genes in different species and treatments in livestock, including cattle [11], pigs [12], ovine [13], goats [14], horse [15], birds [10,18,30,31] and fish [16,17], as well as in plants [32,33]. To date, no other study analyzing reference genes under these experimental conditions (environments of heat stress and comfort, different genetic groups and sex) in chickens are known. These results may be used to guide future studies under similar experimental conditions in G. gallus. Gene expression is the most basic level in which the genotype of an organism leads to phenotype. Furthermore, genes of interest may have a specific pattern of expression under certain experimental condition, giving a wide range of possibilities for functional studies that will allow the identification of suitable reference genes for gene expression normalization.

Conclusion
Our results indicate RG according to genetic group (Peloco, Caneluda and Cobb 500 1 ), sex and environment (heat stress and comfort) effects. For genetic group effect, the recommended RG for Peloco are RPL5, MRPS30 and MRPS27; for Cobb 500 1 are HMBS, RPL5 and MRPS27; for Caneluda are EEF1, HMBS and HPRT1. For environment effect, the recommended genes are HMBS, HPRT1 and LDHA for thermal comfort; MRPS27, RPL5 and MRPS30 for heat stress. For sex effect, there was no difference among the most stable genes between males and females, and the recommended genes for both are MRPS27, MRPS30 and RPL5.
The MRPS27, RPL5 and MRPS30 genes remained stable among all analyzed factors and may be recommended for normalization of gene expression data in RT-qPCR studies of breast muscle tissue of chickens from different genetic groups and sex under heat stress conditions. These three genes must be used as a normalization factor calculated by geometric mean to be used in the normalization of target genes.
Supporting information S1