Selection of Suitable Reference Genes for RT-qPCR Gene Expression Analysis in Siberian Wild Rye (Elymus sibiricus) under Different Experimental Conditions

Elymus sibiricus, which is a perennial and self-pollinated grass, is the typical species of the genus Elymus, which plays an important role in forage production and ecological restoration. No reports have, so far, systematically described the selection of optimal reference genes for reverse transcriptase quantitative real-time polymerase chain reaction (RT-qPCR) analysis in E. sibiricus. The goals of this study were to evaluate the expression stability of 13 candidate reference genes in different experimental conditions, and to determine the appropriate reference genes for gene expression analysis in E. sibiricus. Five methods including Delta Ct (ΔCt), BestKeeper, NormFinder, geNorm, and RefFinder were used to assess the expression stability of 13 potential reference genes. The results of the RefFinder analysis showed that TBP2 and HIS3 were the most stable reference genes in different genotypes. TUA2 and PP2A had the most stable expression in different developmental stages. TBP2 and PP2A were suitable reference genes in different tissues. Under salt stress, ACT2 and TBP2 were identified as the most stable reference genes. ACT2 and TUA2 showed the most stability under heat stress. For cold stress, PP2A and ACT2 presented the highest degree of expression stability. DNAJ and U2AF were considered as the most stable reference genes under osmotic stress. The optimal reference genes were selected to investigate the expression pattern of target gene CSLE6 in different conditions. This study provides suitable reference genes for further gene expression analysis using RT-qPCR in E. sibiricus.


Introduction
Elymus is the largest genus of the tribe Triticeae with worldwide distribution, including approximately 150 species [1]. E. sibiricus (Siberian wild rye), which is the typical species of the genus Elymus, is a perennial, allotetraploid, and self-pollination grass [2,3]. It is one of the most important forages in Northern China due to its high protein content, strong adaptability, superior cold, and drought tolerance [4]. Meanwhile, E. sibiricus were widely used in artificial grassland and ecological governance in recent years [5]. It is noteworthy that severe seed shattering is the main reason for seed yield losses in E. sibiricus. According to previous reports, the degree of seed shattering in E. sibiricus the petri dish for 30 days, and then transplanted to 10 cm (diameter) flowerpots filling with 40% peat soil, 40% vermiculite, and 20% washed sand. The seedlings were grown in the greenhouse under 12 h of a photoperiod at 28 • C/14 • C day and night temperatures and 30% relative air humidity. Routine management was carried out during the whole process of growth.

Treatments and Tissue Collection
Six accessions were used to analyze the expression of reference genes, and three clones with consistent growth of each genotype were selected as three biological replicates. A total of 138 samples were collected from different genotypes, developmental stages, organs, and stress treatments (Table 1).
Specially, for gene expression analysis across different genotypes, leaf samples were collected from six different genotypes (ZN, XH, MQ, HZ, LQ, and LT) at 4-6 leaves stages. Leaf samples were collected from genotype ZN at different developmental stages: seedling, tillering, jointing, heading, and flowering. Tissue samples were collected from root, stem, leaf, and inflorescence at the flowering stage. For stress treatment, 4-6 leaf stage plants were removed from the soil and transferred to Hoagland's solution. For salt treatment, plants were treated by adding 250 mmol/L NaCl to Hoagland's solution for 0, 12, and 24 h. For osmotic stress, plants were treated by adding 20% PEG6000 to Hoagland's solution for 0, 12, and 24 h. For heat stress, plants were moved to a growth chamber at 40 • C for 0, 12, and 24 h. For cold stress, plants were moved to another growth chamber at 3 • C for 0, 12, and 24 h. All samples were frozen in liquid nitrogen and stored in a −80 • C refrigerator for later use.

RNA Extraction
Total RNA was extracted from each tissue using the plant Total RNA Kit (Sangon Biotech, Shanghai, China), according to the manufacturer's instructions. The concentration and purity of extracted RNA were measured by NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). Only the RNA absorbing ratio of 1.8-2.0 at OD260 nm/OD280 nm and 2.0-2.6 at OD260 nm/OD230 nm were used for subsequent analysis. In addition, the integrity of RNA was evaluated by 1% agarose gel electrophoresis.

