Validation of reference genes for real-time quantitative polymerase chain reaction analysis in Lactobacillus plantarum R23 under sulfur dioxide stress conditions

Abstract: Background and Aims: Malolactic fermentation (MLF), accomplished by lactic acid bacteria (LAB), is indispensable for the production of some red wines. As the key additive for wines, sulfur dioxide (SO₂) frequently affects LAB growth and MLF development. We aimed to investigate suitable reference genes in Lactobacillus plantarum R23, an excellent LAB with good tolerance to SO₂, for a real‐time quantitative PCR analysis of SO₂ tolerance. Methods and Results: The expression of eight candidate reference genes under SO₂ stress conditions for L. plantarum R23 was analysed, and geNorm and BestKeeper software packages were employed to assess the reference genes. The expression of mleR1 and mleP1 under SO₂ stress conditions was used to confirm the results. The geNorm analysis indicated that the stability ranking of reference genes was rpoB > rpoC > recA > ldh > cfal > mtlR > asd2 > gpp. Both software packages suggested that four genes should be selected for accurate normalisation, and the validated experiment further confirmed the suitability of the reference genes revealed in this study. Conclusions: The genes rpoB, rpoC, recA and ldh were the most stably expressed reference genes and thus represent the reference gene combination that should be used to obtain the most accurate results for different SO₂ stress conditions in L. plantarum R23. Significance of the Study: This is the first detailed study to evaluate selected reference genes in L. plantarum R23 under SO₂ stress conditions, which should allow the quantification of bacterial gene expression levels under SO₂ stress conditions and also provide a foundation for the more accurate use of real‐time quantitative PCR in L. plantarum.

Lactobacillus plantarum R23; malolactic fermentation; real‐time quantitative PCR; reference genes; SO₂ stress

Real‐time quantitative PCR (RT‐qPCR) is a nucleic acid quantitative analysis method that was developed based on the traditional PCR technique and has the following characteristics: accuracy, high sensitivity, repeatability and high throughput. It has been widely used to analyse gene expression (Gachon et al. [9] , Huggett et al. [13] ). Although RT‐qPCR is a powerful technique for obtaining quantitative gene expression data, artificial differences may arise between samples that could affect the final analysis of gene expression (Bustin [6] , Bohle et al. [4] ). Recent studies have shown that the expression of commonly used reference genes under different conditions may be unstable (Gutierrez et al. [11] ). Stable expression analysis in the absence of a reference gene can have a great effect on the target gene. Small differences in the blind use of a reference gene as an internal control could make it difficult to measure gene expression, with possible errors or even opposite conclusions appearing (Volkov et al. [34] ). Therefore, the appropriate selection of reference genes is particularly important. Selecting the appropriate reference genes is the first step in analysing the expression of the genes of interest.

Housekeeping genes are widely used as reference genes to detect expression changes of the target genes during certain developmental stages or under environmental stress conditions. These housekeeping genes, however, cannot always be stably expressed under all physiological conditions as ideal reference genes (Scharlaken et al. [27] ). During a gene expression analysis, the most common reference genes cannot meet the requirements of accurate quantification. The stably expressed reference genes can be screened based on some statistical analysis software packages, such as geNorm (Vandesompele et al. [33] ), BestKeeper (Pfaffl et al. [24] ) and NormFinder (Andersen et al. [1] ).

The implementation of malolactic fermentation (MLF) is important for red wines because it reduces the acidity, brings biological stability and improves the sensory characteristics of the product (Fahimi et al. [8] , Bartowsky et al. [3] ). Malolactic fermentation is accomplished by lactic acid bacteria (LAB), such as Oenococcus oeni and Lactobacillus plantarum (Genisheva et al. [10] , Krista et al. [16] ). Sulfur dioxide (SO2) is the key additive for wine preservation (Jackowetz and Mira de Orduna [14] ) and is most frequently employed to control LAB growth and MLF development during winemaking because of its antioxidant and selective antimicrobial properties, especially against spoilage microorganisms (Carmen et al. [7] ). Lactobacillus plantarum R23 is an excellent LAB with an SO2 tolerance that can reach 120 mg/L. It was acquired from fermented loquat wine and is a desirable LAB because it typically produces pleasant substances, aromas and flavours (He et al. [12] ), but its SO2 stress control mechanism remains obscure. To analyse the effect of SO2 on the expression of L. plantarum R23 genes through RT‐qPCR, the chosen reference genes should be unaffected by SO2 stress conditions. Adequate reference genes, however, for the analysis of SO2 tolerance gene expression under SO2 stress conditions by RT‐qPCR have not been fully described. Thus, this study aimed to investigate the stability of selected reference genes in L. plantarum R23 by surveying suitable reference genes for the RT‐qPCR analysis of the genes involved in SO2 tolerance. Two software programs, geNorm and BestKeeper, were employed to predict which gene or gene combination would be best suited as a reference for L. plantarum R23 via RT‐qPCR under variable SO2 stress conditions. This study provides the theoretical basis for engineering lactobacillus with high SO2 resistance and provides technical support for the development and application of superior MLF strains in wine.

