Selection and validation of reference genes for gene expression studies in Pseudomonas brassicacearum GS20 using real-time quantitative reverse transcription PCR

Pseudomonas brassicacearum GS20 is an antagonistic strain of bacteria recently isolated from the rhizosphere of Codonopsis pilosula. No validated reference gene has been indentified from P. brassicacearum to use in the normalization of real-time quantitative reverse transcription-PCR data. Therefore, in this study, nine candidate reference genes (recA, gyrA, rpoD, proC, gmk, rho, 16S, ftsz, and secA) were assessed at different growth phases and under various growth conditions. The expression stability of these candidate genes was evaluated using BestKeeper, NormFinder and GeNorm. In general, the results showed rho, rpoD and gyrA were the most suitable reference genes for P. brassicacearum GS20. The relative expression of iron-regulated gene (fhu) was normalized to verify the reliability of the proposed reference genes under iron-replete and iron-limited conditions. The trend in relative expression was consistent with the change in siderophore production under different iron conditions. This study presents reliable reference genes for transcriptional studies in P. brassicacearum GS20 under the chosen experimental conditions.


Introduction
Real-time quantitative reverse transcription PCR (qRT-PCR) is a technique that is widely used to monitor changes in gene expression [1]. The method is regarded to be sensitive, specific, reproducible and accurate [2], even though factors such as the quality of extracted RNA, primer selection and the efficiency of cDNA synthesis can affect the outcome [3]. Generally, these variables are minimized by relative normalization, whereby the expression of the target gene is normalized to the expression of an internal control gene [4]. Thus, choosing an internal control gene is vital to the accurate determination of relative gene expression [5]. Housekeeping genes are usually chosen as internal controls to normalize real-time qRT-PCR data; however, recent reports have claimed that the expression of housekeeping genes may vary according to the genes, cell types and experimental conditions [6]. Confirming the stability of PLOS ONE | https://doi.org/10.1371/journal.pone.0227927 January 27, 2020 1 / 15 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Growth phases: P. brassicacearum GS20 was cultured in LB liquid medium at 30˚C and pH 7.0; bacterial cells were harvested in the lag phase (OD 600 = 0.2), exponential phase (OD 600 = 0.6) and stationary phase (OD 600 = 2.0). The culture time of these three growth phases was 6h, 12h and 36h respectively.
Temperatures: The strain was cultured overnight in LB liquid medium at three different temperatures (25˚C, 30˚C and 37˚C) and the pH was 7.0. The bacterial cells were collected during the exponential phase (OD 600 = 0.6).
pH values: The strain was cultured overnight in LB liquid medium under three pH values (pH 5.0, pH 7.0 and pH 9.0) and at 30˚C. The bacterial cells were collected during the exponential phase (OD 600 = 0.6).
Iron levels: The expression stability of the proposed reference genes was analyzed under different concentrations of iron in LB culture medium. The high iron group was cultured in LB containing ferrous iron (FeSO 4 ); the low iron group was cultured in LB containing an iron chelator (2,2 0 -dipyridyl); the control group was cultured in LB containing MnSO 4 , which provides divalent ions (Mn 2+ ), for the studies of Fur-mediated iron regulation. The final concentrations of FeSO 4 , 2,2'-dipyridyl and MnSO 4 were 100 μM each in LB culture media. P. brassicacearum GS20 was initially cultured for 8 h in LB liquid medium at 30˚C and pH 7.0; then, the bacterial cells were harvested 30 min after the addition of FeSO 4 , 2,2'-dipyridyl or MnSO 4 .
All of the abovementioned culture samples were conducted in triplicate. The bacterial cells were centrifuged (at 4˚C and 3074 g for 10 min) and collected in Eppendorf tubes before extraction of total RNA.

