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.

1. 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 was related to genotypes [3,6], the development of abscission zone [4], and several key genes involved in plant hormones, lignin biosynthesis, and cell wall-degrading enzymes (e.g., ETR, PAL, and CesA) found by transcriptome sequencing [7]. Currently, the regulatory mechanisms associated with seed shattering in E. sibiricus has not been clearly elucidated. As a result, the analysis of expression patterns of these key genes will contribute to a better understanding of regulatory mechanisms regarding seed shattering in E. sibiricus. The appropriate reference gene is crucial for the gene expression analysis in E. sibiricus.

Reverse transcriptase quantitative real-time polymerase chain reaction (RT-qPCR) is a sensitive, precise, specific, rapid, and repeatable method to investigate and quantify the gene expression level [8,9,10]. However, the application of RT-qPCR are commonly affected by several factors, including the quality of RNA (purity and integrity), the efficiency of cDNA (complementary DNA) synthesis, reaction parameters, and primer specificity [11,12]. Therefore, in order to obtain the accurate normalization of RT-qPCR, it is necessary to screen the suitable reference genes for correcting the above-mentioned experimental errors [13,14].

Currently, transcriptome sequencing has been widely used in plant biological research, which greatly facilitates the understanding of important molecular mechanisms in plants [15,16]. Simultaneously, transcriptome data provide important resources for identification and exploration of reference genes. An optimal reference gene should be expressed at a constant level across various experimental conditions such as different genotypes, developmental stages or tissues, and even different abiotic stress. Hence, housekeeping genes are usually selected as the reference genes without validating their performance since their expression levels were thought to be stable [17,18,19]. Nevertheless, a large number of research studies recently found that the common reference genes like glyceraldehydes-3-phosphate dehydrogenase (GAPDH), β-actin (ACT), and 18S ribosomal RNA (18S rRNA) had significantly different expression levels under different tissues or different experimental treatments [20,21,22,23]. Consequently, it is essential to validate the stability of potential reference genes, which can effectively identify whether candidate reference genes are suitable for current experimental conditions. Using improper reference genes will lead to inappropriate interpretations [24]. It is reported that the expression stability of reference genes showed divergence in different species. For instance, TUB was a suitable reference gene for normalization of RT-qPCR in Sedum alfredii [25] and Jatropha curcas [26], but had the worst expression stability in Rhododendron molle [27] and Camellia sinensis [28]. The reason for this phenomenon may be the difference in the primer design and evaluation method, the distinct experimental conditions, and the diversity of gene function among different species [29]. Therefore, screening the reference genes that possess a high degree of expression stability in specific species is crucial to obtain the accurate results of RT-qPCR. Up to now, many suitable reference genes have been reported from different plants, including Arabidopsis thaliana [30], rice [31], Tibetan hulless barley [32], soybean [21], orchardgrass [33], Seashore paspalum [34], Kentucky bluegrass [35], and Setaria viridis [36]. However, there is no report for the selection of suitable reference genes for RT-qPCR analysis in E. sibiricus.

The objective of this study is to screen suitable reference genes with stable expression under various experimental conditions. Thirteen candidate reference genes including actin2 (ACT2), cyclophilin19 (CYP19), DNAJ heat shock N-terminal domain-containing protein (DNAJ), eukaryotic translation initiation factor 3A (eIF-3A), eukaryotic translation initiation factor 3C (eIF-3C), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), histone 3.3 (HIS3), protein phosphatase 2A (PP2A), TATA-binding protein 2 (TBP2), translation elongation factor 2 (TEF2), tubulin α-2 (TUA2), tubulin β-3 (TUB3), and U2 auxiliary factor (U2AF) were selected from the transcriptome data of E. sibiricus [7].

2. Materials and Methods

2.1. Plant Materials and Growth Conditions

In the present study, six wild E. sibiricus accessions (Table S1) were selected based on a previous screening for agronomic traits in 28 E. sibiricus accessions [37]. The materials seeds were germinated in 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.

2.2. 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).

Table 1

Experimental sets in the present study.

