Introduction

There are 27  breeds  of  sheep  in  Iran, 26 of which  are  fat-tailed  breeds. As  Iran is located in arid and semi-arid regions of  the world and experiences periods of feed abundance and scarcity during the year, fat-tail adipose tissue as an energy reserve developed in native sheep  breeds  of  Iran  as an evolutionary  adaptation  that  serves  as a source of energy to increase survival during the shortage in pasture and of feed scarcity. In fat-tailed breeds, deposition  of fat in tail region during feed  abundant season can keep the animal alive during periods of scarcity without considerable elevation of plasma  non-esterified  fatty  acid (NEFA) concentration. To our knowl- edge, there is no study investigating the biological pathways such as gene expres- sion of regulators  and  enzymes  involved  in adipose tissues metabolism of fat-tailed sheep breeds during  periods  of  negative and positive  energy  balance.  Evaluating the effect of negative and positive energy balances on gene expression might reveal  the pathways involved in adipose tissue metabolism, as quantitative real time PCR (RT-qPCR) can produce data with high sensitivity and reproducibility. As  stability of  the  reference  genes   varies   according to type of tissue and physiological stage (Svingen et al., 2015; Kaur et al., 2018), de- fining a suitable set of  reference  genes  is an absolute prerequisite for any RT-qPCR analysis, since a stable reference gene in various physiological and environmental conditions as an internal standard can ac- curately depict changes in expression of target genes. Moreover,  several software programs have been developed to rank the reference candidate genes from the most stable to the least stable one. These soft-

ware programs have been shown  to  rank  the reference genes differently (Najafpanah et al., 2013) due to acquisition  of  differ-  ent algorithms  (Kim et al., 2011).  Hence, the  objective  of  the   current   study   was to evaluate the stability of 8  commonly  used candidate reference genes including Glyceraldehyde 3-phosphate dehydroge- nase  (GAPDH),  Peptidylprolyl  isomerase A (PPIA), Tyrosine 3-monooxygenase/ tryptophan 5-monooxygenase activation protein, Zeta polypeptide (YWHAZ), B-actin, Glucose-6-phosphate dehydroge- nase (G6PDH), RNA polymerase II sub-  unit A (POLR2A), Phosphoglycerate ki- nase 1 (PGK1) and Beta-2-microglobulin (B2M) in adipose and muscle tissues of fat-tailed Lori-Bakhtiari male lambs under periods of negative and positive energy balance by 3 software programs including BestKeeper, NormFinder and geNorm and also consensus ranking of these software programs.

Materials and Methods

Ethics statement

The experiment was done according to  the  recommendations  in  the  Guide  for  the Care and Use  of  Laboratory  Animals of  the  Research  Station   of   Department of Animal Science, University of Tehran, Iran. The protocols were approved by the Animal Care and Use Committee of the University of Tehran Institutional Animal Care and Use Committee.

Animal, housing and sampling

The experiment was carried out at  the Natural Resources & Agricultural Research Farm of  Tehran University, Karaj, Iran. Eighteen Lori-Bakhtiari male lambs with av- erage body-weight of 45.10 ± 3.50 and age

of 5-6 months  were  divided  into  3 groups  of 6 lambs in each treatment according to their body-weight. Lambs were placed  in  individ- ual pens.  The experiment began  after two weeks of  an  adaptation  to  pen  and  lasted  for about  42 days. All lambs  were fed a balanced total mixed ration  (TMR)  formulat- ed by Cornell net carbohydrate and protein system (CNCPS) software  program  1.5  fold  of their maintenance requirement during ad- aptation period. The diet  was consisted of concentrate (44 %) and forage (56 %; alfalfa hay and wheat straw;  Table  S1). The amount of feed  was adjusted weekly according  to lambs body-weight change during the whole experiment. The lambs were  fed  twice  daily  at 8:00 and 17:00  (equal  amount)  and  had  free access to water. At the end of adaptation period, the first group (6 lambs)  was random-  ly selected and weighted after 16 h depriving from feed and slaughtered  to  collect  samples of adipose tissues and longissimus dorsi mus- cle. The remained 2 groups were fed  90,  80 and 70  %  of  their  maintenance  requirement in weeks 1, 2 and 3 of the experiment re- spectively. At the end of week 3, the second group was randomly selected  and  slaughtered to collect samples and the remained group (group  3)  was  fed  ad-libitum  until  the  end  of experiment (day 42) and then was slaugh- tered to collect samples. All samples were immediately frozen in liquid nitrogen, trans- ferred to the laboratory and kept  at  -80  °C until analysis. Lambs were weighed and bled weekly for calculation of changes in body- weight and plasma NEFA concentration.

