Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Normalization for Relative Quantification of mRNA and microRNA in Soybean Exposed to Various Abiotic Stresses

  • Weican Liu ,

    Contributed equally to this work with: Weican Liu, Yu Deng

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • Yu Deng ,

    Contributed equally to this work with: Weican Liu, Yu Deng

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • Yonggang Zhou,

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • Huan Chen,

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • Yuanyuan Dong,

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • Nan Wang,

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • Xiaowei Li,

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • Aysha Jameel,

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • He Yang,

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • Min Zhang,

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • Kai Chen,

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • Fawei Wang ,

    fw-1980@163.com (FW); hyli99@163.com (HL)

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

  • Haiyan Li

    fw-1980@163.com (FW); hyli99@163.com (HL)

    Affiliation College of Life Sciences, Engineering Research Center of the Chinese Ministry of Education for Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

Abstract

Plant microRNAs are small non-coding, endogenic RNA molecule (containing 20–24 nucleotides) produced from miRNA precursors (pri-miRNA and pre-miRNA). Evidence suggests that up and down regulation of the miRNA targets the mRNA genes involved in resistance against biotic and abiotic stresses. Reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) is a powerful technique to analyze variations in mRNA levels. Normalizing the data using reference genes is essential for the analysis of reliable RT-qPCR data. In this study, two groups of candidate reference mRNAs and miRNAs in soybean leaves and roots treated with various abiotic stresses (PEG-simulated drought, salinity, alkalinity, salinity+alkalinity, and abscisic acid) were analyzed by RT-qPCR. We analyzed the most appropriate reference mRNA/miRNAs using the geNorm, NormFinder, and BestKeeper algorithms. According to the results, Act and EF1b were the most suitable reference mRNAs in leaf and root samples, for mRNA and miRNA precursor data normalization. The most suitable reference miRNAs found in leaf and root samples were 166a and 167a for mature miRNA data normalization. Hence the best combinations of reference mRNAs for mRNA and miRNA precursor data normalization were EF1a + Act or EF1b + Act in leaf samples, and EF1a + EF1b or 60s + EF1b in root samples. For mature miRNA data normalization, the most suitable combinations of reference miRNAs were 166a + 167d in leaf samples, and 171a + 156a or 167a + 171a in root samples. We identified potential reference mRNA/miRNAs for accurate RT-qPCR data normalization for mature miRNA, miRNA precursors, and their targeted mRNAs. Our results promote miRNA-based studies on soybean plants exposed to abiotic stress conditions.

Introduction

MicroRNAs are endogenous and non-coding small RNAs that are generally found in plants. The mature miRNAs are approximately 20–24 nucleotides and are derived from miRNA precursors (i.e., pre-miRNA and pri-transcripts) by Dicer-like enzymatic digestion. The miRNAs play an important regulatory role at the post-transcriptional level by targeting mRNA cleavage or translation repression through the RNA-induced silencing complex, where they form complexes with mRNA[1]. Under stress conditions, plant miRNAs are up- or down-regulated by targeting specific stress-associated mRNAs to increase tolerance to adverse environmental conditions[24]. Identification and characterization of miRNAs, miRNA precursors, and their targeted mRNAs are essential for the analysis of the miRNA molecular regulatory mechanism.

RT-qPCR is considered the gold standard for the quantification of mRNA expression because of its accuracy, sensitivity, specificity, reproducibility, and robustness[57]. However, because miRNAs are small, RT-qPCR for miRNA detection is more complicated than traditional RT-qPCR methods for mRNA detection. To address this problem, several RT-qPCR-based procedures have been developed to detect miRNA. First, Schmittgen et al. [8, 9] recommended using the traditional RT-qPCR method to monitor the expression of miRNA precursors, including pri-miRNA and pre-miRNA containing the hairpin sequence. The RNA template was reverse transcribed to cDNA using either gene-specific primers or random hexamers and reverse transcriptase. To quantify samples, sense and antisense primers were designed based on the hairpin sequence of miRNA precursors. The pri-miRNA and pre-miRNA templates were simultaneously amplified. Second, Poly(A)-RT-qPCR[10] were developed to monitor mature miRNAs, in which the RNA templates, including miRNAs, were polyadenylated by polymerase and then reverse transcribed to cDNAs using poly(T) adapters. The designed miRNA-specific forward primer and the universal reverse primer complementary to the poly(T) adapter were used for RT-qPCR. The third method involved stem-loop RT-qPCR[9,1114]. The stem-loop reverse primers were designed with a universal backbone sequence to form a stem-loop structure due to the complementarity between the 5′ and 3′ ends. The specificity of the stem-loop reverse primers for an individual miRNA were based on the fact the last six nucleotides were the reverse complement of the six nucleotides at the 3′ end of the miRNA. The reverse-transcribed product was amplified using the miRNA-specific forward primers and the universal reverse primers. The three RT-qPCR methods to detect miRNAs have their own advantages and disadvantages respectively. Poly(A)-tailing and stem-loop methods are mainly used to detect mature miRNAs. They cannot distinguish between different miRNAs in a family if they produce the same mature miRNA. Therefore, to detect different types of miRNAs, the specific precursors must be analyzed[8]. Mou et al.[15] reported that stem-loop RT-qPCR and poly(A)-tailing RT-qPCR can detect highly abundant miRNAs at similar levels of accuracy and specificity. However, stem-loop RT-qPCR cannot detect low abundant miRNAs, suggesting poly(A)-tailing RT-qPCR may be a better option for miRNA analysis. In this study, we selected poly(A)-tailing RT-qPCR to detect mature miRNAs because the polyadenylated miRNA reverse transcription template enabled universal miRNA detection. This allowed the detection of numerous miRNAs using one template set. Another reason we used poly(A)-tailing RT-qPCR was that stem-loop RT-qPCR was more expensive (e.g., higher costs associated with specific primer synthesis). Moreover, miRNA precursor RT-qPCR provides complementary data for a more comprehensive characterization of miRNA expression.