Sample Sets	Tissue Type	Sampling Dates	Germplasm	Biological Replicates	Total Number of Samples
Different genotypes	Leaf	4–6 leaves stage	ZN, XH, MQ, HZ, LQ, LT	3	18
Different developmental stages	Leaf	Seedling, Tillering, Jointing, Heading, Flowering	ZN	3	15
Different tissues	Root, Stem, Leaf, Tiller bud, Inflorescence	Flowering stage	ZN	3	15
Salt stress (250 mmol/L NaCl)	Leaf, Root	4–6 leaves stage, 0, 12, 24 HAT 1	ZN	3	18
Cold stress (3 °C)	Leaf, Root	4–6 leaves stage, 0, 12, 24 HAT	ZN	3	18
Heat stress (40 °C)	Leaf, Root	4–6 leaves stage, 0, 12, 24 HAT	ZN	3	18
Osmotic stress (20% PEG6000)	Leaf, Root	4–6 leaves stage, 0, 12, 24 HAT	ZN	3	18
Control	Leaf, Root	4–6 leaves stage, 0, 12, 24 HAT	ZN	3	18

Abbreviation:1 HAT, hours after treatment.

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.

2.3. 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.

2.4. 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.

2.5. 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).

Table 2

Description of 13 reference genes and one target gene.

Gene	Gene Description	Arabidopsis Homolog Locus	Amino Acid Identity with E. sibiricus (%)	Rice TIRG Identifier	Amino Acid Identity with E. sibiricus (%)
ACT2	Actin2	At5g09810	98.41	LOC_Os11g06390	98.94
CYP19	Cyclophilin 19	At2g29960	75.12	LOC_Os06g49480	82.73
DNAJ	DNAJ heat shock N-terminal domain-containing protein	At1g76700	66.67	LOC_Os08g41110	86.04
eIF-3A	Eukaryotic translation initiation factor 3A	At4g11420	62.36	LOC_Os01g03070	79.39
eIF-3C	Eukaryotic translation initiation factor 3C	At3g56150	62.18	LOC_Os07g03230	82.43
GAPDH	Glyceraldehyde-3-phosphate dehydrogenase	At3g04120	74.93	LOC_Os08g03290	92.12
HIS3	Histone 3.3	At4g40030	82.93	LOC_Os06g04030	100.00
PP2A	Protein phosphatase 2A	At1g13320	74.28	LOC_Os09g07510	78.36
TBP2	TATA binding protein 2	At1g55520	64.50	LOC_Os03g45410	59.65
TEF2	Translation elongation factor 2	At1g56070	91.10	LOC_Os02g32030	94.90
TUA2	Tubulin α-2	At4g14960	96.67	LOC_Os11g14220	97.56
TUB3	Tubulin β-3	At1g75780	89.04	LOC_Os06g46000	90.13
U2AF	U2 auxiliary factor	At5g42820	74.48	LOC_Os09g31482	69.97
CSLE6	Cellulose synthase-like protein E6	At1g55850	47.92	LOC_Os09g30130	76.64

Table 3

Details of primers and amplicons for 14 genes.

Candidate Reference Gene	Primer Sequence F/R (5′–3′)	Tm (°C)	Amplicon Length/Bp	Efficiency (%)	R2
ACT2	F: CCACGAGACGACCTACAATTCCATC	60.0	206	98.46	0.981
	R: CTCCGATCCAGACACTGTACTTCCT
CYP19	F: GGTGGTGAATCAATCTACGGCACAA	60.4	141	93.73	0.996
	R: GCTCGTGGTTACAGTGGTGATGAAG
DNAJ	F: GCAATGGCGTCAATGGCTTCAC	59.6	209	90.35	0.992
	R: GCATCACTAAGTCTGGACACCTCAG
eIF-3A	F: GAATCAGGCACAAGCACTGGAAGA	59.9	284	95.4	0.984
	R: ACAACCTACGGAACTCGGTGGAT
eIF-3C	F: GAATCATAAGGCTGCTGCGAAGGT	60.0	183	92.84	0.990
	R: AACGGTGGTGGTCCTCTATTGTCA
GAPDH	F: GTTACTGTCTTCGGCGTCAGGAAC	60.2	137	97.05	0.998
	R: ACCTTCTTGGCACCACCCTTCA
