Selection and Validation of Reference Genes for Quantitative Real-Time PCR in White Clover (Trifolium repens L.) Involved in Five Abiotic Stresses

White clover (Trifolium repens L.) is a widely cultivated cool-season perennial forage legume in temperate grassland systems. Many studies have analyzed the gene expression in this grass species using quantitative real-time reverse transcription PCR (qRT-PCR). The selection of stable reference genes for qRT-PCR is crucial. However, there was no detailed study on reference genes in different tissues of white clover under various abiotic stress conditions. Herein, 14 candidate reference genes (ACT7, ACT101, TUA1109, TUB, CYP, 60SrRNA, UBQ, E3, GAPDH1, GAPDH2, PP2A, BAM3, SAMDC, and ABC) were selected and analyzed by four programs (GeNorm, NormFinder, BestKeeper, and RefFinder). Samples were taken from two tissues (leaves and roots) under five different abiotic stresses (drought, salt, heat, cold, and heavy metal stress). Our results showed that 60SrRNA and ACT101 were the two top-ranked genes for all samples. Under various experimental conditions, the most stable gene was different; however, SAMDC, UBQ, 60SrRNA, and ACT101 were always top ranked. The most suitable reference genes should be selected according to different plant tissues and growth conditions. Validation of these reference genes by expression analysis of Cyt-Cu/Zn SOD and CAT confirmed their reliability. Our study will benefit the subsequent research of gene function in this species.


Introduction
Quantitative real-time reverse transcription PCR (qRT-PCR) using cDNA as a template is a sensitive and powerful technique for measuring gene expression level [1]. Quantitative real-time PCR can be used not only to analyze the regulation of gene expression, monitor the expression pattern of mRNA, and quantitatively analyze the transcription level of genes, but also to conduct spatial-temporal analysis of target genes in different tissues under various treatments; as such, qRT-PCR has been widely applied in the field of molecular biology [2][3][4]. However, the quality of gene expression is affected by many experimental factors [5], and the relative quantitative method selected by researchers first needs to correct the expression amount of the target gene in the experimental process [6]. So far, using one or more reference genes is the preferred method of normalization [7], and this is also the simplest method of data processing. Accurate determination of target gene expression levels depends on the selection of stable reference genes [8][9][10].
White clover (Trifolium repens L.) is an important cool-season perennial legume forage which is widely cultivated in temperate grassland systems [11]. It has a high feed value and strong N fixation capacity in soil [12]. However, white clover is susceptible to drought, salt, and heat stress [13].

Plant Materials, Growth Conditions, and Abiotic Stress Treatments
Seeds of white clover cv. Haifa were grown in 20 × 15 × 5 cm pots containing 1 kg silica sand. Each pot was sprinkled with 0.4 g of seeds, supplemented with Hoagland's nutrient solution, and maintained in an environmental growth chamber set to a light intensity of 100 µmol/(m·s) at 23/19 • C (day/night) with a 12-h photoperiod. One-month-old plants were used for all stress experiments. We observed the extent of five different abiotic stresses on plants, characterizing stress as mild, moderate, and severe, and then took samples. For drought stress, the plants were treated with 17% concentration of PEG6000, and samples were taken at 0, 6, 8, and 10 days. For salinity stress, the plants were treated with 250 mmol/L NaCl, and samples were taken at 0, 12, 14, and 16 days. For heat stress, plants were put in an environmental growth chamber set to 38/33 • C (day/night), and samples were taken at 0, 15, 22, and 23 days. Cold stress was imposed at 4 • C in an incubator, and samples were taken at 0, 12, 17, and 22 days. For heavy metal treatment, the nutrient solution was filled with 600 µmol/L Cd 2+ Hoagland's, and samples were taken at 0, 12, 17, and 22 days. Meanwhile, control plants were treated with an equal quantity of Hoagland's nutrient solution. For all controls and treatments, the leaf and root tissues were sampled separately at three different time points with four biological replicates. All the samples were frozen immediately in liquid nitrogen and stored at −80 • C for later use.