Lactobacillus plantarum R23 (NCBI: HQ658056) was acquired from natural fermented loquat wine and was stored at −70°C in de Man, Rogosa, Sharpe (MRS) broth containing 16% glycerol (He et al. [12] ). Isolated colonies from MRS agar streak plates were picked and grown twice at 30°C in MRS broth (pH 6.2) for 18–20 h. Then, the strain was separately inoculated into MRS broth (pH 6.2), and 0, 60 and 120 mg/L SO2 (in the form of sulfuric acid, containing 6% SO2) were added such that the initial concentration of bacteria was approximately 107 CFU/mL. The cultures were incubated at 30°C for 12–14 h. All experiments were in triplicate.

First, a 10 mL culture (bacteria and broth) of the three SO2 treatments was collected in sterile 50 mL glass bottles, and 10 mg/L lysozyme was added separately for cell wall lysis at 20°C for 15 min. Thereafter, the supernatants were transferred to sterile 2 mL tubes and centrifuged at 16 500 × g for 5 min. The broth was discarded, and the bacterial pellet was snap‐frozen in liquid N2 and stored at −80°C until required for RNA isolation.

Approximately 50 mg of bacterial pellet was transferred to a sterile 2 mL Eppendorf tube and was crushed to a fine powder under liquid N2. The RNA was isolated using the RNA isolation kit (Tiangen Biotech, Beijing, China) according to the manufacturer's instructions, and RNase water was used where necessary. To avoid residual DNA contamination, RNA samples were treated with DNase I (Tiangen Biotech). The RNA was eluted into 40 μL of sterile elution buffer and stored at −80°C until required.

The RNA concentration was determined using UV spectrophotometry (NanoDrop ND‐1000; NanoDrop Technologies, Wilmington, DE, USA), and the integrity of the RNA was determined by 1% agarose gel electrophoresis. The RNA of L. plantarum R23 with no SO2 treatment was sent to Shenzhen BGI (Beijing Genomics Institute, Shenzhen, Guangdong Province, China) for Illumina HiSeq 2500 (Illumina, San Diego, CA, USA) high‐throughput sequencing analysis.