HIS3	F: CTACAACTGGAGGAGTGAAGAAGCC	59.6	164	98.85	0.993
	R: GAAGCGGAGGTCAGTCTTGAAGTC
PP2A	F: GTGATAATGAGGCTGAAGTGCGGAT	60.0	140	101.08	0.990
	R: GCGAACATGCTGAGACGAATCTGA
TBP2	F: CGGATGAGGCAGCCAAAGATTGT	59.6	103	98.65	0.993
	R: TGTTCTCGAAGGCAGTGTATGTCTC
TEF2	F: AGAAGTCCTGCCGTACCGTTATGA	60.1	107	98.97	0.995
	R: GCCATCATCAATAGCCTCAGCCAAT
TUA2	F: TGGTGATGTTGTGCCGAAGGATG	59.8	195	100.6	0.996
	R: ACGACACTGGTGGAGTTGGAGAT
TUB3	F: GTGCAGAACAAGAACTCGTCCTACT	59.5	233	99.97	0.998
	R: TCGGTGAACTCCATCTCGTCCAT
U2AF	F: ATCGCTGCTCTCGCATCCATAAC	59.7	281	98.5	0.998
	R: TGCTGCTGCCTGATCTTCCTCT
Target gene
CSLE6	F: AAGGATGGTGGAATGGACAGAGGAT	60.0	149	98.37	0.994
	R: CTTGGACTCGTCTTCGTCGCTTAC

2.6. 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 ddH2O 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.

2.7. 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 Vn/Vn+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) 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 (R2) 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.

3. Results

3.1. 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 (R2) ranged from 0.981 (ACT2) to 0.998 (GAPDH, TUB3 and U2AF).

Figure 1

Primer specificity and amplicon size of 13 candidate reference genes. (A) Melting curves of 13 candidate reference genes exhibiting single peaks. (B) Agarose gel electrophoresis (2%) showing specific amplification products of expected size using Real Time-qPCR. M: 500 bp marker. 1–13: ACT2, CYP19, DNAJ, eIF-3A, eIF-3C, GAPDH, HIS3, PP2A, TBP2, TEF2, TUA2, TUB3, and U2AF.

3.2. 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.

Figure 2

RT-qPCR Cq values for all candidate reference genes in all E. sibiricus samples. Whisker caps, boxes, lines, and square boxes represent maximum/minimum, 25/75 percentiles, median, and mean, respectively.

3.3. 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.

3.3.1. 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), the most stable genes were ACT2 and HIS3 (M = 0.29), while U2AF was the least stable gene (M = 0.99). eIF-3A and PP2A (M = 0.34) showed the most stability while ACT2 (M = 1.47) was identified as the least stable gene in different tissues (Figure 3C). ACT2 and TBP2 (M = 0.44) were the most stable genes under salt stress (Figure 3D), while GAPDH (M = 2.16) was the least stable gene. CYP19 and TUA2 (M = 0.30) had the most stability, whereas GAPDH (M = 1.46) had the worst stability in heat stress samples (Figure 3E). For the cold stress group (Figure 3F), TEF2 and U2AF (M = 0.21) were the most stable genes while eIF-3A (M = 1.59) was the least stable gene. Under osmotic stress (Figure 3G), DNAJ and U2AF (M = 0.38) showed the highest degree of expression stability while GAPDH (M = 1.83) presented the lowest degree of expression stability. Meanwhile, PP2A and TBP2 (M = 0.93) were the most stable genes in all samples and GAPDH (M = 2.00) was the least stable gene. In summary, the most stable gene was different in various experimental sample sets.

Figure 3

Average expression stability values (M) of 13 candidate reference genes under different conditions calculated by geNorm. (A) Different genotypes, (B) different developmental stages, (C) different tissues, (D) salt stress, (E) heat stress, (F) cold stress, (G) osmotic stress, and (H) all samples.

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 V12/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 V6/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 V2/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 V3/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 Vn/Vn+1 were also found in several reports [35,43,46]. A minor variation was found between V2/3 (0.156) and V3/4 (0.153) in samples from osmotic stress, which revealed that two reference genes (DNAJ and U2AF) were sufficient for expression analysis. Similarly, V3/4 (0.223) and V4/5 (0.225) from all samples showed slight variation, which illustrates that the top three genes (PP2A, TBP2, and DNAJ) were needed for normalization.

