Skip to main content


Global analysis of uncapped mRNA changes under drought stress and microRNA-dependent endonucleolytic cleavages in foxtail millet

Article metrics



mRNA degradation plays an important role in the determination of mRNA abundance and can quickly regulate gene expression. The production of uncapped mRNAs, an important mechanism of mRNA degradation, can be initiated by decapping enzymes, endonucleases or small RNAs such as microRNAs (miRNAs). Little is known, however, about the role of uncapped mRNAs in plants under environmental stress.


Using a novel approach called parallel analysis of RNA ends (PARE), we performed a global study of uncapped mRNAs under drought stress in foxtail millet (Setaria italica [L.] P. Beauv.). When both gene degradation (PARE) and gene transcription (RNA-sequencing) data were considered, four types of mRNA decay patterns were identified under drought stress. In addition, 385 miRNA–target interactions were identified in the PARE data using PAREsnip. The PARE analysis also suggested that two miRNA hairpin processing mechanisms—loop-last and loop-first processing—operate in foxtail millet, with both miR319 and miR156 gene families undergoing precise processing via the unusual loop-first mechanism. Finally, we found 11 C4 photosynthesis-related enzymes encoded by drought-responsive genes.


We performed a global analysis of mRNA degradation under drought stress and uncovered diverse drought-response mechanisms in foxtail millet. This information will deepen our understanding of mRNA expression under stressful environmental conditions in gramineous plants. In addition, PARE analysis identified many miRNA targets and revealed miRNA-precursor processing modes in foxtail millet.


Transcript abundance is modulated by transcript synthesis and degradation rates. In recent years, genome-wide profiling methods, such as RNA-sequencing (RNA-seq) [1] and gene-chip analysis [2], have been used to study mRNA expression and to identify genes expressed in specific tissues and developmental processes [3] and in response to environmental stimuli [4]. Such studies are generally designed to capture aspects of steady-state mRNA abundance. In addition to transcriptional regulation, however, an understanding of gene expression networks obviously requires data on mRNA degradation and other patterns of mRNA expression regulation. Recent research has indicated that proper mRNA degradation is a component of cellular homeostasis maintenance and contributes to the precise adjustment of gene expression levels in response to various extracellular stimuli [5]. Several highly conserved pathways for mRNA degradation exist in eukaryotes. One such pathway proceeds in the 3′ to 5′ direction, with mRNA decay beginning with deadenylation catalyzed by mRNA deadenylases [6, 7]. Another mRNA degradation pathway, from 5′ to 3′, is often initiated with cleavage of the monomethyl guanosine (m7G) by the Dcp2 decapping enzyme [8, 9]. Dcp2 mainly cleaves m2,2,7G-capped RNA and m7G-capped RNA, while unmethylated capped RNA is a poor substrate [10]. In addition, some internal cleavages, such as those involving the RNA-induced silencing complex (RISC) directed by microRNAs (miRNAs) or small interfering RNAs, can initiate mRNA degradation [11, 12]. miRNA-mediated degradation mainly occurs in the miRNA:mRNA pairing binding region [13]. Similar mRNA decay mechanisms have been reported in plants, with some studies indicating that mRNA degradation is important for the control of gene expression during growth, development and many physiological transitions [1416]. The most recent study on this topic revealed that RNA decay pathways function mainly via protecting transgenes and endogenous genes from inappropriate posttranscriptional gene silencing to regulate gene expression, and act as safeguards of plant development in Arabidopsis [17].

miRNAs are ~21-nucleotide non-coding RNAs that regulate gene expression by base-pairing to their targets, resulting in gene degradation or translational inhibition in most eukaryotes. In plants, the miRNA is almost completely complementary to its target gene and mediates cleavage of the target at the center of the paired region [1820]. The high-throughput method used to discover miRNA targets in plants, which is based on computational prediction using a set of pre-defined rules [21, 22], leads to a large number of false positives. A more reliable, experimentally based approach involves the use of gene-specific 5′ rapid amplification of cDNA ends (RACE) to validate miRNA-target pairs [23]. Using this method, however, every single predicted target must be independently verified, which is labor-intensive, time-consuming and expensive. Although many miRNA-target pairs have been predicted, only a small fraction has been experimentally confirmed. One encouraging development is parallel analysis of RNA ends (PARE), a novel approach for identification of miRNA targets that is both high-throughput and reliable [24]. Using this method, large-scale miRNA–target pairs have been validated in species such as Arabidopsis [25, 26], rice [27, 28] and grapevine [29].

Most miRNA/miRNA* duplexes released from typical miRNA precursor stem-loops undergo two cycles of cleavage by endonucleases: one at the loop-distal position and the other at the loop-proximal position [30]. In animals, the two steps are spatially separated and completed by two different enzymes; specifically, duplexes are first cleaved by Drosha at the loop-distal position in the nucleus and then by Dicer at the loop-proximal position in the cytoplasm [31, 32]. In plants, the two sequential cleavages are both carried out by a Dicer-like (DCL) enzyme (DCL1) in the nucleus but the DCL1 sequential cleavage site is still poorly understood. Recent research has shown that PARE data can be used to probe patterns of miRNA hairpin processing in plants [33, 34].

Foxtail millet (Setaria italica [L.] P. Beauv.; Poaceae) is an important grass crop species widely planted in China. The genome of the foxtail millet cultivar Yugu1 has been sequenced recently [35]. Because of its trivial size and short life cycle, foxtail millet is an ideal model species for Panicoideae crops and C4 photosynthetic species [3638]. Over the long period of its improvement and domestication, foxtail millet has gradually adapted to semiarid and arid climates. Because of its excellent drought tolerance and water-use efficiency, foxtail millet is an ideal material to investigate the mechanisms of drought tolerance in plants. To deepen our understanding of mRNA degradation and stress response mechanisms in plants, we therefore used the PARE deep-sequencing approach in this study to analyze uncapped mRNA transcripts under drought stress conditions in foxtail millet. PARE was also used to identify potential targets for miRNA-directed cleavage and to reveal multiple novel examples of miRNA precursor processing in foxtail millet.


Overview of PARE-seq data in foxtail millet

To characterize mRNA degradation changes during the drought stress response, we profiled uncapped transcripts using PARE-seq (see Additional file 1 for details) in 14-day-old foxtail millet seedlings. PARE libraries were prepared from seedlings subjected to polyethylene glycol (PEG)-simulated drought (Dd) and control (Dc) conditions, with two independent biological replicates for each group. We generated approximately 39 million and 45 million 50-bp reads from the control and drought-treated seedlings, respectively (Additional file 2). After removing repeats/transposons [39] and known non-coding RNA (rRNA, tRNA, small nuclear RNA and small nucleolar RNA) [40] sequences as described in the Methods, the raw reads were mapped to the Yugu1 reference genome [35] using bowtie. The mapped reads, which accounted for about 68.06–76.35 % of the total reads from different PARE-seq libraries, were used to measure gene degradation levels calculated as reads per million (RPM) (see Methods). Because comparisons of biological replicates showed that their expression values were highly correlated (average R 2 = 0.97; Additional file 3), we took the average RPM of the replicates as the degradation level.

We detected 26,802 genes with uncapped transcripts in at least one of the samples according to the PARE analysis. To explore whether 5′ to 3′ mRNA degradation is associated with the drought response, we compared the degradation levels of these genes and found that 1553 genes exhibited significant changes after drought treatment, with two-thirds of them up-regulated (Fig. 1a). This result demonstrates that the process of mRNA degradation plays an important role in the drought response. For example, four late embryogenesis abundant (LEA) domain-containing protein genes were included in the list of the 10 most significantly up-regulated genes (Table 1). This finding is consistent with previous reports that LEA proteins are associated with cellular tolerance to dehydration induced by salinity, freezing or drying [4143]. In addition, we found that genes for two MYB-family and one WRKY-family transcription factor(s), which are regulators of various plant developmental and physiological processes in response to drought stress [44, 45], were among the 10 most significantly down-regulated genes in our experiments (Table 1).

Fig. 1

