Reference Gene Selection for Quantitative Real-Time Reverse-Transcriptase PCR in Annual Ryegrass (Lolium multiflorum) Subjected to Various Abiotic Stresses

To select the most stable reference genes in annual ryegrass (Lolium multiflorum), we studied annual ryegrass leaf tissues exposed to various abiotic stresses by qRT-PCR and selected 11 candidate reference genes, i.e., 18S rRNA, E2, GAPDH, eIF4A, HIS3, SAMDC, TBP-1, Unigene71, Unigene77, Unigene755, and Unigene14912. We then used GeNorm, NormFinder, and BestKeeper to analyze the expression stability of these 11 genes, and used RefFinder to comprehensively rank genes according to stability. Under different stress conditions, the most suitable reference genes for studies of leaf tissues of annual ryegrass were different. The expression of the eIF4A gene was the most stable under drought stress. Under saline-alkali stress, Unigene14912 has the highest expression stability. Under acidic aluminum stress, SAMDC expression stability was highest. Under heavy metal stress, Unigene71 expression had the highest stability. According to the software analyses, Unigene14912, HIS3, and eIF4A were the most suitable for analyses of abiotic stress in tissues of annual ryegrass. GAPDH was the least suitable reference gene. In conclusion, selecting appropriate reference genes under abiotic stress not only improves the accuracy of annual ryegrass gene expression analyses, but also provides a theoretical reference for the development of reference genes in plants of the genus Lolium.


Introduction
Annual ryegrass (Lolium multiflorum) is cultivated in temperate and subtropical regions worldwide and is used in forage and livestock systems as silage and green fodder owing to its high palatability and digestibility [1,2]. It is one of the most widely cultivated plants in the world, and grows in warm and humid climates. It is mainly used for artificial grass, clipping, and grazing [3]. However, environmental degradation from natural challenges, saline-alkali, acidic aluminum, drought, and heavy metals threatens the production of annual ryegrass. To enable more detailed and in-depth studies of gene expression, it is necessary to identify appropriate reference genes.
Quantitative reverse transcription PCR (qRT-PCR) is widely used to determine the expression of genes in different plant tissues owing to its high sensitivity, specificity, accuracy, and repeatability [4]; it is a powerful technical method to detect gene expression levels. However, the accuracy of qRT-PCR ( Figure 1), indicating that the primers for the 11 candidate reference gene were highly specific and can be used for qRT-PCR analyses. Most of the candidate reference genes in various abiotic stress conditions exhibited melting curves with a single peak, showing that the primers are highly specific. Gene amplification curves for each sample had good repeatability, showing that the qRT-PCR results were accurate and reliable ( Figure 2). However, in some samples, a miscellaneous peak of GAPDH indicated that the primer was not highly stable. Most of the candidate reference genes in various abiotic stress conditions exhibited melting curves with a single peak, showing that the primers are highly specific. Gene amplification curves for each sample had good repeatability, showing that the qRT-PCR results were accurate and reliable ( Figure 2). However, in some samples, a miscellaneous peak of GAPDH indicated that the primer was not highly stable.

Analysis of the Expression Stability of Reference Genes under Different Types of Abiotic Stress
Analysis of the raw expression levels across all samples identified variation among reference genes is shown in Figure 3. We present detailed data in Supplementary Table S1. Ct values for the 11 genes ranged from 11.24 to 37.73, and the majority of these values were between 23.44 and 31.94. 18S rRNA was highly expressed compared to the other genes, reaching threshold fluorescence after only 13.6 amplification cycles, whereas the Ct average of all reference genes within the datasets was approximately 27 cycles. As a result, the 18S rRNA transcripts were around 10,000-fold more abundant than the average for the dataset. Unigene71 and Unigene14912 transcripts were least abundant, with Ct values of 34.03 and 33.75, respectively. The individual reference genes had different expression ranges across all studied RNA samples. However, before we performed the software analysis, we observed that GAPDH has the most discrete data and therefore its ranking should be lowest in subsequent analyses. By contrast, E2, eIF4A, HIS3, TBP-1, and Unigene14912 were expected to perform well in subsequent analyses. The wide expression ranges of the eleven tested reference genes confirmed that it is of the utmost importance to select a reliable reference gene to normalize gene expression under certain conditions.

