Stability analysis of reference genes for RT-qPCR assays involving compatible and incompatible Ralstonia solanacearum-tomato ‘Hawaii 7996’ interactions

Reverse transcription-quantitative PCR (RT-qPCR) is an analytical tool for gene expression quantification. Reference genes are not yet available for gene expression analysis during interactions of Ralstonia solanacearum with ‘Hawaii 7996’ (the most stable source of resistance in tomato). Here, we carried out a multi-algorithm stability analysis of eight candidate reference genes during interactions of ‘Hawaii 7996’ with one incompatible/avirulent and two compatible/virulent (= resistance-breaking) bacterial isolates. Samples were taken at 24- and 96-h post-inoculation (HPI). Analyses were performed using the ∆∆Ct method and expression stability was estimated using BestKeeper, NormFinder, and geNorm algorithms. TIP41 and EF1α (with geNorm), TIP41 and ACT (with NormFinder), and UBI3 and TIP41 (with BestKeeper), were the best combinations for mRNA normalization in incompatible interactions at 24 HPI and 96 HPI. The most stable genes in global compatible and incompatible interactions at 24 HPI and 96 HPI were PDS and TIP41 (with geNorm), TIP41 and ACT (with NormFinder), and UBI3 and PDS/EXP (with BestKeeper). Global analyses on the basis of the three algorithms across 20 R. solanacearum-tomato experimental conditions identified UBI3, TIP41 and ACT as the best choices as reference tomato genes in this important pathosystem.

The Ralstonia species complex comprises a diverse set of strains that induce bacterial wilt (BW) disease in tomato (Solanum lycopersicum L.) and across a large number of botanically unrelated hosts 1,2 . The Ralstonia complex is currently composed of the following species: R. solanacearum (= phylotype II from the Americas); R. pseudosolanacearum (= phylotypes I and III); and three subspecies of R. syzygii (= phylotype IV) [2][3][4][5] . The present classification is supported by both genomic and proteomic data 6,7 . Crop losses due to diseases induced by members of this complex are estimated to be around one billion dollars per year 8

.
Ralstonia species are not able to invade intact root surface (epidermal) tissues 9 . As a result, for these pathogens to penetrate host tissues, the presence of wounds or natural openings are required, such as partially exfoliated cells of the outer parenchyma layers at the emergence sites of secondary roots 10 . Pathogenesis of virulent Ralstonia strains has been described as a three-stage process 11 , including colonization of the root surface (day 1), infection of the root cortex (days 2-3), and infection of the vascular parenchyma followed by invasion of the xylem (day 3 and beyond). Ralstonia species synthesize virulence factors that act in facilitating infection, root xylem penetration/establishment, systemic colonization, and induction of symptoms 12 .
The inbred line S. lycopersicum 'Hawaii 7996' is considered the most stable tomato breeding source, displaying high levels of broad-spectrum resistance/tolerance to a wide array of isolates of R. solanacearum and R. www.nature.com/scientificreports/ pseudosolanacearum, as indicated in worldwide multi-location evaluations [13][14][15] . BW resistance in 'Hawaii 7996' is not associated with prevention of root invasion, but with the ability to restrict bacterial colonization of vascular tissues 10 . Recently, screening assays indicated that the expression of 'Hawaii 7996'-mediated resistance/tolerance may vary according to the bacterial strain employed in trials. Challenging 'Hawaii 7996' against a collection of R. solanacearum and R. pseudosolanacearum isolates revealed that, in many circumstances, resistance is better classified as either phylotype-specific or strain-specific [16][17][18] . In Brazil, a subgroup of R. solanacearum isolates was reported to cause severe breakdown of 'Hawaii 7996'-mediated resistance 14,19 . Two of these isolates (namely 'RS 488' and 'CCRMRs223') were selected based on their ability to induce severe BW symptoms in 'Hawaii 7996' . Conversely, a distinct R. solanacearum isolate ('RS 489') was reported as avirulent/incompatible on 'Hawaii 7996' , but highly virulent/compatible on a set of standard susceptible cultivars 14 . From a breeding standpoint, it is of extreme interest to determine which genes are modulated (either triggered or turned off) in 'Hawaii 7996' by these virulent/compatible (= resistance-breaking) R. solanacearum isolates in comparison with the avirulent/incompatible isolates.
Reverse transcription-quantitative PCR (RT-qPCR) is one of the most reliable analytical tools to measure gene expression across a wide range of experimental conditions 20 . Robust analyses, however, require the employment of a suitable set of reference genes for transcript expression normalization 21 . In accordance with the Minimum Information for Publication of Real-Time Quantitative PCR Experiments (MIQE) guidelines 22 , normalization of expression of each target gene in relation to a reference gene is required to adjust data for variations across samples that often occur during cDNA preparation 23 .
For the interaction Ralstonia-tomato, although expression stability of three genes [viz. phosphoglycerate kinase (PGK), alfa-tubulin (TUB), and ACT ] was previously evaluated in stem tissues of inoculated plants 33 , there is so far no consistent set of reference host genes available for analysis of interactions among 'Hawaii 7996' and compatible and incompatible R. solanacearum isolates. In this study, we carried out multi-algorithm stability analysis of eight candidate reference genes for potential employment in host-pathogen interaction analysis in this important pathosystem. The precise determination of relative changes in gene expression profiles in comparative assays provides crucial and applicable information for the discovery of target genes in classical and molecular tomato breeding programs for BW resistance.

