The mRNA and miRNA transcriptomic landscape of Panax ginseng under the high ambient temperature

Background Ginseng is a popular traditional herbal medicine in north-eastern Asia. It has been used for human health for over thousands of years. With the rise in global temperature, the production of Korean ginseng (Panax ginseng C.A.Meyer) in Korea have migrated from mid to northern parts of the Korean peninsula to escape from the various higher temperature related stresses. Under the high ambient temperature, vegetative growth was accelerated, which resulted in early flowering. This precocious phase change led to yield loss. Despite of its importance as a traditional medicine, biological mechanisms of ginseng has not been well studied and even the genome sequence of ginseng is yet to be determined due to its complex genome structure. Thus, it is challenging to investigate the molecular biology mechanisms at the transcript level. Results To investigate how ginseng responds to the high ambient temperature environment, we performed high throughput RNA sequencing and implemented a bioinformatics pipeline for the integrated analysis of small-RNA and mRNA-seq data without a reference genome. By performing reverse transcriptase (RT) PCR and sanger sequencing of transcripts that were assembled using our pipeline, we validated that their sequences were expressed in our samples. Furthermore, to investigate the interaction between genes and non-coding small RNAs and their regulation status under the high ambient temperature, we identified potential gene regulatory miRNAs. As a result, 100,672 contigs with significant expression level were identified and 6 known, 214 conserved and 60 potential novel miRNAs were predicted to be expressed under the high ambient temperature. Conclusion Collectively, we have found that development, flowering and temperature responsive genes were induced under high ambient temperature, whereas photosynthesis related genes were repressed. Functional miRNAs were down-regulated under the high ambient temperature. Among them are miR156 and miR396 that target flowering (SPL6/9) and growth regulating genes (GRF) respectively. Electronic supplementary material The online version of this article (10.1186/s12918-018-0548-z) contains supplementary material, which is available to authorized users.


Background
Ginseng is a well-known traditional herbal medicine in north-eastern Asia and it has been used for human health for over thousands of years. The medicinal effects of ginseng have been reported as a preventative for various diseases such as cancer, cardiovascular disease, hepatotoxicity and aging related and central nervous system related neurodegenerative diseases [1]. Ginseng is a typical short day and annual plant that favors low light intensity and low temperature under 30°C during the day period. Major locations of Korean ginseng (Panax ginseng C.A.Meyer) production in Korea have migrated from mid to northern parts of the Korean peninsula to escape from the various higher temperature related stresses. Such migration is causing shortage in farmland for cultivating ginseng, which is a common and current problem for worldwide ginseng growers.
High ambient temperature is usually defined as the condition of 3 to 6°C higher than the normal temperatures. Under the high ambient temperature, vegetative growth is accelerated and results in early flowering [2]. In the growth chamber study, it has been reported that seed production was lower under higher temperature [3]. Similarly, the CO 2 exchange rate usually increases in a higher temperature, resulting in a decrease of yield of ginseng root. There are a number of studies of the thermosensory pathway in Arabidopsis, however only few studies focused on food crops and perennial plants like ginseng [4]. Unlike annual plants, perennial plants may memorize some environmental stresses as a irreversible manner, which can transmit to the following years so that some tolerances against the stress that they have been exposed to could be acquired.
To investigate the biological mechanisms in ginseng that are affected under the high ambient temperature, we have cultivated P. ginseng plants in a temperature-gradient green house (TGG). Two samples, G3 and G5, were collected that were cultivated under normal temperature (24°C) and high ambient temperature conditions (27°C) respectively. From the collected samples, we performed small RNA and mRNA high throughput sequencing.