Figure 4

Pairwise variation (V) of 13 candidate reference genes under various experimental conditions calculated by geNorm.

3.3.2. 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.

Figure 5

Expression stability of E. sibiricus potential reference genes calculated using NormFinder. (A) Different genotypes, (B) different developmental stages, (C) different tissues, (D) salt stress, (E) heat stress, (F) cold stress, (G) osmotic stress, and (H) all samples. Low to high expression stability is represented over a spectrum from red to blue, respectively.

3.3.3. 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.

Figure 6

Expression stability of candidate reference genes calculated by BestKeeper. (A) Different genotypes, (B) different developmental stages, (C) different tissues, (D) salt stress, (E) heat stress, (F) cold stress, (G) osmotic stress, and (H) all samples. Low to high expression stability is represented over a spectrum from red to blue, respectively.

3.3.4. 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.

Figure 7

Expression stability for 13 candidate reference genes calculated via Delta Ct. (A) Different genotypes, (B) different developmental stages, (C) different tissues, (D) salt stress, (E) heat stress, (F) cold stress, (G) osmotic stress, and (H) all samples. Low to high expression stability is represented over a spectrum from red to blue, respectively.

3.3.5. 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.

Table 4

Stability ranking of 13 candidate reference genes calculated by Delta Ct, BestKeeper, NormFinder, GeNorm, and RefFinder.

Method	Ranking Order (Better—Good—Average)
	1	2	3	4	5	6	7	8	9	10	11	12	13
