Selection of Reference Genes for Gene Expression Studies Related to Intramuscular Fat Deposition in Capra hircus Skeletal Muscle

Abstract

The identification of suitable reference genes is critical for obtaining reliable results from gene expression studies using quantitative real-time PCR (qPCR) because the expression of reference genes may vary considerably under different experimental conditions. In most cases, however, commonly used reference genes are employed in data normalization without proper validation, which may lead to incorrect data interpretation. Here, we aim to select a set of optimal reference genes for the accurate normalization of gene expression associated with intramuscular fat (IMF) deposition during development. In the present study, eight reference genes (PPIB, HMBS, RPLP0, B2M, YWHAZ, 18S, GAPDH and ACTB) were evaluated by three different algorithms (geNorm, NormFinder and BestKeeper) in two types of muscle tissues (longissimus dorsi muscle and biceps femoris muscle) across different developmental stages. All three algorithms gave similar results. PPIB and HMBS were identified as the most stable reference genes, while the commonly used reference genes 18S and GAPDH were the most variably expressed, with expression varying dramatically across different developmental stages. Furthermore, to reveal the crucial role of appropriate reference genes in obtaining a reliable result, analysis of PPARG expression was performed by normalization to the most and the least stable reference genes. The relative expression levels of PPARG normalized to the most stable reference genes greatly differed from those normalized to the least stable one. Therefore, evaluation of reference genes must be performed for a given experimental condition before the reference genes are used. PPIB and HMBS are the optimal reference genes for analysis of gene expression associated with IMF deposition in skeletal muscle during development.

About the Authors:

Wuzheng Zhu

Affiliation: College of Life Science and Technology, Southwest University for Nationalities, Chengdu 610041, Peoples’ Republic of China

Yaqiu Lin

Affiliation: College of Life Science and Technology, Southwest University for Nationalities, Chengdu 610041, Peoples’ Republic of China

Honghai Liao

Affiliation: College of Life Science and Technology, Southwest University for Nationalities, Chengdu 610041, Peoples’ Republic of China

Yong Wang

* E-mail: [email protected]

Affiliation: College of Life Science and Technology, Southwest University for Nationalities, Chengdu 610041, Peoples’ Republic of China

Introduction

Meat from goats is becoming more widely accepted around the world due to the increasing demand for sustainable foods and the low cholesterol content and high nutritive value of goat meat [1]. IMF content, the meat quality trait with the most economic importance [2], has a positive impact on meat characteristics such as tenderness, flavor and juiciness [3]. Thus, understanding the mechanism underlying IMF accumulation in skeletal muscle is of great importance to meat science. Gene expression analysis is a useful technique because it provides information on the regulation of IMF accumulation.

qPCR, with its high sensitivity, specificity and accuracy as well as ease of use [4], has been widely employed as the method of choice to characterize expression profiles of genes of interest and to verify results from microarray studies. Nevertheless, qPCR requires data normalization to correct for variability in the amount and quality of starting material, RNA stability, content, enzymatic efficiencies, and technique in the qPCR experimental process [5]. To address these problems, several strategies have been implemented for data normalization including quantifying RNA input or number of cells used. However, these methods are questionable because they either do not take into the consideration the imbalance in the abundance of mRNA and rRNA or ignore the efficiency of reverse transcriptase [6,7]. In addition, the amount of RNA maybe insufficient at times and accurate computation of cell number is often infeasible. Another proposed strategy is to use reference genes in data normalization; this strategy is currently considered a robust approach in most cases [8]. Reference genes are supposed to be constitutively expressed, and their expression should remain stable irrespective of experimental conditions such as developmental stage, experimental treatment, physiological state and tissue type [9,10]. However, there is an increasing number of reports demonstrating that the expression of most reference genes–including some commonly used reference genes such as GAPDH, ACTB and 18S–is context dependent [5,11], as their expression may vary significantly depending on the experimental conditions being investigated [12]. No reference genes appear to be generally applicable for all circumstances. Accordingly, the validity of candidate reference genes has to be well established to circumvent problems in data normalization.

To date, several studies have attempted to identify suitable reference genes in ruminants. In a previous study on bovine intramuscular fat, RPLP0 was validated as the most stable reference gene in bovine muscle across bovine breeds and ages [13]. In another study, Pérez et.al. [14] suggested that SF3A1, EEF1A2 and HMBS were the most stable reference genes in bovine muscular tissue. Bonnet et.al. [15] demonstrated that UXT, EIF3K, TBP, TOP2B and CLN3 were the best options for data normalization in bovine muscle, liver, mammary gland and adipose tissue. More recently, Najafpanah MJ et.al. [16] carried out an analysis of nine candidate reference genes in goat; HSP-90 was almost always the most stable reference gene in longissimus muscle, liver and visceral and subcutaneous fat.

However, appropriate reference genes for use in the study of goat muscle, especially skeletal muscle during development, are lacking. Furthermore, only a few reports [17,18] on intramuscular fat deposition have appropriately evaluated reference genes in ruminants, and data on qPCR normalization are scarce in Capra hircus. Given the underlying issues associated with the use of non-validated reference genes for data normalization, in this study, we aimed to identify a set of reference genes that could serve as proper reference genes in Capra hircus skeletal muscle. Evaluation of eight reference genes was performed in skeletal muscle across four developmental stages. The geNorm [19], NormFinder [20] and BestKeeper [21] algorithms were applied to assess the stability of reference gene expression; qPCR procedures were performed following Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines [22].