Challenges in the analysis of transcript abundance without a reference Ginseng genome
Ginseng is a medicinal herb with high demands in east Asia. Existing studies so far mainly focused on the medicinal components of ginseng and their efficacy [5]. Ginseng is a tetraploid plant whose haploid genome is estimated to be 3.3 Gbp [6] and is known to have a genomic complexity involving a large number of major repetitive sequences, making it difficult to determine its genomic sequence. Only recently, genomic and transcriptomic studies using the high throughput sequencing technology were performed to reveal the large composition of major repeats in the genome [7], search for root specific miRNAs [8,9] or understand the biological metabolisms in the Ginseng plant, such as the synthesis of ginsenosides [10] or the biological response to autotoxins [11]. However, to our knowledge, there is no study on how genes and miRNAs respond to the ambient hight temperature in P. ginseng.
The plant itself is largely understudied compared to its medical significance and demands worldwide. The lack of experimentally validated genes, proteins and non-coding RNAs of the P. ginseng genome hinders the analysis of gene association, gene regulation or annotation of pathways or coding/non-coding transcripts. As of now, only 29 miRNAs specific to P. ginseng are deposited in miRBase [12]. Also, databases that archive curated or predicted genomic information of the ginseng plant are not present or not accessible.
The analysis of transcriptomic data mainly relies on aligning reads to a reference sequence, which is not available for the non-model plant species P. ginseng plant. Traditional sequencing was initially used to obtain root expressed sequence tags (ESTs) [13], which are very low in quantity. Studies based on high throughput sequencing data performed functional gene analysis based on contigs assembled from their sampled DNA or RNA libraries. In addition, the P. ginseng varieties differ between the contig sequences submitted to public databases. Thus, the quality of contigs assembled from a small number of samples and the sequence variance across sub-species makes it difficult to utilize them directly for sample or condition specific gene analysis. For genes, only partial gene sequences, such as ESTs [13] and assembled contigs [11] are available.
With the aforementioned challenges in preforming transcriptomic analysis of the P. ginseng plant, we have made several efforts to overcome them by the means of understanding the biological response of the plant under the influence of the high ambient temperature. The major tasks performed are categorized into three parts, which are described in more detail in the "Methods and Results" section in order as described below. First, we have collected a total of 65 paired end RNA-seq samples from the Sequence Read Archive (SRA) for the assembly of transcripts (Additional file 1). Among the assembled contigs, candidate protein coding transcripts were annotated based on detecting long open reading frames (ORF) and searching their closest orthologs in the NR [14] and UniRef90 [15] protein databases. Their expression levels were quantified by the high ambient temperature treated samples. Second, miRNAs were searched and quantified using our small RNA-seq data with the aid of miRBase and a de-novo miRNA detection tool. Third, associations between miRNAs and their candidate target transcripts were identified for functional analysis under the influence of the high ambient temperature. For validation purposes, several transcripts were selected, which sequences were validated by performing RT-PCR.

Materials and experimental design
The goal of our study is to investigate biological mechanisms in ginseng that are affected under the high ambient temperature. In addition, we investigated the regulatory roles of miRNAs on mRNA transcripts under the high ambient temperature.

Plant materials
P. ginseng (cv. Yeonpung) plants were cultured in a TGG located in Eumsung, Chungbuk of Republic of Korea. The micro climate inside of TGG was regulated for 24 h in the year round as described in [16] and [17]. The maximum temperature gradient in TGG was reached at 6°C and the leaf samples for RNA extraction were obtained from two different locations of TGG, G3 and G5 which is near the front entrance and in the middle location of TGG, respectively. At the time of collecting the plant samples, the mean temperature of the front entrance of the TGG was 24°C, which was the same as the outdoor temperature. The mean temperature in the middle location of the TGG was 27°C. The temperature of the TGG was measured throughout the year at several locations within the TGG, which is shown in Additional file 2. The average temperature difference between G3 and G5 was 3°C. Leaves of three years-old ginseng were obtained from three to five different plants in each G3 and G5 that were frozen immediately on site with liquid nitrogen on August 3 2016. Total RNAs of ginseng leaves were extracted using Qiagen RNA extraction kit (Qiagen USA). We used homogenizer to break the cell walls. Extraction of total RNA was repeated for three times where the samples for RNA sequencing were mixed in all three total RNAs repetitively. The sequencing were carried by the Macrogen company in Seoul, Korea.

Methods
To identify the high ambient temperature responsive mRNAs and regulatory miRNAs, our analysis was performed in four parts as shown in Fig. 1. First, contigs were assembled to search and annotate putative protein coding genes. Second, known, conserved and novel miRNAs with valid miRNA precursors were predicted from our assembled contig sequences. Third, negatively correlating expression levels between miRNAs and their target genes were searched. At last, selected genes that were predicted to have functional roles or are responsive to high ambient temperature were validated by reverse transcripted (RT) PCR.
In the first part of performing the transcriptomic analysis of the P. ginseng transcriptome data, we assembled and annotated high quality contigs utilizing 65 paired-end RNA-seq samples collected from SRA. For samplespecific contig sequences, we further corrected the assembled contigs using our RNA-seq samples. Here, gene expression was also quantified using our RNA-seq data using the assembled contigs. In the second part, miRNAs and their precursors were predicted using the assembled contigs. In the last part, we identified miRNA and their target genes that have negative correlating expression levels. The result of the workflow analysis is summarized in Table 1.

