Selection of the reference genes for quantitative gene expression by RT-qPCR in the desert plant Stipagrostis pennata

The desert pioneer plant Stipagrostis pennata plays an important role in sand fixation, wind prevention, and desert ecosystem recovery. An absence of reference genes greatly limits investigations into the regulatory mechanism by which S. pennata adapts to adverse desert environments at the molecular and genetic levels. In this study, eight candidate reference genes were identified from rhizosheath development transcriptome data from S. pennata, and their expression stability in the rhizosheaths at different development stages, in a variety of plant tissues, and under drought stress was evaluated using four procedures, including geNorm, NormFinder, BestKeeper, and RefFinder. The results showed that GAPDH and elF were the most stable reference genes under drought stress and in rhizosheath development, and ARP6 and ALDH were relatively stable in all plant tissues. In addition, elF was the most suitable reference gene for all treatments. Analysis of the consistency between the reverse transcription-quantitative PCR (RT-qPCR) and RNA sequencing data showed that the identified elF and GAPDH reference genes were stable during rhizosheath development. These results provide reliable reference genes for assuring the accuracy of RT-qPCR and offer a foundation for further investigations into the genetic responses of S. pennata to abiotic stress.


Results
Identification of candidate reference genes in different stages of rhizosheath development. The conditions of |log2FoldChange| < 1, q-value ≥ 0.05, and FPKM (Fragments per Kilobase of exon model per Million mapped reads) ≥ 6 were used to screen the RNA-Seq data, and the relatively low coefficient of variation (CV) of FPKM was set as a high standard at all sampling points, generating a total of eight candidate reference genes (Table 1). Detailed information of these genes, including unigene name, gene symbol, homologue locus, and E-value, via comparison with the homologous genes in rice is listed in Supplementary Table S1. The FPKM-based heatmap of the eight genes during different stages of rhizosheath development is provided in Fig. 1 and indicates that GAPDH, α-TUB, TIP41, Histone H3 (HIS-3), elF, and ARP6 exhibited stable expression in different stages of rhizosheath development, whereas ALDH and protein phosphatase 2A (PP2A) showed relative unsteady expression.
Verification of primer specificity and PCR amplification efficiency. Specific primers for the eight candidate reference genes were designed and amplified by RT-qPCR to verify the specificity using cDNA from the roots as a template. Referring to MIQE (Minimum Information for Publication of Quantitative Real-Time PCR Experiments) guidelines, we have provided information about the specific primers of the eight reference genes and RT-qPCR experiments (Supplementary Figs. S1, S2) 18 . According to the slope of the standard curve, the amplification efficiency and R 2 of the RT-qPCR assays were calculated, which indicated that the amplification efficiency ranged from 90.8% (TIP41) to 95.4% (PP2A), and the R 2 value ranged from 0.988 (PP2A) to 0.999 (HIS-3) ( Table 1). These results indicated that all the primers of the eight genes had high specificity and amplification efficiency and were thus suitable for further analysis.
Expression profiles of the eight candidate reference genes in S. pennata. RT-qPCR assays were performed with the designed specific primers for the eight candidate reference genes using the cDNA as template extracted from different tissues and treated materials of S. pennata. Based on the Ct values obtained from RT-qPCR, a box diagram was generated to reflect the differences in the expression levels of the eight candidate genes (Fig. 2). The average Ct values of the reference genes ranged between 16.89 and 29.12, which indicated that the expression levels of these reference genes differed in S. pennata. The results showed that ALDH had the smallest variation, followed by GAPDH, elF, and PP2A, while α-TUB had the largest variation. Therefore, ALDH, GAPDH, elF, and PP2A were identified as relatively stable genes.  Table 2). The geNorm program evaluates the expression stability of reference genes by calculating the M-value, with a high value representing low stability. The results showed that the M-values of the eight reference genes were all lower than 1.5 under drought stress, indicating high expression stability of these genes. In addition, TIP41 and elF (M-value of 0.221) under drought conditions, TIP41 and ALDH (M-value of 0.286) in all plant tissues, TIP41 and GAPDH (M-value of 0.212) in rhizosheath development, and elF and GAPDH (M-value of 0.623) in all samples were identified as the two most stable reference genes ( Fig. 3A-D). The geNorm software was also used to calculate the paired variation value V n /V n+1 (n represents the reference gene number) to determine the optimal number of reference genes for RT-qPCR standardization. The results showed that the V n /V n+1 values of the different sample groups were all lower than 0.15, suggesting that two reference genes were sufficient to complete the RT-qPCR normalization in S. pennata ( Fig. 3E-H).
The BestKeeper procedure was utilized to analyze the expression stability of the eight candidate reference genes, with high R-values and low CV ± SD values denoting better stability 19 . As shown in Table 3  www.nature.com/scientificreports/ CV ± SD = 6.43 ± 1.01) were identified as the most stable reference genes. ARP6 was relatively stable in the tissues (R = 0.973, CV ± SD = 9.73 ± 2.78) and all samples (R = 0.944, CV ± SD = 6.65 ± 1.76). As shown above, evaluation by these three procedures produced different results for stable reference genes from those reported in some studies [20][21][22][23] . To provide a comprehensive assessment of the eight candidate reference genes in different conditions, we performed further analysis by RefFinder 24 according to the geometric mean of the reference genes 20,25 to generate a final comprehensive ranking of the expression stability of the reference genes. The results showed that elF and GAPDH in rhizosheath development, GAPDH and elF under drought stress, and ARP6 and ALDH in the tissues were the most stable genes. For all samples, ARP6 and elF were the two most stable reference genes (Table 4). www.nature.com/scientificreports/ www.nature.com/scientificreports/

