From existing SAGE™ data [[21]; constructed using the field-grown seasonal samples described in the methods section], SAGE™ tags that were mapped correctly, annotated, and had moderate expression profiles across seasons (mean of 5-50 copies per virtual ryegrass shoot cell; i.e. similar expression levels to future target genes) were identified. From the list of the SAGE™ tags that met the above criteria, six genes involved at the pre-transcription stage (histone 3, H3), transcription (YT521-B-like protein family protein, YT521-B), translation (eukaryotic translation initiation factor 4 alpha, eIF4A), and in protein biosynthesis (eukaryotic elongation factor 1 alpha, eEF1A), modification (ubiquitin-conjugating enzyme, E2) and degradation (TAT-binding protein homolog 1, TBP-1) were selected as candidate reference genes.

As eEF1A is often identified as a stable reference gene [11,23-27], a more highly expressed eEF1A gene, eEF1A (h), was also selected from the SAGE™ data for testing, along with a strategy proposed by Martin et al. [25] in which transcripts from multigene eEF1A, eEF1A (m), would be identified. The seasonal expression profile and accession numbers of the candidate reference genes from the SAGE™ data are presented in Table 1.

All of the candidate reference genes were moderately abundant (median crossing point [Cp] values 26-31; Figure 1), with the exception of H3. Expression of H3 was generally very low in the tested samples (Cp values averaging 33.09 ± 2.51) or was not detected at all (n = 110 samples); therefore it was excluded from further analyses. The least variation in gene expression across all 442 tested samples was displayed by E2 (< 6 cycles), while eEF1A (m) was the most variable (8 cycles). Mean ± standard deviation of Cp values of the candidate reference genes from different tissue/cultivar/growth stage/treatment combinations are presented in Additional file 1.

Data were segregated into 10 different datasets for evaluation of gene expression stability using geNorm and NormFinder (Table 2). In geNorm, when all 442 perennial ryegrass samples were included (see Table 2 for sample information), the average expression stability (M) of the moderately expressed eEF1A (s) and multigene eEF1A (m) was least, and that of the more highly expressed eEF1A (h), was greatest. This suggests that expression of eEF1A (s) and eEF1A (m) was most stable and eEF1A (h) was least stable (Figure 2A). This was also the case in the dataset containing the 422 field-grown leaf and stubble samples collected at different growth stages following a range of defoliation treatments (Figure 2B). Average expression stability and ranking of the candidate reference genes in the remaining datasets are presented in Figures 2C-J.

The geNorm algorithm also calculated the pairwise variation V_n/V_n+1, which measured the effect of adding further reference genes on the normalisation factor, thus determining the optimal number of reference genes. Evaluation of all plant samples revealed a large decrease in the pairwise variation with the inclusion of a third, and then fourth reference gene (i.e. the differences between V2/3 and V3/4, and V3/4 and V4/5 in dataset A). When a fifth reference gene was added, the V values dropped below the proposed guideline of 0.15 (Figure 3). Thus, according to geNorm, five reference genes (eEF1A (s), eEF1A (m), TBP-1, E2 and YT521-B) are required for accurate normalisation. Pairwise variation within the field-grown samples collected at different growth stages following a range of defoliation treatments indicate that the four most stably expressed reference genes (eEF1A (s), eEF1A (m), TBP-1, and YT521-B) should be used for normalisation, while in the laboratory-grown samples use of the three most stable genes (eEF1A (s), TBP-1, and eIF4A) is sufficient (Figure 3).

The results of the NormFinder analysis are summarised in Table 3. When all 442 samples were included (dataset A), eEF1A (s) was identified as the most stable gene with an expression stability of 0.331, followed by YT521-B (0.429); eEF1A (h) was least stable (0.581). Out of the ten datasets, eEF1A (s) was most often identified as the most stable reference gene, while YT521-B was most often identified as the least stable reference gene.

NormFinder has the added ability of being able to estimate the variation between sample groups or treatments (as described further in the methods section). This function determines the best combination of two reference genes for normalisation. It also establishes whether normalisation using the two reference genes in combination will be more accurate than just using the most stable gene (i.e. if the stability value of the best two gene combination is lower than that of the most stable gene).