RNA Isolation and cDNA Synthesis
Total RNA was extracted from the leaf and root tissues using the HiPure Universal RNA Kit (Magen Biotech Co., Ltd., China) with RNase-free DNase I (GBC, Beijing, China). RNA purity and concentration were measured with a NanoDrop ND-1000 Micro-Volume UV-Vis Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) at 260/280 nm ratio within the range of 1.8-2.2 and 260/230 nm around 2.0. After 1% agarose gel electrophoresis, the RNA integrity was checked. Following the manufacturer's instructions, 0.5 µg total RNA was used for first-strand cDNA synthesis using PrimeScript RT reagent Kit with gDNA Eraser (Perfect Real Time) (TaKaRa, Japan). The threefold dilution products of the generated cDNA were stored at −80 • C and used for qRT-PCR analyses.

qRT-PCR Analysis
qRT-PCR analyses were performed in 96-well blocks with a BIO-RAD CFX96 Real-Time PCR system (Bio-Rad, Hercules, California, USA) using NovoStart SYBR qPCR Supermix Plus (Novoprotein, China) in a 10-µL volume, containing 5 µL 2× NovoStart SYBR qPCR Supermix Plus, 1 µL diluted cDNA, 0.5 µL of forward primer (10 µmol/L), 0.5 µL of reverse primer (10 µmol/L), and 3 µL of ddH 2 O. The cycling conditions were as recommended by the manufacturer: 1 min at 95 • C, followed by 39 cycles of 95 • C for 5 s, 60 • C for 30 s, and 95 • C for 5 s. To confirm the specificity of the primers, melting curves were included after amplification. At the end of the reaction process, the melting curve was derived by heating the amplicon from 65 to 95 • C. All qRT-PCR analyses were run in technical quadruplicates and biological triplicates.