Materials and Methods

Animal and sample collection

All animal procedures were performed according to protocols approved by the Southwest University for Nationalities Animal Care and Use Committee in Sichuan, China.

Twenty-four male Jianyang big-eared goats, belonging to a meat goat breed developed from crossbreeding Nubian goats with Chinese local breeds, were selected for the experiment. All animals (average BW of 14.45±1.21 kg, mean age of 2 months) were castrated by a licensed veterinarian and then housed in four adjacent pens, with six goats per pen. All goats were raised under standard conditions, fed twice a day (08:30 and 17:00) and given free access to water. Goats were slaughtered at 3, 5, 7 and 12 months of age. At each time point, six goats were slaughtered. Slaughters were performed by exsanguination at a commercial abattoir owned by Sichuan AUSDA Husbandry Development Co., Ltd. longissimus dorsi muscle and biceps femoris musclewere dissected, rinsed with RNase-free water and frozen in liquid nitrogen until RNA isolation.

RNA preparation and cDNA synthesis

Total RNA was extracted using Trizol (Invitrogen, Shanghai, China) according to the manufacturer’s protocol. RNA concentration was measured using a NanoDrop ND-1000 spectrophotometer (Agilent). Purity of the total RNA was determined by A260/280 and A260/230 ratios. RNA quality was assessed via agarose gel electrophoresis.

In our study, the A260/280 and A260/230 ratios ranged from 1.92 to 2.04 and from 2.01 to 2.12, respectively, which indicated the samples were of good quality. For each sample, 1μg of total RNA was reverse transcribed for cDNA synthesis using QuantiTect Reverse Transcription Kit (Qiagen) according to the manufacturer’s protocol. Prior to reverse transcription, genomic DNA was removed from the RNA with a specific gDNA Wipeout Buffer as described in the protocol.

Selection of reference genes and primer design

Eight candidate reference genes belonging to different functional classes, most of which have been evaluated as reference genes in ruminants (goat [23], sheep [24], cattle [25]) and other livestock species (pig [26], chicken [27]) as well as in humans [28], were tested for stability. HMBS, PPIB and RPLP0 were chosen because they had been previously reported to be the best reference gene in bovine muscle tissues [14,17,29]. B2M and YWHAZ were selected because they had been identified as suitable reference genes in porcine longissimus thoracis et lumborum (LTL) muscle and frozen whole blood of goats, respectively [23,30]. Additionally, GAPDH, 18S and ACTB were included because these genes have been utilized as reference genes in numerous studies [31,32,33].

Prior to primer design, a pretest study was conducted with primers identified by literature review for several of the reference genes. However, the performance of some of the primers was not as good as expected, and the optimal annealing temperatures of the primers are not identical. Thus, assays using these primers may work poorly and, furthermore, cannot run simultaneously on the same plate, making the results vulnerable to inter-assay variation. To reduce the inter-assay variation as much as possible for improved qPCR performance and to provide a set of species-specific primers for reference genes in goat, primer pairs were designed by primer5.0 software on the basis of caprine sequence information available in GenBank. Primer details are listed in Tables 1 and 2. Prior to RT-qPCR, conventional PCR and agarose gel electrophoresis were performed to test the gene-specific primers and verify the amplified products. All PCR products were sequenced and then aligned against the goat genome with the BLAST program to verify their identity. The specificity of the primers was evaluated by melting curve analysis. Only primers producing a single band of the expected size whose sequence matched published sequences exactly and exhibiting a unique peak during the dissociation step of the melting curve analysis were used for RT-qPCR.

[Figure omitted. See PDF.]

Table 1. Primers and relative information of reference and target genes.

https://doi.org/10.1371/journal.pone.0121280.t001

[Figure omitted. See PDF.]

Table 2. Description of reference and target genes analyzed in this study.

https://doi.org/10.1371/journal.pone.0121280.t002

Quantitative Real Time RT-PCR (qPCR)

Real-time PCR reactions were performed in 96-well plates with a Bio-Rad CFX96TM Real-time PCR Detection System. Reactions were carried out in a final volume of 25 μl containing 12.5μl 2× QuantiFast SYBR Green PCR Master Mix (Qiagen), 500nM of each primer, 1μl cDNA and 9.5μl DNase/RNase-free water. The cycling program was 95°C for 5 min to activate the polymerase followed by 40 cycles of denaturation at 95°C for 10s, annealing at 60°C for 20s and extension at 72°C for 20s. Melting curve analysis was conducted by heating samples from 65°C to 95°C with continuous fluorescent acquisition. All reactions were performed in triplicate for each cDNA sample. In addition, no template controls were tested for each primer pair. Standard curves were established using a 10-fold dilution series of purified PCR fragments as templates.

Data analysis

Three classic software programs, geNorm, NormFinder and BestKeeper, each with a unique advantage, were used to evaluate the reference genes.

