Identification and selection of reference genes for gene expression analysis by quantitative real-time PCR in Suaeda glauca’s response to salinity

Quantitative real-time polymerase chain reaction (qPCR) using a stable reference gene is widely used for gene expression research. Suaeda glauca L. is a succulent halophyte and medicinal plant that is extensively used for phytoremediation and extraction of medicinal compounds. It thrives under high-salt conditions, which promote the accumulation of high-value secondary metabolites. However, a suitable reference gene has not been identified for gene expression standardization in S. glauca under saline conditions. Here, 10 candidate reference genes, ACT7, ACT11, CCD1, TUA5, UPL1, PP2A, DREB1D, V-H+-ATPase, MPK6, and PHT4;5, were selected from S. glauca transcriptome data. Five statistical algorithms (ΔCq, geNorm, NormFinder, BestKeeper, and RefFinder) were applied to determine the expression stabilities of these genes in 72 samples at different salt concentrations in different tissues. PP2A and TUA5 were the most stable reference genes in different tissues and salt treatments, whereas DREB1D was the least stable. The two reference genes were sufficient to normalize gene expression across all sample sets. The suitability of identified reference genes was validated with MYB and AP2 in germinating seeds of S. glauca exposed to different NaCl concentrations. Our study provides a foundational framework for standardizing qPCR analyses, enabling accurate gene expression profiling in S. glauca.

Suaeda glauca L. is a succulent halophyte that belongs to the Chenopodiaceae family 1 . It is widely distributed throughout Asia, particularly in the coastal areas of China 2 . Suaeda glauca is highly resistant to salt stress and thrives even in the presence of salt contents exceeding 0.48%, although no salt glands or vesicles are present in the leaves 1 . This plant is widely used for phytoremediation of saline soils primarily because of its Na + hyperaccumulation and fast-growing properties as well as its ability to produce large amounts of biomass 3,4 . S. glauca also shows potential for phytoremediation of Cd, Pb, and Mn in contaminated soils or mine tailings 5 and is used as a traditional Chinese medicinal plant to treat fever and alleviate stagnation 6 . The abundance of secondary metabolites makes S. glauca a multi-purpose bio-resource for both ecological and health-related issues. Particularly, the ethyl acetate-soluble fraction of the methanol extract of S. glauca is reportedly hepatoprotective 7 . Additionally, the anti-inflammatory, anti-mutagenic, anti-oxidative, and other bioactive properties of the flavonoid gallic acid 8 lead to protective effects in the liver 9 . However, the content of gallic acid in S. glauca ranges between 0.62 and 0.89% 10 ; this low content of the active agent is the primary limitation of nearly all medicinal plants. Previous studies reported that in medicinal herbs, the contents of secondary metabolites including soluble phenolics, anthocyanins, and flavonoids were increased under salt stress 11,12 . Thus, salinity may be used to boost the synthesis of secondary metabolites of interest in salt-tolerant medicinal plants. However, the specific molecular regulatory mechanism is poorly understood. Furthermore, many reports have indicated that the synthesis and accumulation of plant secondary metabolites are significantly correlated with the expression level of functional genes in their biosynthesis pathways 13,14 . Recent research on increasing the production of active compounds in medicinal plants through recalibration and/or engineering of pivotal genes along the biosynthesis

