Selection of reliable reference genes for quantitative real-time PCR gene expression analysis in Jute (Corchorus capsularis) under stress treatments

To accurately measure gene expression using quantitative reverse transcription PCR (qRT-PCR), reliable reference gene(s) are required for data normalization. Corchorus capsularis, an annual herbaceous fiber crop with predominant biodegradability and renewability, has not been investigated for the stability of reference genes with qRT-PCR. In this study, 11 candidate reference genes were selected and their expression levels were assessed using qRT-PCR. To account for the influence of experimental approach and tissue type, 22 different jute samples were selected from abiotic and biotic stress conditions as well as three different tissue types. The stability of the candidate reference genes was evaluated using geNorm, NormFinder, and BestKeeper programs, and the comprehensive rankings of gene stability were generated by aggregate analysis. For the biotic stress and NaCl stress subsets, ACT7 and RAN were suitable as stable reference genes for gene expression normalization. For the PEG stress subset, UBC, and DnaJ were sufficient for accurate normalization. For the tissues subset, four reference genes TUBβ, UBI, EF1α, and RAN were sufficient for accurate normalization. The selected genes were further validated by comparing expression profiles of WRKY15 in various samples, and two stable reference genes were recommended for accurate normalization of qRT-PCR data. Our results provide researchers with appropriate reference genes for qRT-PCR in C. capsularis, and will facilitate gene expression study under these conditions.


