Selection of transcriptome datasets

Global analysis of transcriptome datasets

In order to analyse these multiple global rice transcriptome data in a comparable way, all arrays were normalised in the same way (materials and methods) and present/absent calls were determined MAS5.0 normalisation. The genome was defined as the 57,302 probesets targeted to Oryza sativa, thus the 81 probesets designed for the bacterial/phage controls were not included. The normalised data from all 373 microarrays (Table 1), representing 136 biological samples were collated together and a probeset was considered to be expressed in a particular tissue/sample if all replicates for every sample showed statistically significant present calls (p < 0.05). This cut-off method has previously been used as a way of present/absent determination [10, 21]. Using this principle, the expression for each probeset across all microarrays could be determined. Nearly eight thousand (7,922) probesets were detected in all 373 microarray samples, thereby fulfilling the first criterion for defining reference genes (Figure 1).

Selection of reference genes

The GC-RMA normalised data for all microarrays with publically available CEL files (331 microarrays; Table 1) was used to calculate the mean, standard deviation (SD) and coefficient of variance (CV; CV = SD/mean) for all 7,922 probesets, where a low CV is indicative of lower variation. This was followed by selection process undertaken to determine which of these genes were suitable as reference genes (Figure 1). Only 151 of the 7,922 probesets were defined as stably expressed across the developmental, stress and/or entire dataset (Figure 1).

In order to visualise the expression of these 151 probesets, the log₂ normalised values were hierarchically clustered and as expected, stable expression profiles were observed across the tissue development, stress and hormone microarray experiments (Figure 2A). Only 2 of these genes, LOC_Os07g02340.1 and LOC_Os03g05290.1, have been previously identified as stably expressed, with the former gene identified in a previous rice study [22], and the latter based on orthology with an Arabidopsis reference gene [7] (Figure 2B, red asterisk). A selection of 12 genes that showed stable expression across the microarrays (Figure 2B) were analysed further by QRT-PCR (Genes 1-12; Table 2). These 12 genes were selected on the basis of their CV and included; 2 transcripts with the lowest CV calculated across the stress microarray set (Genes 1-2), 2 transcripts with the lowest CV across the developmental set (Genes 3-4), 3 transcripts with the lowest CV across the entire microarray set (Genes 5-7) and the remaining 4 genes were randomly selected from the 66 probesets with low CV values (</=0.35; Genes 8-12) from the entire microarray set (Figure 1 and 2B; Table 2).

Closer analysis of these 12 genes reveals that the genes encoding, a 3-phosphoinositide-dependent protein kinase-1 (LOC_Os06g48970.1) and a nucleic acid binding protein (LOC_Os06g11170.1) showed stable, moderate expression levels across the stress microarray set (Genes 1-2 in Table 2; Figure 2B). While the genes encoding a tumor protein homolog (LOC_Os11g43900.1) and translation initiation factor SUI1 (LOC_Os07g34589.1) showed stable expression across the developmental and entire microarray sets respectively (Genes 4, 7 in Table 2; Figure 2B). As would be expected, it can be seen that many of these stably expressed genes are involved in core cellular functions such as mRNA splicing and translation initiation (Genes 1-12 denoted by blue asterisks in Figure 2A; 2B; Table 2).

In order to compare the reference genes defined in this study with the expression of some genes defined as stably expressed in these previous studies [7, 10, 22], 14 genes commonly used reference genes were visualised in the same way across the microarrays (Figure 2C) and the mean, SD and CV for each was also calculated (Genes 13-26; Table 2). It can be seen that there is a large amount of variation in transcript abundance for many of the previously proposed stably expressed genes as well as the typical reference genes, such as those encoding Actin and ubiquitin (Figure 2C; high CVs in Table 2). It is particularly evidenced that beta-tubulin transcript expression is variable under bacterial and parasite infection respectively (Figure 2C). Although the heatmap visualisation of the expression for the nucleotide tract-binding protein (LOC_Os03g25980.1) and TIP41-like protein (LOC_Os03g55270.1) appears unchanging (Figure 3 - top 2 genes), it can be seen that the CVs for both of those genes is over 0.4 indicating a higher level of variation in expression (Table 2).

