Skip to main content
Download PDF
  • Research article
  • Open access
  • Published:

Temporal responses of conserved miRNAs to drought and their associations with drought tolerance and productivity in rice

BMC Genomics volume 21, Article number: 232 (2020) Cite this article

Abstract

Background

Plant miRNAs play crucial roles in responses to drought and developmental processes. It is essential to understand the association of miRNAs with drought-tolerance (DT), as well as their impacts on growth, development, and reproduction (GDP). This will facilitate our utilization of rice miRNAs in breeding.

Results

In this study, we investigated the time course of miRNA responses to a long-term drought among six rice genotypes by high-throughput sequencing. In total, 354 conserved miRNAs were drought responsive, representing obvious genotype- and stage-dependent patterns. The drought-responsive miRNAs (DRMs) formed complex regulatory network via their coexpression and direct/indirect impacts on the rice transcriptome. Based on correlation analyses, 211 DRMs were predicted to be associated with DT and/or GDP. Noticeably, 14.2% DRMs were inversely correlated with DT and GDP. In addition, 9 pairs of mature miRNAs, each derived from the same pre-miRNAs, were predicted to have opposite roles in regulating DT and GDP. This suggests a potential yield penalty if an inappropriate miRNA/pre-miRNA is utilized. miRNAs have profound impacts on the rice transcriptome reflected by great number of correlated drought-responsive genes. By regulating these genes, a miRNA could activate diverse biological processes and metabolic pathways to adapt to drought and have an influence on its GDP.

Conclusion

Based on the temporal pattern of miRNAs in response to drought, we have described the complex network between DRMs. Potential associations of DRMs with DT and/or GDP were disclosed. This knowledge provides valuable information for a better understanding in the roles of miRNAs play in rice DT and/or GDP, which can facilitate our utilization of miRNA in breeding.

Background

MicroRNAs (miRNAs) are a large class of small noncoding RNAs of 20 to 24 nucleotides (nt) in length [36, 43, 44]. The miRNA and its target mRNA can form the miRNA-induced silencing RISC complex, which inhibits the protein of its target genes by either destabilizing the mRNA or by inhibiting its translation [43, 44]. The RISC complex negatively regulates gene expression at the posttranscriptional level. miRNA target transcription factors, many of which are critical regulators in plant growth, development, and reproduction (GDP), and stress responses [36, 38, 49]. Therefore, a miRNA that has great impacts on the transcriptome, is located at the center of complex gene regulatory networks associated with plant GDP and stress-tolerance, [7, 38]. The ability of plants to employ miRNAs to posttranscriptionally inactive or induce the expression of stress-responsive genes provides an advantage compared with regulation by transcription factors alone [49]. It makes miRNAs good targets to improve crop stress-tolerance [7, 38, 49]. However, most miRNAs do not work independently in response to environmental stresses. Their stress responses are tightly coordinated with multiple developmental processes via the complex regulatory network [7, 38] or multihormone responses [29]. Increasing evidence has shown that a miRNA involved in stress tolerance commonly exerts pleiotropic effects on the GDP of plants [38]. It means we should avoid potential negative effects on productivity when developing tolerant cultivars by genetically modifying miRNAs. This requires improved understanding of the association of a miRNA with stress-tolerance and/or GDP.

Drought is a major limiting environmental factor for crops and causes great loss in yield annually. It is essential to develop drought-tolerant crops for food security [12]. Recently, attentions have been focused on the importance of posttranscriptional regulation by miRNAs in drought tolerance (DT) due to their central roles in the regulatory network [7, 38]. With the fast development of next-generation sequencing, drought-responsive miRNAs (DRMs) have been identified in diverse crops, including cotton [48], rapeseed [21], maize [1], tomato [28], and rice (Oryza sativa) [3, 5, 6]. Many DRMs have been characterized as important modulators in DT via regulating the expression of drought-responsive genes [7]. Most miRNAs are induced by drought and downregulate their target mRNAs [7], which are negative factors in the drought response [8, 53]. Conversely, some other miRNAs are downregulated by drought [7], leading to the accumulation of target mRNAs positively contributing to drought adaptation [25, 37].

Rice is one of the most important cereal food for more than half of the global population. Unfortunately, elite rice is very sensitive to drought due to its long-term domestication in irrigated fields [4, 45]. The improvement of DT in rice is thus a primary breeding aim for “green super rice” [31, 52]. For this purpose, the roles played by miRNAs in rice drought-resistance have been widely studied. There have been 604 pre-miRNAs, which encode 738 mature miRNAs, identified in rice and recorded in miRBase (release 22.1, [22]). Hundreds of miRNAs have been determined as DRMs by several genome-wide investigations in different genotypes or tissues [2, 3, 6, 56]. However, there is still a large knowledge gap between the identification of DRMs and the characterization of their associations with DT [17]. According to large number of recommended DRMs, only very low proportions of DRMs have been functionally proven in rice, including miRNA162 [40], miRNA164 [8], miRNA166 [51], miRNA393 [46], and miRNA408 [37]. Among these drought-tolerant miRNAs, miRNA408 [37, 50] and miRNA393 [46] have been reported to have unwanted pleiotropic effects on GDP. For better utilization of miRNAs, it is necessary to understand their associations with DT and/or GDP in rice.

Many former studies have typically investigated a single genotype [3, 6, 56] or two rice genotypes of contrasting DT [5] to identify DRMs. A miRNA that is differentially regulated in response to environmental stress is not necessarily associated with stress tolerance [17]. Therefore, it is essential to study diverse genotypes, which allows us to eliminate bias caused by a limited number of genotypes. Rice adaptation to drought is a progressive process with sequential molecular, physiological, and morphological alterations [9, 35]. However, the time course of miRNA expression and regulations in rice under drought has not been fully understood, but from it we can learn potential associations between miRNAs and physiological/ morphological responses [17]. To understand the potential roles played by miRNAs in rice DT, we investigated the genome-wide expression of miRNAs in six rice genotypes at five time points under drought stress and one time point at recovery. Meanwhile, we also investigate the transcriptomes of six genotypes by RNA-sequencing, from which we can learn the potential impacts of miRNAs on the rice transcriptome. The design of our experiment allows us to address the following questions: (1) How are miRNAs sequentially regulated in response to progressive drought? (2) Do any DRMs associate with drought-tolerance and/or GDP? (3) Which miRNAs are good candidates for improving rice DT? This knowledge can advance our utilization of miRNAs to improve DT without yield penalty in rice.

Results

Alterations of morphological and physiological traits among rice genotypes under drought conditions

The growth, development, and productivity of six rice genotypes were greatly affected by drought, as reflected in reduced plant height, number of seeds per plant, seed weight per plant, and biomass, and delayed heading date (Fig. 1, Additional file: Table S1). Drought also caused the accumulation of H2O2 content (Additional file: Table S2) and dead leaves (Additional file: Table S1) in the rice genotypes. To resist drought, the rice activated mechanisms of osmotic adjustment and ROS scavenging, as reflected in the largely increased osmotic potential (Additional file: Table S3) and total antioxidant capacity (Additional file 1: Table S4) under drought conditions, particularly in later drought time points (D3-D5).

Fig. 1
figure 1

Relative performances (performance under drought (DT) /that under well-watered (CK)) of six rice genotypes. * indicates significant differences between traits measured in DT and those measured in CK

Full size image

Sequence analysis of small RNAs in sequenced samples

A total of ~ 1.068G raw reads were obtained from 66 samples (libraries). After the removal of low-quality reads, adapters, reads shorter than 18 nt, and other contaminating sequences, 792.8 M clean reads (74.3%) were finally retained, including 172.4 M unique reads (Additional file 1: Table S5). Among total clean reads between 18 and 32 nt, 39.5% reads were matched to miRNA (~ 21%), tRNA (~ 5%), and rRNA (~ 12%), respectively (Additional file 2: Figure S1). The distribution of reads in various sizes of small RNAs was not homogeneous. The most abundant were small RNAs of 21 nt (27.4%) and 24 nt (20.1%) in length (Additional file 2: Figure S2). We should also point out that proportions of miRNAs of 21 nt and 24 nt in length had great variations among genotypes, time points, and treatments (Additional file 2: Figure S2).