geNorm defines the standard deviation (SD) of the expression ratio of two reference genes as a pairwise variation under the assumption that the genes tested are not co-regulated and that the expression ratio is identical across all samples. The gene-stability measure (M) is then calculated as the average pairwise variation of a specific gene compared with all of the other reference genes. Genes with the highest M values have the least stable expression. Eventually, all candidate reference genes are ranked by stepwise elimination of the genes with the highest M value. The outstanding advantage of geNorm is that it can identify the number of genes required for dependable normalization.

NormFinder is a model-based approach in which two types of variation are utilized to determine systematic error. This method takes into account both intragroup and intergroup variation; the latter is commonly overlooked in pairwise methods.

In contrast to geNorm and NormFinder, BestKeeper uses raw data instead of relative quantities. The estimate of reference gene stability is calculated based on the SD values of each gene. Then, the BestKeeper index is calculated from the geometric mean of the candidates C_T values for each specific sample. The most stable reference genes are the ones with the lowest SD values and highest coefficients of correlation with the BestKeeper index.

To evaluate the efficacy of the selection of reference genes for normalization, the relative change in the expression of target gene (PPARG) mRNA in skeletal muscle was calculated across different developmental stages using the 2^−ΔΔCT method [34], in which ΔΔCT = [(C_{T PPARG}−C_{T reference gene})_{month x}−(C_{T PPARG}−C_{T reference gene})_{month 12}]. Normalization was carried out with the least stable reference genes and the geometric mean of the most stable reference genes.

All relative quantification data were analyzed by one-way ANOVA with developmental stage as a factor using the PROC ANOVA procedure of the SAS software package, version 9.0. Duncan’s multiple-range test was performed for multiple comparisons when one-way ANOVA results were significant. Differences of P<0.05 were considered significant.

Results

Primer validation

For all primer pairs, the amplification efficiency ranged between 96% and 104% and the coefficient of determination (R²) varied from 0.998 to 1 (Table 2). Single peaks were observed in the analysis of melting curves, and agarose gel electrophoresis revealed a unique band at the expected size in the gel. In addition, the sequences of the amplified DNA fragments matched the reference and target gene sequences available in GenBank. Specific details of primer pair validation are presented in Supporting Information.

Expression profiling of candidate reference genes

The expression of eight candidate reference genes was assayed in samples of longissimus dorsimuscle and biceps femoris muscle. The distribution of C_T values is displayed in Fig. 1. Mean C_T values ranged from 20.19 to 26.56 in longissimus dorsi muscle and 18.48 to 26.41 in biceps femorismuscle and provide an overview of the variation in the expression of the candidate reference genes. Though average C_T values varied, the overall distribution of the C_T values was similar in longissimus dorsi muscle and biceps femoris muscle. This may be because both tissues are types of skeletal muscle.

[Figure omitted. See PDF.]

Fig 1. The distribution of reference gene expression levels in longissimus dorsi muscle (A) and biceps femoris muscle (B).

Means and medians are indicated by asterisks and lines, respectively. The boxes encompass the 25th to the 75th percentiles. Whisker caps denote the maximum and minimum values.

https://doi.org/10.1371/journal.pone.0121280.g001

GeNorm Analysis

The ranking of candidate genes was determined by geNorm software on the basis of their stability. In geNorm analyses, genes with M values below the cut-off value of 1.5 are considered to be stable, while those with M values above the threshold should not be used for data normalization. Fig. 2 shows the M values of all candidate genes during development for longissimus dorsi muscle, biceps femoris muscle and the combined group. All candidate genes evaluated in our study exhibited M values below 1.1, indicating acceptable expression stability. With the lowest M values, PPIB and HMBS had the most stable gene expression during the development of both longissimus dorsi muscle and biceps femoris muscle, whereas 18S exhibited significant variation and was the least stable gene. PPIB and HMBS were also the most stable genes when all samples were analyzed together. Similarly, 18S remained the least stable gene with the highest score in the combined group. In addition, the commonly used reference genes GAPDH and YWHAZ were not stably expressed in any sample. GAPDH was particularly variable, ranking just below 18S in biceps femoris muscle and the combined group. Hence, small differences existed in the geNorm rank of reference genes in the different muscle types.

[Figure omitted. See PDF.]

Fig 2. Gene expression stabilities and rankings of reference genes as calculated by geNorm.

Rank-order of gene expression stability is shown for longissimus dorsi muscle (A-1), biceps femoris muscle (B-1) and the combined group (C-1) according to the average expression stability values (M) for the reference genes from the least stable (left) to the most stable (right). Pairwise variation analysis (V) to determine the optimal number of reference genes for data normalization in longissimus dorsi muscle (A-2), biceps femoris muscle (B-2) and the combined group (C-2).

https://doi.org/10.1371/journal.pone.0121280.g002

The optimal number of reference genes needed for accurate normalization was determined by geNorm on the basis of the average pairwise variation (Vn/n+1) value. Our cut-off value was 0.15. Below this value, inclusion of an additional reference gene offers little further improvement to the data normalization. As depicted in Fig. 2, the value of V2/3 fell below 0.15 for all groups, implying that two reference genes—PPIB and HMBS—would be sufficient for accurate normalization across all samples.

