Identification and validation of reference genes for qRT-PCR analysis in mulberry (Morus alba L.)

Mulberry (Morus alba L.) is an important economic tree species in many countries. Quantitative real time PCR (qRT-PCR) has become a widely used method for gene expression studies in plants. A suitable reference gene is essential to ensure accurate and reliable results for qRT-PCR analyses. However, no reports describing the selection of reference genes have been published for mulberry. In this work, we evaluated the stability of twenty candidate reference genes in different plant tissues and under different stress conditions by qRT-PCR in mulberry using algorithms in two programs—geNorm and NormFinder. The results revealed that TUB2, UBI4, ACTIN3 and RPL4 were ranked as the most stable reference genes in the samples subsets, whereas EF1α4 and TUB3showed the least stability with both algorithms. To further validate the stability of the reference genes, the expression patterns of six genes of mulberry were analyzed by normalization with the selected reference genes. Our study will benefit future analyses of gene expression in mulberry.


Introduction
Mulberry (Morus alba L.) is widely cultivated as an economic tree species in most developing countries in Asia, as it is the sole food source of the domesticated silkworm (Bombyx mori L.). In addition, mulberry has been used as a food product and in herbal medicines for a long time, as it can be used as a pharmaceutical treatment for hyperglycemia, inflammation, hypertension, and fever [1,2]. The mulberry genome has recently been sequenced, and high-throughput sequencing has been used in several mulberry studies [3][4][5]. Genome sequence and large-scale transcriptome data have greatly facilitated molecular studies in mulberry.
Gene expression analysis is widely used in many fields of plant biological research including in development and in responses to stress. Quantitative real-time polymerase chain reaction (qRT-PCR) has become the most popular approach for gene expression studies because of its rapidity, sensitivity, and specificity [6]. However, the reliability of the results obtained from qRT-PCR depends on accurate normalization using stably expressed reference genes [7]. Thus, it is necessary to select appropriate reference genes with minimal variability during the experimental conditions before qRT-PCR analysis. In plants, some reference genes have been reported, but the stability of reference genes in different plant species is inconsistent [8][9][10]. TUB-B, TUB-A, and UBC are the most stable reference genes in celery [9]. NTB and TIP41 are most stable in bamboo [11]. In addition, the expression of reference genes varies under different experimental treatments. In a study of switchgrass, eEF-1α and ACT12 were most stable across all tissue samples, whereas CYP5 and FTHS4 were the most suitable reference genes during seed development under abiotic stress conditions as determined by NormFinder [10]. To date, a systematic study validating reference genes has not been reported in mulberry. Hence, it is necessary to identify the suitable reference genes in various tissues and under different abiotic stress conditions, which will be helpful for the accurate and reliable analysis of gene expression in this plant.
In this study, we examined the stability of twenty candidate reference genes (three ACTIN genes, three TUB genes, three UBI genes, three EF1α genes, two GAPDH genes, two MDH genes, two RPL genes, one CYP gene, and one PP2A gene) for use in qRT-PCR expression studies in mulberry. The expression patterns of these reference genes were quantified in various tissues and under different abiotic stress conditions using both geNorm and NormFinder algorithms [12,13]. In addition, three chalcone synthase genes of mulberry (MaCHS5, MaCHS6, MaCHS7) and three plant hormone related genes (MaERF, MaDELLA, MaJAZ) were used to evaluate the reliability of the most stable reference genes identified in this work.

Amplification specificity and efficiency for each candidate reference gene
To identify suitable reference genes for mulberry, 20 candidate reference genes were selected based on previous reports of control genes used in various plant species [10,14,15]. The gene sequences of ACTIN2 and ACTIN3 were obtained from the National Center for Biotechnology Information (NCBI, USA). The sequences of the eighteen other candidate reference genes were selected from our transcriptome database that was generated via high-throughput Illumina sequencing (S1 Table) [4], and the expression of these eighteen candidate reference genes was relatively stable at all treatment times of three mulberry varieties after infection with Ralstonia solanacearum in our previous study (S2 Table) [16]. The products of these twenty genes are associated with a wide variety of biological functions.
To check the specificity of the primers for these candidate reference genes, we performed qRT-PCR melting curve and agarose gel electrophoresis analyses. The product of each primer pair was a single peak in the melting curve and a single band with the expected size after agarose gel electrophoresis and polyacrylamide gelelectrophoresis (PAGE) electrophoresis (Fig 1  and S1 Fig). The qRT-PCR amplification efficiency for the fifteen candidate reference genes varied from 95.3% (ACTIN3) to 105.7% (EF1α1), and correlation coefficient (R 2 ) values ranged from 0.990 to 0.999 (Table 1).