General description of drought-responsive and recovery-related miRNAs detected in the six rice genotypes

A total of 632 conserved mature miRNAs in miRBase were detected in 66 sequenced samples (Additional file 1: Table S6). Among the expressed miRNAs, 549 miRNAs were available for further analysis (TPM > 0.1 in at least one sample) (Additional file 1: Table S6). During the drought period, 354 miRNAs in 57 families were identified as drought-responsive miRNAs (DRMs). Moreover, 80 differentially expressed miRNAs were detected at the recovery stage and were determined to be recovery-related miRNAs (RRMs) (Additional file 1: Table S6). A considerable proportion (48.6%, 172 out of 354) of DRMs were regulated in a genotype-specific (Additional file 2: Figure S3) or temporal-specific (Additional file 2: Figure S4) manner. There were 78–239 DRMs and 77–116 RRMs among the six genotypes (Table S6). A susceptible genotype S18 (239) and a tolerant genotype S11 (216) had the most DRMs. Meanwhile, a susceptible genotype S24 (78) and a tolerant genotype S28 (87) had the least DRMs. This result indicated the number of DRMs should be not related with rice drought tolerance. However, 107 DRMs could be frequently (frequency ≥ 3) detected among different genotypes and time points (Additional file 2: Figure S5a), suggesting that they have universal roles in rice adaptation to drought. We also detected great variance in number of RRMs among the six genotypes. Interestingly, the three tolerant genotypes S3, S11, and S28 possessed more RRMs (from 6 to 66), while the susceptible ones had less RRMs. Finally, we detected no recovery-specific differentially expressed miRNAs (Additional file 2: Figure S5b). In addition, regulation patterns of most characterized miRNAs (e.g. miR160, miR162, miR393, miR397, and miR408) were consistent with previous studies (Additional file 1: Table S6). However, regulations of some other characterized miRNAs (e.g. miR166, miR172, and miR396) represented great variation among genotypes (Additional file 1: Table S6).

Correlations of expressions among DRMs

Coexpression relationships between DRMs were revealed by their positive or negative correlations (Fig. 2). Pearson correlation coefficients (PCCs) (0.620 ± 0.013, p < 0.001) between miRNAs of the same family (e.g., miRNA169, miRNA395, miRNA818) or PCCs between pairs of miRNAs derived from the same pre-miRNAs (0.453 ± 0.059, p < 0.001) (e.g., miRNA1320-3p/5p, miRNA528-3p/5p, miRNA7695-3p/5p) were significantly higher than the average PCC (0.084 ± 0.007) by both Mann-Whitney and Kolmogorov-Smirnov tests. High PPC values could also be frequently detected between some unrelated miRNAs (e.g. miRNA1862 with miRNA169 and miRNA869, miRNA156 with miRNA169 and miRNA815) (Fig. 2). Additionally, we detected many negatively correlated miRNAs (e.g., miRNA818 with miRNA169 and miRNA166, miRNA395 with miRNA169). These results indicated complicated regulatory networks of miRNAs in response to drought.

Fig. 2
figure 2

A heatmap of the matrix of Pearson’s correlation coefficient among drought-responsive miRNAs based on their expressions. Red and blue frames represent some examples of significantly positive and negative correlations among miRNAs

Full size image

Correlations of miRNAs with GDP and DT traits

Based on their correlations with GDP and DT, miRNAs could be generally grouped into five clusters (Fig. 3). Cluster Ib contained 39 miRNAs. Their expression levels were generally positively correlated with GDP traits, while their expression/regulation levels were negatively correlated with DT traits. Cluster IIa contained 138 miRNAs. Their expression levels were generally negatively correlated with GDP traits, while their expression/regulation levels were positively correlated with DT traits (Fig. 3). miRNAs in cluster Ib and IIa played opposite roles in regulating GDP and DT. Only a few miRNAs were both positively/negatively correlated with GDP and DT traits (mainly distributed in cluster Ia and cluster IIc) (Fig. 3).

Fig. 3
figure 3

A heatmap of Pearson’s correlation coefficient (PCC) between drought-responsive miRNAs (DRMs) and agronomic (in blue) and drought-tolerant (DT) (in red) traits. Five types of DRMs are at right: Type I, a miRNA is significantly correlated (|PCC| > 0.6) with at least one of measure agronomic traits; Type II, a miRNA is significantly correlated (|PCC| > 0.6) with at least one of measured DT traits; Type III, a miRNA is positively or negatively correlated with both agronomic and DT traits; Type IV, a miRNA is oppositely correlated (|PCC| > 0.6) with measured agronomic and DT traits; Type V, a pair of miRNAs are oppositely correlated (|PCC| > 0.6) with measured agronomic and DT traits

Full size image

Four types of miRNA could be further defined by their correlations with GDP and/or DT using threshold of |PCC| ≥ 0.6. There were 74, 68, 21, and 30 miRNAs classified into type I, II, III, and IV, respectively (Fig. 3). The prediction based on the correlation analysis was partially validated by the miRNAs that have been functionally characterized in rice (Additional file 1: Table S7) [14,15,16, 20, 24, 27, 47, 54, 57]. The regulation of a miRNA in response to drought always tended to enhance DT and had negative impacts on GDP (Table 1). Interestingly, DRMs of the same type were more generally highly correlated (Fig. 2, Additional file 1: Table S8). In addition, DRMs of different types, which played similar roles in DT and/or GDP, possessed higher mean PPCs. For example, the mean PPCs among types I-b, III-b, and IV-b, which were positively correlated with GDP traits, ranged from 0.254~0.335 (Additional file 1: Table S8). Similarly, the mean PPCs between types II-b and III-b, which tended to increase DT, were as high as 0.176. Above results indicated DRMs with similar functions worked together to resist drought (Additional file 1: Table S8).

Table 1 miRNAs of different types in responses to drought
Full size table

We also noticed that a pair of mature miRNAs derived from the same pre-miRNA may sometimes have opposite and independent impacts on GDP and DR. In this study, nine pairs of mature miRNAs derived from the same pre-miRNA demonstrated this pattern (defined as type V) (Fig. 3). For example, OsmiR1870-3p and OsmiR1870-5p, were derived from the same stem-loop structure of pre-OsmiR1870. The expressions of OsmiR1870-3p and OsmiR1870-5p were not correlated (PCC = 0.141, p > 0.05) (Fig. 2). The expression of OsmiR1870-3p was negatively correlated with plant height (PCC = -0.802) and biomass (PCC = -0.664) (Fig. 3), indicating a negative role in regulating rice growth and productivity. The expression of OsmiR1870-5p were positively correlated with AOC (PCC = 0.602), relative seed-setting ratio (PCC = 0.826), relative seed weight (PCC = 0.826), and relative biomass (PCC = 0.996) (Fig. 3), indicating its positive role in rice DT.

Time course of the regulation of miRNAs

Based on the regulation of miRNA expression in response to drought (log2FC), 354 DRMs formed five major time course clusters (Fig. 4). Cluster− 1 contained 37 DRMs (5, 3, 3, and 3 for types I, II, III, and IV, respectively), which were gradually downregulated throughout the progress of drought. Cluster-2 also contained 37 DRMs (18, 5, 2, and 3 for types I, II, III, and IV, respectively), which were highly upregulated starting at time point D2, particularly at the later drought time points (D3-D5). This indicated that DRMs in cluster-2 might be associated with DT in the late stage. Cluster-3 contained 18 DRMs (1, 4, 0, and 11 for types I, II, III, and IV, respectively), which were significantly upregulated at the early drought time points (D1 and D2). This indicated that DRMs in cluster-2 might be associated with DT in the early stage. Cluster-4 contained 28 DRMs (2, 8, 3, and 3 for types I, II, III, and IV, respectively), which had significant changes in expression (upregulation and/or downregulation) at one or two time points. Cluster-5, which could be further divided into two sub-clusters (5a and 5b), contained 234 DRMs. In particular, 143 DRMs (31, 30, 8, and 7 for types I, II, III, and IV, respectively) in cluster-5b exhibited greater fold changes than DRMs in the cluster-5a and were gradually upregulated throughout the period of drought, indicating that they have important roles in rice DT.

