Validation of reference genes as internal control for studying viral infections in cereals by quantitative real-time RT-PCR

ANOVA results

GeNorm Results

GeNorm v3.4 software was used to analyse the expression stability of tested genes in various samples, and to rank them accordingly. The GeNorm is a statistical algorithm which determines the gene stability measure (M) of all the investigated genes, based on the geometric averaging of multiple reference genes and mean pairwise variation of a gene from all other reference genes in a given set of samples [12]. It relies on the principle that the expression ratio of two ideal reference genes is identical in all the samples, regardless of experimental condition and cell-type. Genes with the lowest M-values have the most stable expression. We analysed our data in two sets, one with all samples combined and the second according to plant species. When all the samples were combined, the average M-value of GAPDH and TUBB was lowest, and that of TUBA was highest (Table 3, Fig. 2A). The results remained very similar, when the M-value was measured for species series, with least value for GAPDH and TUBB (Table 3, Fig. 2A). According to the GeNorm tool, the TUBA was by far the least reliable reference gene when all samples were combined as one set, whereas for barley and oats the worst reference gene was EIF4A, and for wheat it was TUBA.

gene	GeNorm stability value (M)					NormFinder stability value
	wheat	barley	oat	all	wheat	barley	oat	all
GAPDH	1,09	1,06	1,08	2,0	0,088	0,112	0,035	0,074
TUBB	1,14	1,27	1,12	2,0	0,041	0,085	0,108	0,098
TUBA	2,13	1,54	1,59	4,43	0,785	0,201	0,205	0,289
18S rRNA	1,48	1,27	1,14	2,15	0,076	0,158	0,125	0,248
28S rRNA	1,82	1,71	1,90	2,78	0,378	0,354	0,258	0,487
EIF4a	1,48	2,81	2,59	2,82	0,352	0,489	0,301	0,499
EEF1A	2,08	2,48	2,49	2,87	0,181	0,187	0,355	0,241

The greater expression of stability is indicated by a lower stability value (M). GeNorm stability is based on an estimate of the pairwise variation (M). For NormFinder analysis, samples were grouped into healthy, abiotic stress, and virus-infected groups. The stability is calculated from the intra- and inter-group variation and the best combination of genes is also given.

Theoretically, the high expression stability of a gene indicates that the use of a single reference gene is appropriate. However, for many studies, no single gene may be adequate and may require normalisation with two or more stable reference genes. Therefore, pairwise variations were calculated using GeNorm for each data set to determine the optimal number of reference genes for normalisation. For this, first the normalisation factors (NF) were calculated for the most stable reference genes (with lowest M-value) and then for other genes by stepwise inclusion of the reference gene that remains most stable. Subsequently, pairwise variations of NF_n and NF_n+1 were calculated, reflecting the effect of including additional (n + 1) gene.

For wheat and oats, four genes with the most stable expression were optimal for reliable normalisation with a pairwise variation cut-off value of 0.15, whereas for barley three genes were necessary (Fig. 2B). Below this cut-off value there is no need for inclusion of an additional reference gene. It was apparent that the inclusion of a fourth/fifth gene had no significant effect on pairwise variation. However, an attempt to calculate the optimal number of reference genes for all samples combined failed; even if all the reference genes were taken into account, the pairwise variation exceeded the cut-off value of 0.15 (Fig. 2B). The quantification might not be accurately normalised even if all the tested genes were used in the assay. This was additional evidence of the infeasibility of comparing BYDV between different cereal species by relative quantification, and also further confirmed the suspicion raised by analysis of raw C_q values that the expression of the tested genes differed among species.

NormFinder results

NormFinder is another approach to assess the stability of expression of reference genes [14]. In the present study, NormFinder indicated that the genes with the most stable levels of transcript abundance were GAPDH and TUBB for all the samples taken as one set and for barley and oat samples, and TUBB and 18S rRNA for wheat samples (Table 3, Fig. 2C). These were almost the same genes as those identified by GeNorm, with the exception that 18S rRNA performed better than GAPDH for wheat samples in NormFinder. The least stable genes, according to NormFinder, were EIF4A for all the samples analysed as one set and also for barley, TUBA for wheat samples, and EF1A for oat samples. Therefore, there was a very strong correlation between the results obtained from GeNorm and NormFinder, despite the fact that the methods of calculation are fundamentally different. NormFinder ranks reference genes according to the least estimated intra- and inter-group variation, which is more effective to control the influence of co-regulation of reference genes. NormFinder can account for heterogeneity in the tested samples, such as different treatment groups, and so distinguishes between stability and bias.