Part1: Contig analysis for sample specific gene annotation
For species without a fully annotated reference genome like P. ginseng, de novo assembly must be performed for identifying genes and for further downstream gene expression analysis. Previous studies performed EST sequencing, whole genome sequencing and transcriptome sequencing. None of them are enough to complete the complex Ginseng genome. Furthermore, the subspecies also differ between the studies. Hence the studies performed sample specific contig assembly utilizing only a small number of samples (i.e., a maximum of up to 12 RNA-seq samples) [11]. Since the sequence variation between the existing contigs and EST sequences are large, the RNA-seq alignment are poor, which may not be used for designing primer sequences during sequence validation. In this study, we utilized as many available P. ginseng high throughput samples as possible to assemble high quality contigs, which are then subsequently corrected by building consensus sequences using sample specific RNA-seq reads. In this manner, we were able to mask out sequence variations caused by sub-species differences and perform an effective downstream analysis to reveal biological mechanisms under the high ambient temperature.
The assembly process generated a total of 100,672 contigs. From the contigs, a total of 34,497 protein coding genes were annotated through the sequence similarity search of identified open reading frames (ORF) against the non-redundant protein NR and UniRef90 databases. Differentially expressed gene (DEG) analysis and the gene set enrichment analysis were performed to identify enriched GO (Gene Ontology) terms and pathways.

De novo assembly
From the SRA database, 65 paired end RNA-seq data sampled from P. ginseng root, stem and leaf tissues were collected. Reads from the collected data were pooled together and assembled into the new reference transcriptome sequences. The transcriptome assembly was performed using Trinity (version 2.1.1) [18]. To be able to process the large RNA-seq data sets, we used the Trinity in silico normalization. All parameters were set as default. As a result 62.9 million base pairs and 744,507 contigs were assembled. We further performed scaffolding using CAP3 [19] to acquire a total of 613,448 scaffolded contigs. To compensate for sequence variations caused by sub-species difference, our transcriptome RNA-seq reads were used to build sample specific consensus sequences. As a result, the N50 was 845 bp and the average contig length was 625 bp. The GC content was 36.01%.

Measuring contig expression levels
To measure the expression levels of contigs, we aligned the paired RNA-seq reads from G3 and G5 samples using Bowtie2 [20] to the contigs. For the alignment options, we discarded unpaired and discordant paired reads. Subsequently, we computed FPKM (Fragments per kilobase of exon per million fragments mapped) of each contig for quantifying their expression levels. Contigs with expression level below 1 FPKM were not considered for further analysis.

Annotation of expressed transcripts in contigs
For the functional annotation of contigs, we performed similarity search of the contigs against the NR and UniRef90 non-redundant protein databases. We first extracted putative coding regions from the assembled contig sequences using TransDecoder, which is a Trinity software package. A total of 34,497 protein sequences with a minimum length of 100 amino acids were collected. These peptide sequences were searched using BLASTP with an E-value < 10 −6 against the protein databases. Using the UniProtKB [21] and GO database [22], we further annotated genes with their associated KEGG IDs and GO terms. A total of 26,460 genes were annotated. The list of annotated genes including their searched orthologs in other plant species are provided in the Additional file 3.

Identifying differentially expressed genes
To identify genes that are differentially expressed under the high temperature, we used DEGseq [23], which detects DEGs based on the MA-plot method with a random sampling model. A stringent cut-off value of P < 0.001 was for DEG detection. Multiple testing correction has been applied during the DEG detection process (i.e., q-values are included in the Additional file 4, which uses the Benjamini-Hochberg procedure). A total of 1256 DEGs were identified where 1028 and 228 genes were significantly up and down regulated in G5 compared to G3, respectively. The MA-plot is shown in Fig. 2. As shown, most of the DEGs are induced under the high ambient temperature. We categorized differentially expressed genes into three groups: Photosynthesis related genes, temperature responsive genes and flowering/developmental related genes. The photosynthesis related genes were down  (Table 2). On the other hand, temperature responsive and flowering and growth related genes were induced under the high ambient temperature ( Table 3). The expression level and significance of differentially expressed genes are provided in the Additional file 4.