Results
Primer specificity and amplification efficiency test of candidate reference genes. Standard PCR was used to verify the primer specificity in the amplification of all 13 candidate reference genes. Based on the results, UBC28 and TIM showed multiple bands using cDNA as the template (Supplementary File S2A), whereas EF1α showed multiple bands amplified from both cDNA and genomic DNA (Supplementary File S2). These three genes were therefore discarded from further analysis. The other primers yielded a single and clear band of the expected size and there was no primer-dimer formation (Supplementary File S2). Further Sanger sequencing confirmed that the sequences of the amplified fragments were consistent with that given on NCBI, indicating the specificity of the primer pairs. The presence of a clear and single peak in melting curve analysis further verified the specificity of each primer set (Supplementary File S2). The details of the gene names, accession number, primer sequence, amplification length and efficiency, and correlation coefficient are shown in Table 1. The qPCR efficiency of all 10 candidate reference genes ranged from 87% (MPK6) to 119% (ACT7). The determination coefficients (R 2 ) of the regression equation varied from 0.9763 for MPK6 to 0.9994 for UPL1. Expression profiling of candidate reference genes in different tissues under various salinity stresses. The cycle threshold (Cq) values generated from qPCR for all 10 candidate reference genes were obtained and analyzed, and their variation under each treatment is presented in Fig. 1.
The Cq values of these genes varied from 18.92 (TUA5, Fig. 1B,F) to 31.75 (DREB1D, Fig. 1A,F), suggesting that TUA5 and DREB1D were the highest and least expressed genes, respectively. Candidate reference genes in different treatments showed a somewhat similar expression pattern. However, their expression levels differed under different treatments. While selecting the ideal reference genes among different tissues, PP2A showed the highest expression level (mean Cq of 20.91), and DREB1D exhibited the lowest expression level (mean Cq of 30.67) following treatment with only 1/10 MS liquid medium (Fig. 1A). The mean Cq values for TUA5 and ACT11 were 19.20 and 27.04, and 19.12 and 26.01, for samples obtained from plants exposed to 100 and 300 mM NaCl solution, respectively (Fig. 1B,C). Interestingly, the expression levels of DREB1D and H + -ATPase were markedly increased by NaCl solution treatment ( Fig. 1A-C). ACT7 (mean Cq of 21.10) and TUA5 (mean Cq of 20.13) had the highest expression levels in the seed and seedling samples treated with different NaCl concentrations, respectively. However, DREB1D (mean Cq of 27.22) and PHT4;5 (mean Cq of 28.93) showed the lowest expression levels in the seed and seedling samples treated with different NaCl concentrations, respectively (Fig. 1D,E). TUA5, PP2A, and MPK6 showed relatively high expression levels with a relatively narrow range of ΔCq of 2.94, 3.32, and 4.00, respectively, in all samples, whereas DREB1D, V-H + -ATPase, ACT7, and PHT4;5 showed relatively high variation with ΔCq values of 9.73, 7.82, 7.49, and 7.25, respectively (Fig. 1F). geNorm analysis. In geNorm analysis, a lower M value indicates more stable gene expression and vice versa. All M values of the 10 candidate reference genes in each treatment were lower than 0.7 (Fig. 2), which is much lower than the default limit of 1.5, indicating high expression stability. However, for tissue experiments, the three most stable genes were TUA5, PP2A, and MPK6 with the lowest M values, whereas DREB1D was the most unstable gene ( Fig. 2A-C). In the seed group, TUA5 and PP2A were the two most stable genes at all different concentrations of NaCl treatments (Fig. 2D). In the seedlings group, PP2A and MPK6 were the most stable genes, followed by V-H + -ATPase, whereas ACT7 and DREB1D were the least stable genes (Fig. 2E). Finally, for all sample sets, TUA5 and PP2A showed the lowest M value, and ACT7 was the most unstable gene (Fig. 2F). To determine the optimal number of reference genes required for accurate experimental normalization, pairwise variations of the normalization factor (V n /V n+1 ) were analyzed using geNorm. The V 2 / 3 values of all different experimental groups were below the threshold value of 0.15, indicating that two reference genes are sufficient to normalize gene expression data (Fig. 3).
NormFinder analysis. NormFinder, which is an analysis program based on the variation estimation statistical approach, was used to evaluate the expression stability of the candidate reference genes. A lower average expression stability indicates more stable gene expression. The stability values of the candidate reference genes in each treatment are listed in Table 2. For tissues experiment under 1/10 MS liquid medium, PP2A and TUA5 were the most stable reference genes and DREB1D was the least stable gene, which agreed with the results of geNorm analysis. For tissues experiment under 100 mM NaCl, MPK6 and PP2A were the most stable reference genes, whereas ACT11 was the least stable gene. Similarly, the most sta- Table 1. Primer sequences and amplification characteristics of candidate reference genes of S. glauca used in this study. Amplicon length was obtained using cDNA as a template.
The ranking of gene stability across seedlings under NaCl treatment was: In all samples, TUA5, PP2A, and MPK6 were the most stable reference genes and DREB1D and PHT4;5 were the least ranked reference genes, which agrees with the results of geNorm analysis.      www.nature.com/scientificreports/ BestKeeper analysis. The SD and r value of candidate reference genes in each treatment are listed in Table 3. Similar to the results of geNorm and NormFinder analysis, the three most stable genes were MPK6, PP2A, and TUA5, whereas DREB1D was the least stable gene in tissues (1/10 MS liquid medium). In addition to PP2A, TUA5, and MPK6, V-H + -ATPase was regarded as relatively stable in tissue experiments, both in 100 and 300 mM NaCl treatments. However, the two least stable genes in tissues after 100 mM NaCl treatment were DREB1D and ACT7. DREB1D and PHT4;5 were the two least stable genes in 300 mM NaCl-treated tissues. For seed and seedling experiments, the top three most stable genes were again PP2A, TUA5, and MPK6, whereas DREB1D was the least stable gene for the two experiments. Overall, in the total samples, similar to the geNorm and NormFinder Average expression stability M Least stable genes Most stable genes Most stable genes Most stable genes Most stable genes Most stable genes Most stable genes    www.nature.com/scientificreports/ gradient), and all samples) and third in two experiments (tissues (100 mM NaCl), and seedlings (NaCl gradient)) (Fig. 4). MPK6 was followed by TUA5, consistently ranking third. However, the order of the least stable genes among all treatment samples was DREB1D > ACT7 > PHT4;5 (Fig. 4).