NormFinder Analysis

NormFinder can identify the optimal reference gene with the most stable expression. As illustrated in Table 3, PPIB, having the lowest stability value in all groups, was again identified as the most stable gene, while the second-most stable gene was B2M in longissimus dorsi muscle and HMBS in biceps femoris muscle and the combined group. 18S exhibited the worst stability of the eight tested reference genes in all sample sets. Although ACTB and RPLP0 exhibited moderate variability, the commonly used reference genes GAPDH and YWHAZ tended to rank near the bottom of the candidate genes, which is in agreement with the results generated by geNorm analysis.

[Figure omitted. See PDF.]

Table 3. Ranking and stability values of reference genes calculated by NormFinder.

https://doi.org/10.1371/journal.pone.0121280.t003

BestKeeper Analysis

In BestKeeper, repeated pairwise correlation analysis was performed to determine the optimal reference genes. The premise of BestKeeper is that the SD values of any studied genes should not be higher than 1. The more stable genes have Pearson coefficients of correlation (R) close to1. YWHAZ, 18S and B2M exhibited substantial variation in all sample groups with SD values exceeding 1 (S3 Table). Thus, these genes were excluded and the analysis repeated. GAPDH was also excluded from the combined group because its SD value was 1. Table 4 displays the results for the remaining genes. PPIB and HMBS ranked as the two most stable reference genes in all sample sets with R-values ranging from 0.957 to 0.981. In addition, BestKeeper allows for the analysis of reference genes and target gene at same time. PPARG exhibited considerable variation in both longissimus dorsi muscle and biceps femoris muscle, with SD values of 1.37 and 1.12, respectively (S4 Table), demonstrating that PPARG is expressed variably during skeletal muscle development.

[Figure omitted. See PDF.]

Table 4. Reference genes stability calculated by BestKeeper based on C_T.

https://doi.org/10.1371/journal.pone.0121280.t004

Reference gene validation

To determine whether there is a significant difference between target gene expression normalized to the most stable reference genes and that normalized to the least stable reference genes, we analyzed the expression of PPARG at different developmental stages normalized to GAPDH, 18S and the geometric mean of PPIB and HMBS.

PPARG is an important regulator of many genes involved in adipogenesis [35] and is mainly expressed in intramuscular connective tissues in goat [36]. In vitro, up-regulation of PPARG allows for differentiation of pre-adipocytes into mature adipocytes in bovine adipose tissue, which is essential for promoting adipogenesis [37]. Additionally, some of the genes required for the terminal differentiation of adipocytes are under the control of PPARG [37]. PPARG is also involved in sustaining the adipocyte differentiation program and is an activator of fatty acid synthesis and storage [38], which is crucial to intramuscular fat deposition. Because of its critical role in intramuscular fat development, PPARG was chosen as the target gene for our study.

Our results suggest that the expression of PPARG varies during the development of skeletal muscle, which is consistent with the findings of our BestKeeper analysis. In both longissimus dorsi muscle and biceps femoris muscle, the expression pattern of PPARG across different developmental stages was similar when normalized to either PPIB&HMBS or GAPDH, whereas the expression pattern of PPARG was distinguishably different when normalized to 18S (Fig. 3). In longissimus dorsi muscle, compared with normalization against PPIB&HMBS, the most stably expressed genes, normalization against GAPDH led to a dramatic up-regulation of PPARG expression, with up to 2.8-fold changes observed at the peak of relative expression at 7 months. When PPARG expression was normalized to 18S, the peak of the expression occurred at 5 months, with up to 6.2-fold changes in expression compared with normalization against PPIB&HMBS. Similar results were detected in biceps femoris muscle.

[Figure omitted. See PDF.]

Fig 3. Effect of normalization on PPARG gene expression in skeletal muscles.

The expression of PPARG was normalized to the geometric mean of PPIB&HMBS, 18S or GAPDH expression and is shown relative to expression at 12 months for longissimus dorsi muscle (A) and biceps femoris muscle (B). Error bars indicate SD.

https://doi.org/10.1371/journal.pone.0121280.g003

Discussion

Appropriate reference genes are indispensible for accurate data normalization and thus reliable results in studies of gene expression. However, numerous studies choose reference genes without proper validation, picking genes that have been used for qPCR analyses in other investigations or arbitrarily selecting commonly used reference genes (eg: GAPDH, ACTB and 18S). Unfortunately, the use of such empirical reference genes is problematic, as increased variability may occur under certain experimental conditions [39,40]. Moreover, many published studies have shown that a number of commonly used reference genes, which are subjectively thought to be invariable irrespective of experimental conditions and have been used in published IMF research [41,42,43], are regulated and can vary considerably [44,45,46]. Thus, there is no doubt that reference genes appropriate for a certain context may differ from those suitable for other circumstances and that universal reference genes cannot exist. Thus, proper validation of reference genes should be performed prior to data normalization.

