Expression stability of reference genes in the skeletal muscles of beef cattle

The objectives of this study were to evaluate the relative expression stability of five candidate reference genes in the semimembranosus (SM) and longissimus thoracis (LT) of beef steers, heifers and young bulls. The mRNA levels of Beta-Actin, eukaryotic initiation factor-2B Subunit 2, glyceraldehyde-3phosphate dehydrogenase, peptidylprolyl isomerase-A and succinate dehydrogenase complex-subunit A were quantified using real-time polymerase chain reaction (qPCR). The combined analysis using the geNorm algorithm revealed that glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and eukaryotic initiation factor-2B subunit 2 (EIF2B2) were the most stable gene pairs. However, individual experimental conditions showed that succinate dehydrogenase complex-subunit A was most stably expressed in bulls and heifers SM, and in bulls and steers LT. Glyceraldehyde-3-phosphate dehydrogenase was most stably expressed in bull SM, steer SM, bull LT and heifer LT. The expression stability ranking order differed between experimental conditions, but all genes had low expression variability. Therefore, using the two most stable reference genes, namely GAPDH and EIF2B2, would result in more accurate normalizations for quantitative real-time PCR studies in the SM and LT muscles of beef cattle. The need for prior evaluation of candidate reference genes in different muscles and sex groups of beef cattle is thus emphasized by the present results.


INTRODUCTION
Gene expressions of proteolytic enzymes are routinely analyzed to provide insight into the factors that influence the quality of meat.Variations in the mRNA levels of a specific gene are best studied using quantitative real-time *Corresponding author.E-mail: cmberema@unam.na.Tel: + 264 65 223 5000.Fax: +264 65 223 5294.
Author(s) agree that this article remains permanently open access under the terms of the Creative Commons Attribution License 4.0 International License PCR (qPCR) because the method has proven to be more accurate, reliable and the enhanced sensitivity means that even extremely low levels can be quantified (Bustin et al., 2005).Detected variations are usually due to biological changes within the cell, as well as the compounding experimental errors arising from the qPCR process (Bustin, 2010).These technical variations can be controlled through normalization of qPCR data using a set of housekeeping and reference genes.This is a highly recommended approach because it accounts for differences in RNA quantity and handling and as well as downstream errors that arise during cDNA synthesis (Bustin et al., 2005;Stürzenbaum and Kille, 2001).It is expected that these reference genes are stably expressed under various experimental treatments and physiological conditions (Huggett et al., 2005).This assumption is generally violated due to the fact that the various housekeeping genes used as reference controls can be regulated because of their other functions besides basal cellular metabolism (Thellin et al., 1999).
Indeed, tissue specific differences in expression stability of reference genes have been reported in bovine longissimus muscle (Pérez et al., 2008), liver (Lisowski et al., 2008) as well as in visceral fat tissues, subcutaneous fat and mammary glands (Saremi et al., 2012) and many other experimental conditions (Vorachek et al., 2013).This indicates that there is no consensus regarding the expression stability of a given reference gene.Moreover, information regarding the stability of reference genes in various bovine skeletal muscles and from beef cattle of different sex-age groups is currently limited.Considering that muscle type and sex of the animal will bring about differential expression levels, it is necessary to determine the most stable gene and the number of genes that would be required for accurate normalization in each condition.Thus a systematic analysis of reference genes remains the most suitable strategy towards a robust normalization of qPCR data (Vandesompele et al., 2002).This study therefore aimed to evaluate the expression stability of selected reference genes in the semimembranosus and longissimus thoracis of beef steers, heifers and young bulls.

Study animals and sample collection
This study used 18 Hereford-cross cattle that consisted of 6 young bulls, 6 steers and 6 heifers.The steers were young bulls which were surgically castrated when they were 2 months old.All animals were weaned at 10 months of age and then allowed to feed on high quality grass and grass silage ad libitum.The animal were slaughtered in a licensed commercial abattoir in West Yorkshire, UK, in accordance with the European Community Directive, 86-609-EC.The young bulls, steers and heifers were marketed when they attained 550 kg live weight at the average age of 547, 764, and 889 d, respectively.Meat samples were excised from the semimembranosus (SM) and the longissimus thoracis (LT) of each animal and then submerged in individual tubes containing RNAlater TM reagent to stabilize the RNA.The tubes containing the samples were incubated at 4°C for 24 h and then preserved at -20°C until required for downstream processes that were carried out within one month.

