Selection of suitable reference genes for quantitative real-time PCR gene expression analysis in Salix matsudana under different abiotic stresses

Salix matsudana is a deciduous, rapidly growing willow species commonly cultivated in China, which can tolerate drought, salt, and heavy metal stress conditions. Selection of suitable reference genes for quantitative real-time PCR is important for normalizing the expression of the key genes associated with various stresses. To validate suitable reference genes, we selected 11 candidate reference genes (five traditional housekeeping genes and six novel genes) and analyzed their expression stability in various samples, including different tissues and under different abiotic stress treatments. The expression of these genes was determined using five programs—geNorm, NormFinder, BestKeeper, ΔCt, and RefFinder. The results showed that α-TUB2 (alpha-tubulin 2) and DnaJ (chaperone protein DnaJ 49) were the most stable reference genes across all the tested samples. We measured the expression profiles of the defense response gene SmCAT (catalase) using the two most stable and one least stable reference genes in all samples of S. matsudana. The relative quantification of SmCAT varied greatly according to the different reference genes. We propose that α-TUB2 and DnaJ should be the preferred reference genes for normalization and quantification of transcript levels in future gene expression studies in willow species under various abiotic stress conditions.

Gene expression analysis has been applied to understand different kinds of biological processes 39 . Quantitative real-time polymerase chain reaction (qRT-PCR) is widely used for gene expression analysis due to its high sensitivity, accuracy, specificity, and reproducibility [40][41][42] . However, factors such as sample amount, RNA integrity, reverse transcription efficiency, and cDNA quality can significantly influence the reliability of the gene expression results [43][44][45] . To reduce the influence of these factors, internal reference genes are used to obtain accurate biologically meaningful expression values 46 ; however, unstable reference genes can cause significant biases and misinterpretations of the expression data 47,48 . Actin (ACT) and β-tubulin (β-TUB) have been used as reference genes for qRT-PCR normalization in gene expression analysis in S. matsudana under salt and copper stresses 37,49 ; however, a systematic study to validate reference genes has not been reported for S. matsudana under abiotic stresses. To obtain accurate expression data, it is necessary to select suitable reference genes for each plant species and to verify their stability under the specific experimental conditions of interest.
In this study, we determined the expression profiles of 11 candidate reference genes from S. matsudana in six different tissues and under three kinds of abiotic stresses. The 11 candidate genes were ACT, alpha-tubulin 1 (α-TUB1), alpha-tubulin 2 (α-TUB2), chaperone protein DnaJ 49 (DnaJ), E3 ubiquitin-protein ligase ARI8 (ARI8), F-box family protein (F-box), histone H2A (H2A), heat shock 70 kDa protein (HSP 70), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), membrane-anchored ubiquitin-fold protein (MUB), and β-TUB. The transcriptome data of S. matsudana were used as the source of the potential reference genes (Unpublished data). The stabilities of the 11 reference genes were analyzed using five statistical algorithms-geNorm 43 , NormFinder 44 , BestKeeper 50 , Δ Ct method 51 , and RefFinder, a web-based software 52 . The expression levels of the defense response gene SmCAT (catalase) as a target gene were assayed to verify the selected reference genes. The results will provide suitable reference genes for qRT-PCR normalization for accurate gene expression analysis in S. matsudana under different stress conditions.