Different genotypes
Delta Ct	TBP2	HIS3	eIF-3A	CYP19	PP2A	ACT2	TUA2	TEF2	DNAJ	eIF-3C	TUB3	GAPDH	U2AF
BestKeeper	HIS3	ACT2	CYP19	DNAJ	eIF-3C	TUA2	TBP2	TEF2	TUB3	eIF-3A	PP2A	GAPDH	U2AF
NormFinder	TBP2	eIF-3A	TUA2	PP2A	CYP19	TEF2	HIS3	ACT2	DNAJ	eIF-3C	TUB3	GAPDH	U2AF
GeNorm	DNAJ | eIF-3C		CYP19	ACT2	HIS3	TBP2	TUA2	eIF-3A	PP2A	TEF2	TUB3	GAPDH	U2AF
RefFinder	TBP2	HIS3	CYP19	DNAJ	ACT2	eIF-3A	eIF-3C	TUA2	PP2A	TEF2	TUB3	GAPDH	U2AF
Different developmental stages
Delta Ct	PP2A	TUA2	HIS3	DNAJ	TUB3	CYP19	eIF-3C	ACT2	GAPDH	TEF2	TBP2	eIF-3A	U2AF
BestKeeper	DNAJ	TUA2	HIS3	ACT2	CYP19	eIF-3C	PP2A	TUB3	GAPDH	TBP2	TEF2	eIF-3A	U2AF
NormFinder	PP2A	TUA2	TUB3	HIS3	DNAJ	CYP19	eIF-3C	ACT2	GAPDH	TEF2	TBP2	eIF-3A	U2AF
GeNorm	ACT2 | HIS3		DNAJ	TUA2	PP2A	CYP19	eIF-3C	TUB3	GAPDH	TBP2	TEF2	eIF-3A	U2AF
RefFinder	TUA2	PP2A	HIS3	DNAJ	ACT2	TUB3	CYP19	eIF-3C	GAPDH	TBP2	TEF2	eIF-3A	U2AF
Different tissues
Delta Ct	PP2A	TBP2	eIF-3A	CYP19	TEF2	DNAJ	HIS3	GAPDH	eIF-3C	TUA2	U2AF	TUB3	ACT2
BestKeeper	TBP2	DNAJ	PP2A	eIF-3C	eIF-3A	CYP19	HIS3	TEF2	GAPDH	TUA2	U2AF	TUB3	ACT2
NormFinder	TBP2	CYP19	PP2A	eIF-3A	TEF2	DNAJ	eIF-3C	HIS3	GAPDH	TUA2	U2AF	TUB3	ACT2
GeNorm	eIF-3A | PP2A		CYP19	TBP2	TEF2	DNAJ	HIS3	eIF-3C	GAPDH	TUA2	U2AF	TUB3	ACT2
RefFinder	TBP2	PP2A	eIF-3A	CYP19	DNAJ	TEF2	eIF-3C	HIS3	GAPDH	TUA2	U2AF	TUB3	ACT2
Salt stress
Delta Ct	ACT2	PP2A	TBP2	TUA2	HIS3	TUB3	U2AF	CYP19	TEF2	DNAJ	eIF-3A	eIF-3C	GAPDH
BestKeeper	eIF-3C	DNAJ	U2AF	TUB3	TEF2	PP2A	TBP2	ACT2	TUA2	HIS3	CYP19	eIF-3A	GAPDH
NormFinder	ACT2	PP2A	TBP2	TUA2	HIS3	CYP19	TUB3	TEF2	U2AF	DNAJ	eIF-3A	eIF-3C	GAPDH
GeNorm	ACT2 | TBP2		PP2A	TUB3	U2AF	TEF2	DNAJ	TUA2	HIS3	CYP19	eIF-3A	eIF-3C	GAPDH
RefFinder	ACT2	TBP2	PP2A	TUB3	U2AF	TUA2	DNAJ	eIF-3C	TEF2	HIS3	CYP19	eIF-3A	GAPDH
Heat stress
Delta Ct	ACT2	TUA2	PP2A	CYP19	DNAJ	TBP2	TUB3	U2AF	HIS3	TEF2	eIF-3A	eIF-3C	GAPDH
BestKeeper	eIF-3C	TEF2	U2AF	TUB3	TBP2	PP2A	DNAJ	ACT2	TUA2	CYP19	HIS3	eIF-3A	GAPDH
NormFinder	ACT2	TUA2	PP2A	CYP19	DNAJ	TBP2	TUB3	HIS3	U2AF	TEF2	eIF-3A	eIF-3C	GAPDH
GeNorm	CYP19 | TUA2		ACT2	DNAJ	PP2A	TBP2	TUB3	U2AF	TEF2	HIS3	eIF-3C	eIF-3A	GAPDH
RefFinder	ACT2	TUA2	CYP19	PP2A	DNAJ	TBP2	TUB3	eIF-3C	U2AF	TEF2	HIS3	eIF-3A	GAPDH
Cold stress
Delta Ct	PP2A	ACT2	DNAJ	TBP2	TEF2	TUA2	CYP19	HIS3	U2AF	TUB3	eIF-3C	GAPDH	eIF-3A
BestKeeper	eIF-3C	TUB3	TBP2	DNAJ	ACT2	TUA2	TEF2	U2AF	PP2A	CYP19	HIS3	GAPDH	eIF-3A
NormFinder	PP2A	ACT2	DNAJ	HIS3	TEF2	TBP2	CYP19	U2AF	TUA2	TUB3	eIF-3C	GAPDH	eIF-3A
GeNorm	TEF2 | U2AF		PP2A	DNAJ	TBP2	ACT2	TUA2	CYP19	HIS3	TUB3	eIF-3C	GAPDH	eIF-3A
RefFinder	PP2A	ACT2	DNAJ	TEF2	TBP2	U2AF	eIF-3C	TUB3	TUA2	HIS3	CYP19	GAPDH	eIF-3A
Osmotic stress
Delta Ct	DNAJ	PP2A	ACT2	TUB3	CYP19	U2AF	TBP2	HIS3	TEF2	TUA2	eIF-3C	eIF-3A	GAPDH
BestKeeper	eIF-3C	TEF2	U2AF	DNAJ	TUB3	TBP2	PP2A	ACT2	HIS3	CYP19	TUA2	eIF-3A	GAPDH
NormFinder	ACT2	CYP19	HIS3	PP2A	TUB3	DNAJ	TBP2	TUA2	U2AF	TEF2	eIF-3C	eIF-3A	GAPDH
GeNorm	DNAJ | U2AF		TEF2	TUB3	TBP2	PP2A	ACT2	CYP19	HIS3	TUA2	eIF-3C	eIF-3A	GAPDH
RefFinder	DNAJ	U2AF	ACT2	PP2A	TUB3	TEF2	CYP19	eIF-3C	TBP2	HIS3	TUA2	eIF-3A	GAPDH
All samples
Delta Ct	PP2A	TBP2	CYP19	TEF2	DNAJ	HIS3	TUA2	TUB3	ACT2	U2AF	eIF-3C	eIF-3A	GAPDH
BestKeeper	eIF-3C	DNAJ	U2AF	TBP2	TUB3	TEF2	PP2A	ACT2	CYP19	TUA2	HIS3	eIF-3A	GAPDH
NormFinder	PP2A	CYP19	TBP2	HIS3	TEF2	TUA2	DNAJ	TUB3	ACT2	U2AF	eIF-3C	eIF-3A	GAPDH
GeNorm	PP2A | TBP2		DNAJ	TEF2	TUB3	CYP19	TUA2	HIS3	ACT2	U2AF	eIF-3C	eIF-3A	GAPDH
RefFinder	PP2A	TBP2	DNAJ	CYP19	TEF2	eIF-3C	TUB3	HIS3	TUA2	U2AF	ACT2	eIF-3A	GAPDH
Gene name	ACT2	CYP19	DNAJ	eIF-3A	eIF-3C	GAPDH	HIS3	PP2A	TBP2	TEF2	TUA2	TUB3	U2AF
Number of times the best gene appears	9	1	5	1	6	0	2	11	8	1	2	0	2

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.