In this study, eight candidate reference genes have been evaluated for stability across a number of developmental stages in both longissimus dorsi muscle and biceps femoris muscle. Several algorithms with different working rationales, such as geNorm, NormFinder and BestKeeper, have been developed to facilitate selection of appropriate reference genes. geNorm concentrates on pairwise comparisons of genes with the highest degree of similarity in the expression profile, considering only intergroup variation. Thus, the geNorm algorithm may not work well when co-regulated genes are included. Unlike geNorm and Bestkeeper, which are based on pairwise comparisons, NormFinder takes both intragroup and intergroup variations into account in the stability value it calculates for each candidate gene, which causes it to be less sensitive to genes with correlated expression. A direct evaluation of expression variation is performed for every candidate gene in each sample with NormFinder, testing each gene in a given context [20]. NormFinder, unlike geNorm, does not determine the optimal number of reference genes. BestKeeper enables the user to analyze reference and target genes simultaneously. For increased reliability, more than one type of algorithm should be used to evaluate reference genes because algorithms have different strengths and weakness.

To prevent possible bias, three complementary algorithms (geNorm, NormFinder and BestKeeper) were used in our study. Despite the differences between the algorithms, the results obtained from all three algorithms were consistent. PPIB and HMBS were found to be the most stable genes during the development of longissimus dorsi muscle and biceps femoris muscle, regardless of whether the samples were analyzed separately or together. PPIB was previously validated as a reference gene for studies on IMF deposition in beef [17]. In addition, PPIA, which belongs to the same group of cyclophilin proteins as PPIB [47], has been shown to be the best reference gene for longissimus dorsi muscle development in pig [48]. Thus, identical results were observed in our research and in previous studies, demonstrating that PPIB is an appropriate normalization factor in skeletal muscle. HMBS outperformed all of the candidates aside from PPIB and is regarded as an optimal reference gene in bovine skeletal muscle [14] and multiple goat tissues [49]. The reference genes RPLP0 and B2M, which are the most stable genes under some conditions [15,23,50], displayed only moderate stability in our study, suggesting that the best reference genes for one environment may be unsuitable for another. Our geNorm results indicate that the two reference genes with the most stable expression in all groups (PPIB and HMBS) are sufficient for data normalization. However, if more than two reference genes are used, as Vandesompele recommends [19], RPLP0 may be a good choice. RPLP0 ranked as the third-most stable gene in all groups in our geNorm analysis and was identified as the optimal reference gene for research involving IMF development in cattle [29] and milk somatic cells in goat [23]. 18S and GAPDH were the least stable genes tested, ranking in the bottom half of the list in all analyses. As such, neither are appropriate reference genes for skeletal muscle studies.

To determine the impact normalization to different reference genes has on the expression pattern of target genes, the relative expression of PPARG was evaluated by normalizing it to several reference genes. The expression of PPARG in different muscle types, normalized against GAPDH and PPIB&HMBS, followed a similar pattern during development. However, an obvious overestimation of PPARG expression changes occurred when PPARG was normalized to GAPDH rather than PPIB&HMBS. When 18S was used in the normalization, significant fold changes in expression were detected at 5 and 7 months. Therefore, our comparison revealed substantial differences in the relative expression levels of PPARG when different genes were used for normalization, demonstrating the need for accurate validation of reference genes. Particular attention should be paid to commonly used reference genes, as inappropriate usage could lead to misinterpretation of results. In addition, when PPIB and HMBS were independently utilized for normalization, the relative abundance of PPARG was not altered markedly (data not shown). This maybe because there is little difference in their stability; both genes exhibited stable expression in all analyses by the three algorithms.

Given the differences observed in our comparison, we conclude that normalization against a single “commonly used reference gene” may lead to serious discrepancies in data interpretation. GAPDH was found to be unsuitable for normalization in skeletal muscle development in pig [48]; similar results were reported in goat [16] and cattle [51]. In addition, GAPDH was observed to be up-regulated and down-regulated under various conditions [7,9,52,53]. Likewise, a remarkable change in the expression of 18S occurred during osteoblast differentiation, with a continual increase in expression observed over time [6]. Given the afore mentioned influence of the environment on GAPDH expression and the distinct impact normalization to GAPDH had on the relative expression level of PPARG, we hypothesized that there are great fluctuations in the expression of GAPDH and 18S over the course of development. To test this hypothesis, relative expression of 18S and GAPDH was plotted by normalizing to the geometric average of PPIB&HMBS. As illustrated in Fig. 4, in longissimus dorsi muscle we observed a pronounced decrease in GAPDH expression between 3 and 7 months followed by a sharp increase in expression. The expression pattern was opposite to that of PPARG. The dramatic regulation of GAPDH expression over the course of skeletal muscle development may account for the overestimation of the changes in PPARG expression when GAPDH was used for normalization. An even greater modulation of expression was seen over time when expression was normalized to 18S because the relative expression level of 18S was significantly different at all time points except 3 and 7 months. Similar results were observed in biceps femoris muscle. Thus, these data support our hypothesis that GAPDH and 18S expression varies significantly during skeletal muscle development.

[Figure omitted. See PDF.]

Fig 4. Expression profiling of GAPDH and 18S in skeletal muscles at different developmental stages.

