Testing reference genes for transcript profiling in Uromyces appendiculatus during urediospore infection of common bean

Uromyces appendiculatus is a major pathogen on common bean. Like other rust fungi, it uses effectors to influence its host plant. Effectors are assumed to possess characteristic expression profiles, reflecting their activity during the infection process. In order to determine expression profiles using RT-qPCR, stably expressed reference genes are necessary for normalization. These reference genes need to be tested. Using samples representing seven different developmental stages of the urediospore-based infection process we employed RT-qPCR to measure the expression of 14 candidate reference genes and determined the most suitable ones based on the range of Cq values and comparative calculations using the geNorm and NormFinder algorithms. Among the tested genes RPS14 had the smallest Cq range, followed by Elf1a and Elf3; geNorm rated Tub and UbcE2 best with CytB as a third and NormFinder found UbcE2, Tub and Elf3 as best reference genes. Combining these findings using equal weight for the rankings UbcE2, Elf3 and Tub can be considered the best reference genes. A combination of either two reference genes, UbcE2 and Tub or three reference genes, UbcE2, Tub, and Elf3 is recommended for normalization. However, differences between most genes were relatively small, so all tested genes can be considered suitable for normalization with the exception of RPS9, SDH, Ubc and PDK.


Introduction
Uromyces appendiculatus is a major pathogen on common bean, Phaseolus vulgaris [1]. Common bean is cultivated almost globally both in subtropical and temperate regions. Both pods of P. vulgaris and dry seeds are consumed. With a high protein content and also rich in minerals like iron and zinc common bean is a highly important component of the diet of many people, especially in developing countries where meat may not be easily available [2].
U. appendiculatus is an autoecious rust fungus forming all phases of its infection cycle on one host. This means that in spring basidiospores are formed from teliospores on bean leaf litter, these can infect young bean plants to form pycnia and aecia and aeciospores and later urediospores will further infect the plants. Especially urediospores are produced in vast amounts and can lead to epidemic spread of the disease. Later in the year teliospores are formed to complete the sexual cycle and for overwintering. The urediospore stage of the infection cycle is of highest interest since this is the stage that is most frequently encountered in the field and leads to the highest yield losses. Also subtropical races of U. appendiculatus exist that have lost the ability to form teliospores [1].
U. appendiculatus can be effectively controlled with fungicides and also cultural practices like intercropping or crop rotation together with good sanitation in the field can reduce disease incidence. Since common bean is mostly produced by small scale farmers the efficiency of fungicide application is very limited. Therefore, breeding for resistance is highly important. Several resistance genes have been identified and both cultivar mixes and multi-lines are used for rust control. However, due to the high genetic variability of the pathogen the efficiency of these efforts remains limited [1]. Therefore, broad spectrum and stable resistance is still a highly important objective. Here, studies on rust biology, genetics and physiology can contribute to identify the optimal strategy to establish durable resistance. Studying effectors of rust fungi and their targets is a novel strategy to identify novel sources for resistance. This strategy is also pursued for U. appendiculatus [3][4][5][6]. Part of the characterization of (candidate) effectors is to determine their expression profile. The gold standard for expression profiling is RT-qPCR. For normalization in RT-qPCR internal reference genes are needed. Optimally reference genes should be tested for their stability under the studied conditions [7][8][9]. Optimal reference genes can be identified by comparing between a number of candidate reference genes, for example housekeeping genes. Algorithms that can identify the most stable genes from these comparisons have been published [10,11]. While reference genes for P. vulgaris have been validated, including biotic stress conditions [12], no systematic test for the optimal reference genes has been performed for U. appendiculatus, so far. This might also explain why Cooper et al. [13] used a plant gene as reference for measuring gene expression and fungal biomass in their gene silencing study. Here, a study to identify suitable reference genes for normalization in U. appendiculatus during the phases of infection is presented.
The urediospore infection can be structured into three developmental phases: the penetration phase, the parasitic phase and a sporulation phase [14]. For penetration the spore first attaches firmly to the leaf. Then germination occurs, leading to the formation of germ tubes (gt) and later appressoria (ap) that are positioned over stomata. A penetration hypha and the so called substomatal vesicle are formed, followed by an infection hypha and the haustorial mother cell. The haustorial mother cell is formed upon contact with a plant mesophyll cell and is the basis for plant cell penetration and, as the name implies, for haustorium formation. Up to the haustorial mother cell these infection structures can be generated in vitro on polyethylene sheets that were treated to produce grooves resembling the edges of stomata. Formation of haustoria marks the beginning of the parasitic phase that lasts for several (five to nine) days before sporulation starts. During the parasitic phase the fungus spreads through the plant tissue forming more haustoria. Then hyphae form aggregates that differentiate into uredia producing urediospores. Sporulation can go on for a long time provided the host plant stays fit. However, leaves that are strongly infected-which is always the case upon artificial inoculation -become senescent after a few days of sporulation. For this study samples were evaluated that represent the three phases of the urediospore infection.