The wheat samples differed more from barley and oat samples than these two species differed from each other. This was especially the case for TUBA which performed well for barley and oats but very poorly for wheat. However, there were significant differences in expression of most of the genes, even between oats and barley, and therefore no intra-species comparison by real-time relative RT-qPCR can be recommended.

Choice of reference genes

The results of NormFinder and GeNorm analyses corresponded well. The two tools have not always agreed on the particular order of the individual genes' expressions' stability values; however, the final choice of the best reference genes was almost uniform. The calculation of V-values by GeNorm for the proposed genes (Fig. 2B) is useful for deciding the optimal number to be used in an expression study; pairwise variation between samples is reduced by the inclusion of additional reference genes and therefore indicates the number of genes required to achieve an arbitrarily selected threshold of reference gene stability; a recommended cut-off value is 0.15. If we took into account the results of the pairwise variations, then four genes were necessary for wheat and oats, and three genes for barley. For barley, the three genes with the most stable expression across the samples as determined by GeNorm and NormFinder were GAPDH, 18S rRNA and TUBB. For wheat and oat, the combination of GAPDH, 18S rRNA, TUBB and TUBA (oat) or EIF4A (wheat) should provide safe normalisation. For oats and barley, the poorly performing genes were 28S rRNA, EIF4A and EF1A; these could not be recommended as reference genes for these species. The use of TUBA as a reference gene is not recommended as appropriate for wheat tissues, as with an M-value of 2.13 it appears to be regulated.

Virus titre comparison results

The analyses described above were applied to the quantification of BYDV in the oat samples. Oats were chosen because the differences in symptoms manifest between the susceptible and the tolerant variety were the greatest (Fig. 3). The virus titre was calculated under three different settings. First, the absolute quantification with cloned standards and no normalisation gene was carried out, as it is still a common practise of real-time quantification of viruses among plant virologists. Secondly, 28S rRNA, EIF4A and EF1A were chosen as reference genes that performed very poorly in both the GeNorm and NormFinder analyses; and finally, GAPDH, TUBB, TUBA and 18S rRNA were chosen as reference genes that performed best in both analyses. The results were analysed by one- and two-way ANOVA (Fig. 4).

There were no significant differences between the viral titre in cv. Saul (susceptible) and cv. Atego (tolerant) according to relative quantification, independently of the reference genes used in the assays. However, the trends of the two reference genes groups were slightly different (Fig. 4B). The GAPDH, TUBB, TUBA and 18S rRNA assays resulted in almost identical data from the two groups; however, the 28S rRNA, EIF4A and EF1A results might lead to the assumption that the tolerant cv. Atego had a lower virus load than the susceptible variety. Furthermore, the BYDV titre was determined to be statistically higher in the tolerant compared to the susceptible variety, by absolute quantification. This seems quite unlikely, and contrary to the plant symptoms, indicating that the absolute quantification with no normalisation is inadequate to quantify viruses in plants. These results support the necessity of the correct choice of reference genes for valid experimental data.

A similar study was conducted by Balaji et al. [20], in which the titre of BYDV-PAV was determined in susceptible and resistant wheat lines, using 18S rRNA as a reference gene. In their study, the titre of the virus was measured 1, 2, 4, 6, 8, 10, 12, and 14 d post-infection (dpi); the greatest differences were in the early stages of infection (4 dpi), when the titre varied as much as tenfold, the most susceptible wheat line being that with higher titre. However, later in infection (≥ 12 dpi), the titre in the susceptible line decreased to almost the same as in the resistant line. As the oat samples in this study were collected at 30 dpi, when the symptom differences appeared, the results of the relative quantification seem to be accordant. Nevertheless, for the comparison of such subtle differences, missing or inappropriately chosen normalisation could easily lead to erroneous interpretation of data.

