Validation of reference genes for the normalization of RT-qPCR gene expression in Bacillus siamensis 1021 grown in different culture media.

Background and
Objectives: House-keeping genes are generally selected as reference genes in gene expression analysis. However, some genes may not be stably expressed across all experimental conditions. Thus, this study aimed to validate seven house-keeping genes for gene expression analysis in Bacillus siamensis 1021. Materials and
Methods: Strain 1021 was grown in potato dextrose broth, nutrient broth and mineral salt medium. Reverse-transcription quantitative PCR was used to determine Cq values of seven reference genes including gyrA, gyrB, ssb and dnaB, rpsU, gat_Yqey and udp in these media. Expression stability of these genes was analyzed, using geNorm and Normfinder applications. The target gene ftsZ was used for assessment of the best candidate genes.
Results: Based on geNorm and Normfinder, ssb was the most-stably expressed gene, while udp was the least-stably expressed gene. Pairwise variation indicated the combination of ssb, gyrA, gyrB and gatB_Yqey was suitable for the normalization of ftsZ expression. ftsZ expression in potato dextrose broth and mineral salt medium was higher than that in nutrient broth. In contrast, the normalization against udp resulted in an under- and overestimation of ftsZ expression in potato dextrose broth and mineral salt medium, respectively.
Conclusion: The combination of ssb, gyrA, gyrB and gatB_Yqey was the best candidate for normalization of target gene expression in B. siamensis 1021 in these media. This study emphasized the significance of reference gene validation for gene expression analysis and provided a guideline for future gene expression studies in B. siamensis.

INTRODUCTION

Reverse-transcription quantitative PCR (RT-qP-CR) is a relatively simple and inexpensive technique that is generally used to determine the expression level of target genes. It has been applied for research in various field, e.g., agricultural science, biomedical science, molecular genetics, etc. It is also used for assessing differential gene expression under different conditions in various tissues and organisms. The expression level of a target gene is normalized against reference genes, whose expression is relatively stable across tested conditions (1). In previous studies, house-keeping genes were selected as reference genes because of their constitutive expression. They were also required for basic growth and cellular functions under various conditions (2). However, they may be differentially expressed under different growth conditions or in distinct organisms. For example, in Pseudomonas brassicacearum GS20, gmk expression stability was less affected by changes in temperature, while gyrA was the most suitable reference gene under different iron concentrations in the medium (3). In contrast, the more stably expressed genes in Pseudomonas aeruginosa were proC and rpoD (4). Thus, the minimum information for publication of quantitative real-time PCR experiments (MIQE) indicated that all reference genes must be tested for their expression stabilities across cell and tissue types as well as experimental conditions (5). The use of a single reference gene for normalization must be avoided (5). This is excepted only when a gene was provided with evidence of its invariable gene expression under the experimental design (5).

Various algorithms and applications were available for testing gene expression stability of candidate genes. These included geNorm, Normfinder, Bestkeeper and RefFinder (6). Among these, geNorm and Normfinder were widely used in previous studies. It analyzed the expression stability of a given set of reference genes. The minimum number of reference genes was also determined, based on the calculated geometric mean (7). Normfinder used the model-based approach to analyze the expression variation of candidate genes (8). It calculated the stability value of each candidate gene and determined the best combination of two candidate genes (8). A previous study in Acanthamoeba spp. showed that geNorm and Normfinder consistently indicated that the 18S rRNA and HPRT genes were the most stably expressed reference genes (6). However, geNorm suggested different numbers of reference genes to be used in different Acanthamoeba spp. strains (6).