Validation of identified stable reference genes.
To validate the suitability of the two selected reference genes, PP2A and TUA5, and the least stable reference gene DREB1D, the salt-induced expression level of MYB and AP2 in germinating seeds was normalized via qPCR. MYB (Fig. 5A) and AP2 (Fig. 5B) showed similar expression levels when a single or a combination of reference genes (PP2A and TUA5) was used to normalize the expression. The expression levels of MYB and AP2 were upregulated in all seed samples treated with different concentrations of NaCl. However, when DREB1D (unstable gene) was used for normalization, the relative expression patterns of MYB and AP2 differed as compared to the relative expression values obtained using the two most stable reference genes (PP2A and TUA5) (Fig. 5), highlighting the significance of selecting appropriate internal reference genes under salt stress.

Discussion
Gene expression analysis is a primary method for revealing the mechanism of plants' responses to different conditions, including abiotic stresses 23 . qPCR is an efficient and reliable method for assessing gene expression levels. However, the accuracy of the results is affected by the expression level and stability of internal reference genes. Using inappropriate reference genes can lead to misleading results 42 . Hence, application of steadily expressed internal reference genes is key to obtaining accurate results 20,43 . Commonly used methods for selecting candidate internal reference genes include searching the literature or using housekeeping genes, such as genes involved in cell metabolism, glycolysis, protein degradation, and synthesis 44 . However, previous studies reported that the expression level of internal reference genes obtained by this method may be unstable under different experimental conditions 45 , greatly affecting the accuracy of the results. Compared with the validated reference genes, the expression level of target genes can show 100-fold deviation using inadequately validated reference genes. Thus, using reference genes obtained by traditional methods may lead to erroneous or even contradictory results 21,43 . Therefore, constant expression in different tissues or experimental treatments is a vital feature of an ideal reference gene 46 . Generally, candidate reference genes belong to different types of functional genes, which can significantly reduce the risk of gene expression being affected. S. glauca, a succulent halophyte 1 , is used as a traditional Chinese medicinal herb because of its health-promoting pharmacodynamic effects 6 . In a previous study, ACTIN was used as the reference gene to normalize the salt response gene under salt stress 2 . Similarly, EF1α was used as a reference gene in another Suaeda species (S. salsa) to normalize the gene expression level under salt stress 47 . However, there are no reports on the systematic screening and identification of internal reference genes of S. glauca under salt stress. In this study, after a series of systematic screening steps, 10 reference genes were selected for further screening and validation as potential gene expression normalization genes under salt stress in S. glauca. Further assessment with ΔCq, geNorm, NormFinder, BestKeeper, and RefFinder software identified PP2A and TUA5 as the two most stable reference genes, whereas DREB1D was the most unstable gene in all samples. Notably, these reference genes belong to different functional gene classes. A previous study reported that GAPDH, ACTB, and 18S and 28S rRNA were the most frequently used reference www.nature.com/scientificreports/ genes in articles published by high-impact journals (more than 90% cases) 48 . However, in this study, none of these genes were identified as suitable reference genes for S. glauca. GAPDH was reported as the most appropriate reference gene under different salinity stress conditions for S. aralocaspica, which belongs to the same genus as S. glauca 24 . The constitutively expressed, housekeeping enzyme GAPDH converts glyceraldehyde-3-phosphate into 1,3-bisphosphoglycerate, which plays a key role in glycolysis. Hence, GAPDH is widely used as a reference gene 49 . Additionally, AvGAPDH-2 has been identified as the most stable reference gene in another traditional Chinese herb, A. venetum, under salinity stress 26 . Previous studies reported that GAPDH was one of the best reference genes in different tissue and organ samples of teak 50 and Dendrocalamus latiflorus 34 because of its high stable expression. However, GAPDH was not a suitable reference gene in rice under water-deficient conditions because of the high expression instability 51 . These results indicate that reference genes effective in one species may not be appliable in other species 34 . Although orthologous genes of DREB1D, PHT4;5, and ACT7 have been effectively used in studies of Fagopyrum tataricum 52 , A. hypogaea 29 , and A. venetum 26 , our comprehensive analysis showed that these were the top three most unstable genes in S. glauca. In contrast, PP2A and TUA5 were identified as the two most stable reference genes in the response of S. glauca to salinity. These genes were also chosen as optimal reference genes in C. lanatus 28  RefFinder stability value  RefFinder stability value  RefFinder stability value   www.nature.com/scientificreports/ them under various experimental conditions 55 , which showed that commonly used reference genes exhibited universality in different plants. However, every gene chosen in this method must be verified before it can be used in a new plant. A universally applicable reference gene does not exist 45,56 .