Gene expression data for both longissimus dorsi(A、C) and biceps femoris muscle (B、D)were normalized to the geometric mean of PPIB&HMBS and are relative to expression at 12 months following the 2^−ΔΔCT method [34]. Error bars depict SD. Superscript capital letters indicate significant differences (P<0.05).

https://doi.org/10.1371/journal.pone.0121280.g004

Conclusions

In our study, three different algorithms were used to analyze a panel of reference genes. All three algorithms gave consistent results. PPIB and HMBS proved to be the optimal reference genes in both longissimus dorsi muscle and biceps femoris muscle during development. GAPDH and 18S were the least stable genes, suggesting they are not suitable for expression data normalization in studies of goat muscle. Importantly, proper validation of reference genes is essential for reliable interpretation of data and should be performed prior to data normalization.

Supporting Information

[Figure omitted. See PDF.]

S1 Fig. Melting curves and standard curves of eight reference genes and a target gene.

https://doi.org/10.1371/journal.pone.0121280.s001

(DOCX)

S2 Fig. Agarose gel electrophoresis identification of gene-specific primers of reference genes for qPCR.

https://doi.org/10.1371/journal.pone.0121280.s002

(DOCX)

S1 Table. Sequencing results of PCR products from the amplification of primers of reference and target genes designed for this experiment.

https://doi.org/10.1371/journal.pone.0121280.s003

(DOCX)

S2 Table. Sequencing results of genes using BLASTN from NCBI against nucleotide collection (nr / nt).

https://doi.org/10.1371/journal.pone.0121280.s004

(DOCX)

S3 Table. Reference genes stability calculated by BestKeeper based on C_T.

https://doi.org/10.1371/journal.pone.0121280.s005

(DOCX)

S4 Table. Descriptive statistics of target gene (PPARG) analyzed by BestKeeper in Longissimus dorsi muscle and Biceps femoris muscle, respectively.

https://doi.org/10.1371/journal.pone.0121280.s006

(DOCX)

Author Contributions

Conceived and designed the experiments: WZ YW. Performed the experiments: WZ HL. Analyzed the data: WZ YL HL. Contributed reagents/materials/analysis tools: WZ YL HL. Wrote the paper: WZ YW.

Citation: Zhu W, Lin Y, Liao H, Wang Y (2015) Selection of Reference Genes for Gene Expression Studies Related to Intramuscular Fat Deposition in Capra hircus Skeletal Muscle. PLoS ONE 10(3): e0121280. https://doi.org/10.1371/journal.pone.0121280

References

1. Talpur FN, Bhanger MI, Sherazi ST. Intramuscular fatty acid profile of longissimus dorsi and semitendinosus muscle from Pateri goats fed under traditional feeding system of Sindh, Pakistan. Meat Sci. 2008; 80: 819–822. pmid:22063602

2. Gerbens F, Verburg FJ, Van Moerkerk HT, Engel B, Buist W, Veerkamp JH, et al. Associations of heart and adipocyte fatty acid-binding protein gene expression with intramuscular fat content in pigs. J Anim Sci. 2001; 79: 347–354. pmid:11219443

3. Wood JD, Enser M, Fisher AV, Nute GR, Sheard PR, Richardson RI, et al. Fat deposition, fatty acid composition and meat quality: A review. Meat Sci. 2008; 78: 343–358. pmid:22062452

4. Han X, Lu M, Chen Y, Zhan Z, Cui Q,Wang Y. Selection of reliable reference genes for gene expression studies using real-time PCR in tung tree during seed development. PLoS One. 2012; 7: e43084. pmid:22912794

5. Hruz T, Wyss M, Docquier M, Pfaffl MW, Masanetz S, Borghi L, et al. RefGenes: identification of reliable of reliable and condition specific reference genes for RT-qPCR data normalization. BMC Genomics. 2011; 12: 156. pmid:21418615

6. Stephens AS, Stephens SR, Morrison NA. Internal control genes for quantitative RT-PCR expression analysis in mouse osteoblasts, osteoclasts, macrophages. BMC Res Notes. 2011; 4: 410. pmid:21996334

7. Nygard AB, Jørgensen CB, Cirera S, Fredholm M. Selection of reference gene for gene expression studies in pig tissues using SYBR green qPCR. BMC Mol Biol. 2007; 8: 67. pmid:17697375

8. Modesto P, Peletto S, Pisoni G, Cremonesi P, Castiglioni B, Colussi S, et al. Evaluation of internal reference genes for quantitative expression analysis by real-time reverse transcription-PCR in somatic cells from goat milk. J Dairy Sci. 2013; 96: 7932–7944. pmid:24119819

9. Huth A, Vennemann B, Fracasso T, Lutz-Bonengel S, Vennemann M. Apparent versus true gene expression changes of three hypoxia-related genes in autopsy derived tissue and the importance of normalization. Int J Legal Med. 2013; 127: 335–344. pmid:23108458

10. Chaillou T, Malgoyre A, Banzet S, Chapot R, Koulmann N, Pugnière P, et al. Pitfalls in target mRNA quantification for real-time quantitative RT-PCR in overload-induced skeletal muscle hypertrophy. Physiol Genomics. 2011; 43: 228–235. pmid:21156833