First-Strand cDNA Synthesis
The cDNA synthesis was performed with the M-MLV cDNA synthesis kit (Sangon Biotech, Shanghai, China), according to the manufacturer's instructions. First, one microgram RNA was mixed with 1 µL oligonucleotide dT primer (0.5 µg/µL) in a 1.5 mL centrifuge tube, and then RNase free double distilled water was added to a final volume of 12 µL. Second, the reaction mixture was centrifuged for 3 to 5 s. The tube was incubated at 65 • C for 5 min, which was followed by cooling in ice cubes for 30 sec. Then, the mixture was centrifuged for 3-5 s. Third, the 5× Reaction Buffer 4 µL, RNase Inhibitor (40 U/µL) 1 µL, dNTP Mix (10 mmol/L) 2 µL, and M-MuLV RT (200 U/µL) 1 µL were added to the tube. Then the reaction mixture was centrifuged for 3-5 s. The first-strand cDNA was synthesized by incubating at 42 • C for 45 min. Lastly, reverse transcription was terminated by heating the reaction mixture for 10 min at 70 • C. The cDNA sample was stored at −80 • C and diluted 1:20 using RNase free water before RT-qPCR analysis.

Selection of Candidate Reference Genes and Primer Design
In this study, we selected 13 potential reference genes and one target gene (Table 2). These genes were obtained from our transcriptome database (https://www.ncbi.nlm.nih.gov/sra/SRX2617497) [7] by BLAST search using the sequences of reported Arabidopsis and rice reference genes. The 14 genes were named based on their similarity to known nucleotide sequences with the BLAST score value greater than 600 and identity ranging from 89% to 97% [33]. The primers for RT-qPCR were designed via Primer Premier 6.0 using the following parameters: melting temperature (Tm) at 58-62 • C (optimum Tm of 60 • C), PCR product length at 100-300 bp, GC content at 45-55%, and length of primers at 18-25 bp (Table 3).

RT-qPCR
We conducted RT-qPCR in 96-well blocks with a Bio-Rad CFX96 real-time PCR system (Bio-Rad, Hercules, CA, USA). The final reaction volume was 10 µL, and each reaction mix contained 2xSG Fast qPCR Master Mix (Low Rox) (Sangon Biotech, Shanghai, China) 5 µL, 10 µM Forward Primer 0.2 µL, 10 µM Reverse Primer 0.2 µL, 30 ng/µL cDNA 1 µL, DNF Buffer 1 µL, and ddH 2 O 2.6 µL. The amplification procedure included a denaturation step at 95 • C for 3 min, which was followed by 40 cycles of 95 • C for 3 s and 60 • C for 30 s. After the cycling process, the melting curves of RT-qPCR amplifications were obtained by raising the temperature from 60 • C to 95 • C. We also selected target gene CSLE6 to assess the expression stability of reference genes. Each RT-qPCR reaction was carried out for the independent sample with three technical replicates.

Data Analysis
In the present study, the expression stability of candidate reference genes was evaluated with four algorithms: (geNorm v3.5) (https://genorm.cmgg.be/) [38], (NormFinder v0.953) (https://moma. dk/normfinder-software) [39], (BestKeeper v1.0) (https://www.gene-quantification.de/bestkeeper. html) [40], and Delta Ct [41], and then a comprehensive ranking was obtained by the RefFinder program [42]. The raw RT-qPCR data were obtained by the CFX equipment software, and the Cq (cycle quantification) values were used for further analysis. For geNorm and NormFinder methods, raw Cq values were converted into the relative quantities, according to the formula 2 −∆Cq (∆Cq = each corresponding Cq value-lowest Cq value) [43]. The expression stability (M) of potential reference genes was calculated by the geNorm algorithm based on the average pairwise variation of each gene with all other control genes [38]. Genes with lower M values reflect more expression stability. To determine the optimal number of reference genes for normalization, the pairwise variations (V) were calculated by geNorm. If the V n /V n+1 (n is the number of reference genes) value is below or equal to 0.15, the number of suitable reference genes is equal to n. The degree of variance within and between groups was evaluated via an ANOVA-based model in the NormFinder program, and the gene with the lowest stability value has the most stable expression [39]. For the BestKeeper method, the standard deviation (SD), the coefficient of variation (CV), and correlation coefficient (r) were calculated by Cq values and the more stable reference gene possessed the lower SD value [40]. The Delta Ct algorithm used the standard deviation (SD) to rank the stability of all reference genes, and the reference gene with the lowest SD value showed the most stable performance. The RefFinder is a web-based analysis tool. It can integrate the results from Delta Ct, BestKeeper, NormFinder, and geNorm analysis to generate a comprehensive ranking. The RefFinder program (https://github.com/fulxie/RefFinder) running in my computer works as a local server, and we deployed it to a Php-based server (Apache + PHP) at first. The standard curves were generated by using a series of two-fold diluted cDNA as templates for each primer. The correlation coefficient (R 2 ) and slope were obtained from the linear regression model in Microsoft Excel 2016, and the PCR efficiency (E) was calculated according to the formula: E = (10 −1/slope − 1) × 100% [44]. The regression coefficient (R 2 ) should be greater than 0.98 with the amplification efficiency (E) over 90% but less than 110%. To calculate the relative expression level of the target gene, the 2 −∆∆Ct method [45] was applied, and the significant difference analysis was conducted by the SPSS statistical software v20.0.

Verification of Primer Specificity and PCR Amplification Efficiency
In the present study, we selected 13 potential reference genes based on transcriptome data of E. sibiricus to investigate the expression stability under different conditions. Subsequently, specific primers for RT-qPCR were designed according to the sequences of 13 candidate reference genes. In order to check the specificity of all primers, agarose gel electrophoresis (2%) and melting curve analysis of the amplification products were conducted ( Figure 1). The results of agarose gel electrophoresis showed that each primer amplified a single product of the expected size. Additionally, a single peak for each primer was observed from the melting curve analysis. These results indicated that all of the primers in this experiment were specific. The amplification efficiency (E) of 14 primers varied from 90.35% (DNAJ) to 101.08% (PP2A), and the regression coefficient (R 2 ) ranged from 0.981 (ACT2) to 0.998 (GAPDH, TUB3 and U2AF).

Cq Data Collection and Variations in Reference Genes
In order to analyze the expression level of all potential reference genes under different experimental conditions, the Cq values of all samples were obtained by RT-qPCR ( Figure 2). The Cq values of 13 candidate reference genes ranged from 19.3 (TUA2 expressed in the stem sample) to 39.44 (eIF-3A expressed in the root sample under osmotic stress for 24 h) among all samples. Furthermore, the top three highly expressed genes were HIS3 (mean Cq = 24.45), TUB3 (mean Cq = 25.32), and TUA2 (mean Cq = 25.49). Three lowly expressed genes were eIF-3A (mean Cq = 31.14), U2AF (mean Cq = 30.12), and eIF-3C (mean Cq = 29.57). The coefficient of variation (CV) of the Cq values reflects the expression stability of the reference gene, and the reference gene with a low CV value represents a high degree of stable expression. Therefore, eIF-3C (CV = 5.40%) possessed the most stable expression due to the lowest CV value, and the GAPDH (CV = 18.87%) with the highest CV value had the most unstable expression.

Expression Stability Analysis of the Candidate Reference Genes
To determine the most stable reference gene of E. sibiricus, four common algorithms (geNorm, NormFinder, BestKeeper, and Delta Ct) [38][39][40][41] and an integrated analysis tool (RefFinder) [42] were applied to assess the expression stability of all candidate reference genes. Moreover, the top six stable reference genes from the four algorithms were shown in Figure S1.

GeNorm Analysis
The expression stability of 13 potential reference genes in the different genotypes, developmental stages, tissues, and different abiotic stress treatments was evaluated by the geNorm algorithm ( Figure 3). According to the threshold of M value recommended by the geNorm program, a candidate gene could be used as a reference gene for expression analysis only when its M value is under 1.5. Our results showed that DNAJ and eIF-3C were the most stable genes with the lowest M value (0.03) in different genotypes, while U2AF was the least stable gene with the highest M value (0.52) ( Figure 3A). In different developmental stage samples ( Figure 3B  In order to obtain the optimal number of reference genes in different conditions, pairwise variation (V) was calculated by the geNorm program ( Figure 4). Moreover, 0.15 was used as the threshold value to determine the optimal number of reference genes. For non-stress treatment groups (different genotypes, different developmental stages and different tissues), all pairwise variations except for V 12/13 (0.265) of different tissues were lower than 0.15, which demonstrated that one reference gene was sufficient for gene expression analysis. For the stress groups (salt stress, heat stress, cold stress, and osmotic stress), however, it showed different results. The V 6/7 (0.137) value from salt stress was lower than 0.15, which indicated that the top six reference genes (ACT2, TBP2, PP2A, TUB3, U2AF, and TEF2) were required for expression analysis. The V 2/3 (0.120) value from heat stress samples was lower than 0.15, which suggested that the top two reference genes (CYP19 and TUA2) were sufficient for normalization. The V 3/4 (0.144) value from the cold stress group was lower than 0.15, which indicated that the top three genes (TEF2, U2AF, and PP2A) were required for normalization. Nevertheless, the threshold value of 0.15 should not be regarded as a rigorous standard, and higher cut-off values of V n /V n+1 were also found in several reports [35,43,46]. A minor variation was found between V 2/3 (0.156) and V 3/4 (0.153) in samples from osmotic stress, which revealed that two reference genes (DNAJ and U2AF) were sufficient for expression analysis. Similarly, V 3/4 (0.223) and V 4/5 (0.225) from all samples showed slight variation, which illustrates that the top three genes (PP2A, TBP2, and DNAJ) were needed for normalization.

NormFinder Analysis
The 13 candidate reference genes were ranked based on the expression stability value calculated using the NormFinder program ( Figure 5 and Table S2). The most stable gene possessed the highest expression stability value. Accordingly, TBP2 had the highest ranking in different genotypes and different tissues with the stability values of 0.10 and 0.07, respectively. PP2A was the most stable reference gene in different developmental stages, cold stress, and all samples with the stability values of 0.11, 0.23 and 0.29, respectively. As for the salt stress, heat stress, and osmotic stress, ACT2 was identified as the most stably expressed reference genes with the stability values of 0.15, 0.13, and 0.22, respectively. Although the most stable reference genes according to the NormFinder algorithm were different from that of geNorm in most of the experimental sets, the least stable reference genes in all groups were consistent between NormFinder and geNorm. For example, according to the geNorm analysis, DNAJ and CYP19 were the most stable reference genes in different genotypes and heat stress, whereas their stability rankings were ninth and fourth in the NormFinder analysis, respectively. Furthermore, U2AF, ACT2, GAPDH, and eIF-3A had lower expression stability among all experimental groups by using the geNorm and NormFinder program.

BestKeeper Analysis
The BestKeeper program was used to analyze the expression stability of 13 reference genes according to the values of standard deviation (SD), the coefficient of variation (CV), and correlation coefficient (r). The lower SD value of genes represented the higher expression stability ( Figure 6 and Table S3). Hence, HIS3 was the most stable reference gene in different genotypes, while it ranked fifth and seventh in geNorm and NormFinder analysis, respectively. DNAJ exhibited the best expression stability in different developmental stages, which ranked third and fifth in geNorm and NormFinder analysis, respectively. In different tissues, TBP2 was the most stable reference gene in both BestKeeper and NormFinder analysis, but it ranked fourth in geNorm analysis. For salt stress, heat stress, cold stress, osmotic stress, and all samples, eIF-3C presented the most stable expression with the lowest SD values, while it had low ranking in geNorm and NormFinder analysis. U2AF, ACT2, GAPDH, and eIF-3A were considered to be the least stable genes in all experimental conditions due to their highest SD values, which was consistent with the results from geNorm and NormFinder analysis.

Delta Ct Analysis
The values of standard deviation (SD) was used as the indicator for evaluating the expression stability of reference genes using the Delta Ct method. A gene with the lowest SD value indicated the most stable reference gene (Figure 7 and Table S4). Accordingly, TBP2, PP2A, ACT2, and DNAJ showed the most stable expression in experimental conditions. Meanwhile, U2AF, ACT2, GAPDH, and eIF-3A were the least stable genes in all sample sets, which was consistent with geNorm, NormFinder, and BestKeeper analysis.

RefFinder Analysis
RefFinder, which is a comprehensive analysis tool for expression stability of reference genes, was used to calculate the synthetic rankings of 13 potential reference genes based on four approaches including geNorm, NormFinder, BestKeeper, and Delta Ct ( Figure S2 and Table 4). For different genotypes, TBP2, HIS3, and CYP19 were the top three stable genes according to the RefFinder algorithm analysis. TUA6, PP2A, and HIS3 were considered as the three most stable genes for different developmental stages. In different tissues, TBP2, PP2A, and eIF-3A showed high expression stability. ACT2, TBP2, and PP2A were identified as the most stably expressed reference genes under salt stress. ACT2, TUA2, and CYP19 had the best expression stability under heat stress. PP2A, ACT2, and DNAJ exhibited a high level of expression stability under cold stress. As for osmotic stress, DNAJ, U2AF, and ACT2 were considered as the most stable genes. For all samples, PP2A, TBP2, and DNAJ were the most stable genes. However, U2AF was the least stable reference genes in different genotypes and different developmental stages. ACT2 was identified as the least stable reference genes in different tissues. Under salt stress, heat stress, and osmotic stress and for all samples, GAPDH had the worst expression stability. eIF-3A showed the lowest level of expression stability under cold stress. The results of the above five methods indicated that the most stable genes differed when different experimental sample sets were compared, which illustrates the significance of suitable reference genes that are specific to each experimental condition. Therefore, the top two stable reference genes and the least reference gene in different experimental conditions were selected for validating the expression stability of reference genes, according to the RefFinder analysis.

Discussion
RT-qPCR is an effective technique for gene expression analysis due to its high accuracy, simplicity, specificity, and sensitivity [8][9][10]. It is reported that inappropriate reference genes will lead to opposite conclusions in gene expression analysis [24]. Consequently, evaluating the expression stability of potential reference genes is necessary to quantify the expression level of target genes, and to analyze the expression pattern of genes of interest [47]. Numerous studies suggested that none of the reference genes maintain the consistent expression stability among various experimental conditions, and it is imperative to carry out reference gene screening under specific experimental conditions [11,29,[48][49][50][51]. To our knowledge, there have been no studies regarding the selection of appropriate reference genes in E. sibiricus, and, thus, this study is the first report with respect to the systematic selection and evaluation of reference genes for RT-qPCR normalization in E. sibiricus.
In the present study, 13 candidate reference genes were selected from transcriptome sequencing data to evaluate their expression stability in different experimental conditions [7]. The results showed that all reference gene primers possessed good amplification efficiency (90.35% to 101.08%), regression coefficient (0.981 to 0.998) and Cq values (25 to 30), which illustrates that the RT-qPCR data are suitable for further analysis. Four commonly used methods including geNorm, NormFinder, BestKeeper, and Delta Ct [38][39][40][41] were used to assess the expression stability of 13 candidate reference genes. The results suggested that the least stable reference genes of the same experimental samples were consistent among four algorithms (Table 4), while the most stable reference genes were inconsistent in different conditions ( Figure S1 and Table 4). For example, U2AF and ACT2 were unstable reference genes based on four algorithms in non-stress conditions. GAPDH and eIF-3A also had the worst stable performance under stress conditions. Furthermore, DNAJ was the most stable reference gene found by the geNorm algorithm in different genotypes, while it ranked ninth in NormFinder and Delta Ct analysis. The eIF-3C had the highest expression stability by the BestKeeper algorithm in four stress conditions, but it was considered the unstable reference gene in geNorm, NormFinder, and Delta Ct analysis. The heterogeneous result was also found in previous studies, and that may be caused by different algorithms [29,52]. Fortunately, we obtained an integrated evaluation result of potential reference genes by using RefFinder. RefFinder is regarded as a comprehensive analysis tool, which has extensive recognition in determining the optimal reference genes for gene expression analysis [12,21,[53][54][55]. Lastly, we adopted the rankings derived from the RefFinder method as the ultimate evaluation results of expression stability of candidate reference genes, and selected suitable reference genes for validating the expression stability based on the results of RefFinder.
Previous studies demonstrated that there are no any reference genes with consistent expression stability among different conditions [11,29,[48][49][50][51]. Our results were similar to previous studies. For example, ACT2 exhibited a high degree of expression stability under four stress conditions, whereas ACT2 was the least stable reference gene in different tissues. TUA2 was the most stable reference gene in different developmental stages, but it ranked eleventh under osmotic stress. These results indicated that the importance of selecting appropriate reference genes in different experimental conditions. The plant cell wall provides a guarantee for plants to adapt to various environmental conditions. Cellulose, which is the main component of cell walls in the plant, is synthesized by cellulose synthases (CesA) [56]. A considerable number of cellulose synthase-like (CSL) genes have similarities with CesA gene [57], and, thus, the CSL gene also plays an important role in cellulose synthesis [58]. Previous research studies indicated that the expressions of cellulose synthesis genes are closely related to the adaptation of plants under different stress conditions [56]. In addition, other studies showed that cellulose synthesis genes were associated with shedding or shattering in leaves, flowers, fruits, and seed [7,59,60]. CSLE6, which is a member of the superfamily of genes referred to as glycosyltransferase 2 (GT2), has the function of synthesizing cellulose in the plant cell wall [61]. However, there are few studies regarding the CSLE6 gene. This gene was identified in Arabidopsis (At1g55850, named AtCSLE1) and rice (LOC_Os09g30130, named OsCSLE6), respectively [57,62]. In the present study, to validate the expression stability of reference genes, we selected CSLE6 as a target gene to investigate the expression pattern of CSLE6 in different conditions. Several previous reports suggested that selecting two or more stable reference genes to calculate the relative expression levels of target genes will generate more reliable results [21,25,26,63]. Hence, we selected two stable reference genes and the least stable reference gene in different conditions for the expression pattern analysis of CSLE6. As shown in Figure 8, CSLE6 exhibited distinct expression levels in different experimental conditions. In different genotypes, the expression level of CSLE6 in the ZN (low seed shattering) genotype was higher than the XH genotype (high seed shattering). This is similar to the OsCel9D [60], which is a seed shattering related gene. This gene indicates that CSLE6 may have similar functions in seed shattering. On the one hand, CSLE6 directly hampers the abscission process in seed shattering by changing the cell wall components, such as cellulose content and pectin content [60]. On the other hand, CSLE6 indirectly reduces the seed shattering by affecting the lignin biosynthesis that is closely related to seed shattering [64,65]. In different developmental stages, the order of expression levels of CSLE6 was tillering > flowering > jointing > seedling > heading. This is different from the previous study of E. sibiricus [7], likely due to different tissues and periods. Plants generate a large number of new leaves and cells during the tillering stage, which may explain that CSLE6 had the highest expression in the tillering stage. For different tissues, the order of expression levels of CSLE6 is leaf > root > stem > inflorescence > tiller bud. This was similar to Arabidopsis [62], which the highest expression was found in the leaf sample. This was followed by the root sample. The leaf samples in this analysis were not young leaves (flowering stage) and the cellulose accounts for a high proportion in the cell wall, where the high expression of CSLE6 may play a role in replacing the homogalacturonan (HGA) in the cell wall [62,66]. Furthermore, previous reports revealed that cellulose synthesis genes usually have high expression in tissues with cell division and expansion, such as root and hypocotyl [64,67]. Under salt stress, the expression of CSLE6 in leaf and root samples showed a trend of decreasing first and then increasing, which is different from that in Arabidopsis [68]. For cold stress, the expression level of CSLE6 was down-regulated in leaf and root samples. Many studies indicated that abiotic stress (e.g., salt stress and cold stress) will reduce the cell expansion and inhibit plant growth [69][70][71]. In addition, the expression levels of a great number of cellulose synthesis genes were reduced by cold stress or salt stress in Arabidopsis [62,67], poplar [72], and cotton [73], while these genes showed increased expression in rice [74]. Hence, the cellulose synthesis genes have distinct expression patterns in different experimental conditions and different species because of the multiple ways in which plants adapt to the environment. More importantly, the above results illustrated that the diverse expression patterns of CSLE6 reflected the adaptation mechanisms for E. sibiricus under a variety of environmental conditions.

Conclusions
In summary, the most stable reference genes were different under distinct experimental conditions in this study. To obtain the precise results of gene expression analysis, it is recommended to adopt suitable reference genes in specific experimental conditions. One or more stable reference genes should be selected to investigate the expression pattern of target genes based on the comprehensive evaluation results from RefFinder. The expression pattern of CSLE6 may provide a basis for studying the resistance mechanism in E. sibiricus. Moreover, our study screened several suitable reference genes in specific conditions for E. sibiricus, and offered some guidelines for the selection of reference genes for other plant species.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/6/451/s1. Figure S1: The top six most stable reference genes generated by delta Ct, BestKeeper, NormFinder, and geNorm. The purple, yellow, green, pink, and circles each contain the top six most stable reference genes of delta Ct, BestKeeper, NormFinder, and geNorm, respectively. The genes in the overlap area are the ones confirmed as the top six most stable reference genes by more than one algorithm.