Expression profiles of the candidate reference genes
Expression levels of the 20 candidate reference genes were measured in the 10 different samples with three biological repeats obtained from leaves, stems, roots, green fruit, and red fruit of each plant under normal growth conditions and from leaves from plants under different abiotic stress conditions. The cycle threshold (Ct) values in the qRT-PCR reactions were used to identify the differences in transcript expression levels, with lower Ct values indicating higher transcript abundance and vice versa. We investigated the Ct values of all candidate reference genes within the set of samples, and variation in the expression of these genes is shown in Fig  2. A wide range of expression levels of these genes in all tested samples was observed. Most of the candidate reference genes had Ct values that range from 18 to 24. EF1α1 had the highest mean Ct value at 24.7 and thus had the lowest level of expression among the candidate   (Fig 2). None of the candidate reference genes showed a constant expression level among all samples. Therefore, it was necessary to carry out further analysis to select the most suitable reference genes for normalizing gene expression under particular experimental conditions.

Expression stability analysis
Two commonly used analysis programs, geNorm and NormFinder, were applied to evaluate the expression stability of the candidate reference genes [12,13]. We used this software to analyze gene expression stability across different subsets of samples: 1) all samples, 2) tissue samples, 3) abiotic stress samples, and 4) green fruit and red fruit samples. The ranks of the candidate reference genes are presented in Fig 3 according to their average expression stability value (M) using geNorm. When all the samples were taken together, TUB2 and UBI4 showed the lowest average expression stability value (M = 0.531), and EF1α4 showed the highest value (M = 1.286), indicating that TUB2 and UBI4 had the most stable expression and EF1α4 had the highest level of expression variation among the 20 candidate genes. Pairwise variation was also calculated using geNorm (Fig 4). Pairwise variation analysis showed the optimal number of reference genes required for normalization, with a cut-off value of 0.15 being widely accepted as the criterion for determining a suitable number of reference genes [12]. For all samples, the value of V4/5 was 0.165 and of V5/6 was 0.128. The value of V5/6 was lower than the cut-off value of 0.15, indicating that the use of the five most stable reference genes was required for accurate normalization. Meanwhile, stability of expression for candidate reference genes was evaluated by NormFinder software according to the intra-and intergroup variations for normalization factor calculations. As shown in Table 2, TUB2 and UBI5 with the lowest stability value of 0.019 were identified as the best combination by NormFinder, and the most stable gene was TUB2 (V = 0.027), whereas EF1α4, with the highest stability value of 0.068, was identified as the least stable reference gene.
In the tissue samples subset, which included root, leaf, stem, green fruit, and red fruit of mulberry under normal growth conditions, ACTIN4 and RPL4 were identified as the best pair with an M value of 0.243 by geNorm, whereas TUB3 was the worst gene (Fig 3). Only two reference genes were sufficient for normalization because the pairwise variation value V2/3 was lower than 0.15 for this subset (Fig 4). NormFinder identified ACTIN3 and TUB2 as the best pair for the tissue samples subset, whereas ACTIN3 and EF1α4 were identified as the best and the worst reference genes, respectively ( Table 2).
For abiotic stress samples, TUB2 and GAPDH1 were the best candidates (M = 0.151) for normalization by geNorm, whereas TUB4 was the worst gene (M = 0.881) (Fig 3). Similarly, two control genes were satisfactory for normalization with a V2/3 value of 0.104 (Fig 4). UBI4 and RPL4 were identified as the most suitable pair by NormFinder, whereas TUB2 and EF1α4 were identified as the best and the worst reference genes, respectively (Table 2).
In the green fruit and red fruit sample subset, ACTIN4 and MDH1 were identified as the most stable genes by geNorm (Fig 3), whereas TUB2 and EF1α1 were identified as the most stable genes by NormFinder (Table 2). Two reference genes were sufficient for normalization because the pairwise variation value V2/3 was 0.072 for this subset by geNorm (Fig 4). Both algorithms identified EF1α4 as the least stable gene.