The RNA isolated from each treatment was converted to cDNA using the Thermo Scientific cDNA synthesis kit (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer's instructions. All cDNA was stored at −80°C until required for RT‐qPCR analysis. The reaction mixtures contained 2 μL of 5× reaction buffer, 1 μL RiboLock RNase Inhibitor (20 U/L), 2 μL dNTP (deoxy‐ribonucleoside triphosphate) Mix (10 mmol), 1 μL RevertAid M‐MuLV RT (200 U/μL), 10 μL total RNA (0.2 g/L) and sterile H2O, comprising a final volume of 20 μL. The cDNA synthesis was carried out by incubating the complete reaction mix for 60 min at 42°C, followed by 5 min at 70°C to denature the reverse transcriptase.

The ldh, recA, mtlR, rpoB, rpoC, gpp, cfa1 and asd2 genes were selected as candidate reference genes according to the HiSeq 2500 (Illumina Inc. USA) high‐throughput sequencing results. The PCR primers of these candidate reference genes were designed using the Primer 5.0 software (Table ). Primers were synthesised by Takara Bio (Dalian, Liaoning Province, China). The gene sequences were amplified by conventional PCR amplification; the PCR reaction conditions were 5 min at 95°C for initial denaturation, 30 s at 94°C final denaturation and annealing for 30 s at 56°C, followed by 50 cycles of extension at 72°C for 30 s and then extension at 72°C for 4 min.

Candidate genes and their primers for real‐time quantitative polymerase chain reaction.

Using cDNA of L. plantarum R23 under different SO2 stress as the template, the standard melting curve was obtained according to a concentration gradient of 1/3–1/35 dilutions. All RT‐qPCR reactions were performed on the Applied Biosystems 7500 Real‐Time PCR System (Applied Biosystems, Carlsbad, CA, USA) according to manufacturer's instructions. The RT‐qPCR reaction conditions were 30 s at 95°C for the initial denaturation, followed by 40 cycles of denaturation at 95°C for 5 s and annealing at 60°C for 34 s. Each individual treatment was in triplicate, which means that nine RT‐qPCR reactions were performed for each reference gene.

The amplification efficiency of the candidate reference genes, the slope of the standard curve, the target fragment solution temperature chain and the melting curve analysis were analysed with the 7500 software (V 2.0.6) (Applied Biosystems).

The expression stability of the selected reference genes in L. plantarum R23 under different SO2 stress conditions was determined by measuring the mRNA expression of the reference genes. The transcription level of the selected genes across different experimental manipulations was compared by converting the average Ct value of each duplicate reaction to raw data for subsequent analysis with the geNorm (Vandesompele et al. [33] ) and BestKeeper (Pfaffl et al. [24] ) programs.

The malolactic enzyme gene (mle) plays an important role in regulating MLF and includes mleA, mleP and mleR (Miller et al. [22] ). The selected reference genes were validated by quantifying the gene expression level of mleR1 and mleP1 in L. plantarum R23 in response to 0, 60 and 120 mg/L SO2 using the most stable reference gene(s) and the most unstable gene in the same RT‐qPCR conditions mentioned above. l‐Malic acid (1.5 g/L) was added to the MRS broth, and the pH value was adjusted to 6.2 with KH2PO4 and NaH2PO4. The primer pairs of mleR1 were F: 5′‐CGTTGTGCTAACGGATGTTG‐3′, R: 5′‐TTTGGCTACGTAACGCTGTG‐3′, and those of mleP1 were F: 5′‐CGTTGCCAGGATTACGTTTT‐3′, R: 5′‐AATGGTGGTATTGGCCATGT‐3′. Results were analysed using the comparative critical threshold (ΔΔCT) method (Livak and Schmittgen [18] ).

The agarose gel electrophoresis results of total RNA for L. plantarum R23 under different SO2 stress are shown in Figure . The 23S molecular mass of L. plantarum R23 ranged from 800 to 1200 bp. The 16S molecular mass was approximately 700 bp. The results of repeated experiments demonstrated that all RNA samples had good integrity and that there was no obvious degradation. The RNA concentration of the sample was 120–293 ng/μL, as measured by UV spectrophotometry. The ratio of OD260/OD230 ranged between 1.82 and 2.16. This ratio indicated that the RNA had high purity and could be used for subsequent high‐throughput sequencing and RT‐qPCR analysis.

The standard curve of each reference gene was obtained by PCR amplification using a gradient dilution, and the relevant parameters of all primers are shown in Table . Eight candidate reference genes were used for the data analysis, and the results displayed a relatively wide expression level range, from the lowest mean Ct value (21.71) in rpoC to the highest (26.46) in asd2. The PCR amplification efficiency (E) of the eight genes varied from 90.1 to 98.7%, and all correlation coefficients (R2) were above 0.996, which indicated that the results were accurate and reliable. Simultaneously, the melting curves of different SO2 stress culture samples all appeared only as a single peak, and there were no impurity peaks in the absence of primers and dimers. Furthermore, the melting peak corresponding to the temperature was close to the theoretical temperature of the amplified product (Figure ). This result indicated that the amplification products of the reference gene PCR were a single‐specific product.

Information used for the real‐time quantitative polymerase chain reaction analysis of candidate reference genes.

All the eight candidate reference genes exhibited an excellent efficiency value, and as a result they were examined further to validate their suitability for use as reference genes under SO2 stress conditions.

The gene stability measure ranking of the eight reference genes after determination by geNorm was gpp, asd2, mtlR, cfal, ldh, recA, rpoC and rpoB (from least to most stable). All the eight genes reached high expression stability with a low M‐value (average expression value) below the default limit of M = 1.5 (Figure ). GeNorm analysis also revealed that the pairwise variation value V4/5 was lower than the others (Figure ). The increasing variation in this ratio corresponds to decreasing expression stability due to the inclusion of a relatively unstable fifth gene. Thus, four genes (ldh, recA, rpoC and rpoB) are necessary to normalise RT‐qPCR data for SO2 stress gene expression in L. plantarum R23.

The gene expression variation was calculated for all eight candidate reference genes based on the Ct values and was displayed as the SD and coefficient of variance (CV). The BestKeeper software highlighted rpoB as the reference gene that exhibited the least overall variation from the list of eight candidate genes (Table ), which represented an acceptable 1.08‐fold change in expression. The highest SD values (less stable) to the lowest SD (most stable) are as follows: gpp > asd2 > mtlR > cfal > ldh > recA > rpoC > rpoB. Of the expressed reference genes, the four most stable were identical to those determined using geNorm and had the same stability ranking. From the two programs, ldh, recA, rpoC and rpoB could be considered the best reference genes.

BestKeeper descriptive statistical analysis for reference genes based on threshold cycle values.

1 AM (Ct), arithmetic mean of Ct; Ct, threshold cycle; CV (% Ct), coefficient of variance expressed as a proportion of the Ct level; GM (Ct), geometric mean of Ct; Min (Ct) and Max (Ct), extreme values of Ct; Min (x‐fold) and Max (x‐fold), extreme values of expression levels expressed as absolute x‐fold over or under coefficient; n, number of samples; SD (±Ct), standard deviation of the Ct; SD (±x‐fold), standard deviation of the absolute regulation coefficients.

To further validate the selected housekeeping genes, the relative expression level of mleR1 and mleP1 genes in L. plantarum R23 under different SO2 treatments was investigated using one or two of the most stable reference genes. The result revealed that the expression level of both mleR1 and mleP1 decreased progressively with increasing SO2 concentration, and showed similar changes with a slight difference when using rpoB, rpoC, recA and ldh alone or the combination (rpoB + rpoC, rpoB + ldh) as reference gene (s) for normalisation (Figure ). These changes, however, were completely obscured during normalisation using the least stable reference gene (gpp). These results illustrated the adverse effect of using an unsuitable reference gene for normalisation.

Although RT‐qPCR is a powerful technique for obtaining quantitative gene expression data, it is affected by several factors, such as the quantity of the initial materials, the quality of the RNA, the efficiency of cDNA synthesis, primer performance and the methods used for statistical analysis (Maroufi et al. [21] ). Normalising the expression of the target gene to one or more reference genes is an effective solution to eliminate such effects (Udvardi et al. [31] ). To obtain reliable and accurate results, it is critical to ensure that the expression of the reference genes themselves is not affected by changes in the experimental conditions. Selecting appropriate reference genes depends on experimental conditions. Thus, the accuracy of results from qRT‐PCR analysis largely depends on the reference genes used. Recent studies, however, indicate that the traditional reference genes are not always stably expressed when tested in other species or in a wider range of experimental treatments (Mukesh et al. [23] , Jian et al. [15] ). This means that even the most stable reference gene(s) reported in a species should be validated when used in other species under study or under new experimental conditions. Therefore, it is advisable to validate the expression stability of candidate reference genes under specific experimental conditions prior to their use in RT‐qPCR normalisation, rather than using reference genes published elsewhere (Remans et al. [25] , Lee et al. [17] ).

In this study, the rpoB, rpoC, ldh, recA, mtlR, gpp, cfa1, asd2, 16S r RNA and gapdh genes were initially chosen as candidate reference genes according to Illumina HiSeq 2500 high‐throughput sequencing results of L. plantarum R23, and they had been used in L. plantarum for normalisation of RT‐qPCR data. For example, 16S r RNA (Zhao et al. [35] , Bron et al. [5] ), rpoB and recA (Marco and Kleerebezem [20] , Zhao et al. [35] ) were used as reference genes in L. plantarum. Ideally, the amplification efficiency (E value) for each set of primers should be 100% (or as close to it as possible), implying that there was perfect doubling of the PCR product after every cycle during PCR amplification (Udvardi et al. [31] ). The E value of the 16S r RNA and gapdh genes was found to be less than 80% upon testing; then, further analysis was rejected, and the rpoB, rpoC, ldh, recA, mtlR, gpp, cfa1 and asd2 were selected for the next analysis. In a perfect RT‐qPCR assay, the correlation coefficient (R value) should be greater than 0.9 (Thomas et al. [30] , Saikaly et al. [26] ), and melt‐curve analysis should show amplification of only a single‐specific product in the absence of primer and dimer or non‐specific product formation. In this study, the amplification efficiency of eight candidate genes was in the range from 90.1 to 98.6%, and the correlation coefficient of all the genes was above 0.996 (Table ). These results indicated that the eight genes were suitable for reliable RT‐qPCR analysis.

The gene stability measure M was calculated using the software geNorm. The M value is defined as the average pairwise variation of a certain gene with all other considered control genes, whereas the variation of this certain gene to another is determined as the standard deviation of the log2‐transformed expression ratios. The gene with the highest M value can be considered the least stable one expressed and is eliminated from further recalculation of new M values. The geNorm analysis observed that the stability ranking of reference genes was in the following order: rpoB > rpoC > recA > ldh > cfal > mtlR > asd2 > gpp, and the M value of all genes was below the default limit of M = 1.5. The choice of the reference gene(s) for normalisation is critical to ensure reliability and validity of the RT‐qPCR assay (Suleman and Somai [29] ), but the use of a single reference gene for calibration and standardisation will affect the accuracy of the results (Zhu et al. [36] , Valihrach and Demnerova [32] ). Schmid et al. ([28] ) considered that using two or more than two reference genes contributes to the correction system deviation in a given set of samples or experimental conditions. It is becoming practice to utilise more than one reference gene in order to obtain normalisation data that are more accurate and reliable. In this study, four candidate genes were suggested as the reference genes according to the results of geNorm analysis.

The four most suitable reference genes were equally identified by the BestKeeper and geNorm software. There was a different stability ranking only for the other four genes. From the analysis of these two programs, rpoB, rpoC, recA and ldh can be considered the best reference genes in L. plantarum R23 under SO2 stress conditions due to the following reasons: first, because they were found to be most stable reference gene by each program (the M or SD value was lower than for the others); second, of the eight candidate genes, the four reference genes had the highest expression in L. plantarum (E value above 95%); and third, melt‐curve analysis of the four genes showed amplification of only a single specific product. All these results suggested that the expression of rpoB, rpoC, recA and ldh was not affected by the experimental SO2 stress conditions used in this study and that these were good candidate genes for use as reference genes.

The malolactic enzyme (MLE) is the key enzyme for MLF and has been purified from several species of LAB including Lactobacillus spp., O. oeni and Lactococcus lactis (Bartowsky [2] ). The mle gene from several of these species has been sequenced. Several factors, including malic acid, pH, ethanol, temperature and SO2, are known to influence the MLE activity of LAB (Lonvaud‐Funel [19] , Miller et al. [22] ). To further validate the selected reference genes revealed in our study, the relative expression level of the genes mleR1 and mleP1 mRNA of L. plantarum R23 with different SO2 treatments using reference genes for normalisation were compared. The results showed that normalisation using the most stable reference gene(s) led to interpretation of the expression level of the genes of interest. Using the least stable reference gene, however, the led to misinterpretation of the expression level of the gene of interest (Figure ). These results further illustrate the suitability of the reference genes revealed in the present study.

This is the first detailed study to evaluate selected reference genes in L. plantarum R23 under SO2 stress conditions. By using two software program for data analysis, this study observed that the genes rpoB, rpoC, recA and ldh were the most stably expressed reference genes and thus represent the reference gene combination that should be used to obtain the most accurate results for different SO2 stress conditions in L. plantarum R23.

This study was financially supported by the Natural Science Foundation of Fujian Province, China (no. 2015J01101).

Figure S1. Agarose gel electrophoresis of RNA for Lactobacillus plantarum R23 generated under stress conditions in the culture of 0 mg/L SO2 (A1, A2 and A3), 60 mg/L SO2 (B1, B2 and B3) and 120 mg/L SO2 (C1, C2 and C3). M, marker.

Figure S2. Melting curves generated for eight candidate reference genes: (a) rpoB, (b) rpoC, (c) ldh, (d) recA, (e) mtlR, (f) gpp, (g) cfal and (h) asd2.

PHOTO (COLOR): Expression stability values of reference genes calculated using the geNorm software.

PHOTO (COLOR): Optimal number of reference genes generated for normalisation. Pairwise variation analysis between the normalisation factors NFn and NFn + 1 to determine the optimal number of reference genes for normalisation.

PHOTO (COLOR): Relative quantification of the expression of (a) mleR1 and (b) mleP1 using validated reference genes for normalisation under SO2 treatments of 0 (), 60 () and 120 () mg/L. The genes rpoB, rpoC, recA, ldh alone and the combination (rpoB + rpoC, rpoB + ldh) were the most stable reference genes, and gpp was the most unstable gene. The expression level of treatment with 0 mg/L SO2 was set to 1. Each value represented the mean of three replicates.