11. Serrano M, Moreno-Sánchez N, González C, Marcos-Carcavilla A, Van Poucke M, Calvo JH, et al. Use of Maximum Likelihood-Mixed Models to select stable reference genes: a case of heart stress response in sheep. BMC Mol Biol. 2011; 12: 36. pmid:21849053

12. Brugè F, Venditti E, Tiano L, Littarru GP, Damiani E. Reference gene validation for qPCR on normoxia- and hypoxia- cultured human dermal fibroblasts exposed to USA: is β-actin a reliable normalize for photoaging studies? J Biotechnol. 2011; 156: 153–162. pmid:21963587

13. Wang YH, Byrne KA, Reverter A, Harper GS, Taniguchi M, McWilliam SM, et al. Transcriptional profiling of skeletal muscle tissue from two breeds of cattle. Mamm Genome. 2005; 16: 201–210. pmid:15834637

14. Pérez R, Tupac-Yupanqui I, Dunner S. Evaluation of suitable reference genes for gene expression studies in bovine muscular tissue. BMC Mol Biol. 2008; 9: 79. pmid:18786244

15. Bonnet M, Bernard L, Bes S, Leroux C. Selection of reference genes for quantitative real-time PCR normalization in adipose tissue, muscle, liver and mammary gland from ruminants. Animal. 2013; 7: 1344–1353. pmid:23552195

16. Najafpanah MJ, Sadeghi M, Bakhtiarizadeh MR. Reference genes selection for quantitative real-time PCR using RankAggreg method in different tissues of Capra hircus. PLoS One. 2013; 8: e83041. pmid:24358246

17. da Costa AS, Pires VM, Fontes CM, Mestre Prates JA. Expression of genes controlling fat deposition in two genetically diverse beef cattle breeds fed high or low silage diets. BMC Vet Res. 2013; 9:118. pmid:23767408

18. Lancaster PA, Sharman ED, Horn GW, Krehbiel CR, Starkey JD. Effect of rate of weight gain of steers during the stocker phase. Ⅲ. Gene expression of adipose tissues and skeletal muscle in growing-finishing beef cattle. J Anim Sci. 2014; 92: 1462–1472. pmid:24663156

19. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002; 3: RESEARCH0034.

20. Andersen CL, Jensen JL, Ørntoft TF. Normalization of real-time quantitative reverse transcription-PCR data: a model-based variation estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004; 64: 5245–5250. pmid:15289330

21. Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper—Excel-based tool using pair-wise correlations. Biotechnol Lett. 2004; 26: 509–515. pmid:15127793

22. Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, et al. The MIQE guidelines: Minimum informatio-n for publication of quantitative real-time PCR experiments. Clin Chem. 2009; 55: 611–622. pmid:19246619

23. Jarczak J, Kaba J, Bagnicka E. The validation of housekeeping genes as a reference in quantitative Real Time PCR analysis: Application in the milk somatic cells and frozen whole blood of goats infected with caprine arthritis encephalitis virus. Gene. 2014; 549:280–285. pmid:25068405

24. Peletto S, Bertuzzi S, Campanella C, Modesto P, Maniaci MG, Bellino C, et al. Evaluation of internal reference genes for quantitative expression analysis by real-time PCR in ovine whole blood. Int J Mol Sci. 2011; 12: 7732–7747. pmid:22174628

25. Bai WL, Yin RH, Zhao SJ, Jiang WQ, Yin RL, Ma ZJ, et al. Technical note: Selection of suitable reference genes for studying gene expression in milk somatic cell of yak (Bos grunniens) during the lactation cycle. J Dairy Sci. 2014; 97: 902–910. pmid:24342693

26. Cinar MU, Islam MA, Pröll M, Kocamis H, Tholen E, Tesfaye D, et al. Evaluation of situable reference genes for gene expression studies in porcine PBMCS in response to LPS and LTA. BMC Res Notes. 2013; 6: 56. pmid:23394600

27. Yue H, Lei XW, Yang FL, Li MY, Tang C. Reference gene selection for normalization of PCR analysis in chicken embryo fibroblast infected with H5N1 AIV. Virol Sin. 2010; 25: 425–431. pmid:21221921

28. Leal MF, Belangero PS, Cohen C, Figueiredo EA, Loyola LC, Pochini AC, et al. Identification of suitable reference genes for gene expression studies of shoulder instability. PLoS One. 2014; 9: e105002. pmid:25122470

29. Wang YH, Bower NI, Reverter A, Tan SH, De Jager N, Wang R, et al. Gene expression patterns during intramuscular fat development in cattle. J Anim Sci. 2009; 87: 119–130. pmid:18820161

30. McBryan J, Hamill RM, Davey G, Lawlor P, Mullen AM. Identification of suitable reference genes for gene expression analysis of pork meat quality and analysis of candidate genes associated with the trait drip loss. Meat Sci. 2010; 86: 436–439. pmid:20579813

31. Serrano C, Joy M, Serrano M, Rodellar C, Calvo JH. The effect of feeding system in the expression of genes related with fat metabolism in semitendinous muscle in sheep. Meat Sci. 2011; 89: 91–97. pmid:21543161