There is a lack of consensus regarding how best to perform and interpret RT-qPCR experiments. Many technical defects can affect the accuracy of RT-qPCR analysis. Therefore, the minimum information required to publish RT-qPCR experiments (MIQE) [16] and the Real-Time PCR Data Markup Language (RDML) were published[17]. These standards for RT-qPCR experiments focus on the reliability and consistency of results. According to the MIQE, the selection of suitable reference genes is one of the essential components to ensure the accuracy of a RT-qPCR assay. This is because RT-qPCR results are influenced by factors such as variations in RNA extraction, reverse transcription, and efficiency of amplification. Thus, data normalization of RT-qPCR using reference genes enables comparisons of different mRNA concentrations from various samples [16]. Ideal reference genes should be stably expressed between different plant varieties, tissues, stages of development, and biotic or abiotic stress conditions[18, 19]. However, numerous reports have determined that candidate reference genes may be stably expressed only under certain conditions [2024]. It is necessary to identify and select suitable reference genes through specific experimental procedures [16]. Unfortunately, there is no universally accepted method for selecting reference genes and assessing stability. GeNorm[25], NormFinder[26], and BestKeeper[27] are three algorithms widely used to analyze the stability of candidate reference genes. Using different algorithms to evaluate the stability of reference gene expression may facilitate the selection of reliable reference genes for precise data normalization of RT-qPCR.

Soybean (Glycine max) is an important crop for seed protein and edible oil production. In our earlier study, the miRNAs and transcriptional profiles of genes associated with several stress responses were sequenced using deep sequencing technology [28, 29]. The objective of the previous study was to confirm the accuracy of deep sequencing technology and determine the variations in expression of mature miRNAs, miRNA precursors, and their targeted mRNAs in soybean under various abiotic stress conditions. In the current study, we aimed to identify suitable reference mRNA/miRNAs from the eight traditional candidate reference mRNA genes and eight candidate reference miRNAs. Soybean roots and leaves exposed to various abiotic stress conditions (PEG-simulated drought, salinity, alkalinity, salinity+alkalinity, and abscisic acid) that were used for relative quantification analysis. We used soybean cultivar ‘Williams 82’ of its known genome sequence [30]. Statistical analysis was completed using the geNorm, NormFinder, and BestKeeper algorithms. To the best of our knowledge, this study is the first to use the entire set of reference mRNA/miRNAs for accurate RT-qPCR data normalization for mature miRNAs, miRNA precursors, and their targeted mRNAs in soybean leaf or root tissues exposed to various abiotic stresses.

Materials and Methods

Plant materials and abiotic stress treatments

Soybean seeds (‘Williams 82’) were treated with ethanol for 10 min and then rinsed several times with sterile distilled water. The seeds were cultured in Hoagland nutrient solution and grown at 30°C with a 16 h light/8 h dark photoperiod (80 μmol m−2 s-1 photon flux density) and 50% relative humidity. When the first pair of unifoliate leaves fully opened, we initiated stress treatments as follows: PEG-simulated drought (8% PEG 8000), salinity (120 mM NaCl), alkalinity (100 mM NaHCO3), salinity+alkalinity (70 mM NaCl+50 mM NaHCO3), and ABA (200 μM ABA). Untreated plants were used as controls. The seedlings were incubated for 0, 1, 3, 6, 9, and 12 h, with leaf and root samples collected at each time point. Three biological replicates were performed for each stress treatment. All samples were immediately frozen in liquid nitrogen and stored at −80°C until required. All samples information above can be found in S1 Table.

RNA extraction and quality controls

Total RNA was extracted using RNAiso Plus (Takara, Japan) according to the manufacturer’s instructions. The purified RNA was analyzed using a Thermo Scientific NanoDrop2000 spectrophotometer. A 260 nm/280 nm optical density ratio of 1.8–2.0 indicated high quality RNA. The RNA integrity was checked by 1% agarose gel electrophoresis. Clearly visible RNA bands and a 25S/18S ratio close to 2:1 indicated intact RNA.

Primer designing and their validation

We selected eight candidate reference mRNA genes for data normalization of mRNAs and miRNA precursors during RT-qPCR analysis. We also selected eight candidate reference miRNAs for data normalization of mature miRNAs. The selection of the candidate reference mRNA/miRNAs was based on previous reports [20, 21, 3139]. The mRNA sequences were obtained from the Phytozome 10.3 website [40] (http://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Gmax) and the miRNA sequences were obtained from the miRBase(Release 21) website [41] (http://www.mirbase.org/). Primers were designed using DNAstar software[42]. The mRNA primers were designed possibly across two exons while the miRNA precursor primers were designed at stem-loop regions [8, 43]. Primers were designed to amplify mRNA and miRNA precursor products of 80–200 bases. The possible secondary structures were assessed using Mfold[44] (http://unafold.rna.albany.edu/?q=mfold) and the amplicon specificity of the primers was analyzed using BLAST tools (http://blast.ncbi.nlm.nih.gov/Blast.cgi)[45]. Additionally, the miRNA-specific forward primers were designed using the NCode VILO miRNA cDNA Synthesis Kit (Invitrogen) following the manufacturer’s protocol. The miRNA universal reverse primers were provided in the kit. The designed primer sequences were listed in Table 1. All primers were synthesized at the Suzhou GENEWIZ Biological Technology Services Company, China.

thumbnail
Table 1. Descriptions of mRNAs and miRNAs and RT-qPCR amplification characteristics.

https://doi.org/10.1371/journal.pone.0155606.t001

Robust and precise RT-qPCR assays are usually correlated with high PCR efficiency[16]. The PCR amplification efficiency should be calculated using the formula E = 10(−1/Slope)-1 and the slope of the calibration curve. Calibration curves were produced using mixed cDNA samples as a starting template and five or six replicates of 10-, 4-, and 2-fold serial dilutions to create a gradient of concentrations. The Mx3000P RT-qPCR system automatically generated calibration curves. Only primers with PCR amplification efficiencies close to 100% were used. Furthermore, the amplicon specificity for all candidate reference mRNA/miRNAs was confirmed by the presence of a single peak.

cDNA synthesis and qPCR

The mRNA and miRNA precursor expression levels were determined using standard RT-qPCR [7]. Following treatment with gDNA Eraser, 1 μg total RNA was used as the template for first-strand cDNA synthesis with the PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Japan). Reverse transcription was completed with oligo(dT) and random hexamer primers. The cDNA was diluted 10-fold for mRNA and miRNA precursor analysis. The expression of mature miRNA was measured using Poly(A)-RT-qPCR[10]. Total RNA including miRNA was reverse transcribed using the NCode VILO miRNA cDNA Synthesis Kit according to the manufacturer’s protocol. The miRNA was polyadenylated using ATP and poly-A polymerase, and cDNA was synthesized by reverse transcribing the tailed RNA using the universal RT Primer in a single reaction volume of 20 μl. We diluted the cDNA 10-fold for miRNA analysis stored at −20°C until required.