A C T 7 A C T 1 1 C C D 1 T U A 5 U P L 1 P P 2 A D R E B 1 D V -H + -A T P a s e
Many previous studies reported that a single reference gene was not adequate for normalization. The expression results would generate 3-6.4-fold error using only one internal reference gene 21 . To obtain accurate results, two or more stable reference genes are considered as necessary 26,57 . The geNorm program was used to determine the optimal number of internal reference genes for normalization with the pairwise variation parameters 21 . The result of pairwise variation analysis of all samples in this study revealed that the V2/3 value was below the threshold of 0.15 21 . Therefore, the combination of the most stable two internal reference genes would guarantee an accurate result. Transcription factors MYB 40 and AP2 41 are involved in the salt stress response. A previous study reported that the most responsive transcription factor families were AP2/ERF and MYB of Lotus japonicus treated with NaCl 58 . Hence, these salt-response genes were used to validate the suitability of the candidate reference genes in germinating seeds exposed to different NaCl concentrations. As shown in Fig. 5, both MYB and AP2 were upregulated in all seed samples treated with NaCl when a single or a combination of reference genes (PP2A and TUA5) was used. The results are consistent with the findings that MYB and AP2 are induced by salt stress. However, compared to using a single reference gene, the expression profiles of AP2 and MYB were more accurate when a combination of reference genes (PP2A and TUA5) was used to normalize the expression (Fig. 5). Overall, our results revealed the most relevant reference genes for normalization of relative gene expression of S. glauca tissues and response to salinity stresses.