Reference gene validation in gene expression study
To evaluate the reliability of the reference genes identified by geNorm and NormFinder, the relative expression patterns of three chalcone synthase genes (MaCHS5, MaCHS6, and MaCHS7) in different tissues samples and three plant hormone related genes (MaERF, MaDELLA, and MaJAZ) under different abiotic stress samples were compared by qRT-PCR using different combinations of reference genes for normalization (Fig 5). Mulberry is rich in secondary metabolism flavonoids, and chalcone synthase is one of the key enzymes involved in flavonoid biosynthesis [17]. MaERF, MaDELLA, and MaJAZ were the crucial transcription factors of plant hormone signaling (ethylene, gibberellins, and jasmonic acid), which play important roles in plant responses to a wide range of abiotic stresses [18]. We investigated the relative gene expression of three MaCHS genes in different mulberry tissues samples. Four different combinations of reference genes were used for the validation, including the most stable reference gene ACTIN3, two most stably pairs (ACTIN3/TUB2 and RPL4/ACT4), and the least stably expressed gene EF1α4 in the tissue samples subset (Fig 3 and Table 2). MaCHS5 was highly expressed specifically in green fruit. Considering the expression of root samples as 1, the relative gene expression of MaCHS5 in other four tissues samples was found to be unbiased when three most stable groups (ACTIN3, ACTIN3/TUB2 and RPL4/ ACT4) were used as reference genes for normalization. However, when the least stably expressed gene EF1α4 was used as the reference gene, the relative gene expression of MaCHS5 in three tissues (stem, leave, and green fruit) was about two times comparing with which used three most stable groups (ACTIN3, ACTIN3/TUB2 and RPL4/ACT4) as reference genes ( Fig  5A). MaCHS6 and MaCHS7 were specifically and highly expressed in root and red fruit, respectively. Similar to the MaCHS5 expression trend, the relative expression of MaCHS6 and MaCHS7 were also markedly higher than other reference combinations for three tissues samples (stem, leave, and green fruit) when normalized by EF1α4 (Fig 5A). Based on the reference genes selected for normalization, the expression of target genes can show greater variability. We noted that using the three best genes for normalization-ACTIN3, TUB2, and RPL4-produced similar results to those from the best gene pairs (Fig 5A). Thus, two reference genes were sufficient for normalization in different tissues samples, which further confirmed the conclusion of the pairwise variation analysis by geNorm. Meanwhile, we surveyed relative transcript accumulation of three plant hormone related genes (MaERF, MaDELLA, and MaJAZ) under different abiotic stress samples using different combinations of reference genes. Four reference gene groups were used for the validation, including the most stable reference gene TUB2, two most stably pairs (ACTIN3/TUB2 and GADPH2/MDH2), and one of the least stably expressed gene EF1α4 in abiotic stress samples subset (Fig 3 and Table 2). Considering the expression of drought stress samples as 1, the expression profiles of the three target genes in other four abiotic stress samples showed similar trends when stable reference gene groups (TUB2, TUB2/GAPDH1 and UBI4/RPL4) were used for normalization. When the least stably expressed gene EF1α4 was used as an internal gene, the expression level of the target genes were showed observably lower than other reference combinations for three abiotic stress samples (salt, heat, and wound) (Fig 5B). These results validated the effects of reference genes for normalization in mulberry by the two algorithms.