Fig. 4
figure 4

A heatmap of time-series regulations of drought-responsive miRNAs (DRMs) during drought period. The regulation of a DRM is quantified by Log2 (its expression under drought/ its expression under well-watered condition). Five major clusters (1~5) are generated by hierarchical clustering (Euclidean method)

Full size image

At the recovery stage, 80 RRMs could be divided into 5 types based on their regulatory patterns from the D5 to R stages (Additional file 1: Table S9). In patterns A and B, 52 RRMs maintained their upregulation and downregulation at recovery (Additional file 1: Table S9). In patterns C and D, 24 RRMs had recovered from their regulatory patterns during drought conditions (Additional file 1: Table S9). The regulatory patterns of another four RRMs were varied among the six genotypes (coded as Type E).

Impacts of DRMs and RRMs on the rice transcriptome

Predicted by both TargetFinder and psRobot, 195 DRMs had 612 target genes in total, which contained 202 DRGs (Additional file 1: Table S10). Many of these target DRGs were characterized as with GDP- and/or DT-associated genes (Additional file 1: Table S10). Noticeably, 35.3% (166 out of 470) and 18.1% (85 out of 470) predicted target DRGs were significantly (PCC > 0.242, p < 0.05) and moderately (|PCC| > 0.400) correlated with its regulatory miRNAs (Additional file 1: Table S11). This result indicated that correlations between miRNA and DRGs among genotypes can provide important cues to determine the candidate target gene for a DRM. For example, expressions of a predicted target gene LOC_Os12g40890 was positively correlated (PCC = 0.464) with its regulatory osa-miR408-5p (Additional file 1: Table S12). As expected, the expression of LOC_Os12g40890 was increased in the overexpression lines of pre-OsmiR408 (Additional file 2: Figure S6), which potentially regulated DT in the transgenic lines. However, we still noticed that a proportion of predicted target genes were not correlated with its regulatory miRNA among genotypes. For example, expressions of OsUCL8 (LOC_Os03g50140) and LOC_Os08g37670, which have been verified as the target of osa-miR408-3p [37, 50], was not correlated with osa-miR408-3p among genotypes (Additional file 1: Table S12). It means that target genes were not always significantly correlated with its regulatory miRNA among genotypes.

Besides the direct regulation on target genes, miRNAs could also have profound indirect influences on the rice transcriptome, as revealed by the large number of correlated DRGs. A DRM was highly (|PCC| > 0.6) and moderately (|PCC| > 0.4) correlated with ~ 90.0 (ranging from 0 to 1630) and 868.4 (ranging from 26 to 2895) DRGs on average (Additional file 1: Table S11; Additional file 3). Based on their highly correlated DRGs, these DRMs and RRMs were involved in diverse GO biological processes (Additional file 2: Figure S7) and KEGG pathways (Additional file 2: Figure S8).

Different types of DRMs had their preferential categories of GOBPs (Fig. 5) and KEGG pathways (Additional file 2: Figure S8). For example, DRMs of type I tended to be involved in reproduction, cellular component organization and biogenesis, cell homeostasis, response to external stimulus, response to abiotic stimulus, and response to endogenous stimulus. DRMs of type II tended to be involved in the generation of precursor metabolites and energy, photosynthesis, lipid metabolic process, transport, cell death, biosynthetic process, anatomical structure morphogenesis, and post-embryonic development. DRMs of type III tended to be involved in carbohydrate metabolic process, generation of precursor metabolites and energy, protein metabolic process, and secondary metabolic process. DRMs of type IV tended to be involved in DNA metabolic process and post-embryonic development (Fig. 5). At the recovery stage, RRMs in patterns A and B tended to be involved in reproduction, cellular component organization and biogenesis, nucleic acid metabolic process, and catabolic process. RRMs in patterns C and D tended to be involved in carbohydrate metabolic process, protein modification process, transport, response to stress, cell death, and signal transduction (Fig. 5).

Fig. 5
figure 5

A heatmap of preferential index (PI) of five types of drought-responsive miRNAs (DRMs) and two patterns of recovery-related miRNAs (RRMs) in the categorized Gene Ontology (GO) biological processes. A high PI indicates the type of DRMs or the pattern of RRMs preferentially anticipates the certain biological process

Full size image

The transcriptomic impact of a miRNA and its involvement in certain biological processes can provide additional information for us to understand its role in DT and/or GDP. For instance, OsmiR1870-3p and OsmiR1870-5p were predicted to be associated with GDP and DT, respectively (Fig. 3). OsmiR1870-3p were correlated with 289 DRGs (absolute PCC > 0.4, p < 0.05), while OsmiR1870-5p were correlated with 1598 DRGs (absolute PCC > 0.4, p < 0.05) (Additional file 3, Additional file 2: Figure S9a). DRGs correlated with OsmiR1870-3p and OsmiR1870-5p were rarely overlapped (Additional file 2: Figure S9a). Among the correlated DRGs, there were many characterized genes of GDP and DT (Additional file 1: Table S13). Moreover, the result of GO enrichment indicated that OsmiR1870-5p were involved in response to stimuli (GO:0050896) and response to stress (GO:0006950) (Additional file 2: Figure S9c). This partially supported the predictions of OsmiR1870-5p was associated with DT.

Candidate miRNAs of drought-tolerance

The DRMs from types II, III, and V, which were associated with DT but without potential adverse effects on GDP, could be valuable candidates for the improvement of rice DT. A miRNA has a high probability to be associated with DT if it has direct or indirect impacts on the expression of DT genes. Based on the above considerations, 24 DRMs in 15 families were recommended as good candidates for DT improvement. These DRMs had 2–23 moderately correlated (|PPC| > 0.4) known DT genes (Table 2, Additional file 3). They belonged to time course cluster-2 and cluster-5b, which play roles in DT at different stages. Based on their highly correlated DRGs (|PPC| > 0.6), the candidate miRNAs anticipated many biological processes (Additional file 2: Figure S7) and metabolic pathways (Additional file 2: Figure S8) relevant to DT. For example, osa-miR399k may enhance DT by regulating water transport (GO:0006833) and fluid transport (GO:0042044) (Additional file 2: Figure S7), while osa-miR169 (g and f.1) may regulate DT by thiamine (ko00730) and vitamin B6 (ko00750) metabolism (Additional file 2: Figure S8).

Table 2 Candidate miRNAs for drought tolerance (DT). A correlated DT gene means its absolute Pearson’s correlation coefficient is > 0.4
Full size table

Discussion

Time course of the regulation of miRNAs under drought conditions provides key information to identify miRNAs associated with DT

miRNAs are in the response to drought stress and play essential roles in rice DT [7, 49]. Many previous studies have described responses of miRNAs to drought in different rice genotypes [3, 5, 6] and tissues [2, 18]. The common design of these experiments is to generate drought-responsive miRNAs from a single genotype [2, 6, 56] or a pair of genotypes with contrasting tolerance [3, 5]. The involvement of a limited number of genotypes and time points when investigating DRGs may result in the loss of valuable information [17]. By time course investigation of miRNAs in response to drought among six rice genotypes, we identified greater number of conserved drought-responsive miRNAs than previous studies (354 mature miRNA in 57 families). Genotype-specific expression and regulation of a miRNA in response to drought is common [33, 49]. The expression and regulation of a miRNA could be randomly different between two genotypes. Therefore, the identification of candidate miRNAs simply by generating differentially regulated miRNAs from two genotypes with contrasting drought tolerances may lead to false-positive results. By including a considerable number of samples, we can apply a correlation analysis to build the association of a miRNA with DT/GDP. The correlation-based prediction is reliable, as it is validated by many former functional studies of miRNAs. This correlation-based prediction provides us the key information to identify miRNAs associated with DT.