RNA extraction and cDNA synthesis
Total RNA was extracted from 100 mg of RNAlater TM preserved samples using the RNeasy fibrous tissue extraction kit (Qiagen, Germany) according to the manufacturer's protocol.The eluted RNA was thereafter DNAse-treated using the SIGMA AMP-D1 kit (Sigma-Aldrich, St. Louis, USA) to eliminate any carryover genomic DNA.The NanoDrop ND-1000 spectrophotometer was used to determine the concentration and quality of the eluted RNA where 1 unit at 260 nm is equivalent to 40 ng µl -1 of RNA and the 260/280 nm ratio provided an estimate of the quality.The integrity was assessed by performing electrophoresis on a 2% agarose gel and by exposure to UV light using the SYNGENE system (SynGene Ltd., Cambridge, UK).The components of the Verso cDNA synthesis kit (Abgene Ltd., Epsom, UK), which consisted of 1× cDNA buffer, Verso reverse transcriptase, oligo-dT primers and dNTPs, were prepared as a mastermix and then aliquot into thinwalled PCR tubes.500 ng of total RNA was added to each PCR tube and topped with RNase-free water to make up a 20 µl reaction.Each PCR run also had a tube for the no-enzyme control and a non-template control.All tubes were then placed in an Eppendorf Mastercycler Gradient 5331 thermal cycler (Eppendorf AG, Hamburg, Germany) programmed at 42°C for 30 min and 95°C for 2 min for the reverse transcription step and enzyme deactivation step, respectively.Several replicates were prepared for each sample in order to prepare sufficient cDNA for all selected genes.The cDNA was then stored at -20°C until required for downstream processes.

Quantitative real-time PCR
The primers for and Beta-Actin (ACTB), Eukaryotic initiation factor-2B subunit 2 (EIF2B2), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), peptidylprolyl isomerase-A (PPIA) and succinate dehydrogenase complex-subunit A (SDHA), were obtained from a Bovine HKG gene set designed, optimized and supplied commercially by PrimerDesign Ltd (Southhampton, UK).The genes were selected for evaluation as they have distinct functions and their details are outlined in Table 1.An aliquot was taken from each cDNA tube, pooled and then serially diluted for a 5-point standard curve.Triplicate reactions were run for the standard curve analysis and duplicates for all other samples in a 25 μl amplification reaction mixture containing 5 μl of cDNA, 12.5 μl ABsoluteTM Blue QPCR SYBR® Green master mix (ABgene Ltd., Epsom, UK), 300 nM (final concentration) of gene specific primers and PCR-grade water.For the negative control, 5 μl of PCR-grade water was used instead of the cDNA.
qPCR reactions were performed in the Mastercycler® ep realpex thermal cycler (Eppendorf AG, Hamburg, Germany), using the following conditions: 1 cycle of enzyme activation for 15 min at 95°C, 40 cycles of 15 s denaturation at 95°C, 30 s annealing at 60°C and 30 s elongation at 72°C as well as a melting curve analysis step.The resultant PCR products were subjected to electrophoresis on a 2% agarose gel with 1 × TBE buffer.The amplifications were confirmed specific by the presence of a single band of expected size as visualized by exposure to UV light with the SYNGENE bio-imaging system (SynGene Ltd., Cambridge, UK).

Data analysis using geNorm
The qPCR data in the form of Ct values were exported to Microsoft Excel spreadsheet for subsequent processing and analysis.Amplification efficiency (E) of each reference gene was derived from the slope of the standard curve using the formula, E = (10 (- 1/slope) -1)*100.The Ct values were then transformed into relative quantification values (Q) following the delta-Ct formula, Q =E (ΔCt) , where E is the efficiency of a specific gene and ΔCt is the difference between the lowest intra-gene Ct value and the mean Ct of each technical replicate sample (min Ct -sample Ct).These Qvalues were then used to obtain the gene expression stability measure (M) by calculating pair-wise variation with the geNorm applet.The gene with the lowest M-value is the most stably expressed gene in a given environmental or experimental condition.Furthermore, the geNorm applet calculated the V-score, which indicates the benefit of including a 3 rd , 4 th and 5 th gene in the calculation of a normalization factor (NF).The authors (Vandesompele et al., 2002) stated that a V-score less than 0.15 indicate that using additional reference gene to calculate the normalization factor would not result in significant increase in accuracy.The analysis was carried out to rank all the 5 HKGs under six different conditions resulting from 3 groups of animals (bulls, steers and heifers) and 2 muscles (LT and SM).