Drought stress-responsive genes at the degradation and transcription levels. a Number of genes showing significant degradation level changes after drought treatment. b Venn diagram showing significantly changed genes at the degradation and transcription levels under drought treatment. c Correlation between fold-change values at the degradation and transcription levels after drought stress. Correlation values (R 2) are Pearson’s product–moment correlation coefficients. Up arrow and down arrows represent up- and down-regulation after drought treatment, respectively. d Coverage of parallel analysis of RNA end reads and RNA-sequencing reads on selected genes; each had two replicates. The genes Si021866m and Si013398m belong to type-I, Si036695m and Si022331m to type-II and Si016654m and Si023461m to type-IV. “+PEG” and “−PEG” represent samples with (+PEG, drought-treated) and without (−PEG, control) PEG treatment. The y-axis represents the normalized read depth

Table 1 Function annotations of top 10 up-regulated and down-regulated genes in PARE analysis

Different transcript regulation patterns during the drought response

mRNA synthesis and degradation both affect mRNA abundance. To study the connection between mRNA synthesis and degradation in transcript regulation under drought treatment, we next measured genome-wide gene expression levels using our previously published RNA-seq data [4]. A total of 2824 genes showed significant expression level changes after drought treatment, with almost equal numbers up-regulated and down-regulated. In addition, we found that 34.6 % (1126) of genes displayed significant changes in both expression and degradation levels after drought treatment (Fig. 1b), implying the existence of different transcript regulation patterns during the drought response. To reveal the overall trend of mRNA synthesis and degradation changes, we calculated the fold changes in gene expression and degradation levels for 19,814 genes identified by both RNA-seq and PARE data analysis. We uncovered a positive correlation (Pearson’s correlation: R 2 = 0.56) between the fold-change values of transcription and degradation levels under drought conditions (Fig. 1c). We then classified these genes into nine categories with a log2 fold change of ±1 as the threshold, and found that 86.4 % (17,118 genes, class E in Fig. 1c) of the genes were unchanged. The remaining 2696 genes, comprising the other eight classes (classes A, B, C, D, F, G, H and I in Fig. 1c), were inferred to be involved in the drought response and regulated by either RNA transcription or RNA degradation. We further classified these drought-responsive genes into four types according to their characteristic changes (Table 2) and performed a Gene Ontology (GO) enrichment analysis using the WEGO online tool ( [46] (Additional file 4). We were thus able to aggregate genes with different transcript regulation patterns and view their functional classifications. Notably, only two genes were found to belong to type III (Table 2; Classes A and I in Fig. 1c), in which transcript and uncapped transcript abundance showed opposite trends. As a consequence, no further analysis was performed on type-III genes.

Table 2 Four different transcript regulation patterns in drought stress response

Type-I genes (Table 2; classes C and G in Fig. 1c), which were characterized by transcript and uncapped transcript abundances changing in the same direction after drought treatment, were enriched in catalytic and various oxidation-related enzymes (oxidoreductase, antioxidants and peroxidase) in the molecular function category as well as the biological process subcategories ‘metabolic process’, ‘response to stimulus’ and ‘response to stress’ (Additional file 4). PARE and RNA-seq read coverage is shown for two examples in Fig. 1d. One of these genes, Si021866m, encodes a DNAJ heat shock N-terminal domain-containing protein. Its homologous gene in Arabidopsis thaliana plays an important role in cellular stress sensors [47]. The other example, Si013398m, encodes a transcription factor belonging to the MYB family reported to play crucial roles in plant responses to abiotic stress [44, 48].

Genes in the type-II category (Table 2; classes B and H in Fig. 1c), comprising genes showing significant changes in uncapped transcript abundance after drought treatment but no changes in transcript abundance, were typified by Si036695m, a No Apical Meristem (NAC) transcription factor, and Si022331m, a bZIP transcription factor (Fig. 1d). Genes related to transcription factors, transcription regulators, pigmentation and regulation processes (under the biological and cellular metabolic categories) were overrepresented in this group (Additional file 4).

Type-IV genes (Table 2; classes D and F in Fig. 1c) showed significant changes in transcript abundance while their uncapped transcripts remained unchanged after drought treatment. Genes associated with membranes were enriched in the cellular component category, and various important catalytic activities (hydrolase, lyase and isomerase) in the molecular function category and ‘response to chemical stimulus’ and ‘oxidation reduction’ in the biological process category were heavily represented. Si016654m, a representative of type IV (Fig. 1d), encodes an arginine decarboxylase (ADC) protein. Previous studies have found that ADC is involved in responses to salt, drought and other abiotic stresses [4951], with other investigations revealing that transgenic ADC rice plants show increased biomass under salinity-stress conditions compared with non-transformed controls [52]. Si023461m, another type-IV example, encodes ribulose bisphosphate carboxylase, the crucial enzyme in photosynthesis and photorespiration. Sharkey et al. [53] showed that mild water stress affects ribulose bisphosphate carboxylase activity in intact leaves.

Detailed functional annotations for type-I, −II and -IV genes are given in Additional file 5. Taken together, our data demonstrate that different transcript regulation patterns exist during the drought response and may be correlated with gene function.

Sequence characteristics correlated with different mRNA decay patterns

Sequence characteristics have been reported to contribute to uncapped mRNA abundance [15, 16]. To reveal the characteristics of mRNAs with different decay patterns, we calculated the mRNA lengths, GC contents, minimal folding free energy indexes (MFEIs) [54] of secondary structures, untranslated region (UTR) lengths and intron numbers for type-I, −II and -IV genes (Fig. 2 and Additional file 6). The results showed that the mean lengths of 5′ UTRs, 3′ UTRs and mRNAs of type-I, −II and -IV genes were significantly greater (P < 0.001) than those of all genes and 898 (average number of I, II and IV genes) randomly selected genes (Fig. 2a–c). Moreover, the mean MFEI values of 5′ UTRs, 3′ UTRs and mRNAs of type-I, −II and -IV genes were significantly lower than those of all genes and randomly selected genes (Additional file 6). In contrast, no significant differences were found for the mean GC contents of 5′ UTRs and 3′ UTRs among these gene groups (Additional file 6). Overall, the lengths and MFEIs of UTRs and mRNAs correlated with the structural features of genes involved in the drought response.

Fig. 2

Gene transcript features and sequence motifs contributing to different mRNA decay patterns. a-d A display of mRNA length, 5′ UTR length, 3' UTR length and number of introns for different gene types. “I”, “II”, “IV”, “R” and “A” represent type-I, −II, −IV, randomly selected genes and all genes, respectively. “***” indicates statistically significant difference at P-value < 0.001 (Student’s wilcox-test). e-h Enriched motifs (E-value < 0.001) in the 5' UTRs of type-I (e), type-II (f and g) and type-IV (h) genes

Recent studies have revealed a relationship between introns and mRNA stability [55, 56]. In our study, the average number of introns in type I and IV genes was higher than in all genes and randomly selected genes (P < 0.001), whereas no significant differences were found between type-II genes and the latter gene categories (Fig. 2d). This result indicates that the number of introns has an obvious influence on drought-responsive genes, and this effect is related to mRNA synthesis regulation rather than degradation regulation.

Previous research has shown that enriched motifs in 5′ and 3′ UTR regions can affect mRNA stability [15, 16, 57]. We therefore used an integrated motif-discovery program (MEME) [58, 59] to identify possible motifs located in the 5′ and 3′ UTRs of type-I, −II and -IV gene transcripts. We identified significantly enriched motifs (E-value < 0.001) in the mRNA 5′ UTRs: one in type I (Fig. 2e), one in type IV (Fig. 2h) and two in type II (Fig. 2f and g), implying the existence of some specific regulators or regulatory mechanisms in their 5′ UTRs. Next, we analyzed these enriched motifs using Tomtom tool. The best-matched motif to ‘[AG]G[AG][GA][GA][AG][AG]GA[GA]’ (Fig. 2e) was RNCMPT00044 (P = 2.31528e-05). The protein Poly (rC)-binding protein 2 (PCBP2) was reported to bind the RNCMPT00044 motif [60] and appears to be multifunctional [61, 62]. RNCMPT00019 was the best-matched motif to ‘[AG][GA][AG][GA][GA][AG]AGAA’ (Fig. 2f) (P = 0.002) and RNCMPT00088 was the best-matched motif to both ‘T[CT][TC][TC][CT][TC]TCT[CT]’ (Fig. 2g) (P = 2.16363e-06) and ‘[GA][GA]AG[AG][AG][GA][GA]A[GA]’ (Fig. 2h) (P = 6.99643e-05). The binding protein for both RNCMPT00019 and RNCMPT00088 is serine/arginine-rich splicing factor 10 (SRSF10) [60]. SRSF10 is known to function as a sequence-specific splicing activator [63] and can promote both exon inclusion and exclusion in chicken cells [64]. This suggests that some RNA-binding proteins may play a role in the behavior of these gene classes. In contrast, no significantly enriched motifs were found in the 3′ UTRs of type-I, −II and -IV genes.

Identification of endogenous miRNA cleavage targets

In plants, miRNAs play key roles in many developmental events and regulate their target transcripts through two modes of action: degradation and translation inhibition [6567]. In this study, PARE-seq was used to identify the cleavage sites of targets mediated by miRNA-programmed RISCs [24]. Using the PARE data, we identified 385 putative miRNA-guided cleavages in foxtail millet (Additional file 7); six prominent examples were selected for detailed discussion (Fig. 3). miR160 guides cleavage within the coding regions of Si005991m (993 reads across the Dc PARE libraries, Fig. 3a), Si034525m (322 reads, Fig. 3b) and Si016509m (117 reads, Fig. 3c). The proteins encoded by Si005991m, Si034525m and Si016509m are homologs of A. thaliana auxin response factor 16 (ARF16). A previous study determined that ARF16, targeted by miR160, controls root cap cell formation in A. thaliana [68]. miR169c guides cleavage within the coding region of Si037045m (75 reads across the Dc PARE libraries, Fig. 3d), while nov-sit-miR64 guides cleavage of Si008818m (65 reads, Fig. 3e) and Si035794m (188 reads, Fig. 3f). We found that the Si008818m and Si035794m proteins contain the same characteristic regions, namely QLQ (Gln, Leu and Gln) and WRC (Trp, Arg and Cys) domains, as A. thaliana growth-regulating factor proteins (AtGRFs). AtGRFs are involved in cell expansion in leaf and cotyledon tissues [69]. Examination of PARE-seq and RNA-seq reads mapping to the six miRNA target transcripts (Fig. 3) revealed a prominent cluster of reads at predicted cleavage locations in the PARE libraries, while no such pattern emerged in the RNA-seq library. These results, which reveal that miRNA-mediated degradation is the main pathway of mRNA degradation for some miRNA targets, were visualized using the Integrative Genomics Viewer [70].

Fig. 3

a-f Examples of confidently identified miRNA-directed cleavage. The complementary patterns of miRNA sequences and partial sequences of the target mRNAs are shown in the upper part of the figure and the numbers from parallel analysis of RNA end tags corresponding to cleavage sites are indicated by vertical arrowheads. “D” and “R” represent coverage of parallel analysis of RNA end tags (D) and RNA-sequencing reads (R) on selected miRNA targets. The mapped tags in “D” with the frequency at the position between bases 10 and 11 (from the miRNA 5') of the inset miRNA target alignment are indicated by red vertical arrowheads. Full details of all confidently identified miRNA targets are shown in Additional file 6. The sequences used for this figure came from the control sample. The number of reads mapping to each gene is indicated at the upper right

Insights gained into miRNA precursor metabolism

In plants, most primary miRNAs (pri-miRNAs) that are transcribed from miRNA genes by RNA polymerase II undergo two sequential cleavages by DCL1 to yield an RNA duplex containing the mature miRNA and miRNA* sequences [71, 72]. The locations of PARE tags mapped to the precursor can provide valuable hints to help reveal the details of the DCL1-guided two-step cleavage action on miRNA precursors [33].

To probe the patterns of miRNA hairpin processing in foxtail millet, the PARE reads in our datasets were mapped to 301 annotated [73] foxtail millet miRNA precursors (pre-miRNAs). Of the annotated pre-miRNAs, 114 (37.8 %) had one or more matching PARE tags (Additional file 8). The matching PARE tags were mainly mapped to precursors in one of two places: either the lower or upper cleavage site of the stem-loop 3′ arm. As shown in Fig. 4a, the matching PARE tags were mainly mapped to sit-miR166d, sit-miR166a-2, sit-miR167b-2 and sit-miR529a precursors at the lower cleavage site of the stem-loop 3′ arm; this indicates that these pre-miRNA hairpins could be processed by DCL1 via the classical loop-last mechanism [33] in which the first cleavage of pre-miRNA hairpins occurs at the loop-distal position (Fig. 4b). As shown in Fig. 4c, in contrast, matching PARE tags were mainly mapped to nov-sit-miR14, sit-miR156b-2, sit-miR319-1 and sit-miR535 precursors at the upper cleavage site of the stem-loop 3′ arm. This matching position implies that these pre-miRNAs are processed by DCL1 via an unusual loop-first mode [34, 74] in which the first cleavage of pre-miRNA hairpins occurs by precise processing at loop-proximal sites (Fig. 4d). It is noteworthy that this loop-first sequential processing of pre-miR319 family hairpins (sit-miR319-1 in Fig. 4a and sit-miR319-2 in Additional file 9) is also seen in A. thaliana, Physcomitrella and rice, and can generate two distinct miRNA/miRNA* duplexes [34, 74]. We also discovered that all miR156 family precursors identified in our study (sit-miR156b-2 in Fig. 4a and sit-miR156a-1, sit-miR156b-3 and sit-miR156d-2 in Additional file 9) are processed via the loop-first mechanism, whereas all identified miR166 family precursors (sit-miR166d and sit-miR166a-2 in Fig. 4b and sit-miR166a-1, sit-miR166a-3 and sit-miR166a-5 in Additional file 9) are associated with the loop-last mechanism.

Fig. 4

PARE tags mapping to foxtail millet miRNA hairpins. a and c Examples of “loop-last” (a) and “loop-first” (c) miRNA precursor processing. b and d A diagram of “loop-last” (b) and “loop-first” (d) processing. Regions within the pink and blue bars in (a) and (c) indicate the positions of the miRNA and miRNA* in the precursor, respectively. Two distinct miRNA/miRNA* duplexes were generated from sit-MIR319-1 and the two darker bars in sit-MIR319-1 indicate the miRNA/miRNA* duplexes of nov-sit-miR149. The read count at each position is indicated as a scatter plot

In pre-miRNA hairpin processing, DCL1-mediated cleavage occurs on each strand of the stem region [75]. The resulting 3′ cleavage products with poly (A) tails can be cloned by PARE high-throughput sequencing. If first-step cleavage occurs on both arms simultaneously, only cleavage signals mapped to the 3′ arm of pre-miRNA hairpins will be detected. If first-step cleavage occurs on both arms non-simultaneously, however, cleavage signals mapped to both the 3′ and 5′ arms of pre-miRNA hairpins will be identifiable in the PARE data. Our PARE data showed hardly any cleavage signals mapped to the stem-loop miRNA 5′ arm (Fig. 4 and Additional file 9), indicating that first-step cleavage occurs on both arms, most likely simultaneously, during the processing of most foxtail millet miRNA hairpins.

Distinct mRNA decay patterns among gene functional classes

As a diploid Panicoid crop species, foxtail millet uses C4 pathway photosynthesis. Although the leaves of C4 crops have increased nitrogen and water use efficiencies compared with C3 species [76], virtually nothing is known about how the pathway is regulated under drought stress. In C4 photosynthesis, carbon is shuttled as a C4 acid from the mesophyll to the bundle sheath cells to create a CO2 pump through a series of enzyme catalytic reactions. These enzymes include phosphoenolpyruvate carboxylase (PEPC), malate dehydrogenase (MDH), NADP-malic enzyme (ME) and pyruvate phosphate dikinase (PPDK). Using Phytozome, we found 32 C4 carbon shuttle enzyme (PEPC, MDH, PPDK and ME) genes in foxtail millet. To better understand the effect of drought on C4 enzymes, we analyzed the transcription and degradation of these 32 C4-related genes after PEG-induced drought stress (Fig. 5a). The transcript abundance of nine genes (four MDH genes, four ME genes and one PPDK gene) and the degradation abundance of eight genes (two MDH genes, four ME genes, one PPEC gene and one PPDK gene) showed significant changes after drought treatment. Among these drought-responsive genes, six displayed significant changes at both the transcription and degradation levels after drought treatment. These results suggest that various transcriptional and degradation regulatory mechanisms operate in C4-related genes under drought stress and may function to regulate water use efficiency in foxtail millet.

Fig. 5

Expression pattern of C4 photosynthetic-related genes and miRNA pathway-related regulators under drought stress. a and b Heatmaps showing the degradation level (RPM) and transcription level (RPKM) of genes encoding C4 photosynthesis pathway-related enzymes (a) and miRNA pathway-related regulators (b). Asterisks represent significant changes in transcription level after drought stress. Number signs represent significant changes in degradation level after drought stress. “c-R” and “d-R” represent transcription levels in the control (c-R) and drought-treated (d-R) samples revealed by RNA-seq. “c-D” and “d-D” represent degradation levels in the control (c-D) and drought-treated (d-D) samples revealed by PARE-seq. c The distribution of the ratio of relative uncapped mRNA abundance (RPM) versus total mRNA abundance (FPKM). c: control sample; d: drought-treated sample. “***” indicates statistically significant difference at P-value < 0.001 (Student’s Wilcox-test)

We also analyzed the effect of drought on relevant core regulators in the miRNA pathway, such as AGO, DCL and RNA-dependent RNA polymerase gene family members. Unlike C4-related genes, none of the 22 miRNA pathway-related regulator genes showed significant changes after drought treatment at the transcription or degradation levels (Fig. 5b). Interestingly, we found that the largest numbers of transcripts of these miRNA pathway-related regulators were enriched in uncapped forms (Fig. 5c). The distribution of relative uncapped to total mRNA ratios was found to be significantly biased (P < 0.001), which is consistent with results reported previously in Arabidopsis [15].


Different transcript degradation patterns were revealed during drought stress responses

To examine both gene synthesis and gene degradation, which were revealed respectively by RNA-seq and PARE-seq, genes were divided into four groups according to their change patterns after drought stress (Table 2).

Type-II genes (Table 2; classes B and H in Fig. 1c) could not be detected by analysis of RNA-seq data alone, as the amounts of synthesized mRNA were unchanged by the environmental stress conditions. GO analysis (Additional file 4) and detailed functional annotation (Additional file 5) revealed that many of the type-II genes belonged to diverse families of transcription factors such as WRKYs, MYBs, bZIPs and NACs. Many of these transcription factors play important roles in responses to drought stress [44, 45, 48, 77, 78]. Two significantly enriched motifs were identified in the 5′ UTRs of type-II gene mRNAs (Fig. 2f and g), implying the existence of some specific regulatory mechanism in this gene group. In contrast, the amounts of degraded type-IV gene (Table 2; classes D and F in Fig. 1c) mRNAs remained unchanged after drought treatment, whereas synthesized mRNA quantities showed obvious alterations. In the GO analysis, catalysis-related genes, such as hydrolases, isomerases, lyases, oxidoreductases and peroxidases (Additional file 4), were obviously over-represented among type-IV genes. Peroxidases are widely accepted as ‘stress enzymes’ [79]. Oxidoreductase can decompose H2O2 to water and molecular oxygen and is one of the key enzymes involved in the removal of toxic peroxides [80, 81]. Induction of the activity of these enzymes has been documented under a variety of stressful conditions, such as water stress [8284], chilling [85] and salinity [86], implying that these type-IV genes may serve as an intrinsic defense tool to resist drought stress in foxtail millet.

In type-I genes (Table 2; classes C and G in Fig. 1c), synthesized and degraded mRNA amounts followed similar trends after drought treatment. The GO enrichment analysis (Additional file 4) indicated that these genes were enriched in catalytic and various oxidation-related enzymes, suggesting that these types of genes may be important regulators of reactive oxygen species removal to maintain redox balance under drought stress conditions. The amounts of synthesized and degraded mRNAs showed opposite trends in type-III genes (Table 2; Classes A and I in Fig. 1c) following drought treatment. Unexpectedly, we detected only two type-III genes, for which changes in intact mRNA levels were enhanced by opposing changes in synthesis and degradation. However, in research on Brachypodium distachyon under cold stress, there were 1166 genes in type III (in which changes in transcript and uncapped transcript abundance showed opposite trends after cold treatment), but no obvious functional enrichment among these genes was found in GO analysis [16]. There are three possible reasons for this difference: (i) plants may have different regulation patterns in response to different abiotic stresses (cold and drought), (ii) different plants, e.g., B. distachyon (C3) and foxtail millet (C4), have different adaptation mechanisms in response to environmental stress, and (iii) plants have different regulatory mechanisms in response to different periods of abiotic stress (cold treatment for 24 h and drought treatment for 7 h). Perhaps a strong organismal response was not needed because of the short (7-h) duration of the drought treatment in our study.

miRNA targets

Using the PARE data in this study, we identified 385 putative miRNA-guided cleavages in foxtail millet (Additional file 7). Thus far, the transcripts of eight protein-coding targets of miRNA-mediated cleavage have been confirmed by gene-specific 5′ RACE in foxtail millet [73, 87]. Among these eight miRNA targets, seven (Si005991m, Si016509m, Si034525m, Si016508m, Si001804m, Si006975m and Si025305m) were found in our PARE data. A previous 5′ RACE analysis revealed that the Si016508m gene has two cleavage sites (at bp-positions 1085 and 1082) mediated by different miRNAs [73]. We also identified these two breakpoints in the Si016508m gene in our PARE analysis (Additional file 7). These results indicate that the miRNA-guided cleavages identified by the PARE analysis are genuine. A functional miRNA is expected to regulate target transcripts through two modes of action, either degradation or translation inhibition [13, 66, 67]. Because of the absence of detectable slicing, the PARE analysis was unable to find targets of miRNAs that act by repressing translation [30].

In a previous study, 43 known miRNAs and 212 novel miRNAs were identified in foxtail millet [73]. In our PARE analysis, we confirmed the targets of 80 % (34) of these known miRNAs, but only 34 % (73) of the novel ones (Additional file 7). Compared with known sit-miRNAs, nov-sit-miRNAs have been reported to have relatively lower expression levels and to exhibit higher tissue-specific expression [73]. Thus, we may have identified smaller numbers of targets of known and especially novel miRNAs because we did not analyze many developmental stages or different tissues.

Analysis of miRNA target expression after PEG-induced drought stress revealed that 50 miRNA targets were associated with different decay patterns during drought response (Additional file 10). The expression levels of some miRNA targets have been reported to be significantly altered after drought stress [50, 77, 88, 89]. For example, NAC genes targeted by miR164 negatively regulate drought resistance in rice [77], and nuclear factor Y, the target of miR169, is regulated transcriptionally and post-transcriptionally by miR169 to promote drought resistance in Arabidopsis [89]. In our study, the NAC protein genes Si017567m and Si010553m are also targets of miR164 and involved in response to drought stress. And Si037045m, a nuclear factor Y gene, is also targeted by miR169 and strongly induced by drought stress (Additional file 10). In addition, previous research revealed that the expression levels of sit-miR156, sit-miR160 and sit-miR397 were significantly altered after PEG-induced drought stress in foxtail millet [90]. Here, we found some of their targets, such as Si001804m (targeted by sit-miR156), Si016509m (targeted by sit-miR160) and Si001277m (targeted by sit-miR397), were associated with different decay patterns during the drought stress response (Additional file 10).

Besides their possible involvement in plant drought stress resistance, the miRNA targets identified in this study also play fundamental roles in plant growth and development. For example, GRFs targeted by nov-sit-miR64 have been shown to play an important role in cell expansion in leaf and cotyledon tissues in A. thaliana [69]. SPLs and AP2, targeted by miR156 and miR172, respectively, are responsible for the juvenile to adult transition in Arabidopsis [91].

miRNA precursor metabolism

The critical step of miRNA biogenesis is the precise processing of miRNA/miRNA* duplexes from precursor miRNA hairpins. Our PARE data suggest that both processing mechanisms exist in foxtail millet; some miRNA biogenesis was consistent with loop-last processing (Fig. 4a), whereas the precise processing of some miRNA precursors followed the unusual loop-first mode (Fig. 4b). We also observed that several different small RNAs originated from the same miRNA precursors. One such example is nov-sit-miR110 (conserved with miR159 in Arabidopsis) and nov-sit-miR155 (Additional file 9). This observation is consistent with prior research [92] demonstrating that the miR159 precursor can generate more than one mature miRNA in Arabidopsis. A similar example was found in the pre-miR319 family (Fig. 4 and Additional file 9), where two distinct miRNA/miRNA* duplexes were released from pre-miR319, both through loop-first processing (Fig. 4 and Additional file 9). Previous research has also demonstrated the occurrence of precise loop-first processing of an artificially generated miRNA with an Arabidopsis miR319 backbone [86]. One possible explanation for these observations is that one or more cis-acting hairpin features may direct DCL1 first-cleavage of miRNA hairpins at the loop-proximal position that is spatially separated from the loop-distal regions [34]. Identifying the features that direct both loop-last and loop-first modes of precise plant miRNA hairpin processing is an important goal for future research.

In addition to revealing sites corresponding to the remnants of the DCL-catalyzed cut, a previous study uncovered the existence of other cleavage signals in the middle of either miRNA- or miRNA*-coding loci on pre-miRNAs that are mediated by miRNA or miRNA* self-regulation [25]. In our study, we detected cleavage signals in the middle of miRNA- or miRNA*-coding loci in sit-miR535 (Fig. 4b), nov-sit-miR144-2 and nov-sit-miR120 (Additional file 9) precursors in foxtail millet. This finding is consistent with related observations in rice [28] and Arabidopsis [33] that imply that miRNA- or miRNA*-mediated self-regulation of certain miRNA precursors, although not widespread, exists in various organisms.

Some unexpected cleavage signals were also detected in several miRNA hairpins (Additional file 9C). For example, many different cleavage signals were detected in the region where mature nov-sit-miR121-1-5p and nov-sit-miR121-1-3p were located, while a few cleavage signals were found in other locations in the nov-sit-miR121-1 precursor. In addition, many different cleavage signals were detected in all nov-sit-miR35, nov-sit-miR82 and nov-sit-miR111 precursors. Similar unexpected results have also appeared in previous studies [28, 33] and could not be caused by methodological problems [28]. Although the currently understood model of miRNA biogenesis has been well studied [12, 28, 33, 34, 92], some novel processes not yet revealed must be responsible for these unexpected results and need further study.


In this study, mRNA synthesis and degradation, revealed respectively by RNA-seq and PARE data, were both taken into consideration in a global exploration of gene expression regulation under drought stress. This global analysis revealed new insights into gene expression under drought stress, confirmed many known regulatory mechanisms, and provided a window into many additional potentially novel pathways. Specific degradation patterns, such as miRNA mediation of target degradation and DCL1-mediated miRNA biogenesis were uncovered. Our results will not only deepen the understanding of mRNA degradation under stress conditions, but will also allow further insights into many targets of both known and novel miRNAs. Finally, these findings shed light on miRNA-precursor processing mechanisms in gramineous crops and biofuel grasses, which have close evolutionary relationships with foxtail millet.


Plant material and sequencing

Foxtail millet inbred line Yugu1 [35] was used in this study. Seeds of Yugu1 were germinated on a moist filter containing two layers of damp filter paper and incubated at 28 °C for 24 h. The germinated seedlings were then grown for 9 days in pots filled with a 1:1 mixture of nutrient soil: vermiculite under an illuminating incubator (28 °C day/20 °C night, 14-h photoperiod and 70 % relative humidity). The seedling roots were gently washed, transferred to 1/4 Hoagland’s solution, and allowed to grow for another 5 days. The aerial parts of seedlings (shoots), either subjected to drought stress (20 % PEG 6000 for 7 h) or untreated, were pooled and used for construction of the Dc and Dd libraries, respectively. The drought treatment was performed according to [4]. Total RNA was extracted from the 14-day-old Yugu1 seedlings using Trizol reagent (Invitrogen, USA) following the manufacturer’s protocol. PARE libraries were prepared according to [24], but with modifications described in [93]. Detailed descriptions of each step are provided in Additional file 1. The sequencing libraries were analyzed on a HiSeq 2000 sequencer (Illumina, USA) to produce 50-bp single reads. The generated raw reads have been deposited in NCBI’s SRA database under the accession number SRP061964.

PARE analysis and RNA-seq analysis

The foxtail millet genome sequence was downloaded from the Phytozome database (Sitalica_164_hardmasked.fa) [35, 94]. To analyze the PARE sequencing data, sequences corresponding to known tRNAs, rRNAs, small nucleolar RNAs, small nuclear RNAs and repeats/transposons were first removed [39, 40], and the remaining sequences were mapped to the foxtail millet genome using bowtie [95], with one mismatch allowed. Low-quality reads were excluded from subsequent analyses. PARE reads mapping to multiple positions were divided among the different positions [16]. An in-house developed Perl script was used to identify the total reads for specific transcripts. To compare specific transcripts, the mapped reads were enumerated and normalized against the total count of genome-mapped reads, reported as RPM, for each respective library. To assess the biological replicates, the log2 normalized data (RPM value + 1) were used to calculate correlation coefficients. Because the correlation between biological replicates was high (average R 2 = 0.97), we took the average RPM as the expression quantity. Genes with an adjusted P-value < 0.001, as determined by the DEGseq [96] software, and log2(fold change) ≥ 1 were considered up-regulated, while those with log2(fold change) ≤ −1 were considered down-regulated. Genes with −1 < log2(fold change) < 1 and a P-value > 0.001 were designated as unchanged.

To analyze the RNA-seq data, our previous analysis method was used [4]. After low quality reads were removed, Illumina sequencing reads were mapped to the foxtail millet genome using TopHat v2.0.4 [97], allowing two mismatches for the 100 bp reads. Then, the SAM files generated by TopHat were used as inputs for the Cufflinks software [97] and their expression levels (FPKM, fragments per kilobase of transcript per million fragments) were calculated with the foxtail millet genome annotation file. Genes with log2 (fold change) ≥ 1 and a P-value < 0.001 were considered up-regulated and genes with log2 (fold change) ≤ −1 and a P-value < 0.001 were considered down-regulated. Genes with −1 < log2 (fold change) < 1 and a P-value > 0.001 were designated as unchanged.

Gene ontology annotation

GO analysis was performed using WEGO [46]. The annotation frequency for each gene type was compared with corresponding annotations for the entire foxtail millet genome [35]. The Pearson chi-square test was used for statistical analysis. GO categories that showed significant (α = 0.05) enrichment were analyzed and displayed in the WEGO output histogram.

Analysis of sequence features

To analyze the features of mRNAs with different decay patterns, foxtail millet transcript sequences (release 164) and annotation data were downloaded from Phytozome [98]. The MFE value of the MFEI (MFE/(length of RNA sequence) × 100/(G + C)%) for secondary structures was calculated using RNAfold [99]. All other pattern searches and calculations were performed using custom scripts in Perl and R. Enriched motifs were identified using the MEME software [59] with E-value < 0.001. The Tomtom tool in MEME Suite ( was used to determine whether there were RNA-binding proteins in the analyzed gene classes.

miRNA target identification and miRNA precursor metabolism analysis

The classification of miRNA target categories was performed using the PAREsnip [100] pipeline based on PARE sequences (Dc library), foxtail millet transcripts (release 164) and miRNA sequences as inputs. PARE sequences (Dc library) matching miRNA hairpins according to BLASTN were assessed by comparisons of their respective genomic coordinates (Additional file 8). Foxtail millet miRNA and miRNA hairpin sequences were obtained from a previous publication [73]. Detailed information for the miRNAs discussed in this work is provided in Additional file 11.







Parallel analysis of RNA ends


RNA-induced silencing complex


Rapid amplification of cDNA ends


Dicer-like enzyme


The PARE libraries of drought treatment group


The PARE libraries of control group


Reads per million


Late embryogenesis abundant


Gene ontology


No Apical Meristem


Arginine decarboxylase


Minimal folding free energies index


Untranslated region


RNA-induced silencing complexes


Poly (rC)-binding protein 2


Serine/arginine-rich splicing factor 10


Auxin response factor 16


(Gln, Leu, Gln)


(Trp, Arg, Cys)


A. thaliana growth-regulating factor


Primary miRNA


miRNA precursor


Phosphoenolpyruvate carboxylase


Malate dehydrogenase


NADP-malic enzyme


Pyruvate phosphate dikinase




Dependent RNA polymerase


Nuclear factor Y


  1. 1.

    Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10(1):57–63.

  2. 2.

    Maskos U, Southern EM. Oligonucleotide hybridizations on glass supports: a novel linker for oligonucleotide synthesis and hybridization properties of oligonucleotides synthesised in situ. Nucleic Acids Res. 1992;20(7):1679–84.

  3. 3.

    Chen J, Zeng B, Zhang M, Xie S, Wang G, Hauck A, et al. Dynamic transcriptome landscape of maize embryo and endosperm development. Plant Physiol. 2014;166(1):252–64.

  4. 4.

    Qi X, Xie S, Liu Y, Yi F, Yu J. Genome-wide annotation of genes and noncoding RNAs of foxtail millet in response to simulated drought stress by deep sequencing. Plant Mol Biol. 2013;83(4–5):459–73.

  5. 5.

    Shim J, Karin M. The control of mRNA stability in response to extracellular stimuli. Mol Cells. 2002;14(3):323–31.

  6. 6.

    Mitchell P, Tollervey D. Musing on the structural organization of the exosome complex. Nat Struct Biol. 2000;7(10):843–6.

  7. 7.

    Mitchell P, Petfalski E, Houalla R, Podtelejnikov A, Mann M, Tollervey D. Rrp47p is an exosome-associated protein required for the 3' processing of stable RNAs. Mol Cell Biol. 2003;23(19):6982–92.

  8. 8.

    Coller J, Parker R. Eukaryotic mRNA decapping. Annu Rev Biochem. 2004;73:861–90.

  9. 9.

    Li Y, Kiledjian M. Regulation of mRNA decapping. Wiley Interdiscip Rev RNA. 2010;1(2):253–65.

  10. 10.

    Cohen LS, Mikhli C, Jiao X, Kiledjian M, Kunkel G, Davis RE. Dcp2 Decaps m2,2,7GpppN-capped RNAs, and its activity is sequence and context dependent. Mol Cell Biol. 2005;25(20):8779–91.

  11. 11.

    Chen X. MicroRNA metabolism in plants. Curr Top Microbiol Immunol. 2008;320:117–36.

  12. 12.

    Ramachandran V, Chen X. Small RNA metabolism in Arabidopsis. Trends Plant Sci. 2008;13(7):368–74.

  13. 13.

    Brodersen P, Sakvarelidze-Achard L, Bruun-Rasmussen M, Dunoyer P, Yamamoto YY, Sieburth L, et al. Widespread translational inhibition by plant miRNAs and siRNAs. Science. 2008;320(5880):1185–90.

  14. 14.

    Gutierrez RA, Ewing RM, Cherry JM, Green PJ. Identification of unstable transcripts in Arabidopsis by cDNA microarray analysis: rapid decay is associated with a group of touch- and specific clock-controlled genes. Proc Natl Acad Sci U S A. 2002;99(17):11513–8.

  15. 15.

    Jiao Y, Riechmann JL, Meyerowitz EM. Transcriptome-wide analysis of uncapped mRNAs in Arabidopsis reveals regulation of mRNA degradation. Plant Cell. 2008;20(10):2571–85.

  16. 16.

    Zhang J, Mao Z, Chong K. A global profiling of uncapped mRNAs under cold stress reveals specific decay patterns and endonucleolytic cleavages in Brachypodium distachyon. Genome Biol. 2013;14(8):R92.

  17. 17.

    Zhang X, Zhu Y, Liu X, Hong X, Xu Y, Zhu P, et al. Plant biology. Suppression of endogenous gene silencing by bidirectional cytoplasmic RNA decay in Arabidopsis. Science. 2015;348(6230):120–3.

  18. 18.

    Valencia-Sanchez MA, Liu J, Hannon GJ, Parker R. Control of translation and mRNA degradation by miRNAs and siRNAs. Genes Dev. 2006;20(5):515–24.

  19. 19.

    Haley B, Zamore PD. Kinetic analysis of the RNAi enzyme complex. Nat Struct Mol Biol. 2004;11(7):599–606.

  20. 20.

    Martinez J, Tuschl T. RISC is a 5′ phosphomonoester-producing RNA endonuclease. Genes Dev. 2004;18(9):975–80.

  21. 21.

    Wang XJ, Reyes JL, Chua NH, Gaasterland T. Prediction and identification of Arabidopsis thaliana microRNAs and their mRNA targets. Genome Biol. 2004;5(9):R65.

  22. 22.

    Rhoades MW, Reinhart BJ, Lim LP, Burge CB, Bartel B, Bartel DP. Prediction of plant microRNA targets. Cell. 2002;110(4):513–20.

  23. 23.

    Yekta S, Shih IH, Bartel DP. MicroRNA-directed cleavage of HOXB8 mRNA. Science. 2004;304(5670):594–6.

  24. 24.

    German MA, Luo S, Schroth G, Meyers BC, Green PJ. Construction of Parallel Analysis of RNA Ends (PARE) libraries for the study of cleaved miRNA targets and the RNA degradome. Nat Protoc. 2009;4(3):356–62.

  25. 25.

    German MA, Pillay M, Jeong DH, Hetawal A, Luo S, Janardhanan P, et al. Global identification of microRNA-target RNA pairs by parallel analysis of RNA ends. Nat Biotechnol. 2008;26(8):941–6.

  26. 26.

    Addo-Quaye C, Eshoo TW, Bartel DP, Axtell MJ. Endogenous siRNA and miRNA targets identified by sequencing of the Arabidopsis degradome. Curr Biol. 2008;18(10):758–62.

  27. 27.

    Wu L, Zhang Q, Zhou H, Ni F, Wu X, Qi Y. Rice MicroRNA effector complexes and targets. Plant Cell. 2009;21(11):3421–35.

  28. 28.

    Li YF, Zheng Y, Addo-Quaye C, Zhang L, Saini A, Jagadeeswaran G, et al. Transcriptome-wide identification of microRNA targets in rice. Plant J. 2010;62(5):742–59.

  29. 29.

    Pantaleo V, Szittya G, Moxon S, Miozzi L, Moulton V, Dalmay T, et al. Identification of grapevine microRNAs and their targets using high-throughput sequencing and degradome analysis. Plant J. 2010;62(6):960–76.

  30. 30.

    Voinnet O. Origin, biogenesis, and activity of plant microRNAs. Cell. 2009;136(4):669–87.

  31. 31.

    Kim VN, Han J, Siomi MC. Biogenesis of small RNAs in animals. Nat Rev Mol Cell Biol. 2009;10(2):126–39.

  32. 32.

    Ha M, Kim VN. Regulation of microRNA biogenesis. Nat Rev Mol Cell Biol. 2014;15(8):509–24.

  33. 33.

    Meng Y, Gou L, Chen D, Wu P, Chen M. High-throughput degradome sequencing can be used to gain insights into microRNA precursor metabolism. J Exp Bot. 2010;61(14):3833–7.

  34. 34.

    Addo-Quaye C, Snyder JA, Park YB, Li YF, Sunkar R, Axtell MJ. Sliced microRNA targets and precise loop-first processing of MIR319 hairpins revealed by analysis of the Physcomitrella patens degradome. RNA. 2009;15(12):2112–21.

  35. 35.

    Bennetzen JL, Schmutz J, Wang H, Percifield R, Hawkins J, Pontaroli AC, et al. Reference genome sequence of the model plant Setaria. Nat Biotechnol. 2012;30(6):555–61.

  36. 36.

    Doust AN, Kellogg EA, Devos KM, Bennetzen JL. Foxtail millet: a sequence-driven grass model system. Plant Physiol. 2009;149(1):137–41.

  37. 37.

    Zhang G, Liu X, Quan Z, Cheng S, Xu X, Pan S, et al. Genome sequence of foxtail millet (Setaria italica) provides insights into grass evolution and biofuel potential. Nat Biotechnol. 2012;30(6):549–54.

  38. 38.

    Muthamilarasan M, Prasad M. Advances in Setaria genomics for genetic improvement of cereals and bioenergy grasses. Theor Appl Genet. 2015;128(1):1–14.

  39. 39.

    Repbase Update. Access date: 26 April 2014..

  40. 40.

    The Rfam database. Access date: 26 April 2014.

  41. 41.

    Xiao B, Huang Y, Tang N, Xiong L. Over-expression of a LEA gene in rice improves drought resistance under the field conditions. Theor Appl Genet. 2007;115(1):35–46.

  42. 42.

    Colmenero-Flores JM, Campos F, Garciarrubio A, Covarrubias AA. Characterization of Phaseolus vulgaris cDNA clones responsive to water deficit: identification of a novel late embryogenesis abundant-like protein. Plant Mol Biol. 1997;35(4):393–405.

  43. 43.

    Xu D, Duan X, Wang B, Hong B, Ho T, Wu R. Expression of a late embryogenesis abundant protein gene, HVA1, from barley confers tolerance to water deficit and salt stress in transgenic rice. Plant Physiol. 1996;110(1):249–57.

  44. 44.

    Abe H, Yamaguchi-Shinozaki K, Urao T, Iwasaki T, Hosokawa D, Shinozaki K. Role of arabidopsis MYC and MYB homologs in drought- and abscisic acid-regulated gene expression. Plant Cell. 1997;9(10):1859–68.

  45. 45.

    Ren X, Chen Z, Liu Y, Zhang H, Zhang M, Liu Q, et al. ABO3, a WRKY transcription factor, mediates plant responses to abscisic acid and drought tolerance in Arabidopsis. Plant J. 2010;63(3):417–29.

  46. 46.

    Ye J, Fang L, Zheng H, Zhang Y, Chen J, Zhang Z, et al. WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 2006;34(Web Server issue):W293–7.

  47. 47.

    Rajan VB, D’Silva P. Arabidopsis thaliana J-class heat shock proteins: cellular stress sensors. Funct Integr Genomics. 2009;9(4):433–46.

  48. 48.

    Chen T, Li W, Hu X, Guo J, Liu A, Zhang B. A cotton MYB transcription factor, GbMYB5, is positively involved in plant adaptive response to drought stress. Plant Cell Physiol. 2015;56(5):917–29.

  49. 49.

    Urano K, Yoshiba Y, Nanjo T, Ito T, Yamaguchi-Shinozaki K, Shinozaki K. Arabidopsis stress-inducible gene for arginine decarboxylase AtADC2 is required for accumulation of putrescine in salt tolerance. Biochem Biophys Res Commun. 2004;313(2):369–75.

  50. 50.

    Rabbani MA, Maruyama K, Abe H, Khan MA, Katsura K, Ito Y, et al. Monitoring expression profiles of rice genes under cold, drought, and high-salinity stresses and abscisic acid application using cDNA microarray and RNA gel-blot analyses. Plant Physiol. 2003;133(4):1755–67.

  51. 51.

    Alcazar R, Planas J, Saxena T, Zarza X, Bortolotti C, Cuevas J, et al. Putrescine accumulation confers drought tolerance in transgenic Arabidopsis plants over-expressing the homologous Arginine decarboxylase 2 gene. Plant Physiol Biochem. 2010;48(7):547–52.

  52. 52.

    Roy M, Wu R. Arginine decarboxylase transgene expression and analysis of environmental stress tolerance in transgenic rice. Plant Sci. 2001;160(5):869–75.

  53. 53.

    Sharkey TD, Seemann JR. Mild water stress effects on carbon-reduction-cycle intermediates, ribulose bisphosphate carboxylase activity, and spatial homogeneity of photosynthesis in intact leaves. Plant Physiol. 1989;89(4):1060–5.

  54. 54.

    Zhang BH, Pan XP, Cox SB, Cobb GP, Anderson TA. Evidence that miRNAs are different from other RNAs. Cell Mol Life Sci. 2006;63(2):246–54.

  55. 55.

    Zhao C, Hamilton T. Introns regulate the rate of unstable mRNA decay. J Biol Chem. 2007;282(28):20230–7.

  56. 56.

    Wang HF, Feng L, Niu DK. Relationship between mRNA stability and intron presence. Biochem Biophys Res Commun. 2007;354(1):203–8.

  57. 57.

    Masaru OT. Control of mRNA stability in plants--cis-acting elements which affect mRNA stability. Tanpakushitsu Kakusan Koso. 1993;38(14):2367–71.

  58. 58.

    Bailey TL, Boden M, Buske FA, Frith M, Grant CE, Clementi L, et al. MEME SUITE: tools for motif discovery and searching. Nucleic Acids Res. 2009;37(Web Server issue):W202–8.

  59. 59.

    Bailey TL, Williams N, Misleh C, Li WW. MEME: discovering and analyzing DNA and protein sequence motifs. Nucleic Acids Res. 2006;34(Web Server issue):W369–73.

  60. 60.

    Ray D, Kazan H, Cook KB, Weirauch MT, Najafabadi HS, Li X, et al. A compendium of RNA-binding motifs for decoding gene regulation. Nature. 2013;499(7457):172–7.

  61. 61.

    Fujimura K, Kano F, Murata M. Identification of PCBP2, a facilitator of IRES-mediated translation, as a novel constituent of stress granules and processing bodies. RNA. 2008;14(3):425–31.

  62. 62.

    You F, Sun H, Zhou X, Sun W, Liang S, Zhai Z, et al. PCBP2 mediates degradation of the adaptor MAVS via the HECT ubiquitin ligase AIP4. Nat Immunol. 2009;10(12):1300–8.

  63. 63.

    Kalsotra A, Cooper TA. Functional consequences of developmentally regulated alternative splicing. Nat Rev Genet. 2011;12(10):715–29.

  64. 64.

    Zhou X, Wu W, Li H, Cheng Y, Wei N, Zong J, et al. Transcriptome analysis of alternative splicing events regulated by SRSF10 reveals position-dependent splicing modulation. Nucleic Acids Res. 2014;42(6):4019–30.

  65. 65.

    Llave C, Xie Z, Kasschau KD, Carrington JC. Cleavage of Scarecrow-like mRNA targets directed by a class of Arabidopsis miRNA. Science. 2002;297(5589):2053–6.

  66. 66.

    Jones-Rhoades MW, Bartel DP, Bartel B. MicroRNAS and their regulatory roles in plants. Annu Rev Plant Biol. 2006;57:19–53.

  67. 67.

    Li S, Liu L, Zhuang X, Yu Y, Liu X, Cui X, et al. MicroRNAs inhibit the translation of target mRNAs on the endoplasmic reticulum in Arabidopsis. Cell. 2013;153(3):562–74.

  68. 68.

    Wang JW, Wang LJ, Mao YB, Cai WJ, Xue HW, Chen XY. Control of root cap formation by MicroRNA-targeted auxin response factors in Arabidopsis. Plant Cell. 2005;17(8):2204–16.

  69. 69.

    Kim JH, Choi D, Kende H. The AtGRF family of putative transcription factors is involved in leaf and cotyledon growth in Arabidopsis. Plant J. 2003;36(1):94–104.

  70. 70.

    Thorvaldsdottir H, Robinson JT, Mesirov JP. Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform. 2013;14(2):178–92.

  71. 71.

    Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116(2):281–97.

  72. 72.

    Bologna NG, Schapire AL, Palatnik JF. Processing of plant microRNA precursors. Brief Funct Genomics. 2013;12(1):37–45.

  73. 73.

    Yi F, Xie S, Liu Y, Qi X, Yu J. Genome-wide characterization of microRNA in foxtail millet (Setaria italica). BMC Plant Biol. 2013;13:212.

  74. 74.

    Bologna NG, Mateos JL, Bresso EG, Palatnik JF. A loop-to-base processing mechanism underlies the biogenesis of plant microRNAs miR319 and miR159. EMBO J. 2009;28(23):3646–56.

  75. 75.

    Ji X. The mechanism of RNase III action: how dicer dices. Curr Top Microbiol Immunol. 2008;320:99–116.

  76. 76.

    Aubry S, Brown NJ, Hibberd JM. The role of proteins in C(3) plants prior to their recruitment into the C(4) pathway. J Exp Bot. 2011;62(9):3049–59.

  77. 77.

    Fang Y, Xie K, Xiong L. Conserved miR164-targeted NAC genes negatively regulate drought resistance in rice. J Exp Bot. 2014;65(8):2119–35.

  78. 78.

    Yang X, Wang X, Ji L, Yi Z, Fu C, Ran J, et al. Overexpression of a Miscanthus lutarioriparius NAC gene MlNAC5 confers enhanced drought and cold tolerance in Arabidopsis. Plant Cell Rep. 2015;34(6):943–58.

  79. 79.

    Passaia G, Margis-Pinheiro M. Glutathione peroxidases as redox sensor proteins in plant cells. Plant Sci. 2015;234:22–6.

  80. 80.

    Man H, Wells E, Hussain S, Leipold F, Hart S, Turkenburg JP, et al. Structure, Activity and Stereoselectivity of NADPH-Dependent Oxidoreductases Catalysing the S-Selective Reduction of the Imine Substrate 2-Methylpyrroline. Chembiochem. 2015;16(7):1052–9.

  81. 81.

    Verma S, Dubey RS. Lead toxicity induces lipid peroxidation and alters the activities of antioxidant enzymes in growing rice plants. Plant Sci. 2003;164(4):645–55.

  82. 82.

    Sharma P, Dubey RS. Drought induces oxidative stress and enhances the activities of antioxidant enzymes in growing rice seedlings. Plant Growth Regul. 2005;46(3):209–21.

  83. 83.

    Khanna SM, Taxak PC, Jain PK, Saini R, Srinivasan R. Glycolytic enzyme activities and gene expression in Cicer arietinum exposed to water-deficit stress. Appl Biochem Biotechnol. 2014;173(8):2241–53.

  84. 84.

    Zhang F, Zhong H, Han X, Guo Z, Yang W, Liu Y, et al. Proteomic profile of Aspergillus flavus in response to water activity. Fungal Biol. 2015;119(2–3):114–24.

  85. 85.

    Xu J, Yang J, Duan X, Jiang Y, Zhang P. Increased expression of native cytosolic Cu/Zn superoxide dismutase and ascorbate peroxidase improves tolerance to oxidative and chilling stresses in cassava (Manihot esculenta Crantz). BMC Plant Biol. 2014;14:208.

  86. 86.

    Ranjit SL, Manish P, Penna S. Early osmotic, antioxidant, ionic, and redox responses to salinity in leaves and roots of Indian mustard (Brassica juncea L.). Protoplasma, 2015 [Epub ahead of print].

  87. 87.

    Han J, Xie H, Sun Q, Wang J, Lu M, Wang W, et al. Bioinformatic identification and experimental validation of miRNAs from foxtail millet (Setaria italica). Gene. 2014;546(2):367–77.

  88. 88.

    Mallory AC, Vaucheret H. Functions of microRNAs and related small RNAs in plants. Nat Genet. 2006;38(Suppl):S31–6.

  89. 89.

    Li WX, Oono Y, Zhu J, He XJ, Wu JM, Iida K, et al. The Arabidopsis NFYA5 transcription factor is regulated transcriptionally and posttranscriptionally to promote drought resistance. Plant Cell. 2008;20(8):2238–51.

  90. 90.

    Khan Y, Yadav A, Bonthala VS, Muthamilarasan M, Yadav CB, Prasad M. Comprehensive genome-wide identification and expression profiling of foxtail millet [Setaria italica (L.)] miRNAs in response to abiotic stress and development of miRNA database. Plant Cell Tiss Organ Cult. 2014;118(2):279–92.

  91. 91.

    Wu G, Park MY, Conway SR, Wang JW, Weigel D, Poethig RS. The sequential action of miR156 and miR172 regulates developmental timing in Arabidopsis. Cell. 2009;138(4):750–9.

  92. 92.

    Zhang W, Gao S, Zhou X, Xia J, Chellappan P, Zhou X, et al. Multiple distinct small RNAs originate from the same microRNA precursors. Genome Biol. 2010;11(8):R81.

  93. 93.

    Ma Z, Coruh C, Axtell MJ. Arabidopsis lyrata small RNAs: transient MIRNA and small interfering RNA loci within the Arabidopsis genus. Plant Cell. 2010;22(4):1090–103.

  94. 94.

    Goodstein DM, Shu S, Howson R, Neupane R, Hayes RD, Fazo J, et al. Phytozome: a comparative platform for green plant genomics. Nucleic Acids Res. 2012;40(Database issue):D1178–86.

  95. 95.

    Langmead B. Aligning short sequencing reads with Bowtie. Curr Protoc Bioinformatics, 2010. Chapter 11: p. Unit 11.7. doi:10.1002/0471250953.bi1107s32.

  96. 96.

    Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26(1):136–8.

  97. 97.

    Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7(3):562–78.

  98. 98.

    Phytozome v9.1: Setaria italica.

  99. 99.

    Gruber AR, Lorenz R, Bernhart SH, Neuböck R, Hofacker IL. The Vienna RNA websuite. Nucleic Acids Res. 2008;36(Web Server issue):W70–4. Access date: 10 October 2014.

  100. 100.

    Folkes L, Moxon S, Woolfenden HC, Stocks MB, Szittya G, Dalmay T, et al. PAREsnip: a tool for rapid genome-wide discovery of small RNA/target interactions evidenced through degradome sequencing. Nucleic Acids Res. 2012;40(13):e103.

Download references


We thank Xin Qi for her help in material preparation and Yuwei Liu for useful suggestions on data analysis. This work was supported by the National Basic Research Program of China (Grant No. 2012CB215301) and the National Transgenic Major Program of China (Grant No. 2014ZX08003-002; 2013ZX08003-002).

Author information

Correspondence to Jingjuan Yu.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

JY and FY conceived and designed the research. FY prepared the samples, and wrote the original manuscript. FY and JC performed the analysis. JY revised thoroughly the manuscript and finalize the manuscript. All authors read and approved the final manuscript.

Authors’ information

Not applicable.

Availability of data and materials

Not applicable.

Additional files

Additional file 1.

Overview of the PARE-seq method for the isolation of uncapped mRNA. (PDF 80 kb)

Additional file 2.

Summary of PARE analysis. (XLSX 9 kb)

Additional file 3.

Correlation between PARE biological replicates. We calculated the expression level (RPM) for each replicate separately and compared them to one another. The normalized data of log2 (RPM value + 1) was used to calculate the correlation coefficient and the correlation between the biological replicates was high (average R 2 = 0.97). Dd: the PARE libraries of drought treatment group, Dc: the PARE libraries of control group. (PDF 65 kb)

Additional file 4.

GO functional enrichment analysis for different decay pattern mRNA. Gene Ontology (GO) analysis was performed for the type I, II, IV genes using WEGO which organizes information for cellular component categories, molecular function and biological process. The Pearson Chi-square test was used for statistical analysis. GO categories that show a significant (α =0.05) enrichment were analyzed and displayed here. (PDF 396 kb)

Additional file 5.

Fold change and function annotations of four types of genes with different decay. Annotations were retrieved from phytozome. Note: Columns are blank if no corresponding data is available. (XLSX 252 kb)

Additional file 6.

Gene transcript features to different mRNA decay patterns. Analyzed GC content and MFEI value of mRNA, 5' UTR and 3' UTR for type-I, −II, −IV, randomly selected genes and all genes in foxtail millet. “R” represented randomly selected genes, “A” represented all genes, “***” indicate statistically significant differences at P value < 0.001 (Student’s wilcox-test). (PDF 257 kb)

Additional file 7.

miRNA/target interactions identified in the degradome library by PAREsnip and function annotations of targets. Annotations were retrieved from phytozome. Note: Columns are blank if no corresponding data is available. (XLSX 40 kb)

Additional file 8.

PARE signatures mapping to foxtail millet miRNA hairpins. (XLSX 34 kb)

Additional file 9.

Foxtail millet PARE tags matching miRNA hairpins. (A) “Possible loop-last” processing of miRNA hairpins. (B) “Loop-first” processing of miRNA hairpins. (C) Unexpected cleavage signals in several miRNA hairpins. Regions within the pink and blue bars indicate the positions of the miRNA and miRNA* in the precursor, respectively. The read count at each position is indicated as a scatter plot. (PDF 356 kb)

Additional file 10.

Fold change and function annotations of miRNA targets which are response to drought. Annotations were retrieved from phytozome. Note: Columns are blank if no corresponding data is available. (XLSX 14 kb)

Additional file 11.

The detail information about miRNAs which have discussed in this work. (XLSX 14 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Yi, F., Chen, J. & Yu, J. Global analysis of uncapped mRNA changes under drought stress and microRNA-dependent endonucleolytic cleavages in foxtail millet. BMC Plant Biol 15, 241 (2015) doi:10.1186/s12870-015-0632-0

Download citation


  • Uncapped mRNA
  • Drought stress
  • PARE
  • Foxtail millet
  • miRNA target
  • miRNA precursor
  • C4