Total RNA extraction and cDNA synthesis
Total RNA was extracted from 5 mL of fresh culture using TRIzol reagent (Invitrogen, Carlsbad, CA, DNA). The RNA was dissolved in ribonuclease-free water and its purity and concentration were quantified using a NanoDrop spectrophotometer (ND-2000; NanoDrop Technologies Inc., Wilmington, DE, DNA). RNA integrity was analyzed by 3% agarose gel electrophoresis. Approximately 1 μg RNA was used for real-time PCR in each sample. Complementary DNA (cDNA) was synthesized with a PrimeScript RT reagent kit with genomic DNA (gDNA) Eraser (Perfect Real Time, Takara Biochemicals, Dalian, China), according to the manufacturer's instructions. The kit eliminated gDNA residues before the synthesis of first-strand complementary DNA by the following reaction: 5 × gDNA Eraser Buffer (2.0 μL), gDNA Eraser (1.0 μL), total RNA (� 1 μg), then RNase-free H 2 O was added up to 10 μL; the reaction system was carried out at 42˚C for 2 min. The first-strand cDNA was synthesized in a 20 μL volume reaction containing denatured reaction products (

Real-time-qRT-PCR assay
Prior to the real-time-qRT-PCR assays, the synthesized cDNA samples were diluted with RNase-free water. Real-time-qRT-PCR was performed with TB Green 1 Premix Ex Taq II (Tli RNaseH Plus, Takara Biochemicals, Dalian, China) and conducted on an ABI StepOnePlus Real-Time PCR System (PerkinElmer Applied Biosystems, Foster City, CA, USA). The reaction mixture included 2 μL of diluted cDNA, each primer (final concentration 400 nM), ROX reference dye (final concentration 500 nM) and 2 × TB Green qPCR SuperMix adjusted with nuclease-free water to a final volume of 20 μL. The primer sequences for the nine selected genes are listed in Table 1. The PCR procedure was as follows: 95˚C for 30 s, 40 cycles of 95˚C for 5 s and 60˚C for 33 s. The melting temperature-determining dissociation step was run at 95˚C for 15 s and 60˚C for 1 min and 95˚C for 15 s at the end of the amplification. The realtime-qRT-PCR reactions were carried out in triplicate for each cDNA sample.

Reference gene expression stability analyses
The expression stability of the nine selected reference genes was analyzed using three Microsoft Excel-based tools: BestKeeper [19], GeNorm [20] and NormFinder [21]. For GeNorm and NormFinder, the raw quantification cycle (Cq) values needed to be transformed into relative quantification data by the formula 2 (−ΔCT) . ΔCT included all minus the lowest Cq value among all samples for each reference gene. GeNorm evaluates candidate reference genes by the expression stability measurement (M) value. The lowest M value indicates the most stable expression. Generally, genes with M values �1.5 are considered not suitable as reference genes [20]. NormFinder ranks the stability values of all candidate reference genes, similarly to GeNorm; the reference gene with the lowest stability value is considered to be the most stable. NormFinder takes intragroup and intergroup expression variations into account to find the stability value for each candidate reference gene and this algorithm can avoid misinterpretations caused by co-regulated genes. BestKeeper assesses the stability of reference gene expression by the coefficient of variation [CV (% Cq)] and the standard deviation (SD) of the Cq values [SD (±Cq)]. BestKeeper is used to analyze the raw Cq values without any transformation. BestKeeper regards the reference gene with the lowest SD and CV-values to be the most stable. Genes with SD �1 are considered to be inconsistent by BestKeeper. Moreover, BestKeeper carries out numerous pairwise correlation analyses between each candidate reference gene and calculates the Pearson correlation coefficient (r) and the probability values. This program considers a higher r value to represent greater stability in gene expression.

Siderophore halo production under different iron levels
In our previous study, P. brassicacearum GS20 was shown to be strongly antagonistic to pathogenic fungi and to produce extracellular siderophores on Chrome azurol S (CAS) agar plates. CAS agar plates were prepared according to the method described by Schwyn and Neilands [22]. In order to investigate the production of siderophore halos under different concentrations of iron. FeSO 4 , MnSO 4 , and 2,2'-dipyridyl were added separately to CAS agar plates to final concentrations of 100 μM. P. brassicacearum GS20 was inoculated on these CAS agar plates and cultured at 30˚C for 12 h. The halos on the different CAS agar plates were photographed and measured using a microscope (Leica M205C, Buffalo Grove, IL, USA).

