Evaluation of Reference Genes for Quantitative Real-Time PCR Analysis of the Gene Expression in Laticifers on the Basis of Latex Flow in Rubber Tree (Hevea brasiliensis Muell. Arg.)

Latex exploitation-caused latex flow is effective in enhancing latex regeneration in laticifer cells of rubber tree. It should be suitable for screening appropriate reference gene for analysis of the expression of latex regeneration-related genes by quantitative real-time PCR (qRT-PCR). In the present study, the expression stability of 23 candidate reference genes was evaluated on the basis of latex flow by using geNorm and NormFinder algorithms. Ubiquitin-protein ligase 2a (UBC2a) and ubiquitin-protein ligase 2b (UBC2b) were the two most stable genes among the selected candidate references in rubber tree clones with differential duration of latex flow. The two genes were also high-ranked in previous reference gene screening across different tissues and experimental conditions. By contrast, the transcripts of latex regeneration-related genes fluctuated significantly during latex flow. The results suggest that screening reference gene during latex flow should be an efficient and effective clue for selection of reference genes in qRT-PCR.


INTRODUCTION
Gene expression analysis provides the potential to explore biological processes. Northern blotting, transcriptome sequencing and qRT-PCR are universal methods that provide reliable and accurate gene expression analysis (Kubista et al., 2006;Josefsen and Nielsen, 2011;Lang et al., 2014;Le et al., 2014). Owing to its technical ease, cost-effectiveness, high specificity and accuracy, qRT-PCR has become the preferred method for detecting and measuring gene expression across different samples (Garson et al., 2005;Kubista et al., 2006). The accuracy of qRT-PCR is influenced by several factors such as the RNA extraction yield, reverse transcription and the efficiency of amplification (Huggett et al., 2005). It is necessary to normalize the expression data of the target genes with gene expression data of reference genes (Vandesompele et al., 2002;Huggett et al., 2005). A reference gene is defined as having a stable expression pattern across different tissues and experimental conditions. Generally, several housekeeping genes such as actin (ACT), eukaryotic elongation factor 1-alpha (eEF-1a), 18S rRNA, glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and polyubiquitin (UBQ), have been widely used as reference genes for qRT-PCR in the past decade (Bustin, 2000;Goidin et al., 2001;Kim et al., 2003;Jain et al., 2006). However, accumulating data have shown that the expression levels of these genes are variable across different developmental stages and experimental conditions (Garson et al., 2005). The use of non-validated reference genes could lead to inaccurate results and consequently, erroneous conclusions (Tricarico et al., 2002;Huggett et al., 2005;Remans et al., 2014). The 11 golden rules and Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines are recommended for qRT-PCR performance in the last decade (Udvardi et al., 2008;Bustin et al., 2009Bustin et al., , 2013. Several statistical algorithms, such as geNorm and NormFinder, have been developed for the selection of reference genes for qRT-PCR analysis (Vandesompele et al., 2002;Andersen et al., 2004). The geNorm is a Visual Basic application tool which calculates the expression stability value (M) for all reference genes. The gene with a lower M value suggests its expression is more stable. Additionally, the pairwise variation (Vn/Vn+1) of references can also be estimated by using geNorm. It is recommended that no additional genes are required for normalization when the value of Vn/Vn+1 is below 0.15 (Vandesompele et al., 2002). As another statistical algorithm, the NormFinder ranks the stability of candidate reference genes and calculates not only the candidate reference gene variation but also the variation between subgroup samples in a sample set (Andersen et al., 2004). The lowest variety value indicates the most stable. Application of both algorithms contributes to identify the best stable reference genes for different experimental samples.
The rubber tree (Hevea brasiliensis) is the main source of natural rubber due to the good yield and excellent physical properties of its rubber products (Eschbach and Lacrotte, 1989). Natural rubber biosynthesis takes place in latex, the cytoplasm of highly specialized laticifer cells (Chow et al., 2012;Wang et al., 2013Wang et al., , 2015. In natural rubber production, latex is exploited by successive tappings. Duration of latex flow after tapping is one of crucial factors that determines rubber yield of rubber tree, and easily influenced by multiple factors such as genotypes, age of trees, seasons, etc. (Sethuraj, 1992;Raj et al., 2005;Chao et al., 2015a). Tapping causes rubber tree to be wounded and loose latex (Tang et al., 2010). The loss of latex in turn enhances latex regeneration in laticifer cells (Jacob et al., 1989). Previous studies show that tapping activates the expression of mass of latex metabolism-related genes at transcriptional level, which mainly relay on the polymerase II transcriptional complex that are responsible for eukaryotic mRNA transcription (Hahn, 2004;Hager et al., 2009;Chao et al., 2015a). Such tapping effects should be mainly ascribed to latex flow which results in latex loss from laticifer cells. To screen reference genes on the basis of latex flow may be a more efficient clue than the previously reported methods on the basis of different treatments (Li et al., 2011;Pirrello et al., 2014;Long et al., 2015). In this study, 23 candidate reference genes were evaluated for their expression stability across 18 samples obtained from three stages of latex flow of two rubber tree genotypes with differential duration of latex flow by using geNorm and NormFinder algorithms.

Plant Materials
In this study, we used 11-year-old regularly tapped rubber tree clones CATAS7-33-97 and CATAS8-79 with the same circumference. The trees were grown at the Experimental Station of the Rubber Research Institute of the Chinese Academy of Tropical Agricultural Sciences in Danzhou city, Hainan province. After tapping, the latex samples were respectively collected at 1, 30, and 60 min for CATAS7-33-97, and at 1, 80, and 150 min for CATAS8-79, which represented the early, middle, and late stage of latex flow, respectively. Each batch of latex samples were individual collected from three trees for each clone, and stored at −80 • C until total RNA extraction.

RNA Isolation and cDNA Synthesis
Total RNA was isolated by using the RNAprep Pure Plant Kit (Tiangen, China) and genomic DNA was eliminated according to the manufacturer's instructions. The concentration and quality of RNA were examined using a NanoDrop 2000 (Thermo Scientific Inc., USA). The RNA integrity of the samples was checked using 1.5% agarose gel electrophoresis. Synthesis of cDNA was performed using a RevertAid TMFirst Strand cDNA Synthesis Kit (Fermentas, Canada) following the manufacturer's protocol.

qRT-PCR Analysis
A total of 23 housekeeping genes (18S rRNA, ACT, ACT7a, ACT7b, ADF, ADF4, CYP2, eIF2, eIF3, eIF1Aa, eIF1Ab, FP, PTP, RH2a, RH2b, RH8, ROC3, UBC1, UBC2a, UBC2b, UBC3, UBC4, and YLS8), which were evaluated in previous studies (Li et al., 2011;Chao et al., 2015b;Long et al., 2015), were chosen as candidate reference genes in the present study (Supplementary  Table S1). Reactions were performed using a SYBR PrimeScript RT-PCR Kit (TaKaRa Biotechnology, Japan) on the CFX96 System (Bio-Rad Laboratories Inc., USA) as follow: 95 • C for 3 min followed by 45 cycles at 95 • C for 15 s, 56 • C for 60 s and 72 • C for 30 s, and a melt curve from 55 to 95 • C increasing by 0.5 • C every 30 s. Each real-time PCR reaction was performed in triplicate. The Bio-Rad CFX96 Manager 3.0 software was used for visualizing and analyzing the data, including the cycle threshold (Ct) values, amplification efficiencies (E) and melting curve of reference genes (Supplementary Table S1; Supplementary Figure S1). The primer efficiency of each gene was evaluated using a standard curve generated by qRT-PCR using a fivefold dilution series of cDNAs (Supplementary Table S1).

Data Analysis
Based on the Ct values for the 23 candidate reference genes, a boxplot was drawn using GraphPad Prism 5 software 1 to show the expression variation of each gene. The expression stability of the reference genes was analyzed using the geNorm (version 3.5 2 ) and NormFinder (version 0.953 3 ) software packages. To evaluate the gene expression stability, the Ct values of the candidate reference genes were converted into relative quantities using the formula E − Ct ( Ct = Ct value of each sample-the minimum Ct value) (Li et al., 2016), and imported into the software according to the corresponding manuals of the Norm-Finder and geNorm algorithms (Vandesompele et al., 2002;Andersen et al., 2004).

The Expression of Latex Regeneration-Related Genes during Latex Flow
The HbHMGR1 (Genebank accession number: AB294692) and small rubber particle protein (HbSRPP, Genebank accession number: AEH05971) were rubber biosynthesis-related genes (Sando et al., 2008;Guo et al., 2014). HbRpb11 (Genebank accession number: KU301750) and HbTFIIB (Genebank accession number: KU301749) were the genes homologous to RNA polymerase II transcriptional complex members. They were selected to represent latex regeneration-related genes. Specific primers were designed using Primer Premier 5 software 4 (Supplementary Table S1) and the amplified efficiency was checked (Supplementary Figure S1). The UBC2b, 18S rRNA and FP were respectively used as reference genes to normalize the transcripts of the target genes. The relative expression level of all of the target genes was evaluated using the 2 − Ct method Hu et al., 2016).

Statistical Analysis
Statistical analysis was performed with SPSS Statistics 17.0 5 by analysis of variance (ANOVA) based on Duncan's test (multiple group comparisons) or a t-test (two group comparisons). For Duncan's test, the capital letter or lowercase represented p < 0.01 or p < 0.05, and the same letter indicated no significant difference among groups. For t-test, * * * represented p < 0.001.

Determining the Time Intervals for Latex Collection Based on the Different Durations of Latex Flow in the CATAS8-79 and CATAS7-33-97
Latex is usually obtained by tapping the trunk bark of rubber tree every 2-day interval ( Figure 1A). The rubber tree clones CATAS8-79 and CATAS7-33-97 exhibited a huge difference in the duration of latex flow and latex production (Figures 1B,C). For CATAS8-79, the latex flow lasted for more than 160 min after tapping, which was significantly longer than that of CATAS7-33-97 (∼70 min) ( Figure 1B). As a result, the latex production of CATAS8-79 was notably higher than that of CATAS7-33-97 ( Figure 1C). Based on the differential duration of latex flow between the two clones, latex samples were collected at 1, 30, and 60 min after tapping for CATAS7-33-97 and at 1, 80, and 150 min after tapping for CATAS8-79, which represented the early, middle and late stages of latex flow, respectively.

Experiment Set Classification and Transcript Abundance Evaluation
The qRT-PCR data from 18 samples were divided into three experiment sets to assess the influence of latex flow on the expression of the candidate reference genes. The first set consisted of nine samples from CATAS7-33-97. The second set was comprised of nine samples from CATAS8-79. The third set included both the first and second sample sets.
The Ct values for 23 candidate reference genes were determined across all of the experimental samples using qRT-PCR (Supplementary Table S1; Supplementary Figure S1). The boxplot analysis was performed using GraphPad Prism 5 software (Figure 2). The results indicated that the selected candidate reference genes displayed a wide range of Ct values, ranging from 9 to 26. The Ct values of the majority of the candidate reference genes ranged from 19 to 22. The 18S rRNA gene was the most abundant reference gene with the lowest mean Ct value (9.85), whereas FP was the least abundant reference gene with the highest mean Ct value (25.94) (Figure 2).

Evaluation of Expression Stability Using geNorm
geNorm was used to assess the expression stability (M) of each reference gene based on the average pair-wise variation among all of the tested genes (Vandesompele et al., 2002). The lowest M value indicates the most stable expression, whereas a higher M value is indicative of less stable expression. All of the 23 reference genes had an M value below the geNorm threshold of 1.5 and were ranked according to their M values in three sample sets (Figures 3A-C). For CATAS7-33-97 sample set, UBC2a and eIF3 were the most stable genes with the M values of 0.13 ( Figure 3A). The seven top-ranked genes were UBC2a/eIF3, UBC2b, RH8, UBC3, eIF1Ab and ROC3. For CATAS8-79 sample set, UBC2a and UBC4 were the most stable genes with the M values of 0.17 ( Figure 3B). The seven top-ranked genes were UBC2a/UBC4, YLS8, ADF, ADF4, UBC1 and UBC2b. For total experiment set, RH8 and UBC2b were the most stable genes with the M values of 0.16 ( Figure 3C). The seven top-ranked genes were RH8/UBC2b, UBC2a, UBC4, eIF3, UBC1 and ADF4.
To determine the optimal number of reference genes required for normalization, pairwise variation analysis (Vn/Vn + 1) was performed using the geNorm software ( Figure 3D). The cut-off threshold value for Vn/n + 1 is 0.15 and no additional reference gene is required for normalization when the value below this value (Vandesompele et al., 2002). In the present study, the value The difference in latex production between CATAS7-33-97 and CATAS8-79. Significant difference was indicated by the asterisks above the bars ( * * * p < 0.001).  (Figure 3D), which suggests that two reference genes would be sufficient for gene normalization in each experimental set.

Evaluation of Expression Stability Using NormFinder
NormFinder is an algorithm tool to identify optimally stable reference genes by determining the expression stability (Andersen et al., 2004). Here, we use NormFinder software to further confirm the results obtained using the geNorm program. The data are listed in Table 1. The seven top-ranked genes were UBC2b, RH8, eIF3, UBC3, UBC2a, eIFAa, and ACT7a in the CATAS7-33-97 experiment set ( Table 1), six of which were also present in the top seven genes in the geNorm evaluation ( Figure 3A). The seven top-ranked genes were UBC2a, UBC2b, UBC4, ADF4, YLS8, ADF, and RH8 in the CATAS8-79 experiment set ( Table 1), six of which were also present in the top seven genes in the geNorm evaluation ( Figure 3B). The UBC2a, UBC2b, and RH8 were common between the top seven genes in both experiment sets using both NormFinder and geNorm. Moreover, they were the top three genes in the combined experiment set using both NormFinder ( Table 1) and geNorm ( Figure 3C).

Influence of Latex Flow on the Expression of Latex Regeneration-Related Genes
The HbHMGR1, HbSRPP1, HbTFIIB and HbRpb11 were selected to analyze the expression of latex regeneration-related gene at different stages of latex flow. The relative expression was respectively normalized against reference genes UBC2b, 18S rRNA and FP. The UBC2b was the most stable reference gene while the traditionally used 18S rRNA and the more stable FP reported by Li et al. (2011) were less stable in the present study. The transcript level of all the four genes were significantly fluctuated during latex flow when normalized by any of the three reference genes (Figure 4). The fluctuating pattern of each target gene was, however, different to some extent among the three normalizations (Figure 4). There was no difference in the transcript level of HbHMGR1 between "middle" and "late" stages of latex flow by using UBC2b normalization in either CATAS7-33-97 or CATAS8-79 ( Figure 4A). But it showed significant difference (p < 0.01) by using 18S rRNA normalization in CATAS7-33-97 and difference (p < 0.05) by using FP normalization in CATAS8-79 ( Figure 4A). The transcript level of HbSRPP between "middle" and "late" stages showed significant difference (p < 0.01) in CATAS7-33-97 while had no difference in CATAS8-79 by using UBC2b normalization ( Figure 4B). It had no difference by using FP normalization while showed significant difference (p < 0.01) by using 18S rRNA in either CATAS7-33-97 or CATAS8-79 ( Figure 4B). The relative transcript level of HbRpb11 had no difference between "middle" and "late" stages of latex flow in CATAS7-33-97 by using either UBC2b or FP normalization while significant difference (p < 0.01) by using 18S rRNA normalization ( Figure 4C). It showed difference in CATAS8-79 for 18S rRNA normalization and significant difference for both UBC2b and FP normalization ( Figure 4C) The UBC2b normalized expression pattern of HbTFIIB in either CATAS7-33-97 or CATAS8-79 was similar with the FP normalized but different from the 18S rRNA normalized expression pattern ( Figure 4D).

DISCUSSION
Accurate interpretation of the qRT-PCR results mainly depends on the use of stable reference genes for data normalization to minimize non-biological variations among samples (Huggett et al., 2005). However, increasing evidence shows that the best reference gene significantly depends on the experimental conditions (Garson et al., 2005). Therefore, validating suitable reference genes is necessary to obtain accurate gene expression patterns under specific conditions using qRT-PCR (Ji et al., 2014;Wang et al., 2014;Delporte et al., 2015;Kong et al., 2015;Niu et al., 2015). The rubber tree is one of the most economically important members of the genus Hevea because the milky latex exploited from the tree is the primary source of natural rubber. The latex exploitation is an effective factor that enhances latex regeneration (Jacob et al., 1989). After tapping, the loss of latex from latex flow results in significant changes in turgor pressure of laticifer cells, which may significantly reprogram the transcriptional regulation of latex regeneration-related genes in laticifer cells at the molecular level (Chao et al., 2015a). In the present study, we showed that not only the rubber biosynthesisrelated genes such as HbHMGR1 and HbSRPP significantly fluctuated, but also the primary metabolism-related genes such as HbTFIIB and HbRpb11, the homologs of transcriptional complex members, were notably reprogrammed during latex flow (Figure 4). The results indicate that latex regenerationrelated genes are strongly influenced in the process of latex flow. This factor is largely ignored when other factors such as tapping and hormone applications are considered in screening reference genes (Li et al., 2011;Pirrello et al., 2014;Long et al., 2015).
The expression stability of 22 candidate reference genes was evaluated based on tapping, genotypes, tissues, and plant hormones (Li et al., 2011). The authors suggested that the top five reference genes were UBC2a, UBC2b, UBC4, eIF3, and FP. Thereafter, 20 of the 22 candidate reference genes were evaluated using successive tapping and the top five reference genes were UBC2a, UBC2b, UBC4, UBC3, and ADF (Long et al., 2015). Three genes were identical between the top five genes. The RH2b, UBC2a, RH8 are the top three best references for abiotic and harvesting stress of rubber tree by screening 11 of the 22 candidate genes and simply evaluating based on the standard deviation and coefficient of variance of Ct value that generated by qRT-PCR (Pirrello et al., 2014). In the present study, the selected 23 candidate reference genes included the 22 candidate reference genes (Li et al., 2011) and an ACT gene (Chao et al., 2015b). The expression stability of all the 23 candidate reference genes was evaluated on the basis of latex flow using geNorm and NormFinder algorithms. The top five genes were UBC2a, UBC2b, UBC4, eIF3, and RH8. Four of the five genes are identical to the genes reported by Li et al. (2011) and three of the five genes are identical to those reported by Long et al. (2015). It is interesting that three UBC family members (UBC2a, UBC2b, and UBC4), are common among these experimental groups of rubber trees (Figure 5).
UBC is a large gene family comprising dozens of members, which widely expands in eukaryotes (Sheng et al., 2012). In the present study, two UBCs, the UBC2a and UBC2b, had the best stability based on latex flow in rubber trees. It is interesting that UBC2a and UBC2b are also stable in laticifer cells and bark tissue of rubber trees. For example, UBC2b was also found to be highly stable across multiple experimental conditions, such as hormone treatments, genotypes and latex regeneration between two tapping intervals (Li et al., 2011;Long et al., 2015). UBC2a was the most stable gene in the bark tissues of rubber trees treated with coronatine (data unpublished). The homolog of UBC2a and UBC2b is also stable in many species such as Brachypodium distachyton (Chambers et al., 2012), Trifolium pretense (Mehdi Khanlou and Van Bockstaele, 2012), Lolium perenne FIGURE 5 | Comparing suitable reference genes screened here with previous studies on rubber tree. The top five best references of each experimental group were analyzed here. (Huang et al., 2014), and Withania somnifera (Singh et al., 2015) under different experimental conditions. In plants, UBC2a and UBC2b encode ubiquitin-conjugating enzyme subunits that are necessary for the degradation of ubiquitinated proteins by the eukaryotic conserved 26S protein degradation pathway (Stone, 2014). Protein ubiquitination is mediated by the sequential action of E1 (ubiquitin-activating enzyme, UBA), E2 (ubiquitin-conjugating enzyme, UBC), and E3 (ubiquitinligating enzyme) (E et al., 2015). E2 UBC enzymes accept activated ubiquitin from E1 through a transthiolation reaction, and then transfer ubiquitin to either a substrate directly aided by E3 or a cysteine of ubiquitin protein ligase through a second transthiolation reaction that then transfers ubiquitin to the substrate (Jue et al., 2015). The housekeeping of UBCs may attribute to the fact that UBCs function as an intermediate step of the ubiquitin-proteasome degradation system to control the abundance of cellular proteins required for plant growth and development (Stone and Callis, 2007;Vierstra, 2009). Taken together, screening reference genes on the basis of latex flow should be an efficient and effective clue for selection of reference genes in qRT-PCR. The UBC2a and UBC2b could be used as suitable reference genes for qRT-PCR analysis of latex regeneration-related genes following different treatments.

AUTHOR CONTRIBUTIONS
JC designed and carried out the experiment of this study, and wrote the manuscript. SY, YC participated and analyzed data in the experiment. WT planned the study and participated in the design of the experiment. All authors have read and approved the manuscript in its final form.