Abiotic and biotic stresses impose challenging physiological hurdles to plants. As a response to adverse environmental conditions, plants re-program their cellular activities through multiple gene regulatory mechanisms including post-transcriptional regulation of gene expression. Transcription factors (TFs) and non-coding RNAs are the two major regulatory elements in functional genomics (Deihimi et al., 2012; Mahdi et al., 2013; Panahi et al., 2013; Chiasson et al., 2014).
Small RNAs, particularly microRNAs (miRNAs), have emerged as key post-transcriptional gene regulatory molecules in plants and animals (Trindade et al., 2010; Giacomelli et al., 2012; Sunkar et al., 2012). MiRNAs are 18–22 nucleotides in length and regulate expression of their target genes by degradation or translational repression (Voinnet, 2009; Sunkar et al., 2012). In particular, miRNAs play important roles in plant responses to abiotic and biotic stresses including cold, salt, heat, dehydration, oxidative and mechanical stresses (Jin et al., 2013; Panahi et al., 2013). For example, some of the genes in ARF transcription factor family are targeted by miR160 and miR167 through auxin regulation interference (Gutierrez et al., 2009), and seem to be important for the attenuation of plant growth and development under stress conditions (Sunkar et al., 2012). Another microRNA, miR398, is associated with response to abiotic and biotic stress condition (Zhu et al., 2011) and its expression is up-regulated in response to water deficit in Medicago truncatula (Trindade et al., 2010). MiR403 controls the expression of AGO2 (Allen et al., 2005), which is known to have an antiviral role (Harvey et al., 2011). Srivastava et al. (2013) reported that miR426 is up-regulated in response to arsenic stress in Brassica juncea. In another study, miR842 expression levels were observed to be up-regulated in Arabidopsis (Moldovan et al., 2010). These are handful examples of miRNAs that are differentially expressed in response to various stress conditions.
Sunflower (Helianthus annuus L.) is one of the most important oilseed crops and is resistant to various abiotic stresses, due to its metabolic, physiological and morphological adaptation strategies. This crop is of special interest for its adaptation to high temperatures, limited water availability, high salinity and heavy metal concentrations in the soil (Merah et al., 2012). Using bioinformatics methods, previous studies have identified and reported several miRNAs in sunflower. A number of these were experimentally verified (Barozai et al., 2012; Monavar Feshani et al., 2012). However, their expression patterns under different stress conditions have not been studied.
We have examined the roles of seven of these sunflower miRNAs in response to different abiotic stresses including drought, high temperature, salinity and cadmium. Seven mentioned miRNAs, miR160, miR167, miR172, miR398, miR403, miR426 and miR842, play crucial roles in plant growth and show different expression patterns under both biotic and abiotic stress conditions in several plant species. Simultaneously, we analyzed the expression levels of these miRNAs concurrently with their potential targets to corroborate the putative miRNA-target relationship. Based on our results, several gene regulatory networks were constructed to gain a comprehensive overview of the role of these miRNAs in stress responses.
Materials and Methods
Prediction of miRNA Targets
Most plant miRNAs bind to their targets with perfect or near-perfect complementarity (Jones-Rhoades and Bartel, 2004). We have used this characteristic for predicting potential miRNA targets. Sequences of pre-miRNAs were obtained from the public database, miRBase (Kozomara and Griffiths-Jones, 2011). The number of mismatches at complementary sites between miRNA sequence and potential mRNA target were kept below 4 and no gap at the complementary sites was tolerated. The miRNA targets were predicted using psRNAtarget online sever (http://plantgrn.noble.org/psRNATarget/).
Plant Material, Growth Conditions and Stress Treatments
Sunflower (Helianthus annuus, variety Sirna) seeds were surface sterilized and placed in petri dishes containing two layers of damp sterile filter paper. For germination, the seeds were incubated at room temperature for 4 days.
For applying drought stress treatment, five germinated seeds were sown in round pots (10.5 × 13.5 cm, D × H; five plants/pot) filled with soil (¾ river sand and ¼ soil). Plants were grown in the glasshouse under the controlled condition (16 h of light, 24°C; 8 h of dark, 15°C) and well watered.
For heat, salinity and cadmium stress treatments, five seedlings were grown in a hydroponic system. The system contained aerated Hoagland solutions, pH 5.8 (Hoagland and Arnon, 1950) in plastic pots (1.5 liters volume). One-week-old plantlets were planted in ¼ strength Hoagland's nutrient for 1 more week to immerse the growing roots in the nutrient solution. Plantlets were then transferred to half-strength nutrient solution and grown until 6-leave stage (about 4 weeks). Four weeks old plantlets were then submitted to different stresses as specified in the following sections.
Four weeks old plantlets (6 leaf stage) of sunflower were subjected to water stress by withholding water for 12, 24 and 48 h. Relative water content (RWC) was measured in leaves to determine the plant water status following Faatsky method (Čatský, 1960).
Four weeks old plantlets in ½ strength Hoagland solution were transferred to humidity growth chamber (Memmert, Germany) with 70% relative humidity and maintained at 42 ± 1°C for 1.5, 3 and 6 h (Senthil-Kumar et al., 2003; Mashkina et al., 2010; Mangelsen et al., 2011). We selected stress-related HSP70-related protein gene (Gene Bank ID: AAB57695.1) to investigate the effect of heat stress in root tissue of sunflower.
Plantlets were moved into solution contained either 75 or 150 mM of NaCl (Di Caterina et al., 2007). Leaf and root samples from stressed plants were collected 24 h later alongside control samples. Potassium to sodium ratio was measured as criteria for examining different salinity stress levels.
Plantlets were transferred into containers with either 5 or 20 mg/L of CdCl2 (Niu et al., 2007). Stress and control samples were collected after 24 h from leaf and root tissue and immediately frozen in liquid nitrogen before they stored at −80°C. Cadmium content was measured as an indicative of stress.
Total RNA Isolation from Leaves and Roots
Total RNAs were extracted from sunflower root and leaf tissues using TRIZOL Reagent (Invitrogen, USA) according to the manufacturer's instructions. The extractions were performed separately for two independent biological replicates. Equal amount of total RNA was subsequently pooled for each sample based on their concentration. The quantity and quality of purified RNAs were assessed with a Biophotometer spectrophotometer (Eppendorf, Germany) and the integrity was evaluated by conducting a denaturing agarose gel electrophoresis. All RNA samples were stored at −80°C until further processing.
Stem-loop Reverse Transcription
Stem-loop RT primers for H. annuus miR160, miR167, miR172, miR398, miR403, miR426 and miR842 were designed according to Varkonyi-Gasic et al. (2007) (Table S1). The miRNA specific RT reactions were carried out using prime script RT reagent kit (Cat. #RR037A, TAKARA, Japan). The miRNA stem–loop reverse transcription was accomplished according to Varkonyi-Gasic et al. (2007) with minor modifications at cDNA synthesis, using 200 ng of total RNA sample of treated leaf and root tissues (1 μL), 1 μL of 1 μM miRNA primer RT and 5.5 μL nuclease-free water. The mix was incubated at 65°C for 5 min followed by incubation on ice for 2 min. Two microliters of Primescript buffer (5×) and 0.5 μL Primescript RT enzyme mix were added to each tube and the RT reaction was performed for 30 min at 16°C followed by 60 cycles at 30°C for 30 s, 42°C for 30 s, 50°C for 1 s and terminated by incubation at 85°C for 5 min. Control reactions were performed with no RT primer and no RNA samples.
Quantitative Real-time PCR
To measure and compare the expression levels of the selected H. annuus miRNAs in root and leaf tissues under different stress treatments, RT-qPCR was conducted using SYBR Premix Ex Taq II (Cat. #RR820A, TAKARA) in a Rotor-GeneQ Real-Time PCR system (Qiagen).
For RT-qPCR analysis, 6 μL SYBR Premix Ex Taq II (2×), 1 μL (10 pmol) each of forward and reverse primers, 3 μL nuclease-free water and 1 μL RT stem-loop cDNA product were mixed. Forward primers were specifically designed for each miRNA, and 5′-GTGCAGGGTCCGAGGT-3′ sequence was used as the universal reverse primer (Varkonyi-Gasic et al., 2007) (Table S1). The RT-qPCR reactions were performed using following conditions; initial denaturation at 95°C for 30 s, followed by 40 cycles at 95°C for 5 s and 60°C for 20 s. The melting curves were generated during denaturation step at 95°C followed by the cooling of PCR products at 50°C and the fluorescence signals were collected continuously from 50 to 95°C as the temperature increased at 0.2°C per second. All reactions were repeated four times. For each condition, the RT-qPCR experiments were run as pooled biological duplicates and expression levels were normalized according to the previous studies (Schmittgen and Livak, 2008). Relative fold changes in miRNA expression were calculated using the comparative Ct (2−ΔΔCt) method with 18S rRNA as the endogenous control (Schmittgen and Livak, 2008). The miRNAs clustering analysis, based on their relative expression, and Pearson correlation between expression miRNAs and target genes in each stress were performed using Minitab 16.0. The heat map of expression was visualized using R 3.0 package.
Analysis of miRNAs Target Expression with RT-qPCR
To determine the expression of predicted miRNA targets and possible discovery of novel miRNA target genes under the drought, heat, salinity and cadmium stresses in sunflower, expression levels of miRNA-related target genes were measured with quantitative real-time PCR. ESTs for target quantification analysis were selected based on two criteria: (1) protein found in BLAST search should be a previously published target of the related miRNA in other plant species. (2) The EST should have a possible protein-coding ORF that can allow us to design RT-qPCR primers for the conserved 3′ UTR region. RT-qPCR primers for these selected genes (Table S1) were designed using Primer3 according to following criteria: (1) Primers 3′ self-complementary = 0. (2) Primer annealing temperatures = 62 ± 0–3°C. (3) Product size limit for primer pairs = 150 bp.
Total cDNAs were synthesized from 1 μg RNA using RevertAid™ First Strand cDNA Synthesis Kit (Fermentas) according to the manufacturer's instructions. Further expression analysis was performed for all miRNA targets using the same batch of RNA samples for miRNA RT-qPCR assay. One microliter of this cDNA was amplified with 0.6 μM of specific primers in a total of 10 μl volume using SYBER Premix Ex Taq II (Cat. #RR820A, TAKARA) with Rotor-GeneQ Real-Time PCR system (Qiagen). The quantification was performed using 18S rRNA (Gene ID: 18250984, forward: TTCAGACTGTGAAACTGCGAATGG /reverse: TCATCGCAGCAACGGGCAAA) as a normalizer and four independent PCR results with acceptable efficiency (1.8–2.2) were averaged. Specified RT-qPCR thermal setup was adjusted as follows: pre-denaturation step at 95°C for 1 min, followed by 40 cycles of 95°C for 30 s, 60°C for 1 min. The melting curves were generated as mentioned for miRNAs.
Network Analysis of miRNAs
To construct miRNAs and target genes interaction network, RESNET Plant database of Pathway Studio software v.9 (Elsevier) was used. This database includes new aliases for genes in the model and non-model plant species including barely, corn, tomato, potato and tobacco and collects data through MedScan (text mining tool) to extract functional relationships between miRNAs, proteins, stresses and cellular processes (Nikitin et al., 2003; Alimohammadi et al., 2013; Ebrahimie et al., 2014). In addition to adding the result of this study to prediction database of RESNET database, this database was also updated by MedScan, especially from literature on miRNAs/target genes and drought, heat, salt and cadmium stress conditions before network construction. To predict the interaction networks, the software makes different groups of miRNAs and finds the relations between a protein and its group using algorithms such as Fisher's Exact Test (Alanazi et al., 2013; Ebrahimie et al., 2014).
Network constructed by union selected and physical interaction algorithms were used to make statistical subnetworks based on miRNAs and their target genes. Also, GO (gene ontology) analysis was performed through DAVID Functional Annotation web-tool (https://david.ncifcrf.gov/) and separate tables were produced for biological process, molecular function and cellular component categories. MiRNAs and putative target genes are marked with yellow circles (Figures 6, 7).
In this study, expression levels of seven miRNAs including miR160, miR167, miR172, miR398, miR403, miR426 and miR842 and their targets were investigated in response to four abiotic stress conditions including drought, heat, salt and cadmium stresses. In addition, the impacts of these stresses on the physiological and molecular characteristics of plants were evaluated by measuring the changes in RWC, expression of HSP70-related protein, sodium, potassium and cadmium concentrations in different tissues of plant, grown under stress (Supplementary 1). The expression level of HSP70-related protein up-regulated over time compared to control in root tissue with the highest level at 6 h (Table S2). Moreover, alteration of sodium, potassium and cadmium concentration were remarkably higher in root tissue compared to leaf, in particular in later stage of stress (Table S3A). The average RWC of plants in response to drought stress showed slight decline at 12 and 24 h, and it started to reduce sharply subsequently (Table S3B). The results indicated that plants were significantly affected by stress exposure and were able to initiate stress related responses.
Prediction and Determination of H. annuus miRNA Targets
The mRNA targets of these miRNAs (Table 1) were predicted via psRNAtarget database using selected miRNAs as queries. Targets of miR167 and miR403 were previously reported in the other plant species (Table 1). We predicted and analyzed new targets for miR160, miR398, miR426 and miR842 for the first time in sunflower. Protein sequence analysis of the new putative Helianthus miR160 target, QHG18J04.yg.ab1 showed the presence of ARF, phosphorylase and kinase domains (data not shown) in this protein. Interestingly, COX5b, which was previously introduced as miR398 target in Arabidopsis, was identified as a predicted target for miR172 in sunflower.
Table 1. Putative Helianthus annuus stress-inducible miRNA target genes.
Expression of miRNAs and Putative Targets under Drought Stress
MiRNAs expression levels were significantly down-regulated (P < 0.05) in leaves of plants grown under drought stress with the lowest expression levels were at 48 h (except miR403). MiR172 expression modulation was not significant (P < 0.05) in all period of stress in root tissue. The expression patterns of miR160, miR426 and miR842 were similar in both tissues, except miR160 which showed opposite pattern at 48 h after drought stress in root tissue. The miRNAs, miR167 and miR398 showed similar trend under drought stress in both tissues. They were down-regulated at all-time points within range of 2- to 19-fold change, whereas expression of miR167 was slightly up-regulated at 48 h compared to moderate stress in root tissue. MiR403 showed the opposite pattern in leaf and root tissues at 24 h. Its expression abruptly decreased (24-fold compared to control) in leaf tissue. It, however, exhibited an increasing trend in root tissue while its expression was still lower than control condition (Table 2; Figure 1).
Table 2. Expression of miRNAs and their target genes response to abiotic stress in Helianthus annuus.
Figure 1. Expression patterns of miRNAs and their target genes under drought stress. Sunflower seedling at the six leaf stages grew on soil (¾ river sand and ¼ soil) under normal conditions and leaves and roots of them were harvested as control then they treated by withholding the water for 12, 24 and 48 h. (A) Expression patterns of miRNAs in the leaf. (B) Expression patterns of miRNAs in the root. (C) Expression patterns of target gene in the leaf. (D) Expression patterns of target genes in the root.
The expression of QHG18J04.yg.ab1 gene as target for the miR160, was down- regulated at 12 h and 24 h (2- to 8-fold), and up-regulated at 48 h in both tissues. Except at 48 h in leaf tissue, correlation of miR160 and its target showed similar expression pattern over time and tissue. The expression of ARF6 (target of miR167), COX5b (target of miR172), AGO2 (target of miR403) and TPP2 (target of miR426) displayed similar patterns in both tissues; their expression was induced only in the root tissue at 48 h of drought stress. Interestingly, the coherent trend was observed only at 48 h in the root tissue for miR167, miR403 and miR426 and their respective targets. The transcript of NtGT5b (target of miR398) was constantly decreased in both tissues with 2- to 10-fold under drought conditions where the minimum peak was at 48 h and 24 h after treatment in root and leaf tissue, respectively. The expression of (R)-mandelonitrile lyase (target of miR842) increased at 12 h and decreased at subsequent time points in leaf tissue whereas it showed opposite pattern in root. In general, most targets exhibited their highest expression level at 48 h after stress depending on the tissue (Table 2; Figure 1).
Expression of miRNAs and Putative Targets under Heat Stress
The expression levels of all seven miRNAs in leaf tissue were slightly up-regulated after 1.5 h exposure to heat stress, except for miR172 which had constant expression level. However, expression of miRNAs in leaves indicated a mixed pattern at 3 and 6 h after treatment. The expression levels of miR167 and miR172 were down-regulated within a range of 2- to 3-fold at 3 and 6 h after stress. But in root tissues, their reduction was between 2 and 12 fold at these two time points. MiR398 and miR403 showed similar trend in leaf and root tissues. Their expression was immediately up-regulated in leaf at the initial time point, although, their expressions were decreased at 3 and 6 h after stress. The levels of their expression were higher than control condition. However, in root tissue, they were up-regulated at 1.5 h but down-regulated subsequently compared to control with a 2- to 8-fold change at 3 and 6 h. The miR160, miR426 and miR842 exhibited similar trend under heat stress in leaf tissue. Their expressions were induced at 1.5 h and down-regulated at 3 h and again induced at severe stress. Nevertheless, their expressions were lower compared to the control. In root tissue, except 1.5 h after stress, miR842 and miR426 showed similar pattern within a range of 4- to 38-fold change after stress. The expression of miR160 instantly decreased at 3 h (10-fold change compared to control) and its decrease continued constantly at 6 h (Table 2; Figure 2).
Figure 2. Expression patterns of miRNAs and their target genes under heat stress. Plants treated at 42 ± 1°C for 1.5, 3, and 6 h. (A) Expression patterns of miRNAs in the leaf. (B) Expression patterns of miRNAs in the root. (C) Expression patterns of target gene in the leaf. (D) Expression patterns of target genes in the root.
The expression of QHG18J04.yg.ab1 was significantly down-regulated at 1.5 h and was reversed at 3 h, while its expression was only up-regulated at 6 h in root. The pattern of its transcript had negative correlation with miR160. Transcript of ARF6 was down-regulated at all-time points with a sharp decrease in the root tissue at 3 h. The temporal variation of COX5b was significant only in the root tissue while it decreased sharply at 3 h after heat stress. AGO2 expressed at a significant level at 1.5 h in both tissues where its expression was up-regulated with more than two-fold change, compared to the control. NtGT5b displayed similar pattern in both tissues. The expression level declined at an initial time point and was induced at 3 h after treatment, whereas it showed coherent type with miR398. In leaves, the expression of (R)-mandelonitrile lyase was abruptly down-regulated at 6 h, whereas its expression was declined slightly over time in the root. Expression of TPP2 dropped in leaf tissue, but it decreased in root only at 3 h with 14-fold changes after stress. However, there was a slight increase in its expression at 6 h in the root. The expression was still lower compared with control (Table 2; Figure 2).
Expression of miRNAs and Putative Targets under Salt Stress
The expression of miR167 showed opposite pattern in both tissues. In leaves, salt stress reduced the expression of miR167 by two-fold while the expression level was up-regulated in the root tissues with a two-fold change. Interestingly, miR403 exhibited opposite pattern of expression in leaf and root, whereas its expression had increasing and declining trend in leaf and root, respectively. The expression trends of miR160, miR426 and miR842 were similar in both tissues at severe stress. Their expression showed decreasing pattern in leaf and increasing pattern in root tissue. Interestingly, an abrupt gradient was observed for miR842 in leaf tissue at 150 mM NaCl. None of the tissues disclosed a significant alteration for miR172 expression in response to salt stress. Salt stress induced the miR398 expression in both tissues with the sharp increase observed after 75 mM NaCl treatment in both tissues. In spite of this, there was a considerable reduction in its expression at 150 mM concentration in leaf and root, and the expression level was higher than control condition (Table 2; Figure 3).
Figure 3. Expression patterns of miRNAs and their target genes under salt stress. Plant subjected under NaCl treatment in two concentration, 75 and 150 mM. (A) Expression patterns of miRNAs in the leaf (B) Expression patterns of miRNAs in the root. (C) Expression patterns of target gene in the leaf. (D) Expression patterns of target genes in the root.
Under salt stress, the transcripts of QHG18J04.yg.ab1 slightly increased in the root. Although, there was a decline pattern at 150 mM of NaCl in comparison with an early stage of stress in the leaf tissue, its expression was higher compared to control. The expression of ARF6 was not significant in the root. However, its expression showed down-regulation at both stages of treatment in leaf tissue. Its expression level revealed a slightly increasing trend at 150 mM compared with 75 mM of NaCl. COX5b expression was lower in both tissues at all stages of treatment with sharp reduction at 150 mM of NaCl in the leaf tissue. NtGT5b exhibited similar expression pattern in both tissues where its mRNA level was up-regulated at 75 mM concentration and was abruptly decreased at 150 mM. The AGO2 transcript was highly accumulated at all stages in the root tissue, but in the leaf tissue, was induced at 75 mM of NaCl, and down-regulated subsequently compared with control. The expression pattern of TPP2 was opposite in leaf and root tissue, as decreased in leaf and increased in root tissue with its peak of expression at severe stress. The expression of (R)-mandelonitrile lyase was induced in both tissues with its peak at 150 mM concentration (Table 2; Figure 3).
Expression of miRNAs and Putative Targets under Cadmium Stress
Cadmium treatment resulted in up-regulation of all miRNAs in roots, with fold-changes between 4 and 385. The lowest and highest peaks were for miR842 in 20 mg/L and miR398 in 5 mg/L concentration of cadmium, respectively. The expression of miR426 was induced only at 5 mg/L in leaf, while in root induced at all stages of stress. The opposite pattern for miR160 and miR842 observed only at 20 mg/L in leaf tissue. The level of miR160 was drastically decreased in leaf tissue, but the expression of miR842 was slightly raised. In root tissue, their expression was increased although; the levels of miR842 at 20 mg/L were declined in comparison to 5 mg/L, where its expression was still higher compared to control. MiR403 expression was increased in both tissues whereas expression levels in root tissue were higher than leaf tissue. The highest expression of miR403 was in 5 mg/L concentration with 52-fold change compared to control. The expression pattern of miR172 was induced at all stages of treatment in both tissues with its peak at 20 mg/L concentration in root. Interestingly, the level of miR398 was up-regulated at 5 mg/L, but down-regulated at 20 mg/L with 1.6-fold change in leaf tissue. Similar trend was observed for miR167 and miR398 in root tissue, with the highest peak at 5 mg/L concentration (Table 2; Figure 4).
Figure 4. Expression patterns of miRNAs and their target genes under cadmium stress. Plants treated with cadmium in 5 and 20 mg/L concentration. (A) Expression patterns of miRNAs in the leaf. (B) Expression patterns of miRNAs in the root. (C) Expression patterns of target gene in the leaf. (D) Expression patterns of target genes in the root.
Transcript of QHG18J04.yg.ab1 showed up-regulation trend in both tissues where in root, it was over 10-fold higher than leaf. The target of miR167 and miR403 exhibited similar trend in both tissues as their expressions were significantly induced in the root with their highest peak at 20 mg/L of CdCl2. The expression of COX5b and TPP2 was declined in leaf and was increased in root tissue in both concentrations of cadmium with its peak in 20 mg/L. (R)-mandelonitrile lyase highly accumulated in the root with abrupt increase at 20 mg/L. In contrast, its expression was decreased in leaf tissue after treatment. Expression of NtGT5b revealed temporal variation where it was down-regulated at 5 mg/L and up-regulated at 20 mg/L in both tissues. Interestingly, miR172, miR398, miR426 and miR842 showed inverse correlation with their targets in the leaf tissue at 5 mg/L (Table 2; Figure 4).
Cluster Analysis of miRNAs in H. annuus Root and Leaf Tissues based on their Expression Patterns
Cluster analysis was performed to further analyze the pattern of expression of these miRNAs under different stress conditions. A comparison of expression profiles of these seven miRNAs in both tissues was conducted to find their tissue-specific expression patterns, regardless of their response to different stress conditions. Three clusters were formed, when the expression pattern of the seven miRNAs were analyzed in response to stress (Figure 5). Cluster I included miR167 and miR398. The expression pattern of leaf and root was opposite, but their trend was similar in each tissue over time and stress. In root tissue, these miRNAs revealed slight alteration in three treatments; the highest peak was in 5 mg/L cadmium stress. Cluster II contained miR172 and miR403 families and showed a fluctuated pattern in leaves and roots in all of stress conditions. Median value of expression of these miRNAs exhibited increasing pattern in four stresses. Due to similarity in expression patterns of miR160, miR426 and miR842, they grouped in the cluster III. They displayed similar trend in each tissue during stress. The lowest peak of miR842 was in severe stress through drought and salt condition in the leaf tissue, while it was at 3 h after heat treatment in the root tissue.
Figure 5. Clustering of miRNAs expression profiles. Heat map diagram of miRNA expression prepared with two-way unsupervised hierarchical clustering of miRNAs expression under different abiotic stress. miRNAs are given in the rows and each columns represent a sample. The miRNA clustering tree is shown on the left (cluster I, II and III). Abbreviations: L, Leaf; R, Root; h, Hour; D, Drought stress; S, Salt Stress; Cd, Cadmium stress.
Relationship of miRNAs and Target Gene in Stress Response
Analysis of Pearson correlation showed positive and negative correlation between expression of miRNAs and target genes against stress in both tissues. The results indicated significant negative correlation between miR172, miR398 and miR403 and their putative targets in leaf tissue under cadmium stress. On the contrary, miR398 presented significant negative correlation in root tissue after heat treatment. Weak correlation was observed between some miRNAs and their candidate targets in both tissues against specific stress such as miR160 in drought stress and miR167 after heat treatment (Table S4). Some miRNAs displayed significantly positive correlation, which indicated other mechanisms are involved in target gene regulation.
Interaction between miRNAs and TFs-Mediated Gene Regulatory Subnetworks
The regulatory subnetworks that are constructed for each miRNA elucidated some of the intermingled miRNA and TF relationships as well as miRNA-miRNA relationships and the involvement of other miRNA families in miRNA specific post-transcriptional regulation pathways (Figures 6, 7). References of interaction relations between miRNAs and genes in subnetworks and their correlated references are listed in Table S5.
Figure 6. miRNA-mediated gene regulatory sub networks in response to abiotic stress. miRNA-target-abiotic stress interaction are shown. (A) Feedback loop between miR160, its target and abiotic stress. (B) Feedback loop between miR167, its target and abiotic stress and abiotic stress. (C) Feedback loop between miR172, its target and abiotic stress. (D) Feedback loop between miR398, its target and abiotic stress.
Figure 7. miRNA-mediated gene regulatory sub networks in response to abiotic stress. miRNA-target-abiotic stress interaction are shown. (A) Feedback loop between miR403, its target and abiotic stress. (B) Feedback loop between miR842, its target and abiotic stress.
GO analysis is a robust approach in understanding underlying molecular mechanisms of different cellular events and offer a reliable tool for GO based gene selection (Fruzangohar et al., 2013). GO classification showed that the target genes are involved in auxin signaling, RNA mediated silencing, ethylene signaling, DNA methylation and response to abiotic stress.
Subnetwork of miR426 did not display any connectivity with other participants in the network. As miR426 exhibits stress and tissue specific pattern, this miRNA is a suitable candidate for further studies in the context of miRNA mediated gene regulatory pathway. MiR842 showed connectivity with MBB18_8 and (R)-mandelonitrile lyase that drought, heat, salt and cadmium affected its expressions (Figure 7B).
MiR160 centered network exposed interaction between QHG18J04.yg.ab1 and three TFs belonging to the ARF family (ARF16, ARF17, and ARF18), where four abiotic stress affect expression of miR160 and QHG18J04.yg.ab1 (Figure 6A). Interestingly, in this network, two other ARF family members; ARF6 and ARF8 regulate JA (Jasmonic acid) pathway and are in turn regulated by miR167. In addition, sub network of miR167 includes regulation of TAS1C, which encodes a ta-siRNA (Figure 6B). Expression of TFs involved in ethylene signaling pathway and flower and seed development is regulated by the members of the miR172 family. But miR172 itself is regulated by miR156 and DDL (a gene silencer) (Figure 6C). MiR398 revealed cross talk with SPL, with central roles in many cellular processes such as auxin metabolism, root growth and cytokinesis; where SPL is regulated by miR156, miR157 and miR159. Therefore, these three miRNAs are speculated to regulate the expression of miR398 indirectly. MiR156 have excelled as a key regulator in the regulation of other miRNAs. MiR398 causes mRNA cleavage and disease resistance, which regulates expression of CSDs gene and NtGT5b (Figure 6D). Subnetwork of miR403 revealed that SE and miR402 regulate AGO2. It is known that SE has role in miRNA biogenesis, RNA splicing and chromatin modification. MiR402 displayed connectivity with DML1 and DML3, which are involved in DNA methylation and transcriptional control (Figure 7A).
The role of seven miRNAs in response to several abiotic stresses was studied in H. annuus leaf and root tissues by RT-qPCR. Changes in RWC, expression of HSP70-related protein, sodium, potassium, and cadmium concentration revealed obviously various mechanisms stimulation to cope with stress. The measured alterations of metabolic and physiological status in H. annuus provide additional support for modification reactions toward stress conditions.
The temporal and spatial expression profiles of seven miRNAs in both tissues at altered time points were in agreement with previous studies (Sunkar et al., 2012; Zandkarimi et al., 2015; Zhang, 2015). For instance, high-throughput sequencing of Brassica napus showed up-regulation of miR172 in the root tissue under cadmium excess, up-regulation of miR167 from 2 to 24 h after exposure in 300 mM NaCl in Arabidopsis and down-regulation of miR398 under drought stress in cotton which treated with PEG.
The different expression patterns of miR167 observed in photosynthetic and non-photosynthetic tissues suggest that miR167 may have a role in tissue-specific adaptation to stress. MiR160 and miR167 regulate ARF families and play roles in the auxin signaling pathway (Khraiwesh et al., 2012), and were differentially expressed under stresses compared to control conditions. In addition, ARF domain was found in the protein sequence of QHG18J04.yg.ab1 (target of miR160). Therefore, this gene may be involved in auxin signaling pathway and have a role in promoting plant tolerance to stress. The alteration of miR160 and miR167 in H. annuus at early stage of treatments suggests that they are responsive to early stress response. Reduced expression of miR160 and miR167 in leaves under severe stress may alter basal level of auxin and consequently restrict plant growth through the antagonism between abscisic acid and auxin because of escape of tension. In turn, this could lead to attenuation of plant growth and development under stress to cope with the imposed stress (Sunkar et al., 2012). Even though they have common targets, in our study, miR160 and miR167 did not group together in one cluster. On the other hand, these miRNAs showed various expression patterns in other plant species under different stress (Khraiwesh et al., 2012; Sunkar et al., 2012). Besides, under normal conditions, their expression pattern was also variable in different plant genus and species (Zeng et al., 2010).
Some studies have suggested roles for miR172 during cadmium, drought, cold and heat stress conditions (Sunkar et al., 2012; Zhou et al., 2012). In line with the previous studies, miR172 revealed temporal up- or down-regulation in leaf and root in response to stress conditions except for the salt stress. Expression of miR172 showed inverse correlation with its target COX5b, only in leaf tissue under cadmium stress. The post-transcriptional silencing of COX5b by miR172 may reflect the providence of Helianthus plant from energy loss via avoiding respiration under excessive cadmium concentrations (Sunkar et al., 2012).
We found that miR398 is especially up-regulated in leaf tissues during heat stress in line with the previous reports (Guan et al., 2013). Interestingly, miR398 was down-regulated in root tissue while expression of HSP70-related was up-regulated, which may indicate that miR398 has a tissue specific mode of action and localization during heat stress condition. Differential expression of miR398 in response to drought, salt, heat and cadmium stress have been shown in several species such as wild emmer wheat, Medicago, Nicotiana, Brassica and Arabidopsis (Trindade et al., 2010; Frazier et al., 2011; Kantar et al., 2011; Zhou et al., 2012; Guan et al., 2013). In contrast, under severe drought and cadmium stress conditions, miR398 was down-regulated in leaf tissue indicating a tissue-specific and stress-specific response orchestrated by this miRNA. In this study, the probable target of miR398 was NtGT5b which is a microsomal enzyme responsible for glucuronidation reactions with a role in the storage of secondary metabolites and plants defense against stress (Miners et al., 2002; Ko et al., 2006). On the other hand, new target genes were predicted for miR398 and miR167 in Phaseolus vulgaris and Malus hupehensis in vegetative phase which involved in monogalactosyl diacylglycerol synthase, acyltransferase and dioxygenase, gluconeogenesis pathway and glycolytic process (Heyndrickx and Vandepoele, 2012; Han et al., 2014; Xing et al., 2014). Furthermore, in Nicotiana tabacum in response to TiO2 nanoparticles, these miRNAs were grouped in one cluster (Frazier et al., 2014). All of these results have led to the conclusion that miR398 may be involved in the sugar biosynthesis pathway, associated with reduction in energy consumption for photosynthesis, and increase the tolerance of plant in abiotic stress conditions.
Surprisingly, the pattern of miR403 was varying in each stress condition. Its expression was declined and was induced in both tissues in drought and cadmium stress, respectively. In contrast, miR403 displayed increased expression in the leaf and decreased expression in the root during heat and salt stress. In the previous study, abundance of miR403 was high in heat and salt libraries in Raphanus sativus (Wang et al., 2014; Sun et al., 2015). As a result, this discrepancy suggested that miR403 was potentially expressed in stress-, tissue-, and species- specific manner during abiotic stress. Also, the expression of AGO2 at all stages of treatment in both tissues was aberrant. In plants, AGOs are involved in various small RNA pathways from post-transcriptional gene silencing to epigenetic silencing phenomena such as RNA-directed DNA methylation (RdDM) pathway in Arabidopsis (Schraivogel and Meister, 2014). AGO1 and AGO2 proteins were regulated by miR403 (Figure 7A). In addition, AGO2 has been shown to have an antiviral role (Harvey et al., 2011). Therefore, it is possible that this miRNA is involved in the regulation of miRNA-mediated RNA cleavage carried out by other miRNAs during stress conditions. Furthermore, AGO1 is also regulated by miR172 family (Ronemus et al., 2006), which may designate a crucial role for miR172 in general miRNA mediated gene silencing pathways (Figure 6C). This result, based on altered expression of miR403 and its target under different abiotic stress and its subnetwork (Figure 7A) which showed DML1 and DML3, involved in DNA methylation, may pave the ways for intricate control mechanisms for drought, heat, salt and cadmium stress tolerance in sunflower. Interestingly, miR172 and miR403 were grouped in one cluster (Figure 5) and they showed common targets; AGO1 and AGO2, which may suggest a general role for these miRNAs in small RNA pathway and DNA methylation.
In this study, miR426 and miR842 show differential expression under abiotic stress. However, their expression exhibited stress-dependent manner during stress which was declined in root tissue under heat and drought stress and was reversed in salt and cadmium treatment. Also, inverse correlation of these miRNAs with their target was temporary in some stages. Their possible targets, TPP2 and (R)-mandelonitrile lyase, were induced in response to some abiotic stresses and are reported to be involved in defense mechanism, oxidation-reduction process, cyanide biosynthetic process and alcohol metabolic process, which is indicated as a biocatalyst in organic chemistry (Hu and Poulton, 1999; Shima et al., 2007). MiR842 revealed no significant change after waterlogging conditions (Moldovan et al., 2010), whereas it was repressed after ABA treatment of roots in Arabidopsis (Jia and Rock, 2013) which was similar with the expression of
Spring frost is an important environmental stress that threatens the production of Prunus trees. However, little information is available regarding molecular response of these plants to the frost stress. Using high throughput sequencing, this study was conducted to identify differentially expressed miRNAs, both the conserved and the non-conserved ones, in the reproductive tissues of almond tolerant H genotype under cold stress. Analysis of 50 to 58 million raw reads led to identification of 174 unique conserved and 59 novel microRNAs (miRNAs). Differential expression pattern analysis showed that 50 miRNA families were expressed differentially in one or both of almond reproductive tissues (anther and ovary). Out of these 50 miRNA families, 12 and 15 displayed up-regulation and down-regulation, respectively. The distribution of conserved miRNA families indicated that miR482f harbor the highest number of members. Confirmation of miRNAs expression patterns by quantitative real- time PCR (qPCR) was performed in cold tolerant (H genotype) alongside a sensitive variety (Sh12 genotype). Our analysis revealed differential expression for 9 miRNAs in anther and 3 miRNAs in ovary between these two varieties. Target prediction of miRNAs followed by differential expression analysis resulted in identification of 83 target genes, mostly transcription factors. This study comprehensively catalogued expressed miRNAs under different temperatures in two reproductive tissues (anther and ovary). Results of current study and the previous RNA-seq study, which was conducted in the same tissues by our group, provide a unique opportunity to understand the molecular basis of responses of almond to cold stress. The results can also enhance the possibility for gene manipulation to develop cold tolerant plants.
Almond is a perennial species belonging to Prunoideae, a subfamily of Rosaceae. Almond is predominantly self-incompatible. It has sixteen (2n = 16) small chromosomes and a small diploid genome of approximately 300 Mbp . Spring frost is one of the most important environmental factors, which limits almond’s production. During winter, buds are dormant and hard. But they swell and expand into blossoms when the temperature rises. As a result, they become less resistant to the cold. Frost can damage 10% to 90% of the almond tree in different phonological stages of flowering . Following our previous study on transcriptome of almond reproductive tissues in response to spring frost , we have decided to further expand our understanding of the molecular basis of almond responses to cold stress by conducting small RNA-Seq analysis in the same tissues. As miRNAs are mega bio-regulators in plants  and are involved in different stress responses [5, 6], we thought that changes in gene expression in response to cold would also reflect a similar pattern in miRNAs expression pattern.
Plants involvement of miRNAs and their target genes in response to stress conditions is well documented in several studies such as cold signaling regulation [7, 8] and salt and nutrient deficiency [9–11]. For instance, Barakat et al.  have identified 157 conserved and 230 non-conserved miRNAs sequences in Prunus persica in response to cold stress. In Medicago truncatula, 283 and 293 miRNAs were identified from control and drought stress libraries, respectively . Zhang et al.  reported 106 conserved miRNAs from tea leaves treated with cold (Camellia sinensis). Sun and colleagues  identified 136 known and 68 novel miRNAs in Raphanus sativus in a salinity stress experiment. Among these miRNAs, several were characterized as stress responsive. Many studies have shown that miR156/157, miR159/319, miR165/166, miR169, miR172, miR393, miR394, miR396, miR397 and miR398 are up-regulated in response to cold stress [8, 16, 17]. Some of the studies, moreover, have shown the down-regulation of miRNAs such as miR156g-j, miR475a, b and miR476a under cold stress in Populus . Although information regarding miRNAs expression in many species is building up quickly, we do not have any reports regarding the number of miRNAs and their expression pattern in plants such as almond (Prunus. dulcis (Mill.) D.A.Webb).
Different approaches including forward genetics, direct cloning and bioinformatics prediction followed by experimental validation can identify miRNAs [6, 19]. However, these approaches have limitations such as time, cost and the versatility of the methods. In recent years, next generation sequencing platforms like Genome Analyzer (Illumina Inc.) or Genome Sequencer TM FLX (454 Life Science TM and Roche Applied Science) have been used to detect small RNA molecules and determine the expression levels of both conserved and novel miRNAs [19–22]. In many plant species, this method has been successfully applied to discover miRNAs responses to different abiotic stress. Some of the examples of such applications are the identification of drought stress responsive miRNAs in Populus trichocarpa , salt stress responsive miRNAs in Saccharum spontaneum  and heavy metal stress responsive miRNAs in Medicago truncatula .
Until now, there is no report on the responses of almond’s miRNAs to cold stress, especially in the reproductive tissues, where the cold stress would have the most devastating effects. As the whole genome sequencing of almond is not available, the results of our two studies could be merged to uncover the molecular responses of reproductive tissues of almond to cold stress.
Construction and high-throughput sequencing of small RNAs libraries in almond
cDNA libraries of small RNAs were constructed from anther and ovary under cold stress and control conditions in order to identify miRNAs in almond. Our sample collection consisted of anther samples under stress (HSA) and under control (HCA); and ovary samples under stress (HSO) and under control (HCO). All samples were obtained from H gentoptype of almond, which is beleived to be tolerant to cold stress. We have obtained 50 to 58 million raw reads, using Illumina sequencing platform with average read length of 51 nucleotides and the quality of 39 for all samples. GC contents for all samples were 52%, while it was 51% in the HSA.
HCA samples yielded the highest number of small RNA reads of about 58,664,982, while HSO samples contained 50,245,161 reads. 85.74% (12.77% unique reads) from HCA, 78.33% (13.45% unique reads) from HSA, 91.03% (16.09% unique reads) from HCO and 88.85% (17.77% unique reads) from HSO were obtained after filtering. The adaptors and all contaminanted small RNAs including tRNA, snRNA, snoRNA and rRNA were removed and residual reads were queried against known miRNAs in miRBase (version 20). 2,137,645, 1,626,754, 2,767,683 and 2,197,689 reads from HCA, HSA, HCO and HSO libraries were identified as conserved miRNAs. The rest of the sequences were considered as non-conserved sRNA sequences.
Identification of conserved miRNAs in almond
After performing a blast against miRBase database (version 20), sequence similarity releaved the presence of 131, 122, 123 and 119 miRNAs in the HCA, HCO, HSA, and HSO libraries, where 94 miRNAs were unique. Out of these 94 uniquly identified conserved miRNA families in four librarires, 26 were highly conserved and 68 were conserved in some plant species (S1 Table). The highly conserved miRNA families contained 1–10 members in our dataset. miR482f family had the highest number of members (10 members). This was followed by miRNA 171 family, having 8 members (Fig 1). On the other hand, among the conserved miRNAs family, miR6267, miR6291, miR7122 with 3 members had the highest number of members, while the most of the other conserved miRNA families conatined only one member (Fig 2).
The distribution of highly conserved miRNAs family size in prunus dulcis.
The distribution of conserved miRNAs family in prunus dulcis.
A comparison between miRNA family members in these four libraries showed that 32 miRNA families had tissues-specific expression pattern. Among these miRNAs, 21 were expressed only in the anther and 11 in the ovary (Fig 3A). Interstingley, 10 miRNAs found to be expressed solely under cold-stress conditions. Nineteen miRNAs were specifically expressed under control conditions and therfeore, their expression was not detectable under stress conditions (Fig 3B). Also the size distribution of miRNA families between the four libraries revealed 62 common miRNAs, though the number of members were different in each library (Fig 4A and 4B).
Common and specific miRNA families identified by high throughput sequencing.
The distribution of known miRNA families in HCA, HSA, HCO and HSO libraries.
Identification of novel miRNAs in almond
In our study, we found many reads without any hits in the RNAs datasets. They were 289,302, 245,497, 151,928 and 177,393 reads in the HCA, HSA, HCO and HSO libraries, respectively. Therefore, these reads were mined to identify non-conserved miRNAs. Following miRNA annotation criteria set by Meyer and collaborators , we successfully identified 59 novel miRNAs sequences. The stem-loop structures of these miRNA precursors were obtained by Mfold  and are presented in Fig 5. The mauture length of these miRNAs varied between 18 and 22 nt. The most abundant length was 21 nt. The precursor lengths were found to be between 45 to 147 nt. The minimum free energy (MFE) varied from -71.8 to -8.5 kcal/mol (S2 Table). The MFE index (MFEI) for these precursers, as the most important criterion for prediction of stem-loop structure of miRNA , were calculated according to Zhang et al. . Notabley, the mature sequences of 34 novel miRNAs were located on 3' arm of their precursers and the rest were located on 5' arm.
The structure of candidate novel miRNAs precursor.
Expression pattern analysis of known miRNAs
Expression analysis of conserved miRNAs indicated that miR159, miR166, miR482 and miR1511 families had the highest expression in the four libraries. However, distinct expression variation was observed within the members of some families. For instance, in HCA library miR482f had 109668 copies, whereas the read number for miR482b-5p was only 10. The expression pattern of miRNAs in anther indicated that 53 miRNAs had higher expression in control conditions compared to ovary. In contrast, 48 miRNAs were up-regulated under stress conditions in anther compared to ovary. Moreover, 41 miRNAs showed similar pattern in both tissues exposed to cold stress. Therefore, we suggest that these miRNAs should be considered as cold stress responsive miRNAs. From these miRNAs, 22 were down-regulated and 19 were up- regulated. Among differntially expressed miRNAs, miR482d-3p was up-regulated in anther, while it was down-regulated in ovary. In contrast, miR172a-5p and miR1511-3p were up- regulated in ovary and down-regulated in anther (S3 Table). Comparative analysis of expression pattern in these two reproductive tissues revealed 11 and 9 up-regulated miRNAs in anther and ovary, respectively (Fig 6A). Additionally, 20 and 13 down-regulated miRNAs were found in these tissues respectively (Fig 6B). Futhermore, the analysis identified 5 up-regulated and 11 down-regulated miRNAs in anther with no significant changes in ovary. However, we found 8 repressed miRNAs in anther, while their expression was iduced in ovary tissues.
Relative differentially expressed conserved miRNA in reproductive tissues.
Expression pattern validation of miRNAs in reproductive tissues of almond
To validate the expression patterns detected for miRNAs in the high throughput sequencing, 16 differentially expressed miRNAs were selected and analyzed with qPCR. In addition, to further investigate the miRNA responses to cold stress, we included a rather cold-sensitive genotype (Sh12) in our analysis. Expression pattern comparison between small RNA sequencing data and qPCR results indicated very similar patterns of expression for 11 and 7 miRNAs in anther and ovary samples under -2°C treatment, respectively. In anther, furthermore, miR162, miR166d, miR168, miR171a, miR398a-3p, miR403, miR482f, miR6285 and miR8123-5p were expressed differentially in H genotype and Sh12 varieties. However, in ovary, miR403, miR1511-3p and miR7122a-3p displayed differential expression in these two genotypes. On the other hand, miR7122a-3p in anther and miR162, miR166d, miR398b, miR398a-3p, miR160a, miR319a, miR168, miR6285 and miR8123-5p in ovary had similar expression pattern in both tolerant and sensitive varieties. Unsurprisingly, some miRNAs displayed different mode of expression in anther and ovary samples of these two almond verities. For example, miR162, miR477a-3p and miR482f in H genotype and miR168, miR403, miR6285 and miR8123-5p in Sh12 showed different expression between these studied tissues. (Figs 7, 8A and 8B). Among all 16 differentially expressed miRNAs studied by qPCR, cold stress repressed the expression of 10 miRNAs in anther tissue of H genotype and had no effects on the expression of miR160a, miR319a, miR394b, miR398b, miR1511-3p and miR7122a-5p. In contrast, miR162, miR171a, miR319a, miR394b, miR398b were up-regulated in Sh12 under cold stress. In addition, it was found that -2°C cold treatment induced the expression of miR398b, miR477a-3p, miR394b, miR7122a-3p, miR162 and miR166d in both ovary and anther tissues of Sh12 variety. Different expression pattern were detected for some miRNAs including miR1511-3p, miR398a-3p, miR166d and miR7122a-5p in these varieties treated with two temperatures. For instance the expression of miR1511-3p, miR398a-3p and miR166d was lower in 0°C compared to that of -2°C. Transcript of miR394b and miR398b was not altered in the anther tissue of H genotype in -2°C. But the two miRNAs were down-regulated in 0°C. Among the tested miRNAs, miR477a-3p, which increased 8.6- fold in the ovary tissue of H genotype under 0°C, and miR166d, with 8.27- fold increase in Sh12 under -2°C, showed the most changes in their expression pattern. In anther tissue of H genotype, miR166d and miR160a had the greatest reduction of expression under 0 and -2°C, respectively. MiR319a was the only miRNA, which was up-regulated in both varieties under both stress conditions. In contrast, miR160a was the only miRNA that was down-regulated under cold stress in both varieties. In other studies this miRNA was introduced as a cold-induced miRNA .
qPCR analysis for validation of differentially expressed miRNAs expression in reproductive tissues of almond.
qPCR analysis of cold stress-responsive miRNAs and their targets expression profiles.
In this study, among the predicted targets for cold-responsive miRNAs (S4 Table), the expression pattern of 6 miRNA's targets namely DUF1639 (miR398b target), ABC transporter and AAE7 (miR477a-3p targets), FBX (miR394b target), BHLG041 (miR7122a-3p target), DCL1 (miR162 target) and ATHB-15 and ABA-insensitive 5 (miR166d targets) were examined in the lab. Negative correlation between some of these miRNAs and related targets was detected (Table 1). For instance, cold-stress repressed miR477a-3p expression was negatively correlated with the conversely up-regulated expression of both ABC transporter and AAE7 genes in anther tissue of H genotype. On the other hand, the up-regulation of miR394b was followed by the down-regulation of the expression of FBX gene in Sh12. Among the studied targets, ABC transporter gene, one of the miR477a-3p targets, showed the highest number of fold changes in anther tissues of cold-tolerant genotype H.
Expression profile between miRNA and corresponding targets.
Identification of miRNA-guided cleavage of target mRNA using RLM-RACE
In order to corroborate whether FBX and DCL1 genes are direct targets of repression by miRNA394b and miRNA162, through miRNA-directed cleavage of their mRNA, we used the modified 5′ RACE technique for mapping the cleavage sites of miRNA on these two genes transcripts. We used RNA from plants treated for 3 h at 0°C. Cleavage of the miR394b target (FBX), and miR162 target (DCL1) was confirmed by 5′ RACE assay. PCR RACE products were purified and sequenced. The results indicated that the site of cleavage is located downstream of the second nucleotide from the 5′ end of the FBX mRNA (Fig 9A) and the second nucleotide from the 3′ end of the DCL1 mRNA (Fig 9B). These findings confirmed that FBX gene and DCL1 gene were directly targeted for cleavage by miR394b and miR162, respectively.
Experimental validation of predicted miRNA targets. Bottom strands depict target mRNA sequences.
Like other abiotic stresses, cold stress threatens many plant species. Most organisms have an optimum temperature for appropriate growth and development. Thermal fluctuations beyond the threshold level lead to stress . Plants respond to cold stress through specific mechanisms that involve reprogramming gene expression. The important role of miRNAs in response to different environmental stresses has been verified in previous studies [21, 15, 31–35]. Several works have introduced cold stress responsive miRNAs in different plant species such as Arabidopsis [7, 16], tea , tomato [36, 37], rice , wheat [39, 40], peach , Populus  and Brachypodium . No study has been done on the miRNAs regulatory response to cold stress in almond. This study, therefore, was conducted to identify the cold responsive miRNAs and their targets in almond.
Identification of conserved and novel miRNAs in almond
High throughput sequencing analysis of sRNA libraries in almond's reproductive tissues under two thermal treatments (10 and -2°C) led to identification of 94 conserved miRNA families and 59 novel miRNAs. Size distribution of conserved miRNAs showed that miR482, miR171, miR159 and miR396 families have the highest number of miRNAs. According to, the study of Barakat et al. , however, miR395 and miR169 were represented most frequently in peach. The comparison between known miRNA families in some Rosaceae species revealed that 28 miRNAs were specific to Prunus tree species. This result can potentially show the roles of these miRNAs in initiation of specific characterization in this species.
The distribution of miRNA family size in some Rosaceae species
The number of memebers in miRNA families of four Rosaceae tree species, including Prunus persica, Malus domestica, Prunus mumeand, and Pyrus bretschneideri, were compared with those of Prunus dulcis. Almond and peach shared 70 known conserved miRNA families, indicating a remarkable homology. Approximately, all highly conserved miRNAs are present in these species, except miR 535 and miR530 that are only presnt in apple and miR 397 and miR 530 that are specific to Japanese apricot. Among the miRNA families stduied, miR156 (with 31 members) was the largest family in Malus domestica. The size distribution of other groups of known miRNAs showed that some of them are shared only among almond, peach and apricot. This group of shared miRNAs includes miR6257, miR6258, miR6259, miR6261, miR6262, miR6263, miR6264, miR6266, miR6267, miR6269, miR6270, miR6271, miR6274, miR6279, miR6282, miR6283, miR6284, miR6285, miR6287, miR6288, miR6289, miR6290, miR6291, miR6292, miR6293, miR6294, miR6295 and miR6297. That is why Gao et al.  considered these miRNAs as unique to Prunus plants. Interstingley, 10 miRNA families including miR414, miR862, miR2655, miR3434, miR4995, miR5368, miR5539, miR6173, miR6207 and miR6234 did not have any orthologous in other species. As a result, they were considered to be specific to almond (Fig 10).
Comparison of known miRNAs across almond (pdu) and other tree species of Rosaceae family.
Identification of cold-responsive miRNAs in almond
The miRNAs that were identified with differential expression in the current study have different roles in plant molecular biology. Various studies have shown that miR156, miR163, miR169, miR172, miR398 and miR399 play important role in flowering time. They are also introduced as temperature-responsive miRNAs . MiRNAs, including miR172 and miR156, affect the development of the plant by controlling AP2 (floral patterning genes) and SPL (squamosal promoter binding like). These miRNAs have regulatory roles on transition from juvenile to the reproductive phase . MiR156-SPL and miR172- SVP (short vegetative phase) make up the regulatory circuit and are known as flowering-time regulators in response to ambient temperature. Using these regulatory mechanisms, plants can adapt their development with temperature fluctuation . MiR159 targets the gibberelic acid MYB (GAMBY) gene and regulates seed germination and the formation of the anther . Some of the identified miRNAs have regulatory roles on other miRNA biosynthesis. For example, miR162 effects the regulation of DCL1 (dicer-like 1) , and miR171, which in turn controls SCR gene (Scarecrow-like protein). In its own turn, SCR protein, in cooperation with SHR, activates the expression of miR165a and miR166b . Notably, some of the miRNAs are mainly involved in regulating stress-responsive genes. For instance, it was revealed that miR160 is a cold-induced miRNA . Or it was found that miR169, the main regulator of NFYA1 (Nuclear transcription factor), displayed different expression pattern in response to drought, salinity, low temperature, etc. In addition, miR7122a-3p gets involved in response to environmental tensions of plants by targeting BHLG041 gene as part of replication elements . In study of Alisoltani et al , this miRNA was reported as a positive regulator under cold stress. In this investigation, miR7122a-3p/HOS1 was introduced the only post-transcriptional regulator in gene regulatory network in almond in response to cold stress. MiR166 regulates several genes such as ATHB, REV and ABA insensitive 5, is involved in stress-response of plants to different conditions. In soybean, for example, miR166b was up-regulated under drought, salinity, and alkalinity , while miR166a-5P and miR166f were down-regulated under abiotic stress . This variable regulatory pattern revealed different function of miRNAs from the same family .
By using high throughput sequencing, this study identified 12 up-regulated miRNAs (miR159, miR164a/b/c, miR171a/c/g/f, miR394a/b, miR482f and miR7122a-5p) and 15 down-regulated miRNAs (miR156a/b, miR162, miR166a/b/c/d/e, miR396b, miR398a-3p, miR6262 and miR8123-5p) that respond to cold stress in both anther and ovary tissues of almond.
The induction of miR169, miR393 and miR394 and miR398 in almond under cold stress is in agreement with findings in Arabidopsis [7, 8, 16]. The repression of miR482 in the anther tissues of H genotype was consistent with the findings about tomato . On the other hand, the induction of miR166 and miR396 in Arabidopsis and the induction of miR477 in Populus are in contrast to the results about almond [7, 41]. The differences in the expression of these miRNAs among different studies may be due to differences in species, cultivars (genotypes), and the extent and duration of cold stress in different studies .
Among the miRNAs, which are involved in cold stress response, miR319a is known as positive regulator by high throughput sequencing and experimental confirmation. Up-regulation of miR319a under both cold stress treatments in the two studied varieties was observed. Other studies have also noted the over-expression of miR319 under cold stress treatments [29, 51, 52]. The members of this miRNA family are usually up-regulated when plants are exposed to different stresses . MiR319a is the first plant miRNA that was identified by forward genetics. This miRNA regulates TCP family, especially TCP4. This transcription factor is involved in the development of leaf and flower, mostly the petal and stamen [29, 53–55]. The other predicted target (GAMYB) regulates the growth of the plant. Zhou et al.  has demonstrated that the enhanced abiotic stress tolerance in Agrostis stolonifera is the result of the overexpression of osa-miR319a and consequently the significant down-regulation of its target gene. Therefore, previous studies and our own results imply that miR319a could be the main target for engineering tolerant plants to abiotic stress.
It is well known that the genotypes of different plant species affect their response to the same abiotic stress. To test this point, the expression pattern of 16 candidate miRNAs was investigated in a sensitive and a tolerant variety under cold stress treatments. Among the studied miRNAs, miR162, miR166d, miR168, miR171a, miR398a-3p, miR403, miR482f, miR6285 and miR8123a-5p in anther and miR403, miR1511-3p and miR7122a-3p in ovary appeared with differential expression in H genotype and Sh12. With this result on the expression pattern of miRNAs and their targets, positive and negative regulators under cold stress could be identified in almond. Such results could then be further integrated into more sophisticated genetic engineering approaches to develop tolerant plants. Among the miRNAs, miR6285 and miR8123a-3p have more potential to be considered as specific cold-responsive miRNAs in almond. Interestingly, these two miRNAs were also identified as specific miRNAs in Prunus species.
Assessment of cold stress response mediated by regulatory roles of studied miRNAs
When plants encounter different kinds of stresses, responsive genes are expressed. These stress-responsive genes could be categorized into two groups. The first group includes stress tolerance inducer proteins such as sugar transporter, which is predicted for miR394b in this study. The second group mostly contains transcription factors, which regulate signal transduction pathways. Function analysis of the second group can provide useful information on the regulatory gene networks involved in cold stress response . In this investigation the majority of the predicted targets for the altered miRNAs belong to the second group. For example it was found that miR172a-5p and miR169a regulates two transcription factors, namely APETALA2 protein and nuclear transcription factor Y (subunit A-1).
Furthermore, the expression patterns of target genes of 6 miRNAs were studied. One of these genes, ATHB (homeodomain leucine zipper) was predicted to be miR166d target. In a study by Kim et al.  on hot pepper, the induction of ATHB and CaCBFIB (as dehydration-responsive element) under low temperature stress (4°C) and their interaction during cold stress have been discovered. Kim et al.  consequently identified ATHB gene as an important regulator of dehydration-responsive element. Different studies have shown the overexpression of dehydration-responsive elements followed by the induction of cold-responsive genes, which leads to the inhibition of cellular dehydration and increase in freezing-stress tolerance . These results indicate that exposure of plants to cold stress down-regulates miR166 and, subsequently, increase the expression of ATHB gene. In the almond the negative regulatory role of miR166d on ATHB, which led to overexpression of this gene, was only observed in the anther of the H genotype under 0°C, while in other treatments such negative correlation was not detected.
The other cold-induced miRNA in almond was miR398b; and DUF1639 was predicted to be its target. In addition to DUF1639, miR398 controls the expression of at least four target genes. Many of DUFs (Domains of unknown function) are highly conserved in plants, indicating their important biological roles. For example Bischoff et al.  showed TBL role (a DUF protein) in synthesis and deposition of secondary wall cellulose. Moreover, Mawlong et al.  suggested the important roles of DUF sequence, which is present in AP2/ERF-N22 (2) (Apetela type 2 transcription factor/ ethylene responsive factor) in stress tolerance, plant growth and its development. It is believed that DUFs or PUFs (proteins of unknown function) are required in special conditions  such as cold stress. Different studies displayed the importance of DUFs in plants’ response to stress. For instance one DUF named Esk1 has been recognized with negative regulatory function in cold acclimation . In another study, mutation in ESKIMO1 (a DUF protein) induced freezing tolerance in Arabidopsis .
In this study the induction of miR394b in cold stress (0°C and -2°C) caused repression of FBX gene expression in both the reproductive tissues of Sh12 variety. In other words, the negative regulatory role of miR394b on the production of F-box protein was only observed in Sh12. The overexpression of miR394b in response to cold stress in almond was similar to those of Arabidopsis, tobacco, rice and poplar [16; 63, 64]. But the over-expression of an F-box protein gene in rice caused a reduction in abiotic stress tolerance and the growth of the root. From these results one can confer that the mechanism of plants’ response to abiotic stresses is through the induction of miR394b and the repression of F-box protein. Further, the role of miR394 in direct response to multiple stresses has been proved . In a study by Huang and collaborates , miR394 has been shown as a versatile miRNA that play important role in various stresses. A large gene family in plants encodes F-box proteins. In Arabidopsis, 700 different F-box proteins were recognized, but the function of many of them has not been investigated . These proteins are involved in a large number of biological processes such as self-incompatibility, leaf senescence, branching and response to different biotic and abiotic stress . The proteins are also involved in abiotic stress response through ubiquitin pathway . One of the reasons a negative correlation has not been observed between the expression of miR394b and its target in the H genotype might be the lag in the response to the stress. Other reasons could be that the expression of other targets such as STP13, SBT might be affected by miR394b and not FBX.
In our study, ABC transporter C family member 8-like and acetate/butyrate CoA ligase AAE7 were identified as targets of miR477a-3p. The up-regulation of these genes in the cold tolerant genotype (H) was compared to that of the cold sensitive almond varieties (Sh12) in both tissues. The ABCC multidrug resistance associated proteins (ABCC-MRP), a subclass of ABC transporters, is involved in multiple physiological processes such as cellular homeostasis and metal detoxification . The genes belonging to this subfamily were identified in various species . For instance they have 120 members in Arabidopsis thaliana and rice (Oryza sativa) . Most ABC proteins, as membrane proteins, mediate MgATP-energized trans-membrane transport and/or regulate other transporters. These proteins are recognized as the key manager for plants under abiotic and biotic stresses . The identified roles of this group of proteins include polar auxin transport, lipid catabolism, xenobiotic detoxification through transporting toxic compounds from cytosol into the vacuole  disease resistance, and stomatal function . The induction of this ABC transporter in cold-tolerant tomato cultivar has been reported . In our study, the significant up-regulation of ABC transporter in the H genotype, and not in the Sh12 genotype, revealed the involvement of this gene in cold tolerance, possibly through different processes such as detoxification mechanisms.
Acetate/butyrateCoA ligase AAEs, the other miR477a-3p target, participates in a variety of anabolic and catabolic pathways. Most AAEs share conserved structural elements such as a motif, which is necessary for binding to ATP and adenylate formation . AAE7 show high expression level in many tissues including root, leaf, stem, flowers and the developing seeds . Allen and colleagues  found the relationship between glutamine level and AAE7 accumulation. In their study, mutations in AAE7 or ACN1 (acetate non-utilizing 1) gene led to nearly 50% decrease in 13C-labelling of glutamine, which is a major carbon sink in seedlings. Glutamine (Gln) is a major nitrogen source in plants and is a central intermediate for carbon-nitrogen assembly. Glutamine synthetases (GS) was employed to maintain a sufficient Gln supply. This enzyme has cytosolic GS1 and plastidic GS2 for Gln production. Ji  showed that GS1 is responsive to various environmental stresses and also the involvement of GS1s in Gln homeostasis. For instance the activity of GS1 in tomato was increased after exposure to different NaCl concentrations . Two-dimensional (2D) gel electrophoresis analysis of proteins revealed that GS1 accumulates in response to cold treatment . According to these findings and the induction of AAE7 under cold stress in the current study, it can be concluded that miR477a-3p, via its regulatory role on AAE7, could prevents the reduction of Gln supply, which is essential for growth and developments.
Over the last decade, miRNAs have emerged as playing a major role in plants’ responses to different stress conditions. This is especially true by employing advanced NGS technologies, which have significantly increased the number of identified miRNAs in several plant species. In this study, we have outlined the current status of miRNA research in almond as an economically important fruit tree. Although miRNA research in fruit trees has mostly been descriptive to date, miRNAs offer great potential for the trait improvement of these crops. Hundreds of diverse miRNA genes including highly conserved, non-conserved and several species-specific miRNAs have been identified in peach, apple, citrus and grapevine species. Putative targets for the identified miRNAs have been described in the literature and include homologues of well-characterized targets from model plant species. An important goal of this study was to identify target mRNAs for species-specific miRNAs, and to elucidate their associated functions. Although several experimental approaches have been used to demonstrate function of miRNAs in model species such as Arabidopsis, rice and poplar, these methods have not been readily applied to fruit trees and almond to date.
One of the areas of plant biology that could benefit from NGS study is the investigation of gene expression under stress. Small RNAs contents of plant tissues could be easily sequenced and analyzed to identify important miRNAs under stress. Since spring frost stress is the most important environmental constraints for almond, we have analyzed cold-stress-related miRNAs in almond at two levels (0°C and -2°C) of cold exposure by stem_loop qPCR. Furthermore, we retrieved the expression profile of predicted target genes from two almond genotypes, which were obtained from a public database. Different expression profiles of target genes were observed among varieties. Comparative analysis of the gene expression of the miRNAs and their related targets, between cold tolerant and sensitive varieties grown under the same condition, allowed us to identify responsive miRNAs that are associated with cold tolerance. It further allowed us to make distinction between different varieties of the same species. In this study miR162, miR166d, miR168, miR171a, miR398a-3p, miR403, miR482f, miR6285, miR8123-5p, miR403, miR1511-3p and miR7122a-3p were introduced as differentially expressed miRNAs between H and Sh12 genotype. Further studies such as overexpression and/or RNAi strategies are needed to elucidate the precise role of these genes in almond.
Materials and Methods
Plant materials and stress condition
Prunus dulcis Mill genotype H, a cold tolerant and late blooming plant (REF), was used for sRNA sequencing. Sh12, which is also a cold sensitive genotype, was used for expression pattern confirmation. Branches were collected in the popcorn stage during early spring. The cut end of branches were placed in a 5% sucrose solution and then subjected to 0°C for 3h, -1°C for 2 h and -2 for 1 h, consecutively. The temperatures were chosen based on previous studies of Miranda et al.  and Kodad et al. . As control treatment, few branches were kept at 10°C. Collected samples (ovary and anther) were immediately placed in liquid nitrogen and later stored at -80°C.
RNA isolation and sequencing
Total RNA was isolated from 200 mg of anther and ovary samples using the method described by Rubio-Piña and Zapata-Pérez . In this protocol, after samples were powdered, 1 ml extraction buffer containing 2% (w/v) CTAB; 0.1 M Tris-HCl (pH 8); 1.4 M NaCl; 20 mM EDTA (pH 8); 2% (w/v) PVP and 70 μl of β-mercaptoethanol was added to each sample. Then they were incubated at 65°C for 10 min. After incubation 800 μl chloroform was added and centrifuged at 10,000 rpm for 10 min at 4°C. Two other centrifugations were done after 800 μl phenol/chloroform (1:1) and chloroform/isoamyl alcohol (24:1) were added in equal volume to the supernatant. Finally, the last step was to add LiCl (8M) and to do centrifugation. Formed RNA pellets were washed and after a short time for drying, dissolved with DEPC water. Extracted RNA were qualified and quantified after DNAaseI treatment by spectrophotometer and agarose gel electrophoresis. cDNA library construction and sequencing (Illumina HiSeq 2000 platform) were perfomed by Macrogen company (Korea). Data was produced in FASTQ format and was obtained from the company’s server. The raw sRNA data has been deposited in the sequence reads archive (SRA), NCBI, and could be accessed using Bioproject number PRJNA245549 and Biosample accession numbers SAMN04317490, SAMN04317478, SAMN04317221 and SAMN04317155.
Identification of conserved miRNAs in almond
For identification and analysis of our reads we have followed the procedure undertaken by Barakat et al.  with some modifications as follows. After quality control using SolexaQA software , sequences between 18 and 24 were extracted using an in-house perl script. Using BLASTN program, sRNA sequences that passed the size filter were then queried against rfam (http://rfam.sanger.ac.uk/), chloroplast and mitochondrial genomes (http://gobase.bcm.umontreal.ca/), rRNA (http://lowelab.ucsc.edu/GtRNAdb/), tRNA (http://www.psb.ugent.be/rRNA/), snoRNA (http://bioinf.scri.sari.ac.uk/cgi-bin/plant_snorna/home/) and all contaminating rRNA, tRNA and snoRNA sequences were removed. To identify conserved miRNAs, unique sRNA sequences were blasted against known miRNAs (http://www.miRbase.org/index.shtml) using default parameters, allowing up to two mismatches.
Identification of non-conserved, almond-specific miRNAs
Sequences with no similarity to known miRNA sequences were used to search for non-conserved miRNAs in almond. The genome of almond has not been sequencesd yet. In order to obtain the possible precurseors for miRNAs, therefore, we had to map small RNA sequences according to the almond RNA-seq assemblies provided by Mousavi et al. . Bowtie software (version 2) was used to map the reads on the contigs . Regions harboring sRNA hits (150 at each side) were then obtained using a shell script. Next, Vienna RNA package was used for predicting miRNAs secondary structure . The structures were then checked for miRNA features using MiRCheck program. Sequences that passed MiRCheck were sorted for redundancy. The sRNAs that passed the MiRCheck filter were then manually folded using mFold and assessed for miRNAs structure.
Differential expression analysis of miRNAs related to cold stress
The normalization of known miRNAs frequency in four libraries was calculated to earn transcripts per million (TPM) according to the formula:
If the abundance of normalized values was zero, it was modified to 0.01. The p-value and fold-change (log2 ratio (treat/control)) were calculated for differential expression analysis. The p-value was obtained according to Gao et al.  and Li et al. . The criteria for significant up/down-regulation was that, log2 ratio (treat/control) > 0.5 and < −0.5, and the p-value of < 0.05.
Target prediction for miRNA sequences
The target genes of miRNAs were identified using psRNATarget (http://plantgrn.noble.org/psRNATarget/). This server is widely used for plant miRNA target prediction . In order to identify almond specific target genes, which are under control by conserved miRNAs, we have used almond contigs assembled by Mousavi et al. . Appropriated BLASTN hits with the highest identity were selected as predicted targets.
Confirmation of miRNAs and their potential target genes by qPCR
Sixteen known miRNAs with the distinct difference in expression in cold stress vs. control condition were selected for further confirmation. In addition to miRNAs, eight predicted target genes of six differentially expressed miRNAs (miR162, miR166d, miR394b, miR398b, miR477a-3p and miR7122a-5p) were selected to assess reverse expression pattern between miRNAs and their target genes. Relative expression of candidated miRNAs in the cold tolerant genotype (H) was compared with the sensitive variety (Sh12) in various thermal treatments (10 and -2°C). The first strand cDNA synthesis for each miRNA was performed according to Varkonye-Gasic et al.  method. In this protocol 2μg of extracted RNA was reverse transcribed to cDNA using M-MuLV Reverse Transcriptase and stem-loop RT-PCR miRNA primers. All these primers were designed with miRNA primer designer software. To evaluate the expression of target genes by qPCR, the first strand cDNA Was synthesized using oligo-dT primers. The qPCR reaction in a total volume of 12 μl was performed on a Roter-Gene Q (Qiagen). Eeach reaction contained 2 μl of cDNA, 6 μl SYBR Premix Ex Taq (Takara) and, 0.83μM of forward and reversed primers. Three technical replicates were included for each biological repeats. The PCR amplification condition was performed according to Mousavi et al. . miRNA specific primers, used in the current study, are in Table 2. For data normalization of miRNAs, 18S rRNA  was used, while PduAct1  was applied as the internal control for target genes. The analysis was performed based on 2−ΔΔCT method .
The list of miRNA specific primers for confirmation by qPCR.
Validation of direct target genes for selected miRNAs
2μg of extracted RNA from ovary tissue of H genotype under cold stress (-1°C) treatment was ligated to a 5′RACE DNA adaptor. Random nine-mer primers were then used for cDNA synthesis. Nested PCR was performed for amplification of a cDNA fragment, which contains the cleavage site of the targets. The first PCR reaction was done using the 5′RACE outer primer and a gene-specific outer primer. The second PCR was performed using the outer PCR reaction product as template, the 5′RACE inner primer and gene-specific inner primer (Table 3). PCR products were purified using the Gene JET™ Gel Extraction kit according to the manufacturer's recommendations (Fermentase) and were sequenced by Macrogen Company (Korea) through Sanger sequencing.
The lists of primers for confirmation of predicted miRNA target genes.
The authors wish to thank the Shahrekord University for providing facilities. We are also grateful to the ‘Seed and Plant Improvement Institute (SPII)’, Karaj, Iran, and Dr Ali Imani for providing plant material used in this work. We would like also to express our gratitude to Parisa Shiran (PhD candidate at the University of Melbourne) for editing the manuscript.
The authors received no specific funding for this work.
The authors confirm that all data underlying the findings are fully available without restriction. The data have been uploaded at SRA part of NCBI http://www.ncbi.nlm.nih.gov/sra [accession number PRJNA245549].
1. Martínez-Gómez P, Sánchez-Pérez R, Dicenta F, Howad W, Arús P, Gradziel TM Almond In: Kole C, Editor. Fruits and Nuts. Springer-Verlag; Berlin Heidelberg; 2007. pp. 229–242.
2. Miranda C, Santesteban G, Royo B, De Pamplona U. Variability in the relationship between frost temperature and injury level for some cultivated Prunus species. Hort Science. 2005; 40(2): 357–361.
3. Mousavi S, Alisoltani A, Shiran B, Fallahi H, Ebrahimie E Imani A, Houshmand S. De novo transcriptome assembly and comparative analysis of differentially expressed genes in Prunus dulcis Mill. in response to freezing stress. PLoS ONE. 2014. 14; 9(8): e104541 doi: 10.1371/journal.pone.0104541[PMC free article][PubMed]
4. Baloch IA, Barozai MY, Din M. MicroRNAs: the mega regulators in eukaryotic genomes. Pure Appl Bio. 2013; 2(3): 83–88.
5. Bej S, Basak J. MicroRNAs: The Potential Biomarkers in Plant Stress Response. American Journal of Plant Sciences. 2014; 5: 748–759.
6. Zhang J, Xu Y, Huan Q, Chong K. Deep sequencing of Brachypodium small RNAs at the global genome level identifies microRNAs involved in cold stress response. BMC Genomics. 2009; 10: 449 doi: 10.1186/1471-2164-10-449[PMC free article][PubMed]
7. Liu HH, Tian X, Li YJ, Wu CA, Zheng CC. Microarray-based analysis of stress-regulated microRNAs in Arabidopsis thaliana. RNA. 2008; 14: 836–843. doi: 10.1261/rna.895308[PMC free article][PubMed]
8. Zhou X, Wang G, Sutoh K, Zhu JK, Zhang W. Identification of cold-inducible microRNAs in plants by transcriptome analysis. Biochim Biophys Acta. 2008; 1779: 780–788. doi: 10.1016/j.bbagrm.2008.04.005[PubMed]
9. Chiou TJ. The role of microRNAs in sensing nutrient stress. Plant Cell Environ. 2007; 30(3): 323–332. [PubMed]
10. Kawashima CJ, Yoshimoto N, Maruyama-Nakashita A, Tsuchiya YN, Saito K, Takahashi H, et al. Sulphur starvation induces the expression of microRNA-395 and one of its target genes but in different cell types. Plant J. 2009; 57: 313–321. doi: