Selection and Validation of Reference Genes for RT-qPCR Analysis in Spinacia oleracea under Abiotic Stress

Reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) is an accurate and convenient method for mRNA quantification. Selection of optimal reference gene(s) is an important step in RT-qPCR experiments. However, the stability of housekeeping genes in spinach (Spinacia oleracea) under various abiotic stresses is unclear. Evaluating the stability of candidate genes and determining the optimal gene(s) for normalization of gene expression in spinach are necessary to investigate the gene expression patterns during development and stress response. In this study, ten housekeeping genes, 18S ribosomal RNA (18S rRNA), actin, ADP ribosylation factor (ARF), cytochrome c oxidase subunit 5C (COX), cyclophilin (CYP), elongation factor 1-alpha (EF1α), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), histone H3 (H3), 50S ribosomal protein L2 (RPL2), and tubulin alpha chain (TUBα) from spinach, were selected as candidates in roots, stems, leaves, flowers, and seedlings in response to high temperature, CdCl2, NaCl, NaHCO3, and Na2CO3 stresses. The expression of these genes was quantified by RT-qPCR and evaluated by NormFinder, BestKeeper, and geNorm. 18S rRNA, actin, ARF, COX, CYP, EF1α, GAPDH, H3, and RPL2 were detected as optimal reference genes for gene expression analysis of different organs and stress responses. The results were further confirmed by the expression pattern normalized with different reference genes of two heat-responsive genes. Here, we optimized the detection method of the gene expression pattern in spinach. Our results provide the optimal candidate reference genes which were crucial for RT-qPCR analysis.


Introduction
Reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) is an accurate and convenient method of quantifying mRNA levels of gene expression [1]. Selection of appropriate reference genes is crucial for validating accurate gene expression [2,3]. Improper reference genes used in data processing may lead to inaccurate and even wrong results [4]. The commonly used reference genes in RT-qPCR analysis are the housekeeping genes because they are usually expressed steadily at mRNA levels in any organs under various conditions [5,6]. Although 18S ribosomal RNA (18S rRNA), actin, and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) are usually taken as reference genes [5], their mRNAs are not always stable in any cases [7][8][9][10]. Therefore, selection and optimization of reference genes are important steps in RT-qPCR experiments.
Some statistical algorithms, such as BestKeeper, NormFinder, and geNorm, are widely used for analyses of reference genes for RT-qPCR. BestKeeper analyzes the stability by calculating the standard deviation (SD) of the quantification cycle (Cq) values of candidate reference genes [11]. NormFinder compares the variation within and between sample groups of candidate genes and calculates the stability value of each gene based on the 2 −ΔCq of genes [12]. geNorm evaluates the stability of candidate genes through calculating the stability value based on the geometric mean of 2 −ΔCq of genes and mean pairwise variation in sample groups, as well as providing the optimal numbers of reference genes under each condition [13].
Spinach (Spinacia oleracea) belongs to the Amaranth family and is rich in carotenoid, vitamins, and minerals. It is a favorite vegetable all over the world. In our previous study, the draft of the spinach genome has been sequenced, with 25,495 encoding genes predicted [14]. Moreover, the patterns of gene expression and protein abundance have been reported in spinach in response to diverse stresses (e.g., heat, salinity, heavy metal, and virus) using molecular genetic, transcriptomic, and proteomic approaches [15][16][17][18][19][20][21][22][23].
To date, there are only a few reports on gene characterization and function analyses in spinach in response to stresses. A recent study reported that the S. oleracea heat shock 70 (SoHSC70) was induced by heat stress, and overexpression of SoHSC70 enhanced the heat tolerance in spinach calli [24]. Besides, spinach cytochrome P450 85A1 (SoCYP85A1) was upregulated in response to Phytophthora nicotianae infection, and 35S-promoted SoCYP85A1 overexpression conferred resistance to P. nicotianae pathogen inoculation in tobacco [25]. In addition, the spinach nonsymbiotic hemoglobin family gene (SoHb) was induced by several stress treatments (i.e., polyethylene glycol, NaCl, H 2 O 2 , salicylic acid, and nitric oxide) but suppressed by a nitric oxide scavenger, nitrate reductase inhibitor, and nitric oxide synthase inhibitor. Overexpression of SoHb in Arabidopsis resulted in the decreases in nitric oxide level and sensitivity to nitrate stress [26]. In these studies, some housekeeping genes in spinach were used as reference genes in RT-qPCR experiments for normalization analysis. 18S rRNA was used as a reference gene to detect the expression patterns of chilling-/drought-responsive SoCAP85 (85 kD cold acclimation protein) [27], drought-/salt-/oxidative stress-responsive SoHb [26], 13 heat-responsive genes (including SoHSFB2b and Sob-ZIP9) [22], 15 nitrate transport and assimilation-related genes [28], and anthocyanin biosynthesis-related genes in various spinach germplasms [29]. In addition, actin, GAPDH, and ubiquitin 5 (UBQ5) were used to normalize the expression levels of several genes in response to various stresses, such as drought [30], biotic stress [25,31], gibberellic acid (GA) treatment, and gender-specific condition [32]. However, only one study evaluated the stability of these reference genes [27]. Five commonly used housekeeping genes (i.e., GAPDH, actin, 16S rRNA, tubulin alpha chain (TUBα), and 18S rRNA) were evaluated between partially/fully hydrated versus dry seeds of spinach under chilling, desiccation, and optimum conditions. Among them, 18S rRNA appeared to be most stable, but still fluctuated under several treatments [27].
In this study, ten candidate genes (i.e., 18S rRNA, actin, ARF, COX, CYP, EF1α, GAPDH, H3, RPL2, and TUBα) were selected for reference genes in spinach. Three Excel programs (i.e., BestKeeper, geNorm, and NormFinder) were used to evaluate the stability of these candidate genes in different organs, as well as stresses of heat, heavy metal, NaCl, Na 2 CO 3 , and NaHCO 3 . Optimal reference genes for each condition were verified. In addition, the two stable reference genes, ARF and actin, and the commonly used reference gene TUBα were selected to normalize the mRNA levels of two representative heat-responsive genes (SobZIP9 and SoHSFB2b) on references.