Plant material
Common bean (P. vulgaris) cultivar Primel (Vilmorin, La Ménitré, France) was grown in the greenhouse with a 16 h light / 8 h dark regime at 22˚C / 18˚C. The plants were inoculated when they were 21 d old; at this stage they were forming the second trifoliate leaf.

Fungal isolate, plant inoculation, and in vitro production of early infection structures
Urediospores of U. appendiculatus (SWBR1, laboratory collection, Phytopathology, University of Hohenheim, Stuttgart, Germany) were suspended in water (0.05%) together with 0.08% milk powder and 0.01% Tween22 (w/v) and stirred for 30 min at room temperature in the dark. This suspension was sprayed onto P. vulgaris plants using a chromatography vaporizer. The plants were then kept at 95% humidity in the dark at 15˚C for 24 h and then put back to the greenhouse with 16 h light / 8 h dark at 15˚C / 13˚C.
To obtain germ tubes (gt) and appressoria (ap), urediospores were finely dispersed onto polyethylene (PE) sheets using fine gauze material for sifting. For producing ap, the PE sheets were scratched using a brass brush beforehand. After dispersing the spores, the PE sheets were sprayed with water using a chromatography vaporizer and then kept at 20˚C and 95% humidity in the dark. For gt generation the structures were harvested after 4 h, for ap after 9 h. Formation of both gt and ap was checked by microscopy. For later RNA preparation the structures were scraped from the PE sheets, dried by vacuum filtration and stored in 2 ml tubes at -70˚C after freezing in liquid nitrogen.

RNA preparation and cDNA synthesis
Harvesting of structures grown in vitro is described above. Tissue from soybean leaves either uninfected or infected at different stages was gained by cutting out pieces with a cork borer. A cork borer with 12 mm diameter was used; every sample was composed of four pieces cut out from different leaves giving an amount of roughly 100 mg. The samples were collected in plastic tubes, frozen in liquid nitrogen and kept at -70˚C until RNA preparation.
While still frozen, samples were homogenized with a FastPrep 1 -24 MP (Biomedicals GmbH, Eschwege, Germany) set to a shaking speed of 4.5 m/s. For this, two 2 mm steel beads were added to the sample and homogenization was run three times for 20 s; in between homogenization rounds the samples were cooled again in liquid nitrogen. The Agilent Plant RNA Isolation Mini Kit (Agilent Technologies, Santa Clara, CA, USA) was used for RNA preparation. The preparation was done according to the manufacturer's instructions. Because the mini prefiltration columns tended to clog otherwise, an additional centrifugation step to remove cell debris was added after adding the extraction buffer. RNA quality and concentration were checked photometrically by measuring OD 260 and the OD 260 /OD 280 ratio. To test for RNA integrity, also DNA electrophoresis was performed with 1% agarose gels under denaturing conditions. The isolated RNA samples were stored at -70˚C.
To remove residual DNA RNA was treated with DNaseI (Thermo Fischer Scientific, Waltham, MA). Following the recommendations of the manufacturer, 1 μg of total RNA was treated in a 10 μl reaction. cDNA synthesis was performed with the RevertAid RT Reverse Transcription Kit (Thermo Fischer Scientific), also following the manual. To the 11 μl RNA resulting from DNAseI treatment after addition of EDTA the reverse transcriptase and corresponding buffers were added together with random hexamer primers to give a 20 μl reaction mix that was incubated at 25˚C for 10 min, 42˚C for 60 min, and heat inactivated for 10 min at 70˚C. For production of larger amounts of cDNA, reactions were scaled up by proportionally increasing all amounts and volumes. Storage of the cDNA was at -20˚C.