Validation of reference genes in quantitative RT-PCR in tissue and stress samples

In order to confirm stable expression of the reference genes identified in this study primers were designed to 26 genes,12 stably expressed genes identified in this study and 14 previously suggested reference genes (Table 1, Additional file 1, Table S1). The stability of transcript abundance of these genes was analysed by QRT-PCR across 15 different samples from a variety of developmental (dry seed, imbibed seed, leaf and roots from young and old plants) and stress treated tissues (shoots from cold treated and heat treated young seedlings over time; Materials and methods). High quality total RNA was isolated from these samples and reverse transcribed to generate cDNA. The same cDNA pool from each of the samples was used to measure the transcript abundance by QRT-PCR, with melt curve analysis for each gene confirming primer specificity.

The geNORM v3.5 software was used to analyse the expression stability for the reference genes analysed by QRT-PCR from the 12 tissue samples (Additional file 1, Table S1) [13]. This software allows calculation of a gene stability measure (M) value for all the genes analysed, where genes with the lowest M value shown the most stable expression (Figure 3A). Authors of the geNORM software suggest using the 3 most stable genes (3 lowest M values) as the most appropriate reference genes [13]. It can be seen that even when commonly or previously suggested reference genes and the novel reference genes from this study are analysed together, all 3 of the most stable genes are the novel reference genes identified in this study (Figure 3A). It is important to note that this M value is only calculated based on data from the limited number of samples that were analysed by QRT-PCR, thus not representing the wide variety of tissues/treatments analysed by microarrays. Therefore, in order to visualise the variation in expression across in the microarrays in parallel, the CV values for each gene was also plotted with the M values, where a lower CV value indicates greater stability. In this way, the most stable genes were identified as those with both low M and CV values. In this combined analysis, the 12 genes chosen all outperformed previously used reference genes, particularly in terms of having a lower CV (Figure 3A), the genes indicated with a black diamond all had lower CV values as indicated by the bar graph, with a gene encoding a nucleic acid binding protein (LOC_Os06g11170.1) apparently the most stable (Figure 3A).

To further test the stability of the reference genes defined in this study, the expression of the 12 novel reference genes defined in this study were analysed independently by geNORM for the samples analysed by QRT-PCR (Figure 3B and 3C). Overall, it can be seen that the most stable genes had low M values as well and low CV values, indicating stable expression (Figure 3B). Furthermore, the geNORM pair-wise analysis to determine the number of control genes recommended for use in normalisation [13], revealed that 2 or even one gene is stable enough for accurate normalisation, however 2 genes is recommended for more robust normalisation (V < 0.15; Additional file 2, Figure S1) [13]. Using QRT-PCR analysis, we also compared the expression of these 12 reference genes to 3 heat or cold responsive genes including, an Apetala type transcription factor (AP2), a heat shock responsive factor (HSF-82) and alternative oxidase (AOX) over time under i) cold or ii) heat conditions respectively (Figure 3C). It can be observed that under cold treatment, all 12 reference genes show very stable expression over time (Figure 3Ci). Similarly, despite slight variation of some genes under heat conditions, it is evidenced that overall, these genes are also stably expressed over time following heat treatment (Figure 3Cii). In addition, the observed induction of AP2 and HSF-82 under cold and heat treatment, confirmed the success of the respective treatments (Figure 3C). Furthermore, comparison of this induction (at 3 h) to the induction observed from the analogous microarray data, showed that normalisation of the QRT-PCR data using the reference genes defined in this study resulted in comparable increases to those seen using the microarray data (Figure 3C).

Comparison to previous studies and other expression platforms