We also describe time course regulations of DRMs of different types during drought periods. These results can provide additional clues to understanding the role a miRNA plays in DT and/or GDP. For instance, the DRMs in cluster− 1 and cluster-5b were regulated in a drought-dependent manner, with regulation levels that gradually increased along with drought progression. GDP-associated DRMs in this cluster were passively regulated to slow down rice growth, while DT-associated DRMs in this cluster had minor but additive impacts on DT. The DRMs in cluster-2 were largely upregulated starting at D2, indicating their crucial roles during later drought periods, such as the strong activation of osmotic adjustment. In contrast, the DRMs in cluster-3 play a role in early drought. It is noteworthy that cluster-3 contains a high proportion of DRMs that are inversely correlated with both DT and GDP (type IV) and potentially have pleiotropic effects on GDP. Their activations limited in the early drought may bring less penalties on GDP.

Complex regulatory networks of miRNAs in rice DT and/or GDP

First, we detected a complex regulatory network of DRMs based on the correlation analysis. We found that DRMs of the same type worked together to activate the tolerant mechanisms and/or to inhibit rice GDP. Second, a miRNA could have profound impacts on the plant transcriptome [7, 38]. In our study, we observed that a large number of DRGs are highly correlated with DRMs. Which make them be involved in diverse GO biological processes and/or KEGG metabolic pathways. The involvement of DRMs in certain GO biological processes and KEGG pathways could provide additional information to understand the role of a miRNA plays in DT and/or GDP. For example, DT-associated DRMs tended to be involved in energy metabolism (e.g., photosynthesis, generation of precursor metabolites and energy). This result highlights the point that DT in a crop requires the ability to maintain normal GDP under drought conditions, rather than merely to survive [32]. The GDP-associated DRMs (miRNAs of type I) tended to be involved in stress responses (e.g., tropism, cell homeostasis, response to abiotic stimulus). This indicates that the miRNA-mediated regulation of GDP is relevant to environmental adaptation. We also detected a few (21 out of 354) miRNAs belonging to type III. They tended to be involved in carbohydrate metabolic processes (e.g., trehalose metabolic process, trehalose biosynthetic process). These metabolic or biological processes have been reported to be closely related to better productivity under drought conditions [10, 34]. Therefore, DRMs of type III are valuable candidates for improving DT.

Tradeoff between DT and GDP by miRNAs

The activation of DT mechanisms is always coordinated with multiple developmental processes [11, 35]. Given the central role of miRNAs in the regulatory network [7, 38], a DT-associated miRNA can have pleiotropic effects on GDP [38, 46]. It represents as inversely correlations of a miRNA with traits of GDP and DT, potentially results in tradeoffs between DT and GDP. Based on our prediction, approximately ~ 14.2% of miRNAs (30 out of 211) can have potential negative pleiotropic effects on GDP. For instance, OsmiRNA439(a-i) negatively regulates OsGDCH (LOC_Os10g37180, PCC = -0.624) [26] while it positively regulates OsPIN5b (LOC_Os08g41720, PCC = 0.653) [30]. Moreover, OsmiRNA439(a-i) positively regulates OsNAC5 (LOC_Os11g08210, PCC = 0.569) [19] and OsALDH10A5 (LOC_Os04g39020, PCC = 0.582) [39], which act as positive regulators in DT. As a result, the upregulation of OsmiRNA439(a-i) under drought stress leads to inhibition of GDP and enhancement of DT by these DRGs.

Such a tradeoff could also be caused by a pair of mature miRNAs from the same pre-miRNA. These miRNAs commonly posttranscriptionally regulate different target genes [6], which may sometimes have opposite roles in regulating GDP and DT, as indicated by previous studies [37, 50] and by this study. This arrangement occurs at a frequency of 21.9% (nine pairs from type V out of 41 pairs in total) across the rice genome. We provide an example of OsmiR1870-3p and OsmiR1870-5p in this study. OsmiR1870-3p negatively regulates GDP while OsmiR1870-5p positively regulates DT via their different impacts on the rice transcriptome. Pre-OsmiRNA408, which forms two mature miRNAs (OsmiRNA408-3p and OsmiRNA408-5p), is another example. As reported by two independent studies, overexpression of pre-OsmiRNA408 in rice decreases rice drought tolerance [37] while it positively regulates grain yield [50]. It has been reported that OsmiRNA408-3p positively regulates grain yield by suppressing UCL8 (LOC_Os03g50140) [50]. We speculate that OsmiRNA408-5p negatively regulates drought tolerance by suppressing its predicted target gene LOC_Os12g40890 (OsIAA30). Different roles of two mature miRNAs from the same pre-miRNA in DT and GDP raise a challenge for their utilization. We cannot simply manipulate the pre-miRNA, which is regularly applied in crop improvement [37, 50], to improve a given agronomic trait or obtain stress tolerance.

Valuable candidates for the improvement of drought tolerance

While a miRNA can have systematic impacts on DT, it can sometimes have many pleiotropic effects [38, 46]. We should avoid unwanted pleiotropic effects on GDP by a DT gene when applying it in breeding strategies. In this study, we recommend some candidates of type II, which have low probabilities of unwanted pleiotropic effects on GDP. We also recommend four miRNAs of type III, which can be advantageous in both DT and GDP. Among all candidates, miRNA398a/b in Medicago truncatula [42] has been reported to be associated with DT. Three candidates, miRNA169g [55], miRNA398a [3], and miRNA5505 [49], have been recommended by other studies on rice. Given the considerable number of samples and the correlation-based analyses, we also identified many other potential candidates associated with DT. These miRNAs have great prospects in developing DT in rice.

Conclusion

In this study, we detected 344 drought responsive miRNAs among six rice genotypes in total and described their temporal regulations along a progressive drought. These miRNAs can be divided into five categories based on their associations with drought tolerance and productivity. About 15% miRNAs possessed adverse impacts on rice drought tolerance and productivity, which limited their applications in breeding drought tolerant rice cultivars. Meanwhile, we also identified 21 drought responsive miRNAs that can have both advantages in drought tolerance and productivity once they were appropriately activated. The information of potential associations of miRNAs between rice drought-tolerance and productivity based on the correlation analysis can facilitate our utilization of miRNA in breeding.

Methods

Plant materials