The RT-qPCR was completed in 96-well blocks using an Applied Mx3000P Real-Time Thermocycler (Stratagene). The reaction mixtures were prepared using the SYBR Premix Ex Taq II kit (Takara, Japan). The primer sequences were listed in Table 1. The RT-qPCR program consisted of an initial step at 95°C for 2 min to activate the Taq DNA polymerase, followed by 40 cycles of 94°C for 5 s and 62°C for 20s. We verified the amplicons specificity by melting curve analysis from 62°C to 95°C. Each reaction was completed with three technical replicates. However, three biological replicates have been pooled together after cDNA synthesis for expression stability assessment of candidate reference mRNA/miRNAs. All candidate reference mRNA genes or candidate reference miRNAs were quantified using the same batch of cDNA to minimize experimental variation. The no-template controls were included in each sample batch.

Analysis of gene expression stability

We initially examined the expression stability of eight candidate reference mRNA genes and eight candidate reference miRNAs using boxplot analysis[46]. Boxplots of the quantification cycle (Cq) values for each candidate reference mRNA/miRNA were assessed using MATLAB software. The expression variation range of Cq values (ΔCq) was calculated using the formula ΔCq = Cqmax − Cqmin. Then we used geNorm[25] analysis (https://genorm.cmgg.be/). The sample with the highest expression level (i.e., with the lowest Cq value) for each gene was used as a control with an expression level of 1. The relative expression levels for the other samples were calculated using the formula 2−(average Cq samples − average Cqmin). The obtained data were analyzed to determine the gene stability value (M) for each candidate reference mRNA/miRNA. Additionally, the average pairwise variation (Vn/n+1) was calculated to determine the optimal number of reference mRNA/miRNAs. Next we analyzed the samples using NormFinder [26] (http://moma.dk/normfinder-software). The raw average Cq values were converted to linear scale expression quantities using a standard curve. The expression stability (S) values and best gene combination between inter-group were assessed by NormFinder. We also completed BestKeeper[27] analysis (http://gene-quantification.com/bestkeeper.html), which was based on average raw Cq values of the candidate reference mRNA/miRNAs. The pair-wise correlation assessments of all candidate mRNA/miRNAs pairs were conducted and geometric expression means were calculated. BestKeeper algorithm identified the most stably expressed gene based on the highest correlation coefficient. Finally, the consensus rankings of the candidate reference mRNA/miRNAs were determined according to the geometric rank means from three analyses (i.e., geNorm, NormFinder, and BestKeeper).

Validation of gene expression stability

To confirm the identities of the most stable candidate reference mRNA/miRNAs and the best reference combinations, the relative expression levels of Growth Regulating Factor 9 (GRF9), precursor gma-miR396a (Pre-396a), and mature gma-miR396a-5p (396a) in leaf and root samples treated for 3 h under various abiotic stresses(PEG-simulated drought, salinity, alkalinity, salinity+alkalinity, and abscisic acid)were determined using RT-qPCR. Each reaction was completed with three technical and three biological replicates. The primer details are listed in Table 1. The relative expression levels were calculated using the 2−ΔΔCt formula[5]. The significance of the expression level differences was determined using one-way analysis of variance in SPSS 13.0.

Results

RNA quality and characteristics of RT-qPCR amplified products

RNA quality, primer specificity, and amplification efficiency are key issues related to RT-qPCR according to the MIQE[16] guidelines. The ribosomal RNA bands observed during 1% agarose gel electrophoresis were clearly visible and the 25S:18S ratio was approximately 2:1 (S1A Fig). The 260 nm/280 nm optical density ratio ranged from 1.8 to 2.0 (S1B Fig), indicating the RNA quality was appropriate for RT-qPCR assays.

The candidate reference mRNA/miRNAs were selected based on their amplification efficiency and specificity. The gene-related information and RT-qPCR amplified product characteristics were listed in Table 1. The primer pair locations on transcript sequences are indicated in S2 File. Primer specificity was confirmed based on the dissociation curves of the amplicons consisted of a single peak (S1 File). The amplification efficiency for the tested gene primers varied from 94.5 to 108.8% (Table 1), and the correlation coefficients of standard curves (S1 File) were between 0.991 and 1.000 (Table 1).

Candidate reference mRNA/miRNA expression profiles

The boxplots of Cq values for each candidate reference mRNA/miRNA were assessed in sets of samples respectively (Fig 1). The sets of samples were listed in S1 Table. The Cq values for each candidate reference mRNA/miRNA were ranged in 18–30 cycles, which corresponded to a wide dynamic normalization range (Fig 1).

thumbnail
Fig 1. Boxplots of the quantification cycle (Cq) values of candidate reference mRNAs and candidate reference miRNAs.

Boxplots of Cq values for each candidate reference mRNA/miRNAs were assessed in sets of samples, respectively(S1 Table). (A) Candidate reference mRNA genes were analyzed in leaf and root combined samples (n = 52). (B) Candidate reference mRNA genes were analyzed in leaf samples (n = 26). (C) Candidate reference mRNA genes were analyzed in root samples (n = 26). (D) Candidate reference miRNAs were analyzed in leaf and root combined samples (n = 52). (E) Candidate reference miRNAs were analyzed in leaf samples (n = 26). (F) Candidate reference miRNAs were analyzed in root samples (n = 26). The box indicates the 25th and 75th percentiles. The line across the box corresponds to the median, and whiskers represent the maximum and minimum values. The plus sign indicates the maximum or minimum outlier.

https://doi.org/10.1371/journal.pone.0155606.g001

Because the cDNA templates for the different samples were reverse transcribed from the same amount of total RNA, we assumed that ΔCq was related to the expression stability of candidate reference mRNA/miRNAs. A narrow range of Cq values indicated more stable expression. The ΔCq for each candidate reference mRNA gene in leaf and root combined samples ranged from 2.35 to 3.57 cycles (Fig 1A), while the ΔCq for each candidate reference miRNA in leaf and root combined samples ranged from 1.785 to 6.67 cycles (Fig 1D). The ΔCq values indicated a wider inter-tissue range. The ΔCq for each candidate reference mRNA gene in leaf samples ranged from 0.95 (EF1a) to 1.75 (TuA) cycles (Fig 1B), while in root samples it ranged from 1.185 (EF1a) to 2.29 (TuB) cycles (Fig 1C). The ΔCq for each candidate reference miRNA in leaf samples ranged from 0.69 (156a) to 1.955 (393a) cycles (Fig 1E), while in root samples it ranged from 0.785 (156a) to 2.15 (172a) cycles (Fig 1F). The ΔCq values exhibited a narrow intra-tissue range. All results indicated that the intra-tissue candidate reference mRNA/miRNAs expression was more stable than the inter-tissue expression. However, comparisons of ΔCq were insufficient, which required further assessments using three statistical algorithms to validate the expression stability of the candidate reference mRNA/miRNAs.

GeNorm assessment of expression stability

The geNorm algorithm eliminates the least stable reference gene and then recalculates the average M values for the remaining candidate reference genes. The genes with smaller M values were the more stably expressed genes [25]. Additionally, the optimal number of reference genes was determined based on pairwise variation (Vn/n+1). If Vn/n+1 > 0.15, an additional (n + 1) reference gene was required, but if Vn/n+1 ≤ 0.15, an additional (n + 1) reference gene was not required[25]. The assessment of candidate reference mRNA/miRNAs expression stability was completed in sets of samples, respectively (same as in Fig 1A to 1F) using geNorm. In stepwise analysis, we observed that the M values of the most stable candidate reference mRNAs and candidate reference miRNAs in leaf and root combined samples were bigger than those of leaf or root samples with M values of EF1a/Act (0.33) > EF1a/Act (0.17) or EF1a/EF1b (0.16) (Fig 2A, 2B and 2C) and 172a/393a (0.35) > 166a/167a (0.17) or 156a/171a (0.21) (Fig 2D, 2E and 2F). Meanwhile, the Vn/n+1 values of the candidate reference mRNAs and candidate reference miRNAs for leaf and root combined samples all were greater than those of leaf or root samples (Fig 3), the data demonstrated that the best combinations of reference mRNA/miRNAs was better in leaf or root tissues than in leaf and root combined samples. More specifically, the Vn/n+1 values of the candidate reference miRNAs for leaf and root combined samples all were higher than the threshold value of 0.15 (Fig 3), it indicated we could not get suitable combination of reference miRNAs for leaf and root combined samples. This results was consistent with the results of the ΔCq analysis. The intra-tissue candidate reference mRNA/miRNAs expression was more stable than the inter-tissue expression. Therefore, we subsequently analyzed the gene expression stability of candidate reference mRNA/miRNAs in leaf and root samples.

thumbnail
Fig 2. Ranking of the candidate reference mRNAs and candidate reference miRNAs by geNorm.

The ranking of candidate reference mRNA/ miRNAs was based on average expression stability values (M). A lower M value indicates a more stable expression level. The six tested groups are the same as those indicated in Fig 1 (Group A to F).

https://doi.org/10.1371/journal.pone.0155606.g002

thumbnail
Fig 3. Determination of the optimal number of reference mRNAs and candidate reference miRNAs for normalization according to pairwise variation (Vn/n+1) using geNorm.

Vn/n+1 > 0.15 indicated an additional (n + 1) reference was required while Vn/n+1 ≤ 0.15 indicated an additional reference was not required. The six tested groups are the same as those indicated in Fig 1 (Group A to F).

https://doi.org/10.1371/journal.pone.0155606.g003

The geNorm gene expression stability results for candidate reference mRNA/miRNAs in leaf and root samples are provided in Table 2 and Fig 2. Regarding candidate reference mRNA genes for mRNA and miRNA precursor data normalization, the most stable reference mRNA (i.e., with the lowest M values) in leaf and root samples were EF1a/Act (M value: 0.17) and EF1a/EF1b (M value: 0.16), respectively. The most unstable reference mRNA (i.e., with the highest M values) were Cyp (M value: 0.36) and TuB (M value: 0.43), respectively. In terms of candidate reference miRNAs for mature miRNA data normalization, the most stable reference miRNA in leaf and root samples were 166a/167a (M value: 0.17) and 156a/171a (M value: 0.21), respectively. The most unstable reference miRNA were 1520d (M value: 0.56) and 172a (M value: 0.46), respectively. Additionally, to determine the optimal number of candidate reference mRNA genes for data normalization in leaf and root samples, the V2/3 values were 0.064 and 0.085, respectively (Fig 3B and 3C); for candidate reference miRNAs in leaf and root samples, the V2/3 values were 0.073 and 0.084, respectively (Fig 3E and 3F). To determine the optimal number of candidate reference mRNA/miRNAs for data normalization in leaf and root samples all V2/3 values were lower than the cut-off value of 0.15 (Fig 3B, 3C, 3E and 3F), which meant a third reference mRNA/miRNA was unnecessary. Therefore, the optimal number for both the reference mRNA/miRNAs was 2. The best reference mRNA combinations for mRNA and miRNA precursor data normalization in leaf and root samples were EF1a + Act and EF1a + EF1b, respectively. The best reference miRNA combinations for mature miRNA data normalization in leaf and root samples were 166a + 167a and 171a + 156a, respectively.

thumbnail
Table 2. Ranking of candidate reference mRNAs and candidate reference miRNAs using the geNorm, NormFinder, and BestKeeper algorithms.

https://doi.org/10.1371/journal.pone.0155606.t002

NormFinder assessment of expression stability

NormFinder was an algorithm to identify the optimal reference gene for data normalization among a set of candidate reference genes [26]. It ranked the expression stability of the candidate reference genes according to their S values. The lowest S values represented the most stable expression levels. The S values and rank order of the candidate reference mRNA/miRNAs are provided in Table 2 and S3 File. Regarding candidate reference mRNA genes for mRNA and miRNA precursor data normalization, the most stably expressed reference mRNA were TuB (S value: 0.093) in leaf samples and EF1b (S value: 0.106) in root samples. The most unstable reference mRNA were Cyp (S value: 0.198) in leaf samples and TuB (S value: 0.240) in root samples. TuB was the most stable gene in leaf samples, but the least stable gene in root samples. In terms of candidate reference miRNAs for mature miRNA data normalization, the most stable reference miRNA were 166a (S value: 0.021) in leaf samples and 167a (S value: 0.067) in root samples. The most unstable reference miRNA were 397a (S value: 0.218) in leaf samples and 172a (S value: 0.225) in root samples.

NormFinder uses a solid statistical framework to estimate the overall variability in expression of the candidate genes and the variations among tested sample subgroups [26]. Using NormFinder, we calculated intra-group and inter-group variations among the five groups exposed to different abiotic stresses (Table 2 and S3 File). We determined that EF1b + Act (S value: 0.063) was the best mRNA combination in leaf samples and 60s + EF1B (S value: 0.081) in root samples for mRNA and miRNA precursor data normalizations. In leaf and root samples, the best miRNA combinations for mature miRNA data normalization were 166a + 167a (S value: 0.020) and 167a + 171a (S value: 0.067), respectively.

BestKeeper assessment of expression stability

The BestKeeper algorithm evaluates gene expression stability based on the standard deviation (SD), coefficient of variation (CV), and Pearson correlation coefficient (r). Estimating the reference gene expression stability based on the SD and CV values, reference genes with SD values higher than 1 should be excluded. The BestKeeper Index is calculated as the geometric mean of each candidate gene. This index is used along with a pair-wise correlation analysis of all pairs of candidate genes to identify the optimal reference genes[27]. A high r value indicates the reference gene pairs have very similar expression patterns, which makes them stable reference genes. The BestKeeper analysis results are provided in Table 2 and S4 File. Each candidate reference mRNA/miRNA had low variation values with SD values below 1. The r values suggested that the most suitable reference mRNAs for mRNA and miRNA precursor data normalization were Act (r value: 0.95) in leaf samples and EF1b (r value: 0.87) in root samples. We determined that the most suitable reference miRNAs for mature miRNA normalization was 166a (r value: 0.95) in leaf samples and 167a (r value: 0.95) in root samples. Cyp (r value: 0.65) and TuB (r value: 0.74) were the most unsuitable reference mRNAs for mRNA and miRNA precursor data normalization in leaf and root samples, respectively. 397a was the most unsuitable reference miRNA for mature miRNA normalization in leaf (r value: 0.31) and root (r value: 0.68) samples.

Comprehensive assessment of expression stability

The summarized data in Table 2 indicated that the geNorm, NormFinder, and BestKeeper algorithms produced similar rank orders for the candidate reference mRNA/miRNAs expression stability. However, the results lacked consistency, which may be due to differences in calculation methods among the three algorithms.

We also determined the consensus rank of each candidate reference mRNA/miRNAs using geometric means of the rankings determined by the geNorm, NormFinder, and BestKeeper algorithms (Table 2). Based on a comprehensive assessment, we identified Act and EF1b as the most suitable reference mRNA for mRNA and miRNA precursor data normalization in leaf and root samples, respectively, while the most unsuitable candidate reference mRNA were Cyp and TuB in leaf and root samples, respectively. We determined that 166a and 167a were the most suitable reference miRNA for mature miRNA data normalization in leaf and root samples, respectively, while the most unsuitable candidate reference miRNA were 397a and 172a in leaf and root samples, respectively.

Validation of reference mRNA/miRNA for relative quantification

To validate the accuracy of the selected reference mRNA/miRMA in analyses of relative quantification, we compared the relative expression levels of Pre-396a, 396a and GRF9, calculated using the best combination indicated by geNorm and NormFinder and each candidate reference mRNA/miRMA. As indicated in Fig 4, the relative quantities of Pre-396a, 396a, and GRF9 calculated using data normalized by the best combination and the most stable reference mRNA/miRNA were similar. Except for ABA treatment, all abiotic stresses produced the expected expression tendencies for each gene in leaf and root samples. For example, we observed up-regulated Pre-396a and 396a expression (Fig 4A and 4B) and down-regulated GRF9 expression (Fig 4C) in leaf samples. In contrast, the expression of Pre-396a and 396a were down-regulated (Fig 4D and 4E) while the expression of GRF9 was up-regulated in leaf samples (Fig 4F). As reported previously, miR396 expression is affected by various environmental stresses [28,4751]. Our results indicated that there was a positive correlation between Pre-396a and 396a expression levels, which is likely because 396a is derived from Pre-396a. The gene expression levels of 396a and GRF9 were inversely correlated, probably because GRF genes are miR396 targets in plants [5153]. However, following ABA-treatment, there were no obvious changes to GRF9, Pre-396a, and 396a expression levels. This was consistent with the results of previous studies that concluded miR396 gene expression was insensitive to ABA[54]. All of our results were in agreement with those of other reports, which increased our confidence in the reliability and viability of the selected reference mRNA/miRNAs. However, one way analysis of variance detected significant differences between the results generated using recommended mRNA/miRNA combinations or the most stable reference mRNA/miRNA and the results produced using the other candidate reference mRNA/miRNAs. Significant differences were observed for the results from the same experimental conditions. For example, in some instances, TuA as the reference gene for mRNA and miRNA precursor normalization in leaf tissue samples; 172a or 393a as the reference miRNA for mature miRNA normalization in leaf tissue samples; TuB as the reference gene for mRNA and miRNA precursor normalization in root tissue samples; 172a as the reference miRNA for mature miRNA normalization in root tissue samples. These results indicated that the use of inappropriate reference genes may lead to inaccurate results.

thumbnail
Fig 4. Validation of reference mRNAs and reference miRNAs for relative quantification under various stress conditions.

Comparison of Pre-396a, 396a, and GRF9 relative expression levels determined using data normalized with the best reference mRNA/miRNA combination (according to geNorm and NormFinder) and each candidate reference mRNA/miRNA (consensus ranking from best to worst). (A) Pre-396a normalized by the best reference mRNA combination and each candidate reference mRNA in leaf samples; (B) Pre-396a normalized by the best reference mRNA combination and each candidate reference mRNA in root samples; (C) Mature miR396 normalized by the best reference miRNA combination and each candidate reference miRNA in leaf samples; (D) Mature miR396 normalized by the best reference miRNA combination and each candidate reference miRNA in root samples; (E) GRF9 normalized by the best reference mRNA combination and each candidate reference mRNA in leaf samples; (F) GRF9 normalized by the best reference mRNA combination and each candidate reference mRNA in root samples. Samples were treated for 3 h with various stresses (PEG-simulated drought, salinity, alkalinity, salinity+alkalinity, and ABA). The error bars represent the standard deviations (SD) of three biological replicates. P < 0.05 (instead of *) indicated a significant difference according to one way analysis of variance between the data normalized using the recommend combinations or the best reference mRNA/miRNA and the data normalized using the other candidate ones.

https://doi.org/10.1371/journal.pone.0155606.g004

Discussion

RT-qPCR has become the preferred method for determining gene expression levels because it is more sensitive, simpler, and less time-consuming than approaches such as northern blotting, microarrays, and deep sequencing[57]. Normalization using reference genes eliminates sampling differences to enable the identification of gene-specific variations and is an essential component of a reliable RT-qPCR assay[16]. For an ideal reference gene, the gene expression variations among tested samples should be minimal or none [18, 19]. However, such genes are difficult to identify because plant gene expression is affected by environmental conditions[2024]. In accordance with the MIQE, the reference genes for RT-qPCR assays should be identified and selected using specific experimental protocols [16].

Reference genes in soybean have been evaluated after exposure to some abiotic stresses in other studies. However, the current study differs from previous investigations because one of our objectives was to identify more stable reference genes for data normalization of samples exposed to specific experimental conditions. As Kulcheski et al. (2010)[36] reported, miR156b and miR1520d are stable reference miRNA in different soybean tissues and genotypes treated with abiotic or biotic stresses. In our study, miR1520d expression was unstable in leaf and root samples. Our observations indicated that 166a in leaf samples and 167a in root samples are the most suitable reference miRNA for mature miRNA data normalization. This inconsistency in results may be due to differences in cultivars, developmental stages, tissues, abiotic stress conditions, and treatment methods causing changes to miRNA expression levels. As Le et al. (2012) [35] reported, the most stable reference genes in soybean root samples treated with four abiotic stresses [dehydration, salinity, cold, and ABA (n = 9)] were ABC, 60s, and EF1b. According to our findings, the three most stable reference genes in soybean root samples treated with various abiotic stress conditions (PEG-simulated drought, salinity, alkalinity, salinity+alkalinit, and abscisic acid) were EF1b, EF1a, and 60s. The abiotic stresses used in these two studies were similar, but we used different treatment methods. Additionally, we investigated the effects of alkalinity and salinity+alkalinity, which differed from other related studies. Similar to other reports, we identified 60s and EF1b as suitable reference genes, which inspired confidence when they normalized in root samples treated with abiotic stresses. Ma et al. (2013)[37] reported that EF1b and UKN2 were the two most reliable reference genes for soybean leaf tissue (‘Jidou 7’ and ‘Nannong 1138–2’) treated with various stresses (i.e., salt, drought, darkness, and soybean mosaic virus infection). However, it is contradictory, the results of Ma et al. (2013)[37] also showed EF1b and UKN2 expression levels are not stable enough in soybean leaf tissue exposed to drought stress. In our study, we identified Act as the most stable reference gene in soybean leaf samples exposed to PEG-simulated drought, salinity, alkalinity, salinity +alkalinity, and ABA.

In conclusion, we evaluated the expression stability of eight candidate reference mRNAs and eight candidate reference miRNAs for use in the normalization of expression level data for mature miRNAs, miRNA precursors, and mRNAs in soybean leaf and root samples treated with various abiotic stresses. The results of ΔCq analysis (Fig 1) and geNorm anlysis (Figs 2 and 3) indicated the intra-tissue candidate reference mRNA/miRNA expression was more stable than the inter-tissue expression. Especially, we found candidate reference miRNAs unlike tanditional housekeeping gene, there weren’t reference miRNA or reference miRNA combination stable enough for data normalization in leaf and root combined samples. This may be because of the differential expression of candidate reference mRNA/miRNAs between leaf and root samples. So we suggest to select reference mRNA/miRNA in leaf or root samples, not in leaf and root combined samples. In addition, it is proved that candidate reference miRNAs have some specific function[5558], this may lead to expression variation related to its own function under some conditions. The selected reference mRNA/miRNAs may express differently in different experimental conditions, so it is recommended to keep the conditions used in this study. In contrast the stability of candidate reference mRNA/miRNAs has to be re-evaluated to adopt their own experimental conditions. Generally, We assessed the expression stability of candidate reference mRNA/miRNAs in our study according to the consensus rank orders determined by geNorm, NormFinder, and BestKeeper algorithms in leaf or root samples (Table 2). We recommended that Act is the most suitable reference gene in leaf samples and EF1b in root samples for mRNA and miRNA precursor normalization. The most suitable reference is l66a in leaf samples and 167a in root samples for mature miRNA normalization. Furthermore, the use of multiple combinations of reference for data normalization can improve the reliability of RT-qPCR results, but it is time-consuming and more expensive. We also suggest to consider the accuracy and cost of reference combinations to be used. The best combination of reference genes for mRNA and miRNA precursor normalization were EF1a + Act (geNorm) and EF1b + Act (NormFinder) for leaf samples, and EF1a + EF1b (geNorm) and 60s + EF1b (NormFinder) for root samples. The best combination of reference for mature miRNA normalization were 166a + 167a (geNorm or NormFinder) in leaf samples, and 171a + 156a (geNorm) and 167a + 171a (NormFinder) for root samples.

Supporting Information

S1 Fig. Quality assessment of random RNA samples.

(A) 1% agarose gel electrophoresis was used to check RNA integrity. The ribosomal RNA bands are clearly visible and the 25S:18S ratio was approximately 2:1, which indicated that the RNA was intact. (B) RNA purity was determined using a NanoDrop spectrophotometer. A 260 nm/280 nm optical density ratio ranging from 1.8 to 2.0 indicated high quality RNA. (A and B) Samples were randomly selected from the various abiotic stress samples used for RT-qPCR (see S1 Table).

https://doi.org/10.1371/journal.pone.0155606.s001

(TIF)

S1 File. Primer pair amplification specificities for RT-qPCR.

https://doi.org/10.1371/journal.pone.0155606.s002

(DOC)

S2 File. Primer pair annealing locations on their respective transcripts.

https://doi.org/10.1371/journal.pone.0155606.s003

(DOC)

S1 Table. Description of the samples used for RT-qPCR.

https://doi.org/10.1371/journal.pone.0155606.s006

(DOC)

Author Contributions

Conceived and designed the experiments: HL WL FW. Performed the experiments: WL YD YZ. Analyzed the data: WL. Contributed reagents/materials/analysis tools: HC XL NW YD HY MZ KC. Wrote the paper: HL WL FW. Grammar and spelling: AJ.

References

  1. 1. Voinnet O.Origin, biogenesis, and activity of plant microRNAs. Cell. 2009; 136(4):669–687. pmid:19239888
  2. 2. Sunkar R, Li YF, Jagadeeswaran G.Functions of microRNAs in plant stress responses. Trends Plant Sci. 2012;17(4):196–203. pmid:22365280
  3. 3. Kumar R.Role of microRNAs in biotic and abiotic stress responses in crop plants. Appl Biochem Biotechnol. 2014;174(1):93–115. pmid:24869742
  4. 4. Zhou M, Luo H. MicroRNA-mediated gene regulation: potential applications for plant genetic engineering. Plant Mol Biol. 2013;83(1–2):59–75. pmid:23771582
  5. 5. Derveaux S V J, Hellemans J.How to do successful gene expression analysis using real-time PCR. Methods. 2010;50: 227–230. pmid:19969088
  6. 6. Heid CA, Stevens J, Livak KJ,Williams PM. Real time quantitative PCR. Genome Res. 1996;6: 986–994. pmid:8908518
  7. 7. Nolan T H R, Bustin SA. Quantification of mRNA using real-time RT-PCR. Nat Protoc. 2006;1: 1559–1582. pmid:17406449
  8. 8. Schmittgen TD, Jiang J, Liu Q, Yang L. A high-throughput method to monitor the expression of microRNA precursors. Nucleic Acids Res. 2004;32(4):e43. pmid:14985473
  9. 9. Schmittgen TD, Lee EJ, Jiang J, Sarkar A, Yang L, Elton TS, et al. Real-time PCR quantification of precursor and mature microRNA. Methods. 2008;44(1):31–38. pmid:18158130
  10. 10. Shi R, Chiang VL. Facile means for quantifying microRNA expression by real-time PCR. Biotechniques. 2005;39(4):519–525. pmid:16235564
  11. 11. Varkonyi-Gasic E, Wu R, Wood M, Walton EF, Hellens RP. Protocol: a highly sensitive RT-PCR method for detection and quantification of microRNAs. Plant Methods. 2007;3:12. pmid:17931426
  12. 12. Varkonyi-Gasic E, Hellens RP. Quantitative stem-loop RT-PCR for detection of microRNAs. Methods Mol Biol. 2011;744:145–157. pmid:21533691
  13. 13. Lao K, Xu NL, Yeung V, Chen C, Livak KJ, Straus NA. Multiplexing RT-PCR for the detection of multiple miRNA species in small samples. Biochem Biophys Res Commun. 2006;343(1):85–89. pmid:16529715
  14. 14. Chen C, Ridzon DA, Broomer AJ, Zhou Z, Lee DH, Nguyen JT,et al. Real-time quantification of microRNAs by stem-loop RT-PCR. Nucleic Acids Res. 2005;33(20):e179. pmid:16314309
  15. 15. Mou G, Wang K, Xu D, Zhou G. Evaluation of three RT-qPCR-based miRNA detection methods using seven rice miRNAs. Biosci Biotechnol Biochem. 2013;77(6):1349–1353. pmid:23748783
  16. 16. Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, et al. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009;55(4):611–622. pmid:19246619
  17. 17. Lefever S, Hellemans J, Pattyn F, Przybylski DR, Taylor C, Geurts R, et al. RDML: structured language and reporting guidelines for real-time quantitative PCR data. Nucleic Acids Res. 2009;37(7):2065–2069. pmid:19223324
  18. 18. Kozera B, Rapacz M.Reference genes in real-time PCR. J Appl Genet. 2013;54(4):391–406. pmid:24078518
  19. 19. Dundas J, Ling M. Reference genes for measuring mRNA expression. Theory Biosci. 2012;131(4):215–223. pmid:22588998
  20. 20. Machado RD, Christoff AP, Loss-Morais G, Margis-Pinheiro M, Margis R, Korbes AP. Comprehensive selection of reference genes for quantitative gene expression analysis during seed development in Brassica napus. Plant Cell Rep. 2015;34(7):1139–1149. pmid:25721200
  21. 21. Lin YL, Lai ZX. Evaluation of suitable reference genes for normalization of microRNA expression by real-time reverse transcription PCR analysis during longan somatic embryogenesis. Plant Physiol Biochem. 2013;66:20–25. pmid:23454294
  22. 22. Kozera Bartłomiej, Rapacz Marcin. Reference genes in real-time PCR. J Appl Genetics. 2013;54:391–406.
  23. 23. Dekkers BJ, Willems L, Bassel GW, van Bolderen-Veldkamp RP, Ligterink W, Hilhorst HW, et al. Identification of reference genes for RT-qPCR expression analysis in Arabidopsis and tomato seeds. Plant Cell Physiol. 2012;53(1):28–37.
  24. 24. Huang H, Cheng F, Wang R, Zhang D, Yang L. Evaluation of four endogenous reference genes and their real-time PCR assays for common wheat quantification in GMOs detection. PLoS One.2013;8(9):e75850. pmid:24098735
  25. 25. Vandesompele J, De Preter K, Pattyn F, Poppe B, Van Roy N, De Paepe A, et al. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002;3(7): 0034.
  26. 26. Andersen CL, Jensen JL, Orntoft TF. Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets. Cancer Res. 2004;64(15):5245–5250. pmid:15289330
  27. 27. Pfaffl MW, Tichopad A, Prgomet C, Neuvians TP. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper—Excel-based tool using pair-wise correlations. Biotechnol Lett.2004;26(6):509–515. pmid:15127793
  28. 28. Li H, Dong Y, Yin H, Wang N, Yang J, Liu X, et al. Characterization of the stress associated microRNAs in Glycine max by deep sequencing. BMC Plant Biol. 2011;11:170. pmid:22112171
  29. 29. Fan XD, Wang JQ, Yang N, Dong YY, Liu L, Wang FW, et al. Gene expression profiling of soybean leaves and roots under salt, saline-alkali and drought stress by high-throughput Illumina sequencing. Gene. 2013;512(2):392–402. pmid:23063936
  30. 30. Schmutz J, Cannon SB, Schlueter J, Ma J, Mitros T, Nelson W, et al. Genome sequence of the palaeopolyploid soybean. Nature. 2010;463(7278):178–183. pmid:20075913
  31. 31. Wang Y, Yu K, Poysa V, Shi C, Zhou Y.Selection of reference genes for normalization of qRT-PCR analysis of differentially expressed genes in soybean exposed to cadmium. Mol Biol Rep. 2012;39(2):1585–1594. pmid:21625860
  32. 32. Jian B, Liu B, Bi Y, Hou W, Wu C, Han T. Validation of internal control for gene expression study in soybean by quantitative real-time PCR. BMC Mol Biol.2008;9:59. pmid:18573215
  33. 33. Hu R, Fan C, Li H, Zhang Q, Fu YF. Evaluation of putative reference genes for gene expression normalization in soybean by quantitative real-time RT-PCR. BMC Mol Biol.2009;10:93. pmid:19785741
  34. 34. Li Q, Fan CM, Zhang XM, Fu YF. Validation of reference genes for real-time quantitative PCR normalization in soybean developmental and germinating seeds. Plant Cell Rep. 2012;31(10):1789–1798. pmid:22588479
  35. 35. Le DT, Aldrich DL, Valliyodan B, Watanabe Y, Ha CV, Nishiyama R, et al. Evaluation of candidate reference genes for normalization of quantitative RT-PCR in soybean tissues under various abiotic stress conditions. PLoS One. 2012;7(9):e46487. pmid:23029532
  36. 36. Kulcheski FR, Marcelino-Guimaraes FC, Nepomuceno AL, Abdelnoor RV, Margis R. The use of microRNAs as reference genes for quantitative polymerase chain reaction in soybean. Anal Biochem, 2010;406(2):185–192. pmid:20670612
  37. 37. Ma S, Niu H, Liu C, Zhang J, Hou C, Wang D. Expression stabilities of candidate reference genes for RT-qPCR under different stress conditions in soybean. PLoS One. 2013;8(10):e75271. pmid:24124481
  38. 38. Feng H, Huang X, Zhang Q, Wei G, Wang X, Kang Z. Selection of suitable inner reference genes for relative quantification expression of microRNA in wheat. Plant Physiol Biochem. 2012;51:116–122. pmid:22153247
  39. 39. Ferdous J, Li Y, Reid N, Langridge P, Shi BJ, et al. Identification of reference genes for quantitative expression analysis of microRNAs and mRNAs in barley under various stress conditions. PLoS One. 2015;10(3):e0118503. pmid:25793505
  40. 40. Goodstein DM, Shu S, Howson R, Neupane R, Hayes RD, Tricker PJ. Phytozome: a comparative platform for green plant genomics. Nucleic Acids Res. 2012;40:D1178–1186. pmid:22110026
  41. 41. Griffiths-Jones S, Grocock RJ, van Dongen S, Bateman A, Enright AJ.miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res.2006;34:D140–144. pmid:16381832
  42. 42. Plasterer TN. PRIMERSELECT. Primer and probe design. Methods Mol Biol. 1997;70:291–302. pmid:9089623
  43. 43. Thornton B, Basu C. Rapid and simple method of qPCR primer design. Methods Mol Biol. 2015;1275:173–179. pmid:25697660
  44. 44. Zuker M. Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Res. 2003;31(13):3406–3415. pmid:12824337
  45. 45. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;15(3):403–410.
  46. 46. Royeen CB. The boxplot: a screening test for research data. Am J Occup Ther. 1986; 40(8):569–571. pmid:3752225
  47. 47. Wang M, Wang Q, Zhang B. Response of miRNAs and their targets to salt and drought stresses in cotton (Gossypium hirsutum L.). Gene. 2013;530(1):26–32. pmid:23948080
  48. 48. Kantar M, Lucas SJ, Budak H. miRNA expression patterns of Triticum dicoccoides in response to shock drought stress. Planta. 2011;233(3):471–484. pmid:21069383
  49. 49. Omidbakhshfard MA, Proost S, Fujikura U, Mueller-Roeber B. Growth-Regulating Factors (GRFs): A Small Transcription Factor Family with Important Functions in Plant Biology. Mol Plant. 2015;8(7):998–1010. pmid:25620770
  50. 50. Gao P, Bai X, Yang L, Lv D, Li Y, Cai H, et al. Over-expression of osa-MIR396c decreases salt and alkali stress tolerance. Planta. 2010;231(5):991–1001. pmid:20135324
  51. 51. Liu D, Song Y, Chen Z, Yu D. Ectopic expression of miR396 suppresses GRF target gene expression and alters leaf growth in Arabidopsis. Physiol Plant. 2009;136(2):223–236. pmid:19453503
  52. 52. Wang L, Gu X, Xu D, Wang W, Wang H, Zeng M, et al. miR396-targeted AtGRF transcription factors are required for coordination of cell division and differentiation during leaf development in Arabidopsis. J Exp Bot. 2011; 62(2):761–773. pmid:21036927
  53. 53. Casadevall R, Rodriguez RE, Debernardi JM, Palatnik JF, Casati P. Repression of growth regulating factors by the microRNA396 inhibits cell proliferation by UV-B radiation in Arabidopsis leaves. Plant Cell. 2013;25(9):3570–3583. pmid:24076976
  54. 54. Xing L, Zhang D, Zhao C, Li Y, Ma J, An N, et al. Shoot bending promotes flower bud formation by miRNA-mediated regulation in apple (Malus domestica Borkh.). Plant Biotechnol J. 2016 Feb;14(2):749–70. pmid:26133232
  55. 55. Kitazumi A, Kawahara Y, Onda TS, De Koeyer D, de los Reyes BG. Implications of miR166 and miR159 induction to the basal response mechanisms of an andigena potato (Solanum tuberosum subsp. andigena) to salinity stress, predicted from network models in Arabidopsis. Genome. 2015, 58(1):13–24. pmid:25955479
  56. 56. Liu N, Wu S, Van Houten J, Wang Y, Ding B, Fei Z, et al. Down-regulation of AUXIN RESPONSE FACTORS 6 and 8 by microRNA 167 leads to floral development defects and female sterility in tomato. J Exp Bot. 2014, 65(9):2507–2520. pmid:24723401
  57. 57. Morea EG, da Silva EM, ES GF, Valente GT, Barrera Rojas CH, Vincentz M, et al. Functional and evolutionary analyses of the miR156 and miR529 families in land plants. BMC Plant Biol. 2016, 16(1):40.
  58. 58. Ma Z, Hu X, Cai W, Huang W, Zhou X, Luo Q, et al. Arabidopsis miR171-targeted scarecrow-like proteins bind to GT cis-elements and mediate gibberellin-regulated chlorophyll biosynthesis under light conditions. PLoS Genet. 2014, 10(8):e1004519. pmid:25101599