3.4. Validation of the Stability of Reference Genes

In order to detect the performance of expression stability among reference genes, we selected two reference genes with a high degree of stability and one unstable gene to analyze the expression patterns of Cellulose synthase-like protein E6 gene (CSLE6) under different experimental conditions (Figure 8). The reference genes were selected for expression normalization in different experimental conditions as follows: TBP2, HIS3, and U2AF for different genotypes, TUA2, PP2A, and U2AF for different developmental stages, TBP2, PP2A, and ACT2 for different tissues, ACT2, TBP2, and GAPDH for salt stress, ACT2, TUA2, and GAPDH for heat stress, PP2A, ACT2, and eIF-3A for cold stress, DNAJ, U2AF, and GAPDH for osmotic stress. Meanwhile, the 2−ΔΔCt method was used to calculate the relative expression level of target gene CSLE6.

Figure 8

The relative expression level of target gene CSLE6 in various experimental conditions. The most two stable reference genes and the most unstable reference genes in different conditions were selected for expression normalization. (A) Different genotypes, (B) different developmental stages, (C) different tissues, (D) salt stress, (E) heat stress, (F) cold stress, and (G) osmotic stress. Different letters in the same sample represent a significant difference among three reference genes at the 0.05 level. Error bars indicate standard deviation.

As shown in Figure 8, the expression patterns of CSLE6 generated by two stable reference genes were similar. In contrast, the relative expression level of CSLE6 showed abnormal trends when the unstable reference genes were selected for expression normalization. Moreover, there was no significant difference (p < 0.05) in the relative expression level of CSLE6 between two stable reference genes in most of the sample sets, while the expression trend of CSLE6 derived from unstable reference genes was significantly different from that of stable reference genes. For example, the relative expression level of CSLE6 in HZ (Figure 8A) and tillering (Figure 8B) samples were underestimated, and abnormally up-regulated expressions were found in the stem sample (Figure 8C), the root and leaf samples under salt stress at 24 h (Figure 8D), the root sample under heat stress at 12 h (Figure 8E), the leaf sample under cold stress (Figure 8F), and osmotic stress (Figure 8G) at 24 h. Furthermore, the expression pattern analysis of CSLE6 exhibited that CSLE6 had a high expression level in the HZ genotype sample, the tillering sample, and the leaf sample with non-stress treatment. The target gene also had high expression in the leaf sample under salt stress at 24 h, and in the root sample under heat and osmotic stress at 12 h. The expression level of CSLE6 was down-regulated in leaf and root samples under cold stress.

4. 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.

5. 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.