Materials and Methods
Plant materials and stress treatments. Cuttings (approximately 10 cm long) from annual branches of S. matsudana were grown in hydroponics. Plants were supplemented with water containing 1/4 strength Hoagland 53 solution on alternate days under normal conditions (25 °C, 16 h light/8 h dark). After 45 days of culture, groups of S. matsudana seedlings were subjected to different abiotic stresses in solutions containing 1/4 strength Hoagland solution at pH 6.0 as follows: drought (15% PEG 6000), salt (100 mM NaCl), and heavy metal (100 μ M CdCl 2 ). Untreated seedlings were used as the control. The roots of the treated plants were sampled at 0 h, 12 h, 24 h, 48 h, and 72 h. Tissues from the root, xylem, bark, stem, leaf, and flower were collected from the untreated plants. All the samples from three biological replicates were carefully harvested, immediately frozen in liquid nitrogen, and stored at − 80 °C until total RNA extraction.
Total RNA isolation and cDNA synthesis. Total RNA from each sample was isolated from approximately 0.1 g fresh root using a total RNA kit (NORGEN, Thorold, Canada) and treated with DNase I (TaKaRa, Dalian, China) to remove any genomic DNA contamination. The RNA concentration of each sample was determined using a NanoDrop-2000 spectrophotometer (Thermo, Wilmington, USA). Samples with a 260/280 ratio of 1.9-2.1 and a 260/230 ratio ≥ 2.0 were chosen to determine the quality and purity of the RNA preparations. The integrity of the purified RNA was checked by 1.0% (p/v) agarose gel electrophoresis. Subsequently, first-strand cDNA was synthesized in a 20-μ L reaction mixture in an Invitrogen SuperScript First Strand Synthesis System (Invitrogen, Carlsbad, USA) following the manufacturer's instructions, and stored at − 20 °C until use.
Screening of candidate reference genes and primer design. We identified 11 candidate reference genes and one target gene (Table 1) from the S. matsudana transcriptome data. Primers were designed based on the sequences the 11 genes using Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/primer3/) with the following criteria: GC content 45-65%, optimal Tm 58-61 °C, primer length 18-22 bp, and amplicon length 120-220 bp ( Table 1). The specificity of each selected primer pair was observed via standard RT-PCR using Premix Ex Taq (TaKaRa, Dalian, China), and each gene was verified by 2% agarose gel electrophoresis and sequenced to ensure its reliability.
qRT-PCR. qRT-PCR amplification was performed in 96-well plates with a Applied Biosystems 7300 Real-Time PCR System (Applied Biosystems, CA, USA) using SYBR ® Premix Ex Taq ™ (TaKaRa, Dalian, China). PCR reactions were prepared in 20 μ L volumes containing: 2 μ L of 50-fold diluted synthesized cDNA, 10 μ L 2 × SYBR Premix Ex Taq ™ , 0.8 μ L of each primer, 0.4 μ L ROX reference dye (50×), and 6.8 μ L ddH 2 O. The reactions comprised an initial step of 95 °C for 30 s, followed by 40 denaturation cycles at 95 °C for 5 s and primer annealing at 60 °C for 31 s. Next, the melting curves ranging from 60 °C to 95 °C were evaluated in each reaction to check the specificity of the amplicons. Biological triplicates of all the samples were used for the qRT-PCR analysis, and three technical replicates were analyzed for each biological sample. The threshold cycle (Ct) was measured automatically.
Statistical analysis to determine the expression stability of the candidate reference genes.
Standard curves were generated in Microsoft Excel 2013 to calculate the gene-specific PCR efficiency and the correlation coefficient from 5-fold series dilution of a mixed cDNA (flower, bark, and stem) template for each primer pair. The amplification plots, melting curves and sequencing peaks were shown in Figure S1a,b,c. The PCR amplification efficiency (E) and the correlation coefficient were calculated using the slope of the standard curve according to the equation E = [5 −1/slope − 1] × 100. Stabilities of the 11 selected reference genes were evaluated by four algorithms-geNorm, NormFinder, BestKeeper, and the Δ Ct method. Finally, RefFinder (http://www. Scientific RepoRts | 7:40290 | DOI: 10.1038/srep40290 fulxie.0fees.us), a comprehensive evalution platform integrating the four above algorithms, ranked the overall stabilities of these 11 candidate genes. Pairwise variations based on the geNorm calculation were used to determine the optimal number of candidate reference genes for accurate normalization.
Expression normalization of SmCAT gene based on different reference genes. The defense response gene SmCAT was selected as the target gene to measure the stabilities of the candidate reference genes by quantifying SmCAT expression levels in all the tested samples. SmCAT gene expression levels were normalized with the two most stable candidate reference genes (α-TUB2 and DnaJ), as well as one of the least stable reference genes (β-TUB).

Results
qRT-PCR data for the candidate reference genes. The 11 selected candidate reference genes ( Table 1) are orthologs of genes in Salix purpurea, for which the whole genome has been sequenced. The specificity and accuracy of the primers designed for the selected genes were determined by 2% agarose gel electrophoresis ( Figure S2a), and further confirmed by a single peak in the melting-curve analysis ( Figure S2b). The primer sequences, amplicon length, correlation coefficient, and PCR amplification efficiency are shown in Table 1. Furthermore, the qRT-PCR products were sequenced (File S1) to determine the accuracy of the 11 genes.
To evaluate the stability of the 11 candidate reference genes at the transcript level under the three abiotic stress conditions, the gene expression levels were determined by the average Ct values, which varied from 17 to 30 (Fig. 1). According to the average Ct values of all the samples, α-TUB1 was the most abundantly expressed gene, followed by DnaJ, α-TUB2, and F-box, while H2A was the least abundantly expressed gene, followed by β-TUB, ACT and MUB.
Analysis of gene expression stability. Expression stabilities of the 11 candidate reference genes were determined using geNorm, NormFinder, Δ Ct, and BestKeeper, and their overall stabilities were ranked by RefFinder across all the stress treatments and tissue samples. geNorm analysis. The stabilities of the 11 candidate reference genes of S. matsudana calculated using geNorm were ranked in the different tissues and abiotic stress treatments according to their M values, as shown in Fig. 2. The lowest M value indicates the most stable reference gene, and the highest M value indicates the least stable one. DnaJ and ARI8 had the highest expression stabilities in the six tissues, and all the genes had M values below the threshold of 1.5 (Fig. 2a). The top two most stable genes were DnaJ and α-TUB2 for drought and heavy metal stresses, and α-TUB2 and MUB for salt stress (Fig. 2b,c,d). When the stabilities from all the samples were combined, DnaJ and α-TUB2 were determined to be the most stable reference genes in all the samples (Fig. 2e), while β-TUB had the less stability.
The pairwise variation (V n /V n+1 ) between two sequential normalization factors NF n and NF n+1 was calculated by the geNorm algorithm to determine the optimal number of reference genes for accurate normalization. A cutoff value of 0.15 is the recommended threshold indicating that an additional reference gene will make no remarkable contribution to the normalization. The V 2/3 values in the tissues, salt, and heavy metal were less than 0.15 (Fig. 3), which suggested that the top two reference genes were sufficient for accurate normalization. For the  drought stress samples V 4/5 was 0.123, indicating that the top four reference genes (DnaJ, α-TUB2, MUB, and ACT) were needed for accurate normalization. For the total samples V 3/4 was 0.138, showing that three reference genes (DnaJ, α-TUB2, and MUB) were required. Table 2, DnaJ was the most stable gene (lowest stability value) in the salt and drought subsets calculated using NormFinder. For the heavy metal samples, α-TUB2 was the most stable gene, while ARI8 was the most stable gene in the different tissues. When all samples were taken together to determine the stability of reference genes, the three most stable genes were α-TUB2, ARI8, and DnaJ. ΔCt analysis. The 11 candidate reference genes from the most to least stable expression, as calculated by the Δ Ct method, are listed in Table 3. α-TUB2 was the most stable reference gene in the drought, heavy metal, and total samples subsets. MUB and ARI8 were the most stable genes for the salt subset and different tissues, respectively, and were considered the ideal reference genes.

NormFinder analysis. As shown in
BestKeeper analysis. BestKeeper determined the stabilities of the candidate reference genes based on their standard deviation (SD). Genes with SD > 1 was considered unacceptable reference genes. The genes are listed from most to least stable in Table 4. DnaJ was the most stable gene in the tissue and drought subsets, while GAPDH and α-TUB2 were the most stable genes in the heavy metal and salt subsets.
RefFinder analysis. To acquire reliable results for the expression stabilities of the 11 candidate reference genes of S. matsudana, the rankings of the four algorithms were integrated by RefFinder and the results are shown in Table 5. The 11 genes were ranked from the most to least stable expression by RefFinder (Fig. 4). The expression of α-TUB2 was ranked the most stable under the salt and heavy metal stress treatments, and the expression of DnaJ was ranked the most stable under the drought stress treatment. Overall, the best reference gene for accurate transcript normalization in all of the samples was α-TUB2, which had the lowest Geomean (geometric mean) of the ranking values.   Reference gene validation. To validate the performance of the best ranked candidate reference genes, the expression patterns of SmCAT (catalase) were analyzed (Fig. 5). CAT as abiotic stress inducible genes, are up-regulated by drought 54 , salt 55 , and Cd 56 treatments. The CAT with low affinity towards H 2 O 2 but with a high processing rate 57 , can operate through a complex networking machinery to avoid damage caused by ROS 58 . In this study, we used the most stable reference genes (α-TUB2 and DnaJ) and the least stable gene (β-TUB) as internal controls for normalization of SmCAT according to the RefFinder rankings. The expression profiles of SmCAT were determined in different tissues and under drought, salt, and heavy metal stresses. When the stable reference genes α-TUB2 and DnaJ were used for normalization, SmCAT exhibited similar expression trends. However, when the least stable reference gene β-TUB was used for normalization, the expression patterns of SmCAT were different from that obtained using the two stable reference genes.

Discussion
Abiotic stress conditions including drought, salt, and heavy metals bring great losses to forestry productivity and increase the risk of environment. To guarantee sustainable forestry productivity and decrease the risk of environment, it is imperative to understand the regulation and function of the key genes under different abiotic stresses.
To study gene expression variations and determine gene regulation patterns, suitable reference genes are prerequisite to accurately determine the expression levels of target genes. qRT-PCR is a reliable and accurate technique for measuring gene expression levels. Some suitable reference genes under abiotic stresses, such as GAPDH 59,60 and DnaJ 10 , have been detected in plants; however, the number of reference genes evaluated is limited, especially for woody plants.
S. matsudana is an important afforestation and greening material in China that can adapt to harsh environments including drought, salt, and heavy metal. A good understanding of the molecular mechanisms related to abiotic stress responses in woody plants will not only help in improving forestry productivity but also help to decrease the risk of environment. A few studies have explored the ability of S. matsudana to withstand different abiotic stresses; however, the study of reference genes in willows has lagged behind that of other major plant species. To address this problem, we analyzed the expression of 11 candidate reference genes, five traditional reference genes (ACT, α-TUB1, α-TUB2, β-TUB, and GAPDH) and six new genes (DnaJ, ARI8, MUB, HSP70,  F-box, and H2A)   qRT-PCR methods. The best and worst candidate reference genes were further verified by expression profiling of the defense response gene SmCAT. We used five different statistical algorithms to determine the stabilities of candidate reference gene(s) under various abiotic stress conditions in S. matsudana. The results listed in Table 5 showed that, for the most parts, geNorm, NormFinder, Δ Ct, and RefFinder consistently ranked the same genes as the most stable candidate reference genes. The BestKeeper algorithm is different from the other algorithms, which explains why the BestKeeper results showed the least correlation with the others 61 . Therefore, we selected the reference gene(s) determined by geNorm, NormFinder, Δ Ct, and RefFinder.
α-TUB2 and DnaJ were the two most stable reference genes in all the sample sets according to the four algorithms. α-TUB2 encoding a cytoskeleton structure protein 62 and DnaJ encoding a cellular chaperone have the ability to repair heat-induced protein machinery damage 63,64 . Our results are in agreement with several previous studies, which showed that α-TUB2 and DnaJ were established as the most stable reference genes in plants under abiotic stresses; for example, in Syntrichia caninervis under drought, salt, and heavy metal stresses 65 , Corchorus capsularis under drought stress 10 , Buchloe dactyloides under salt stress 66 , and Platycladus orientalis under salt stress 67 . Normalization with multiple reference genes is an effective way to avoid erroneous data that may result from using a single reference gene 68 . In this study, two top ranked reference genes, DnaJ and α-TUB2 under heavy metal stress and α-TUB2 and MUB under salt stress, were appropriate for gene expression normalization, Meanwhile. Four reference genes, DnaJ, α-TUB2, MUB, and ACT under drought stress, were needed for accurate normalization. Two reference genes were found to be sufficient to analyze the expression of target genes in sorghum 62 , jute 10 , and moss 65 .
Significant differences were revealed in the expression patterns of the target gene SmCAT when was normalized with the two most stable genes (α-TUB2 and DnaJ) compared with one of the least stable genes (β-TUB) (Fig. 5), The results emphasize the importance of using stable reference genes for normalization. Our findings indicated that α-TUB2 and DnaJ either singly or in combination are suitable for normalization of gene expression in S. matsudana under different abiotic stresses. Consequently, we recommend α-TUB2 and DnaJ as the most   suitable reference genes for normalization of qRT-PCR expression data in S. matsudana under diverse abiotic stress conditions. To the best of our knowledge, this is the first report on the identification and validation of suitable reference genes for qRT-PCR analysis in S. matsudana under abiotic stresses.

Conclusion
To validate suitable reference genes for gene expression normalization in S. matsudana under drought, salt, and heavy metal stresses, we selected 11 candidate reference genes using four systematic statistical algorithms (geNorm, NormFinder, Δ Ct, and BestKeeper). The obtained results were compared and ranked using RefFinder. Based on the gene stability analysis, we identified α-TUB2 and DnaJ as the most stable reference genes for normalization of gene expression under drought, salt, and heavy metal stress conditions. Furthermore, the expression profiles of SmCAT validated α-TUB2 and DnaJ could be used as suitable reference genes. The reference genes identified in this study will facilitate accurate and consistent expression analysis of stress tolerance genes in willows and woody plants under various abiotic stress conditions for functional genomic studies.