Gene set enrichment functional analysis
The functional analysis of DEGs were performed by the GO and the pathway enrichment analysis. For both enrichment analysis, we performed Fisher's exact test with P < 0.05. For the GO enrichment test, we used genes with annotated GO terms as the background set. A total of 45,075 unique GO terms were present in the background set. The list of DEGs with annotated GO terms were used as the query set, where 3,966 unique GO terms were present. 73 GO terms showed a significant enrichment. The two most highly enriched GO terms were ' ATP binding' related GO terms (i.e., ' ATP binding' and 'Chloroplast thylakoid membrane'). Light harvesting related GO terms, such as 'Protein-chromophore linkage' , 'Photosystem I/II' , 'Photosynthesis' , 'Light harvesting in photosystem I' and 'Pigment binding' were next highly enriched. The list of enriched GO terms are provided in Additional file 5(a).
For the pathway enrichment test, we used genes with annotated KEGG pathway IDs (i.e., KO IDs) as a background set. A total of 2,415 unique KO IDs were present in the background set. The list of DEGs with annotated KO IDs were used as the query set, where 260 unique KO IDs were present. 22 pathways showed a significant enrichment of DEGs. There were a significant number of DEGs in two photosynthesis pathways (i.e., 'Photosynthesis' and 'Photosynthesis -antenna proteins'), and the 'Plant hormone signal transduction' pathway. A majority of genes, including DEGs that mapped to photosynthesis pathways, were significantly down regulated. Within the 'Plant hormone signal transduction' pathway, the genes along its sub-pathways related to Auxin (i.e., AUX1, TIR1, ARF etc.), Giberrellin (i.e., GID1, GID2, DELLA, TF) and Abscics acid (i.e. PYR/PYL, PP2C and SnRK2), were significantly up regulated. The three sub-pathways are shown to induce cell enlargement, plant growth, stem growth or germination and stomatal closure and seed dormancy respectively. The list of enriched pathways is shown in the Additional file 5(b).
The first cDNA strand was synthesized by RT-PCR which was carried out with RT-PCR kit of Invitrogen using both gene specific primers. In addition, we used a gene specific primer for upper primer and oligo (dT) for down primer to amplify 3 -end regions and to compare with contigs we constructed as described.
Among the six sequences, the contig sequence aligning to R3H showed a clear band during the validation. As a result, two sequences with length of 665 and 657 nt were detected (Additional file 6 -Sequence 1 and Sequence 2). The sequences were BLASTN searched in    well aligns to the Phytoclock Like 1 (PCL1) protein sequence, which encodes a MYB family transcription factor with a single MYB DNA-binding domain (type SHAQKYF) that is unique to plants and is essential for circadian rhythms, specifically for transcriptional regulation within the circadian clock. It is also known as the symbol LUX, which is required for normal rhythmic expression of multiple clock outputs in both constant light and darkness. It is co-regulated with TOC1 and seems to be repressed by CCA1 and LHY by direct binding of these proteins to the evening element in the LUX promoter [25]. In A. thaliana, PCL1 is also known to encode a novel GARP protein that is essential for the circadian clock [26].
While not discussed with further detail, three other genes were successfully synthesized, which were the contigs that well aligned to MET1, DDM1 and Chlorophyll A/B binding protein encoding gene sequences (Additional file 6 -Sequence 3, 4 and Sequence 5).

Part2: Sample specific small non-coding RNA analysis
To investigate how small RNAs interact with genes and regulate gene expression level under the high ambient temperature, we searched and predicted miRNAs from the small RNA-seq data. Among the annotated ginseng miRNAs in miRBase, 6 miRNAs were searched with significant expression levels. A total of 214 miRNAs conserved in other species were detected. In addition, 60 novel miRNAs with associated template precursors were predicted from our assembled contigs. Interestingly, functional miRNAs showed down regulation under the high ambient temperature. Consistently, their target genes were up regulated.