Validation of Fur-regulated gene (fhu) expression
We analyzed the genome sequence of P. brassicacearum GS20 and found the iron-regulated gene (fhu). In order to verify the reliability of the proposed reference genes in the relative expression normalization of fhu gene under iron-replete and iron-limited conditions, the fold changes of fhu were calculated by the 2 -ΔΔct method as the description of Livak and Schmittgen [23]. The relative expression results of fhu gene were statistically analysed under different iron contents. Turkey tests of one-way analysis of variance (ANOVA) were used determine the significance of any differences (P < 0.01).

Real-time qRT-PCR amplification efficiency and specificity of candidate reference genes
Based on previous reports on suitable internal control genes for Pseudomonas species, nine candidate reference genes were chosen for this study. Gene names and descriptions, primer sequences, PCR efficiencies, and amplicon sizes are listed in Table 1. The gene proC should be pointed out because its Cq values were >30 even when undiluted cDNA was used as the amplification template. Furthermore, we changed its primers twice and the results were the same. Therefore, we excluded proC from further study. Primers designed for the amplification were specific and the amplicon sizes were as expected (S1 Fig). All dissociation-curves revealed only single peaks, which indicated that there was no formation of primer-dimers or non-specific PCR products (S2 Fig). Amplification efficiency reflects the amplification rate of a template during the PCR reaction and could be determined by the slope of calibration curve. The acceptable range is generally between 90% and 110% [24]. For all of the candidate genes in this study, the PCR efficiencies extended from 96.78% (secA) to 107.24% (rho), which were in the acceptable range. Moreover, all correlation coefficients (R 2 ) ( Table 1) were reasonable and verified the quality of the primers designed for the real-time-qRT-PCR analysis.

Expression of candidate reference genes under different conditions
The candidate reference gene expression levels could be assessed by their Cq values in all samples. As shown in Fig 1, the mean Cq values ranged from 11 (16S) to 28 (gyrA). Low Cq values designate high expression levels and the opposite indicate low expression levels [25]. The 16S was the most abundantly expressed gene in all samples (11.27 ± 1.03, mean Cq ± standard deviation), while gyrA was the least abundantly expressed gene (28.17 ± 1.35).

Expression stability analysis by BestKeeper, NormFinder and GeNorm
The expression stabilities of the eight selected candidate genes were evaluated under different conditions and measured by three statistical algorithms: BestKeeper, NormFinder, and GeNorm.
BestKeeper performs repeated pairwise correlation and regression analyses for a given gene with all other genes [26]. The results of the Bestkeeper analysis for the eight genes were shown in Table 2. BestKeeper estimated the standard deviation (SD) of the Cq values [SD (± Cq)] and SD of the absolute regulation coefficients [SD (± x-fold)] and then considered whether a reference gene was acceptable. Generally, suitable reference genes should have values of SD [± Cq] < 1 and SD [± x-fold] <2 [27]. According to this standard, the candidate reference genes rho, 16S, gyr, and gmk passed this filter and gmk had the lowest SD [± Cq] (0.59) and SD [± x-fold] (1.50), while ftsz had both the highest SD [± Cq] (1.31) and SD [± x-fold] (2.49). Among the genes gmk, 16S, gyrA and rho, 16S exhibited the highest coefficient of variation (CV) value (7.77). BestKeeper performed pairwise comparisons and calculated coefficients of correlation (r). The r value is related to the stability of gene expression and the ideal r value is 1, which indicates the highest gene expression stability. Among rho, 16S, gyrA and gmk, the r value of rho (0.924) was the closest to 1, so it was the most stable gene. The gmk was the least stable gene and its r value was 0.783.
NormFinder also measures candidate gene stability. Based on this algorithm, the most stable candidate genes were usually those with the lowest stability values. The stability values of the eight genes in all samples are shown in Table 3. Under different pH values, ftsz (0.062) was the most stable gene, while recA (0.651) was the least stable. The ftsz (0.318) was ranked highest, whereas gmk (0.087) was the most stably expressed under different temperatures. At different growth phases, rho was the most stable gene with an M value of 0.024 and recA was  Selection and validation of reference genes for Pseudomonas brassicacearum GS20 using real-time-qPCR than for the other conditions. The M values of 16S and recA were the lowest (0.04), while the M value of secA was the highest (0.33). In analyzing all of the conditions together, gyrA and rpoD ranked as the most stable genes, with an M value of 0.43, while recA had the highest value for M at 1.28 and, therefore, showed the least stable expression under all conditions. GeNorm was also used to determine a suitable number of reference genes required for the precise normalization of target gene expression by pairwise variation (V value). The V values of the eight reference genes under different conditions are shown in Fig 3. For different pH values, temperatures, and iron levels, the V 2/3 value was lower than 0.15, which was below the proposed V value threshold of 0.15. This meant that the first two genes were enough to be used for the validation. At different grow phases, all the V values were more than 0.15. In analyzing all of the different conditions together, the V 5/6 value was 0.14; hence, the first five genes (gyrA, rpoD, gmk, rho, and ftsz) were enough to provide accurate validation and there was no need to introduce 16S to improve it. In fact, the V value is only used to guide the selection of the optimal number of reference genes and the proposed value (0.15) should not be taken as a strict cut-off [20]. Usually, the three best reference genes are used for a valid normalization strategy, which can result in a more reliable normalization than the use of only one reference gene.
Based on the algorithms of BestKeeper, NormFinder and GeNorm, the expression stability of the eight reference genes was not exactly in accordance. Therefore, we confirmed the rank of the eight reference genes using the geometric mean of the rankings from the three programs Selection and validation of reference genes for Pseudomonas brassicacearum GS20 using real-time-qPCR (Table 4). The result showed that rho was the most stable gene, followed by gyrA and rpoD and, finally, recA. Therefore, rho, gyrA, and rpoD were selected as the most stable internal reference genes for accurate validation in the subsequent study.