Bacillus siamensis was first described as a Gram-positive, rod-shaped bacterium isolated from salted crabs in Thailand (9). Based on a previous genomic study, it was found closely related to Bacillus velezensis and Bacillus amyloliquefaciens (10). Several strains of B. velezensis and B. amyloliquefaciens were demonstrated for their antagonistic and biocontrol activities against various phytopathogens and plant diseases (11–13). Similar to its relatives, strains of B. siamensis were recently proposed as potential biocontrol agents against many pathogenic fungi, e.g. Alternaria alternata, Rhizoctonia solani, Verticillium lateritium, Macrophomina phaseolina (14–16). The antagonistic activities of several B. siamensis strains were derived from the production of nonribosomal peptides (NRPs) (14, 17). This class of lipopeptides was synthesized by modular nonribosomal peptide synthetases (NRPSs) (18). Previous studies in other Bacillus species suggested that NRP production correlated with expression of NRPS genes (19–20). Several regulatory genes of NRP biosynthesis were also determined (21–22). Similar to NRPS genes, expression levels of these regulatory genes affected NRP biosynthesis in B. subtilis ZK0 (23). Thus, gene expression analysis could provide a better understanding on the regulation of NRP production. B. siamensis 1021 was previously isolated and identified, based on the phylogenetic analysis of the 16S rRNA gene (24). The production of antibiotic compounds of this strain against Pyricularia oryzae, the causative agent of rice blast, depended on culture media (24). However, little was known about the regulation and expression of NRPS genes in B. siamensis. For accurate quantification of NRPS genes in strain 1021, reference genes must be validated. Thus, the aim of the present study was to identify the most stably expressed reference genes among seven reference genes, including gyrA, gyrB, ssb, dnaB, rpsU, gat_Yqey and udp. The reference genes could be used for the future analysis of NRPS gene expression in this strain. This study also provided a guideline for gene expression analysis in other B. siamensis strains.

MATERIALS AND METHODS

Bacterial strain and culture.

B. siamensis strain 1021 was grown in nutrient agar (NA; Sisco Research Laboratories, India) at 30°C for 24 hours to obtain single colonies. A single colony was used for inoculation in 10 ml of nutrient broth (NB; Sisco Research Laboratories, India) and grown at 30°C for 24 hours, on a rotary shaker at 160 rpm. The light absorbance of the cell suspension was determined at the 600 nm wavelength (OD600 ). Cell concentration was adjusted to 0.1 OD600, and 100 μL of this suspension was used for inoculation in 50 ml of NB, potato dextrose broth (PDB; Himedia, India) and mineral salt medium (MSM) (17). The bacterial cultures were grown at 30°C for 60 hours on a rotary shaker at 160 rpm. The cultures were sampled at 8, 12, 24, 36, 48 and 60 hours post inoculation for growth determination. Cell suspension was serially diluted and grown on NA at 30°C for 18 hours for colony enumeration.

RNA extraction and cDNA synthesis.

For RNA extraction, 2 mL of cell suspension in PDB, NB and MSM were collected at 48 hours post inoculation. The cells were harvested by centrifugation at 13,500 rpm for 10 minutes. Favorprep Tissue Total RNA mini kit (Favorgen, Taiwan) was used for RNA extraction, according to the manufacturer’s protocol. RNA solution was treated with RNase-free DNase I (Geneaid, Taiwan) and purified with RNA Cleanup kit (Geneaid, Taiwan), according to the manufacturer’s protocols. RNA quality and concentration was determined using NanoDrop spectrophotometer (ThermoFisher Scientific, USA). Purified RNA solution was used for agarose gel electrophoresis to determine the absence of contaminating genomic DNA. cDNA was synthesized from one μg of RNA and random hexamer primers, using iScript Select cDNA Synthesis kit (Bio-Rad, USA). The temperature profile followed the manufacturer’s protocol.

Reference gene selection and primer design.

Whole-genome sequence of strain 1021 (assembly accession number: GCF_014534735.1) was obtained from the GenBank database. The annotated nucleotide sequences of gyrA¸grB, dnaB, rpsU, ssb, gatB_Yqey, udp and ftsZ genes were obtained. These genes functioned in DNA replication (gyrA, gyrB, ssb and dnaB), translation (rpsU and gat_Yqey) and the epimerization of UDP-N-acetylglucosamine (udp). Primers were manually designed and used for qRT-PCR (Table 1).

Table 1.

Primers used in this study.

Gene	Accession number	Primer sequence (5’-3’)	Tm (°C)	Product size (bp)	r2 0.998	PCR efficiency
gyrA	MBD0406835	F: CGATTAACCGGTCTGGAGCG	62.5	71		99.3%
		R: AGCTCGGCGATAAGCGCAAC	62.5		0.998
gyrB	MBD0406836	F: GCCGGCGGTAAATTTGACGG	62.5	72		99.4%
		R: TACGACAGACGCCCCTACAC	62.5		0.995
dnaB	MBD0406901	F: TTCAAGCGACAACCGCCAGC	62.5	64		104.5%
		R: CCGTGCCAGCGATTTTAGCG	62.5		0.998
rpsU	MBD0407962	F: GCTCTTCGTCGCTTCAAACGC	63.2	66		98.6%
		R: TTCGCGCTTTCTTGCTTCTTGC	62.1		0.998