INTRODUCTION
With the increase of global environmental awareness, more and more people are actively purchasing goods made from ecologically friendly materials. Unfortunately, many of these materials are intrinsically unrecyclable, including many of the predominantly used polymers. However, many natural fibers such as jute and kenaf possess properties that are comparable to more traditional composites, including stiffness, flexibility, impact resistance, and elasticity (Sydenstricker et al., 2003). The overlap between natural fiber properties and those of traditional, reinforced composites are underscored by their environmentally friendly nature. For instance, fiber components are biodegradable, renewable, and result in low energy consumption. Collectively, these characteristics lend natural fibers to being logical substitutes for non-renewable synthetic fibers (Oksman et al., 2003;Corrales et al., 2007).
Jute (Corchorus capsularis L.) is an annual herbaceous fiber crop. It is found predominantly in Southeast Asia and is the second cheapest and most commercially available fiber crop, whereby it provides a biodegradable and renewable lignocellulose fiber. Jute is a promising and featured fiber material among many value-added industrial products, due in large part to its high luster, moisture absorption properties, ability for rapid water loss, and easy breakdown. Furthermore, jute has much potential for use in the industrial production of packaging materials (Sydenstricker et al., 2003). With the recent and increasing interest in jute, work has been conducted to better understand its physiological and biochemical properties (Corrales et al., 2007;del Río et al., 2009;Defoirdt et al., 2010). At a molecular level, several studies have also been done to examine relevant molecular markers (e.g., SSR, ISSR, RAPD, and AFLP) (Qi et al., 2003a,b;Basu et al., 2004;Roy et al., 2006;Mir et al., 2008Mir et al., , 2009, ESTs, stress response factors (Alam et al., 2010), and a transformation system (Chattopadhyay et al., 2011;Zhang et al., 2013). However, many of these previous experiments used single or multiple gene expression analysis via classical RT-PCR and/or Northern blot (Alam et al., 2010). Only a handful of studies have used quantitative reverse transcription PCR (qRT-PCR) to determine the expression pattern of functional jute genes (Zhang et al., 2013. As of this writing, qRT-PCR provides the most efficient, sensitive, low cost, and reproducible method for accurate and rapid detection and quantification of mRNA transcription levels for a given gene of interest (Bustin, 2002). However, several factors including RNA integrity, reverse transcription efficiency, cDNA quality, sample amount, and/or extraneous tissue and cell activities can significantly influence the accuracy of gene expression (Bustin, 2002;Huggett et al., 2005). To lessen these problems, one or more reference genes are required to account for the variance between samples and/or reactions. To this end, the transcription level of an ideal reference gene should remain constant across different tissues, treatments, and developmental stages (Gutierrez et al., 2008). A number of commonly used housekeeping genes (e.g., β-actin, elongation factor 1α, 18S ribosomal RNA, and polyubiquitin) have been used, but results indicated that their expression levels are somewhat unstable. Thus, their use as internal reference genes should be taken with some caution under a given set of experimental conditions (Gutierrez et al., 2008). Taken together, these inconsistencies highlight the need to evaluate the stability of candidate reference genes under the relevant experimental conditions prior to using for gene expression normalization by qRT-PCR.
In the present study, 11 genes were selected as candidate reference genes for evaluation expression stability in jute: 18S ribosomal RNA (18S rRNA), actin (ACT), actin 7 (ACT7), chaperone protein dnaJ (DnaJ), eukaryotic elongation factor 1-alpha (EF1α), ras-related small GTP-binding protein (RAN), alpha-tubulin (TUBα), beta-tubulin (TUBβ), ubiquitinconjugating enzyme like protein (UBC), ubiquitin extension protein (Imai et al., 2014) and ubiquitin (UBQ). We then sought to reveal which reference gene(s) were the best option for gene expression qRT-PCR analysis of jute under various experimental treatments using three available statistical algorithms, geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), and BestKeeper (Pfaffl et al., 2004). Finally, we validated the expression levels of the transcription factor CcWRKY15 using the selected reference gene(s). The results from this work will facilitate future studies on gene expression as well as foster a better understanding of how novel genes function in the molecular mechanisms of jute biological and/or physiological processes.

Plant Materials and Treatments
Jute (Corchorus capsularis L.) variety Huangma 108 was used for all experimental treatment groups. To ensure disease-free materials, seeds were rinsed under running water for 10 min and sterilized with 5% sodium hypochlorite for 10 min. They were then washed three times with sterile water before germinated on filter paper that had been saturated with water in complete darkness at 28 • C. After 3 days, seedlings were grown in the greenhouse in 1/4 Hoagland solution under a 16/8-h light/dark cycle at 30/26 • C (day/night). Seedlings were assessed at the 3-5 leaf stage and the most consistent were used for (i) the abiotic and biotic stress groups or for (ii) harvesting of different plant tissues (e.g., root, stem, and leaf).
For the salinity and drought treatments, seedlings were subjected to 200 mM sodium chloride and 15% (w/v) PEG6000, respectively, and harvested at 0, 2, 4, 6, 8, 12, and 24 h. For the biotic treatment, seedlings were inoculated with 2 ml (10 6 spores per ml) of Colletotrichum siamense spores in suspension previously described by Ma et al. (2009) and then sampled at 0, 1, 6, 12, 24, 48, 72, and 96 h. Tissues from the roots, stems, and leaves were collected from plants in the 3-5 leaf stage that had grown well under the experimental conditions. All samples were harvested from three replicate plants, giving a total of 25 samples comprised of 14 abiotic and eight biotic stress treatment samples and three tissue-specific samples. Samples were immediately frozen in liquid nitrogen and stored at −80 • C until RNA extraction.

Total RNA Extraction and cDNA Synthesis
Total RNA was extracted from approximately 100 mg fresh leaf tissue using the OMEGA isolation kit (R6827-01, USA) according to the manufacturer's instructions. The genomic DNA contamination was eliminated using RNase-free DNase I (TaKaRa, Japan) and RNA sample quality was then determined using the NanoDrop 2000 spectrophotometer (NanoDrop, Thermo Scientific). RNA integrity was assessed by 2.0% agarose gel electrophoresis. Finally, RNA samples with an A 260 /A 280 ratio between 1.9 and 2.1 and an A 260 /A 230 ratio greater than 2.0 were used for further analysis. Subsequently, the first-strand cDNA was synthesized from 1 µg total RNA in a volume of 20 µl using the PrimeScript R RT reagent kit (TaKaRa, Japan) by following the manufacturer's protocol. All final cDNA samples were diluted 10-fold for subsequent quantitative real-time RT-PCR reactions. Samples were stored at −20 • C until use.
To check the specificity of the amplicon, all primer pairs were initially tested via standard RT-PCR using the Premix Ex Taq (TaKaRa, Japan) and the amplification product of each gene was verified by electrophoresis on a 2% agarose gel. Realtime amplification reactions were performed with the Applied Biosystems 7500 Real-Time PCR System using SYBR R Premix Ex Taq ™ (TaKaRa, Japan). Reactions were prepared in a 20 µl volume containing the following: 2 µl cDNA template, 0.4 µl of each amplification primer, 0.4 µl ROX Reference Dye II, 10 µl 2x SYBR Premix Ex Taq ™ , and 6.8 µl dH 2 O. Amplifications were performed with an initial 30 s step of 95 • C followed by 40 denaturation cycles at 95 • C for 5 s and primer annealing at 60 • C for 34 s. The melting curve ranged from 60 to 95 • C and temperature was increased in increments of 0.2 • C every 10 s for all PCR products. ABI Prism Dissociation Curve Analysis Software was used to confirm the occurrence of specific amplification peaks. All qRT-PCR reactions were carried out in triplicate with template-free negative controls being performed in parallel.

Statistical Analyses of Gene Expression Stability
To select a suitable reference gene, three publicly available software packages, geNorm (version 3.5), NormFinder and BestKeeper, were used to analyze the stability of each reference gene. All analyses using these packages occurred according to the manufacturers' instructions. For geNorm and NormFinder algorithms, the raw Ct values from the qRT-PCR were transformed into the relative expression levels using the following formula: E − Ct , where Ct equaled each corresponding Ct value minus the minimum Ct value. Then, the relative expression values were imported into geNorm and NormFinder to analysis gene expression stability. For BestKeeper analysis, the Ct value was used as input data to calculate the coefficient of variation (CV) and the standard deviation (SD). The comprehensive rankings of the best reference genes were obtained by integrating the results of three algorithms. To validate the reliability of the selected reference genes, we analyzed the relative expression levels of the transcription factor CcWRKY15 in all tested samples. Additionally, a standard curve was generated from a 10-fold dilution of cDNA in a qRT-PCR assay using Microsoft Excel 2003. The PCR efficiency (E) and the regression coefficient (R 2 ) were calculated using the slope of the standard curve according to the equation E = [10 −(1/slope) − 1] × 100%. All other multiple comparisons were performed using the statistical analysis software SPSS 22.0 (SPSS Inc., USA).

qRT-PCR Data of Candidate Reference Genes
A total of 11 candidate reference genes, including 18S rRNA, ACT, ACT7, DnaJ, EF1α, RAN, TUBα, TUBβ, UBC, UBI, and UBQ, were identified and assessed under abiotic stress (NaCl stress and drought stress), biotic stress (Colletotrichum siamense) and in different tissues in this study. For each gene, the specificity of the designed primers was verified using agarose gel electrophoresis and the subsequent presence of a single band with the expected size ( Figure S1), and further confirmed by the presence of a single peak in the melting curve analysis, which was done prior to performing qRT-PCR ( Figure S2). As described in Table 1, the amplicon size ranged from 130 to 219 bp. The PCR efficiency (E) was greater than 90% and varied from 90.1% (18S rRNA) to 113.7% (ACT), and the regression coefficient (R 2 ) ranged from 0.988 (TUBα) to 0.999 (DnaJ and TUBβ) ( Table 1).
To evaluate stability of the reference genes across all experimental samples, the transcript abundances of the 11 candidate reference genes were detected by their mean Ct values. The Ct values of these candidates varied from 15.22 to 31.28, with the majority falling between 20.13 and 24.40 (Figure 1). Across all samples, 18S rRNA was the most abundantly expressed gene, with the lowest average Ct ± SD (16.85 ± 0.78), followed by EF1α (18.28 ± 0.94), UBQ (20.13 ± 1.54), UBC (21.20 ± 1.13), UBI (21.44 ± 0.62), DnaJ (22.22 ± 0.75), ACT7 (22.24 ± 0.92), and TUBβ (22.40 ± 1.07). TUBα was found to have the lowest level of expression of any of the genes tested, with a mean Ct ± SD of 28.43 ± 2.15, followed by RAN (26.37 ± 1.29) and ACT (25.75 ± 1.31) ( Table 1). Small co-variance (CV) of the Ct value indicates that a given gene is more stably expressed. Among these 11 candidate reference genes, UBQ showed the greatest variation with CV value of 7.65%, whereas both UBI (2.89%) and DnaJ (3.38%) showed the least variation in their expression levels across all tested samples. The ranking of gene stability by CV was as follows: (Table 1). Collectively, these results indicate that the transcript levels of the candidate reference genes varied across different experimental samples. Thus, it is essential to screen the most appropriate reference genes in jute in order to normalize gene expression analysis.

Stability Analysis of Reference Genes by geNorm
A geNorm-based analysis was carried out to determine which candidate reference gene(s) would be optimal in each of the tested samples sets. As shown in Figure 2, genes were ranked according to their M values. Since a lower M value indicates increased stability, RAN, and ACT7 were determined to be the most stable reference genes in total samples. In contrast, TUBα and ACT were the least stable reference genes. For each subset of the treatment, the top two reference genes for qRT-PCR normalization were ACT7 and RAN in biotic stress subset, TUBβ, and UBI in different tissues, EF1α and RAN for NaCl stress, and FIGURE 2 | Expression stability of 11 candidate genes in jute as calculated by geNorm. Mean expression stability (M) was calculated following stepwise exclusion of the least stable gene in biotic stress samples, tissue samples, NaCl-and PEG-treated samples and all samples. The least stable genes are on the left and the most stable genes on the right.
Frontiers in Plant Science | www.frontiersin.org EF1α and UBC for drought stress. In general, the most stable genes across all experimental samples were RAN and ACT7, while TUBα, ACT, and 18S rRNA were the least stable (Figure 2).
The pairwise variation (V n ) between normalization factors (NF n ) calculated by the geNorm algorithm also determines the optimal number of reference genes for accurate normalization. A cut-off value of V n/n+1 < 0.15 (Vandesompele et al., 2002) indicates that an additional reference gene makes no significant contribution to the normalization. As depicted in Figure 3, the V 2/3 values in biotic stress, NaCl stress, and PEG stress subsets were below 0.15, which indicated that two reference genes (ACT7 and RAN for biotic subset; EF1α and RAN for NaCl subset; and EF1α and UBC for drought subset) were sufficient for accurate normalization. In the tissue subset, four reference genes (TUBβ, UBI, EF1α, and RAN) were needed for accurate normalization, as the V 4/5 value was lower than 0.15. When the total samples were taken into account, the V3/4 value (0.132) was lower than the cutoff value of 0.15, which indicated that three genes (ACT7, RAN, and DnaJ) were suitable for all samples in this study (Figure 3).

Stability Analysis of Reference Genes by Normfinder
The NormFinder approach was used to determine the stability of reference genes based on inter-and intra-group variance in expression. As shown in Table 2, similar results were generated by NormFinder, which predicted ACT7 and RAN to be the two most stably expressed normalization factors in both biotic and total subsets. For the NaCl and PEG treated samples, ACT7 and RAN performed as the best reference genes by the NormFinder analysis, while ACT7 was ranked the second in NaCl subset, and ACT7 and RAN the fifth and third in PEG subset by the geNorm analysis. For the tissue subset, the four most stable genes, TUBβ, EF1α, RAN, and UBI determined by NormFinder were ranked the first, second, third and first by geNorm algorithm, respectively (Table 2; Figure 2). FIGURE 3 | Determination of the optimal number of reference genes for normalization by pairwise variation (V) using geNorm. The pairwise variation (V n /V n+1 ) was calculated between normalization factors NF n and NF n+1 by geNorm to determine the optimal number of reference genes for qRT-PCR data normalization.

Stability Analysis of Reference Genes by Bestkeeper
The BestKeeper program was used to evaluate the stabilities of reference genes based on the coefficient of variance (CV) and the standard deviation (SD) of the average Ct values. The most stable genes were identified as those that exhibit the lowest CV and SD (CV ± SD), and genes with SD greater than 1 were considered unacceptable and should be excluded (Chang et al., 2012;Xiao et al., 2014). In the biotic stress subset, 18S rRNA (2.96 ± 0.51) and UBI (2.39 ± 0.52) had lowest CV ± SD values, and showed remarkably stable expression. In the NaCland PEG-treated subset, TUBα had the lowest CV ± SD values of 0.79 ± 0.24 and 0.96 ± 0.29, respectively, and showed the most stable expression. In the total samples subset, UBI (2.32 ± 0.50) and 18S rRNA (3.76 ± 0.63) were identified as the two best reference genes for normalization. These results are inconsistent with those acquired from the geNorm and NormFinder methods (Figure 2; Tables 2, 3). In the tissue samples subset, the most stable reference genes identified by BestKeeper were UBI (1.98 ± 0.42) and TUBβ (2.03 ± 0.43). This result is consistent with the results obtained from geNorm and NormFinder analyses (Figure 2; Tables 2, 3).

Comprehensive Stability Analysis of Reference Genes
To acquire a consensus result of the best reference genes, three algorithms rankings of the stability were integrated, generating a comprehensive ranking according to the geometric mean of three rankings (Xiao et al., 2014). The comprehensive rankings were shown in Table 4: RAN and ACT7 were ranked as the top two stable reference genes in the biotic stress subset, NaCl stress subset and total samples subset; UBC and DnaJ were the two most stable genes in the PEG stress subset; TUBβ and UBI were the most stable genes, followed by EF1α and RAN in the tissue subset ( Table 4). Taken the number of reference genes to use suggested by geNorm and the comprehensive rankings into consideration, the most stable and least stable combination of reference genes in each subset was shown in Table 5.

Reference Genes Validation
To validate the selected reference genes, the relative expression levels of the target gene, CcWRKY15 under different experimental conditions were evaluated using qRT-PCR.
CcWRKY15 is a homolog of AtWRKY15, which is known to be a central regulator in the response to oxidative stress and pathogenic infection (Vanderauwera et al., 2012). Thus, it would be expected to have similar expression patterns as AtWRKY15 under abiotic and biotic stress conditions. However, a substantial divergence can occur in its relative transcript abundance when normalized to different kinds of reference genes. We therefore used the most stable reference genes found in each subset (ACT7 and RAN for Biotic stress and NaCl stress; UBC and DnaJ for PEG stress) either singly or in combination, and the least stable reference gene (ACT or TUBα), to perform a qRT-PCR analysis. Results showed that in accordance with the behavior of AtWRKY15 previously described in A. thaliana, the CcWRKY15 functions as a negative regulator and was induced by salt-stress, oxidative-stress, and pathogenic infections (C. siamense) in jute (Figure 4). In addition, we also examined the reference gene (EF1α and UBC for abiotic stress; TUBβ for fungal stress) selected by Ferdou et al. (2015) in NaCl stress, PEG stress, and fungal stress subsets. The results showed that EF1α could sever as a stable reference gene for normalization but UBC unstable under NaCl stress condition ( Figure 4A). For the PEG stress subset, the relative expression folds of CcWRKY15 normalized by EF1α were slightly decreased compared to the stable genes UBC and DnaJ (Figure 4B). For the biotic stress subset, reference gene TUBβ was not as stable as that described by Ferdou et al. (2015), conversely, similar to the expression pattern of the worst gene TUBα (Figure 4C). Our tissue type analysis revealed that the transcript abundance of CcWRKY15 was the highest in the leaf, followed by the stem and then the root (p < 0.01) ( Figure 4D). Contrastingly, when the least stable reference gene was used as normalization factor, the expression level of CcWRKY15 was significantly overestimated (p < 0.01). For example, the relative expression folds of CcWRKY15 were approximately eight-fold (p < 0.01) higher than that of ACT7, RAN, or their combination at 4 h   under NaCl stress, when ACT was used as normalization factor ( Figure 4A).

DISCUSSION
There has been a surge in interest for environment-friendly materials and jute fiber-used as reinforcement componenthas been widely used in the textile, papermaking, automotive, and aerospace industries. Given its popularity, the need for the application of new genomic tools has become increasingly important (Sydenstricker et al., 2003;Corrales et al., 2007). New technologies such as qRT-PCR now make it possible to understand the molecular mechanisms underpinning the commercially important traits of jute. qRT-PCR is an essential tool that can be used in studies of target gene expression patterns. Although past studies have used 18S rRNA as a reference gene for examining the functions of UDP-glucose pyrophosphorylase and caffeoyl-CoA 3-O-methyltransferase in jute (C. capsularis) (Zhang et al., 2013, and in the most recent study were 7 reference genes tested under abiotic and biotic stress in the other jute species Corchorus olitorius (Ferdou et al., 2015), little information is available on the systematic exploration and validation of a set of suitable reference genes in jute (C. capsularis). In addition, it has been reported that ribosomal RNA genes are not adequate reference genes due to their high transcription abundance. Ultimately, this could lead to experimental error when normalizing genes with weak expression (Jain et al., 2006;Niu et al., 2014). Comparing our results to that of Ferdou et al. (2015), the stability values of three out of seven reference genes (EF1α and UBC for abiotic stress; TUBβ for biotic stress) were significantly lower in our experiment. For instance, when UBC was used as reference gene, CcWRKY15 had a 5.0-fold higher expression value at 24 h under NaCl stress condition ( Figure 4A); when TUBβ as an internal reference, CcWRKY15 had a 6.0 times higher at 24 h after inoculation ( Figure 4C). Additionally, TUBβ was found to be one of the least stable genes in our study. These expression stability differences might be a result of different species/cultivars analyzed in the compared experiment conditions. It is consistent with the previous studies that the stability of reference genes is not only cultivar/species specific but may also be tissue specific and influenced by the different experimental treatments (Nicot et al., 2005;Štajner et al., 2013;Zhuang et al., 2015). This is also the primary reason for validating the stability of reference genes for a specific genotype and/or experimental treatment.
In this study, we used three publicly available programs, geNorm (Vandesompele et al., 2002), NormFinder (Andersen et al., 2004), and BestKeeper (Pfaffl et al., 2004), to evaluate the expression stability of 11 candidate reference genes in a total of 25 jute samples taken from different tissues and experimental treatment groups. We found different rankings for the selected genes after comparison to the ranking of the candidates generated by the three algorithms (Figure 2; Tables 2, 3). This apparent divergence probably reflects the discrepancies in the three statistical algorithms to calculate stability. NormFinder takes the inter-and intra-group variations into account, and combines them into a stability value, and finally ranks the top genes with minimal inter-and intra-group variation. In contrast, geNorm identifies two reference genes with the highest degree of similarity in expression profile and the lowest intra-group variation (Andersen et al., 2004;Jian et al., 2008;Cruz et al., 2009;Liu et al., 2012). As for BestKeeper, this program determines the stability ranking of the reference genes based on the coefficient of variance (CV) and the standard deviation (SD) values. The most stable genes are identified as those that exhibit the lowest CV ± SD values (Chang et al., 2012;Xiao et al., 2014). Similar methods have been used in previous studies of different species, such as Salicornia europaea Xiao et al., 2014), Populus euphratica , and Cynodon dactylon . Thus, referring to the previous studies (Štajner et al., 2013;Xiao et al., 2014), the integrated results were obtained from three programs, leading to a more comprehensive ranking and better accuracy for each gene (Table 4).
Of the top three reference genes defined by three algorithms in the total samples subset, ACT7 and RAN were found to be the best candidates for the biotic stress and NaCl stress subsets. The strong performance of ACT7 in jute (C. capsularis) was consistent with the results obtained in the developmental stage series of G. max (Jian et al., 2008), during grape berry development of V. vinifera (Reid et al., 2006), and across different tissues and cold-treated samples of Platycladus orientalis (Chang et al., 2012). However, this gene performed poorly in studies of O. sativa (Jain et al., 2006), S. officinarum (Ling et al., 2014), and S. tuberosum (Nicot et al., 2005), suggesting that the expression levels of reference genes are variable among different species. RAN was the other most stable reference gene found in the present study. It has also been shown to be the best performer across different tissue and hormone treatment samples of M. acuminata  as well as in the leaf or plant growth regulator treatment samples of C. melo (Kong et al., 2014). DnaJ showed relatively stable expression and ranked third across all samples (Figure 2; Tables 3, 4). This gene has also been recommended as the best reference gene in different tissues and NaCl-treated samples of P. orientalis (Chang et al., 2012).
However, the data presented show that the candidate reference genes EF1α, UBC, UBI, UBQ, and TUBβ were moderately expressed and had variable rankings in this study. For example, EF1α was ranked first by geNorm analysis and was stably expressed in NaCl-and PEG-treated samples, but was shown to be the fourth in NaCl stress subset analyzed by NormFinder; the fifth in PEG stress subset by NormFinder and Comprehensive ranking. Previously, EF1α showed stable expression during abiotic and biotic stress in O. sativa (Jain et al., 2006), L. chinensis , S. tuberosum (Nicot et al., 2005), and C. sativus (Wan et al., 2010), but was shown as an unsatisfactory reference gene in T. aestivum (Paolacci et al., 2009), G. max (Jian et al., 2008), N. benthamiana , and P. orientalis (Chang et al., 2012). UBC, an ubiquitin-conjugating enzyme gene, was ranked first in the PEG stress subset, which was also the optimal reference gene under NaCl-and PEG-treated conditions in P. orientalis (Chang et al., 2012) and C. olitorius (Ferdou et al., 2015); however, it was not the best choice for normalization in the different tissues of bamboo (Fan et al., 2013). UBQ were relatively weakly expressed across all the experimental samples according to the three programs and comprehensive analyses, but were the optimal reference genes for developmental stages and under NaCl-and PEG-treated conditions in P. orientalis (Chang et al., 2012). As previously noted, TUBβ was the most stable reference gene across various developmental stages of G. max (Jian et al., 2008), leaf senescence system of H. annuus (Fernandez et al., 2011), and among different tissues and PEGtreated samples of P. orientalis (Chang et al., 2012). In this study, similarly, all of the algorithms ranked TUBβ in the first position, indicating that TUBβ was the optimal choice as an internal control gene for different tissues investigated in jute (C. capsularis). UBI, another stable reference gene, showed relatively small variation in tissues subset, which was also the most stable reference gene in different tissues of C. sativus (Wan et al., 2010). In different tissues subset, in addition to expression of TUBβ and UBI, the expression of EF1α and RAN was also stable and therefore they were considered as the suitable reference genes (Figure 2; Table 4).
Interestingly, the commonly used reference genes TUBα, 18S rRNA, and ACT performed poorly and were not suitable for most of the experimental conditions. Several studies have shown similar results. For example, TUBα was found to be unstable as reference gene in the developmental stage in tomato (Expósito-Rodríguez et al., 2008) and in different flax tissue (Huis et al., 2010). 18S rRNA was considered as the least reliable reference gene in viral-infected N. benthamiana  and under different conditions in C. sativus (Wan et al., 2010). ACT was proved unsuitable for normalization in different flax tissues (Huis et al., 2010), during abiotic and biotic stress in S. tuberosum (Nicot et al., 2005), and in viral-infected N. benthamiana . Taken together, these findings indicate that large numbers of experimental data on gene expression should be acquired to investigate the transcript stability of commonly used reference genes under different experimental conditions.
To further validate the feasibility of the reference genes screened in this study, we analyzed the transcription profiles of the WRKY domain gene CcWRKY15, a homolog AtWRKY15 of A. thaliana. AtWRKY15 has been shown to be a key regulator of plant growth and salt/osmotic stress responses in A. thaliana (Vanderauwera et al., 2012). In this study, the expression of CcWRKY15 was normalized using the most stable reference genes in each subset both singly and combined as well as a least stable gene as an internal control (Figure 4). Our results showed that expression of CcWRKY15 was negatively induced by NaCland PEG-treated stress (Figures 4A,B) and was significantly increased after 24 h of inoculation treatment ( Figure 4C) (p < 0.01). By comparing the expression pattern of CcWRKY15 with that reported in A. thaliana (Vanderauwera et al., 2012), a supported result was found. Therefore, the results obtained from this study are credible. Moreover, we compared the results from the jute species C. olitorius selected by Ferdou et al. (2015) under the same conditions, the results showed great differences between them. These results underscore the fact that inappropriate utilization of reference genes without validation may generate bias in the analysis and lead to misinterpretation of qRT-PCR data.

CONCLUSION
We present here a systematic attempt to validate a set of candidate reference genes for the normalization of gene expression using qRT-PCR in jute (C. capsularis) under abiotic (salt and drought) and biotic (Colletotrichum siamense) stress conditions as well as across different tissue types. The expression stability of the 11 candidates was analyzed by the three applications (geNorm, NormFinder, and BestKeeper), and their results were furtherly integrated into a comprehensive ranking based on the geometric mean. For gene expression study under biotic stress and NaCl stress, we recommend ACT7 and RAN to normalize the qRT-PCR data. For gene expression study under PEG stress, UBC, and DnaJ are the two most suitable reference genes in C. capsularis. For the study of gene expression in the different tissues, TUBβ, UBI, EF1α, and RAN are recommended as the best reference genes for normalization. In addition, the two least stable reference genes 18S rRNA and ACT should be carefully used for normalization. Furthermore, the feasibility of the reference genes screened was further confirmed by comparing the expression pattern between CcWRKY15 and AtWRKT15, and the selected reference genes can perform significantly better in gene expression normalization. In particular, the reference genes selected in current study will facilitate the future work on gene expression studies in C. capsularis.