32. Waters SM, Kelly JP, O'Boyle P, Moloney AP, Kenny DA. Effect of level and duration of dietary n-3 polyunsaturated fatty acid supplementation on the transcriptional regulation of Delta9-desaturase in muscle of beef cattle. J Anim Sci. 2009; 87: 244–252. pmid:18791145

33. Zhao JP, Bao J, Wang XJ, Jiao HC, Song ZG, Lin H. Altered gene and protein expression of glucose transporter1 underlies dexamethasone inhibition of insulin-stimulated glucose uptake in chicken muscles. J Anim Sci. 2012; 90: 4337–4345. pmid:22859751

34. Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative C(T) method. Nat Protoc. 2008; 3: 1101–1108. pmid:18546601

35. Rosen ED, MacDougald OA. Adipocyte differentiation from the inside out. Nat Rev Mol Cell Biol. 2006; 7: 885–896. pmid:17139329

36. Li L, Jiang J, Wang L, Zhong T, Chen B, Zhan S, et al. Expression patterns of peroxisome proliferator- activated receptor gamma 1 versus gamma 2, and their association with intramuscular fat in goat tissues. Gene. 2013; 528: 195–200. pmid:23895798

37. Graugnard DE, Piantoni P, Bionaz M, Berger LL, Faulkner DB, Loor JJ. Adipogenic and energy metabolism gene networks in longissimus lumborum during rapid post-weaning growth in Angus and Angus × Simmental cattle fed high-starch or low-starch diets. BMC Genomics. 2009; 10:142. pmid:19335898

38. Varga T, Czimmerer Z, Nagy L. PPARs are a unique set of fatty acid regulated transcription factors controlling both lipid metabolism and inflammation. Biochim Biophys Acta. 2011; 1812: 1007–1022. pmid:21382489

39. Hu Q, Guo W, Gao Y, Tang R, Li D. Reference gene selection for real-time RT-PCR normalization in rice field eel (Monopterus albus) during gonad development. Fish Physiol Biochem. 2014; 40: 1721–1730. pmid:25079246

40. Zinzow-Kramer WM, Horton BM, Maney DL. Evaluation of reference genes for quantitative real-time PCR in the brain, pituitary, and gonads of songbirds. Horm Behav. 2014; 66: 267–275. pmid:24780145

41. Mannen H. Identification and utilization of genes associated with beef qualities. Anim Sci J. 2011; 82: 1–7. pmid:21269353

42. Wang L, Li L, Jiang J, Wang Y, Zhong T, Chen Y, et al. Molecualar characterization and different expression patterns of the FABP gene family during goat skeletal muscle development. Mol Biol Rep. 2015; 42: 201–207. pmid:25245957

43. Zhang HB, Wang ZS, Peng QH, Tan C, Zou HW. Effects of different levels of protein supplementary diet on gene expressions related to intramuscular deposition in early-weaned yaks. Anim Sci J. 2014; 85: 411–419. pmid:24450936

44. Thomas KC, Zheng XF, Garces Suarez F, Raftery JM, Quinlan KG, Yang N, et al. Evidence based selection of commonly used RT-qPCR reference genes for the analysis of mouse skeletal muscle. PLoS One. 2014; 9: e88653. pmid:24523926

45. Zhao DJ, Guo K, Kang L. Identification of condition-specific reference genes from microarray data for locusts exposed to hypobaric hypoxia. FEBS Open Bio. 2012; 2: 235–240. pmid:23650605

46. Kim I, Yang D, Tang X, Carroll JL. Reference gene validation for qPCR in rat carotid body during postnatal development. BMC Res Notes. 2011; 4: 440. pmid:22023793

47. Naoumov NV. Cyclophilin inhibition as potential therapy for liver diseases. J Hepatol. 2014; 61: 1166–1174. pmid:25048953

48. Feng X, Xiong Y, Qian H, Lei M, Xu D, Ren Z. Selection of reference genes for gene expression studies in porcine skeletal muscle using SYBR green qPCR. J Biotechnol. 2010; 150: 288–293. pmid:20887758

49. Zhang Y, Zhang XD, Liu X, Li YS, Ding JP, Zhang XR, et al. Reference gene screening for analyzing gene expression across goat tissue. Asian-Australas J Anim Sci. 2013; 26: 1665–1671. pmid:25049756

50. Vossaert L, O'Leary T, Van Neste C, Heindryckx B, Vandesompele J, De Sutter P, et al. Reference loci for RT-qPCR analysis of differentiating human embryonic stem cells. BMC Mol Biol. 2013; 14:21. pmid:24028740

51. Saremi B, Sauerwein H, Dänicke S, Mielenz M. Technical note: identification of reference genes for gene expression studies in different bovine tissues focusing on different fat depots. J Dairy Sci. 2012; 95: 3131–3138. pmid:22612949

52. Montero-Melendez T, Perretti M. Gapdh gene expression is modulated by inflammatory arthritis and is not suitable for qPCR normalization. Inflammation. 2014; 37: 1059–1069. pmid:24493325

53. Deloffre LA, Andrade A, Filipe AI, Canario AV. Reference genes to quantify gene expression during oogenesis in a teleost fish. Gene. 2012; 506: 69–75. pmid:22750322