ssb	MBD0407252	F: GGATACAACGAAGGAAACAGCG	62.1	61		99.3%
		R: GATTATCATTTTGGCCTCCGCC	62.1		0.993
gatB_Yqey	MBD0407961	F: CGAGCTCGAAAGCGGACATG	62.5	83		101.8%
		R: GTTAATCACGCTGCCGTCAGC	63.2		0.999
udp	MBD0407318	F: AATCCCGCAGTACGTGAGGC	62.5	55		98.6%
		R: GCACTCTGTCCGAATCAGCG	62.5		0.996
ftsZ	MBD0406176	F: GCCGAAGCCGCTAAAAAGGC	62.5	70		108.5%
		R: AAACACCTTGCGCTCCGTCG	62.5

Reverse-transcription quantitative PCR (RT-qPCR).

Equal volumes of all cDNA samples were pooled to obtain the 1× cDNA solution. This solution was used for preparation of five-fold serial dilutions (1/5×, 1/25× and 1/125×). All RT-qPCR reactions were performed, using Luna Universal qPCR Master Mix (New England Biolabs, USA), according to the manufacturer’s protocol. Detection of amplification products was based on SYBR Green fluorescence. Two μl of cDNA solutions of each concentration were used as the template. The temperature profile included one cycle of 95°C for 60 seconds and 40 cycles of 95°C for 15 seconds and 60°C for 30 seconds. Melt-curve analysis was done from 60°C to 95°C with a 0.5°C increment. The reactions were done on CFX96 Touch Real-Time PCR Detection System (Bio-Rad, USA). Cq (quantitation cycles) values were determined from three technical replicates and used for generating the standard curves with the logarithmic scale of the cDNA concentrations. r2 and PCR efficiency (E) (25) of the primers for each gene were determined. E was determined based on the slopes of the standard curves: E (%) = (10−(1/slope) − 1) × 100% (25).

Expression data and reference gene stability.

cDNA samples obtained from PDB, NB and MSM were diluted five folds. Cq values of the seven candidate genes were determined, using RT-qPCR. A single threshold baseline was set, using the CFS Manager software (Bio-Rad, USA). The Cq values in each sample was converted to relative quantities (RQ), using the formula RQ = E-ΔCq (6). ΔCq was determined as ΔCq = Cqmin – Cqsample, where Cqmin was the lowest Cq value for each gene and Cqsample was the Cq value of each sample. The RQ values were subsequently used for calculation of gene expression stability, using the geNorm (7) and Normfinder (8). Pairwise variation Vn/n+1 was performed, as described previously (7).

Validation of reference genes.

Relative ftsZ gene expression of strain 1021 in each medium was normalized against the selected reference genes, using the ΔΔCq method (26). Bacteria grown in NB medium was used as the reference group, because NB is the general medium for growing bacteria.

Statistical analysis.

All experiments were performed in triplicates. Statistical analysis was carried out, using IBM SPSS Statistics Version 28.0. Analysis of variance (ANOVA), followed by multiple comparisons with Tukey’s test, was used to determine statistically significant differences (P < 0.05).

RESULTS

Differential growth of strain 1021 in PDB, NB and MSM.

Strain 1021 grew differently in PDB, NB and MSM media (Fig. 1). Growth of strain 1021 in these media was comparable up to 24 hours post inoculation. At 36 hours, bacterial growth in PDB (2.3 + 0.14 ×108 CFU/mL) became considerably higher than that observed in the other two media. Additionally, growth in MSM medium (1.1 + 0.14 × 108 CFU/ml) became significantly higher than that in NB (5.8 × 107 CFU/mL). At 48 hours, growth in PDB slightly increased to 2.5 + 0.17 × 108 CFU/mL and was approximately four and two folds of that in NB (5.8 + 0.06 ×107 CFU/mL) and MSM (1.08 + 0.11 × 108 CFU/mL), respectively.

Fig. 1.

Growth of strain 1021 in PDB, NB and MSM media at 30°C on a rotary shaker at 160 rpm. The cultures were grown for 60 hours. Error bars indicate standard deviations. Different letters indicate statistically significant differences between samples at the same growth period.

Reference gene analysis for RT-qPCR.