Real time PCR parameters
All Real Time PCR reactions were run on a CFX96™ Real-Time-PCR System (Bio-Rad Laboratories, Hercules, CA, USA). FrameStar 1 96 Well Skirted PCR Plates with white wells and black frame (4ti-0961, 4titude, Brooks Automation, Chelmsford, MA, USA) were used. Plate design, data logging, and basic data analysis including the determination of Cq values as well as the calculation of primer efficiencies from dilution series was performed with the Bio-Rad CFX Manager 2.1 (Bio-Rad Laboratories). For Cq determination, thresholds were first used as determined automatically by the CFX Manager, for full comparability 3000 rfu as a rounded estimated average of automatic thresholds was used.
20 μl Real Time PCR reactions were performed using the SensiFAST™ No-ROX mix (Bioline Reagents Ltd., London, UK) with primers at a concentration of 0.4 μM and 2 μl cDNA template (directly from reactions described above) or alternatively water in the NTC or RNA (1 μg / 20 μl in the same buffers as the cDNA) in the NRT controls respectively.
A two step PCR protocol was used with an initial denaturation step of 95˚C for 5 min, followed by 40 cycles with 5 s at 95˚C and 15 s at 60˚C for annealing and elongation. A melt curve analysis was also run starting at 65˚C to 95˚C in steps of 0.5˚C with a duration of 5 s and measuring of fluorescence after every step. For correct readout, the instrument was set to the BR white plate setting and the SYBR/FAM only scan mode was used.

Data analysis
For determination of the most stable candidate reference genes the GenEx software package (GenEx 6.0.1.612, MultiD Analysis AB, Göteborg, Sweden) was used. Cq values exported from the CFX Manager were first averaged for the replicates, and then corrected for efficiency. Both the NormFinder and geNorm algorithms were run as implemented in GenEx. Graphs in this publication showing Cq values are showing Cq values after efficiency correction. The GenEx tool Data editor was also used for normalization and calculation of relative expression data. All graphs were created in MS Excel (Microsoft Corporation, Redmont, WA, USA).

Primers
Primers (Table 1) were designed using the webtool Primer3 [15,16] with the following settings: Primer size range 18-23 nt with an optimum of 20 nt, optimal melting temperature (AT+GC) 60˚C with a range from 58-62˚C, GC content 50% with a range from 30-70%, product size 75-150 bp. Other settings were kept at default. The primers were cross-checked with GeneRunner (version 6.5.52 Beta, generunner.net) and occasionally shifted by a few bases to have a G or C at the 3' end. Primers were ordered in HPCL quality from biomers.net GmbH (Ulm, Germany). Primers were tested for efficiency using cDNA from RNA from spores in four concentrations, undiluted, 1:10, 1:100, 1:1000. NTC and NRT controls were also included and the primers were also tested against cDNA from RNA prepared from non-infected plant material. When primers were found unsuitable due to formation of artifacts (in NTC), alternative primers proposed by Primer3 were used or, as a third alternative, primers were designed directly using GeneRunner.

Infection stages sampled
Samples were taken to represent the major phases of urediospore infection, the penetration phase, the parasitic phase with formation of haustoria and the sporulation phase. During the penetration phase the amount of fungal RNA in the infected leaf is too low for reproducible amplification with RT-qPCR. Therefore, we prepared RNA just from urediospores (sp), from germinated spores or germ tubes (gt) and from infection structures differentiated on structured PE sheets for 9 h when more than 50% of all structures showed appressoria (ap). Incubating for a longer time might have led to a higher percentage of structures reaching the appressorial stage. However, our earlier study on reference genes for Phakopsora pachyrhizi [17] showed, that waiting too long also leads to degradation of RNA. Collectively sp, gt, and ap are also named in vitro stages or samples as opposed to in planta samples.
Infected plants were sampled at 3 days post inoculation (dpi), 5 dpi, 7 dpi, and 14 dpi. 3 dpi and 5 dpi represent the parasitic phase, 7 dpi marks the beginning of the sporulation phase and 14 dpi the later sporulation phase with leaf senescence. 3 dpi was considered the earliest time point when amplification with RT-qPCR from infected plant samples was possible.