Stability Ranking of Candidate Reference Genes for Different Abiotic Stresses
In this study, GeNorm, BestKeeper, NormFinder, and RefFinder were used to analyze the expression of relatively stable reference genes in different abiotic stress conditions. The results showed that reference gene exhibited expression differences depending on the abiotic stress.

GeNorm Analysis
GeNorm was used to rank the reference genes with respect to expression stability by calculating average pairwise expression ratios (M1), where values below 1.5 indicate stable expression. In this study, Unigene77 and Unigene755, with the same M1 value of 0.962, were the most stably expressed in all samples, while GAPDH, with an M1 value of 2.257, was the least stably expressed gene ( Figure  4E). For the samples subjected to acidic aluminum stress, drought stress, heavy metal stress, and saline-alkali stress, the pairs SAMDC and Unigene14912, E2 and TBP-1, E2 and Unigene71, and Unigene755 and Unigene14912 ranked the highest, respectively, suggesting that these gene pairs were stable and could be used as references ( Figure 4A-D).
In some gene expression analysis of study, sometimes to express quantitatively the high precision requirements, or in some less obvious purpose gene expression quantity change, using a single reference genes cannot accord with the requirement of the experiment. Therefore, it has become a new criterion in quantitative PCR research to correct the expression of the target gene with two or more stable reference genes [21]. Accordingly, this study analyzed the pairwise variation of each reference gene ( Figure 5), and the results showed that the pairwise variation of the two to eleven reference gene was less than 0.15. Vandesompele et al. [17] suggested that when the paired variation V value was less than 0.15, there was no need to add multiple reference genes for the correction. Therefore, we can see that the optimal number of reference genes is two. Figure 3. qRT-PCR Ct values for all candidate reference genes in annual ryegrass leaf samples under various abiotic stress conditions. Variation is displayed as medians (lines), 25th to 75th percentiles (boxes), and ranges (whiskers).

Stability Ranking of Candidate Reference Genes for Different Abiotic Stresses
In this study, GeNorm, BestKeeper, NormFinder, and RefFinder were used to analyze the expression of relatively stable reference genes in different abiotic stress conditions. The results showed that reference gene exhibited expression differences depending on the abiotic stress.

GeNorm Analysis
GeNorm was used to rank the reference genes with respect to expression stability by calculating average pairwise expression ratios (M 1 ), where values below 1.5 indicate stable expression. In this study, Unigene77 and Unigene755, with the same M 1 value of 0.962, were the most stably expressed in all samples, while GAPDH, with an M 1 value of 2.257, was the least stably expressed gene ( Figure 4E). For the samples subjected to acidic aluminum stress, drought stress, heavy metal stress, and saline-alkali stress, the pairs SAMDC and Unigene14912, E2 and TBP-1, E2 and Unigene71, and Unigene755 and Unigene14912 ranked the highest, respectively, suggesting that these gene pairs were stable and could be used as references ( Figure 4A-D).
In some gene expression analysis of study, sometimes to express quantitatively the high precision requirements, or in some less obvious purpose gene expression quantity change, using a single reference genes cannot accord with the requirement of the experiment. Therefore, it has become a new criterion in quantitative PCR research to correct the expression of the target gene with two or more stable reference genes [21]. Accordingly, this study analyzed the pairwise variation of each reference gene ( Figure 5), and the results showed that the pairwise variation of the two to eleven reference gene was less than 0.15. Vandesompele et al. [17] suggested that when the paired variation V value was less than 0.15, there was no need to add multiple reference genes for the correction. Therefore, we can see that the optimal number of reference genes is two.