Validation of Fur-regulated gene variation under different iron levels
The genes rho, gyrA and rpoD were selected as the most stable reference genes to validate Furregulated gene variation under iron-replete and -limited conditions. As shown in Fig 4, fhu  Selection and validation of reference genes for Pseudomonas brassicacearum GS20 using real-time-qPCR expression was normalized by rho, gyrA, rpoD and their geometric means, respectively; the results of the four kinds of normalization were consistent. They all showed that, when compared with control group, fhu was significantly upregulated (P < 0.01), and its expression was increased 1.64-2.20-fold in the iron-limited group; while fhu was significantly downregulated (P < 0.01), and its expression was decreased 0.40-0.61-fold in the iron-replete group. The fhu expression differences were in accordance with the variations in siderophore yields under exposure to different levels of iron. We inoculated P. brassicacearum GS20 onto CAS plates containing different levels of iron; an orange halo zone surrounding the colonies signified the production of siderophores (Fig 5). The diameter of the halo zone was much smaller on the iron-replete plates (5,774.894 μm) than the iron-limited plates (13,546.008 μm).

Discussion
Quantification of gene expression contributes to the study of cellular responses to changes in conditions and provides an understanding of regulatory mechanisms. To obtain reliable gene expression results using real-time qRT-PCR, stable internal control genes should be selected and used for the accurate normalization of target gene expression [25]. In our study, P. brassicacearum GS20 had recently been isolated and showed strong antifungal properties, which highlighted its potential as a biocontrol strain. There are no previous reports of validated reference genes for P. brassicacearum. There were several studies relating to internal reference genes for Pseudomonas species, but the results were inconsistent. Savli et al. (2003) studied the expression stability of six housekeeping genes (ampC, fabD, proC, pbp-2, rpoD, and rpoS) in P. aeruginosa and found that proC and rpoD comprised the most stable pair [14]. Chang et al. (2009) investigated the expression stability of eight reference genes (rpoN, rpoD, dbhA, phaF, 16S, gst, lexA, and atkA) in P. putida mt-2 during its degradation of p-xylene and revealed that rpoN, rpoD, 16S, and atkA were suitable reference genes with highly Selection and validation of reference genes for Pseudomonas brassicacearum GS20 using real-time-qPCR stable expression [18]. Alqarni et al. (2013) analyzed the expression stability of 13 housekeeping genes (rpoS, proC, recA, rpsL, rho, oprL, anr, tipA, nadB, fabD, ampC, algD, and gyrA) under carbon starvation of P. aeruginosa and found rpoS to be the only stably expressed gene [28]. Inconsistencies in reported reference genes attributed to the lack of a universal reference gene for use in gene expression studies and, therefore, the expression stability of reference Selection and validation of reference genes for Pseudomonas brassicacearum GS20 using real-time-qPCR genes could vary under different experimental conditions [20]. Ideally, the expression stability of each candidate reference gene should be verified under each condition [29]. Other studies used 16S rRNA genes as internal references for validation [30,31], but transcript levels of 16S may depend on the state of the bacteria [32]. Also, rRNA genes are abundant and this could affect the accurate quantification of rare and less abundant mRNA transcripts [33]. Moreover, due to their high abundance, the cDNA used needs to be more dilute prior to real-time qRT-PCR reactions; thus, dilution errors may compromise the accuracy of the analyses. Therefore, it was essential to identify suitable internal reference genes for P. brassicacearum GS20.
In this study, we selected nine candidate reference genes (recA, gyrA, rpoD, proC, gmk, rho, 16S, ftsz and secA) to assess their expression stabilities. The programs GeNorm, NormFinder, and BestKeeper were used to analyze data from all the samples. The gene proC was excluded due to its low expression quantity, which was different to the result obtained by Savli et al [14]. This may be due to proC expression differences in the two strains. Based on our analysis of the results calculated by these three programs, we found that the most stable gene was not identical in all conditions but the least stable gene was similar in some conditions. For example, secA and recA performed as the least stable genes in most conditions (different growth phases, pH values and iron levels), which was in accordance with other studies [34,35]. So, if the tested samples are from very different treatments, it is important to choose the specific reference gene for validation under each condition [36]. Among the eight genes, the most stable gene was rho, as analyzed by NormFinder and BestKeeper, while GeNorm determined gyrA and rpoD to be the most stable. The pairwise variation determined using GeNorm showed that the top five reference genes (gyrA, rpoD, rho, gmk and ftsz) could be used as the minimum number required to carry out an accurate normalization of the real-time qRT-PCR data. Based on the different algorithms, the results obtained using GeNorm, NormFinder, and BestKeeper yielded slight differences. Others have reported gene stability ranking differences among these three algorithms [35]. Therefore, taking the three algorithms into consideration, we confirmed the ranking of the eight reference genes' expression stabilities by calculating the geometric mean of the individual rankings from each of the different algorithms. The result suggested that rho ranked first, followed by gyrA, rpoD, gmk, ftsZ, 16S, secA and recA. Using only one reference gene for validation may affect the accuracy, so deferring to the suggestion provided be the GeNorm result, we selected the three most stable genes (rho, rpoD and gyrA) for the use in validation studies.
To verify the reliability of the selected reference genes, we applied them individually and together to normalize the relative expression of iron-related genes under different iron conditions. They could all objectively reflect the changes in fhu expression, even though the foldchange values calculated by the selected reference genes showed some minor difference, which may have been caused by their minor stability differences. By analyzing the stability results, we found the results for fold-change were consistent with the results for stability ranking determined by NormFinder. The ranking of gene stability according to NormFinder was gyrA, rpoD, and rho. Iron-limited conditions could upregulate fhu gene expression in P. brassicacearum GS20 and stimulate the strain to produce more siderophores and obtain iron from the surroundings, thus maintaining normal metabolism. The production of siderophores is an important aspect in antagonizing phytopathogenic fungi, especially under iron-limited conditions [13]. Plant growth-promoting fluorescent pseudomonads could be beneficial as they produce extracellular siderophores under iron-limited conditions and thus utilize environmental iron efficiently and reduce the amount of iron available to other microorganisms, including pathogenic fungi, so as to inhibit their growth [37]. Apart from siderophore production, volatile organic compounds released by P. brassicacearum GS20 may also inhibit the growth of phytopathogenic fungi (unpublished). Therefore, study of the related antagonistic mechanisms is required. The internal reference genes selected in this study will be valuable to further research in this area.