Validation of the reference genes by RT-qPCR.
To verify the expression stability of the identified reference genes, two genes, namely CL9727 and U3887, which are closely related to rhizosheath development, were selected for RT-qPCR to analyze the similarity between the RT-qPCR results and the RNA-Seq data. The two stable genes of elF and GAPDH and two unstable genes of ALDH and α-TUB screened by the comprehensive assessment ranking under the condition of rhizosheath development were used as reference genes for normalization to calculate the relative expression levels of CL9727 and U3887. The results indicated that the expression levels of CL9727 and U3887 detected by RT-qPCR were highly consistent with the RNA-Seq data. When the unstable reference genes of ALDH and α-TUB were selected for standardization, the expression levels of CL9727 and U3887 detected by RT-qPCR were significantly different from the RNA-Seq data (Fig. 4). Additionally, correlation analysis of the relative expression levels of CL9727 and U3887 normalized by elF and GAPDH between RT-qPCR and RNA-Seq was performed, indicating a strong positive correlation between the RT-qPCR results and RNA-Seq data (R 2 = 0.6061-0.9257), once again validating the stability and reliability of the two reference genes.

Discussion and conclusion
Reverse transcription qPCR has become a common method for gene expression analysis due to its high sensitivity, good repeatability, high specificity, and high throughput, using reliable internal reference genes for the correction of target gene data to obtain accurate results [25][26][27][28] . Due to the lack of reference gene information in S. pennata, the reliability and accuracy of related gene expression detection significantly limits further gene function and genetic studies. Housekeeping genes such as GAPDH and TUB have typically been used as internal controls for the normalization of RT-qPCR, obtaining high effectiveness and reliability in other plant species 12,13 . However, there is great variability in the expression stability of some reference genes as a result of differences among conditions, species, tissues, growth and development stages, and experimental treatments 12,[14][15][16][17] . In this study, a total of eight candidate reference genes, including ARP6, elF, PP2A, ALDH, GAPDH, TIP41, HIS-3, and α-TUB, were screened for suitability in S. pennata ( Table 1). The algorithms of geNorm, NormFinder, BestKeeper, and RefFinder are commonly utilized to assess the expression stability of candidate reference genes to identify the optimal reference genes, such as in soybean 4 , potato 29 , Plukenetia volubilis L. 30 , and banana 31 . Using these four programs, the eight screened reference genes showed varying expression stability under different conditions, with the different programs obtaining different results. The most stably expressed genes in the present study included GAPDH by NormFinder and BestKeeper under drought stress, elF and TIP41 by geNorm under drought stress, ARP6 by NormFinder and BestKeeper in all tissues, ALDH and TIP41 by geNorm in all tissues, elF by NormFinder and geNorm in rhizosheath development,  (Tables 2, 3; Fig. 3). The 18S rRNA and 25S rRNA genes have been widely used as internal controls for the normalization of RT-qPCR in rice under stress conditions. In NaCl-and mannitol-treated rice seedlings, 18S rRNA and 25S rRNA showed the most stable expression levels of all reference genes 14 . It was also found that the 26S rRNA gene in Arabidopsis thaliana and other herbaceous plants was stably expressed under 10% PEG treatment 32 . In this study, GAPDH and elF were identified as the most stable reference genes under PEG treatment (Table 4). Different reference genes have been used in different plants and tissues 14,[33][34][35] . The polyubiquitin genes UBQ4 and UBQ10 demonstrated the most stable expression in different tissues of Brachypodium distachyon 17 . GAPDH showed the best expression stability in different tissues and organs of Saccharum sp. 36 . GAPDH and EF1α exhibited better expression stability and were selected as optimal internal reference genes during fruit development in Lycium barbarum L. EF1α and TUA were used as two internal reference genes for gene expression analysis during fruit development in Amomum villosum Lour 37 . Our comprehensive analysis identified ARP6 and ALDH in the tissues, elF and GAPDH in rhizosheath development, and ARP6 and elF in all samples as suitable reference genes in S. pennata (Table 4). In Vitis amurensis Rupr. berries, the expression of GAPDH was more stable than in Triticum aestivum L. 38,39 , while ACT and UBI exhibited better expression stability in wheat but worse expression stability in Solanum lycopersicum L. 39,40 . Comparative analysis of the correlation between the RNA-Seq data and RT-qPCR results using the screened stable genes of elF and GAPDH in rhizosheath development, which was done to quantify the target gene expression of CL9727 and U3887, indicated their high stability (Figs. 4 and 5), validating that these screened genes constituted reliable and effective internal reference genes for the normalization of RT-qPCR in S. pennata.
In conclusion, we screened eight candidate reference genes for RT-qPCR normalization based on transcriptome datasets for S. pennata. Using the three procedures of geNorm, NormFinder, and BestKeeper, the comprehensive assessment analysis by RefFinder, as well as comparison and correlation analysis between the RT-qPCR results and RNA-Seq data, we identified reliable and suitable internal reference genes for RT-qPCR normalization under various conditions, thereby providing a foundation for further investigations of the genetic functions and regulatory mechanisms at the molecular level in S. pennata.   RT-qPCR analysis. The RT-qPCR was performed on a LightCycler ® 480 real-time PCR system (Roche Diagnostics, Mannheim, Germany) using the SYBR Green-based PCR assay. The total reaction volume was 20 µL, which included 10 µL 2× SuperReal PreMix Plus (SYBR Green, TANGEN BIOTECH, Beijing, China), 0.5 µL each of 10 µM forward and reverse gene-specific primer, 3 µL template (first-strand cDNA), and 6 µL ddH 2 O 2 .