Total RNA extraction, clean-up, and cDNA synthesis

Total RNA was extracted according  to  the method of Chomczynski and Sacchi (2006) using YTzol reagent (Yekta Tajhiz Azma Co., Tehran, Iran) and treated with

RNase-free DNase I in order to remove the remnant genomic DNA from the samples (TaKaRa, Shuzo, Kyoto, Japan). The RNA abundance was estimated by nanodrop spec- trophotometry at  260  nm,  and  the  purity was checked by determining the absorption ratio at 260/280 nm. The quality of  the extracted RNA was assessed by electropho- resis at 1% agarose-gel that contained ethid- ium bromide. The first-strand  complementa- ry DNA (cDNA) was synthesized from  100 ng of  total RNA by   cDNA  synthesis  kit (M- MuLV Reverse Transcriptase, Cinaclon Co, Tehran, Iran, Cat No; PR911658),   an oligo (dT) primer and random hexamers ac- cording to manufacturer’s instructions. The process of cDNA synthesis was initiated by annealing of the primers at 37 ºC for 1 min followed by cDNA synthesis at 42 ºC for

60  min  and  terminated  by  inactivation  of the reverse  transcriptase  enzyme  at  85  ºC for 5 min. The  synthesized  cDNA was  kept at -20 ºC to be used later.

Primer design