Candidate reference genes
In our study on reference genes for P. pachyrhizi [17] we tested 15 candidate reference genes carefully chosen from other studies testing or using reference genes for rust fungi [18][19][20][21]. Since this study was successful this previous choice was regarded as sensible, useful, and good to be used again. So, for this present study on U. appendiculatus, it was decided to test the respective homologous genes. The homologs were identified by BLASTing the predicted protein sequences from P. pachyrhizi against the predicted protein sequences of the U. appendiculatus transcriptome [3]. This approach allowed the identification of the genes and sequences listed in Table 2. For RPS11 (40S ribosomal protein S11) no adequate homologous sequence could be identified.

Expression of the candidate reference genes in different samples
For the selected genes primers were designed and as recommended by Pfaffl [22] amplification efficiencies were determined (Table 1). For GAPDH, PDK, RPS14, and Tub, the primers that were chosen first produced artifacts in the NTC. So for these genes alternative primers were chosen; for GAPDH and RPS14 this was a modified version of alternative primers suggested by Primer3; for PDK both primers and for Tub the new reverse primer were freely chosen using GeneRunner. RT-qPCR was run for all candidate reference genes and all samples. For an overall crude impression on expression of the candidate reference genes the efficiency corrected means of the Cq values of the technical replicates were plotted for all samples (Fig 1). In the plot that shows the Cq values for all samples separately (Fig 1a), it initially seems that the expression of all candidate reference genes is not stable between the samples. This apparent instability is not contradictory to the use of these potential reference genes, however. The differences in the Cq values can be explained by general changes in the percentage of fungal mRNA in the total RNA prepared from the different samples: Cq values for the gt samples are generally lower than those for the sp samples. On average the difference between the Cq values for the two stages is three, meaning that gene expression in germinated spores is eight times higher than in the spores. Since the metabolism is resting in the spores, this difference is to be expected. By contrast, the differences between the gt samples and the ap samples are very small, as should be expected for stably expressed genes. This is, however, in contrast to our results for P. pachyrhizi where Cq values for the appressorial samples were much higher than for the germinated spore samples [17]. As already mentioned, we attributed these higher Cq values to  degradation of mRNA in those specimen that were already dead for a while at the time they were collected. The fact that Cq values for the chosen housekeeping genes are now stable between the gt and ap samples is clearly a better and correct result and it can be concluded that the decision to shorten the incubation time for appressoria generation was the right choice.
There is a second shift in Cq values from an average of 19.2 in the ap samples to an average of 27 in the 3 dpi samples. The reason for this shift is that in the in planta samples the fungal RNA is only a fraction of the total RNA-fungal RNA is "diluted" by the RNA of the host plant. Towards the later stages the average Cq values for the samples are again decreasing. This is explained by the increasing fungal biomass in the leaf over time as the fungus is growing. Since the graphs for the different candidate reference genes are mostly parallel, it appears that these genes nicely represent overall mRNA abundance in the different samples. Together with all other samples also the samples taken from non-infected plants were tested. Almost all primer pairs gave some amplification with this sample. The Cq values obtained from this control were all >30 (S1 File) and the melting curves showed peaks at similar temperatures as the melting curves obtained during testing the primer efficiencies. This might indicate that the primers also have some affinity to plant RNA, presumably the P. vulgaris homologs of the candidate reference genes. The possibility that the non-infected plant RNA sample actually contained some U. appendiculatus contamination might be another reason for these results. However, since the RNA was also tested with other U. appendiculatus primers with no amplification, this eventuality can be ruled out. Subsequently these Cq values were compared with those for the other samples. Almost all were much higher, so it can be supposed that the products from the amplification of plant RNA do not interfere with the analysis of the amplification of the fungal genes. Two exceptions were noted: RPS9 and SDH that both have relatively high Cq values came close to the Cq values of the control in the 3 dpi sample. Therefore, RPS9 and SDH were excluded from all further analyses including in planta samples.