Seven reference genes, gyrA, gyrB, dnaB, rpsU, ssb, gatB_Yqey and udp, were selected for the analysis of their gene expression stability in the three media. Based on the standard curves, the r2 and PCR efficiency values of the primers ranged from 0.993 to 0.999 and from 98.6% to 104.5%, respectively (Table 1). From the melt-curve analysis, a single peak and a single Tm value was obtained from each primer pair (Fig. 2). The Tm values ranged from 80.0°C to 82.5°C. This indicated there was only one type of PCR product amplified in each PCR reaction.

Fig. 2.

Melt-curves of the amplification product from the candidate genes gyrA, gyrB, dnaB, rpsU, ssb, gatB_Yqey and udp and the target gene ftsZ. A single peak indicated the presence of one type of amplification product in each reaction.

Cq values of the house-keeping genes of strain 1021 grown in PDB, NB and MSM were determined for the calculation of gene expression stability, using geNorm and Normfinder programs (Table 2). Based on geNorm, the expression stability values (M) of each reference gene ranged from 0.672 (ssb) to 1.314 (udp) (Table 2). Low M values indicated high levels of gene expression stability, while high M values signified low levels of gene expression stability (15). Thus, ssb was the most stably expressed gene, while udp was the least stably expressed gene. Stepwise exclusion of the least stably expressed reference genes was performed by the program geNorm. The combination of gyrA and gyrB displayed the lowest M value (Fig. 3). To determine the number of required reference genes, pairwise variation (Vn/n+1) was calculated. The variation level at 0.15 was previously proposed as the cut-off level where lower variation levels indicated the inclusion of an additional reference gene was not required (7). We found that V4/5 (0.13) and V5/6 (0.14) were below the cut-off level (Fig. 4). This suggested that the combination of the four most-stably expressed genes, ssb, gyrA, gyrB and gatB_Yqey, was the best candidates for normalization of target gene expression. Consistently, Normfinder (8) suggested that ssb (stability value = 0.086) and udp (0.740) were the most- and the least-stably expressed genes (Table 2). Additionally, the combination of ssb and gyrB was suggested by Normfinder with the stability value at 0.128.

Table 2.

Expression stability of seven reference genes obtained from geNorm and Normfinder.

Gene	Expression stability
	geNorm	Normfinder
ssb	0.672	0.086
gyrA	0.796	0.373
gatB_Yqey	0.828	0.381
rpsU	0.908	0.476
gyrB	0.709	0.258
dnaB	1.053	0.545
udp	1.314	0.740

Fig. 3.

Ranking of the reference genes based on the average expression stability (M) after stepwise exclusion of the least-stably expressed gene by geNorm.

Fig. 4.

Pairwise variation (Vn/n+1) was calculated with geNorm. The analysis determined the optimum number of reference genes required for accurate normalization of target gene expression. Vn/n+1 values lower than the cut-off value (0.15) indicated the addition of additional reference genes would not yield a significant improvement of normalization.

ftsZ expression analysis in strain 1021.

The r2 and PCR efficiency of the ftsZ primers were 0.996 and 108.5%, respectively (Table 1). The melt-curve analysis revealed a single peak and a single Cq value which suggested there was only one PCR product from the ftsZ primers (Fig. 2). Relative gene expression of ftsZ in PDB, NB and MSM was normalized against the combination of four genes (ssb-gyrA-gyrB-gatB_Yqey) (Fig. 5). ftsZ expression of strain 1021 in PDB (2.40 + 0.37) was comparable to that in MSM (2.37 + 0.65), whereas the expression level in NB (1.00 + 0.10) was significantly lower. A similar result was observed with the normalization against the combination of two genes (ssb-gyrB) that was suggested by Normfinder. In contrast, normalization against udp, the least-stably expressed gene, resulted in a different expression pattern. ftsZ expression in MSM (3.58 + 1.73) was overestimated, while an underestimation of ftsZ expression was observed in PDB (1.02 + 0.16). ftsZ expression in PDB was insignificantly different from that in NB.

Fig. 5.

Relative expression levels of the ftsZ gene of strain 1021 grown in PDB, NB and MSM, when normalized against different reference genes. Different letters indicate statistically significant differences (P < 005) between media with the same reference genes.

DISCUSSION