Six rice genotypes were selected to investigate the genome-wide responses of miRNAs to long-term progressive drought, as well as their morphological and physiological performances. The materials were collected from International Rice Research Institute and conserved in the seed bank of Shanghai Agrobiological Gene Center (http://seed.sagc.org.cn/). The six rice genotypes have great differences in drought tolerance, which makes them a good system to study the association of miRNAs with drought tolerance. Based on our pre-evaluation and estimations in this study, genotype S3, S11, and S28 possess better drought tolerance, while genotype S9, S18, and S24 are susceptible to drought. Meanwhile, they are important breeding materials for water-saving and drought-resistant rice [31]. The knowledge of expression patterns of a drought-tolerant miRNA among the six genotypes can provide informative cues for us to utilize it by using the appropriate genotype in breeding. We made a minor adjustment in their dates of germination and transplanting to ensure they can meet the drought-stress before heading (Additional file 1: Table S14).

Field experiment

The field experiment was conducted in a drought-resistance facility at Baihe Experimental Station, Shanghai, in 2014. The canopy of the facility was normally opened, and could be closed on rainy days to enable continuous drought. The depth of the soil layer in the drought-treated field was limited to 30 cm. With this design, the development of roots below this depth was equal between both cultivars and therefore the differences in drought-avoidance could be largely mitigated [32]. Moreover, planting rice in the experimental field rather than in pots led to more homogenous levels [32]. A similar field nearby was well-watered during the experiment as the control treatment (21.4% soil-water content). Plants for each genotype were planted in a plot of 10 rows × 10 hills at 18 cm intervals. Three replicates were set up for each genotype in the drought-treated and well-watered fields. The field arrangement was followed with a single factor randomized block design. The dates for germination and transplanting for the six genotypes were adjusted to make their heading during drought. We started the drought treatment on the 16th of July and continued the artificial drought for 38 days. We measured the soil water content at a depth of 30 cm in the drought-treated field every 3–5 days to monitor the progress of the drought. The drought-treated field was re-watered on the afternoon of the 22nd of August until the SWC decreased to ~ 10% (Additional file 2: Figure S10).

Physiological traits were measured using three replicates of leaf samples. Each replicate contained the three top leaves of the main tillers from each plot. A further three replicates of leaf samples were collected for miRNA and RNA sequencing. We harvested the leaf samples at six time points: the 24th of July (D1, tillering stage), the 29th of July (D2, booting stage), the 5th of August (D3, booting stage), the 11th of August (D4, heading stage), the 22nd of August (D5, heading stage), and the 23rd of August (R, recovery). All leaf samples were strictly collected between 13:00 and 14:00 for each time point and then put into liquid nitrogen immediately.

Osmolality were measured to reflect the osmotic adjustment of the leaf samples. We conducted the measurement of the osmolality via the Vapro™ vapor pressure osmometer (Wescor Model 5600). We measured the total antioxidant capacity (AOC), which reflects the capacity for ROS scavenging, using the total antioxidant capacity assay kit (product#A015, Nanjing Jiancheng Bioengineering Institute, Jiangsu, China). We measured the content of H2O2 and dead leaf ratio to reflect the physiological damage caused by drought stress. The content of H2O2 was measured by an H2O2 assay kits (product#A064, Nanjing Jiancheng Bioengineering Institute, Jiangsu, China). Measurements were taken from D1 to D5 during the drought period. The dead leaf ratio was estimated from five plants per plot on the 23rd of August. We also measured a number of important agronomic traits to reflect the growth, development, and reproduction of each genotype. Plant height and number of tillers were measured from five individual plants in each plot. The number of seeds, seed-setting ratio, biomass, and grain weight for each genotype were measured after harvest from eight plants per plot. Relative performances (the value of a trait in drought/that in CK) were calculated for these agronomic traits to reflect their DT. We applied independent t-tests to detect any significant differences in measured physiological and agronomic traits between the treatments at each time point.

RNA extraction, library construction, and Illumina HiSeq sequencing

Total RNA was extracted from the leaf tissue using PureLink® Plant RNA Reagent (Life Technologies) according to the manufacturer’s instructions. Genomic DNA was removed using DNase I (Takara). Then, RNA quality was determined by a 2100 Bioanalyser (Agilent) and quantified using an ND-2000 (NanoDrop Technologies). Only high-quality RNA samples (OD260/280 = 1.8~2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5, 28S:18S ≥ 1.0, > 10 μg) were used to construct a sequencing library. The concentration of RNA and the purity of the samples were estimated using a NanoDrop Spectrophotometer (Thermo Fisher Scientific) and Qubit Fluorometer (Thermo Fisher Scientific), respectively. An aliquot of the sample was also separated on an Agilent RNA Bioanalyzer chip to check for integrity. Libraries for sRNA sequencing were prepared with a TruSeq Small RNA Sample Preparation Guide (Illumina, USA) and sequenced on an Illumina NextSeq500 platform. Before library construction, three replicates of RNA samples from each genotype from one time point were mixed together. Total RNA libraries were constructed following the specifications in the TruSeq® RNA Sample Preparation v2 Guide (Illumina) and sequenced on an Illumina HiSeq 2500 platform. Three replicates of RNA samples for S3, S11, and S28 were mixed together to construct the libraries, while three replicates of RNA samples for S9, S18, and S24 were independently constructed from the libraries. High-throughput sequencing was conducted at Shanghai Majorbio Biopharm Technology Co., Ltd. (Shanghai, China). The raw data of miRNA (Additional file 1: Table S5) and mRNA (Additional file 1: Table S15) for each sample by high-throughput sequencing will be submitted to the NCBI Sequence Read Archive (SRA) under the number PRJNA609211 before publication.

Data analysis

Identification and quantification of conserved miRNAs

The raw reads were trimmed and quality controlled by Fastx-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). After the identical sequences were merged, all clean reads were separately aligned to the reference genome (Nipponbare, msu7.0) (http://rice.plantbiology.msu.edu/) using Bowtie 2 [23]. The reads that could be aligned to the reference genome were then compared to the database (http://www.mirbase.org/, Release 21, June 2014) to identify conserved miRNAs. We used miRDeep2 (https://www.mdc-berlin.de/content/mirdeep2-documentation) to count the clean reads and quantify the expression level of a detected miRNA by Transcripts per Million (TPM). A miRNA was used for further analysis if the TPM was > 0.10 in at least one sample. The expression data of the available miRNAs from each sample are provided in Additional file 4. The drought-responsive miRNA (DRMs) and recovery-related miRNAs (RRMs) for each genotype were determined using p-value< 0.05 and |log2FC| ≥ 1 during the drought period and at the recovery stage, respectively. We selected six DRMs to validate their expression in 10–15 randomly chosen samples by qPCR (Additional file 2: Figure S11). The primers for the miRNAs and the reference U6 are listed in Additional file 1: Table S16.

Determination of differentially expressed genes between samples from CK and treated fields

We used SeqPrep to strip adaptors and/or merge paired reads with overlap into single reads (https://github.com/jstjohn/SeqPrep) and used Sickle to remove low-quality reads (https://github.com/najoshi/sickle). We then assembled the clean data using the software Cufflinks and mapped them to the reference genome (Nipponbare, msu7.0) and to mitochondrial and chloroplast genomes (http://rice.plantbiology.msu.edu/) via Tophat with no more than two base mismatches allowed in the alignment [41]. We determined the gene expression levels with the Fragment Per Kilobase of exon per Million fragments mapped (FPKM) method via the widely used software program Cuffdiff [41]. The gene expression quantified by FPKM can be well validated by qPCR in our previous study using the same samples (Additional file 2: Figure S12).

We determined differentially expressed genes (DEGs) between samples from drought and CK fields for each genotype at six time points. Since S9, S18, and S24 had three biological replicates, we determined their DEGs via a false discovery rate (FDR) < 0.05 and a logarithm two-fold change |log2FC| ≥ 1. Given the mixed nature of the cDNA libraries of S3, S11, and S28, we determined their DEGs with a p-value< 0.05 and |log2FC| ≥ 1. Moreover, as S9, S18, and S24 had three biological replicates, their FPKMs were averaged as a unit in further analysis. The expression data (FPKM) for all DEGs are provided in Additional file 5.

Studying impacts of miRNAs on the rice transcriptome, GDP, and DT by correlation analyses

To construct the regulatory network of the DRMs, correlation analyses (Pearson’s correlation coefficient, PCC) were applied to detect any significant correlations among DRMs based on their expressions. We predicted target genes for each miRNA by TargetFinder and psRobot to investigate the direct impacts of miRNAs on gene expression, particularly of genes known to be associated with GDP and/or DT (Additional file 2: Table S17). Detailed information can be found in the China Rice Data Center (http://www.ricedata.cn/gene/). Correlation analysis was further applied between the expression of DRMs and DEGs (from D1-D5 and R) to study the indirect impacts of miRNAs on the rice transcriptome. |PCC| > 0.6 and |PCC| > 0.4 was set as the threshold value to indicate a high and moderate correlation between a DRM and a DEG, respectively.

To predict the association of a DRM with GDP, correlation analyses were conducted between DRM expression and the measured agronomic and physiological traits among samples during drought (D1 to D5). We also conducted correlation analyses between the average fold change of a DRM during drought and four DT traits (relative no. of panicle, relative seed-setting ratio, relative biomass, and relative grain weight) to predict the association of the DRM with DT. To be more cautious, |PCC| > 0.6 was set as the threshold value to indicate the potential association of a DRM to a given trait. Based on the results of the correlation analysis, four types of DRMs were defined: (1) Type I, a DRM highly correlated with at least one GDP trait; (2) Type II, a DRM highly correlated with at least one DT trait; (3) Type III, a DRM positively or negatively correlated with both GDP and DT traits; and (4) Type IV, a DRM inversely correlated with at least one GDP trait and a DT trait.

We conducted enrichment analyses of Gene Ontology (GO) and KEGG pathways for each DRM based on its highly correlated DRGs. If a DRM had < 30 highly correlated DRGs, it was not included in the GO and KEGG enrichment analyses. The enriched GO biological processes (GOBPs) were further categorized into classifications by the web tool “GO Terms Classifications Counter” (http://www.animalgenome.org/cgi-bin/util/gotreei) [13]. We defined a preferential index (PI) to determine the preferential GO classifications for each type of DRM. PI is calculated as \( \frac{\mathrm{Ng}/\mathrm{Nm}}{\mathrm{Ntg}/\mathrm{Ntm}} \), where Ng indicated no. of enriched GOBPs of a classification by one type of DRMs; Nm indicated no. of DRMs in this type; Ntg indicated no. of enriched GOBPs of a classification by total DRMs, and Ntm indicated no. of total DRMs.

Temporal gene regulation and relevant biological processes in tolerant and susceptible groups

Time course analyses on the regulation of miRNAs from D1 to D5 were conducted via hierarchical clustering (Euclidean method). Fold changes in the expression of each miRNA in the six genotypes from a single time point were averaged in this analysis. Based on this method, miRNAs could be divided into five major clusters based on the time course regulation of the miRNAs. We further defined five patterns of miRNA regulation at recovery (D5 to R): (A) upregulated at both D5 and recovery; (B) downregulated at both D5 and recovery; (C) upregulated at D5 while downregulated at recovery; (D) downregulated at D5 while upregulated at recovery; (E) regulation levels varied among genotypes. The preferences in GO classifications were also analyzed by the PI for RRMs in different patterns.

Validation of a predicted target gene of osa-miR408-5p by qPCR

The fragment including the precursor of osa-miR408 was amplified from genomic DNA of upland rice cv. IRAT109 and placed to pCAMBIA1300-EGFP to control GFP expression. Then it was placed respectively to pCAMBIA1300-EGFP for fusion expression under the control of the CaMV35S promoter. And so, for the miR408 promoter driving GUS expression vector. All the constructs were transformed into the japonica rice cv. Zhonghua− 11 by the Agrobacterium-mediated (stain EHA105) transformation method42. Three transgenic lines (OX− 1~3) were used to test the expression of osa-mR408-5p and its target gene LOC_Os12g40890 by qPCR. U6 and Actin were used as the reference gene for miRNA and mRNA, respectively (Additional file 1: Table S16).

Availability of data and materials

All data supporting the conclusions of this article are provided within the article and its Additional files 1, 2, 3, 4 and 5. Raw data of sequenced sample will be submitted to the NCBI Sequence Read Archive (SRA) under the number PRJNA609211 before publication.

Abbreviations

DRM:

Drought-responsive miRNA

DT:

Drought tolerance

GDP:

Growth, development, and reproduction

GO:

Gene Ontology

GOBP:

Biological processes in Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

PCC:

Pearson correlation coefficient

PI:

Preferential index

RRM:

Recovery-related miRNA

References

  1. Aravind J, Rinku S, Pooja B, Shikha M, Kaliyugam S, Mallikarjuna MG, Kumar A, Rao AR, Nepolean T. Identification, characterization, and functional validation of drought-responsive microRNAs in subtropical maize inbreds. Front Plant Sci. 2017;8:941.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Bakhshi B, Fard ME, Nikpay N, Ebrahimi MA, Bihamta MR, Mardi M, Salekdeh GH. MicroRNA signatures of drought signaling in rice root. PLoS One. 2016;11:e0156814.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  3. Balyan S, Kumar M, Mutum RD, Raghuvanshi U, Agarwal P, Mathur S, Raghuvanshi S. Identification of miRNA-mediated drought responsive multi-tiered regulatory network in drought tolerant rice, Nagina 22. Sci Rep. 2016;7:15446.

    Article  CAS  Google Scholar 

  4. Bernier J, Atlin GN, Serraj R, Kumar A, Spaner D. Breeding upland rice for drought resistance. J Sci Food Agric. 2008;88:927–39.

    Article  CAS  Google Scholar 

  5. Cheah BH, Nadarajah K, Divate MD, Wickneswari R. Identification of four functionally important microRNA families with contrasting differential expression profiles between drought-tolerant and susceptible rice leaf at vegetative stage. BMC Genomics. 2015;16:692.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  6. Chung PJ, Jung H, Jeong DH, Ha SH, Choi YD, Kim JK. Transcriptome profiling of drought responsive noncoding RNAs and their target genes in rice. BMC Genomics. 2016;17:563.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. Ding YF, Tao YL, Zhu C. Emerging roles of microRNAs in the mediation of drought stress response in plants. J Exp Bot. 2013;64:3077–86.

    Article  CAS  PubMed  Google Scholar 

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

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Fang YJ, Xiong LZ. General mechanisms of drought response and their application in drought resistance improvement in plants. Cell Mol Life Sci. 2014;72:673.

    Article  PubMed  CAS  Google Scholar 

  10. Garg AK, Kim JK, Owens TG, Ranwala AP, Choi YD, Kochian LV, et al. Trehalose accumulation in rice plants confers high tolerance levels to different abiotic stresses. Proc Natl Acad Sci U S A. 2002;99:15898–903.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Harb A, Krishnan A, Ambavaram MMR, Pereira A. Molecular and physiological analysis of drought stress in Arabidopsis reveals early responses leading to acclimation in plant growth. Plant Physiol. 2010;154:1254–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Hu H, Xiong L. Genetic engineering and breeding of drought-resistant crops. Annu Rev Plant Biol. 2014;65:715–41.

    Article  CAS  PubMed  Google Scholar 

  13. Hu ZL, Bao J, Reecy JM. CateGOrizer: a web-based program to batch analyze gene ontology classification categories. Online J Bioinform. 2008;9:108–12.

    Google Scholar 

  14. Huang J, Li ZY, Zhao DZ. Deregulation of the OsmiR160 target gene OsARF18 causes growth and developmental defects with an alteration of auxin signaling in rice. Sci Rep. 2016;6:29938.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Guo S, Xu Y, Liu H, Mao Z, Zhang C, Ma Y, Zhang Q, Meng Z, Chong K. The interaction between OsMADS57 and OsTB1 modulates rice tillering via DWARF14. Nat Commun. 2013;4:1566.

    Article  PubMed  CAS  Google Scholar 

  16. Gao F, Wang K, Liu Y, Chen YP, Chen P, Shi ZY, Luo J, Jiang DQ, Fan FF, Zhu YG, Li SQ. Blocking miR396 increases rice yield by shaping inflorescence architecture. Nat Plants. 2015;15:196.

    Google Scholar 

  17. Jeong DH, Green PJ. The role of rice microRNAs in abiotic stress responses. J Plant Biol. 2013;56:187–97.

    Article  CAS  Google Scholar 

  18. Jeong DH, Park S, Zhai JX, Gurazada SGR, Paoli ED, Meyers BC, Green PJ. Massive analysis of rice small RNAs: mechanistic implications of regulated microRNAs and variants for differential target RNA cleavage. Plant Cell. 2011;23:4185–207.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Jeong JS, Kim YS, Redillas M, Jang G, Jung H, Bang SW, Choi YD, Ha SH, Reuzeau C, Kim JK. OsNAC5 overexpression enlarges root diameter in rice plants leading to enhanced drought tolerance and increased grain yield in the field. Plant Biotechnol J. 2013;11:101–14.

    Article  CAS  PubMed  Google Scholar 

  20. Jiao YQ, Wang YH, Xue DW, Wang J, Yan MX, Liu GF, Dong GJ, Zeng DL, Lu ZF, Zhu XD, Qian Q, Li JY. Regulation of OsSPL14 by OsmiR156 defines ideal plant architecture in rice. Nat Genet. 2010;42:541–5.

    Article  CAS  PubMed  Google Scholar 

  21. Jian H, Wang J, Wang T, Wei L, Li J, Liu L. Identification of rapeseed microRNAs involved in early stage seed germination under salt and drought stresses. Front Plant Sci. 2016;7:658.

    PubMed  PubMed Central  Google Scholar 

  22. Kozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 2019;47:D155–62.

    Article  CAS  PubMed  Google Scholar 

  23. Langmead B, Salzberg S. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Lee YS, Lee DY, Cho LH, An GH. Rice miR172 induces flowering by suppressing OsIDS1 and SNB, two AP2 genes that negatively regulate expression of Ehd1 and florigens. Rice. 2014;7:31.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Li WX, Oono Y, Zhu J, He XJ, Wu JM, Iida K, Lu XY, Cui X, Jin H, Zhu JK. The Arabidopsis NFYA5 transcription factor is regulated transcriptionally and posttranscriptionally to promote drought resistance. Plant Cell. 2008;20:2238–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Lin HC, Karki S, Coe RA, Bagha S, Khoshravesh R, Balahadia CP, Sagun JR, Tapia R, Israel WK, Montecillo F, de Luna A, Danila FR, Lazaro A, Realubit CM, Acoba MG, Sage TL, von Caemmerer S, Furbank RT, Cousins AB, Hibberd JM, Quick WP, Covshoff S. Targeted Knockdown of GDCH in Rice Leads to a Photorespiratory-Deficient Phenotype Useful as a Building Block for C4 Rice. Plant Cell Physiol. 2016;57:919–32.

    Article  CAS  PubMed  Google Scholar 

  27. Liu H, Jia SH, Shen DF, Liu J, Li J, Zhao HP, Han SC, Wang YD. Four AUXIN RESPONSE FACTOR genes downregulated by microRNA167 are associated with growth and development in Oryza sativa. Funct Plant Biol. 2012;39:736–44.

    Article  CAS  PubMed  Google Scholar 

  28. Liu MM, Yu HY, Zhao GJ, Huang QF, Lu YE, Ouyang B. Identification of drought-responsive microRNAs in tomato using high-throughput sequencing. Identification of drought-responsive microRNAs in tomato using high-throughput sequencing. Funct Integr Genomics. 2018;18:67–78.

    Article  CAS  PubMed  Google Scholar 

  29. Liu Q, Zhang YC, Wang CY, Luo YC, Huang QJ, Chen SY, Zhou H, Qu LH, Chen YQ. Expression analysis of phytohormone-regulated microRNAs in rice, implying their regulation roles in plant hormone signaling. FEBS Lett. 2009;583:723–8.

    Article  CAS  PubMed  Google Scholar 

  30. Lu GW, Coneva V, Casaretto JA, Ying S, Mahmood K, Liu F, Nambara E, Bi YM, Rothstein SJ. OsPIN5b modulates rice (Oryza sativa) plant architecture and yield by changing auxin homeostasis, transport and distribution. Plant J. 2015;83:913–25.

    Article  CAS  PubMed  Google Scholar 

  31. Luo LJ. Breeding for water-saving and drought-resistance rice (WDR) in China. J Exp Bot. 2010;61:3509–17.

    Article  CAS  PubMed  Google Scholar 

  32. Ma X, Xia H, Liu Y, Wei H, Zheng X, Song C, Chen L, Liu H, Luo LJ. Transcriptomic and metabolomic studies disclose key metabolism pathways contributing to well-maintained photosynthesis under the drought and the consequent drought tolerance in rice. Front Plant Sci. 2016;7:1886.

    PubMed  PubMed Central  Google Scholar 

  33. Mutum RD, Balyan SC, Kansal S, Agarwal P, Kumar S, Kumar M, Raghuvanshi S. Evolution of variety-specific regulatory schema for expression of osa-miR408 in indica rice varieties under drought stress. FEBS J. 2013;280:1717–30.

    Article  CAS  PubMed  Google Scholar 

  34. Nuccio M, Wu J, Mowers R, Zhou H, Meghji M, Primavesi LF, et al. Expression of trehalose-6-phosphate phosphatase in maize ears improves yield in well-watered and drought conditions. Nat Biotechnol. 2015;33:862–9.

    Article  CAS  PubMed  Google Scholar 

  35. Pandey V, Shukla A. Acclimation and tolerance strategies of rice (Oryza sativa L.) under drought stress. Rice Sci. 2015;22:3.

    Article  Google Scholar 

  36. Sun G. MicroRNAs and their diverse functions in plants. Plant Mol Biol. 2012;80:17–36.

    Article  CAS  PubMed  Google Scholar 

  37. Sun MZ, Yang JK, Cai XX, Shen Y, Cui N, Zhu YM, Jia BW, Sun XL. The opposite roles of OsmiR408 in cold and drought stress responses in Oryza sativa. Mol Breed. 2018;38:12.

    Article  CAS  Google Scholar 

  38. Tang JY, Chu CC. MicroRNAs in crop improvement: fine-tuners for complex traits. Nat Plants. 2017;3:17077.

    Article  PubMed  Google Scholar 

  39. Tang W, Sun JQ, Liu J, Liu FF, Yan J, Gou XJ, Lu BR, Liu YS. RNAi-directed downregulation of betaine aldehyde dehydrogenase 1 (OsBADH1) results in decreased stress tolerance and increased oxidative markers without affecting glycine betaine biosynthesis in rice (Oryza sativa). Plant Mol Biol. 2014;86:443-54.

    Article  CAS  PubMed  Google Scholar 

  40. Tian CJ, Zuo ZL, Qiu JL. Identification and characterization of ABA-responsive microRNAs. J Genet Genomics. 2015;42:393–402.

    Article  CAS  PubMed  Google Scholar 

  41. 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:562–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Trindade I, Capitão C, Dalmay T, Fevereiro MP, dos Santos DM. miR398 and miR408 are up-regulated in response to water deficit in Medicago truncatula. Planta. 2010;231:705–16.

    Article  PubMed  CAS  Google Scholar 

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

    Article  CAS  PubMed  Google Scholar 

  44. Wang J, Mei J, Ren G. Plant microRNAs: biogenesis, homeostasis, and degradation. Front Plant Sci. 2019;10:360.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Xia H, Luo Z, Xiong J, Ma XS, Lou QJ, Wei HB, Qiu J, Yang H, Liu GL, Fan LJ, Chen L, Luo LJ. Bi-directional selection in upland rice leads to its adaptive differentiation from lowland rice in drought resistance and productivity. Mol Plant. 2019;12:170–84.

    Article  CAS  PubMed  Google Scholar 

  46. Xia K, Wang R, Ou X, Fang Z, Tian C, Duan J, Wang Y, Zhang M. OsTIR1 and OsAFB2 downregulation via OsmiR393 overexpression leads to more tillers, early flowering and less tolerance to salt and drought in rice. PLoS One. 2012;7(1):e30039.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Xia K, Ou X, Tang H, Wang R, Wu P, Jia YX, Wei XY, Xu XL, Kang SH, Kim SK, Zhang M. Rice microRNA osa-miR1848 targets the obtusifoliol 14a-demethylase gene OsCYP51G3 and mediates the biosynthesis of phytosterols and brassinosteroids during development and in response to stress. New Phytol. 2015;208:790–802.

    Article  CAS  PubMed  Google Scholar 

  48. Xie FL, Wang QL, Sun RR, Zhang BH. Deep sequencing reveals important roles of microRNAs in response to drought and salinity stress in cotton. J Exp Bot. 2015;66:789–804.

    Article  CAS  PubMed  Google Scholar 

  49. Zhang BH. MicroRNA: a new target for improving plant tolerance to abiotic stress. J Exp Bot. 2015;66:1749–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Zhang JP, Yu Y, Feng YZ, Zhou YF, Zhang F, Yang YW, Lei MQ, Zhang YC, Chen YQ. miRNA MiR408 regulates grain yield and photosynthesis via a phytocyanin protein. Plant Physiol. 2017;175:1175–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Zhang JS, Zhang H, Srivastava AK, Pan YJ, Bai JJ, Fang JJ, Shi HZ, Zhu JK. Knockdown of rice microRNA166 confers drought resistance by causing leaf rolling and altering stem xylem development. Plant Physiol. 2018;176:2082–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Zhang Q. Strategies for developing green super rice. Proc Natl Acad Sci U S A. 2007;104:16402–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Zhang XH, Zou Z, Gong PJ, Zhang JH, Ziaf K, Li HX, Xiao FM, Ye ZB. Over-expression of microRNA169 confers enhanced drought tolerance to tomato. Biotechnol Lett. 2011;33:403–9.

    Article  CAS  PubMed  Google Scholar 

  54. Zhang YC, Yu Y, Wang CY, Li ZY, Liu Q, Xu J, Liao JY, Wang XJ, Qu LH, Chen F, Xin P, Yan C, Chu J, Li HQ, Chen YQ. Overexpression of microRNA OsmiR397 improves rice yield by increasing grain size and promoting panicle branching. Nat Biotechnol. 2013;31:848–52.

    Article  CAS  PubMed  Google Scholar 

  55. Zhao B, Liang RQ, Ge LF, Lia W, Xiao HS, Lin HX, Ruan KC, Jina YX. Identification of drought-induced microRNAs in rice. Biochem Biophys Res Commun. 2007;354:585–90.

    Article  CAS  PubMed  Google Scholar 

  56. Zhou LG, Liu YH, Liu ZC, Kong DY, Duan M, Luo LJ. Genome-wide identification and analysis of drought-responsive microRNAs in Oryza sativa. J Exp Bot. 2010;61:4157–68.

    Article  CAS  PubMed  Google Scholar 

  57. Zhu QH, Upadhyaya NM, Gubler F, Helliwell CC. Over-expression of miR172 causes loss of spikelet determinacy and floral organ abnormalities in rice (Oryza sativa). BMC Plant Biol. 2009;9:149.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

Download references

Acknowledgements

We thank Mr. Guo Yuntao, Mr. Guo Youbing, and Miss Shi Caiping for their kindly help in data analyses.

Funding

This work is supported by Shanghai Agriculture Applied Technology Development Program, China (Grant No. 2017-02-08-00-08-F00071), Shanghai Natural Science Foundation (17ZR1425500), and the Youth Talent Development Project supported by the Shanghai Municipal Agriculture Commission 2017 (Grant No. 1–33). The funders had no role in experimental design, data collection, data analysis, and manuscript preparation.

Author information

Author notes

  1. Xia Hui and Yu Shunwu contributed equally to this work.

Authors and Affiliations

  1. Shanghai Agrobiological Gene Center, Shanghai, China

    Hui Xia, Shunwu Yu, Deyan Kong, Jie Xiong, Xiaosong Ma, Liang Chen & Lijun Luo

Authors
  1. Hui Xia
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Shunwu Yu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Deyan Kong
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Jie Xiong
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Xiaosong Ma
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Liang Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Lijun Luo
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

HX, LC, and LJL designed this study. HX, SWY, DYK, XSM, and JX conducted field and molecular experiments. HX, SWY, and JX analyzed the data. HX, XSM, and LJL collected, developed, and prepared plant materials in this study. HX, SWY, and LJL wrote the manuscript. All authors have carefully read this manuscript and agree with all aspects of the work.

Corresponding authors

Correspondence to Hui Xia or Lijun Luo.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Additional file 1: Table S1.

Agronomic traits measure among six rice genotypes in drought-treated and well-watered fields. * and ** indicate significant differences at p < 0.05 and p < 0.01 by independent t test between samples from drought and well-watered fields. NA is short for “not available”. Table S2. Content of H2O2 (mmol/g FW) measured among six rice genotypes in drought and well-watered fields. * indicates significant differences at p < 0.05 by independent t test between samples from drought and well-watered fields. Table S3. Osmotic potential measured among six rice genotypes in drought and well-watered fields. * and ** indicate significant differences at p < 0.05 and p < 0.01 by independent t test between samples from drought and well-watered fields. Table S4. Total anti-oxidant capacity (U/g FW) measured among six rice genotypes in drought and well-watered fields. * and ** indicate significant differences at p < 0.05 and p < 0.01 by independent t test between samples from drought and well-watered fields. Table S5. Basic information for sequenced samples. Table S6. Number of expressed, available, drought-responsive miRNAs (DRMs), and recovery-related miRNAs (RRMs). Table S7. Validation of expression patterns and correlation-based predictions by former functionally characterized miRNAs in rice. Table S8. Number of recovery-related miRNAs in five regulation patterns from timepoint D5 to recovery stage (R). Table S9. Number of recovery-related miRNAs (RRMs) in five regulation patterns from timepoint D5 to recovery stage (R). “Varied” indicates the regulation of a miRNA varies among genotypes. Table S10. Target genes predicted by TargetFinder and psRobot for drought-responsive miRNAs. GDP is short for growth, development, and reproduction. DT is short for drought tolerance. NA is short for not assigned. Table S11. Number of predicted target genes of differentially expressed miRNAs in drought responsive genes (DRGs) with different Pearson’s Correlation Coefficient (PCC) values. Table S12. Target genes predicted for miR408-3p and miR408-5p and their PPC values in correlations between miRNA and mRNA. Table S13. Pearson’s correlation coefficients (PCCs) of genes related to plant height (PH) and drought-tolerance (DT) with OsmiR1870-3p and OsmiR1870-5p. Table S14. Basic information for six rice genotypes. Table S15. Information for RNA-sequenced samples. Table S16. Primers of six miRNAs and the reference (U6) for qPCR. Table S17. List of GDP (growth, development, and reproduction)- and DT (drought-tolerance)-associated drought-responsive genes. These DT-associated genes have been functionally characterized. These GDP-associated genes are related to plant height, number of tillers, grain yield, and biomass by former studies.

Additional file 2: Figure S1.

Distribution of the clean reads in small RNA libraries. Figure S2. Proportions of small RNAs in different length among genotypes (a), time points (b), and different treatments (c). Figure S3. A heatmap of the frequency of genotypes for drought-responsive miRNAs during drought (D1-D5). Figure S4. A heatmap of the frequency of time points (D1-D5) for drought-responsive miRNAs among six genotypes. Figure S5. Summary of drought-responsive miRNAs (DRMs) and recovery-related miRNAs (RRMs). a. Frequencies of a DRM among genotypes and time-points. b. Venn diagram of DRMs and RRMs. Figure S6. Expressions of osa-mR408-5p and its target gene LOC_Os12g40890 quantified by qPCR in transgenic lines overexpressing pre-miRNA408. Figure S7. A heatmap describing the involvements of DRMs in biological processes based on the Gene Ontology enrichment using their highly correlated drought-responsive genes (|PCC > 0.6|). The blue (1) color indicates significant enrichment (FDR < 0.05), while the grey (0) color indicates no significant enrichment (p > 0.05). Figure S8. A heatmap describing the involvements of DRMs in metabolic pathways based on the KEGG enrichment using their highly correlated drought-responsive genes (|PCC > 0.6|). The blue (1) color indicates the significant enrichment (p < 0.05), while the grey (0) color indicates no significant enrichment (p > 0.05). Figure S9. Impacts of OsmiR1870-3p and OsmiR1870-5p on rice transcriptome. a. Venn diagram of positively (PCC > 0.4) and negatively (PCC < -0.4) correlated drought-responsive genes (DRGs) for OsmiR1870-3p and OsmiR1870-5p. b. GO enrichment by correlated DRGs of OsmiR1870-3p. c. GO enrichment by correlated DRGs of OsmiR1870-5p. PCC, Pearson’s correlation coefficient. GO terms in red, orange, and yellow indicate p < 0.001. p < 0.01, and p < 0.05 in the enrichment analyses, respectively. Figure S10. Soil-water content measured during drought. Figure S11. Fold changes (drought/ well-watered) of six drought-responsive miRNAs quantified by high-throughput sequencing is well validated by qPCR. Figure S12. Correlations between expressions of six drought responsive genes quantified by RNA-seq and qPCR. The data can be also found in the reference Ma et al. [29].

Additional file 3.

Matrix of Pearson’s correlation coefficients (PCCs) between drought-responsive miRNAs and drought-responsive genes.

Additional file 4.

Expressions (by TPM) of detected miRNAs.

Additional file 5.

Expressions (by FPKM) of drought-responsive genes.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Xia, H., Yu, S., Kong, D. et al. Temporal responses of conserved miRNAs to drought and their associations with drought tolerance and productivity in rice. BMC Genomics 21, 232 (2020). https://doi.org/10.1186/s12864-020-6646-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-020-6646-5

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us