The nucleotide sequence of 8 candidate reference genes belonging to  the  sheep (Ovis aries) was obtained from public databases (GenBank, National Center for Biotechnology Information). Primer pairs were designed according to these  sequenc- es (optimal Tm at 61 °C and GC between 45-50%)  using  primer3Plus  (Untergasser  et al., 2007) online software programs and the  suitability  of   primers   was   evaluat- ed by OligoAnalyzer 3.1 (http://eu.idtdna. com/analyzer/applications/oligoanalyzer/) and OligoCalc (Kibbe, 2007). The speci- ficity of designed primers was examined through PrimerBLAST software of NCBI database (Ye et al., 2012).  The  sequence and some other characteristics of designed primers are presented in Table 1.

Table 1. The sequence and characteristics of primers used for evaluation of expression of reference genes

Real-time RT-PCR

The real-time quantitative PCR was per- formed using SYBR Green I  technology  on an iQ5 System (BioRad, USA). The reac- tions consisted of 10 µL SYBR Green PCR Master Mix  (SYBR biopars,  GUASNR, Iran), 10 pmol (1 µL) of each specific for- ward and reverse primer, 3 µL  of  cDNA,  and 5 µL nuclease free water, for a final volume of 20  µL. Real-time quantitative  PCR was performed for samples with 6 bio- logical replicates. The PCR temperature cy- cling program was an initial  denaturation  at 95 ºC for 15 min  followed  by 40 cycles  at  95  ºC  (denaturation,  15  sec),  62  ºC (anneal-

ing, 30 sec), and 72 ºC (elongation, 30 sec), followed by a final extension at 72 ºC for 5 min. The  amplified  DNA  was  incubated  at 4 ºC, and 5.5 μL of PCR amplified product was purified using  horizontal  electrophore- sis  in  a  2%  agarose   gel  and  visualized   by

ethidium bromide to confirm  the  specificity  of amplified fragments. The efficiency of RT-PCR was assessed for  each  gene  based on the slope of a linear  regression  model. The bulk of  each  cDNA  sample was  used as  a  PCR  template  to   produce   a   graph of the cycle threshold (Ct) in  a range of 10-fold dilution series. The corresponding RT-PCR efficiencies  were  calculated  based on  the  slope  of  the  standard  curve using the following equation: (E=10 -1/slope-1) (Radonić et al., 2004). A melt-curve analysis was conducted for each amplification be- tween 55-95 ºC to  ascertain  that  non-specif- ic products were not amplified. Three soft- ware programs of NormFinder, geNorm and BestKeeper were used to rank the candidate reference genes according to their stability. The arithmetic mean of the reference genes ranks by 3 software programs  was  calculat- ed as consensus ranking.

Statistical analysis

Data were analyzed by GLM  proce-  dure of SAS software (SAS 2002) to evaluate the  difference  in  Ct  value  of  the candidate reference genes and SAS MIXED procedure was used to analyze bodyweight and plasma NEFA concen- tration during periods of negative and positive energy balance. The difference between treatments was considered to be significant if P<0.05.

Results

Induction of Negative energy balance

As it is shown in Figure 1, feeding 90,

80 and 70 % of maintenance requirement respectively in weeks 1, 2 and 3 of negative energy balance period significantly reduced body-weight and  increased plasma NEFA concentration, hence successfully induced negative energy balance in Lori-Bakhtiari lambs. The lambs experienced the  most severe negative energy balance and body- weight loss during the third week of the restricted feeding and by elimination of the restricted  feeding,  started  to  gain  weight and plasma NEFA concentration was  re- turned (decreased) to the basal level.

Figure 1. The bodyweight and plasma NEFA concentration changes during negative and positive energy balances

Reference genes stability

The mean Ct value with standard devi- ation and also Ct distribution of candidate reference  genes  are  presented  in Figures 2 to 7. In mesenteric  adipose  tissue,  lowest and highest Ct values were observed with B-actin and PGK1 respectively (Figure 2). The range of the Ct value distribution was parallel to the  standard deviation of the reference genes as genes with the lowest standard deviation showed Ct values dis- tributed  in  a   narrower   range   (Figure  3). In fat-tail adipose tissue,  the   lowest  Ct  value was observed with B2M and PPIA, whereas POLR2A and YWHAZ showed the highest Ct value ( Figure 4).  The  Ct  value  of G6PDH showed a narrower distribution range compared to other  candidate  refer- ence genes in fat-tail adipose tissue  (Figure 5). In longissimus dorsi muscle, G6PDH showed the lowest standard deviation (Fig- ure 6) and also had the narrowest distribu- tion (Figure 7). There was  no significant change in Ct  value  of  the  reference  genes in mesenteric adipose tissue (Figure 8),  except for B-actin and G6PDH which in- creased  as  the  experiment  progressed (from

19.32 and 21.65 at the beginning of the experiment to 21.62  and  24.19  at  the  end of the experiment respectively  for  B-actin and  G6PDH).  Induction  of  negative  energy

balance increased the Ct value of all can- didate reference genes and shifting to posi- tive energy balance reduced their Ct value, however the difference was significant only for G6PDH (P<0.02) and PGK1 (P<0.05).

The negative energy balance caused a sig- nificant enhancement in Ct value of GAP- DH, B-actin, B2M  and  PGK1  followed  by a reduction in response to positive energy balance, whereas the Ct value of  G6PDH  was reduced as a consequence of negative energy  balance,  however,  the   difference was not significant.

G6PDH  was  the   most  stable  gene  in the mesenteric adipose tissue defined by NormFinder and  geNorm software pro- grams, whereas by BESTKEEPER  soft- ware, it ranked 3  and  B-actin  was  defined as the most stable gene (Table 2). Gene expression of POLR2A, PGK1 and B2M showed the least stability in mesenteric ad- ipose tissue calculated by NormFinder and geNorm software programs and POLR2A  was replaced by YWHAZ when  BestKeep- er was used. Arithmetic mean of  the  rank- ing by 3 software programs showed that GAPDH, B-actin and PPIA were the most stable and POLR2A, PGK1 and B2M were the least stable genes in mesenteric adipose tissue during negative and positive energy balances.

Figure 2. The Ct value with standard deviation of mesenteric adipose tissue

 Figure 3. The distribution of Ct value of mesenteric adipose tissue

   Figure 4. The Ct value with standard deviation of fat-tail adipose tissue

   Figure 5. The distribution of Ct value of fat-tail adipose tissue

Figure 6. The Ct value with standard deviation of longissimus dorsi muscle tissue

   Figure 7. The distribution of Ct value of longissimus dorsi muscle tissue

Figure 8. The Ct value of reference genes in mesenteric adipose tissue in different energy balances

Table 2. The candidate genes ranked by different software programs and the consensus ranking in mesenteric depot

Ranking of 8 candidate  reference  genes in  fat-tail  adipose  tissue   by   NormFind- er and geNorm software  programs  was quite similar except for B-actin and B2M which   were   exchanged   between   ranks of

5 and 8 (Table 3). PPIA, PGK1 and YWHAZ  were  the  most   stable   refer- ence genes defined by NormFinder and geNorm,  whereas  Best  Keeper  calculat-  ed G6PDH, YWHAZ and POLR2A  as genes  with  least  variability  with  YWHAZ

as the only similarity among 3 software programs. PGK1 was defined as the sec-  ond stable reference gene by NormFinder and geNorm  software  programs,  whereas  it  showed  the  least   stability   calculated  by BestKeeper software program.  Average of the ranking by 3 software programs showed that PPIA, YWHAZ  and  POL- R2A were the most and B-actin, GAPDH and B2M were the least stable genes in fat-tail adipose tissue.

Table 3. The candidate genes ranked by different software programs and the consensus ranking in fat-tail depot

When mesenteric and fat-tail adipose tis- sues were considered together,  GAPDH, PPIA  and  YWHAZ  were  considered  as the most stable genes during negative and positive energy  balances  by  NormFinder  and geNorm software programs, while by using Best Keeper, G6PDH, B-actin and

POLR2A were defined as genes with least variability (Table 4). PPIA which was de- fined as the second stable gene by Norm- Finder  and  geNorm  software  programs, was considered as a  gene  with  low  stabil- ity (ranked sixth) by BestKeeper software program. In addition, G6PDH which was

considered as the best reference gene with least variability by BestKeeper, was among

genes with the lowest stability defined by other software programs.

Table 4. . The candidate genes ranked by different software programs and the consensus ranking in adipose tissue

For longissimus dorsi muscle, G6PDH, POLR2A and YWHAZ were defined  as  the most stable genes by BestKeeper software program, while by NormFinder and geNorm  software  programs,  B-ac-  tin,  PGK1  and  YWHAZ  were   defined as the most stable reference genes with YWHAZ  as  the  only  similarity  (Table 5). G6PDH was considered as the most stable reference gene in muscle tissue by BestKeeper program, whereas by Norm-

Finder  and  geNorm,  it  was  considered  as a gene  with  low  stability.  In  addi-  tion, GAPDH was considered as a gene with low stability by all 3 software pro- grams. Consensus ranking of all software programs defined B-actin, YWHAZ and PGK1 as the  3  most  and  B2M,  G6P-  DH and GAPDH as the 3 least stable reference genes in muscle tissue during periods of negative and positive energy balance.

Table 5. The candidate genes ranked by different software programs and the consensus ranking in skeletal muscle

genes by different software programs  and  also the consensus ranking were affected

ranking of the software  programs  showed that in mesenteric adipose tissue, GAPDH

was among the  3 most stable  reference genes in all  periods of different energy balance, whereas B-actin was not  consid- ered as a stable reference gene in negative energy  balance.  In  fat-tail  adipose  tissue, the stability of reference genes was con- siderably affected by energy balance, how- ever, B-actin and PGK1 were among the

3 most stable genes in both negative and subsequent positive energy balances. When mesenteric and fat-tail adipose tissues were considered together, GAPDH  was  among the 3 most stable genes in all  periods  and then B-actin and PPIA were the most fre- quent genes selected as the 3 most stable reference genes.

Table 6. Three most stable reference genes in various tissues during neutral, negative and positive energy balances.

Discussion

To date, there is no  study  investigating the underlying mechanisms controlling ad- ipose tissue metabolism including gene ex- pression of regulators and enzymes in fat- tailed sheep breeds profoundly. RT-PCR/ quantitative PCR is a sensitive and reliable analysis for investigation of biological pathways involved in tissue  metabolism (Fan et al.,2013) which needs some stable reference genes as an internal normal- ization factor to depict changes in target genes expression. Researchers choose the reference genes from other carried out researches even in closed species  which does not  seem  suitable  as  there  is  not  any reference gene to be stable in all en- vironmental and physiological conditions and also nutritional treatments. Negative energy balance is the consequence of in- creased demands, reduced intake or both which force the body to use its energy reserve to continue the vital metabolic pathways. Enhanced release of free fatty

acids as a consequence of stimulated lip- olysis in response to negative energy  bal- ance leads to increased plasma concentra-  tion of NEFA. In the current study, plasma NEFA concentration increased more than 4 fold at the end of week 3 compared to the beginning of the experiment which demon- strates  that  feed  restriction  was  influential to induce negative energy balance.

The rankings of the candidate reference genes by NormFinder and geNorm soft- ware programs were similar in all studied tissues with some negligible differences, whereas ranking by BestKeeper software program was totally different  from  those  of NormFinder and geNorm. For example, G6PDH  was  defined  as  the  most   sta-  ble reference gene in fat-tail and muscle tissues by BestKeeper software program, whereas by NormFinder and geNorm soft- ware programs it was considered  as  the least stable reference gene.  The  difference in ranking of the reference genes by dif- ferent software programs is in agreement with Najafpanah et al. (2013) and Kaur et

al. (2018) who reported different  ranking  of candidate reference genes by Norm- Finder and geNorm software programs compared to BestKeeper in different  tis- sues. The geNorm program ranks the can- didate reference genes according to mean pairwise variation in a special  candidate gene compared to all other candidate ref- erence genes and represents it as M value and subsequently by  stepwise  elimination of gene with highest M value (Vande- sompele et al., 2007). NormFinder uses an algorithm  rooted  in  mathematical   model of gene expression and a solid statistical framework to estimates both variation of inter and intra-group and provide a stabili-  ty value for each candidate reference gene (Mallona et al., 2010), whereas BestKeeper determines  the  optimal  reference  genes   by repeated pairwise correlation analysis (Pfaffl 2001). These differences  in mathe- matical algorithms used by various soft- ware programs can explain the observed difference in ranking  of  the  reference genes by different software programs.

Consensus ranking which is the average  of candidate reference genes rank calcu- lated  by  3  software  programs,  showed that  the  3  most   stable   genes   defined   for adipose tissues (GAPDH, PPIA and YWHAZ) are different  from those  de-  fined for longissimus dorsi muscle (B-ac-  tin, YWHAZ and PGK1), except for YWHAZ which was among the 3 most stable  reference  genes  in  both  adipose  and muscle tissues. In the studies of Na- jafpanah et al. (2013) and Bonnet et al. (2013), a significant difference  in  stability of candidate reference genes among var- ious tissues was reported in caprine and bovine respectively. Moreover,  in mesen- teric adipose tissue, GAPDH, B-actin and

PPIA was the first 3 most stable reference genes defined by consensus ranking of 3 software programs, whereas in fat-tail adi- pose tissue, PPIA, YWHAZ and POLR2A were the most stable defined reference  genes. In the study  reported  by  Zhang et  al. (2016), adipose tissues from different depots in rat had different most stable ref- erence gene defined by software  programs of NormFinder, geNorm and BestKeeper. The results of current  study  demonstrate that the stability of  reference  genes  not only varies among different   tissues,  but also various  depots  of  a  special  tissue such  as  adipose  tissue  can  be   influen- tial on reference gene  stability.  Moreover,  as it was shown in Fig  3,  Ct  value  of  some reference genes including GAPDH, B-actin and PGK1 in longissimus dorsi muscle, PGK1 and G6PDH in fat-tail ad- ipose tissue and B-actin and G6PDH in mesenteric  adipose  tissue  were   affected  by induction of negative energy balance. These variation in Ct of  reference  genes  can explain the difference in selection of different 3 most stable genes in different periods of energy balance. Gene expres-  sion of candidate reference genes was affected by physiological stage of dairy cows (Macabelli et al., 2014; Jatav et al., 2016). The results indicate that  the  stabil- ity of a reference gene can be affected by physiological status of the animal.

The results of the current study demon- strate that the stability of the  reference  genes varied between mesenteric and fat-  tail adipose tissues and the level  of  en-  ergy balance affects  the stability of the reference genes. Therefore, the ideal  way  for normalization of data related to RT- PCR/quantitative PCR is to define reference genes separately for different  tissues and