RT-qPCR is an important tool for studying gene expression and regulation. Reliable and stably expressed reference genes are crucial for an accurate estimation of target gene expression levels. The MIQE guideline suggested all reference genes used in an expression analysis must be validated for their expression stability in the experimental design (5). Previously, the expression stability of eleven genes were determined in Bacillus cereus ATCC 14579 and Bacillus thuringiensis subsp. konkukian 97-27. Based on geNorm and Normfinder programs, the study reported that the gatB_Yqey gene could be used as a single reference gene for normalization of gene expression throughout the bacterial life cycle in brain heart infusion medium (27). In contrast, another study of B. subtilis used the 16S rRNA gene as the reference gene for quantifying expression of genes involved in bacterial adaptation to phosphate-limiting conditions (28). However, a survey on the use of 16S rRNA as the reference gene indicated its poor expression stability in many validation studies of various bacterial species (29). The survey also found that genes belonging to DNA replication were experimentally validated in a higher proportion of previous studies (29). Additionally, up to the time of the preparation of this manuscript, no reference genes were not yet experimentally validated in other B. siamensis strains. Thus, for gene expression analysis of B. siamensis 1021, we aimed to examine the expression stability of various house-keeping genes that were involved in DNA replication (ssb, gyrA, gyrB, dnaB) and other cellular processes (rpsU, gatB_Yqey, udp).

Based on the expression stability levels obtained from geNorm and Normfinder, we found that the ssb, gyrA, gyrB and gatB_Yqey genes were the four most stably expressed genes. The V4/5 variation level (0.13) was below the cut-off level (7). It also suggested four genes were required for the normalization of target gene expression in strain 1021. The udp gene was the least stably expressed in this study, based on geNorm and Normfinder. In contrast, udp was the second stably expressed gene in B. cereus ATCC 14579 and B. thuringiensis subsp. konkukian 97-27 at various time point of their growth in brain heart infusion medium (27). This distinction emphasized the significance of reference gene analysis prior to a gene expression study.

Strain 1021 grew differently in PDB, NB and MSM. A higher growth was observed in PDB and MSM, compared to that in NB medium, after 24 hours post inoculation. This was rather counter-intuitive, since NB is often used for growing bacteria. However, several previous studies used PDB for growing strains of B. pumilus, Bacillus natto and B. subtilis for bioactive compound production (30–32). MSM was also included in this study because it was used for the optimization of secondary metabolite production by B. siamensis SCSIO 05746 (17). FtsZ is an initiating protein for the formation of the Z ring at the site of cell division in bacteria and directly involved in bacterial growth (33). In Escherichia coli, its expression oscillated during the cell cycle and reached the maximum during DNA replication (33). Thus, this gene was selected as the target gene for validation of the selected reference genes.

Based on the pairwise variation analysis by geNorm, the combination of ssb, gyrA, gyrB and gatB_Yqey were tested for the normalization of ftsZ expression. Based on this strategy, ftsZ expression was considerably higher in PDB and MSM, compared to NB. The normalization against the combination of two reference genes suggested by Normfinder (ssb-gyrB) also yielded a similar result. This ftsZ expression pattern positively correlated with strain 1021 growth in PDB, NB and MSM. It also suggested that the combination of the four reference genes was a suitable normalization strategy in this experimental setting. This was supported by previous studies which showed that the level of ftsZ expression changed during cell cycle. In Prochlorococcus, ftsZ mRNA level reached the maximum during DNA replication. Its protein level also peaked at the time of cell division (34). Another study demonstrated the growth inhibitory effect of bacilysin, obtained from Bacillus amyloliquefaciens FZB42, on Microcystis aeruginosa (35). It showed that ftsZ transcription in M. aeruginosa was reduced during the prolonged treatment of bacilysin. Our expression analysis suggested PDB and MSM had a stimulating effect on both strain 1021 growth and its ftsZ expression. It also supported our proposal of the combinations of the reference genes ssb, gyrA, gyrB and gatB_Yqey for normalization of target gene expression in this strain. These reference genes were potential candidates for future gene expression analysis in strain 1021 and other B. siamensis strains.

CONCLUSION

Seven reference genes were tested for their expression stability in Bacillus siamensis 1021, grown in PDB, NB and MSM. The combination of ssb, gyrA, gyrB and gatB_Yqey was found as the best candidate for the normalization strategy in this experimental design. This was supported by the correlation between strain 1021 growth and ftsZ expression in these media. In contrast, we found that the use of the least stably expressed gene udp resulted in a different pattern of ftsZ expression level. Our study indicated the significance of reference gene analysis and provided the candidates for further gene expression studies in B. siamensis 1021.