Plant Materials and Growth Condition.
A heat-resistant sibling inbred line of spinach, Sp75, was used in this study. Spinach plants were placed in a growth chamber with a temperature regime of 22/18°C, 10/14 h day/night cycle, and a relative humidity of 60%. The top third and fourth leaves, stems, roots, male flowers, and female flowers of plants with uniform growth were sampled at 50 days after planting. Seedlings were sampled at 10 days after planting. These samples were flash-frozen in liquid nitrogen and stored at -80°C for further experiment. Three biological replicates were taken for each organ, and at least three plants were used for each replicate.

Stress Treatment.
Fifty-day-old spinach seedlings with uniform growth were used for stress treatment. For heat treatments, the plants were moved into a chamber (37°C) at 0, 1, 2, 4, 6, 8, 12, 24, and 48 h before they were sampled [15]. For the treatments with heavy metal, salt, and alkali, the seedlings were watered with 200 μmol/L CdCl 2 [51], 200 mmol/L NaCl [52], 200 mmol/L NaHCO 3 [53], and 100 mmol/L Na 2 CO 3 [54,55], and the treatment times are 0, 1, 3, 6, 12, 24, and 48 h [56]. The top third and fourth leaves as well as roots were sampled for RNA extraction. Three biological replicates were taken for each time point of all stress treatments, and at least three plants were used for each replicate.
2.3. RNA Isolation and First-Strand cDNA Synthesis. Plant samples were ground in liquid nitrogen with a mortar and a pestle. Total RNA was isolated from 100 mg sample powder with TRIzol™ LS Reagent (Invitrogen, USA). The RNA samples with 260/280 ratios ranging from 1.8 to 2.1 were used for the following experiment. Total RNA was also examined by electrophoresis with 1% agarose gel to ensure the integrity. One microgram of total RNA was used for first-strand cDNA synthesis in a 20 μL total volume with a mixture of oligo dT primer and Random 6-mer in PrimeScript™ RT reagent (Takara, Japan) according to the manufacturer's instructions.