Discussion
qRT-PCR is the most commonly used method in gene expression studies because of its accuracy, sensitivity, and efficiency [19,20]. It is very important to select a suitable reference gene in qRT-PCR analysis to correct for any errors in RNA quantity and reverse transcription efficiency, and ultimately to determine the real expression of a target gene [21]. A suitable reference gene should maintain invariable expression across different experimental samples [22].
To our knowledge, this is the first published report to examine reference gene selection in mulberry, although some reference genes have recently been used in mulberry studies. ACTIN3 (NCBI number of HQ163775) and ACTIN4 (NCBI number of KT793030) has been used as reference genes in qRT-PCR analyses in mulberry [16,23,24]. EF-1α and UBI-as basic components of metabolic processes in organisms-have been used as control genes in expression analyses of the stress response in mulberry [25,26]. RPL15 (Ribosomal protein L15) was used as a reference gene in the qRT-PCR analysis of different tissues and stress treatments in mulberry [27,28]. In this study, we found that the same reference gene ACTIN3 (HQ163775) was one of the most stable reference genes among most subset samples, suggesting ACTIN3 is a good housekeeping gene in mulberry. However, the stability of ACTIN4 (KT793030) and UBI genes showed not very well (Table 2, Fig 3).
geNorm and NormFinder are two commonly used programs to evaluate and identify suitable reference genes [12,13]. In our study, the evaluation of the candidate reference genes from the two programs showed similar stability, although there were some differences between them (Table 2, Fig 3). The differences may be due to the distinct algorithms of geNorm and NormFinder [12,13]. The reference gene TUB2 had stable expression patterns in all four sample subsets in this study ( Table 2, Fig 3), and thus the reference gene can be considered good housekeeping genes in mulberry. In contrast, TUB3 and EF1α4 had varying levels of expression across the subsets ( Table 2, Fig 3) and should not be used as reference genes. The three tubulin gene family members (TUB2, TUB3, and TUB4) varied in their stability under the same experimental conditions. This phenomenon is consistent with findings from other plants. For example, GAPDHs are stably expressed in several plant species-such as Dioscorea opposita and Panax ginseng [29,30]-but are not stably expressed in watermelon (Citrullus lanatus) and celery (Apium graveolens) [9,15].
We found that TUB2, UBI4, ACTIN3 and RPL4 were ranked as the most stable reference genes in the samples subsets examined from the geNorm and NormFinder evaluation (Table 2, Fig 3). ACTIN and TUB, which encode the basic skeletal components of plant organelles, are the most commonly used reference genes and are usually stably expressed in tissues and organs in plants [31,32]. ACTIN7 is the most stably expressed gene in pear (Pyrus L.) and tung (Vernicia fordii) [33,34] but is not a suitable reference gene in Arabidopsis thaliana and Rhododendron [35,36]. TUBs are considered to be the best reference genes in water lily and soybean [37,38] but are not stably expressed in Rhododendron [35]. GAPDH, EF-1α, and UBI are involved in biochemical metabolic processes in organisms and are traditional reference genes in plants [31,32]. EF1α genes have stable expression patterns in cucumber, potato, and soybean plants [39][40][41]. However, EF1-α genes are not regarded as suitable reference genes in Arabidopsis and bamboo [11,42]. These results suggested that reference genes used in qRT-PCR should be validated in different species and under different conditions. TUB2, UBI4, ACTIN3 and RPL4 were ranked as the most stable reference genes in the samples subsets.
In conclusion, this is the first study to perform a systematic evaluation of mulberry to validate candidate reference genes for qRT-PCR normalization in different plant tissues and under different stress conditions. Twenty candidate reference genes were assessed. TUB2, UBI4, ACTIN3 and RPL4 were ranked as the most stable reference genes in different tissues and under different stress conditions of mulberry. Moreover, ACTIN3/TUB2, RPL4/ACT4, TUB2/GAPDH1 and UBI4/RPL4 were identified as the most stably expressed pairs in our study. These results provide useful guidelines for qRT-PCR data normalization of gene expression analysis in mulberry.