Stability Ranking of Candidate Reference Genes
Three different software programs (GeNorm [24], NormFinder [4], BestKeeper [26]) and the RefFinder (https://www.heartcure.com.au/reffinder/) web tool were used to evaluate the stability of 14 candidate reference genes under various stress conditions. Expression levels of the candidate reference genes were determined by quantification cycle (Cq) values. The three software programs of statistical analyses were conducted according to the manufacturer's instructions. RefFinder was Plants 2020, 9, 996 4 of 14 used to calculate the final rank of the 14 candidate reference genes. Results from CFX manager were exported into Microsoft Excel 2003 and transformed to create input files for each target according to the requirements of each software. For GeNorm and NormFinder, the Cq values were transformed into relative quantities using the formula 2 −∆Cq , in which ∆Cq = corresponding Cq value-minimum Cq value [29]. GeNorm identifies reference genes with good stability by calculating the M value, with a smaller M value indicating a better stability of the reference gene [26]. The program considers M values below 1.5 to indicate stable expression. The software can also calculate the V values, and the optimal number of reference genes for target gene expression normalization was decided by pairwise variation (V n /V n+1 ). A V n /V n+1 cutoff value of 0.15 indicates that an additional reference gene is not necessary [24,30]. The stability value calculated by NormFinder determines inter-and intragroup variation, and the lowest value indicates the highest stability. BestKeeper analyzes the stability of the candidate reference genes based on untransformed Cq values. The reference genes are considered to be the most stable when they exhibit the lowest CV and standard deviation (CV ± SD). Finally, RefFinder assigns an appropriate weight to each gene and calculates the geometric mean of their weights for the overall final ranking based on the rankings from each program.

Validation of Reference Genes by Expression Analysis of Cyt-Cu/Zn SOD and CAT Under Abiotic Stresses
To validate the reliability of the reference genes from software programs analysis, two target genes, namely Cyt-Cu/Zn SOD and CAT, were selected to analyze the expression patterns using the two most stable reference genes and the least stable reference gene. The results were calculated using the 2 −∆∆Cq method [31]. Three technical replicates were performed for each sample. The expression level of Cyt-Cu/Zn SOD in white clover was determined with forward primer 5 -AACTGTGTACCACGAGGACTTC-3 and reverse primer 5 -AGACTAACAGGTGCTAACAACG-3 , while the expression level of CAT was determined with forward primer 5 -AACAGGACGGGAATAGCACG-3 and reverse primer 5 -ACCA GGTTCAGACACGGAGACA-3 .

Verification of PCR Amplicons, Primer Specificity, and Gene-Specific PCR Amplification Efficiency
The amplicon sizes of 14 reference genes were checked by testing each primer pair using electrophoresis on a 2% agarose gel. The PCR products showed that the 80-199 bp fragments of ACT7, ACT101, TUA1109, TUB, CYP, 60SrRNA, UBQ, E3, GAPDH1, GAPDH2, PP2A, BAM3, SAMDC, and ABC were clearly amplified, consistent with the expected fragment sizes, and no impurities or primer dimers were observed ( Figure 1). Plants 2020, 9, x FOR PEER REVIEW 5 of 16 ABC were clearly amplified, consistent with the expected fragment sizes, and no impurities or primer dimers were observed ( Figure 1). Each melting curve of 14 candidate reference genes under various abiotic stresses only exhibited a single peak, and the gene amplification curves had good repeatability (Figure 2), showing that the primers were highly specific for later qRT-PCR and the results were reliable. Each melting curve of 14 candidate reference genes under various abiotic stresses only exhibited a single peak, and the gene amplification curves had good repeatability (Figure 2), showing that the primers were highly specific for later qRT-PCR and the results were reliable.

GeNorm Analysis
In GeNorm analysis, the M values were calculated to rank the average expression stability of 14 candidate reference genes ( Figure 4). Previous studies had confirmed that an M value below the threshold of 1.5 was considered as indicative of a suitable reference gene. A lower M value represents a higher degree of expression stability for the candidate reference gene. In this study, ACT101 and 60SrRNA were the most stable genes for all samples, including different tissues under various abiotic stresses, while GAPDH1 was the least stably expressed gene. For drought stress, E3 and SAMDC were the most stable genes in leaf samples, while ACT101 and PP2A were the most stable genes in root samples. Similarly, the ACT101 and SAMDC genes ranked the highest in terms of stability for leaf samples under salt stress, while ACT7 and TUA1109 were the most stable for root samples. For heat stress, the ACT7 and ACT101 genes were the most stable in leaf samples, while the TUB and GAPDH1 genes were the most stable in root samples. ACT7 and E3 were the most stable genes in cold-treated roots and leaves. For heavy metal stress, UBQ and PP2A were the most stable genes in leaves, while E3 and GAPDH2 were the most stable in roots.

GeNorm Analysis
In GeNorm analysis, the M values were calculated to rank the average expression stability of 14 candidate reference genes ( Figure 4). Previous studies had confirmed that an M value below the threshold of 1.5 was considered as indicative of a suitable reference gene. A lower M value represents a higher degree of expression stability for the candidate reference gene. In this study, ACT101 and 60SrRNA were the most stable genes for all samples, including different tissues under various abiotic stresses, while GAPDH1 was the least stably expressed gene. For drought stress, E3 and SAMDC were the most stable genes in leaf samples, while ACT101 and PP2A were the most stable genes in root samples. Similarly, the ACT101 and SAMDC genes ranked the highest in terms of stability for leaf samples under salt stress, while ACT7 and TUA1109 were the most stable for root samples. For heat stress, the ACT7 and ACT101 genes were the most stable in leaf samples, while the TUB and GAPDH1 genes were the most stable in root samples. ACT7 and E3 were the most stable genes in cold-treated roots and leaves. For heavy metal stress, UBQ and PP2A were the most stable genes in leaves, while E3 and GAPDH2 were the most stable in roots.
GeNorm procedure was also used for determining the optimal number of reference genes required for qRT-PCR. The optimal number of reference genes for target gene expression normalization was decided by pairwise variation (V n /V n+1 ). A 0.15 V n /V n+1 cutoff value indicates that an additional reference gene is not necessary. In this study, except for the case of all samples, the V 2 /V 3 values were less than 0.15 in leaves and roots under abiotic stress ( Figure 5), indicating that the combination of two reference genes was suitable. However, when all samples were analyzed together to determine the optimal number of reference genes, the pairwise variation of V 2 /V 3 was higher than 0.15, and the V 5 /V 6 was just 0.15, indicating that five reference genes should be used for gene expression studies in white clover including various stress conditions. Thus, it was more convenient to select optimal reference genes according to different experimental conditions. GeNorm procedure was also used for determining the optimal number of reference genes required for qRT-PCR. The optimal number of reference genes for target gene expression V2/V3 values were less than 0.15 in leaves and roots under abiotic stress ( Figure 5), indicating that the combination of two reference genes was suitable. However, when all samples were analyzed together to determine the optimal number of reference genes, the pairwise variation of V2/V3 was higher than 0.15, and the V5/V6 was just 0.15, indicating that five reference genes should be used for gene expression studies in white clover including various stress conditions. Thus, it was more convenient to select optimal reference genes according to different experimental conditions.

NormFinder Analysis
The stability value calculated by NormFinder determines inter-and intragroup variation, and the lowest value means the most stable. The expression stability calculated by NormFinder for each gene showed that E3 was the top ranked gene in leaves under drought stress, while UBQ was the most stable reference gene in roots under drought, heat, and cold stress and for all samples. The 60S gene ranked highest in leaves under salt and heat stress. Meanwhile, SAMDC was the best in roots under salt stress, and TUA1109 was the top ranked gene in leaves under cold stress. For heavy metal stress, UBQ and ACT101 were the most stable reference genes in leaves and roots, respectively ( Table  2).

NormFinder Analysis
The stability value calculated by NormFinder determines inter-and intragroup variation, and the lowest value means the most stable. The expression stability calculated by NormFinder for each gene showed that E3 was the top ranked gene in leaves under drought stress, while UBQ was the most stable reference gene in roots under drought, heat, and cold stress and for all samples. The 60S gene ranked highest in leaves under salt and heat stress. Meanwhile, SAMDC was the best in roots under salt stress, and TUA1109 was the top ranked gene in leaves under cold stress. For heavy metal stress, UBQ and ACT101 were the most stable reference genes in leaves and roots, respectively ( Table 2).

BestKeeper Analysis
BestKeeper software was used to synchronously analyze the untransformed Cq values, which reflect the stability of the candidate reference genes. The reference genes are considered to be the most stable when they exhibit the lowest CV ± SD. The results indicated that TUB and GAPDH1 were the two most stable genes in our study. GAPDH1 ranked the highest in leaves under drought, salt, heat, and cold stress and in roots under drought and salt stress. Meanwhile, TUB ranked the highest in leaves under heavy metal stress and roots under heat, cold, and heavy metal stress. For all samples, TUB was also the most stable reference gene (Table 3).

RefFinder Analysis
Finally, RefFinder was used to assign an appropriate weight to each gene and calculate the geometric mean of their weights for the overall final ranking, based on the rankings from each program ( Table 4). For all samples and drought-treated leaves, the two top-ranked genes determined by RefFinder method were the same as those determined by GeNorm. Both UBQ and 60S were the most suitable reference genes in drought-treated samples. UBQ and SAMDC were identified as the most stable reference genes in leaf samples under salt stress, and SAMDC and 60S were identified as the most stable reference genes in root samples under salt stress. For heat stress, the most stable combinations were 60S plus ACT7 in leaves and UBQ and BAM3 in roots. TUA1109 and ACT7 were the most stable genes in cold-treated leaf samples, while TUB and UBQ were the most stable genes in cold-treated root samples. For heavy metal stress, the most stable combinations were UBQ plus CYP in leaves and ACT101 plus E3 in roots.

Validation of the Reference Genes Identified from this Study
The relative expression levels of Cyt-Cu/Zn superoxide dismutase (SOD) and catalase (CAT) genes were used to validate the performance of the identified reference genes in this study. SOD, which catalyzes superoxide to H 2 O 2 and O 2 , initiates the defense system by removing superoxide, and it can be classified into three distinct groups by their metal cofactors: Cu/Zn, Mn, and Fe [32]. Cu/Zn SOD is present in the cytosol and chloroplasts [33]. Transgenic tobacco and cotton overexpressing chloroplastic Cu/Zn SOD and chloroplast-targeted MnSOD showed enhanced photosynthetic rates under chilling stress [34,35]. CAT reacts with H 2 O 2 directly to form H 2 O and O 2 . In most species, SOD and CAT activities are relatively sensitive in response to various abiotic stresses [36].
The relative expression levels of Cyt-Cu/Zn SOD and CAT genes were normalized using the two most stable reference genes and the least stable reference gene in white clover at different times. As shown in Figure 6, the normalized expression level of SOD in roots increased at 12 days and then decreased under cold treatment when using the two most stable genes (TUB and UBQ), while the expression level at 12 days was extremely low when BAM3 was used as a reference gene. In response to drought stress, the expression levels of CAT in leaves were similar at 6 and 8 days. However, the relative expression decreased at 10 days when the two most stable genes (SAMDC and E3) were adopted. Meanwhile, a 10-to 11-fold higher expression level occurred at 10 days when using ACT7 as a reference gene.
Plants 2020, 9, x FOR PEER REVIEW 13 of 16 adopted. Meanwhile, a 10-to 11-fold higher expression level occurred at 10 days when using ACT7 as a reference gene.
(a) (b) Figure 6. Relative expression levels of target genes. (a) Relative expression levels of Cyt-Cu/Zn SOD under cold stress using the two most stable reference genes and the least stable reference gene for normalization in white clover root tissues at different times; (b) relative expression levels of CAT under drought stress using the two most stable reference genes and the least stable reference gene for normalization in white clover leaf tissues at different times.

Discussion
In the process of plant growth and development, it is inevitable to face a lot of adversities. Conventional breeding techniques usually take a long time for selecting valuable genes with stable expression, but transgenic technology could improve the efficiency greatly. When studying the molecular mechanisms of stress resistance in plants and cloning stress resistance genes, real-time quantitative PCR (qRT-PCR) is needed to analyze the expression of a target gene. The accurate determination of relative gene expression mainly depends on the reference genes [37]. Therefore, the selection of suitable reference genes can reduce the experimental error [38]. Previous studies have demonstrated that there is no "universal" reference gene applicable for various experimental conditions [39]. Thus, it is necessary to select matched reference genes for quantitative real-time PCR in white clover involved in various abiotic stresses.
In this study, we observed the growth process of plants under five abiotic stresses. Leaf and root tissues under mild, moderate, and severe stress were sampled. Fourteen frequently used reference genes were picked out, and their stability was analyzed by four software programs (GeNorm, NormFinder, BestKeeperand RefFinder) under five experimental conditions (drought, salt, cold, heat, and heavy metal stress). Notably, the different algorithms evaluating the expression stability of reference genes selected different stable genes due to their different mathematical calculations [40]. Furthermore, RefFinder was used to integrate and generate the comprehensive ranking of the candidate reference genes based on the geometric mean of the weights of every gene calculated by each program [27]. The RefFinder results show us the overall ranking order, which has been widely used to select suitable reference genes in previous studies. The GeNorm results showed that it was better to select two reference genes in most experimental conditions. Finally, we concluded that the top two reference genes as ranked by the RefFinder program should be selected. In order to validate these selected candidate reference genes in white clover, the relative expression levels of Cyt-Cu/Zn SOD and CAT genes were normalized using the two most stable reference genes and the least stable reference gene. The validation results suggested that using inappropriate reference genes may significantly increase the error of target gene expression and make the results unreliable.
Furthermore, we determined that there was no single reference gene that exhibits a constant expression level in all samples of various tissues and under different experimental conditions; this was consistent with previous research [41]. Rafael Narancio [15] determined that EF1a, followed by ACT11 and UBQ, was the most stably expressed gene across organs and treatments in white clover. From our results, ACT and UBQ also showed a high stability across most experimental conditions.

Discussion
In the process of plant growth and development, it is inevitable to face a lot of adversities. Conventional breeding techniques usually take a long time for selecting valuable genes with stable expression, but transgenic technology could improve the efficiency greatly. When studying the molecular mechanisms of stress resistance in plants and cloning stress resistance genes, real-time quantitative PCR (qRT-PCR) is needed to analyze the expression of a target gene. The accurate determination of relative gene expression mainly depends on the reference genes [37]. Therefore, the selection of suitable reference genes can reduce the experimental error [38]. Previous studies have demonstrated that there is no "universal" reference gene applicable for various experimental conditions [39]. Thus, it is necessary to select matched reference genes for quantitative real-time PCR in white clover involved in various abiotic stresses.
In this study, we observed the growth process of plants under five abiotic stresses. Leaf and root tissues under mild, moderate, and severe stress were sampled. Fourteen frequently used reference genes were picked out, and their stability was analyzed by four software programs (GeNorm, NormFinder, BestKeeperand RefFinder) under five experimental conditions (drought, salt, cold, heat, and heavy metal stress). Notably, the different algorithms evaluating the expression stability of reference genes selected different stable genes due to their different mathematical calculations [40]. Furthermore, RefFinder was used to integrate and generate the comprehensive ranking of the candidate reference genes based on the geometric mean of the weights of every gene calculated by each program [27]. The RefFinder results show us the overall ranking order, which has been widely used to select suitable reference genes in previous studies. The GeNorm results showed that it was better to select two reference genes in most experimental conditions. Finally, we concluded that the top two reference genes as ranked by the RefFinder program should be selected. In order to validate these selected candidate reference genes in white clover, the relative expression levels of Cyt-Cu/Zn SOD and CAT genes were normalized using the two most stable reference genes and the least stable reference gene. The validation results suggested that using inappropriate reference genes may significantly increase the error of target gene expression and make the results unreliable.
Furthermore, we determined that there was no single reference gene that exhibits a constant expression level in all samples of various tissues and under different experimental conditions; this was consistent with previous research [41]. Rafael Narancio [15] determined that EF1a, followed by ACT11 and UBQ, was the most stably expressed gene across organs and treatments in white clover. From our results, ACT and UBQ also showed a high stability across most experimental conditions. However, under certain conditions, the most stable genes may be different from other species. Therefore, it is necessary to choose the most suitable reference gene for a more accurate result according to different experimental conditions, as the expression of a target gene is bound to change when using different reference genes.
In conclusion, 14 candidate reference genes were first selected in white clover. However, the optimal reference genes for different tissues (leaves and roots) under different experimental conditions are not identical. For all samples, 60SrRNA and ACT101 were the two top-ranked genes. Under drought stress, SAMDC and UBQ were identified as the most stable reference genes in the leaf and root samples, respectively. UBQ, SAMDC, and 60SrRNA were suggested as suitable reference genes in salt stress. For heat stress, the most stable gene was 60SrRNA in leaves, while UBQ was the most stable in roots. TUA was the most stable gene for cold-treated leaf samples, while TUB was the most stable gene for cold-treated root samples. For heavy metal stress, UBQ was the most stable gene in leaves, while ACT101 was the most stable gene in roots. For the first time, we analyzed the most stable reference genes for different tissues in white clover under five different abiotic stresses, providing the most suitable reference for gene expression analysis in later research. This is of great significance and will be helpful in exploring the potential molecular mechanisms of the abiotic stress response in white clover.