The dataset containing the 422 field-grown leaf and stubble samples collected at different growth stages following a range of defoliation treatments was the only dataset that contained sufficient replication to allow this full analysis (see Table 2 for replicate information). Analysing this dataset with or without the treatment groups identified did not affect the ranking of the genes, although the stability values were reduced with the inclusion of treatment groups (Table 3). With the treatment groups identified in this dataset, NormFinder selected eEF1A (s) as the most stably expressed single gene, with a stability value of 0.039. The best combination of two genes, eEF1A (s) and YT521-B, further reduced the NormFinder stability value to 0.030.

The expression levels of a target gene, chloroplast translational elongation factor Tu (EF-Tu, GenBank dbEST accession number GR517729), were used as an example to show the effect of using different reference genes for normalisation. The EF-Tu expression was normalised using three different strategies: 1) geometric average of the four most stably expressed reference genes selected by geNorm, 2) geometric average of the two most stably expressed reference genes selected by NormFinder, and 3) the least stably expressed gene according to both geNorm and NormFinder used alone.

In perennial ryegrass leaf tissue, there was no effect (P > 0.1) of defoliation frequency or severity, or any interaction between the defoliation treatments. There was, however, a significant (P < 0.001) interaction between leaf regrowth stage and the normalisation strategy used. Normalisation using the least stable reference gene (eEF1A (h)) led to over-estimation of the target gene following defoliation (0-leaf stage) and at the 1-leaf stage of regrowth compared with the geNorm strategy, and at the 1-leaf stage of regrowth compared with NormFinder (Figure 4). Although the geNorm and NormFinder strategies did differ in their estimation of the target gene at the 1- and 3-leaf stages of regrowth, the trend in transcript abundance throughout regrowth remained the same, in contrast with that displayed following normalisation using eEF1A (h).

Frozen perennial ryegrass tissues (callus, etiolated seedlings, inflorescence, leaf, stubble and root) from all conditions were ground independently in liquid nitrogen. Total RNA was extracted using the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol. Residual genomic DNA was removed by on-column DNAse I digestion, using the RNase-free DNase set (Qiagen), and mRNA was purified from total RNA using Dynabeads^® Oligo (dT)₂₅ (Invitrogen Dynal AS, Oslo, Norway). The mRNA concentration and purity were determined using a Nanodrop ND-1000 spectrophotometer (Nanodrop Technologies Inc., Wilmington, DE, USA); each mRNA sample was assayed twice and an average value determined.

Absence of genomic DNA contamination was confirmed by performing qRT-PCR on 0.1 ng of each of the mRNA samples using primers designed for the small subunit of ribulose-1,5-bisphosphate carboxylase gene (forward and reverse primers in 5'→3' direction are GAGGAGTCCGGCAAGGCATAA and TATGCTTTTACATGTAGCCGGTTC, respectively).

Messenger RNA (10 ng) was reverse transcribed to produce cDNA using the Transcriptor First Strand cDNA Synthesis Kit (Roche Diagnostics, Mannheim, Germany) with anchored-oligo (dT)₁₈ primers in total reaction volumes of 20 μl. The two cDNA samples from each season (one pre-grazed and one post-grazed) were bulked together, resulting in four seasonal samples (autumn, winter, spring and summer). All cDNA samples were diluted 100-fold with PCR-grade water.

