- Research article
- Open Access
Validation of reference genes for quantitative expression analysis by real-time RT-PCR in Saccharomyces cerevisiae
BMC Molecular Biology volume 10, Article number: 99 (2009)
Real-time RT-PCR is the recommended method for quantitative gene expression analysis. A compulsory step is the selection of good reference genes for normalization. A few genes often referred to as HouseKeeping Genes (HSK), such as ACT1, RDN18 or PDA1 are among the most commonly used, as their expression is assumed to remain unchanged over a wide range of conditions. Since this assumption is very unlikely, a geometric averaging of multiple, carefully selected internal control genes is now strongly recommended for normalization to avoid this problem of expression variation of single reference genes. The aim of this work was to search for a set of reference genes for reliable gene expression analysis in Saccharomyces cerevisiae.
From public microarray datasets, we selected potential reference genes whose expression remained apparently invariable during long-term growth on glucose. Using the algorithm geNorm, ALG9, TAF10, TFC1 and UBC6 turned out to be genes whose expression remained stable, independent of the growth conditions and the strain backgrounds tested in this study. We then showed that the geometric averaging of any subset of three genes among the six most stable genes resulted in very similar normalized data, which contrasted with inconsistent results among various biological samples when the normalization was performed with ACT1. Normalization with multiple selected genes was therefore applied to transcriptional analysis of genes involved in glycogen metabolism. We determined an induction ratio of 100-fold for GPH1 and 20-fold for GSY2 between the exponential phase and the diauxic shift on glucose. There was no induction of these two genes at this transition phase on galactose, although in both cases, the kinetics of glycogen accumulation was similar. In contrast, SGA1 expression was independent of the carbon source and increased by 3-fold in stationary phase.
In this work, we provided a set of genes that are suitable reference genes for quantitative gene expression analysis by real-time RT-PCR in yeast biological samples covering a large panel of physiological states. In contrast, we invalidated and discourage the use of ACT1 as well as other commonly used reference genes (PDA1, TDH3, RDN18, etc) as internal controls for quantitative gene expression analysis in yeast.
Real-time PCR technology has recently reached a level of sensitivity, accuracy and practical simplicity allowing its use as a routine bioinstrumentation for pathogen detection, single nucleotide polymorphism and gene expression analysis [1–4]. In particular for the latter application, several controls are needed to ensure the integrity of each step along the process  and therefore, to obtain reliable and accurate results. This process includes RNA extraction (yield, integrity, DNA contamination), efficiency of the reverse transcription and PCR steps, amount of RNA added into the reaction, etc. While the quantitative RT-PCR is technically robust, the normalization procedure to correct sample-to-sample variation remains a critical and challenging problem of this method [1, 4, 6, 7]. Several procedures have been suggested based on physical parameters, such as volume or cell number, but these methods are either impractical or unreliable due to the heterogeneity of biological samples. Some authors favour an internal control strategy, which uses an alien RNA molecule that is artificially incorporated into the biological sample . As an example, Liu & Slininger proposed a set of universal external RNA calibrators for microbial mRNA expression analysis . In spite of these initiatives, the most common practice is to normalize to either total RNA amount, ribosomal RNA or to a single internal reference gene termed HouseKeeping gene (HSK). Several mathematical models have been developed that calculate the "relative" mRNA expression changes of a target gene with respect to an HSK. The "2ΔΔCt" approach  is the most popular application in quantitative RT-PCR, but it assumes optimal and identical PCR efficiencies of target and reference genes. Violation of this rule results in a systematic bias that either underestimates or overestimates the initial copy numbers. This problem can be bypassed by adjusting for PCR efficiency, which can be estimated using many approaches [10–12] that can be separated into three groups : serial dilutions, individual graph analysis based on the rate of fluorescence accumulation within the exponential region, or mathematical model fitting. Whatever the method employed for determining the PCR efficiency, accurate relative quantification implies that the expression of the reference gene is perfectly stable in the sample set.
It is empirically assumed that housekeeping genes fulfil the criterion of unregulated expression independent of the experimental condition. However, some evidence shows that these genes are regulated to some extent, reinforcing the idea that there is no universal reference gene whose expression level remains constant whatever the conditions . Since even small variations of an internal control could lead to non-reliable expression data, it is critical to validate that the expression of reference genes is stable prior to their use for normalization in real time RT-PCR analysis. To overcome the "circular problem" of evaluating the expression stability of a candidate gene if no reliable measure is available to normalize the candidate , Vandesompele and colleagues  developed a statistical algorithm termed geNorm. Their strategy relied on (i) a careful selection of a set of genes that display minimal variation across different biological conditions, and (ii) normalization of the genes of interest to the geometric mean of a minimal, albeit optimal number of the selected genes. The strength of using geometric averaging is in smoothing the individual variation of the expression value of a single reference gene, which can lead to large errors of normalized data in samples of interest . Other statistical algorithms were also proposed, e.g. Bestkeeper , which allows including up to ten genes of interest in the analysis, or Normfinder  that is apparently less sensitive toward coregulation of the candidate reference genes.
Real-time RT-PCR is commonly used to validate microarray-generated data [16, 17]. ACT1 and RDN18 are among the most frequently used reference genes in S. cerevisiae studies, because the expression of these genes has been considered relatively stable under the conditions investigated. However, only two recent papers showed the stability of ACT1 expression and some other standard reference genes to normalize the expression of genes involved in central carbon metabolism during short-term glucose pulse , or during the rehydration process in active dry yeast . With the notable exception of these works, we could not find any study dedicated to the selection and validation of suitable reference genes in S. cerevisaie, contrary to other fungal models such as the pathogenic yeast Candida albicans , and the fungi Metarhizium anisopliae  and Aspergillus niger . Therefore, the purpose of the present work was to identify a robust set of reference genes for growth phase-related mRNA profiling in the yeast Saccharomyces cerevisiae. From public microarray datasets, we selected a set of potential reference genes that exhibited minimal variation among various conditions. The most stable subset of internal controls, which gave rise to a robust normalization factor, was then applied to quantify expression of genes involved in glycogen metabolism in response to changing growth conditions, and in a mutant defective in TPS1 which encodes the trehalose-6P synthase subunit .
Cell samples were regularly harvested from the yeast cultures, and only samples from key physiological states were selected and used for mRNA quantification by real time RT-PCR assays (Figure 1 and Additional files 1 and 2). These physiological states were defined from the macrokinetic growth parameters and reserve carbohydrates profiles as follows: the exponential - respiro fermentative- phase (EP); the diauxic shift (DS), which corresponds to the time when the sugar has just been exhausted from the medium while glycogen shows a transient peak of accumulation; the post-diauxic or purely respiratory phase (PDS), which corresponds to the re-assimilation of fermentation products and to a second phase of glycogen accumulation; and finally, the stationary phase (SP) when cells are starved for carbon nutrient. Contrary to the wild type strain behaviour, the tps1 mutant significantly mobilized glycogen during the PDS phase, and consequently more samples were taken up during this period to better characterize the physiological state of this mutant in this growth phase (see also additional file 2 for details of the sampling). In total, we carried out six independent yeast cultures (Figure 2). The wild type KT strain was grown on glucose (sample set B) as the basic and reference growth condition . It was also cultivated on galactose (set C), which was used as the growth control condition for tps1 mutant since this mutant strain cannot grow on glucose . Finally, cultures on galactose of the CEN.PK strain and the corresponding tps1 mutant were made in duplicates (sets D & E for the wild type CEN.PK and sets H & I for the tps1 mutant).
Expression level and stability of candidate reference genes
As stated in the introduction, accurate normalization requires reference genes whose expression changes are negligible under the investigated conditions. Candidate genes were therefore identified using public microarray datasets from De Risi et al  and Gasch et al , because the culture conditions reported in these studies were the closest to our experimental setup. We selected eight potential reference genes based on the stability of their expression during growth on glucose (genes highlighted in bold in Table 1), taking care that these genes belong to different functional categories to minimize the risk of coregulation. The remaining genes listed in Table 1, i.e. ACT1, PDA1, RDN18, IPP1 and TDH3, were also included in the list since they are traditionally used as single reference genes in expression studies by Northern blots or real time RT-PCR.
Transcription profiling using real-time RT-PCR assays was then performed with these 13 candidate genes, in samples from the 6 cultures. We first analyzed transcript abundance of these genes in the different samples by direct comparison of their cycle threshold (Ct), assuming equal Ct for equal transcript number since all RT-PCR reactions were performed with equal quantity of total RNA. As can be seen in Figure 3, most of the selected genes presented Ct values that spanned from 20 to 30 cycles, while Ct values from RDN18 and TDH3 were clearly lower. For RDN18, these values centered around 8 cycles with a very low dispersion. The glycolytic gene TDH3 was also highly expressed as indicated by Ct values around 17 cycles, but it exhibited rather high dispersion over the growth phases and culture conditions as indicated by large whiskers of the box and many outliers. The Ct of the remaining selected genes showed a reasonable dispersion, with expression levels of ALG9, KRE11, TAF10, TFC1 and UBC6 exhibiting smaller variation than that of ACT1, HEM2, IPP1 or PDA1.
These raw Ct data were then analyzed using geNorm to identify the most suitable candidate genes. For each independent culture (e.g. sample set B, Figure 2 left panel), or pool of several cultures (e.g. set A that combined samples from cultures B & C), the 13 genes were ranked according to their gene expression stability measure "M" (Figure 2, right panel, and additional file 3). All genes presented an M value below 1.5, which is the default limit for acceptable expression stability as defined by Vandesompele et al . Another advantage of geNorm is to provide the optimal number of reference genes required for accurate normalization. This number is obtained by calculating the Pairwise variation values (V(n/n+1)) between each combination of sequential normalization factors (NF) (Figure 2, right column). Vandesompele and coworkers  recommended a cut-off value at 0.15, below which the inclusion of an additional gene does not result in a significant improvement of the normalization. According to this criterion, TAF10 and UBC6 turned out to be sufficient as internal controls to normalize expression levels from samples taken from growth on glucose (set A, V2/3 = 0.098), whereas 5 genes were required for normalizing gene expression from the data sets F and G (V5/6 = 0.150 and 0.141, respectively). According to the recommendation of Vandesompele et al , we always used a minimum of three of the most stably expressed genes to calculate the normalization factor. From this analysis, ACT1, IPP1 and TDH3 were excluded from the set of selected genes for normalization as they always ranked among the worst candidates. In contrast, ALG9, TAF10, TFC1, UBC6 and to a lesser extent KRE11 turned out to be the most stable genes in the culture conditions tested in this study (Figure 2 right panel, additional file 3).
To further support this conclusion and the suitability of this set of genes to serve as a reference in a broader panel of experimental conditions, we examined gene stability of these candidate reference genes by using data from the entire microarray datasets from the SGD server (approx. 30 experiments), which altogether correspond to several hundred different experimental conditions. As can be seen in Figure 4, genes like ALG9, TAF10, TFC1, UBC6 presented a significantly higher number of experiments with a log2 ratio close to zero as compared to ACT1. This microarray survey analysis indicated that our initially selected genes exhibited very little expression change over a wide range of experimental conditions. Therefore, this set of genes, ALG9, TAF10, TFC1 and UBC6 should be preferentially used to calculate normalization factors in quantitative RT-PCR expression analysis in the yeast S. cerevisiae.
Impact of reference gene selection on expression ratio values
The strength of GeNorm to select the most suitable reference genes was demonstrated by comparing normalized data calculated from different subsets of potential HSKs. As is shown in Figure 5, the use of a normalization factor based on the geometric mean of expression levels of UBC6, TAF10 and ALG9 (NF(UBC 6, TAF 10, ALG 9)) yielded expected expression patterns of the glycogen metabolic genes GSY2 and GPH1 during growth of the KT strain on glucose (see last results section for more details). Ratios of expression values were almost identical using the following 3 best genes based on geNorm classification (NF(TFC 1, KRE 11, FRP 2), compare grey and hatched bars) or applying the normalization factor calculated from the six best genes together (NF(6 best), not plotted). These results showed that, at least in this condition, any subset of three genes among the most stably expressed candidates was sufficient to calculate robust NF and to normalize expression of Genes of Interest (GOIs).
The advantage of using validated genes for normalization was further analyzed comparing expression results after normalization (NF(UBC 6, TAF 10, ALG 9), etc) to those obtained by using non-validated reference genes, e.g. ACT1 (NF(ACT 1)) or a combination of ACT1, PDA1 and IPP1 (NF(ACT 1, PDA 1, IPP 1)). As it could be expected from the coregulation of these three genes during growth on glucose (Figure 5), identical expression data were obtained with NF(ACT 1)and NF(ACT 1, PDA 1, IPP 1)as normalization factors (compare empty and black dots for RDN18 and GOIs normalized expression data). In contrast, a strong discrepancy between normalization to ACT1 and normalization to validated genes was observed in biological samples collected in the post-diauxic phase (#3 and #4) and in the stationary phase (#5). This strong deviation could be explained by the drop of ACT1 mRNA as well as that of transcripts of other HSK genes (IPP1 and PDA1) during these growth phases (Figure 5). To better visualize the advantage of using normalization to validated genes, data from sample set B were reported on a scatter plot (Figure 6), comparing data normalized to NF(TFC 1, KRE 11, FRP 2)and NF(ACT 1), respectively, to those normalized to NF(UBC 6, TAF 10, ALG 9)(Figure 6). A regression coefficient close to one (R2 = 0.997) was calculated for NF(TFC 1, KRE 11, FRP 2)versus NF(UBC 6, TAF 10, ALG 9). In contrast, the coefficient was extremely low for NF(ACT 1)versus NF(UBC 6, TAF 10, ALG 9)(R2 = 0.598), mainly due to expression data from PDS and SP samples.
Similar analyses were carried out by using different biological situations, as for instance, in a biological set combining samples from cultures on glucose and galactose to analyze the influence of the carbon source (Figure 7, set A), or from wild type and tps1 mutants to analyze the impact of the mutation on the expression data (Figure 8, set K). Again, discrepancies in ratios of expression values calculated by normalization to the "best reference genes" and to ACT1, respectively, were evident in samples collected from yeast cultures entering the diauxic shift. The difference was even more pronounced with late stationary phase samples, as the difference could reach almost 10-fold between the two procedures (see Figure 7 &8). This discrepancy was visualized in the scatter plots presented in Figure 9, which report a larger range of ratio values than those in Figures 7 &8. As expected, data from (NF(TAF 10, FRP 2, ALG 9)) versus (NF(UBC 6, TFC 1, KRE 11)) aligned with a good regression coefficient (panel A, R2 = 0.916), while NF(ACT 1) did not correlate at all with NF(UBC 6, TFC 1, KRE 11) as shown by the worst regression coefficient of the study (panel B, R2 = 0.124). Altogether, these results demonstrated the benefit of using multiple selected genes instead of a single, non-validated gene (e.g. ACT1) for accurate and reliable data normalization.
Application to quantitative expression analysis of genes involved in glycogen metabolism
To test the robustness of this subset of selected reference genes, we analyzed the transcriptional regulation of genes involved in glycogen metabolism in this yeast. It has been reported that large variations of reserve carbohydrate content are associated with coordinated transcriptional regulation of the cognate genes in response to changing growth conditions or under various genetic contexts [23, 24, 28–33]. When examining raw Ct values from GPH1 (glycogen phosphorylase) and GSY2 (glycogen synthase) in the complete dataset (Figure 3), the very long whiskers of the boxes and numerous outliers confirmed the high variability of the expression of these two genes. On the contrary, the SGA1 gene encoding the vacuolar amylo-1,4 -1,6 glucosidase  exhibited much lower dispersion of Ct values, indicating smaller expression change than GSY2 and GPH1 under our growth conditions. Using the 3 best reference genes for data normalization (NF(UBC 6, TAF 10, ALG 9)), we confirmed the induction of GSY2 and GPH1 between the exponential phase (S#1, calibrator sample) and the entry into the diauxic shift (S#2), and found remarkable expression ratios close to 20 for GSY2 and almost 100 for GPH1 (Figure 5). Moreover, this normalization procedure allowed us to show that the expression of these two genes dropped immediately after the diauxic shift, to return to the initial level in stationary phase for GSY2, or close to it for GPH1, while ACT1 normalization indicated stable and high expression of these two genes all along these growth phases. In contrast to GSY2 and GPH1, the expression of SGA1 showed a modest increase when cells entered the diauxic shift to reach a 3-fold activation in the stationary phase. We also analyzed for the first time transcriptional patterns in galactose-grown cells (Figure 7). Unexpectedly, the expression of GPH1 and GSY2 was already very high in the exponential phase (S#6) as compared to cultures on glucose (S#1), and it did not further increase as cells entered the diauxic shift on this carbon source, whereas glycogen accumulated with a kinetic almost similar to that on glucose (see additional files 1 &2). The expression of these two genes then dropped during the post diauxic phase to reach levels even lower than on glucose in the stationary phase. In contrast to GPH1 and GSY2, SGA1 expression was not affected by the carbon source. Finally, expression patterns of these three genes on galactose were the same in the CEN.PK genetic background (Figure 8) as in the KT strain (Figure 7).
The loss of TPS1 function had a strong impact on the glycogen accumulation pattern on galactose, as it caused hyper-accumulation of the polymer at the end of the exponential phase, and also promoted its rapid and sustained degradation during post-diauxic and stationary phases (Figure 1, and in additional file 2). Therefore, to get a preliminary idea on how this mutation could alter the glycogen kinetics at the transcriptional level, we quantified the expression of GSY2, GPH1 and SGA1 during growth on galactose of the WT and its tps1 derivative strain. It is shown in Figure 8 that the expression pattern of GPH1 or GSY2 was identical in both the wild type and in the mutant strains. The only noticeable difference was for SGA1 gene whose expression was already increased in the exponential phase in a tps1 mutant (S#14) as compared to the wild type strain (S#10, calibrator).
In the yeast Saccharomyces cerevisiae, the microarray datasets available on the Saccharomyces Genome Database website now represent a vast treasure-trove of reference genes suitable for gene expression normalization. We therefore used the Gene Expression Connection tool  to search for a set of stably expressed genes in growth dynamics [26, 27]. Other approaches have recently been proposed for selection of internal controls, with statistical analysis of large microarray datasets [36, 37]. Nevertheless, as stated by the authors [36, 37], these in silico searches for stable internal controls must be accompanied by lab-bench work to verify that selected candidate genes are reliable for normalization in a specific experimental context. This is what we actually performed in the present work. Out of 13 genes analyzed in this study, i.e. 8 functionally unrelated genes selected from the microarray datasets together with 5 standard reference genes (e.g. ACT1, PDA1 etc), we identified ALG9, TAF10, TFC1, UBC6 and to a lesser extent KRE11, as the most stable genes in our experimental conditions. Another very important result from this study was the observation that geometric averaging of any subset of three genes among the six most stable genes led to very similar normalization factors, therefore highlighting the robustness of our gene selection. This conclusion was further supported by the weak expression change of this subset of genes as revealed by a survey of microarray datasets from the SGD server, i.e. in approximately 30 large scale transcriptomic studies, which altogether correspond to several hundreds of different samples. Therefore, ALG9, TAF10, TFC1 and UBC6 genes are the most pertinent reference genes, not only for growth phase related mRNA profiling in S. cerevisiae, but more generally for quantitative gene expression analysis with samples that cover a large panel of physiological and metabolic states.
Probably because of the tedious experimental validation of a suitable set of reference genes, the common practice by many authors was to use ACT1, PDA1, TDH3 or RDN18 as a single reference gene for normalization in S. cerevisiae expression studies, assuming a steady state level of expression for these so-called housekeeping genes. However, the validity of ACT1 as an internal standard has been already questioned  and PDA1 was preferred to ACT1 for normalization in Northern blot analysis since stationary phase samples showed a more significant drop of ACT1 mRNA. Still, by the use of Northern blot analysis, it was indeed demonstrated that ACT1 is a representative gene whose transcription is typically repressed following the shift from logarithmic growth to stationary phase [39, 40]. As preliminary molecular clues on how down regulation of this gene occurs, these authors showed that topoisomerase I has a regulatory role in the transcriptional repression of most of the genes following the diauxic shift and in the stationary phase . The second work reported a major role of the RNA polymerase II subunit RPB4, which permitted appropriate transcriptional responses during stress, including nutrient stress that accompanies entry into stationary phase . Global transcriptome analysis  also showed that ACT1 mRNA levels did decrease significantly during growth in a glucose rich medium, a pattern that was confirmed by using absolute quantification by real time RT-PCR . In our study, using the geometric averaging of multiple selected reference genes for relative quantification, we also found a significant drop of ACT1 transcripts and of other frequently used genes like PDA1 (E1 alpha subunit of the pyruvate dehydrogenase complex) and IPP1 (cytoplasmic inorganic pyrophosphatase) when cells proceed from the exponential phase of growth to the stationary phase. Therefore, as already mentioned by Monje-Cajas et al , the use of ACT1 and related transcripts would seriously over-estimate (approx. 10-fold) expression levels of genes of interest in stationary phase, leading to erroneous conclusions. Moreover, the genes encoding glycolytic enzymes, for example TDH3 (glyceraldehyde-3-phosphate dehydrogenase), were amongst the first yeast genes to be isolated. Because of their high expression levels, their promoters have been widely used to construct yeast expression vectors [42–47] and as model systems to study transcription . Nevertheless, as details of the organisation of glycolytic promoters have emerged, it has become clear that these "simple HouseKeeping genes" actually have sophisticated molecular mechanisms controlling their expression . As reviewed in this latter reference, some glycolytic enzymes were induced by glucose while others such as ENO1 and TDH3 appeared to be constitutively expressed irrespective of the carbon source (glucose or other sugars versus non-fermentable carbon sources). Our results, i.e. the large dispersion of Ct values from TDH3 in our experimental conditions (Figure 3) and the fact that this gene ranked amongst the worst candidate reference genes (Additional File 3), somehow contradicted this conclusion and do not encourage the use of TDH3 as a constitutive, internal reference gene for reliable gene normalization. Finally, the gene RDN18 encoding the 18srRNA could have been accepted for normalization, but the very strong expression that is several orders of magnitude above the mean expression of tested reference genes and GOIs, together with the fact that RDN18 does not encode mRNA precluded its use as an internal reference gene.
As an application of this normalization procedure with carefully selected reference genes, we reinvestigated in a quantitative manner the transcriptional response of genes implicated in glycogen metabolism, since this physiological event is an interesting hallmark during long term yeast cultures with significant transcriptional remodeling and huge variations of the polysaccharide content [23, 24]. This study confirmed known transcriptional induction patterns of genes encoding glycogen phosphorylase (GPH1) and glycogen synthase (GSY2) that paralleled the accumulation of glycogen between the exponential phase and the diauxic shift on glucose [23, 24], with a notable induction ratio of approx. 100-fold for GPH1 and 20-fold for GSY2. For the first time, this study also provided expression data of these genes on galactose. Despite similar growth pattern and glycogen accumulation kinetic as compared to glucose, the growth on galactose radically changed the expression pattern of these two genes. Transcripts levels were already very high in the exponential phase as compared to glucose, and we could not observe any upregulation at the entry into the diauxic-shift. Therefore, the glycogen accumulation that was observed at the end of the exponential phase during growth on galactose probably came from a concomitant activation of Gsy2p and the inactivation of Gph1p by protein phosphatases mediated dephosphorylation events . With respect to the vacuolar amylo-glucosidase (SGA1), the large-scale transcriptional study from Gash and co-workers showed some transcriptional activation of this gene during long term yeast cultures on glucose [27, 49]. Wang and co-workers  also provided indirect evidence of transcriptional activation of SGA1 in stationary phase, showing a significant role of SGA1 gene deletion on glycogen accumulation pattern in this phase of growth. Their result suggested that SGA1 may not be strictly sporulation specific, but could be activated under starvation, like in the late stationary phase. Our study supported this assertion and showed that the expression of SGA1 increased during the post-diauxic growth, with maximal 3-fold activation in stationary phase in wild type strains. On the other hand, the causal relationship between higher expression of SGA1 and faster glycogen degradation in tps1 derivative strains as compared to the WT remains to be investigated.
A set of putative internal control genes for real-time RT-PCR analysis was selected from public microarray datasets. Using geNorm software, we validated that ALG9, TAF10, TFC1, UBC6 and to a lesser extent KRE11, turned out to be the most stable genes under all conditions investigated. Geometric averaging of their expression data was then applied to smooth individual variation and to calculate robust normalization factors, which allowed for the demonstration that the use of a single reference gene like ACT1 could lead to erroneous expression data. This set of reference genes was carefully tested in a context of large heterogeneity of samples (different physiological states during long term S. cerevisiae cultures, different carbon sources and genetic contexts) and applied to explore quantitatively the transcriptional regulation of genes involved in glycogen metabolism in this yeast. This study brought new insights into the transcriptional control of GSY2, GPH1 and SGA1 during long term growth on glucose and galactose, suggesting a potential role of SGA1 in the management of glycogen storage in tps1 cells. To summarize, this work provides a set of pertinent reference genes that should be used for validation of expression data from microarrays experiments and more generally for reliable real time RT- PCR analysis in yeast.
Yeast strains, growth conditions and sampling
The Saccharomyces cerevisiae strains, CEN.PK113-7D  and its tps1 derivative , and JLP48-3B (KT1112 context ), were grown at 30°C in a synthetic minimal medium containing 0.17% (w/v) yeast nitrogen base (DIFCO), 0.5% (w/v) ammonium sulfate and 2% (w/v) galactose (YNGal) or glucose (YNGlu). The use of prototroph strains did avoid amino acid complementation of the medium. The pH was adjusted to 5.0 with succinic acid and sodium hydroxide. Cell growth was followed by measurement of OD (600 nm) during at least 10 days. For real independency of duplicates, shake-flasks cultures were performed neither simultaneously, nor from the same inoculums. The residual extracellular carbon source was quantified by HPLC. Intracellular glycogen and trehalose were measured as described in . Yeast samples for real-time PCR analysis (approx. 108 cells) were centrifuged (3,000 rpm, 4°C, 3 min), and the cell pellets were immediately frozen in liquid nitrogen and stored at -80°C until RNA extraction.
Total RNA extraction
Frozen cells were mechanically disrupted using a ball mill (Mikro-Dismembrator S; B. Braun Biotech International). Total RNA was extracted using the RNeasy mini kit (Qiagen). To eliminate genomic DNA contamination, an additional DNase treatment was performed according to the RNeasy kit instruction with the RNase-free DNase set (Qiagen). The extracted RNA was quantified using the Bioanalyzer 2100 with the RNA 6000 Nano LabChip kit (Agilent) and the ND-1000 UV-visible light spectrophotometer (NanoDrop Technologies). As another preliminary quality control assay, the absence of contaminant genomic DNA in RNA preparations was verified using RNA as a template in real-time PCR assays (minus RT control, i.e. RNA not reverse-transcribed to cDNA).
Oligonucleotides for real-time PCR (Table 1) were designed using Beacon Designer 2.0 software (PREMIER Biosoft International), which included a BLAST analysis against S. cerevisiae Genome sequence for specificity confidence, and analysis using the Mfold server to avoid positioning on risky secondary structures.
One microgram of total RNA was reverse-transcribed into cDNA in a 20 μL reaction mixture using the iScript cDNA synthesis kit (Bio-Rad). The cDNA levels were then analyzed using the MyIQ real-time PCR system from Bio-Rad. Each sample was tested in duplicate in a 96-well plate (Bio-Rad, CA). The reaction mix (25 μL final volume) consisted of 12.5 μL of iQ SYBR Green Supermix (Bio-Rad), 2.5 μL of each primer (250 nM final concentration), 2.5 μL of H2O, and 5 μL of a 1/10 dilution of the cDNA preparation. The absence of genomic DNA in RNA samples was checked by real-time PCR before cDNA synthesis (minus RT control). A blank (No Template Control) was also incorporated in each assay. The thermocycling program consisted of one hold at 95°C for 4 min, followed by 40 cycles of 10 s at 95°C and 45 s at 56°C. After completion of these cycles, melting-curve data were then collected to verify PCR specificity, contamination and the absence of primer dimers.
The PCR efficiency of each primer pair (Eff) was evaluated by the dilution series method using a mix of sample cDNAs as the template. Briefly, it was determined from standard curves using the formula 10(-1/slope). For the calculations, the base of the exponential amplification function was used (e.g. 1.94 means 94% amplification efficiency). Relative expression levels were determined with efficiency correction , which considers differences in primer pair amplification efficiencies between target and reference genes, and results in a more reliable estimation of the "real expression ratio" than the 2ΔΔCt method . Expression data and associated technical errors on duplicates were calculated using the gene expression module of the BIORAD iQ5 software, which follows models and error propagation rules outlined in the geNorm manual.
Data analysis using geNorm
The stability of mRNA expression of tested reference genes was evaluated by using the geNorm VBA applet for Microsoft Excel . This program calculates the gene expression stability measure "M" for a potential reference gene as the average pair-wise variation for that gene with all other tested genes. Then it ranks genes considering that those with the lowest M value have the most stable expression. Finally, it determines the optimal number of genes for an accurate normalization by calculating the pair-wise variation (V(n/n+1)) between two normalization factors, namely NFn (normalization factor based on the geometric mean of the n most stable genes) and NFn+1: if Vn/n+1 is superior to 0.15 as a cut-off value, one could consider that the (n+1)th gene has a significant effect on normalization quality and should preferably be included for calculation of a reliable normalization factor. Authors of geNorm nevertheless recommend the minimal use of the three most stable internal control genes for calculation of the normalization factor (NF3), and stepwise inclusion of more control genes until the (n+1)th gene has no significant contribution to the newly calculated normalization factor.
Bustin SA, Benes V, Nolan T, Pfaffl MW: Quantitative real-time RT-PCR--a perspective. J Mol Endocrinol. 2005, 34 (3): 597-601.
Kubista M, Andrade JM, Bengtsson M, Forootan A, Jonak J, Lind K, Sindelka R, Sjoback R, Sjogreen B, Strombom L, et al: The real-time polymerase chain reaction. Mol Aspects Med. 2006, 27 (2-3): 95-125.
VanGuilder HD, Vrana KE, Freeman WM: Twenty-five years of quantitative PCR for gene expression analysis. Biotechniques. 2008, 44 (5): 619-626.
Wong ML, Medrano JF: Real-time PCR for mRNA quantitation. Biotechniques. 2005, 39 (1): 75-85.
Nolan T, Hands RE, Bustin SA: Quantification of mRNA using real-time RT-PCR. Nat Protoc. 2006, 1 (3): 1559-1582.
Ho-Pun-Cheung A, Cellier D, Lopez-Crapez E: [Considerations for normalisation of RT-qPCR in oncology]. Ann Biol Clin (Paris). 2008, 66 (2): 121-129.
Huggett J, Dheda K, Bustin S, Zumla A: Real-time RT-PCR normalisation; strategies and considerations. Genes Immun. 2005, 6 (4): 279-284.
Liu ZL, Slininger PJ: Universal external RNA controls for microbial gene expression analysis using microarray and qRT-PCR. J Microbiol Methods. 2007, 68 (3): 486-496.
Livak KJ, Schmittgen TD: Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods. 2001, 25 (4): 402-408.
Cikos S, Bukovska A, Koppel J: Relative quantification of mRNA: comparison of methods currently used for real-time PCR data analysis. BMC Mol Biol. 2007, 8: 113-
Pfaffl MW: Quantification strategies in real-time PCR. A-Z of quantitative PCR. Edited by: Bustin S. 2003, 5: 912-La jolla, CA, USA: International University Line (IUL)
Rebrikov DV, Trofimov D: [Real-time PCR: approaches to data analysis (a review)]. Prikl Biokhim Mikrobiol. 2006, 42 (5): 520-528.
Andersen CL, Jensen JL, Orntoft TF: Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004, 64 (15): 5245-5250.
Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, Speleman F: Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002, 3 (7): RESEARCH0034-
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 (6): 509-515.
Dallas PB, Gottardo NG, Firth MJ, Beesley AH, Hoffmann K, Terry PA, Freitas JR, Boag JM, Cummings AJ, Kees UR: Gene expression levels assessed by oligonucleotide microarray analysis and quantitative real-time RT-PCR -- how well do they correlate?. BMC Genomics. 2005, 6 (1): 59-
Provenzano M, Mocellin S: Complementary techniques: validation of gene expression data by quantitative real time PCR. Adv Exp Med Biol. 2007, 593: 66-73.
Stahlberg A, Elbing K, Andrade-Garda JM, Sjogreen B, Forootan A, Kubista M: Multiway real-time PCR gene expression profiling in yeast Saccharomyces cerevisiae reveals altered transcriptional response of ADH-genes to glucose stimuli. BMC Genomics. 2008, 9: 170-
Vaudano E, Costantini A, Cersosimo M, Del Prete V, Garcia-Moruno E: Application of real-time RT-PCR to study gene expression in active dry yeast (ADY) during the rehydration phase. Int J Food Microbiol. 2009, 129 (1): 30-36.
Nailis H, Coenye T, Van Nieuwerburgh F, Deforce D, Nelis HJ: Development and evaluation of different normalization strategies for gene expression studies in Candida albicans biofilms by real-time PCR. BMC Mol Biol. 2006, 7: 25-
Fang W, Bidochka MJ: Expression of genes involved in germination, conidiogenesis and pathogenesis in Metarhizium anisopliae using quantitative real-time RT-PCR. Mycol Res. 2006, 110 (Pt 10): 1165-1171.
Bohle K, Jungebloud A, Gocke Y, Dalpiaz A, Cordes C, Horn H, Hempel DC: Selection of reference genes for normalisation of specific gene quantification data of Aspergillus niger. J Biotechnol. 2007, 132 (4): 353-358.
Francois J, Parrou JL: Reserve carbohydrates metabolism in the yeast Saccharomyces cerevisiae. FEMS Microbiol Rev. 2001, 25 (1): 125-145.
Parrou JL, Enjalbert B, Plourde L, Bauche A, Gonzalez B, Francois J: Dynamic responses of reserve carbohydrate metabolism under carbon and nitrogen limitations in Saccharomyces cerevisiae. Yeast. 1999, 15 (3): 191-203.
Gancedo C, Flores CL: The importance of a functional trehalose biosynthetic pathway for the life of yeasts and fungi. FEMS Yeast Res. 2004, 4 (4-5): 351-359.
De Risi JL, Iyer VR, Brown PO: Exploring the metabolic and genetic control of gene expression on a genomic scale. Science. 1997, 278 (5338): 680-
Gasch AP, Spellman PT, Kao CM, Carmel-Harel O, Eisen MB, Storz G, Botstein D, Brown PO: Genomic expression programs in the response of yeast cells to environmental changes. Mol Biol Cell. 2000, 11 (12): 4241-4257.
Hwang PK, Tugendreich S, Fletterick RJ: Molecular analysis of GPH1, the gene encoding glycogen phosphorylase in Saccharomyces cerevisiae. Mol Cell Biol. 1989, 9 (4): 1659-1666.
Ni HT, LaPorte DC: Response of a yeast glycogen synthase gene to stress. Mol Microbiol. 1995, 16 (6): 1197-1205.
Parrou JL, Enjalbert B, Francois J: STRE- and cAMP-independent transcriptional induction of Saccharomyces cerevisiae GSY2 encoding glycogen synthase during diauxic growth on glucose. Yeast. 1999, 15 (14): 1471-1484.
Rowen DW, Meinke M, LaPorte DC: GLC3 and GHA1 of Saccharomyces cerevisiae are allelic and encode the glycogen branching enzyme. Mol Cell Biol. 1992, 12 (1): 22-29.
Teste MA, Enjalbert B, Parrou JL, Francois JM: The Saccharomyces cerevisiae YPR184w gene encodes the glycogen debranching enzyme. FEMS Microbiol Lett. 2000, 193 (1): 105-110.
Unnikrishnan I, Miller S, Meinke M, LaPorte DC: Multiple positive and negative elements involved in the regulation of expression of GSY1 in Saccharomyces cerevisiae. J Biol Chem. 2003, 278 (29): 26450-26457.
Clancy MJ, Smith LM, Magee PT: Developmental regulation of a sporulation-specific enzyme activity in Saccharomyces cerevisiae. Mol Cell Biol. 1982, 2 (2): 171-178.
Ball CA, Jin H, Sherlock G, Weng S, Matese JC, Andrada R, Binkley G, Dolinski K, Dwight SS, Harris MA, et al: Saccharomyces Genome Database provides tools to survey gene expression and functional analysis data. Nucleic Acids Res. 2001, 29 (1): 80-81.
Faccioli P, Ciceri GP, Provero P, Stanca AM, Morcia C, Terzi V: A combined strategy of "in silico" transcriptome analysis and web search engine optimization allows an agile identification of reference genes suitable for normalization in gene expression studies. Plant Mol Biol. 2007, 63 (5): 679-688.
Lee S, Jo M, Lee J, Koh SS, Kim S: Identification of novel universal housekeeping genes by statistical analysis of microarray data. J Biochem Mol Biol. 2007, 40 (2): 226-231.
Wenzel TJ, Teunissen AW, de Steensma HY: PDA1 mRNA: a standard for quantitation of mRNA in Saccharomyces cerevisiae superior to ACT1 mRNA. Nucleic Acids Res. 1995, 23 (5): 883-884.
Choder M: A general topoisomerase I-dependent transcriptional repression in the stationary phase in yeast. Genes Dev. 1991, 5 (12A): 2315-2326.
Choder M, Young RA: A portion of RNA polymerase II molecules has a component essential for stress responses and stress survival. Mol Cell Biol. 1993, 13 (11): 6984-6991.
Monje-Casas F, Michan C, Pueyo C: Absolute transcript levels of thioredoxin- and glutathione-dependent redox systems in Saccharomyces cerevisiae: response to stress and modulation with growth. Biochem J. 2004, 383 (Pt 1): 139-147.
Rosenberg S, Coit D, Tekamp-Olson P: Glyceraldehyde-3-phosphate dehydrogenase-derived expression cassettes for constitutive synthesis of heterologous proteins. Methods Enzymol. 1990, 185: 341-351.
Kingsman SM, Cousens D, Stanway CA, Chambers A, Wilson M, Kingsman AJ: High-efficiency yeast expression vectors based on the promoter of the phosphoglycerate kinase gene. Methods Enzymol. 1990, 185: 329-341.
Hinnen A, Buxton F, Chaudhuri B, Heim J, Hottiger T, Meyhack B, Pohlig G: Gene expression in recombinant yeast. Bioprocess Technol. 1995, 22: 121-193.
Mumberg D, Muller R, Funk M: Yeast vectors for the controlled expression of heterologous proteins in different genetic backgrounds. Gene. 1995, 156 (1): 119-122.
Nacken V, Achstetter T, Degryse E: Probing the limits of expression levels by varying promoter strength and plasmid copy number in Saccharomyces cerevisiae. Gene. 1996, 175 (1-2): 253-260.
Graham IR, Chambers A: Constitutive expression vectors: PGK. Methods Mol Biol. 1997, 62: 159-169.
Chambers A, Packham EA, Graham IR: Control of glycolytic gene expression in the budding yeast (Saccharomyces cerevisiae). Current Genetics. 1995, 29 (1): 1-
Wang Z, Wilson WA, Fujino MA, Roach PJ: Antagonistic controls of autophagy and glycogen accumulation by Snf1p, the yeast homolog of AMP-activated protein kinase, and the cyclin-dependent kinase Pho85p. Mol Cell Biol. 2001, 21 (17): 5742-5752.
van Dijken JP, Bauer J, Brambilla L, Duboc P, Francois JM, Gancedo C, Giuseppin ML, Heijnen JJ, Hoare M, Lange HC, et al: An interlaboratory comparison of physiological and genetic properties of four Saccharomyces cerevisiae strains. Enzyme Microb Technol. 2000, 26 (9-10): 706-714.
Guillou V, Plourde-Owobi L, Parrou JL, Goma G, Francois J: Role of reserve carbohydrates in the growth dynamics of Saccharomyces cerevisiae. FEMS Yeast Res. 2004, 4 (8): 773-787.
Parrou JL, Francois J: A simplified procedure for a rapid and reliable assay of both glycogen and trehalose in whole yeast cells. Anal Biochem. 1997, 248 (1): 186-188.
Pfaffl MW: A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001, 29 (9): e45-
This work was partially supported by ANR Blanc n° 05-2-42128 and by Genopole grant (2005-2007) to JMF. We are indebted to Dr. S. Sokol for his help with statistical analysis, Dr. Thomas Walther for critical reading, and Christine Rettew (Drexel University) for copyediting the manuscript.
JLP, MAT and JMF conceived the study. MAT carried out SGD microarray datasets analysis, primer design and their validation; MAT and MD carried out the cultures, samples treatment and qPCR analysis. JLP and JMF wrote the manuscript and all authors read and approved the final version.
Marie-Ange Teste, Manon Duquenne contributed equally to this work.
Electronic supplementary material
Additional file 1: Growth curve and glycogen content of WT strain on glucose and galactose. Growth (cells) and glycogen content during cultures of KT strain on glucose (set B from Figure 2) and galactose (set C from Figure 2). Cell samples (red dots) analyzed by real-time RT-PCR and sample numbering (S# followed by red numbers in the blue area). Cells (OD600), Glycogen (μg eq.glucose/OD unit). (PDF 41 KB)
Additional file 2: Growth curve and glycogen content of WT and tps1 strains. Growth (cells) and glycogen content during cultures of CEN.PK strains on galactose, WT (set D from Figure 2) and tps1 (set H from Figure 2). Legend as in Additional file 1. (PDF 54 KB)
Additional file 3: Ranking of reference genes according to their expression stability. Compiled data from all sample sets. For each set (A to M), genes are ranked from the least stable (left) to the most stable (right). The two most stable genes cannot be ranked in order. Gene expression stability value (Upper panel) as a function of gene name (lower panel). (PDF 16 KB)
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.