Identification of the best reference genes
To find the most suitable reference genes the range of Cq values over all samples was considered for all genes. To provide graphical representations of these ranges the Cq values are represented in box plots in Fig 1b to 1d. These box plots represent the Cq ranges for the candidate reference genes either for just the in vitro samples (Fig 1b), the in planta samples (Fig 1c) or all samples together (Fig 1d). The Cq range was narrowest for Act when just the in vitro samples were considered (Fig 1b); when the in planta samples or all samples were considered RPS14 had the smallest range of Cq values (Fig 1c and 1d) making it the gene with the least difference in mRNA abundance over all tested samples. Ubc, PDK, RPS9, SDH, and Cytb, on the other hand, had a very wide range of Cq values. These Cq ranges as a first aspect for deciding, which is the best reference gene, were used for a ranking (Table 3). This ranking was included into an overall ranking to represent the aspect of the Cq range.
Next, the stability of expression was tested with the algorithms geNorm and NormFinder. Both algorithms calculate differences between the Cq values of the different genes and compare these differences between the different samples provided. When these differences are mostly constant between the different samples the gene expression is deemed stable, the more variations in the differences are noted, the less stable is the expression of the gene. While geNorm repeatedly compares each gene against an average of all other genes and removes the gene with the worst performance (the highest M value) from the group after every round of calculations, the NormFinder algorithm is model based [11] and as implemented in GenEx provides a standard deviation (SD) for every gene. The functioning of geNorm determines that there is always a pair of genes performing best and the authors state that at least three genes should be used for normalization [10]. NormFinder, on the other hand, can rank all genes by SD and also provides an accumulated standard deviation for a group of genes that indicates what should be the optimal number of genes used for normalization. Like for the Cq ranges these analyses were performed both separately for the in vitro samples and the in planta samples and for all samples together. Fig 2 shows the results for the in vitro samples. The geNorm algorithm predicts Elf3 and SDH as the best genes with Tub third (Fig 2a), while NormFinder finds SDH and Tub best with PDH third (Fig 2b). Altogether all values are quite low so most of the genes seem to be suitable although there is a group of seven genes including SDH, Tub, PDH, PDK, Elf3, RPS14, GAPDH, and UbcE2 that are markedly better than the other six. These two groups are the same for both algorithms indicating an overall similarity of the results from both algorithms. From accumulated SD it appears that two reference genes should be optimal (Fig 2c) but there is so little difference between the values that even one gene could be sufficient. Based on these results, either SDH, Tub, or Elf3 could be used as reference genes for normalization between in vitro samples.
As mentioned already, the genes SDH and RPS9 were excluded from the analysis of the in planta samples. Here geNorm identified Elf1a and Elf3 as the best genes with RPS14 third (Fig  3a). NormFinder calculated Act and PDH as best and UbcE2 third. Two genes seem to be optimal for normalization (Fig 3b and 3c). For these stages the rankings are differing between the two algorithms. Only Ubc and PDK are consistently identified as the least stable genes. Since the values are all relatively low, the different rankings can be interpreted such as there seem to be small differences between the different genes and have similar suitability as reference genes.
An analysis of the genes over all stages (also without SDH and RPS9) resulted in Tub and UbcE2 as best genes predicted by geNorm with Cytb as a third and UbcE2, Tub, and PDH as best with NormFinder (Fig 4a and 4b). For this analysis six reference genes are calculated as optimal (Fig 4). Since there is almost no difference in the accumulated standard deviation between the different combinations except between using only one reference gene or two, however, it seems that two reference genes should be sufficient again.