BestKeeper Analysis
BestKeeper was used to analyze the expression stability of genes by determining which genes had the lowest CV ± SD ( Note: expression stability and ranking of 11 candidate reference genes calculated with BestKeeper under all sample conditions, acidic aluminum, drought, heavy metal, and saline-alkali stress. The stability of eleven reference genes was identified based on the coefficient of variation (CV) and standard deviation (SD).

Figure 5.
Pairwise variation (V) measure of the candidate reference genes. When V n /V n+1 < 0.15, then the optimal number of reference genes is N.

BestKeeper Analysis
BestKeeper was used to analyze the expression stability of genes by determining which genes had the lowest CV ± SD (Table 1). Our results indicated that Unigene71 had the lowest CV values under acidic aluminum and heavy metal stress; Unigene14912 showed the lowest CV values for drought stress, saline-alkali stress, and all samples. However, GAPDH showed the highest CV values of 14.04 ± 3.44, 11.84 ± 2.90, 14.26 ± 3.31 and 12.12 ± 2.82 when plant materials were subjected to acidic aluminum, drought, heavy metal, and saline-alkali stress, respectively, and had the highest CV values for all samples (Table 1).

NormFinder Analysis
NormFinder was used to rank gene expression stability based on comparisons of the average pairwise variation of a gene to all other genes (M 2 ). The stability values of reference genes were calculated by NormFinder and are shown in Table 2. Under different conditions, the lowest M 2 values were assigned to different candidate reference genes. Unigene14912 ranked highest for all samples with an M 2 value of 0.889. For acidic aluminum stress, 18S rRNA was the most-stable gene with an M 2 value of 0.545. Under drought, heavy metal, and saline-alkali stress, eIF4A, Unigene71, and TBP-1 were the most suitable genes, with M 2 values of 0.273, 0.26 and 0.455, respectively.

RefFinder Analysis
Finally, RefFinder was used produce a comprehensive ranking of the most stable candidate reference genes under each stress condition ( Table 3).

Validation of the Usefulness of the Reference Genes Identified from This Study
To verify that reference genes influence qRT-PCR results for annual ryegrass leaves, the expression patterns of delta-1-pyrroline-5-carboxylate (P5CS1) and Cyt-Cu/Zn superoxide dismutase (Cyt-Cu/Zn SOD) were analyzed in varied abiotic stresses by using five different reference genes identified in previous results including three most stable reference genes, the combination of reference gene and the worst reference gene. P5CS1 is a critical gene for proline biosynthesis rate-limiting enzymes regulating the balance of osmotic and preventing the cell turgor pressure increases to adapt severe environment [22][23][24]. SOD is a kind of metalloenzyme commonly found in plants. When plants are under heavy metal stress, SOD is the first line of defense against oxidative damage [25]. SOD will work for hydrogen peroxide and oxygen molecules by disproportionation on superoxide anion, which is efficient removal of oxygen free radicals in plants and protects cells from oxidative damage [26]. Advanced plants have Cu/Zn SOD, which exists in the cytoplasm and chloroplasts [27].
We utilized the 2 −∆∆Ct method to calculate the expression of P5CS1 and Cyt-Cu/Zn SOD under different stresses. As shown in Figure 6A,B, P5CS1 had same expression patterns using eIF4A + TBP-1, eIF4A, E2 and TBP-1 for normalization under drought stress. The similar results were obtained in Cyt-Cu/Zn SOD expression patterns. Specially, using TBP-1 for normalization, the expression of P5CS1 was induced with a 2-fold increase on the first day, 4-fold increase on the 3rd day, and 9-fold increase on the 6th day.
Under acidic aluminum stress, P5CS1 and Cyt-Cu/Zn SOD had same expression patterns using HIS3 + Unigene14912, HIS3, and Unigene14912 for normalization ( Figure 6C,D). Using Unigene14912 for normalization, the expression of P5CS1 was induced with a 3-fold increase on the first day, 5-fold increase on the 3rd day, and 7-fold increase on the 6th day, while the expression of Cyt-Cu/Zn SOD was induced with a 1.5-fold increase on the first day, 3-fold increase on the 3rd day, and 6-fold increase on the 6th day.
The HIS3 had the best efficient for analyzing the expression patterns of P5CS1 and Cyt-Cu/Zn SOD under saline-alkali stress. The expression of P5CS1 was induced with a 3-fold increase on the first day, a 5-fold increase on the 3rd day, and a 35-fold increase on the 6th day ( Figure 6E,F). The similar results were obtained for Cyt-Cu/Zn SOD.
As the results show ( Figure 6G,H) that P5CS1 had same expression patterns using Unigene71, HIS3, and E2 for normalization under heavy metal stress. The same results were obtained in Cyt-Cu/Zn SOD that both two target genes were up-regulated. Of note, using Unigene71 for normalization, the expression of P5CS1 was induced with a 5-fold increase on the first day, 8-fold increase on the 3rd day, and 29-fold increase on the 6th day.

Validation of the Usefulness of the Reference Genes Identified from This Study
To verify that reference genes influence qRT-PCR results for annual ryegrass leaves, the expression patterns of delta-1-pyrroline-5-carboxylate (P5CS1) and Cyt-Cu/Zn superoxide dismutase (Cyt-Cu/Zn SOD) were analyzed in varied abiotic stresses by using five different reference genes identified in previous results including three most stable reference genes, the combination of reference gene and the worst reference gene. P5CS1 is a critical gene for proline biosynthesis ratelimiting enzymes regulating the balance of osmotic and preventing the cell turgor pressure increases to adapt severe environment [22][23][24]. SOD is a kind of metalloenzyme commonly found in plants. When plants are under heavy metal stress, SOD is the first line of defense against oxidative damage [25]. SOD will work for hydrogen peroxide and oxygen molecules by disproportionation on superoxide anion, which is efficient removal of oxygen free radicals in plants and protects cells from oxidative damage [26]. Advanced plants have Cu/Zn SOD, which exists in the cytoplasm and chloroplasts [27].
We utilized the 2 −ΔΔCt method to calculate the expression of P5CS1 and Cyt-Cu/Zn SOD under different stresses. As shown in Figure 6A,B, P5CS1 had same expression patterns using eIF4A + TBP-1, eIF4A, E2 and TBP-1 for normalization under drought stress. The similar results were obtained in Cyt-Cu/Zn SOD expression patterns. Specially, using TBP-1 for normalization, the expression of P5CS1 was induced with a 2-fold increase on the first day, 4-fold increase on the 3rd day, and 9-fold increase on the 6th day.
Under acidic aluminum stress, P5CS1 and Cyt-Cu/Zn SOD had same expression patterns using HIS3 + Unigene14912, HIS3, and Unigene14912 for normalization ( Figure 6C,D). Using Unigene14912 for normalization, the expression of P5CS1 was induced with a 3-fold increase on the first day, 5-fold increase on the 3rd day, and 7-fold increase on the 6th day, while the expression of Cyt-Cu/Zn SOD was induced with a 1.5-fold increase on the first day, 3-fold increase on the 3rd day, and 6-fold increase on the 6th day.
The HIS3 had the best efficient for analyzing the expression patterns of P5CS1 and Cyt-Cu/Zn SOD under saline-alkali stress. The expression of P5CS1 was induced with a 3-fold increase on the first day, a 5-fold increase on the 3rd day, and a 35-fold increase on the 6th day ( Figure 6E,F). The similar results were obtained for Cyt-Cu/Zn SOD.
As the results show ( Figure 6G,H) that P5CS1 had same expression patterns using Unigene71, HIS3, and E2 for normalization under heavy metal stress. The same results were obtained in Cyt-Cu/Zn SOD that both two target genes were up-regulated. Of note, using Unigene71 for normalization, the expression of P5CS1 was induced with a 5-fold increase on the first day, 8-fold increase on the 3rd day, and 29-fold increase on the 6th day.

Discussion
qRT-PCR has become an important technique for the analysis of gene expression owing to its sensitivity, accuracy, and high-throughput nature [28]. The accuracy of qRT-PCR depends on the stability of reference genes used for data normalization [29]. Under different conditions, changes in reference gene expression may lead to incorrect analyses of real-time data. Hence, it is important to identify suitable reference genes for gene expression analyses by qRT-PCR for each condition, rather than relying on reference genes from published data for different plant species [8].
With rapid developments in biotechnology, traditional reference genes have not met the needs of researchers. Accordingly, more stable reference genes are needed for gene expression analyses to detect very small differences in gene expression. In this study, seven traditional reference genes and

Discussion
qRT-PCR has become an important technique for the analysis of gene expression owing to its sensitivity, accuracy, and high-throughput nature [28]. The accuracy of qRT-PCR depends on the stability of reference genes used for data normalization [29]. Under different conditions, changes in reference gene expression may lead to incorrect analyses of real-time data. Hence, it is important to identify suitable reference genes for gene expression analyses by qRT-PCR for each condition, rather than relying on reference genes from published data for different plant species [8].
With rapid developments in biotechnology, traditional reference genes have not met the needs of researchers. Accordingly, more stable reference genes are needed for gene expression analyses to detect very small differences in gene expression. In this study, seven traditional reference genes and four transcriptome Unigenes were analyzed to assess traditional reference genes and explore unknown and stable genes. We compared three different statistical applets, GeNorm, BestKeeper, and NormFinder to evaluate the 11 reference genes in annual ryegrass. The three applets yielded different rankings of stable reference genes depending on treatment (Table 3). It is worth noting that Unigene14912 is the most stable reference gene for all samples and under saline-alkali stress. Additionally, Unigene71 was the most stable of the reference gene under heavy metal stress. By contrast, the traditional reference genes SAMDC and eIF4A only performed well under acidic aluminum and drought stress, respectively. In this study, two new stable reference genes, Unigene14912 and Unigene71, were successfully identified, providing more potential reference genes for future studies.
BestKeeper and NormFinder identified Unigene14912 as the most stable reference gene, followed by eIF4A for all samples. GeNorm identified Unigene77 and Unigene755 as the best reference genes, followed by HIS3. GAPDH ranked last using all three methods. The most stable reference genes varied among the three methods used in this study, but the identification of unstable genes was consistent, and the inconsistent results may be explained by the different models associated with each statistical algorithm [30,31].
The results of this study were consistent with results obtained for other species. Three methods showed that under abiotic stress, annual ryegrass eIF4A was highly stable, similar to the results for closely related species, e.g., Lolium perenne under abiotic stress [32], darnel (Lolium temulentum) [33], and the Kentucky bluegrass [9], indicating that it is the most suitable for reference gene. GAPDH exhibits similar results across species, such as papaya (Carica papaya) [34], rice (Oryza sativa) [35], and tobacco (Nicotiana tabacum) [36], in which it was not suitable for use as a reference gene, but is involved in other biochemical reactions during abiotic stress.
The transcript abundances of Unigenes were generally low. In order to improve reliability of the result, the combination of traditional and novel reference genes are recommended. In the validation experiment under saline-alkali and heavy metal stress, the standard deviation of the combination of the reference gene for normalization is smaller than single reference gene, which indicating that the combination of the reference gene was more accurate and powerful in gene expression analysis. The most suitable reference gene combination was HIS3 + TBP-1 and Unigene71 + HIS3 under saline-alkali stress and heavy metal stress, respectively.
The material we use is a single variety of leaves. Although it has been able to meet the research needs, it is not perfect. In future studies, we will explore between different varieties, different sites, and different developmental stages, reference genes perform excellent.

Materials
This experimental subject is "Chuannong NO.1", developed by the Sichuan Agricultural University (Chengdu, China). Seeds were sown in 25 × 19 × 6 cm pots and the matrix was silica sand. They were placed in an incubator with a light intensity of 100 µmol/(m·s), 25 • C and day/night cycles of 16 and 8 h. When the leaf counts of the plants reached 5-6, stress treatments were initiated. Under drought stress, plants were treated with 15% PEG6000 for 6 days. For acidic aluminum treatment, plants were treated with 0.1 mol/L AlCl 3 (pH 2.9) for 6 days. For heavy metal treatment, the nutrient solution was filled with 200 mg/L Cr 6+ Hoagland's for 6 days. For saline-alkali stress, plants were treated with 0.1 mol/L NaHCO 3 (pH 8.2) for 6 days. We take a sample of several seedlings (upper leaves) in the same pot. Take three mixed samples at four time points (0, 1, 3, 6 days) under four treatments (drought, saline-alkali, acidic aluminum, heavy metal). There are three biological repeats for each sampling point. We treat the before treatment (0 day) as a control group in each treatment, as described by Huang L et al. [37]. All samples frozen in liquid nitrogen, and stored at −80 • C for later use. The 48 annual ryegrass samples during the study and which datasets they were included in Supplementary Table S2.

Reagents and Instruments
We used the following reagents. RNA extraction reagent from the RNA Simple Total RNA Kit (TianGen Biotech Co., Ltd., Beijing, China), reverse transcriptase reagent from the iScript cDNA Synthesis Kit (Bio-Rad Laboratories Inc., Hercules, CA, USA), and quantitative fluorescence PCR kitswere purchased from TAKARA (Shiga, Japan). Primers were obtained from Sangon Biological Engineering (Shanghai), Beijing branch (China). A fluorescence quantitative instrument was used for fluorescence quantitative PCR (CFX-96; Bio-Rad).

RNA Extraction and First-Strand cDNA Synthesis
The RNA simple Total RNA Kit (TianGen Biotech Co., Ltd., Beijing, China) was used to extract RNA from annual ryegrass material and RNase-free DNase I (GBC, Beijing, China) was used to remove DNA. The Ultramicro Spectrophotometer (NanoVue Plus Spectrophotometer, Wilmington, DE, USA) was used to detect RNA purity and concentration at A260 nm/A280 nm and A260 nm/A230 nm. After 1% agarose gel electrophoresis, the iScript cDNA Synthesis Kit (Bio-Rad Laboratories Inc.) was used to reverse transcribe total RNA samples into cDNA, and samples were diluted 1:5 with ultrapure water and stored at −80 • C for later use.

qRT-PCR Amplification Procedure
In a fluorescence quantitative PCR tube (TLS-0851; Bio-Rad), 2 µL of cDNA (30 ng/µL), 1.5 µL ofreverse primer (10 µmol/L), 1.5 µL of forward primer (10 µmol/L), 10 µL 2 × SYBR Premix Ex Taq (5 U/µL), and 5 µL of ddH 2 O were added for a total volume of 20 µL. The procedure was repeated 3 times for each biological sample and 4 times for each technology. The amplification procedure was as follows: 95 • C for 30 s, followed by 95 • C for 5 s, and 64 • C for 30 s, repeated 40 times. Then, an extension phase was performed from 60 • C to 95 • C, where the temperature for each cycle increased by 0.5 • C for 5 s to obtain Tm and fluorescent signals for the melting curve. All qRT-PCRs were run in quadruplicate for technical samples and triplicate for biological samples.

Data Analysis
qRT-PCR analyses and GeNorm, NormFinder, and BestKeeper were used to analyze the stability of 11 candidate reference genes. Finally, RefFinder was used to sequence the 11 reference genes and select the gene with the highest stability under different stress conditions. The top three ranked genes for heavy metal stress and lowest ranked gene were used to calculate P5CS1 and Cyt-Cu/Zn SOD gene expression to evaluate the effectiveness of the reference genes.

Conclusions
In this study, the expression levels of 11 candidate reference genes under various types of stress were analyzed in annual ryegrass. For different abiotic stresses, the most stable reference genes differed. The eIF4Aexpression was most stable under drought stress; the most appropriate reference gene is HIS3+TBP-1 under saline-alkali stress; the expression of HIS3 is best suited for the reference gene under acidic aluminum stress; Unigene71 + HIS3 is the optimal choice of reference gene under heavy metal stress. Unigene14912 was the best gene when considering all samples. The four treatments differed with respect to the most stable gene, suggesting that no single gene can be consistently used for data normalization across different conditions [40], but we found that HIS3, Unigene14912, E2 and TBP-1 were stable in two or more treatments. Based on the analysis of GeNorm software for pairwise variation, we concluded that the combination of Unigene14912 and HIS3 should be selected as the reference gene. The study results indicate genes that can be used as suitable candidate reference genes under different stresses in annual ryegrass, and provide a basis for Lolium expression analyses.
Supplementary Materials: The supplementary materials are available online, Table S1: The Ct values was collected from each treatment and time point by triplicate for biological, Table S2: The 48 annual ryegrass samples during the study and which datasets they were included in.