Plant, virus and vector material

Laboratory isolates of BYDV-PAV and BYDV-PAS maintained on barley seedlings and transmitted by aphids Rhopalosiphum padi and Sitobion avenae were used for the study. For each of wheat (Triticum aestivum L.), barley (Hordeum vulgare L.) and oats (Avena sativa L.), 60 plants were sown into individual pots. For each species, 30 plants were of the BYDV-susceptible and 30 of the BYDV-tolerant variety: wheat cvs. Anza (tolerant) and Vlada (susceptible); oat cvs. Atego (tolerant) and Saul (susceptible); and barley cvs. Wysor (tolerant) and Finesse (susceptible). They were maintained in greenhouses with controlled conditions of 16-h-light period and 20°C. When the plants were at the growth stage where two leaves were unfolded, 20 plants per species (10 per variety) were separated and infected with BYDV-PAS strain by vector transmission. For each plant, approximately five viruliferous R. padi aphids were transferred onto and allowed to feed on the plant for 2 d. After that, the plants were treated with insecticide and kept in the greenhouses for one more week. On the last day of the week period, 20 additional plants of each species were separated and abiotic stress conditions were created for them. The plants were moved into a chamber with 4°C constant temperature and darkness and were kept there for 24 h. Then, all 180 plants were harvested and the whole aboveground biomass of individual plants was used for further analysis.

For the case study, oat plants from a field experiment testing the resistance of different varieties to BYDV were used: cvs. Atego and Saul. Atego did not manifest many symptoms, while Saul symptomatically reacted as very sensitive to the infection. The leaf samples from 20 randomly chosen plants were collected one month after the inoculation with BYDV-PAV and BYDV-PAS. For each plant, the most symptomatic and the least symptomatic leaves were chosen.

RNA purification and cDNA synthesis

The whole plants were ground with a mortar and a pestle in liquid nitrogen, and 100 mg of the material was then used for RNA isolation. The RNA was isolated with a Spectrum™ Plant Total RNA Kit (Sigma Aldrich) according to the manufacturer's instructions. The concentration and purity of the isolated RNA was then spectrophotometrically measured. Of total RNA, 5 μg was diluted to a total volume of 25 μL and digested with DNase I (DNA-free™; Ambion) according to the manufacturer's instructions. cDNA was produced using the M-MLV Reverse Transcriptase (Promega) according to the manufacturer's recommendations for oligo(dT)20 or random hexamers primed cDNA-synthesis. cDNA synthesis was performed using 1 μg of RNA, at 40°C. Finally, cDNA was diluted 1:5 before use in qPCR.

Standards RNA preparation

A specific BYDV nucleotide sequence (294 bp) amplified by RT-PCR was inserted into the vector pGem-T (Promega) and cloned into E. coli JM-109. The plasmid was linearised at the RsaI site and used as a target in an in vitro transcription reaction performed with Megascript T7 kit (Ambion) followed by DNase I (Ambion). The amount of RNA was quantified by spectrophotometry (Nanophotometer, Implen). The μg of single stranded RNA was converted to ρmol using the average molecular weight of a ribonucleotide (340 Da) and the number of bases of the transcript (Nb). The following mathematical formula was applied: ρmol of ssRNA = μg (of ssRNA) × (106 ρg/1 μg) × (1 ρmol/340 ρg) × (1/Nb). Avogadro's constant was used to estimate the number of transcripts (6.023 × 1023 molecules/mol). Thereafter, tenfold serial dilutions of the transcripts were prepared.

Real-Time PCR conditions

Results analysis

The relative gene expression ratios were calculated by a mathematical model, which includes an efficiency correction for real-time PCR efficiency of the individual transcripts [8]. The amplification efficiency was established for each of the targets from serial dilutions of cereal leaves within range 0.80-1.0. The absolute viral gene's quantification values were transformed into relative values by simple proportion. Two-way ANOVA was performed on data, followed by the Bonferroni post-tests, where appropriate. Data in graphs or tables are presented as means with their standard errors of the means. P > 0.05 was considered as non-significant, P < 0.01 very significant and P < 0.001 extremely significant. Reference gene selective analysis was performed with the NormFinder [14] and GeNorm [12] tools.