Plant material and treatments
Cuttings (about 10 cm in length and 1 cm in diameter) of mulberry (M. atropurpurea) cultivar "Tang 10" (bred by our unit) were rooted and planted in plastic pots (50% sand, 50% peat moss). One cutting with a sprout was planted in each pot. Plants were cultivated in a greenhouse (day/night temperature, 25/20˚C; relative humidity, 40-70%) at the Sericulture and Agri-Food Research Institute, Guangdong Academy of Agricultural Sciences, Guangzhou, China.
When they were~90 days old, the seedlings were used for stress treatments and tissue collection. Tissue samples of leaves, stems, and roots were collected from 3-month-old plants growing under well-watered conditions. Tissue samples of green fruit and red fruit were collected from the fruiting stage of adult mulberry cultivar "Tang 10". For drought treatment, plants were subjected to dry conditions with no irrigation for 10 days, at which point slight wilting occurred. Salt stress treatment consisted of watering with 200 mM NaCl for 5 days. For cold and heat shock treatments, plants were transferred to 4˚C and 42˚C, respectively, for 48 h. Wounding stress treatment consisted of mechanical wounding of a leaf three to five times with a sharp knife. The wounded area represented~10% of the leaf surface. Samples were collected from the last fully expanded leaf 24 h after the conclusion of each stress treatment. The samples from the five tissues and five abiotic stress treatments were collected from three replicate plants, giving a total of 30 samples comprised of 15 tissue-specific samples and 15 abiotic stress treatment samples. Samples were immediately frozen in liquid nitrogen and stored at -80˚C until RNA extraction.

RNA extraction and cDNA synthesis
Total RNA was isolated using the Mini-BEST Plant RNA Extraction kit (TaKaRa, Japan) with the addition of an on-column DNase I digestion according to the manufacturer's instructions. First-strand cDNA was synthesized from 1 μg of total RNA in a total volume of 20 mL per reaction using the PrimeScript RT reagent kit with gDNA Eraser (TaKaRa, Japan) following the manufacturer's protocol. The cDNA products were diluted 10-fold with nuclease-free water before being used in the qRT-PCR assays.

Primer design and qRT-PCR
Specific primers for qRT-PCR were designed using the Primer 5.0 software (PE Applied Biosystems, Foster City, USA) with primer lengths of 21-26 bp and amplicon lengths of 100-250 bp (Table 1 and S1 File). To determine the efficiency of each primer pair, a mixture of cDNA from the 30 samples was used to perform qRT-PCR reactions (see below). Five-point standard curves of a fivefold dilution series (1:1, 1:5, 1:25, 1:125, and 1:625) of the pooled cDNA were used. Agarose gel electrophoresis and polyacrylamide gelelectrophoresis (PAGE) electrophoresis of the amplification products of each candidate reference gene were analyzed.
qRT-PCR was carried out on a LightCycler480 System (Roche) using the SYBR Premix Ex Taq II kit (TaKaRa, Japan). Reactions were performed using a total volume of 20 μL, which contained 1 μL of cDNA template (corresponding to 5 ng of the starting amount of RNA), 0.2 mM each primer, and 10 μL 2× SYBR Premix Ex Taq II. The PCR cycling conditions were as follows: 94˚C for 30 s, followed by 40 cycles of 94˚C for 10 s, 55-62˚C for 10 s, and 72˚C for 10 s in a 96-well reaction plate, and the annealing temperature was based on the Tm value of primers. The melting curve was recorded after 40 cycles to verify primer specificity by heating from 65˚C to 95˚C. Each qRT-PCR reaction was performed in triplicate (technical replicates) on samples from three individual plants (biological replicates).

Statistical analysis
Two software programs, geNorm 3.5 and NormFinder 0.953, were used to assess the expression stability of each candidate reference gene according to their user manuals [12,13]. The Ct values were converted into relative expression values, which were calculated in Microsoft Excel 2007 using the highest expression value as the calibrator, and then they were imported into the geNorm and NormFinder software. All other statistical analyses were performed with Microsoft Excel 2007.

Normalization of the verified genes
Specific primers for the three MaCHS genes (MaCHS5, MaCHS6, and MaCHS7) and three plant hormone related genes (MaERF, MaDELLA, and MaJAZ) for qRT-PCR were designed using the Primer 5.0 software (S3 Table). The relative expression of the target gene was calculated using the 2 −ΔΔCt method [43] with normalization using the genes indicated. Measurement of the verified genes expression was made with three biological replicates with three technical repeats per sample.  (DOCX) S1 File. Primers location and the alignment of their gene sequences with homologous genes from other plants. (PPTX)