Results
Transcriptome RPKM values. Alleles of ACT , APT, TUB2, EF1α, EXP, TIP41, PDS, and UBI3 genes in tomato were selected as candidate reference genes for RT-qPCR validation. RPKM values for these genes are summarized in Supplementary Table S1, which shows expression values with overall low standard deviations across all conditions. Calculation of the Cq of samples. Variation of Cq values across treatments for each gene was verified, with a Box-plot representation summarizing expression ranges and average Cq values for each gene under the two pre-defined conditions: interaction with compatible ( Fig. 1A) and incompatible (Fig. 1B) R. solanacearum isolates. Mean Cq values for the candidate reference genes ranged between 22.50 and 29.00. Results indicated that none of the selected genes were uniformly expressed across all samples. As such, it was necessary to evaluate the candidate genes for normalization using different sets of interaction samples (analyses #01 to #20; Table 1). The set of analyses involving all contrasts among compatible and incompatible interactions (at 24 and 96 HPI) was named as the "Global Interaction"-GI (= analysis #01 in Table 1), whereas the set of analyses involving samples of the avirulent/incompatible isolate (at 24 and 96 HPI) was named as "Global Incompatible Interaction"-GII (= analysis #02 in Table 1). The results of these two analyses are summarized in Table 2. All additional specific contrasts are presented in Supplementary Table S2. Analysis of gene expression stability data of the candidate reference genes. Expression stability and ranking of the candidate reference genes in 'Hawaii 7996' was determined using the algorithms geNorm, NormFinder, and BestKeeper. On the basis of geNorm analyses, all candidate reference genes showed high stability values (M < 0.5) in GI (#01) and GII (#02) ( Table 2, Fig. 2). In GI, PDS, TIP41, and UBI3 genes showed the greatest stability (M = 0.205 for all three genes), with EF1α the least stable of the tested genes (M = 0.49) ( Table 2, Fig. 2A). In GII, TIP41, EF1α, and ACT genes displayed greatest stability (M = 0.032, 0.034 and 0.040, respectively), with APT the least stable (M = 0.456) ( Table 2, Fig. 2A). PDS and TIP41 were ranked as the Best Combination of Two Genes (BCTG) for GI. EF1α and TIP41 were the BCTG for the GII (Table 2). In the specific interactions (analyses #03 to #20; Table 1, Supplementary Table S2) EXP was the most stable reference gene in analyses #07, #09, #12 and #19; ACT in #04, #05, #13 and #14; PDS in #10, #15, #17 and #20; TUB2 in #03 and #11; TIP41 in #06 and #18; UBI3 in #08 and APT in #16.
Paired variation analysis (V) indicated that the combination of V2/3 genes met the normalization parameter requirements, with values of V < 0.15 (baseline) in GI and GII (Fig. 2), as well as in the remaining set of analyses (from #03 to #20, data not shown). V2/3 data indicated that any combination of two genes among the three most stable (identified according to their M values) was suitable for employment as reference gene pairs for expression www.nature.com/scientificreports/ normalization of the target genes. The inclusion of additional reference genes does not contribute significantly to the normalization of gene expression data (Fig. 2B).
Analyses using the algorithm BestKeeper revealed all candidate genes employed in GI and GII displaying standard deviation (± CP) values lower than 1.0 (Table 2), with the exception of EF1α, with standard deviation (± CP) values of 1.974 and 1.445. For GI, greatest expression stability was observed with UBI3, PDS, EXP, and ACT . For GII, the most stable candidate genes were ranked as UBI3, followed by TIP41, PDS and ACT (Table 2). With this algorithm, an overlap of candidate reference genes was observed when compared to ranking positions based on geNorm and NormFinder data, especially for GI. For the specific set of analyses involving 'Hawaii 7996'-R. solanacearum interactions (#03 to #20 in Table 1) UBI3 was the gene with the lowest standard deviation [± CP] across most of the conditions (#04, #06, #07, #08, #09, #10, #11, #12, #13, #14, #15, #17, #18 and #19, Table S2). EF1α displayed a value higher than 1 across all the tested conditions. Interestingly, in the subset #11, none of the genes analyzed displayed adequate levels of stability according to the algorithm BestKeeper (all standard deviations > 1.0; Table S2). Subsets #13, #14 and #16 also showed numerous candidate genes with CP values above 1.0.
Overall analysis across the 20 specific conditions and on the basis of the three algorithms indicated that UBI3, TIP41 and ACT can be considered as the most suitable choices as reference genes for normalization for the interaction R. solanacearum-tomato.

Discussion
A number of previous studies have developed reference genes for normalization of RT-qPCR gene expression data in tomato, with the main genes employed to date for expression normalization comprising TUB4, TUB2, UBQ, GADPH, 18S RNA, ACT , EF1α, and APT1. These genes were evaluated under distinct experimental conditions, which included combinations of tomato genotypes, tissue types, developmental stages, distinct abiotic stresses, as well as mycorrhizal fungus-root interactions and response to fungal, viral, and bacterial pathogens [24][25][26][27][28][29][30][31] . In addition, alternative genes have also been employed as references in RT-qPCR analyses in tomato, including EXP, TIP41 26 and PDS 28 . SGN-U314153 (CAC ) and SGN-U346908 ("Expressed") genes have also been reported as an optimum combination of reference genes to normalize gene expression data during fruit development 32 ,   www.nature.com/scientificreports/ with the 2A catalytic subunit (PP2Acs) and SGN-U321250 (TIP41) the most appropriate genes for postharvest fruits subjected to electricity-induced stress 34 , and PP2Acs and actin (ACT2) identified as the most appropriate for seed tissues 35 . In a previous study of the interaction between tomato and the bacterial pathogen Xanthomonas campestris pv. vesicatoria, out of a total of 11 evaluated genes, the most suitable reference genes were identified as those coding for PHD finger family proteins and the U6 snRNA-associated protein LSm7 36 .
In this study, we evaluated eight candidate reference genes [viz. ACT , APT, TUB2, EF1α, EXP, PDS, TIP41, and UBI3] for potential RT-qPCR analysis of gene expression during interactions among R. solanacearum isolates and S. lycopersicum 'Hawaii 7996' , a genotype considered thus far the best breeding source of resistance to this group of pathogens. In addition, we were able to dissect interaction expression profiles that represented sampling of the whole plant, i.e., a bulk of leaf, root, and stem tissues, at 24 HPI and 96 HPI. Moreover, the reference genes described here were found to be appropriate for investigation of host gene expression during interaction with R. solanacearum isolates with contrasting virulence profiles (compatible versus incompatible against this resistant tomato line). Hence, our experimental conditions were distinct and more comprehensive than the single previous study dealing with the tomato-R. solanacearum pathosystem 33 . In this former study, RNA expression stability of the housekeeping genes PGK, TUB, and ACT in stem tissues in the tomato cultivar 'King Kong 2' (which displays only moderate levels of resistance) was evaluated in conjunction with a silicon priming treatment conducted prior to infection with a highly virulent R. solanacearum strain (race 1, biovar 3).
Our results were consistent, and at least one pair of selected reference genes were appropriate for gene expression analysis in each of the 20 possible sample combinations. Sampling time, as well as pathogen isolate, were considered in reference gene development, with ranking of stably expressed reference genes consistently varying across the treatments. These observations reinforce the importance of validating reference genes for specific sets of experimental conditions, especially in the case of the tomato-R. solanacearum pathosystem. Primer sets  (Table 1) and the Global Incompatible Interaction #02 (see also Table 1)  www.nature.com/scientificreports/ developed here were found to be appropriate for investigation of host gene expression at 24 and 96 HPI, which are crucial time points in terms of R. solanacearum infection and pathogenesis. Previous histological investigation of tomato infection with compatible Ralstonia strains reported initial colonization of root surfaces and intercellular cortical spaces at 24 HPI, followed by cortical infection of vascular parenchyma and xylem from 72 HPI onwards 11 . Normalization of RT-qPCR experiments ensures that results are both statistically significant and biologically meaningful 32 , with inappropriate reference genes potentially resulting in misinterpretation of gene expression data and erroneous conclusions 37 . Furthermore, in specific cases for bacterial pathogens, the most suitable reference genes have been shown to be dependent upon the experimental conditions employed for calibrating RT-qPCR analyses. For this reason, reference genes were validated for specific virulent and avirulent 'Hawaii 7996'-R. solanacearum interactions at different times after inoculation. Results revealed that no single algorithm was appropriate to exclusively determine candidate reference gene expression stability in 'Hawaii 7996' in these interactions. Additionally, as each algorithm for analysis of expression stability is based on a different statistical method, differences in gene stability ranking often occur when analyzing datasets. As such, a standard procedure is to consider methods as complementary, with outcomes equally weighted and combined into a 'consensus' ranking. Here, the most stable genes varied according to the software employed, with the exception of the TIP41 gene. Overall, the housekeeping and 'classical' genes performed as well as the 'nonclassical' genes. A certain degree of overlap in the selection of the two most stable reference genes was observed across the algorithms.
Although expression stability analyses in the same experimental treatment revealed differences when employing algorithms geNorm and NormFinder, the most highly stable pairs of genes, as well as the same single gene, were identified by both algorithms in 80% of the analyses. For example, when considering GI and GII, TIP41 was consistently identified as the most stable gene across both experiments with both algorithms, although the second gene was quite divergent. Identical stability ranking was similarly reported 28 , for two out of ten reference genes examined for normalization of gene expression by RT-qPCR in the tomato-begomovirus interaction using geNorm, NormFinder, and BestKeeper algorithms. In this scenario, the employment of two or more algorithms should be endorsed as a standard methodological procedure in order to obtain more precise assessments of reference gene stability 26,38 . BestKeeper is an algorithm that identifies the most stable reference gene based on pairwise correlation analysis of candidate genes, calculating the geometric mean of the most stable. As a default, reference genes displaying a variation in amount of starting material by a factor of two or more are considered unstable 39 . This equates, in a PCR reaction with an amplification rate of two (100% reaction efficiency), to genes whose Cq values display standard deviation (± CP) values > 1. Given that in biological samples the reaction efficiency is rarely 100% 40 , the default setting of BestKeeper might be inappropriate for in vivo samples, limiting applicability to restricted experimental conditions 41 . Examples of discrepancies in rankings of gene stability data with BestKeeper in relation to the algorithms GeNorm and NormFinder has been documented in recent analyses of plant gene expression 42 . In our analysis, BestKeeper determined CP < 1 values for most genes, which were ranked according to their stability. However, samples extracted at 24 HPI appear to show higher CP values, as observed in groups #3, #5, #7, #11 and #16, with no stable gene determined for group #11 (Table S2). The  Table 1) and the Global Incompatible Interaction #02 (see also www.nature.com/scientificreports/ groups where samples extracted at 24 and 96 HPI were analyzed together also showed a tendency for higher CP values, as observed in groups #13 and #14. This may be related to the initial level of gene expression in the plant biological response, with RPKM values from in silico data showing similar expression levels in these genes to the mock 0 HPI sample, decreasing at 96 HPI (Table S1). Across all specific sets of R. solanacearum-tomato 'Hawaii 7996' interactions, TIP41, EXP, UBI3 and ACT outperformed all other genes in terms of expression stability. Interestingly, in contrast to our study, the ACT gene (actin 4; Solyc04g011500) was previously reported as displaying unstable expression 33 , in particular during the early infection phase. Here, we employed a different ACT gene (= Actin 51; SolycSolyc11g005330.2.1), using a distinct primer pair, and sampling both aerial and root tissues. This ACT allele displayed a stable level of mRNA expression across all specific sets of 'Hawaii 7996'-R. solanacearum interactions, highlighting the importance of developing reference genes for normalization with specific alleles that are suitable for both plant genotype and imposed experimental conditions. In addition, it is recognized that the information from paired gene variation is required to perform accurate normalization assays 28 . Hence, according to our results, any of the combinations of pairs of genes PDS/TIP41, PDS/UBI3, or TIP41/UBI3 may be employed as reference genes for GI. Similarly, TIP41/EF1α, TIP41/ACT , or EF1α /ACT were found to be the most appropriate combinations of reference genes for GII.
We employed RNAseq transcriptome data to double-check the stability of the candidate reference genes for the R. solanacearum-tomato 'Hawaii 7996' pathosystem based on RPKM analysis. We found that the transcriptome per se functioned as a reliable predictor of reference gene stability, with all the reference genes tested via RT-qPCR, including the PDS gene (which showed low RPKM values), showing stable expression across the treatments.
In summary, this study is the first comprehensive investigation into validation of appropriate reference genes for one of the most economically and biologically relevant tomato pathosystems. Stable gene combinations were identified which are appropriate for accurate examination of target gene expression via RT-qPCR. The emergence of novel R. solanacearum isolates that are able to break down 'Hawaii 7996'-derived resistance (which is one of the few effective sources of resistance in tomato breeding programs) is a major threat that will require advanced biotechnological approaches to develop genetic solutions for control of these bacterial variants. In this scenario, well-designed transcriptomic and RT-qPCR assays will play a major role in molecular breeding programs. The employment of these sets of genes with stable expression will have crucial methodological relevance for the identification of key genes of interest for development of tomato cultivars with stable and durable BW resistance.

Methods
Plant material and inoculation assays with contrasting R. solanacearum isolates. Whole seedlings (all aerial organs plus root system) of the tomato inbred line 'Hawaii 7996' were employed in the bioassays. The use of plants in the present study complies with international, national and/or institutional guidelines. All plant material was provided by Embrapa. Three R. solanacearum isolates with contrasting virulence/compatibility patterns in relation to 'Hawaii 7996' were employed in the study. The isolates CCRMRs223 (from tomato, Brazil, sequevar IIA-63) 19 and 'RS 488' (from tomato, Brazil, sequevar IIB-1) are able to induce severe wilt symptoms in 'Hawaii 7996' 14 . Contrasting results were obtained with the isolate 'RS 489' (from tomato, Brazil, sequevar IIA-50). 'RS 489' is known to be avirulent/incompatible on 'Hawaii 7996' 14 , but is highly virulent/compatible against a set of susceptible tomato cultivars/lines (viz. 'Caline IPA-6' , 'LS-391' , and 'hybrid TY-2006') 14,18,19 . Seedlings (30 days after sowing, displaying four fully-expanded leaves) were separately inoculated with the R. solanacearum isolates via root spraying, with bacterial suspensions adjusted to 1 × 10 8 CFU mL −1 . Subsequently, inoculated plants were transplanted to 500 mL pots filled with a sterile substrate mixture. Symptoms were observed in susceptible cultivars 5-7 days after inoculation, demonstrating the adequacy of the inoculation procedure.
Experimental design and biological samples. Experiments were carried out under greenhouse conditions at CNPH (Brasília, DF, Brazil) in a completely randomized design, with three biological replicates of ten plants each. A mix of root and entire aerial plant parts were collected from ten tomato seedlings for each treatment. Plant material was flash frozen in liquid nitrogen and stored at − 80 °C until RNA extraction. The biological treatments consisted of one tomato line/three bacterial isolates and sampling at 24 and 96 HPI. Sterile water-treated plants were collected at 0 HPI, representing the mock-inoculated control treatment (Table 1) Table 3. For a subgroup of these selected genes, primers were previously validated in distinct pathosystems involving tomato, as described in the literature (Table 3). Additional in silico amplification assays using these primer pairs were performed using sequence information derived from the transcriptome data in the present study.
Evaluation of amplification efficiency. RT-qPCR reactions were performed using the Universal SYBR  www.nature.com/scientificreports/ primers (0.2 µM each), 1:20 cDNA (2.0 µL) and adjustment with nuclease-free water to a final volume of 10 µL. Amplification conditions were as follows: initial hold at 52 °C for 2 min, denaturation at 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. Each gene was analyzed in three biological replicates, with technical triplicates included per repetition. Primer specificity was also double-checked by confirming the presence of a single peak on melting curves in RT-qPCR analysis (Supplementary Figure S1), together with agarose gel analysis of all amplicons obtained via RT-PCR. Supplementary Figure S2 shows an example of specificity for the tested gene UBI3.
Stability analysis and validation of candidate reference genes. Twenty comparisons were analyzed amongst the seven biological samples from the 'Hawaii 7996'-R. solanacearum pathosystem (Table 1). These sets comprised one incompatible (one avirulent R. solanacearum isolate and 'Hawaii 7996') and two compatible (viz. two virulent/compatible R. solanacearum isolates and 'Hawaii 7996') interactions, as well as the mock-inoculated samples at 24 and 96 HPI, with three replicates (Table 1). Analyses were based on the values of the quantification cycle (Cq), calculated according to the efficiency of the reaction for each primer, estimated using LinReg software 45 . The Cq value represents the number of amplification cycles during the exponential PCR amplification phase that are required to reach a default threshold value for detection. The Cq average and variation for each gene were analyzed with Box-plot graph in the program SigmaPlot v. 12.5 (Systat Software, San Jose, CA, USA). Validation of the candidate reference genes was performed using three different analytical algorithms: BestKeeper 46 and NormFinder 47 in Excel, and geNorm using qBASE software, version 3.0 (Biogazelle, Zwijnaarde, Belgium), according to the default parameters of the software and as previously described 48 .
Performance of the primers for each candidate reference gene. Expression stability was examined across 'Hawaii 7996' samples, at different times after inoculation and with compatible and incompatible R. solanacearum isolates. Melting curve analysis was performed by RT-qPCR after 40 cycles of amplification. The specificity of the primer pairs for each candidate gene was confirmed through both visualization of single amplicons after standard PCR assays, together with RT-qPCR dissociation curve analysis and observation of single peaks, following calculations with LinReg PCR (Supplementary Figure S1). RT-qPCR amplification efficiencies from the standard curves ranged from 86 to 101% (Table 3).

Data availability
All data generated and analyzed during the study is included in the published article and Supplementary Information files.