Materials and methods
Amplifications were performed using the following program: initial denaturation at 94 °C for 2 min, followed by a cycling procedure of 30 s denaturation at 94 °C, 30 s annealing at 55 °C, 30 s extension at 72 °C, and then a final extension at 72 °C for 10 min. The RT-qPCR analysis was tested in three biological replicates. Relative gene expression levels were calculated using the 2 −ΔΔCt method 41 . Determination and validation of the expression stability of the reference genes. The geNorm 25 , NormFinder 42 , BestKeeper 43 , and RefFinder 44 tools are used to screen stable reference genes for data analysis based on raw quantification cycle (Cq) values. The geNorm program selects the stable reference gene by calculating the M-value of each reference gene and determines the number of optimal reference genes according to the V n /V n+1 value. The default value of V is 0.15. If V n /V n+1 < 0.15, the number of optimal reference genes is n, if V n / V n+1 > 0.15, the number of optimal internal reference genes is n + 1. The calculation principle of NormFinder is similar to that of GeNorm, and the most suitable internal parameter gene is selected according to the stability value, with the most suitable gene being the lowest stability value. The correlation coefficient (R), standard deviation (SD), and coefficient of variation (CV) of pairing between each gene can be calculated by BestKeeper, with a larger R value or smaller SD and CV values denoting better stability of the reference gene. We used the online RefFinder software (https:// github. com/ fulxie/ RefFi nder) for the comparative analysis of ΔCt. The geometric mean of the Ct values of all candidate reference genes was analyzed to rank the expression stability.
Validation of the expression stability of the reference genes. The CL9729 and U3887 genes in the transcriptome data showed a close connection with rhizosheath development. In order to verify the analysis results of the candidate genes, two stably expressed (GAPDH and elF) and two unstably expressed (ALDH and α-TUB) reference genes were selected for RT-qPCR validation for these two target genes (CL9729 and U3887).