Primer Design and Evaluation.
The primer pairs of each reference gene were designed according to their sequences by using the online program Primer3Plus (Table 1) (http://www .primer3plus.com/cgi-bin/dev/primer3plus.cgi) [60]. The mixed cDNA from all samples was used as a template in primer evaluation. The PCR amplification products of each primer pair were checked by 2% agarose gel electrophoresis and sequencing, and then, the specificity of each pair of primers was evaluated by melting curve analysis followed by the amplification in RT-qPCR. Standard curves of each primer pair were established using a 5-fold dilution series ([1/1], [1/5] [61]. The amplification efficiencies (E) of these primer pairs were calculated by the slope of standard curves (E = 10 − 1/slope), and the correlation coefficients (R 2 ) were acquired from the standard curves as well [62].
2.6. Real-Time Quantitative PCR. RT-qPCR analysis was performed in 0.2 mL tubes with the Applied Biosystems 7500 Real-Time PCR System (ABI, USA). Each reaction contained 1 μL cDNA (5-time diluted), 10 μL AceQ Universal SYBR qPCR Master Mix (Vazyme, China), 0.5 μL of primer (500 μmol/L), and 8 μL deionized water. The PCR was carried out as the following program: predegeneration at 95°C for 3 min; 40 cycles of degeneration at 95°C for 15 s, annealing at 55°C for 15 s, and extension at 72°C for 30 s; and melting curve analysis at 65°C-95°C. RT-qPCR of each cDNA sample was carried out three times as technical replicates.

Stability Evaluation of Candidate Reference Genes.
The Cq values of 10 candidate genes were put in an Excel sheet, and a boxplot of these Cq values was generated. The stability of these 10 genes was evaluated by three statistical Excel macro programs, including BestKeeper [11], NormFinder [12], and geNorm [13]. The Cq values of each gene in roots, stems, leaves, flowers, and seedlings as well as leaves and roots in response to various stresses (i.e., high temperature, NaCl, Na 2 CO 3 , NaHCO 3 , and CdCl 2 ) were used in this evaluation. In BestKeeper, the standard deviation (SD) of the Cq values of each candidate gene was calculated, and the gene with the lowest SD was taken as the most stable gene [11,63]. In NormFinder, the expression stability (M1) was calculated by Cq values obtained by RT-qPCR of candidate genes and ranked in each sample set. The lowest M1 indicates that the gene is most stable [12]. In geNorm, the stability of candidate genes was evaluated by relative expression levels (Q) transformed from the Cq values for each sample according to the formula of Q = 2 −ΔCq ðΔCq = Cq value of each sample − the minimum Cq value in each setÞ. This formula works under the precondition that the efficiency of primers should range from 90% to 105% [64,65]. An average expression stability value (M2) of each candidate was calculated to demonstrate their stability. The gene with the lowest M2 value was regarded as the most stable expression. Besides, the geNorm software also determines the optimal number of reference genes required for RT-qPCR data normalization under each condition by pairwise variation (V n /V n+1 ) between the normalization factors NF n and NF n+1 . If V n / V n+1 < 0:15, the first n is the optimal number of genes required for this condition [13].
2.8. Normalization of SobZIP9 and SoHSFB2b. The expression pattern of two heat-responsive genes SobZIP9 and SoHSFB2b [22] in spinach was detected in response to heat stress (37°C). The mRNA levels of SobZIP9 and SoHSFB2b were then normalized by the most stable genes suggested by NormFinder and BestKeeper, as well as TUBα, a commonly used reference gene. Primer pairs for SobZIP9 (qSoHSFB2b-F: TCTTTCCACACTCGCTCTGT, qSoHSFB2b-R: CGGA TTACAAGAAGGCAGGC) and SoHSFB2b (qSobZIP9-F: TGCTGGAAACCCTAGGACTG, qSobZIP9-R: CTTCTG GTGCTTCTAGGCCT) [22] were used in this experiment. A linear ANOVA was used for evaluation of the variation of each gene in response to heat stress.

Results
3.1. Analysis of Primer Specificity and PCR Amplification. Ten candidates were chosen as reference genes for RT-qPCR analysis. They are 18S rRNA, actin, ARF, COX, CYP, EF1α, GAPDH, H3, RPL2, and TUBα. To evaluate the specificity of designed primers for the ten genes, the analyses of PCR, gel electrophoresis, and melting curves were performed. The gel electrophoresis showed a single band with expected size of each pair of primers (Figure 1), and the melting curves of each primer pair exhibited a single peak (Supplementary Figure S1), indicating the specificity of these primer pairs of candidate genes. The target amplicons were sequenced, and the results were consistent with their gene sequences in SpinachBase [14,30,58]. The standard curves indicated that the RT-qPCR amplification efficiency of candidate genes ranged from 92.53% (ARF) to 102.80% (EF1α), and the correlation coefficients varied from 0.991 (COX) to 0.999 (actin, CYP, and GAPDH) ( Table 1 and Supplementary Figure S2). Thus, these primers are specific for their respective genes and can be used in RT-qPCR analysis.

Expression Profiles of Ten Candidate Genes in Spinach.
Analysis of the expression levels of ten candidates was performed in young seedlings, roots, stems, leaves, male flowers, and female flowers, as well as leaves and roots in response to stresses of heat, heavy metal, NaCl, Na 2 CO 3 , and NaHCO 3 . Cq values of ten candidates in various samples obtained by RT-qPCR were shown in a boxplot ( Figure 2). The Cq values varied from 10.19 (18S rRNA in organs/seedlings) to 35.79 (TUBα in heat treatment) indicating that these candidate genes present different expression levels. The Cq value range reveals variability among the candidate genes. 18S rRNA shows the minimal range of Cq values in organs/seedlings (2.08, Figure 1(a)) and under NaCl treatment (0.86, Figure 1(d)), while actin (3.07, Figure 1(c)), ARF (1.64, Figure 1(c)), GAPDH (3.09, Figure 1(e)), and CYP (1.38, Figure 1(f)) show the minimal range of Cq values under heat, heavy metal (CdCl 2 ), NaHCO 3 , and Na 2 CO 3 treatments, respectively. These minimal ranges of Cq values indicated that these genes are more stable than others in each condition. However, further analyses are needed, because the comparison of the range of Cq values along is deficient to reveal the stability of the candidate genes. Thus, statistical macro programs were then used in this study.

Expression Stability of Ten Candidate
Genes. Three statistical Excel macro programs (i.e., BestKeeper, NormFinder, and geNorm) were used to evaluate the stability of ten candidate genes, in order to find optimal reference genes in spinach   BioMed Research International for RT-qPCR normalization ( Table 2, Figures 3 and 4). The stability of candidate genes was evaluated in 6 data sets, organs/seedlings, leaves, and roots under heat, heavy metal, NaCl, NaHCO 3 , and Na 2 CO 3 , respectively. The stability of candidate genes was analyzed with Best-Keeper (Table 2), which can calculate the standard deviation (SD) on the basis of the Cq values of all candidate reference genes [11]. The reference genes exhibiting the lowest standard deviation (SD) were taken as the most stable genes [63]. In the samples of different organs, 18S rRNA    (Table 2). Importantly, the TUBα was identified as the least stable gene in organs (SD = 1:58), as well as under treatments of heat (SD = 1:67), NaCl (SD = 1:24), and NaHCO 3 (SD = 2:13).
The expression stability of ten genes was also analyzed using the NormFinder software (Table 2), which can provide a stability value (M1) for each gene by comparing the variation within and between sample groups. For the organs/seedlings, there were six sample groups, including young seedlings, roots, stems, leaves, male flowers, and female flowers. For the heat-treated samples of leaves and roots, there were nine sample groups, which were samples under heat (37°C) stresses for 0, 1, 2, 4, 6, 8, 12, 24, and 48 h, respec-tively. In addition, for the samples under heavy metal, NaCl, Na 2 CO 3 , and NaHCO 3 treatment, there were six sample groups for each treatment (i.e., 0, 1, 3, 6, 12, 24, and 48 h) and totally 24 sample groups. The variation of each candidate gene in each case was compared, respectively. The gene stability was evaluated by M1 values.
The gene with lower M1 indicated that it was more stable [12].     Figure 3: Average expression stability values (M2) in different data sets obtained from the software geNorm. The stability values (M2) and stability rank were obtained from geNorm. A lower M2 value suggests higher stability: (a) stability rank obtained in organs including leaves, stems, roots, flowers, and seedlings, (b) stability rank obtained in the leaves and roots under heat stress, (c) stability rank obtained in the leaves and roots under heavy metal stress, (d) stability rank obtained in the leaves and roots under NaCl, (e) stability rank obtained in the leaves and roots under NaHCO 3 , and (f) stability rank obtained in the leaves and roots under Na 2 CO 3 . Ten candidate reference genes refer to Figure 1. Three biological replicates were taken for different organs/seedlings, as well as leaves and roots at each time point of all stress treatments. 7 BioMed Research International and GAPDH (M1 = 0:523) were the least stable genes under NaHCO 3 and Na 2 CO 3 treatments, respectively (Table 2).
In addition, geNorm was used for gene stability analysis. geNorm evaluates the stability of candidate genes based on the geometric mean of these genes and mean pairwise variation in sample groups [13]. As mentioned above, there were six, nine, and 24 sample groups for organs/seedlings, heat stress, and other stresses, respectively. The stability value (M2) ( Figure 3) and pairwise variation (V n /V n+1 ) (Figure 4) in the results given by geNorm revealed the gene stability and optimal gene numbers in certain case. The genes with the lowest M2 values were regarded as the most stable ones in each case ( Figure 3). Besides, when pairwise variation (V n /V n+1 ) was less than 0.15, the minimum value of n was the optimal number of genes required for such condition (Figure 4). And the suitable gene pair included the genes from rank 1 to rank n in each condition [13]. According to this validation method, we recommended several reference genes for RT-qPCR analyses in organs and stress response. In different organs, EF1α/RPL2 (V 2 /V 3 = 0:132) was a suitable gene pair for mRNA level normalization (Figures 3(a) and 4). EF1α, RPL2, and COX (V 3 /V 4 = 0:125) were identified as an appropriate gene set under heat treatment (Figures 3(b) and 4). Besides, the gene pair of 18S rRNA/actin (V 2 /V 3 = 0:116) was recommended for heavy metal stress (200 μM CdCl 2 ) (Figures 3(c) and 4), while EF1α/H3/ARF (V 3 /V 4 = 0:106) were suggested as reference genes under NaCl (200 mM) (Figures 3(d) and 4) treatment. In addition, CYP/H3/RPL2/actin (V 4 /V 5 = 0:131) and ARF/H3/RPL2/EF1α (V 4 /V 5 = 0:129) were selected as reference gene sets under 200 mM NaHCO 3 (Figures 3(e) and 4) and 100 mM Na 2 CO 3 (Figures 3(f) and 4), respectively.

Validation Test of Candidate Genes.
To prove the feasibility of reference genes for RT-qPCR in spinach, mRNA levels of two heat-responsive genes, spinach basic regionleucine zipper 9 (SobZIP9) and heat stress transcription factor 2b (SoHSFB2b) [22], were detected in spinach under heat treatment and normalized by heat-stable actin and ARF, as well as commonly used TUBα. Under heat stress, the expression level of SobZIP9 was reduced in spinach after 2-12 h heat treatment when normalized by actin and also reduced after 2 h heat treatment when normalized by ARF. However, when normalized by TUBα, it exhibited an increase under heat treatment ( Figure 5(a)). Besides, SoHSFB2b was increased about 10-fold after 1 h heat stress normalized by actin and ARF, but it exhibited 24-time increase when being normalized by TUBα ( Figure 5(b)). This implies that the SoHSFB2b expression level was overestimated when normalized by TUBα. In addition, when normalized by actin and ARF, the levels of heat-reduced SobZIP9 and heat-induced SoHSFB2b were consistent with those in spinach under 35°C for 30 min and 5 h based on results from a previous transcriptomic study [22], which showed that actin and ARF are the appropriate reference genes in spinach for heat response analysis.

Discussion
Diverse PCR approaches have been applied in the evaluation of gene expression levels, such as semiquantitative PCR, RT-qPCR, and digital PCR (dPCR). About ten years ago, digital PCR (dPCR) was developed as a novel technology for mRNA quantitation at a single molecular level [66,67]. The dPCR can detect absolute copy numbers of certain gene expression without reference genes and standard curve [68], and it is more sensitive than qPCR for detecting the smaller copy number variation of genes [69]. However, dPCR is labourintensive and expensive and has relatively low throughput  Figure 4: Pairwise variation (V) analysis of the candidate reference genes. The geNorm software was used to analyze the pairwise variation (V n /V n+1 ) between the normalization factors (NF) NF n and NF n+1 in order to determine the optimal number of candidate reference genes required for RT-qPCR data normalization. If V n /V n+1 < 0:15 (gray dotted line), the minimum value of n is the optimal number of genes required. 8 BioMed Research International when compared with qPCR [67], which leads to being not popularly used for the evaluation of gene expression in plant molecular labs. RT-qPCR is a convenient and accurate method to detect the mRNA levels of certain genes. The accuracy requires one or more stable reference genes for calculating the gene expression levels using the 2 −ΔΔCq method [64]. Transcriptomic results suggest that there is always more than one gene stable in each set of samples, but none of them remains stable under all the conditions [5]. To date, 18S rRNA [22,28,29], actin [25], GAPDH [27,31], and ubiquitin [32] were used as reference genes in previous studies in spinach, although their stability was not validated yet. One study in spinach seeds evaluated the stabilities of GAPDH, actin, 16s rRNA, TUBα, and 18S rRNA, and 18S rRNA was taken as a relative stable reference gene [27]. In this report, we have determined 18S rRNA as the optimal reference gene in roots, stems, leaves, flowers, and seedlings, as well as in leaves and roots in response to salinity stress.
TUB (TUBα or TUBβ) was commonly used as a reference gene in some plant species, such as Raphanus sativus, Baphicacanthus cusia, and Hibiscus cannabinus [33,34,57]. However, in this report, TUBα was regarded as the most unstable reference gene in spinach in many cases. This is similar to the situation in Achyranthes bidentata and Amaranthaceae. In A. bidentata under hormone treatments (e.g., indole-3-butytric acid and methyl jasmonate), NaCl, and drought, TUB was also determined to be an unstable reference gene when compared with other candidate references (e.g., 18S rRNA, actin, APT1, EF1α, GAPDH, TUBβ, UBC, and UBQ) [35]. Besides, EF1α was regarded as the optimal reference gene in spinach under NaCl treatment. Similarly, EF1α was also the optimum gene in A. bidentate under NaCl and drought stresses [35]. It should be noted that EF1α was considered relatively unstable in several other species, such as Silene vulgaris and Papaver rhoeas [37,57]. This implies that homologous genes in different plant species may not always exhibit similar expression levels under different conditions.
The recommended reference genes from NormFinder and BestKeeper look different ( Table 2), because of their different mathematical models [11,12]. BestKeeper uses Cq to calculate SD of candidate genes [11], but NormFinder uses 2 −ΔCq to calculate the expression stability [12]. In our results, the expression patterns of heat-responsive genes (SobZIP9 or SoHSFB2b) are consistent when using the SoARF and SoActin as reference genes, which are recommended from Norm-Finder and BestKeeper, respectively ( Figure 5). Importantly, the heat-responsive patterns of SobZIP9 or SoHSFB2b normalized by SoARF and SoActin are also similar with previous reports [22], but these heat-responsive patterns are absolute opposite or exaggerated when using TUBα as a reference gene that is not recommended by BestKeeper and NormFinder. This indicates that the optimal genes obtained in response to heat stress by NormFinder and BestKeeper are reliable for normalization of gene expression.
Due to different expression levels of genes, the Cq values of some reference genes are high (Cq > 30, Figure 2), or the ΔCq between the reference gene and target genes is relatively large. In this condition, even if the amplification efficiency of both primers of the reference gene and target genes is close to 100%, the error of gene expression level is probably enlarged, due to using the 2 −ΔΔCq method [64,70]. Thus, in this case, it is better to normalize the mRNA levels using the values of primer amplification efficiency rather than 100%.
In this study, optimal reference genes of each sample set were determined in spinach. 18S rRNA and ARF were validated as internal reference genes in different organs. Actin and ARF, instead of 18S rRNA, were the most suitable genes in leaves and roots under heat treatment. For reference genes under various abiotic stresses, EF1α and ARF were suitable for CdCl 2 , 18S rRNA and COX for NaCl, RPL2 and actin   Figure 5: Expression pattern of SobZIP9 and SoHSFB2b normalized by different genes in response to heat stress. Expression patterns of two heat-induced genes, basic region-leucine zipper 9 (bZIP9, a) and heat stress transcription factor 2 b (HSFB2b, b), were normalized with actin, ADP ribosylation factor (ARF), and tubulin alpha chain (TUBα). Three biological replicates were taken for leaves at each time point of heat treatments. Data represent average ± SD, and * indicates significant difference (p < 0:05).
9 BioMed Research International for NaHCO 3 , and ARF and CYP for Na 2 CO 3 . These results demonstrate that different reference genes should be used under different conditions [9]. Taken together, actin should be a reference gene for evaluating gene expression across all organs under various stress conditions, because actin has relatively better rank and M/SD values among all cases ( Table 2).

Conclusion
It is important that optimal genes should be used for certain conditions when RT-qPCR is conducted to determine the normalized gene expression pattern using the 2 −ΔΔCq method. Commonly used housekeeping genes in plant species may not be suitable under all the conditions or in certain species. In this study, the optimal genes were determined for gene expression normalization in spinach organs and in response to stresses of heat, CdCl 2 , NaCl, NaHCO 3 , and Na 2 CO 3 . The data provide a list of useful reference genes for future studies of gene expression patterns in spinach.

Data Availability
The data used to support the findings of this study are included in this published article and its supplementary information files.