Primer pairs were designed to amplify a large portion of the 3' untranslated region (3'UTR) of the candidate reference genes and the target gene using Primer3 software [http://primer3.sourceforge.net/; accessed 2007/2008; [56]] and are described in Table 4. To ensure maximum specificity and efficiency during PCR amplification a stringent set of criteria was used for primer design [57]. This included predicted melting temperatures (T_m) of 58-60°C, primer lengths of 19-24 nucleotides, guanine-cytosine contents of 36-58%, and PCR amplicon lengths of 111-168 base pairs. All primers were custom-ordered from a commercial supplier (Invitrogen, Auckland, New Zealand).

The qRT-PCR were performed in 384-well plates with a LightCycler^® 480 real-time PCR instrument (Roche Diagnostics) using the LightCycler^® 480 SYBR Green I Master kit. The reaction set-up was performed on the epMotion^® 5075LH automated liquid handling system (Eppendorf, Hamburg, Germany). Reactions were performed in triplicate and contained 5 μl SYBR Green I Master, 2 μl PCR-grade water, 2 μl cDNA, and 0.5 μl of each of the 10 μM forward and reverse gene-specific primers in a final volume of 10 μl. In addition, each plate contained no-template controls and two calibrator samples required for normalisation of the target gene (one leaf and one stubble sample collected at the 3-leaf stage, i.e. immediately before defoliation).

The reactions were incubated at 95°C for 5 min to activate the FastStart Taq DNA polymerase, followed by 45 cycles of 95°C for 10 sec, 60°C for 10 sec, and 72°C for 8 sec. The specificity of the PCR reaction was confirmed with a heat dissociation protocol (from 60°C to 95°C) following the final PCR cycle. This ensured the resulting fluorescence originated from a single PCR product, and did not represent primer dimers formed during PCR or a non-specific product. Amplification of a single product of expected size was verified by gel electrophoresis on a 1.5% agarose gel (Sigma-Aldrich, St Louis, MO, USA) and ethidium bromide staining.

LightCycler^® 480 software (version 1.5; Roche Diagnostics) was used to collect the fluorescence data. PCR efficiencies were calculated using the equation E = 10^-1/slope on a standard curve generated using a tenfold dilution series of one sample (leaf and stubble) over three dilution points that were measured in triplicate.

The mean, standard deviation (SD), and coefficient of variation (CV) of the raw triplicate qRT-PCR values within each plate were determined. Samples whose CV were greater than 1.5% were inspected; a reaction was considered an outlier if one of the triplicate reactions deviated by more than 1 SD from the mean and it was excluded from analysis. Samples were repeated if exclusion of one of the reactions still did not result in a CV <1.5%.

Two publicly available software tools, geNorm [17] and NormFinder [22] were used to evaluate gene expression stability. Both tools require the transformation of Cp values to linear scale expression quantities. Using the LightCycler^® software, Cp values were converted into quantities via the standard curve with the Absolute Quantification Fit Points method, and both measures were exported into Microsoft Excel.

To ensure that data from different plates were comparable, the quantities for each gene were then normalised to the quantity of the 1/100 dilution from the standard curve dilution series that was run on each plate. For example, the dilution series on the first plate for eIF4A resulted in an average quantity of 0.0102 for the triplicate 1/100 dilution. Following absolute quantification, the average quantities on the second, third and fourth plates for the triplicate 1/100 dilution were 0.0102, 0.0098 and 0.0108. Normalisation factors for each plate were calculated by dividing 0.0102 by the average quantity for each plate, resulting in normalisation factors of 1.0020, 1.0404, and 0.9501 for the second, third and fourth plates, respectively. The quantities for each of the samples on each plate were then multiplied by the calculated normalisation factor.

The quantities were then imported into the two software tools, geNorm (version 3.5) and NormFinder, which were used as described by Vandesompele et al. [17] and Andersen et al. [22], respectively. NormFinder has the added ability of being able to estimate the variation between sample groups (i.e. treatments). In most of our datasets there was not enough samples per treatment to fully utilise this function; therefore, only a simple analysis of gene stability was carried out which identified the most stable gene in each dataset. The dataset containing the 422 field-grown leaf and stubble samples collected at different growth stages following a range of defoliation treatments did however contain sufficient replication to allow a full analysis, enabling the combination of the two best reference genes to be determined.

Normalised ratios of the target gene EF-Tu in perennial ryegrass leaf tissue collected following different defoliation management were calculated from the LightCycler^® Relative Quantification Software (Roche Diagnostics) using the formula:

This formula provides an efficiency-corrected relative quantification, normalised to a calibrator sample (perennial ryegrass leaf at the 3-leaf stage), where TS is the concentration of the target gene in a sample, RS is the concentration of the reference gene in a sample, TC is the concentration of the target gene in the calibrator sample and RC is the concentration of the reference gene in the calibrator sample.

The EF-Tu expression was normalised using three different strategies: 1) geometric average of the four most stably expressed reference genes selected by geNorm, 2) geometric average of the two most stably expressed reference genes selected by NormFinder, and 3) the least stably expressed gene according to both geNorm and NormFinder used alone.

The normalised ratios of the target gene in leaf tissue were log₁₀-transformed before statistical analysis. Results were analysed using mixed models with a compound symmetry covariance structure for the repeated measurements through time in GenStat 11.1 [58]. The statistical model included defoliation frequency and severity, time (leaf regrowth stage), the normalisation strategy, and all possible interactions between defoliation frequency, severity, time, and strategy. Following statistical analysis, averages were back-transformed, and results expressed as the mRNA transcript abundance relative to the leaf calibrator sample.