Landscape of small RNA transcripts in P. ginseng
We found that the landscape of small RNA transcripts significantly differed in G3 and G5 as shown in Fig. 3. In G3, small RNAs of 21 nt were most abundant while it was not the case in G5. Interestingly, longer small RNAs were expressed more in G5, which is a quite strong trajectory. The expression level of miRNAs was quantified by the miRNA read counts normalized by each sample size. As of now, our scientific community has limited knowledge on the biological mechanisms of P. ginseng and the landscape of small RNA transcripts needs in-depth investigation in the future. We discussed the regulatory roles of miRNAs in the "Discussion" section.

Predicting miRNAs
A total of 329 miRNAs were detected from the small RNA-seq samples. Eight known, 243 conserved and 78 novel miRNA transcripts were annotated as shown in Fig. 4.

Identification of known miRNAs
Currently, 29 miRNAs, with precursor sequences, of the P. ginseng plant are annotated in miRBase [12], which were archived by the previous study [27]. Among them, five miRNAs were commonly expressed in G3 and G5 samples. All five miRNAs were up-regulated in G5. pgi-miR6140a was induced by more than 2-fold as shown in Fig. 5a. Furthermore, we tested if miRDeep2 [28] reported the known miRNAs, and their precursor sequences by providing the annotation information of the known miR-NAs from miRBase. Among the five miRNAs, four miR-NAs with associated precursors were detected. a b Fig. 3 Features of the small RNA population by transcript length. a The number of unique transcripts were the greatest for 21 nt transcripts. G3 and G5 both showed a peak for 24 nt transcripts. The population of unique reads increased with the length of transcripts in G5, whereas it declined in G3 for transcripts longer than 25 nt. b The expression level of small RNA transcripts similarly showed a spike at 21 nt. Also, the expression level increased with the length of small RNA transcripts specific to G5

Identification of novel miRNAs
Using miRDeep2, we searched for novel miRNA precursors and associated mature miRNA sequences from the assembled contigs. First, miRDeep2 aligned the small RNA-seq reads to the contigs without allowing any mismatches. Then it predicted RNA secondary structure from the RNAfold result, which predicts and scores RNA secondary structures of potential miRNA precursors. Potential miRNA precursors were determined by the following three criteria: (1) contigs that fold into an unbifurcated hairpin, (2) contigs that consist of mature (or core miRNA region), loop and star regions, and (3) contigs where 60% of reads map to the mature sequence region. As a result, we identified 101 and 82 novel miRNA precursors in G3 and G5, respectively. In total, 139 novel miRNA precursors were identified where 44 were commonly present in G3 and G5.

Identification of conserved miRNAs
To identify conserved miRNAs in our G3 and G5 samples, we aligned the small RNA reads against the miR-Base miRNA precursors allowing up to 2 mismatches. The short reads were aligned to the hairpins in miRBase that are annotated from more than 206 species. Species with less than 28 annotated miRNA hairpins were filtered out from the alignment, which accounts for 0.1% of the total annotated hairpins in miRBase. The organisms were sorted by the number of aligned reads to the associated organism's hairpin sequence. The top 16 organisms were all plant species. Among them, the P. ginseng had the highest hits, which was very surprising in respect to the small number of annotated miRNA hairpins (i.e., 32 hairpin sequences). The following plants followed in descending order: L. usitatissimum, M. esculenta, T. cacao, R. communis. Similar to the known and novel miRNAs, the conserved miRNAs also showed a strong down-regulation pattern in G5. These conserved miRNAs in multiple organisms may share common functions, which can be valuable knowledge in P. ginseng. The fold-changes of miRNAs are shown in Fig. 6.

Part3: mRNA-miRNA integrated analysis
To search for gene regulatory miRNAs, we identified target contigs of predicted miRNAs using psRNATarget [29]. psRNATarget is a small RNA target analysis tools for plant species. For the miRNA sequence input, we provided the a b c known, conserved and novel miRNAs that are present in our small RNA-seq data. For the miRNA target candidates, we provided the protein coding contig sequences. Target contigs were selected as candidates when contigs with target sites have significant alignment with a miRNA and have strong target-site accessibility based on unpaired energy between nucleotides of the target site and the miRNA. Finally, we selected miRNA-target gene pairs only when the expression level were negatively correlated. As a result, a total of 592 miRNA and their target contig pairs with negative correlation were identified. Among them, a majority of 583 miRNAs were down regulated while their target genes were up regulated. Only 12 pairs were present where miRNAs showed mRNA suppression. We further investigated the miRNA target regions in the contigs. The target site statistics are grouped by genomic Fig. 6 Expression fold change of identified miRNAs. The expression level fold change of the identified 156 miRNAs are shown in log 2 scale. The miRNAs are sorted in ascending order of the log 2 fold change regions such as 3 UTR, ORF or 5 UTR region. From our assembled contigs, non-ORF regions that are located to the left of the starting ORF are annotated as 5 UTR regions. Similarly, the non-ORF regions located at the right of the last ORF region are annotated as 3 UTR regions. In summary, miRNA target sites seem to be distributed somewhat evenly across different regions within the gene body as shown in Fig. 7.