A large-scale study of reference genes in Arabidopsis revealed superior reference genes using Affymetrix microarray data [10]. Using the Inparanoid orthologue output [23] for Arabidopsis and rice, it was seen that only 15 rice orthologues of the 30 novel Arabidopsis reference genes were also expressed across all the microarrays in this study and 3 of these were randomly selected for further analysis by QRT-PCR (Genes 24-26; Table 2). Notably, only 1 gene (LOC_Os03g05290.1) encoding an aquaporin TIP protein, was seen to be stably expressed i.e. one of the 151 stably expressed in this study (red asterisk only; Figure 2A). It may be noted that the overall CV values are higher in this study compared to the CV values calculated in the Arabidopsis study [10]. The main reason for this is likely to be due to significant differences in the variability of the input data from both studies. That is, the Arabidopsis reference gene study used microarray data generated from only 7 studies using a large number of microarrays each e.g. 237 microarrays in the single developmental study [10], whilst this study involved analysis of microarrays from 20 studies carried out in different laboratories, using between 4 and 60 microarrays in each.

Previous studies in rice have examined reference genes using QRT-PCR analysis, however these only involved analysis of a small number of commonly used reference genes such as Actin, Actin1, alpha and beta tubulin, polyubiquitin, ubiquitin, GAPDH and elongation factor 1 in up to 25 samples, under a limited range of conditions [7, 20]. Analysis of these genes in the context of this study (Genes 13-20; Table 2) revealed that some of these were not detected as expressed in one or more tissue/stress microarray experiments, notably, this included Actin1 (LOC_Os05g36290.1; Gene 14 in Table 2) which was not expressed in all 3 biological replicates of the semi apical meristem (GSE6901) (Figure 2C). Similarly, a recent study in rice defined a set of 248 stably expressed genes across 40 developmental tissues that were analysed using Yale/BGI oligonucleotide microarrays [22]. Only 61 of these genes were found to be expressed across all the microarrays analysed in this study, nevertheless 3 of these were randomly selected for further analysis by QRT-PCR (Genes 21-23; Table 2). Notably, one of the 61 genes (LOC_Os07g02340.1) encoding an "expressed protein" was also found to fulfil all the criteria outlined in Figure 1, and showed stable expression across all the samples analysed in the present study (Gene 8 in Table 2; denoted by red and blue asterisk in Figure 2A and 2B).

In order to test the robustness of expression stability for the 12 reference genes identified in this study, two different approaches were undertaken. Firstly the expression patterns of these 12 genes were examined on other expression platforms, specifically the BGI/Yale oligonucleotide and Agilent microarray platforms. Overall a stable expression pattern was observed for all genes examined, with the most stable expression particularly evidenced for LOC_Os11g43900.1, LOC_Os03g46770.1 and LOC_Os07g02340.1 using the Yale oligonucleotide microarrays (Figure 4A). Notably, the latter gene was also grouped within the 248 stably expressed genes defined previously identified [22], thus complementing the identification of this gene in the presented study. Similarly, the 12 reference genes identified in this study were also examined for changes in expression following infection with hemibiotrophic fungus Magnaporthe oryzae [24]. In this study, Agilent Arrays (G4138A) were used for global transcriptomic analysis following infection [24]. The expression of all 12 genes were not found to significantly differ (Students t-test, p < 0.01) following infection (Figure 4B). However, given that this experiment involved stress treatment; AP-2, HSF-82 and AOX expression were also examined following infection and it was observed that AOX was significantly up-regulated (p < 0.01) following infection (Figure 4B). AOX is a known stress responsive gene [25]. Thus the reference genes defined are stable even under biotic stress stimulation, in addition to the abiotic treatments carried out as described above.

Analyses of all publically available rice microarrays

To compile the entire publically available Affymetrix rice microarray (as at 1^st August 2009), all experiments containing CEL files were downloaded from the Gene Expression Omnibus within the National Centre for Biotechnology Information database or from the MIAME ArrayExpress database http://www.ebi.ac.uk/arrayexpress/. The GSE or EXP numbers for the respective rice studies are shown in Table 1. There was a total of 373 microarrays for which there was either MAS5.0 data available, thus all of these were used for present/absent determination in defining the list of 7,922 probesets expressed in all microarrays. However of these 373 microarrays, 7 had no biological replicates and 35 did not have available CEL files, thus the remaining 331 microarrays were used to carry out further normalisation (GC-RMA) and calculation of mean, standard deviation and coefficient of variance (CV). This allowed analysis of 117 tissues/conditions, with a minimum of 2 biological replicates. The 117 included 41 organ/developmental tissues, 65 samples within abiotic and biotic stress experiments and 11 samples within hormone treatment experiments.