Reference genes for Uromyces appendiculatus
Optimally reference genes should be stably expressed and this should be reflected in a constant Cq or a narrow range of Cq values, when a constant amount of RNA/cDNA is used as template for PCR. Therefore, reference genes might be selected for the narrowest Cq range. There can, however, be reasons for changes in Cq that do not diminish the suitability of reference genes. Such reasons are the change from dormant spore to metabolically active germinating spore and the difference in measurement of only fungal RNA and a mixture of fungal and plant RNA. Since the algorithms geNorm and NormFinder consider only the differences between Cq values between the genes, these parallel shifts in Cq are irrelevant. To select reference genes that are both stable by the means of geNorm and Normfinder and also do not have a large Cq range, a combined score was used. For this, the rankings by the three parameters were added up and a new ranking was performed based on the sums (Table 3). In this overall ranking SDH was the best gene for in vitro samples, followed by Elf3a and Tub. For in planta samples it was UbcE2, followed by PDH and Elf3a together with Cytb. For all sampled stages UbcE2, Elf3, and Tub seem to be most suitable.
This study provides reference gene combinations for three sets of samples, in vitro stages, in planta stages and all samples together. Of these the all samples set is probably the most relevant, since studies, that aim to identify genes involved in the establishment and maintenance of the biotrophic interaction, i. e. effectors, are searching for genes that are up-regulated in in planta samples and have to compare in vitro and in planta samples. Studies using only the in vitro samples could be performed to specifically search for genes that are highly expressed in the appressorium while searches for genes that are involved in spore formation could be limited to the in planta stages.
The proposed reference genes are relatively similar for all three sets of samples, so the reference gene combinations can be exchanged. An exception is SDH, which should not be used for in planta samples and also UbcE2 may not be the first choice for a study focusing entirely on the in vitro samples. Otherwise a combination of either two reference genes, UbcE2 and Tub or three reference genes, UbcE2, Tub, and Elf3 is recommended for normalization in RT-qPCR for gene expression profiling during the stages of urediospore infection of U. appendiculatus on P. vulgaris.
In our study for reference genes in P. pachyrhizi [17] we also included a sample infected with Bean Pod Mottle Virus (BPMV) to check the stability of the reference genes in a host induced gene silencing (HIGS) setup. Since we do not perform the corresponding experiments for U. appendiculatus, such a sample was not included here. Cooper et al. [13], however, could demonstrate the effectiveness of HIGS against U. appendiculatus. In their study they do not yet report measuring mRNA levels of the silenced gene during silencing. Based on our results for P. pachyrhizi, where no notable changes in Cq were observed for any of the candidate reference genes it can be assumed that the reference gene combinations that are best for either the in planta samples or all samples could also be used for assessing the success of HIGS experiments.

Validation of the reference genes
To validate the reference gene combination UbcE2, Elf3, and Tub, these were tested for transcript profiling of effector candidates. We recently published a study identifying candidate effectors in U. appendiculatus [4]. Part of that study was to establish the expression patterns of these genes. For this analysis we used Act and CytB to normalize. This selection of reference genes was based on earlier results from this present study using the same strategy for selection of the best reference genes as presented above. In this earlier analysis the genes GAPDH, PDK, RPS14, and Tub were not included because at this time suitable primers for these genes had not yet been identified.
For testing the two effector candidates Uaca_9 and Uaca_12 that showed distinctive expression patterns in our earlier study were chosen. Uaca_9 showed almost no changes in expression whereas Uaca_12 showed strongly induced expression in the in planta samples.
In the experiment I now again found almost constant expression for Uaca_9 and strong induction for Uaca_12 (Fig 5). This can be accepted as additional evidence that the combination of reference genes UbcE2, Elf3, and Tub is suitable for normalization over all included samples. On the other hand, this analysis also validates the results of the old experiment. Even though in the final analysis Act and CytB were not among the best reference genes they are still among those that can be regarded as suitable for normalization. The similarity of expression patterns shown in Fig 5 and those in [4] corroborates this assumption. Therefore, it can be stated that suitable reference genes were identified and the expression patterns presented in [4] are valid.

Conclusions
This study shows that Cq values of housekeeping genes vary between different samples representing the development of U. appendiculatus on P. vulgaris. Part of these variations is due to the fact that the U. appendiculatus RNA is either pure in the in vitro samples or mixed with P. vulgaris RNA in the in planta samples. But there are also definite changes in gene expression between the resting urediospore and the germ tube samples. Since the resting urediospore stage is metabolically inactive and differs in gene expression from all other stages it is not advisable to use this stage as a reference sample; the germ tube stage should be used instead as the sample relative to which expression should be measured.
The Cq values are relatively stable between gt and the ap sample. Compared to our results for P. pachyrhizi [17] this is a considerable improvement. It is, therefore, an important consideration to incubate the in vitro structures on PE sheets for 9 h to observe appressoria but not much longer. Table 3 represents all rankings of the candidate reference genes. Using a combined ranking based on Cq ranges, geNorm, and NormFinder a combination of either two reference genes, UbcE2 and Tub or three reference genes, UbcE2, Tub, and Elf3 is recommended for normalization.