Discussion
From hight-throughput sequencing data collected from the P. ginseng plant, we assembled transcript contigs and identified three types of gene sets that showed response to the high ambient temperature: 1) Heat stress associated genes, 2) Photosynthesis and respiration associated genes and 3) Flowering associated genes. Genes related to development and flowering were significantly induced under the high ambient temperature. Among all genes, the Flowering promoting factor 1 (FPF1) was the most up-regulated gene with 29-fold change in G5 compared to G3. The Phytochrome interacting factor 4 (PIF4) was also present in our contigs list, which was also upregulated by 2.2-fold in G5. Interestingly, the target genes of PIF4, FT [30], was also up-regulated by 1.94-fold in G5. Brassinazole-resistant1 (BZR1) is known to be a cointeracting gene that interacts with PIF4 [31], which was also up-regulated by 2.7-fold.
A strong gene expression trend was observed between temperature (i.e., heat and cold) responsive and photosynthesis related genes. 14 low and cold temperature responsive genes were expressed, of which five were detected as up-regulated DEGs. The remaining genes were not DEGs but were all induced by some level. For example, the Low temperature induced 78 (LTI78) protein encoding gene was induced by more than 4-fold in G5. LTI78 is reported to accumulate in all genotypes in response to low temperature and drought and it is always present when plants were freezing tolerant in the A. thaliana organism [32]. Cold and circadian rhythm related genes, such as CDSP1, CCR1 and CCR2, are induced by 4 to 8-folds in G5. A total of 46 heat related genes were expressed, of which 40 were induced. Among the induced heat associated genes, 19 were DEGs. Only two down-regulated DEGs were detected. Most of the genes were Heat shock proteins where Hsp70 and Hsp90 were the most abundant ones. On the contrary, a majority of photosynthesis related genes were repressed. Among the 21 photosynthesis related genes, 16 were down-regulated, of which 3 were DEGs. Only one up-regulated DEG was detected. Phototropin, Photosystem I and II were the most abundantly found genes in the photosynthesis related gene set. Collectively, most of the temperature responsive genes were induced under the high ambient temperature, whereas photosynthesis related genes were down-regulated. The relation between temperature response and photosynthesis activity related genes showed opposite expression trends.
miRNAs that are associated to stress response and floral development are observed to be down-regulated under the high ambient temperature. Three stress responsive miRNAs conserved in other plant organisms (i.e., miR394, miR399a and miR8175) were expressed in large quantity that were all down-regulated in G5. They are known to be associated with salt and drought response in plants. The development associated miRNAs were also downregulated. The miR156 is reported to be induced under stressful conditions in plants that adopted stress tolerant phenotypes [33]. SPL6 and SPL9 are among the targets of miR156, and when suppressed by miR156, the flowering time is delayed until a suitable conditions is met [34]. We observed that miR156 was repressed in G5 by 3-fold, whereas its target genes, SPL6 and SPL9, were induced correspondingly, which may imply that the flowering mechanism has been promoted under high ambient temperature. The development associated miR-NAs miR159, miR319, miR390a and miR396. miR396 were down-regulated. They are known to antagonize the expression of the growth regulating factor (GRF) gene, which is a transcription factor. All of the GRF genes were induced, including the GRF2 to GRF12 genes. Collectively, the growth of the plant may have been accelerated by several biological mechanisms that includes the repression of miRNAs that target growth promoting genes [35].
In previous studies, the expression levels of small interfering RNAs (siRNAs) are shown to be dramatically reduced with the rising temperature [36]. However, the a b direction of expression pattern of miRNAs in response to stress are more dynamic. In the study [37], the miR169 was shown to be down-regulated in response to drought in A. thaliana, but up-regulated in rice.