RNA quality and PCR performance
The RNA extract was of acceptable quality as indicated by the high (> 1.9) A260/A280 absorbance ratios.A single band of expected size was observed on the agarose gel indicating a lack of primer dimer formations and the specificity of the primers.Moreover, the efficiencies of amplification (Table 1) derived from the 5-point standard curve of each gene were greater than 90 % with correlation coefficients (R 2 ) above 0.99.The mean Ct values and standard deviations of each candidate gene tested are as shown in Table 2. Ct values in range of 15 To 24.31 were observed with EIF2B2 being the least expressed while GAPDH appears to be highly transcribed.The most variation was associated with EIF2B2 expression in Bull LT (24.31 ± 0.73) and the least

Variations in gene expression stability
The overall gene expression stability rankings of the 5 reference genes are illustrated in Figure 1.This shows that GAPDH and EIF2B2 were the most stable gene pair.
However, the M-values of the remaining genes were below the 1.5 mark suggested by Vandesompele et al. (2002), an indication that in the present study PPIA, ACTB and SDHA, were also stably expressed, and in that order.On the other hand, the analysis of individual conditions revealed a different pattern of expression stability (Table 3).The expression stability of the reference genes were generally reversed between the SM and LT muscles of steers and heifers.In the SM of steers, GAPDH was the most stable gene together with EIF2B2 but it was more regulated in the LT.In contrast, the GAPDH was more regulated in the SM of heifers as compared to the LT.EIF2B2 also had the highest Mvalues in the SM and LT of bulls, and the LT of heifers, but it was found to be stably expressed as shown by the lowest M-values in the SM of steers and heifers.In contrast, SDHA was ranked highest in 4 experimental conditions except in the steer SM where it appeared to be variably expressed (Table 3).ACTB and PPIA displayed intermediate stability under most conditions but were both found to be stably expressed in the LT of heifers and steers, respectively.Therefore, the order of expression stability differed between the 2 muscles and 3 sex groups (Table 3).

Optimal number of reference genes for normalization
The optimum number of reference genes required for accurate normalization of gene expression data for SM and LT muscles in beef cattle were illustrated in Figure 2. In this study, all experimental conditions revealed a V 2/3 less than 0.15, indicating that the addition of more than 2 HKGs in the calculation of a normalization factor would not offer significant benefits (Figure 2).These results indicated that the geometric average of at least 2 most stable genes would produce an accurate normalization factor and reliable results for the analysis of gene expression in the LT and the SM muscle.The inclusion of a third reference gene reduced the pairwise variation further in the combined analysis and in the LT analysis but not in the SM (Figure 2).Thus the results indicate that the inclusion of a 4 th relatively less stable gene would have no significant beneficial effect on the normalization factors as compared to using 3.