All raw intensity CEL files were imported into Avadis 4.3 (Strand Genomics) and the standard MAS5.0 normalisation was first carried out in order to determine present/absent/marginal calls for each probeset. For all 331 microarrays with available CEL files (and carried out for biological replicates), GC-RMA normalisation was carried out. The mean expression, SD and CV (=SD/mean) was then calculated for each of the 7,922 probesets across the developmental set, stress set and entire dataset (which included the hormone experiments). On the basis of CV cut-offs, the list of 151 probesets was generated (Figure 1). The averaged log GC-RMA normalised values for these 151 probesets, across the developmental tissues, stress and hormone treatment experiments were hierarchically clustered using average linkage on Euclidean distance. The clustering analysis and heatmap generation was carried out using Partek Genomics Suite, version 6.3 (Partek). For the Agilent microarray comparison, data was retrieved under the accession GSE8518 from the Gene Expression Omnibus within the National Centre for Biotechnology Information database.

Analysis of orthologues

The InParanoid: Eukaryotic Orthologue Groups database (version 7.0) was used to analyse all orthologues between rice and Arabidopsis [23]. The orthologous group file was downloaded for the whole-genome comparison of rice versus Arabidopsis. This produced information for orthologues identified by TIGR identifiers for rice and AGIs for Arabidopsis.

Stress treatments, tissue collection and RNA isolation

In order to analyse the expression of all the genes in Table 2, a selection of tissues were collected across different developmental stages/tissues and under different stress conditions in wild type rice, cv. Amaroo. In order to analyse different developmental tissues; embryos were extracted from dry seed, seeds imbibed for 24 h with oxygen gas (24 h A), seeds imbibed for 24 h in the absence of oxygen gas i.e. in the presence of nitrogen gas (24 h N), seeds imbibed for 24 h under nitrogen gas and switched to oxygen gas for 3 h (27 NA), leaf and root tissues from 2-week old seedlings and 3 month old plants. Furthermore, to examine the effects of abiotic stress, 2-week old seedlings were transferred to 4°C and 42°C for cold and heat treatment respectively over a 9 h time course, whilst the controls remained at a constant temperature of 30°C. All 15 tissue samples were analysed using three biological replicates, the RNA was isolated using the Qiagen RNeasy Plant RNA isolation kit and DNase treated using both the Qiagen on-column DNase digestion as well as the Ambion Turbo DNase treatment exactly as carried out in Howell et al., 2009 [27].

QRT-PCR analysis

Details of the primer sequences and amplicon lengths for each of the genes are shown in Additional file 1, Table S1. The transcript abundance for each gene was analysed using the SYBR green I master (Roche, Sydney) with the Roche LC480. Each sample was analysed in biological triplicate, using individual plants and treatments to test for reproducibility. Following RNA isolation each of the samples was quantitated using a Nanodrop spectrophotometer. This provided the following information for each sample: concentration (ng/μl), the absorbance (A) in nm at 230, 260 and 280, the A₂₃₀/A₂₆₀ and A₂₆₀/A₂₈₀ ratios. Using this information the RNA yield and purity was calculated to ensure that they all had no significant impurities between samples that may affect reverse transcription and/or amplification during QRT-PCR. 1 μg of total RNA was reverse transcribed using the Bio-Rad^® (Sydney) iScript reverse transcription kit, according manufacturer's instructions. In parallel for each sample, another 1 μg of RNA was used for the same reverse transcription reaction, with the exception of the addition of the reverse transcriptase enzyme (no RT samples). Following this, the Qiagen^® PCR purification kit was used according to manufacturer's instructions on all samples (RT and "no RT" samples). This purified cDNA was diluted 1 in 10 with nuclease-free water and 1 μl was used for QRT-PCR analysis. For the no RT samples, no dilution was carried out and 1 μl was used in the same manner as the diluted cDNA for QRT-PCR analysis, this enables the detection of any genomic DNA contamination.