Methods
Plant material and salinity stress treatments. Suaeda glauca seeds were collected from a saline field in Dongying City, Shandong Province, China. All experimental procedures including investigation and collection were in accordance with local and national regulations. The voucher specimen of S. glauca has been deposited in Naturalis Biodiversity Center (herbarium ID L.1678537_544178685). Professor Jie Song (College of Life Science, Shandong Normal University, China) undertook the formal identification of S.glauca used in your study. We were allowed to collect the plant by the local government. Seeds of S.glauca were grown in a greenhouse with a temperature range of 25 °C ± 1 °C and 16-/8-h light/dark photoperiod cycles. The plants were grown in vermiculite and watered using 1/10 Murashige and Skoog (MS) liquid medium. Considering that S. glauca is autochthonous of saline landscapes and that the salt content and types of salts present in nature may differ from those of artificial setup, three groups of experiments were designed for tissue collection. A total of 60 four-week-old seedlings were sorted into three groups of 20 plants. One group was irrigated with 1/10 MS liquid medium, whereas the other two were irrigated with either 100 or 300 mM NaCl dissolved in 1/10 MS liquid medium. For each group, the following samples (roots, stems, leaves, and shoots) collected at different growth stages (seedlings, adult plants, and senescent plants), and followers were used for expression analyses of candidate reference genes screening. Samples for each tissue were collected from at least five individual plants and quickly frozen in liquid nitrogen and stored at − 80 °C throughout the sampling period. A series of NaCl solutions (0, 50, 100, 150, and 200 mM NaCl dissolved in 1/10 MS liquid medium) was used to treat the germinating seeds. Seed samples were collected at 1, 2, and 4 days after treatment. Furthermore, two-week-old seedlings were irrigated with a series of NaCl solutions (0, 50, 100, 150, and 200, and 400 mM NaCl dissolved in   36 . Before qPCR, the PCR products of all primer pairs were electrophoresed on a 2% agarose gel stained with ethidium bromide to verify successful amplification and primer specificity. Primer specificity was validated via melting curve analyses following amplification using qPCR. qPCR conditions and amplification efficiency test. The 7500 Real-Time PCR system (Applied Biosystems, Foster City, CA, USA) was used to perform qPCR under the following conditions: 50 °C for 2 min and initial denaturation at 95 °C for 5 min, followed by 40 cycles at 95 °C for 10 s and 60 °C for 35 s. The dissociation curve was analyzed to determine the specificity of primer amplification, with the dissociation temperature ranging from 60 to 95 °C, by increasing the temperature stepwise by 1 °C/min 26 . The volume of the qPCR system was 20 μL. The reaction mixture consisted of 10 μL 2 × SYBR green mix (Q321, Vazyme), 0.4 μL ROX Reference Dye, 1 μL cDNA, and 0.4 μL of each forward and reverse primer (10 mM). Each sample was evaluated with three biological replicates. To identify the specific PCR amplification efficiency of the reference gene primer pairs, a tenfold dilution series (10-1000-fold dilution) of a cDNA isovolumetric mixing pool was used as templates for qPCR. The primer amplification efficiency (E) and correlation coefficient (R 2 ) were calculated from the raw Cq values using the equation: E = 10 -1/slope -1 28 , where the slope is derived from the regression equation, calculated using Excel linear regression parameter (Microsoft Corp., Redmond, WA, USA).
Stability assessment of candidate genes. The statistical algorithm software programs geNorm 21 , NormFinder 22 , BestKeeper 37 , and RefFinder 38 were used to analyze the expression stability of candidate reference genes. The raw Cq value was used as input data for further analysis with RefFinder and BestKeeper. 2 −ΔCq value and 2 Cq value were used as input data for further analysis with geNorm and NormFinder, respectively. The parameters calculated using geNorm were the average gene expression stability measure (M) and average pairwise variation (V). A lower M value of the reference gene indicated more stable expression 21 . The parameter calculated by NormFinder was the stability value, which is related to the systematic error of each candidate gene. Smaller values indicate more stable gene expression 22 . The parameters calculated using BestKeeper were the SD values of Cq and coefficient of correlation (r) value. A gene with an SD value below 1.0 and r value close to one indicate more stable gene expression 37 . RefFinder generated a final comprehensive ranking of the most consistent and appropriate candidate reference genes according to the weighted geometric mean of the results of geNorm, NormFinder, and BestKeeper 38,39 . The lowest geometric mean indicated the most stable gene.
Validation of reference genes. To validate the reliability of the selected reference genes, the relative expression levels of two transcription factors MYB 40 and AP2 41 , which are involved in the salt stress response, were used as stress indicator genes to validate the suitability of candidate reference genes using qPCR at germinating seed samples.

Conclusion
We systematically evaluated the expression stability of different potential reference genes for qPCR in S. glauca. Thirteen typical reference genes were selected to identify the most stable reference genes for qPCR normalization in S. glauca. Based on the primer specificity results derived from the PCR results, UBC28, TIM, and EF1α were discarded because multiple bands were amplified from the cDNA and gDNA templates. ACT7, ACT11, CCD1, TUA5, UPL1, PP2A, DREB1D, V-H + -ATPase, MPK6, and PHT4;5 were further analyzed in different tissues and under various salt concentration treatments. PP2A and TUA5 were identified as the best reference