ORIGINAL RESEARCH article
Validation of Suitable Reference Genes for Assessing Gene Expression of MicroRNAs in Lonicera japonica
- 1State Key Laboratory Breeding Base of Dao-di Herbs, National Resource Center for Chinese Materia Medica, China Academy of Chinese Medical Sciences, Beijing, China
- 2Beijing Institute of Genomics, Chinese Academy of Sciences, Beijing, China
MicroRNAs (miRNAs), which play crucial regulatory roles in plant secondary metabolism and responses to the environment, could be developed as promising biomarkers for different varieties and production areas of herbal medicines. However, limited information is available for miRNAs from Lonicera japonica, which is widely used in East Asian countries owing to various pharmaceutically active secondary metabolites. Selection of suitable reference genes for quantification of target miRNA expression through quantitative real-time (qRT)-PCR is important for elucidating the molecular mechanisms of secondary metabolic regulation in different tissues and varieties of L. japonica. For precise normalization of gene expression data in L. japonica, 16 candidate miRNAs were examined in three tissues, as well as 21 cultivated varieties collected from 16 production areas, using GeNorm, NormFinder, and RefFinder algorithms. Our results revealed combination of u534122 and u3868172 as the best reference genes across all samples. Their specificity was confirmed by detecting the cycling threshold (Ct) value ranges in different varieties of L. japonica collected from diverse production areas, suggesting the use of these two reference miRNAs is sufficient for accurate transcript normalization with different tissues, varieties, and production areas. To our knowledge, this is the first report on validation of reference miRNAs in honeysuckle (Lonicera spp.). Restuls from this study can further facilitate discovery of functional regulatory miRNAs in different varieties of L. japonica.
Lonicera japonica Thunb. (FLJ), commonly known as Japanese Honeysuckle or “Ren Dong,” has been widely used in traditional medicine in East Asian countries for thousands of years (Jiang et al., 2009). Its dry flower buds, leaves, and stems may be used to prevent and to treat a variety of illnesses, such as hand-foot-and-mouth disease, H1N1 influenza, pancreatic cancer, and some severe acute respiratory syndromes (Kwak et al., 2005; Kang et al., 2007; Lin et al., 2016). However, different parts of this plant have different medicinal properties: the flower buds have anti-inflammatory and anticancer properties, the leaves have antioxidant and tyrosinase-inhibitory activities, and the stems have xanthine oxidase-inhibitory and nitrite-scavenging activities (Vukovic et al., 2012). Several secondary metabolites have been identified and isolated from this plant, including biflavonoids, quercetin, phenolic acids, and dicaffeoylquinic acid (Fu et al., 2013), which are thought to be responsible for these various biological and pharmacological activities (Liu et al., 2016; Zhang B. et al., 2016). Quantitative differences, comparing the chemical composition of different parts of L. japonica, have been investigated: (1) the percentage of oxygenated monoterpenes and chlorogenic acid decreased from flower > leaf > stem; (2) the level of hexadecanoic acid increased from flower < leaf < stem; and (3) the content of luteoloside decreased from leaf > flower > stem (Vukovic et al., 2012; He, 2013; Jiang et al., 2013). In contrast, the content of active compounds also differs significantly among the varieties and production areas of L. japonica, leading to the different pharmacological activities and medicinal qualities (Yuan et al., 2012; Shi et al., 2016; Wu et al., 2016). L. japonica Thunb. var. chinensis (Watts; rFLJ), a Chinese endemic variety, has different active compound contents from those of FLJ (Yuan et al., 2012; Supplementary Figure S1). Both varieties of L. japonica are widely planted in different production areas of China, which can also obviously affect the content of the bioactive components (Yuan et al., 2012; Li et al., 2013; Tan et al., 2016). Previous investigations revealed that 87 EST-SSRs were significantly different between L. japonica and L. japonica var. chinensis in diverse production areas (Jiang et al., 2012). Moreover, random amplified polymorphic DNA (RAPD) markers could identify five varieties of L. japonica collected from different geographic locations of China (Fu et al., 2013). However, only limited information is available on the biosynthetic enzymes involved in the synthesis of active compounds and their regulation in different organs and varieties of L. japonica collected from diverse production areas. Therefore, it is necessary to focus on elucidating the molecular regulatory mechanisms of secondary metabolism in flower buds, stems, and leaves from both varieties of L. japonica in various production areas.
MicroRNAs (miRNAs), which are approximately 19 to 24 nucleotides (nt) in length, are a class of small endogenous non-coding RNAs (sRNA) that originate from their precursors (short stem-loop structures) within primary transcripts (pri-miRNA; Bartel, 2004; Barakat et al., 2007). The mature miRNAs play a key role in negative control of gene expression mainly by inhibiting the translation of target genes or by cleaving their target mRNAs, and thereby exert control over a series of key developmental and stress resistance processes (Lu et al., 2008; Bartel, 2009). To date, 8496 mature miRNAs and 6992 precursor miRNAs (pre-miRNAs) have been identified in Viridiplantae (miRBase1, April 5, 2016), among which six frequently used medicinal herbs have been investigated, including Panax ginseng, Digitalis purpurea, Rehmannia glutinosa, Salvia sclarea, Citrus reticulate, and C. sinensis (Supplementary Table S1). Increasing evidence indicates that miRNAs play important regulatory roles in diverse biological and metabolic processes such as plant secondary metabolism and environmental stress resistance (Shi et al., 2010; Wu et al., 2012a ; Carlile and Werner, 2014). Thus, miRNAs have the potential to serve as biomarkers for differences in medicinal plant varieties and cultivation environmental conditions. Moreover, this class of critical regulators could be a useful tool to determine the regulation mechanisms of secondary metabolism in both varieties of L. japonica in response to multiple environmental factors. However, only limited honeysuckle miRNAs have been reported to date (Zhou et al., 2014; Hirschi et al., 2015; Yang et al., 2015).
In miRNA expression studies, the most common choices of reference genes are ribosomal RNAs (e.g., 5S RNA) and small nuclear RNAs (e.g., U6; Reid et al., 2006; Lagesen et al., 2007). However, the expression levels of these genes may vary across different samples or because of external factors (Kou et al., 2012; Wang et al., 2014; Wu et al., 2014). Moreover, a commonly used reference gene in an organism may not be appropriate in other organisms or under different conditions. Thus, it will be necessary to evaluate the possibility that endogenous miRNAs in plant tissues may serve as potential internal reference genes. In recent years, there have been extensive reports on the expression of miRNAs. For example, mi167 and mi159 are the two most suitable reference genes in wheat (Hao et al., 2012), and miR169, miR171/170, and miR172 are stably expressed in lettuce under abiotic stresses (Borowski et al., 2014). In addition, the combination of miR156b and miR1520d was reported as a normalizer of soybean miRNA rather than the more commonly used protein-coding genes (Kulcheski et al., 2010). However, to our knowledge, no systematic study has been carried out in different organs and varieties of L. japonica for validation of reference miRNAs.
Here, we evaluated the stability of a set of 16 candidate reference genes, including two general reference genes (u437272 from the 5S ribosomal RNA gene and u1325500 from the U6 small nuclear RNA gene), to identify the most suitable reference miRNA for qRT-PCR data normalization in different tissues and varieties of L. japonica obtained from different production areas. The optimum pairs of internal control genes were identified using GeNorm, NormFinder, and BestKeeper statistical algorithms. Our results provide an important basis for further studies on molecular regulatory mechanisms of secondary metabolism in flower buds, leaves, and stems and investigation of the potential biomarkers for differences in varieties of L. japonica.
Materials and Methods
Honeysuckle Sample Preparation
Honeysuckle samples from different geographic locations were collected from April to September 2013 (Table 1). Flower buds, young leaves, and stems were sampled before flowering for gene expression analysis. All samples were immediately frozen in liquid nitrogen and stored at -80°C. Three cultivars of honeysuckle (BJ-YTH, SD-YTH, and SD-DMH) were used for sequencing of small RNAs and characterization by RNA-seq. Each sample tissue of all honeysuckle cultivars consisted of four plants for the expression stability analysis and ranking of candidate reference genes. In addition, eight samples of each honeysuckle cultivar were used to validate the reference miRNAs. All the experiments were repeated in triplicate.
Extraction of miRNA, cDNA Preparation, and qRT-PCR
Total RNA was isolated from different tissues of L. japonica using a TRIzol kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. Isolated RNA integrity was tested by agarose gel electrophoresis. The quality and concentration of the extracted RNA were evaluated using a NanoDrop 2000 spectrophotometer (Thermo Scientific, San Jose, CA, USA).
Total RNA was thawed on ice, and miRNA was extracted using a miRcute miRNA Isolation Kit (TIANGEN, Beijing, China) following the manufacturer’s instructions. The quantity and purity of extracted miRNA were estimated by examining both the absorbance at 260 nm and the 260/280 nm ratio. Subsequently, 0.5 μg of RNA per sample were reverse-transcribed to cDNA with a miRcute miRNA First-Strand cDNA Synthesis Kit (TIANGEN, Beijing, China). The above cDNA was then diluted 1: 10 in water, after which 2 μl was used for quantitative PCR, using the miRcute miRNA qPCR Detection kit (TIANGEN, Beijing, China), performed on an ABI 7500 real-time quantitative PCR instrument (Applied Biosystems, Carlsbad, CA, USA). Primers for mature reference miRNAs were obtained from Sangon (Sangon Biotech, Shanghai, China). The reaction was performed at 94°C for 2 min, followed by 40 cycles of 94°C for 20 s, and 60°C for 34 s. Melting curves were performed immediately after the completion of the qRT-PCR and fluorescence was measured from 55 to 95°C. PCR efficiency (E; Tichopad et al., 2003) for each miRNA was determined (Table 2).
Selection of Candidate Reference Genes and Primer Design
Using high-throughput sequencing of the miRNA transcriptome of L. japonica, we identified a set of fourteen miRNAs that were equally abundant in three small RNA-seq libraries (Accession number SRR3567675; Supplementary Table S2). The Arabidopsis thaliana 5S rRNA gene (AJ307353.2) and Solanum tuberosum U6 small nuclear RNA gene (Z17301.1) were used as queries in Blast searches of the NCBI nucleic acid database to search for homologous sequences in the L. japonica miRNAs database, using identification standards of e-value ≤1e-15, score ≥100. u437272 from the 5S ribosomal RNA gene and u1325500 from the U6 small nuclear RNA gene were selected to serve as general reference genes. Specific poly(A) primers were designed for miRNA cDNA synthesis as previously described (Rui and Chiang, 2005). miRNAs were polyadenylated and reverse-transcribed with a poly(A) adapter into cDNAs for real-time PCR using the miRNA-specific forward primer and the sequence complementary to the poly(A) adapter as the reverse primer. Specific forward primers for each full miRNA sequence were designed (Table 2) and the reverse primer sequence was universal (TIANGEN, Beijing, China) for qRT-PCR amplifications.
Analysis of Stability of Reference miRNAs
To assess the stability of candidate reference miRNAs, statistical algorithms including GeNorm, NormFinder, and BestKeeper were utilized. The statistical software package GeNorm is based on pair-wise comparisons and calculates the M-value of all candidate genes (Vandesompele et al., 2002). A valuable feature provided by GeNorm output is the V-value that reflects the minimal number of genes required for a robust normalization factor (NF; Teste et al., 2009). NormFinder is an ANOVA-based model that provides direct measures of intra- and inter-group variation in the expression of each gene (Tong et al., 2009; Liu et al., 2012). BestKeeper determines the “optimal” housekeeping genes using pair-wise correlation analysis of all pairs of candidate genes, as well as the target gene (Pfaﬄ et al., 2004; Takle et al., 2007). Data analyses were performed with SPSS 19.0 statistical software (SPSS Inc., Chicago, IL, USA). Significance in relative expression levels of miRNAs were determined by one-way ANOVA. P-values below 0.05 were considered statistically significant for all tests.
Selection of Candidate Reference Genes, Primer Specificity, and Amplification Efficiency in L. japonica
For the evaluation of potential endogenous normalizers in L. japonica, 16 candidate reference miRNAs were selected from small RNA libraries, which were chosen in preference according to the following criteria: the mean variation (RPKM value) in miRNA expression levels was <1.1-fold in three miRNA sequencing libraries of different varieties of L. japonica and no significant differences in miRNA expression levels (P > 0.05). Those precursors of selected miRNAs were compared with miRBase miRNA sequences2 to search for homologs of microRNA sequences (Table 2). Additionally, a plant microRNA database (PMRD) was used to predict the secondary structures and examine sequence similarity of the selected reference miRNAs in this study (Figure 1). As expected, lj-miR167a and lj-miR171b were the orthologous sequences of miR167a and miR171b, which have been described in the literature as reference miRNAs in Lotus japonicas and maize (Zhang et al., 2009; De Luis et al., 2012). Interestingly, except for lj-miR167a and lj-miR171b, orthologous sequences of selected miRNAs were quite different but had similar secondary structure sequences, which suggested that the selected reference genes in this study might have unknown functions deserving further investigation.
FIGURE 1. Secondary structure prediction of 16 candidate reference miRNAs and their similar sequences using the plant microRNA database (PMRD). (A) Secondary structures of lj-miR167a and its similar sequence ptc-MIR167g. (B) Secondary structures of lj-miR171b and its similar sequence vvi-MIR171d. (C) Secondary structures of u30297 and its similar sequence mtr-MIR2111q. (D) Secondary structures of u1846379 and its similar sequence mtr-MIR169g. (E) Secondary structures of u1760353 and its similar sequence ghr-MIR399c. (F) Secondary structures of u3464767 and its similar sequence ppt-MIR1031b. (G) Secondary structures of u312335 and its similar sequence tae-MIR1127. (H) Secondary structures of u4339213 and its similar sequence ptc-MIRfl2395. (I) Secondary structures of u3817076 and its similar sequence sit-MIR142-npr. (J) Secondary structures of u2100564 and its similar sequence osa-MIR816. (K) Secondary structures of u821189 and its similar sequence tae-MIR1127. (L) Secondary structures of u534122 and its similar sequence gma-MIR5371. (M) Secondary structures of u3868172 and its similar sequence osa-MIRf11793-akr. (N) Secondary structures of u4631289 and its similar sequence csi-MIR827. (O) Secondary structures of u1325500 and its similar sequence vvi-MIR166c. (P) Secondary structures of u437272 and its similar sequence osa-MIRf11621-akr.
The performance of each primer pair was tested by qRT-PCR. The amplification efficiencies for the 16 primers ranged from 90.1 to 106.9%, and R2 ranged from 0.986 to 0.999, according to the slopes of the standard curves (Table 2; Supplementary Figure S2). Melting curve analysis confirmed the presence of a single PCR product with a single peak from all samples, indicating that all the primers used are highly specific and efficient for qRT-PCR amplification (Supplementary Figure S3).
Expression Profiling of Candidate Reference miRNAs
Expression levels of all 16 candidate reference genes were determined by qRT-PCR. The Ct values showed differential transcript levels in three tissue samples examined, with lower Ct values suggesting higher transcript abundance and vice versa. The range of Ct values indicated the stability of reference gene expression; the wider the range of Ct values for a gene, the more unstable the gene’s expression (Hong et al., 2008; Gao et al., 2012). As shown in Figure 2, the mean Ct value of all genes ranged from 21 to 31. Among all the genes, u1846379 and u4631289 were the two most abundant with mean Ct values <24, and u1325500 always had the lowest level of expression with mean Ct value >30, irrespective of the samples analyzed. The Ct value could change in different tissues. In general, u2100564 accumulated in leaves far more than in other tissues (P < 0.05); however, u3464767 had a lower level of expression in leaves compared with stem and flower buds (P < 0.05). In addition, the range of Ct values for u30297, u3464767, and u4339213 was above six cycles in flower buds, a range wider than that of the other genes tested, which indicated these miRNAs were unstable and inappropriate for reference genes. Notably, u3817076, u821189, u534122, and u3868172 showed a comparatively stable and narrower range of expression variation in all sample sets, indicating that they could possibly serve as reliable reference genes. In contrast, Ct values of u30297, u4339213, and u437272 were highly variable in all sample sets. This variation may be mainly the result of the different expression levels of these miRNAs in different tissues or the diversity in RNA yield and quality of different samples.
FIGURE 2. Transcript abundance of 16 candidate genes analyzed by GeNorm in three Lonicera japonica datasets. (A) Flower buds; (B) leaves; (C) stem. The solid line represents the median value and the boxes are 25th and 75th percentiles. The average is indicated by the point in the box. Whisker caps represent the maximum and minimum Ct values.
Expression Stability Analysis and Ranking of Candidate Reference miRNAs
Expression stability of the candidate reference genes was determined by GeNorm, NormFinder, and RefFinder which evaluated the stability ranking of each reference gene using the Ct values across all the experimental sets and tissue samples (Kozera and Rapacz, 2013; Luo et al., 2014) (Figure 3; Supplementary Tables S3 and S4). The stability ranking analyses done by each of the three software packages are detailed below.
FIGURE 3. Ranking of candidate reference genes using GeNorm analysis of qRT-PCR of candidate reference genes. The x-axis from left to right indicates the ranking of the genes according to their expression stability: (A) flower buds; (B) leaves; (C) stems; (D) all samples.
The GeNorm software, which has been used for reference gene selection in many plants, was developed for ranking the candidate genes and determination of the number of reference gene that should be applied. This method is based on the principle that the expression ratio of two ideal internal reference genes is identical in all samples regardless of the cell types or treatments. The pairwise variation of a particular gene with all the other reference genes is calculated and defined as the M-value, where a lower M-value represents a slight variation. The gene expression stability of all 16 candidate reference genes was evaluated using the GeNorm statistical algorithm. According to GeNorm analysis, the cut-off range of the stability value (M) is <1.5, so the gene with the lowest M-value is considered to be the most stable reference gene in terms of gene expression and vice versa. All tested datasets reached a high expression stability with relatively low M-values far below the default limit of M < 1.5. Gene expression stability of the 16 candidate genes was analyzed on four different datasets: (1) the M-values of lj-mir167a (M = 0.21) and lj-mir171b (M = 0.25) were the least followed by u1760353 and u821189 in flower buds (Figure 3A); (2) u1760353 (M = 0.31) and lj-mir171b (M = 0.39) were the most stable miRNAs in leaves followed by u3868172 and u534122 (Figure 3B); (3) u534122 (M = 0.32) and u3868172 (M = 0.41) were the most stable miRNAs in stem samples followed by u4339213 and u3817076 (Figure 3C); (4) the gene pair u534122 and u3868172 was the most stable combination of genes in all tissue datasets (Figure 3D). Although the M-values of lj-mir167a and lj-mir171b were the least in flower buds, they were not as stable in stems and leaves as u534122 and u3868172. Taking all of the above datasets together, u534122 and u3868172 ranked as the highest stable reference genes, indicating that they have the potential to be widely used for multiple purposes.
The GeNorm algorithm was also used to determine the optimal number of reference genes, as inclusion of two or more internal control genes for transcript normalization in qRT-PCR experiments is considered better for obtaining precise and consistent results. GeNorm determines the pairwise variation (Vn/Vn+1) between sequential NF (NFn and NFn+1) in each sample set. A threshold of 0.15 was proposed to indicate that inclusion of an additional reference gene is unnecessary. Thus, according to the stability value plot (Figure 4), the top “n” reference genes are sufficient as internal controls. Taking all tissue samples together, the V2/3 values were 0.1134, 0.1012, 0.1362, and 0.1322 using four different statistical algorithms, respectively (Figure 4), indicating that the top two reference miRNAs would be adequate for normalization, with a third miRNA being unnecessary. Thus, u534122 and u3868172 should be included as reference genes for accurate gene expression normalization in different tissues of L. japonica.
FIGURE 4. Pairwise variation (V) analysis to determination the optimal number of reference genes for normalization using GeNorm analysis. GeNorm calculates pairwise variation (Vn/Vn+1) for the normalization factors NFn and NFn+1 to determine (V < 0.15) the optimal number of reference genes.
NormFinder is a model-based algorithm to estimate expression variation of candidate reference genes. This strategy enables estimation of not only the overall variation of the genes but also the variation among subgroups of the sample set. This program ranked the candidate miRNAs based on intra- and inter-group variation, with a lower stability value indicating more stable expression. As shown in Supplementary Table S3, NormFinder ranked the 16 candidate miRNAs in all samples from lowest to highest stability values as follows: u3868172, u534122, u312335, lj-mir171b, lj-mir167a, u3817076, u821189, u1760353, u4631289, and u30297. Similar to the results from GeNorm, NormFinder also identified u534122 and u3868172 as the most stably expressed miRNAs in all samples, with u3464767 being the least stably expressed.
BestKeeper analyzes gene expression variation for candidate genes by calculating standard deviation (SD) and a pairwise correlation, and the lowest SD value indicates the most stable reference miRNA expression. As shown in Supplementary Table S4, u534122 and u3868172 were the most stably expressed genes with the lowest SDs in all samples (0.77 and 0.84, respectively). Thus, the results obtained through the BestKeeper program were consistent with those obtained through the GeNorm and NormFinder analyses.
Three ordered lists were generated according to the corresponding statistical values: M-value by GeNorm, stability value by NormFinder, and the correlation coefficient value of BestKeeper. According to the above results, u534122 and u3868172 represent the most adequate normalization candidate miRNAs tested. Only small differences between GeNorm and NormFinder were found. The combination of lj-miR171b and u1760353 was most stable in leaves using GeNorm, while lj-miR171b and u2100564 were the first and second most stable using NormFinder, respectively. In stem samples, the combination of u534122 and u3868172 was most stable using GeNorm, while u534122 and u4339213 were the first and second most stable using NormFinder, respectively. However, our result demonstrated that, in general, u534122 and u3868172 are the two most reliable reference genes.
Validation and Analysis of Candidate Reference miRNAs
To validate further the stability of u3868172 and u534122, we investigated their expression levels in different varieties of L. japonica (n = 8) collected from 21 geographic locations of China. This collection was also the basis for studying regulation mechanisms of secondary metabolism in both varieties of L. japonica in response to multiple environmental factors (Table 2). The expression levels of u3868172 and u534122 were also analyzed in three organs of these samples. The Ct value variance ranges of u3868172 were 1.33, 1.26, and 1.45 in flower buds, leaves, and stems (Figure 5A), respectively. Similarly, u534122 also showed relatively stable expression with Ct value variance ranges of 1.57 in flower buds, 1.62 in leaves, and 1.83 in stems (Figure 5B). Furthermore, the least expression variation in the different samples was shown by the combination of u3868172 and u534122, with Ct variance value ranges of 0.862 in flower buds, 0.878 in leaves, and 0.865 in stems (Figure 5C). In addition, within different samples, u3868172/u534122 had standard deviations of less than 1 (Figure 5C). Moreover, the expression levels of u3868172 and u534122 were also analyzed in three organs in all samples to validate their usage in miRNA regulation in these tissues (Figure 6). The Ct value variance range of u3868172/u534122 was the narrowest among the candidate miRNAs (Figures 5 and 6), which indicated that the combination u3868172/u534122 could be recommended as the optimal reference for miRNAs in L. japonica.
FIGURE 5. Expression levels of selected candidate reference miRNAs among different varieties. qRT-PCR analysis was conducted in L. japonica derived from different varieties in different production areas (n = 8). (A) Expression levels of u3868172 among different varieties. (B) Expression levels of u534122 among different varieties. (C) Expression levels of u3868172/u534122 among different varieties. Values are given as real-time PCR cycle threshold numbers (Ct values). The solid line represents the median value, and the boxes are 25th and 75th percentiles. The average is indicated by the point in the box. Whisker caps represent the maximum and minimum Ct values. (a = HN-DMH, b = FQ-DMH, c = CQ-SDYZ, d = SD-YTLZ, e = SD-DMH, f = GX-SDYZ, g = AH-SXYZ, h = HUB-HBYZ, i = CQ-SXYZ, j = SX-JHSH, k = NX-SDYZ, l = HB-DMH, m = YN-YT, n = JS-YT, o = BJ-YTLB, p = GS-YT, q = JS-HYH, r = BJ-YTH, s = HB-HYH, t = SD-YTH, u = GS-YTH).
FIGURE 6. Expression levels of selected candidate reference miRNAs among different tissues. qRT-PCR analysis was conducted in L. aponica derived from all samples (n = 21 × 8).
The quality of herbal medicine has been very difficult to control and to evaluate primarily because of the complexity and incomplete knowledge of the active medicinal compounds. L. japonica, widely used in traditional medicine, contains many active compound, whose contents differed significantly among the organs, varieties, and production areas of the herbal medicines. Thus, recent focus on variation in quality of medicinal material has been directed toward understanding the molecular regulatory mechanisms of secondary metabolism through transcriptomics or functional genomics approaches. miRNA genes are an important class of fine-tuning regulators, playing a major role in a wide range of developmental, biological, and metabolic processes in plants, including metabolism, stress response, vegetative phase changes, organogenesis, and signal transduction (Mathiyalagan et al., 2013). To date, a large number of miRNA families have been described in plants, some of which are directly associated with secondary metabolism synthesis (Ekström et al., 2015). Increasingly, miRNAs have been identified and characterized in some important medicinal herbs (Song et al., 2009; Wu et al., 2012b) and different tissues of medicinal plants (Xu et al., 2014; Khaldun et al., 2015). Further, the prediction of phenylpropanoid and terpenoid biosynthetic pathway genes targeting miRNAs in various plant species (Shi et al., 2010; Mathiyalagan et al., 2013), as well as evidence of flower development (Song et al., 2015; Niu et al., 2016) and stress resistance (Navarro et al., 2006; Lu et al., 2008; Zhang Y. et al., 2016), miRNAs produce insight into the study of miRNAs in different tissues and varieties of L. japonica.
qRT-PCR has emerged as a useful technique for validation of transcriptomics data owing to its precision, accuracy, convenience, speed, and sensitivity. However, accurate normalization of gene expression remains a major criterion for precise qRT-PCR analysis, as normalization helps in adjusting variation introduced at various steps of qRT-PCR arising from sample-to-sample variation, variation in RNA integrity, PCR efficiency, and cDNA sample loading (Shivhare and Lata, 2016). To date, numerous RNA species, including rRNA, snRNA, and synthetic miRNAs as spike-in controls, have been used as normalizers (Peltier and Latham, 2008). However, few studies have focused on selecting suitable reference miRNAs for different varieties and different tissues in medicinal plants. Thus, the present study performed comprehensive analysis of candidate reference genes for miRNA gene expression studies in different tissues of L. japonica varieties. Fourteen L. japonica miRNA candidate genes were selected for the comprehensive expression stability analysis from three small RNA-seq libraries. The performance of each candidate gene was tested in different samples (flower buds, leaves, and stems) to help determine suitable sets of reference genes for expression profiles studies in L. japonica.
In this research, we provide evidence that miRNAs could have high expression stability in medicinal plants, which has also been recently demonstrated in other plant species (Kulcheski et al., 2010; Hao et al., 2012; Borowski et al., 2014). The commonly used GeNorm, NormFinder, and BestKeeper algorithms identified the best suitable reference gene sets for each group of tissue samples in L. japonica. Although the results from the three algorithms did not rank genes in exactly the same order, these ranking discrepancies were previously reported and reflect differences in the calculation methods of these three approaches. u3868172 and u534122 appeared to be the two most optimal reference genes after analysis with GeNorm; a similar result was confirmed through NormFinder and BestKeeper analyses. After verification of the expression of u3868172 and u534122, the individual expressions of u3868172 and u534122 showed statistical differences within breeding groups, indicating that the selection of a single reference miRNA may be not adequate. Further, the combination of u3868172/u534122 was more suitable as the inner reference gene and was not affected by differences in breeding or plant organs. Additionally, secondary structures of u3868172 and u534122 precursors were predicted by PMRD, whose target genes could be investigated by available L. japonica transcriptome data in the future. Therefore, the group of miRNA reference genes in this study was validated not only as biomarkers for L. japonica but also for performing normalization of the expression of a new miRNAs specific to L. japonica.
Finally, we hope that our work will facilitate the exploration of regulatory mechanisms related to secondary metabolism for production of useful compounds in L. japonica.
The work presented here was carried out in collaboration between all authors. Completed the main experiment and draft the manuscript: YW and JL. Defined the research theme: YY and LH. Completed sequencing the small RNA libraries data: XW and GW. Completed the Data analysis: JZ. Completed validation of candidate reference miRNAs: SL and TC. Revised the manuscript: CJ and LZ.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
This work was supported by the Natural Science Foundation of China (81373959).
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fpls.2016.01101
FIGURE S1 | Some Lonicera japonica samples derived from different varieties in different producing areas (A) FLJ in Beijing; (B) FLJ in Beijing; (C) FLJ in Jiangsu (D) FLJ in GS; (E) FLJ in SD; (F) FLJ in Hubei; (G) FLJ in Ningxia; (H) FLJ in Chongqing; (I) rFLJ in Beijing; (J) rFLJ in Hebei; (K) rFLJ in Jiangsu; (L) rFLJ in GS.
FIGURE S2 | Linear correlation between Ct values and log concentration of miRNAs in a serially diluted miRNA plasmid standards. miRNAs expression levels were measured in 10-fold serial dilutions.
FIGURE S3 | Melting curves for sixteen candidate reference genes. Temperature is displayed in the x-axis, and the fluorescence signal is displayed in the y-axis.
Arenas-Huertero, C., Perez, B., Rabanal, F., Blanco-Melo, D., De la Rosa, C., Estrada-Navarrete, G., et al. (2009). Conserved and novel miRNAs in the legume Phaseolus vulgaris in response to stress. Plant Mol. Biol. 70, 385–401. doi: 10.1007/s11103-009-9480-3
Barakat, A., Wall, K., Leebens-Mack, J., Wang, Y. J., Carlson, J. E., and dePamphilis, C. W. (2007). Large-scale identification of microRNAs from a basal eudicot (Eschscholzia californica) and conservation in flowering plants. Plant J. 51, 991–1003. doi: 10.1111/j.1365-313X.2007.03197.x
Borowski, J. M., Galli, V., Messias, R. D. S., Perin, E. C., Buss, J. H., Silva, S. D. D. A. E., et al. (2014). Selection of candidate reference genes for real-time PCR studies in lettuce under abiotic stresses. Planta 239, 1187–1200. doi: 10.1007/s00425-014-2041-2
De Luis, A., Markmann, K., Coqnat, V., Holt, D. B., Charpentier, M., Parniske, M., et al. (2012). Two microRNAs linked to nodule infection and nitrogen-fixing ability in the legume Lotus japonicus. Plant Physiol. 160, 2137–2154. doi: 10.1104/pp.112. 204883
Devers, E. A., Branscheid, A., May, P., and Krajinski, F. (2011). Stars and symbiosis: microRNA- and microRNA∗-mediated transcript cleavage involved in arbuscular mycorrhizal symbiosis. Plant Physiol. 156, 1990–2010. doi: 10.1104/pp.111.172627
Ekström, L., Skilving, I., Ovesjö, M. L., Aklillu, E., Nylén, H., Rane, A., et al. (2015). miRNA-27b levels are associated with CYP3A activity in vitro and in vivo. Pharmacol. Res. Perspect. 3:e00192. doi: 10.1002/prp2.192
Fu, J., Yang, L., Khan, M. A., and Mei, Z. (2013). Genetic characterization and authentication of Lonicera japonica Thunb. by using improved RAPD analysis. Mol. Biol. Rep. 40, 5993–5999. doi: 10.1007/s11033-013-2703-3
Gao, Z. H., Wei, J. H., Yang, Y., Zhang, Z., and Zhao, W. T. (2012). Selection and validation of reference genes for studying stress-related agarwood formation of Aquilaria sinensis. Plant Cell Rep. 31, 1759–1768. doi: 10.1007/s00299-012-1289-x
Hao, F., Huang, X., Zhang, Q., Wei, G., Wang, X., and Kang, Z. (2012). Selection of suitable inner reference genes for relative quantification expression of microRNA in wheat. Plant Physiol. Biochem. 51, 116–122. doi: 10.1016/j.plaphy.2011.10.010
Hong, S. Y., Seo, P. J., Yang, M. S., Xiang, F., and Park, C. M. (2008). Exploring valid reference genes for gene expression studies in Brachypodium distachyon by real-time PCR. BMC Plant Biol. 8:112. doi: 10.1186/1471-2229-8-112
Jiang, J., Lv, M., Liang, Y., Ma, Z., and Cao, J. (2014). Identification of novel and conserved miRNAs involved in pollen development in Brassica campestris ssp. chinensis by high-throughput sequencing and degradome analysis. BMC Genomics 15:146. doi: 10.1186/1471-2164-15-146
Jiang, K., Pi, Y., Hou, R., Zeng, H., Huang, Z., Zhang, Z., et al. (2009). Molecular cloning and expression profiling of the first specific jasmonate biosynthetic pathway gene allene oxide synthase from Lonicera japonica. Mol. Biol. Rep. 36, 487–493. doi: 10.1007/s11033-007-9205-0
Kang, H. N., Chung, H. J., Jeon, E. J., Mi, K. P., Yong, H. Y., Liu, J. R., et al. (2007). In vitro, biosynthesis of strictosidine using Lonicera japonica, leaf extracts and recombinant yeast. J. Plant Biol. 50, 315–320. doi: 10.1007/BF03030660
Khaldun, A. B. M., Huang, W., Liao, S., Lv, H., and Wang, Y. (2015). Identification of microRNAs and target genes in the fruit and shoot tip of Lycium chinense: a traditional Chinese medicinal plant. PLoS ONE 10:e0116334. doi: 10.1371/journal.pone.0116334
Kou, S. J., Wu, X. M., Liu, Z., Liu, Y. L., Xu, Q., and Guo, W. W. (2012). Selection and validation of suitable reference genes for miRNA expression normalization by quantitative RT-PCR in citrus somatic embryogenic and adult tissues. Plant Cell Rep. 31, 2151–2163. doi: 10.1007/s00299-012-1325-x
Kulcheski, F. R., Nepomuceno, A. L., Abdelnoor, R. V., and Margis, R. (2010). The use of microRNAs as reference genes for quantitative polymerase chain reaction in soybean. Anal. Biochem. 406, 185–192. doi: 10.1016/j.ab.2010.07.020
Kwak, W. J., Cho, Y. B., Han, C. K., Shin, H. J., Ryu, K. H., Yoo, H., et al. (2005). Extraction and purification method of active constituents from stem of i Lonicera japonica /i Thunb., its usage for anti-inflammatory and analgesic drug. EP. E1536811. Patent US20060014240. Washington, DC: US Patent and Trademark office.
Lelandais-Brière, C., Nayaa, L., Sallet, E., Calenge, F., Frugier, F., Hartmann, C., et al. (2009). Genome-wide Medicago truncatula small RNA analysis revealed novel microRNAs and isoforms differentially regulated in roots and nodules. Plant Cell 21, 2780–2796. doi: 10.1105/tpc. 109.068130
Li, W., Cheng, Z., Wang, Y., and Qu, H. (2013). Quality control of Lonicerae japonicae Flos using near infrared spectroscopy and chemometrics. J. Pharm. Biomed. Anal. 72, 33–39. doi: 10.1016/j.jpba.2012.09.012
Lin, L., Wang, P., Du, Z., Wang, W., Cong, Q., Zheng, C., et al. (2016). Structural elucidation of a pectin from flowers of Lonicera japonica and its antipancreatic cancer activity. Int. J. Biol. Macromol. 88, 130–137. doi: 10.1016/j.ijbiomac.2016.03.025
Liu, D., Shi, L., Han, C., Yu, J., Li, D., and Zhang, Y. (2012). Validation of reference genes for gene expression studies in virus-infected Nicotiana benthamiana using quantitative real-time PCR. PLoS ONE 7:e46451. doi: 10.1371/journal.pone.0046451
Liu, Z., Cheng, Z., He, Q., Lin, B., Gao, P., Li, L., et al. (2016). Secondary metabolites from the flower buds of Lonicera japonica and their in vitro anti-diabetic activities. Fitoterapia 110, 44–51. doi: 10.1016/j.fitote.2016.02.011
Luo, X., Shi, T., Sun, H., Song, J., Ni, Z., and Gao, Z. (2014). Selection of suitable inner reference genes for normalisation of microRNA expression response to abiotic stresses by RT-qPCR in leaves, flowers and young stems of peach. Sci. Hortic. 165, 281–287. doi: 10.1016/j.scienta.2013.10.030
Mathiyalagan, R., Subramaniyam, S., Natarajan, S., Kim, Y. J., Sun, M. S., Kim, S. Y., et al. (2013). Insilico profiling of microRNAs in Korean ginseng (Panax ginseng Meyer). J. Ginseng Res. 37, 227–247. doi: 10.5142/jgr.2013.37.227
Moldovan, D., Spriggs, A., Yang, J., Pogson, B. J., Dennis, E. S., and Wilson, I. W. (2010). Hypoxia-responsive microRNAs and trans-acting small interfering RNAs in Arabidopsis. J. Exp. Bot. 61, 165–177. doi: 10.1093/jxb/erp296
Navarro, L., Dunoyer, P., Jay, F., Arnold, B., Dharmasiri, N., Estelle, M., et al. (2006). A plant miRNA contributes to antibacterial resistance by repressing auxin signaling. Science 21, 436–439. doi: 10.1126/science.aae0382
Niu, Q., Li, J., Cai, D., Qian, M., Jia, H., Bai, S., et al. (2016). Dormancy-associated MADS-box genes and microRNAs jointly control dormancy transition in pear (Pyrus pyrifolia white pear group) flower bud. J. Exp. Bot. 67, 239–257. doi: 10.1093/jxb/erv 454
Peltier, H., and Latham, G. (2008). Normalization of microrna expression levels in quantitative RT-PCR assays: identification of suitable reference RNA targets in normal and cancerous human solid tissues. RNA 14, 844–852. doi: 10.1261/rna.939908
Pfaﬄ, M. W., Tichopad, A., Prgomet, C., and Neuvians, T. P. (2004). Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: bestkeeper-excel-based tool using pair-wise correlations. Biotechnol. Lett. 26, 509–515. doi: 10.1023/B:BILE.0000019559.84305.47
Reid, K. E., Olsson, N., and Schlosser, J. (2006). An optimized grapevine RNA isolation procedure and statistical determination of reference genes for real-time RT-PCR during berry development. BMC Plant Biol. 6:27. doi: 10.1186/1471-2229-6-27
Shi, R., Yang, C., Lu, S., Sederoff, R., and Chiang, V. L. (2010). Specific down-regulation of PAL genes by artificial microRNAs in Populus trichocarpa. Planta 232, 1281–1288. doi: 10.1007/s00425-010-1253-3
Shi, Z., Liu, Z., Liu, C., Wu, M., Su, H., Ma, X., et al. (2016). Spectrum-effect relationships between chemical fingerprints and antibacterial effects of Lonicerae japonicae Flos and Lonicerae Flos base on UPLC and microcalorimetry. Front. Pharmacol. 7:12. doi: 10.3389/fpha.2016.00012
Shivhare, R., and Lata, C. (2016). Selection of suitable reference genes for assessing gene expression in pearl millet under different abiotic stresses and their combinations. Sci. Rep. 6:23036. doi: 10.1038/srep23036
Song, J. H., Yang, J., Pan, F., and Jin, B. (2015). Differential expression of microRNAs may regulate pollen development in Brassica oleracea. Genet. Mol. Res. 14, 15024–15034. doi: 10.4238/2015.November.24.10
Takle, G. W., Toth, I. K., and Brurberg, M. B. (2007). Evaluation of reference genes for real-time RT-PCR expression studies in the plant pathogen Pectobacterium atrosepticum. BMD Plant Biol. 7:50. doi: 10.1186/1471-2229-7-50
Tan, T., Lai, C. J. S., Yang, H. O., He, M. Z., and Feng, Y. (2016). Ionic liquid-based ultrasound-assisted extraction and aqueous two-phase system for analysis of caffeoylquinic acids from Flos Lonicerae Japonicae. J. Pharm. Biomed. Anal. 120, 134–141. doi: 10.1016/j.jpba.2015.12.019
Teste, M. A., Duquenne, M., François, J. M., and Parrou, J. L. (2009). Validation of reference genes for quantitative expression analysis by real-time RT-PCR in Saccharomyces cerevisiae. BMC Mol. Biol. 10:20. doi: 10.1186/1471-2199-10-99
Tong, Z., Gao, Z., Fei, W., Zhou, J., and Zhen, Z. (2009). Selection of reliable reference genes for gene expression studies in peach using real-time PCR. BMC Mol. Biol. 10:71. doi: 10.1186/1471-2199-10-71
Vandesompele, J., Preter, K. D., Pattyn, F., Poppe, B., Roy, N. V., Paepe, A. D., et al. (2002). Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 3:research0034. doi: 10.1186/gb-2002-3-7-research0034
Wang, T., Lu, J., Xu, Z., and Zhang, Q. (2014). Selection of suitable reference genes for miRNA expression normalization by qRT-PCR during flower development and different genotypes of Prunus mume. Sci. Hortic. 169, 130–137. doi: 10.1016/j.scienta.2014.02.006
Wu, B., Li, Y., Yan, H., Ma, Y., Luo, H., Yuan, L., et al. (2012a). Comprehensive transcriptome analysis reveals novel genes involved in cardiac glycoside biosynthesis and mlncRNAs associated with secondary metabolism and stress response in Digitalis purpurea. BMC Genomics 13:15. doi: 10.1186/1471-2164-13-15
Wu, B., Wang, M., Ma, Y., Yuan, L., and Lu, S. (2012b). High-throughput sequencing and characterization of the small RNA transcriptome reveal features of novel and conserved microRNAs in Panax ginseng. PLoS ONE 7:e44385. doi: 10.1371/journal.pone.0044385
Wu, J., Su, S., Fu, L., Zhang, Y., Chai, L., and Yi, H. (2014). Selection of reliable reference genes for gene expression studies using quantitative real-time PCR in navel orange fruit development and pummelo floral organs. Sci. Hortic. 176, 180–188. doi: 10.1016/j.scienta.2014.06.040
Wu, J., Wang, X. C., Liu, Y., Du, H., Shu, Q. Y., Su, S., et al. (2016). Flavone synthases from Lonicera japonica and L. macranthoides reveal differential flavone accumulation. Sci. Rep. 6:19245. doi: 10.1038/srep19245
Xu, X., Jiang, Q., Ma, X., Ying, Q., Shen, B., Qian, Y., et al. (2014). Deep sequencing identifies tissue-specific microRNAs and their target genes involving in the biosynthesis of tanshinones in Salvia miltiorrhiza. PLoS ONE 9:e111679. doi: 10.1371/journal.pone.0111679
Xue, W., Wang, Z., Du, M., Liu, Y., and Liu, J. Y. (2013). Genome-wide analysis of small RNAs reveals eight fiber elongation-related and 257 novel microRNAs in elongating cotton fiber cells. BMC Genomics 14:629. doi: 10.1186/1471-2164-14-629
Yuan, Y., Song, L., Li, M., Liu, G., Chu, Y., Ma, L., et al. (2012). Genetic variation and metabolic pathway intricacy govern the active compound content and quality of the Chinese medicinal plant Lonicera japonica thunb. BMC Genomics 13:195. doi: 10.1186/1471-2164-13-195
Zhai, J., Jeong, D. H., De Paoli, E., Park, S., Rosen, B. D., Li, Y., et al. (2011). MicroRNAs as master regulators of the plant NB-LRR defense gene family via the production of phased, trans-acting siRNAs. Genes Dev. 25, 2540–2553. doi: 10.1101/gad.177527.111
Zhang, B., Nan, T., Zhan, Z., Kang, L., Yang, J., Yuan, Y., et al. (2016). Development of a monoclonal antibody-based enzyme-linked immunosorbent assay for luteoloside detection in Flos Lonicerae Japonicae. Anal. Bioanal. Chem. doi: 10.1007/s00216-016-9396-0 [Epub ahead of print].
Zhang, Y., Wang, Y., Xie, F., Li, C., Zhang, B., Nichols, R. L., et al. (2016). Identification and characterization of microRNAs in the plant parasitic root-knot nematode Meloidogyne incognita using deep sequencing. Funct. Integr. Genomics 16, 127–142. doi: 10.1007/s10142-015-0472-x
Zhang, B. H., Pan, X. P., Wang, Q. L., Cobb, G. P., and Anderson, T. A. (2005). Identification and characterization of new plant microRNAs using EST analysis. Cell Res. 15, 336–360. doi: 10.1038/sj.cr.7290302
Zhang, L., Chia, J. M., Kumari, S., Stein, J. C., Liu, Z., Narechania, A., et al. (2009). A genome-wide characterization of microrna genes in maize. PloS Genet. 5:e1000716. doi: 10.1371/journal.pgen.1000716
Zhao, Y. T., Wang, M., Fu, S. X., Yang, W. C., Qi, C. K., and Wang, X. J. (2012). Small RNA profiling in two Brassica napus cultivars identifies microRNAs with oil production- and development-correlated expression and new small RNA classes. Plant Physiol. 158, 813–823. doi: 10.1104/pp.111.187666
Keywords: Lonicera japonica, microRNAs, qRT-PCR, reference gene, normalization, validation, different varieties
Citation: Wang Y, Liu J, Wang X, Liu S, Wang G, Zhou J, Yuan Y, Chen T, Jiang C, Zha L and Huang L (2016) Validation of Suitable Reference Genes for Assessing Gene Expression of MicroRNAs in Lonicera japonica. Front. Plant Sci. 7:1101. doi: 10.3389/fpls.2016.01101
Received: 08 April 2016; Accepted: 11 July 2016;
Published: 26 July 2016.
Edited by:Mohammad Anwar Hossain, Bangladesh Agricultural University, Bangladesh
Reviewed by:Shri Ram Yadav, Indian Institute of Technology, Roorkee, India
Ai-Sheng Xiong, Nanjing Agricultural University, China
Zhensheng Kang, Northwest A&F University, China
Copyright © 2016 Wang, Liu, Wang, Liu, Wang, Zhou, Yuan, Chen, Jiang, Zha and Huang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
†These authors have contributed equally to this work.