DISCUSSION
The present study was conducted to validate the expression stability of a set of HKGs in the skeletal muscles of cattle from a commercial beef herd that comprised of different gender groups, specifically, young bulls, heifers and steers.The longissimus thoracis and the semimembranosus were chosen as economically important muscles in the beef industry with relatively distinct meat quality characteristics.The analysis revealed that the expression stability of the 5 reference genes investigated here varied across all experimental conditions.This variability may be attributed to differences in gene regulation and metabolic activities in the two muscles in interaction with gender.Although the expression of many reference genes in bovine have been evaluated in various tissues and under different experimental treatments, the information on their expression stability in skeletal muscle is limited (Bahar et al., 2007;Bonnet et al., 2013;Pérez et al., 2008).
According to Goossens et al. (2005), as well as De Ketelaere et al. ( 2006), SDHA was found to be one of the most stable reference genes in bovine embryos and in polymorphonuclear leukocyte cells of lactating cows, respectively.In contrast, SDHA was ranked highest in term of expression stability in Bovine neutrophils.However, the higher variability in the overall analysis in the present study is probably due to the variable expression levels of SDHA in various muscle fibre types, a reflection of differential demand for the action of the encoded SDHA enzyme in the metabolism of glucose derivatives in various samples.
GAPDH and ACTB are generally thought of as the most stable genes and are frequently used as reference genes for normalization in more than 90% of the studies (Suzuki et al., 2000).However, other studies have shown that these too can be affected by various experimental conditions and requires prior validation (Bionaz and Loor, 2007;Bustin, 2000;Chapman and Waldenström, 2015;Glare et al., 2002;Stamova et al., 2009).The present study found GAPDH to be a suitable endogenous in the SM and LT of young bulls, SM of steers and the LT of heifers.This is in agreement with Bahar et al. (2007) who found GAPDH as the most stable gene in the Longissimus muscle of Limousine cross heifers.Moreover, GAPDH, and ACTB were found to be stably expressed in the livers of cows subjected to different physiological and dietary treatments (Janovick-Guretzky et al., 2007).
In contrast, the expression stability of GAPDH in steer LT of the present study was ranked lowest.Similar results were demonstrated by analysis of the Longissimus muscles of pigs (Erkens et al., 2006;Nygard et al., 2007) as well as in cattle (Pérez et al., 2008), where GAPDH was found to be highly regulated.Indeed, a review by Chapman and Waldenström, (2015), showed that there are only a few studies that found GAPDH and ACTB to be optimal for normalization of gene expression data.GAPDH expression is regulated by a variety of factors including calcium (Chao, et al., 1990) and insulin levels (Nasrin et al., 1990). DeGraff et al. (2010) further demonstrated that androgen treatment regulates the expression of insulin-like growth factor binding protein (IGFBP-2) and GAPDH.These findings suggest a possible sex induced regulation of GAPDH for the steer muscles in the present study.
With regards to ACTB, the gene encodes an important component that is required for muscle contractions and shape.The SM, being a leg muscle undergoes relatively more frequent contractions than the LT.The present analysis indicated a stably expressed ACTB in the LT but a regulated gene in the SM.In the same view, ACTB was ranked fourth most stable gene in bovine Longissimus muscle (Pérez et al., 2008), but highly variable in bovine embryos (Goossens et al., 2005).Therefore, ACTB may not be a suitable reference gene in active muscles and tissues were it appears to be more regulated.EIF2B2 was also found to be stably expressed in the skeletal muscles of steers but slightly regulated in muscles of young bulls.Differences in protein synthesis efficiencies between bulls, steers and heifers as well as between different muscles have been documented.In general, bulls have an advantage over steers and heifers as they gain much faster and produce leaner carcasses (Seideman et al., 1982;Venkata et al., 2015), suggesting enhanced protein synthesis in growing bulls.The EIF2B2 protein may, depending on the circumstances, slow or enhance protein synthesis (Abbott and Proud, 2004) and this could be associated with the differences in the stability of EIF2B2 expression observed in the present study.

Conclusions
Taking into consideration the performance of all the 5 reference genes studied in the skeletal muscles of beef cattle, more accurate normalization factors can be calculated using the geometric means of GAPDH, PPIA and EIF2B2.The present study also suggests that ACTB is not a suitable reference gene especially in more active muscles such as the SM.Thus the expression stability was influenced by muscle type and the sex of the animal.These results contribute to a list of HKGs that may be useful in gene expression studies in Semimembranosus and Longissimus of beef cattle.Further studies are therefore warranted in other economically important muscles of livestock.

Figure 1 .
Figure 1.A combined expression stability analysis of reference genes in bovine skeletal muscles according to geNorm.Higher M-values indicate lower expression stability.

Figure 2 .
Figure 2. Evaluation of the optimum number of reference genes for accurate normalization determined by geNorm software.The Vn/Vn+1 (V) ratio indicates the benefits in the accuracy of normalization for each additional gene (n) included in the calculation of the normalization factor.

Table 1 .
GeneBank accession numbers, functions and amplification information of the 5 bovine reference genes.

Table 2 .
Mean Ct values and standard deviations of the candidate reference genes tested.

Table 3 .
Expression stability ranking of reference genes in the longissimus and semimembranosus muscles of beef cattle.