Genome-wide identification and analysis of drought-responsive microRNAs in Oryza sativa

Abstract

In addition to regulating growth and development, the most important function of microRNAs (miRNAs) in plants is the regulation of a variety of cellular processes underlying plant adaptation to environmental stresses. To gain a deep understanding of the mechanism of drought tolerance in rice, genome-wide profiling and analysis of miRNAs was carried out in drought-challenged rice across a wide range of developmental stages, from tillering to inflorescence formation, using a microarray platform. Among the 30 miRNAs identified as significantly down- or up-regulated under the drought stress, 11 down-regulated miRNAs (miR170, miR172, miR397, miR408, miR529, miR896, miR1030, miR1035, miR1050, miR1088, and miR1126) and eight up-regulated miRNAs (miR395, miR474, miR845, miR851, miR854, miR901, miR903, and miR1125) were revealed for the first time to be induced by drought stress in plants, and nine (miR156, miR168, miR170, miR171, miR172, miR319, miR396, miR397, and miR408) showed opposite expression to that observed in drought-stressed Arabidopsis. The most conserved down-regulated miRNAs were ath-miR170, the miR171 family, and ath-miR396, and the most conserved up-regulated miRNAs were ptc-miR474 and ath-miR854a. The identification of differentially expressed novel plant miRNAs and their target genes, and the analysis of cis-elements provides molecular evidence for the possible involvement of miRNAs in the process of drought response and/or tolerance in rice.

Drought stress, microarray, miRNA, Oryza sativa, target genes

Introduction

Micro RNAs (miRNAs) are a class of non-protein-coding, single-strand molecules of 20–27 nucleotides in length, that have been reported to regulate gene expression at the post-transcriptional level in plants, animals, nematodes, and fungi (Bartel, 2004; Denli et al., 2004; Dugas and Bartel, 2004; Meister et al., 2004; Bentwich, 2005; Bhattacharyya et al., 2006; Lytle et al., 2007; Ventura et al., 2008). Since the discovery of the first miRNA, lin-4 RNA, which controls the timing of Caenorhabditis elegans larval development (Lee, 1993; Wightman et al., 1993), the number of identified miRNAs in plants and animals has expanded rapidly (Ambros, 2004; Sunkar et al., 2005; Jones-Rhoades et al., 2006; Zhang et al., 2006; Pilcher et al., 2007; Xie et al., 2007; Sunkar and Jagadeeswaran, 2008). These molecules regulate gene expression by binding to and modifying the translation of target mRNAs that contain sequences at least partially complementary to the mature miRNAs. The regulation of the mode of action of miRNAs differs between plants and animals; whereas miRNAs guide translational repression of animal target genes in most cases, in plants miRNA-guided target mRNA cleavage (post-transcriptional repression) seems to be predominant (Llave et al., 2002; Palatnik et al., 2003). In fact it is no longer true that mRNA cleavage is the predominant mechanism in plants since translation inhibition has recently been shown to be widespread (Brodersen and Voinnet, 2009).

In plants, the miRNAs revealed so far mainly affect plant development and physiological processes such as leaf and flower differentiation, flowering time, auxin response, and floral identity (Sunkar et al., 2005; Mallory and Vaucheret, 2006). Among the roles of miRNAs in regulating growth and development processes in plants, the most currently proposed role is their hypersensitivity to biotic and abiotic stresses as well as to diverse physiological processes including nutrient deficiency in Arabidopsis thaliana (Sunkar and Zhu, 2004; Sunkar et al., 2006). Mallory and Vaucheret (2006) have shown that the expression of miR395 and miR399 was induced in low-sulphate and low-phosphate conditions, respectively.

More recently, other miRNAs such as miR168, miR171, and miR396 were found to be responsive to high salinity, mannitol, and cold stress in A. thaliana, thus supporting the hypothesis of a role for miRNAs in the adaptive response to abiotic stress (Liu et al., 2008). This notwithstanding, the total number of miRNAs involved in plant response to drought stress and a more systematic expression analysis for known miRNAs in plants exposed to field-like, water deficit stress conditions have not been reported thus far. Therefore, efforts to identify novel, drought stress-regulated miRNAs in rice and confirmation of their function may help to unravel the mechanism of drought stress tolerance in rice and other crops.

Usually, gene expression studies in stressed plants have deployed different techniques such as low-throughput northern blotting, promoter–reporter gene fusion, real-time quantitative reverse transcription-PCR (RT-PCR), in situ hybridization, high-throughput differential display RT-PCR (DDRT-PCR), and cDNA microarrays. Microarray technology, which allows thousands of genes or miRNAs to be monitored simultaneously, has revolutionized gene expression studies. The comparison of the expression profiles facilitates the identification and cloning of stress-related genes and miRNAs. In this study, miRNA array analyis was applied to detect the expression profile of drought stress-regulated miRNAs in rice (Oryza sativa L.) grown in soil pots. A group of drought-inducible miRNAs were identified and confirmed by detecting the expression patterns of these miRNAs and their deduced targeted mRNAs.

Materials and methods

Plant materials and drought stress treatment

Twenty pre-germinated seeds of rice cv. IRAT109 were planted in soil in a 1/5000a Wagner's pot. Four pots of rice plants were grown in parallel; one was subject to stress treatment at the tillering stage, one was stress challenged at the inflorescence-forming stage, and the other two pots were grown under normal condition; that is, without suffering drought stress caused by withholding of water as in the treatments. At growth stage 6.1 which is characterized by the number of leaves on the main stem of a plant, that is equivalent to the tillering stage in paddy field production, water was withheld to develop drought stress, while control pots received normal watering. Typical stress symptoms, such as leaf rolling, appeared on the eighth day after water withholding (8 DAW) when the absolute water content of the soil decreased to <4.5%. Total RNA was then extracted on 10, 12, and 14 DAW using TRNzol Total RNA Reagent (Tiangen, China) according to the manufacturer's instructions. The growth status of the stressed rice is shown in Fig. 1A. The second drought stress was then imposed when the flag leaf appeared. Stress symptoms started to emerge on 4 DAW and total RNA was extracted on 5 and 6 DAW (Fig. 1B). RNA was extracted from the plants grown in the control pots in the above two phases. These RNA samples were profiled by a miRNA array platform.

Fig. 1.

Growing status of rice plants under drought stress (left) as compared with well-watered conditions (right): (A) 8, 10, 12, and 14 days after water withholding (DAW) at the tillering stage; (B) 0, 5, and 6 DAW at the inflorescence-forming stage.

miRNA microarray analysis

A large number of unique mature miRNAs, including 623 plant miRNAs listed in miRBASE release 10.0 (http://www.sanger.ac.uk/Software/Rfam/mirna/), were used in the preparation of a miRNA chip, and multiple redundant regions were included. Each region was further comprised of a miRNA probe region, which was introduced to detect transcripts of individual miRNAs. In addition, 52 short oligos with no homology to any of the existing miRNA sequence were included as negative controls. Among the negative controls, PUC2PM-20B and PUC2MM-20B were probes that perfectly match and single-base match, respectively, to a 20-mer RNA positive control sequence that was spiked into the RNA samples before labelling. Probes were duplicated five times on a microarray chip.

Microarray assay was performed using a service provider (LC Sciences). The assay started from 2 μg to 5 μg of total RNA sample, which was size fractionated using a YM-100 Microcon centrifugal filter (from Millipore) and the small RNAs (<300 nucleotides) isolated were 3′-extended with a poly(A) tail using poly(A) polymerase. An oligonucleotide tag was then ligated to the poly(A) tail for later fluorescent dye staining; two different tags were used for the two RNA samples in dual-sample experiments. Hybridization was performed overnight on a μParaflo™ microfluidic chip using a micro-circulation pump (Atactic Technologies) (Gao et al., 2004). On the microfluidic chip, each detection probe consisted of a chemically modified nucleotide coding segment complementary to target miRNA or other RNA (control) and a spacer segment of polyethylene glycol (PEG) to extend the coding segment away from the substrate. The detection probes were made by in situ synthesis using PGR (photogenerated reagent) chemistry. The hybridization melting temperatures were balanced by chemical modifications of the detection probes. Hybridization used 100 μl of 6× SSPE buffer (0.90 M NaCl, 60 mM Na₂HPO₄, 6 mM EDTA, pH 6.8) containing 25% formamide at 34 °C. After hybridization, detection used fluorescence labelling employing tag-specific Cy3 and Cy5 dyes. Hybridization images were collected using a laser scanner (GenePix 4000B, Molecular Device) and digitized using Array-Pro image analysis software (Media Cybernetics). Data were analysed by first subtracting the background and then normalizing the signals using a LOWESS filter (locally weighted regression) (Bolstad et al., 2003). Array experiments with the sample and the control were repeated twice. Chip experiments and array analysis were carried out by LC Sciences (Houston, Texas, USA).

Statistical analysis of microarray profiles was carried out adopting a threshold of fold change >2.0 (the value in log2 scale) and P-value <0.001 [analysis of variance (ANOVA)] to identify differentially expressed miRNAs under drought stress.

Expression validation of miRNA using real-time qPCR

Mature miRNA-specific PCR forward primers (sense DNA oligo identical to the entire mature miRNA sequences) were designed according to miRNA sequences, and the NCode™ miRNA First-Strand cDNA Synthesis Kits and qRT-PCR Kits (Invitrogen) were used for miRNA quantitative RT-PCR (qRT-PCR) analysis. Briefly, the same RNA samples used for microarray profiling were also employed, and all the miRNAs in a sample were polyadenylated using poly(A) polymerase and ATP. Following polyadenylation, SuperScript™ III reverse transcriptase and a specially designed Universal RT Primer were used to synthesize cDNA from the tailed miRNA population. The first-strand cDNA was prepared for analysis by qPCR using SYBR Green or SYBR GreenER™ detection reagents, the Universal qPCR Primer provided in the kit, and the forward primer that targets the specific miRNA sequence of interest.

A set of 45 miRNAs with expression levels with at least 1.5-fold change were selected from the array analyses (differentially expressed miRNAs as indicated in Table 1). PCRs were performed using the 7000 Sequence Detection System (Applied Biosystems) with the following cycling parameters: 95 °C for 10 min, 40 cycles of denaturing at 95 °C for 30 s, annealing at 60 °C for 30 s, and extension at 72 °C for 30 s, followed by a disassociation stage (melting curve analysis). Threshold value was empirically determined based on the observed linear amplification phase of all primer sets. Sample cycle threshold (Ct) values were standardized for each template based on the actin gene control primer reaction, and the 2^–ΔΔCT method was used to analyse the relative changes in gene expression from real-time qPCR experiments (Livak and Schmittgen, 2001).

Table 1.

Drought stress-induced miRNAsa identified by microarray analysis

Sampling stage	Down-regulated miRNAs	Fold change	Up-regulated miRNAs	Fold change
Tillering stage
10 DAW	osa-miR171i	–2.18	ppt-miR896	4.06
	ath-miR396a	–2.29	tae-miR1125	2.98
	ath-miR172a	–2.86	ptc-miR474a	2.70
	ath-miR172e	–2.96	ptc-miR474c	2.69
	ath-miR172c	–2.98	ptc-miR474b	2.60
	ptc-miR172g	–3.19	ath-miR854a	2.46
	ptc-miR172i	–3.21
	osa-miR172c	–3.28
	sbi-miR172b	–3.29
	ath-miR396b	–3.32
	sbi-miR172a	–3.37
12 DAW	ath-miR168a	–2.08	ptc-miR474b	5.91
	ath-miR171a	–2.19	ptc-miR474a	5.78
	ptc-miR171c	–2.20	ptc-miR474c	5.72
	osa-miR168a	–2.28	gma-miR319c	3.86
	sof-miR168b	–2.32	osa-miR319a	3.71
	ppt-miR896	–2.42	ppt-miR395	2.44
	osa-miR168b	–2.51	pta-miR159b	2.42
	ptc-miR397b	–2.54	pta-miR319	2.26
	osa-miR397b	–2.58	ath-miR319a	2.06
	ath-miR397a	–2.58
	ath-miR170	–2.62
	mtr-miR171	–2.63
	sof-miR408e	–2.76
	ath-miR396b	–2.74
14 DAW	ath-miR396a	–2.24	ath-miR851-3p	2.99
	ppt-miR1050	–2.28	ath-miR845a	2.95
	smo-miR396	–2.49	ptc-miR474a	2.93
	ath-miR396b	–2.71	ptc-miR474b	2.86
	ppt-miR1035	–3.30	ptc-miR474c	2.84
	ptc-miR396g	–3.87	ppt-miR1026a	2.52
	tae-miR1126	–4.72	ppt-miR901	2.40
			ath-miR854a	2.36
			ppt-miR903	2.34
Reproductive stage
5 DAW	ppt-miR1030j	–2.05	ptc-miR474c	2.76
	ppt-miR1030a	–2.11	ptc-miR474b	2.71
	ath-miR172b	–2.82	ptc-miR474a	2.64
	ppt-miR896	–2.82	ath-miR854a	2.24
			gma-miR319c	2.02
6 DAW	ppt-miR1030j	–2.06	ptc-miR169s	2.74
	ath-miR156h	–2.10	ath-miR169h	2.71
	gma-miR319c	–2.18	ptc-miR169o	2.60
	ppt-miR1030a	–2.21	ptc-miR169v	2.59
	ath-miR396b	–2.40	ath-miR169a	2.57
	ppt-miR529g	–2.49	ath-miR169d	2.46
	ppt-miR529a	–2.51	osa-miR169f	2.32
	ppt-miR529d	–2.58	osa-miR169e	2.31
	ppt-miR896	–2.70	ath-miR169b	2.30
	smo-miR156c	–2.72	smo-miR171b	2.00
	smo-miR1088-5p	–3.45
	ath-miR172b	–4.45

a

miRNAs with fold change over ±2 and P-value <0.001.

Prediction of putative target genes of drought-induced miRNAs

miRU (http://bioinfo3.noble.org/miRNA/miRU.htm), a plant miRNA target prediction server, was used to predict putative miRNA target genes. It reports all potential sequences complementary to an inquiring miRNA sequence with mismatches no more than a specified value for each mismatch type. Each mismatch is converted to a weighted score according to the criterion in Jones-Rhoades and Bartel (2004): that is, a G:U scores 0.5, an indel 2.0, and any other mismatch 1.0. Furthermore, any mismatch other than the G:U wobble in position 2–7 of the 5′ end of an miRNA is penalized by an extra 0.5 to reflect the importance of complementarity to the miRNA 5′ end for target site function (Mallory et al., 2004; Brennecke et al., 2005; Lim et al., 2005). With default settings, the minimal weighed score <3.0 was applied in the prediction.

Accession numbers

The microarray data of rice microRNAs have been deposited in the NCBI/GEO database, accession numbers GLP8015 and GSE14246.

Results

Identification and expression patterns of drought stress-induced miRNAs in rice

The results indicated that 16 miRNAs, namely miR156, miR159, miR168, miR170, miR171, miR172, miR319, miR396, miR397, miR408, miR529, miR896, miR1030, miR1035, miR1050, miR1088, and miR1126, were significantly down-regulated in response to drought stress (Table 1). Conversely, 14 miRNAs, namely miR159, miR169, miR171, miR319, miR395, miR474, miR845, miR851, miR854, miR896, miR901, miR903, miR1026, and miR1125, were significantly up-regulated under drought stress conditions.

Although some miRNA gene families, such as miR171, miR319, and miR896, were identified in both down- and up-regulated groups in this study, the specific miRNA member within each miRNA gene family was different. As shown in Table 1, smo-miR171b was up-regulated while the other four miR171 genes were down-regulated. However, in miR319 and miR896, miRNA genes were both down- and up-regulated depending on the developmental stage and organ considered. Taking gma-miR319c as an example, it was up-regulated on 12 DAW during the tillering stage and on 5 DAW during the inflorescence-forming stage, but it was down-regulated on 6 DAW. As a result, most miRNAs were identified with a distinct expression direction under drought stress while only the expression direction of miR171, miR319, and miR896 was found to overlap.

When the fold change criterion was taken as less than –1.0 or more than +1.0 (the same value on the log2 scale), a total of 74 significantly down-regulated miRNAs were identified during the tillering stage, of which seven miRNAs, belonging to the miR170, miR171, and miR396 families, were down-regulated on 10, 12, and 14 DAW during tillering (Fig. 2A). miR171 and miR396 were also down-regulated by drought stress, rather than up-regulated by abiotic stress as reported in Arabidopsis (Liu et al., 2008). Similarly, the most abundant down-regulated miRNAs during the inflorescence-forming stage were the miR166, miR172, miR896, and miR1030 gene families (Fig. 2B). Additionally, a total of 52 significantly up-regulated miRNAs were identified at the tillering stage, of which seven up-regulated miRNAs belonging to the miR474, miR854, and miR894 gene families were identified in all three of the water-withholding regimes (Fig. 2C). Among the drought-induced up-regulated miRNAs, miR854 was the only one to be up-regulated during both the tillering and inflorescence-forming stages (Fig. 2D). Combining the expression profiles during the two stages, miR172 and miR854 were found to be down- and up-regulated, respectively, throughout both the tillering and inflorescence-forming stages.

Fig. 2.

The Venn diagram illustrates the numbers of common and unique differential miRNAs induced by drought stress. (A and B) Down-regulated miRNAs identified with a 1-fold change in expression level at the tillering and inflorescence-forming stages, respectively, in rice. (C and D) Up-regulated miRNAs identified with a 1-fold change in expression level at the tillering and inflorescence-forming stages, respectively. (This figure is available in colour at JXB online.)

By BLAST analysis, it was found that eight of the 30 (26%) differential miRNAs were conserved in the rice genome. These are miR159, miR169, miR170, miR171, miR172, miR319, miR396, and miR397. Three other miRNAs, miR156, miR168, and miR408, were found in the rice genome with only 1–2 mismatches. Thus, 36% of the significantly detectable miRNAs match the rice genome. Considering that members of the same miRNA families differ from each other by only 1–3 nucleotides, cross-hybridization probably occurred in the microarray experiments between members of the same miRNA family, or between homologous sequences across species. However, hybridization in microarray experiments is highly strict and sensitive, and especially for such a short fragment like a miRNA the cross-hybridization rate should not be nearly 60%, indicating that it is not simple to make this conclusion, and the mechanism needs to be uncovered in further studies.

Confirmation of differentially expressed miRNAs by real-time qPCR

Real-time qPCR has been used for independent validation of the differentially expressed miRNAs determined by miRNA microarray in animal tissues (Lu et al., 2005). This method has been proven to be quantitative and sensitive, and is specific enough to discriminate a single nucleotide difference between miRNAs (Wang and Li, 2005). Among drought-stress induced down-regulated miRNAs, 22 selected miRNAs (detected by microarrays on 10, 12, and 14 DAW during the tillering stage and on 5 and 6 DAW during the inflorescence-forming stage) showed decreased expression intensities (<1.0 in the relative expression level calculated with the 2^–ΔΔCT value) (Fig. 3A). Similarly, 23 selected miRNAs showed an increased relative expression level (2^–ΔΔCT value >1.0) in Fig. 3B. Although the specific expression value by qRT-PCR was not exactly identical to the fold change calculated by the microarray results, the differential expression trend or direction was the same by both methods. The qRT-PCR results ultimately reflected consistency with the microarray data.

Fig. 3.

Confirmation by qRT-PCR of differentially expressed miRNAs obtained from miRNA array results. Actin was used as the reference level of miRNA expression. (A) Down-regulated miRNAs. (B) Up-regulated miRNAs. Plants at 10, 12, and 14 DAW were sampled at the tillering stage, and those at 5 and 6 DAW were sampled at the inflorescence-forming stage. (This figure is available in colour at JXB online.)

Target gene prediction and functional classification of drought stress-inducible miRNAs

The miRNA sequences were searched against rice genomic sequences on the plant miRNA potential target finder (http://bioinfo3.noble.org/miRNA/miRU.htm) to predict plant miRNA target genes. Out of the 165 differentially expressed miRNAs with fold change ≥1.5 or less than or equal to –1.5, 139 (84%) could be searched for their potential target genes. The number of potential targets in rice varies from two to 18 targets per miRNA family. Thus a score cut-off of ≤3 and mismatch <3 nucleotides was chosen to minimize the number of non-authentic targets. The main annotated targets and functional annotation of differential miRNAs with fold change exceeding ±2 are listed in Table 2. The results of the search indicated that the functions of the target genes of the drought stress-induced miRNAs identified are in keeping with previous reports, including the morphological differentiation and development of shoot organs such as root, leaf, and floral organs, hormone signal responses, miRNA regulation, and abiotic stress responses. Except for the deduced and experimentally confirmed functions of target genes, most predicted targets remain unknown regarding their functions in plant growth and development (Fig. 4).

Table 2.

Predicted target genes of differentially expressed miRNAs and their functional annotation

MiRNAsa	Identified direction of expression	Cluster	Target genes	Function
miR156	Down-regulated	I	cDNA teosinte glume architecture 1, cDNA squamosa promoter-binding-like protein 9	Evolution of maize from teosinte, flowering time, affecting sterility. NaCl stress response (Jones-Rhoades and Bartel, 2004; Sunkar and Zhu, 2004)
miR159	Up-regulated	I	MYB-like DNA-binding domain-containing protein, transcription factor GAMYB	Flowering time and flower organ development. NaCl stress response (Jones-Rhoades and Bartel, 2004; Sunkar and Zhu, 2004)
miR168	Down-regulated	II	ARGONAUTE1 protein (AGO1) family protein	Plant development. NaCl, and mannitol stress response (Sunkar et al., 2005; Gray-Mitsumune and Matton, 2006)
miR169	Up-regulated	I	CCAAT-binding transcription factor (CBF-B/NF-YA) family protein	Drought stress response (Jones-Rhoades and Bartel, 2004; Sunkar and Zhu, 2004)
miR170/miR171	Up/down-regulated	I	Scarecrow-like transcription factor, sodium-inducible calcium-binding protein (ACP1)	Floral development. NaCl, mannitol, and 4 °C stress response (Jones-Rhoades and Bartel, 2004; Sunkar and Zhu, 2004)
miR172	Down-regulated	I	cDNA floral homeotic protein APETALA2,bZIP ranscription factore family protein	Flowering time.floral organ identity. 4 °C stress response (Jones-Rhoades and Bartel, 2004; Sunkar and Zhu, 2004)
miR319	Up/down-regulated	I	MYB family transcription factor	Morphogenesis of shoot lateral organs. NaCl and 4 °C stress response (Jones-Rhoades and Bartel, 2004; Sunkar and Zhu, 2004)
miR395	Up-regulated	II	cDNA low affinity sulphate transporter	Response to sulphate nutrition stress in plant (Sunkar et al., 2005; Gray-Mitsumune and Matton, 2006)
miR396	Down-regulated	II	cDNA NB-ARC domain-containing protein,cysteine proteinase, root hair initiation protein (RHL1), disease resistance protein	NaCl, mannitol, and 4 °C stress responses; GRF transcription factors (Sunkar et al., 2005; Gray-Mitsumune and Matton, 2006)
miR397	Down-regulated	II	cDNA L-ascorbate oxidase precursor, monocopper oxidase-like protein SKS1 precursor	4 °C stress response (Sunkar et al., 2005; Gray-Mitsumune and Matton, 2006)
miR408	Down-regulated	II	cDNA chemocyanin precursor,cDNA phosphatidylinositol 3- and 4-kinase family protein	Pollen tube growth or guidance (Sunkar et al., 2005; Gray-Mitsumune and Matton, 2006)
miR474	Up-regulated	III	PPR, protein kinase, kinesin, leucine-rich repeat	Unknown
miR529	Down-regulated	I	cDNA squamosa promoter-binding-like protein 9, zinc finger family protein, U-box domain-containing protein	Possible involvment in flower development (Jones-Rhoades and Bartel, 2004; Sunkar and Zhu, 2004)
miR845	Up-regulated	III	cDNA retrotransposon protein	Unknown
miR851	Up-regulated	III	cDNA ATP-binding protein, putative, expressed	Unknown
miR854	Up-regulated	III	cDNA transferase, putative, expressed	Unknown
miR896	Up/down-regulated	III	cDNA transposon protein, putative, unclassified	Unknown
miR901	Up-regulated	III	cDNA phosphoglycerate kinase, putative, expressed	Unknown
miR903	Up-regulated	III	cDNA protein-binding protein, putative, expressed	Unknown
miR1026	Up-regulated	III	cDNA DNA-directed RNA polymerase 1B	Unknown
miR1030	Down-regulated	III	cDNA hypothetical protein	Unknown
miR1035	Down-regulated	III	cDNA integral membrane protein	Unknown
miR1050	Down-regulated	I	cDNA TATA-binding protein-associated factor 2N	Unknown
miR1088	Down-regulated	III	cDNA expressed protein, putative	Unknown
miR1125	Up-regulated	III	cDNA expressed protein	Unknown
miR1126	Down-regulated	III	cDNA expressed protein	Unknown

a

miRNAs with fold change over±2 and P-value <0.001.

Fig. 4.

Functional classification of differentially expressed miRNAs induced by drought stress at the tillering and inflorescence-forming differentiation stages. (A) Drought-induced up-regulated miRNAs at each sampling stage. (B) Drought-induced down-regulated miRNAs at each sampling stage. (This figure is available in colour at JXB online.)

According to the characteristics of predicted targets and their functions, the differential miRNAs with fold change >2 in down- and/or up-regulation in this study could be classified into three clusters. The first cluster includes eight miRNAs, namely miR156, miR159, miR169, miR171, miR172, miR319, miR529, and miR1050 (Table 2). Their target genes include different transcription factor families such as SBP, CBF, HAP2, ACP1, MYB, bZIP, and zinc finger, which further regulate gene expression and signal transduction that probably play roles in stress responses (Jones-R4hoades and Bartel, 2004; Sunkar and Zhu, 2004). Most of these eight miRNAs in the first cluster were down-regulated by drought stress in this study. Thus, it is presumed that the decrease of these miRNAs would lead to the increased expression of some corresponding transcription factors which, in turn, activate a set of coding genes to play defensive roles against abiotic stresses (experiments are being undertaken in rice and Arabidopsis plants). The second cluster includes five miRNAs, miR168, miR395, miR396, miR397, and miR408, which are directly involved in the process of stress response or stress tolerance (Sunkar et al., 2005; Gray-Mitsumune and Matton, 2006). The last cluster includes 13 miRNAs (miR474, miR845, miR851, miR854, miR896, miR901, miR903, miR1026, miR1030, miR1035, miR1088, miR1125, and miR1126) whose predicted target genes are basically putative expressed cDNAs of unknown function in plant growth (Table 2).

Except for the classification of miRNAs mentioned above, many miRNAs, such as miR156, miR159, miR168, miR171, miR172, and miR319, have been shown to be involved in the regulation of the morphological development process as well as in the stress response or stress tolerance process (Williams et al., 2005; Barakat et al., 2007; Liu et al., 2008). Furthermore, a few of the miRNAs, such as miR389 and miR400, reported to be induced by high salinity, drought, or cold stress in Arabidopsis (Sunkar and Zhu, 2004), were not identified in rice in the present study. The phosphate deficiency-responsive miR399 in Arabidopsis (Miura et al., 2005) was identified in this study, but the fold change was <1 in log2 value. Some miRNAs (e.g. miR171, miR319, and miR896) were identified to be differentially expressed in both down- and up-regulated directions at different sampling dates during the tillering and inflorescence-forming stages, which indicated the complexity of miRNA involvement in regulation during drought stress.

Pathway analysis of some drought stress-induced miRNAs

For the novel drought stress-induced miRNAs that were first reported in the present study, or whose expression profile was opposite to that in Arabidopsis, the search for target genes and bioinformatics analysis on metabolism was carried out through KEGG (http://www.genome.ad.jp/kegg/). The up-regulated miR395 and miR854 are involved in sulphur metabolism and ubiquinone biosynthesis, respectively. The coding protein β-fructofuranosidase, encoded by genes presumed through computational analysis to be target genes down-regulated by miR171, miR397, and miR408, takes part in starch and sucrose metabolism (Supplementary Fig. S1 available at JXB online). Down-regulated miR397 also plays a role in the reductive carboxylate cycle (CO₂ fixation). CO₂ fixation and starch synthesis are important biochemical processes in plant growth. Maintaining a reasonable rate of synthesis of carbon–hydrogen compounds helps to enhance stress tolerance under drought-prone conditions.

Analysis of cis-elements found in the promoter regions of drought induced miRNAs

Cis- and trans-elements in promoter regions have been extensively reported to be involved in expression of stress-induced genes (Jaglo-Ottosen et al., 1998; Kasuga et al., 1999; Zhang et al., 2005; Liu et al., 2008). The 1000 bp upstream promoter regions of 18 drought-induced osa-miRNA gene families were analysed to search for known stress-responsive cis-elements by searching the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) in order to reveal further precise promoter motifs and the relationship between a drought-induced miRNA, its target gene, and the relevant promoter (Lescot et al., 2002). Among the putative cis-elements shown in Table 3, several are known stress-responsive elements, such as the cis-acting element involved in abscisic acid (ABA) responsiveness (CE3), the ABA-responsive element (motif Iib), the TATA-box as the core element for start of transcription, and the GC motif involved in anoxic-specific induction. Some other regulatory elements, namely the TATC-box and the P-box for gibberellin responsiveness, Skn-1 and GCN4 for endosperm expression, the A-box for α-amylase promoters, previously reported to be stress-relevant in the analysis of regulatory sequences in the promoters of zinc finger proteins Zat12 and Zat7, and transcription factor MRKY25 (Rizhsky et al., 2004), were also identified in this study. Therefore, these elements can tentatively be related to the drought stress response or tolerance process. The G-box for the light stress response, although not clearly confirmed to have a role in drought stress responsiveness, was proven to be involved in the complex biochemical systems to perceive and respond to light of different wavelengths for photosynthetic organisms (Martínez-Hernández et al., 2002).

Table 3.

Types of drought stress-relevant cis-elements in the upstream regions of miRNAs in rice

cis-Element type	Element funtion	osa-miRNAs with cis-elements found in the 1000 bp upstream promoter regionsa
ATGCAAAT motif	cis-Acting regulatory element associated with the TGAGTCA motif	miR169(1)
P-box	Gibberellin-responsive element	miR156(1), miR159(1), miR160(1), miR162(1), miR166(1), miR171(1), miR399(1)
Skn-1 motif	cis-Acting regulatory element required for endosperm expression	miR156(1), miR159(5), miR160(1), miR164(2), miR166(2), miR169(1), miR171(2), miR172(4), miR319(2), miR396(2), miR397(1), miR398(2), miR399(1), miR530(1), miR531(1)
TATA-box	Core promoter element around –30 of the transcription start	miR164(1), miR167(4), miR169(1), miR397(2), miR398(1), miR806(1)
GC-motif	Enhancer-like element involved in anoxic specific inducibility	miR159(1), miR171(1), miR172(1), miR531(1)
GCN4-motif	cis-Regulatory element involved in endosperm expression	mi159(1), miR166(2), miR169(1), miR171(1), miR399(1)
TATC-box	cis-Acting element involved in gibberellin responsiveness	miR171(1)
A-box	Aequence conserved in α-amylase promoters	miR160(3), miR164(3)
TATCCAT/C-motif		miR160(3), miR164(1), miR167(1), miR397(2), miR399(1), miR806(1)
G-box	cis-Acting regulatory element involved in light responsiveness	miR164(1), miR169(1), miR806(1)
Motif Iib	Abscisic acid-responsive element	miR166(1)
CE3	cis-Acting element involved in ABA and VP1 responsiveness	miR399(1)

a

Identified differential osa-miRNAs in the microarray result; number in parentheses indicates replicates of the cis-element which is the average number of cis-elements from members of each miRNA family.

Among the promoters of the 18 miRNAs analysed in this study, 15 have the Skn-1 motif required for endosperm expression and seven have a P-box motif for the gibberellin-responsive element. The promoters of six miRNAs had a TATA-box, the core element for start of transcription. Based on the number of cis-element types found in the promoter region, miR169 and miR171 have five kinds of cis-elements in the promoter region, followed by miR159, miR160, miR164, miR166, miR399, and miR806 with fewer and fewer cis-elements, respectively (Table 3). Taken together, enrichment of stress response-relevant cis-elements could provide additional evidence that the above-mentioned eight rice miRNAs are very likely to be involved in the response to drought stress.

Discussion

In the past few years, many studies have been devoted to transcript regulation in response to high salinity, cold, drought, and ABA treatment in rice through cDNA microarray analysis (Kawasaki et al., 2001; Rabbani et al., 2003; Oh et al., 2005). As a result, a large number of abiotic stress-inducible genes have been identified, contributing to the elucidation of the molecular mechanism of the stress response. Expression products of induced genes were reported to function both in the initial stress response and in establishing stress tolerance in plants (Nakashima et al., 2007; Shinozaki and Yamaguchi-Shinozaki, 2007). MiRNAs are a large new class of small non-coding RNAs, which play important roles in post-transcriptional gene regulation by targeting mRNAs for cleavage or repressing translation. However, reports on microarray analysis with all plant miRNAs for drought-stressed rice grown in soil through tillering to the inflorescence-forming stage are not available so far. The microarray analysis of the miRNA expression profile under real drought stress in this study helped to better understand the function of plant miRNAs.

Here, samples from rice leaves at tillering and inflorescence-forming stages were used to identify miRNAs whose expression was significantly changed. Thirty-nine miRNA members with >2-fold down-regulation from 16 miRNA families, and 27 miRNA members with >2-fold up-regulation from 14 miRNA families were identified during the growing stages. The number of miRNAs involved in the drought stress response in rice in this study was larger than that identified by PEG-induced dehydration stress in rice seedlings, which identified only one miRNA (miR169g) (Zhao et al., 2007). The present results demonstrated that the differential expression of several miRNAs is involved in drought stress response or tolerance in rice.

Among the 30 miRNA families that were identified to be significantly down- or up-regulated under drought stress, 11 down-regulated families (miR170, miR172, miR397, miR408, miR529, miR896, miR1030, miR1035, miR1050, miR1088, and miR1126) and eight up-regulated families (miR395, miR474, miR845, miR851, miR854, miR901, miR903, and miR1125) were first revealed to be induced by drought stress in plants, and nine (miR156, miR168, miR170, miR171, miR172, miR319, miR396, miR397, and miR408) showed the opposite expression direction to that observed in drought-stressed Arabidopsis (Liu et al., 2008). The up-regulated miR169 under water-withholding drought stress in this study showed the same trend as in rice seedlings under PEG-monitored drought stress (Zhao et al., 2007). Regardless of the drought-induced expression directions of miRNAs, it could be inferred from the expression profiles that these 30 miRNAs were most probably involved in certain physiological or developmental processes in both Arabidopsis and rice plants under drought stress.

To explore the cross-talk among the differentially expressed miRNAs from the five sampling regimes at the tillering and inflorescence-forming stage, the most consistent down-regulated miRNAs with fold change less than –1.0 were miR170, the miR171 family, and miR396. The role of miR170 was totally unknown and the present study for the first time showed that it was inducible by drought stress in plants. MiR171 and miR396 were found to target transcription factors involved in regulation of gene expression and signal transduction (Sunkar and Zhu, 2004; Axtel et al., 2007; Yang et at., 2007), but there is still no direct experimental evidence such as overexpression or inhibited expression of these miRNAs to verify their roles in the abiotic stress response. Similarly, the most conserved up-regulated miRNAs were miR474 and miR854a. The expression result in Populus trichocarpa implied the potential role of ptc-miR474 in the post-transcriptional regulation of organelle gene expression and RNA processing (Lu et al., 2005). Array analysis in submerged maize roots indicated that ptc-miR474 was reduced at the early submergence phase and induced after 24 h of submergence (Zhang et al., 2008). The results of this study also demonstrated the same trend for this miRNA. However, the specific role of ptc-miR474 remains to be verified in future experiments. Bioinformatics analysis indicated that miR854 has multiple target sites in the 3'-untranslated region of oligouridylate-binding protein1b (UBP1b), a gene encoding a heterogeneous nuclear RNA-binding protein (hnRNP) involved in the regulation of pre-mRNA maturation at different levels, including pre-mRNA splicing (Lambermon et al., 2000). Further study indicated that plants and animals share miRNAs of the miR854 family, suggesting a common origin of these miRNAs as regulators of basal transcriptional mechanisms (Arteaga-Vázquez et al., 2006). Here miR854 was confirmed as a drought stress-induced up-regulated miRNA in rice plants, which was a new finding, but the mechanism of its role remains unknown.

Conclusion

Among the large number of plant miRNAs basically identified in silico, only for a very limited number of miRNAs has it been possible to elucidate their role in post-transcriptional gene regulation and in physiological processes. Here the value of microarray analysis to detect drought stress-regulated miRNAs in rice at tillering and inflorescence-forming stages has been shown. Thirty plant miRNA gene families were identified to be significantly down- or up-regulated by drought stress in rice. Especially important is the finding that 19 miRNAs (details as stated in the Results and Discussion) were first identified to be drought stress inducible, and nine miRNAs (miR156, miR168, miR170, miR171, miR172, miR319, miR396, miR397, and miR408) with an opposite expression regulation to that in drought-stressed Arabidopsis. These observations may facilitate uncovering of the mechanism of action of miRNAs in regulating plant response to abiotic stress, especially in water deficiency-induced stress. Other miRNAs, such as miR170, miR171, miR408, miR474, miR854, and miR1030, could be further investigated to study their roles in stress signalling and in the establishment of drought tolerance in plants. Target genes and cis-element analysis of the target gene promoter regions provide further evidence for the possible involvement of these differential miRNAs in the process of drought response and/or tolerance. The identification of these differentially expressed plant miRNAs would help to reveal the molecular mechanism of drought stress response and tolerance establishment in rice, and in particular provide an expression basis for further investigation of the potential role of individual miRNA under drought stress.

This work was supported by grants 05PJ14085 and 09DJ1400501 from the Science and Technology Committee of Shanghai Municipal Administration, the National Key Basic Research Program (grants 863-2006AA100101-4, 2007CB815701, 2009AAQ04002), the National Natural Science Foundation (grant 30830071) and grants 948-2006-R1 and 2009ZX08009-071B from the Ministry of Agriculture in China. The authors are grateful to Professor W. Zhang, Washington University, St Louis, for his advice and revision during preparation of the manuscript. The authors thank LC Sciences, Houston, Texas, USA for plant miRNA microarray analysis, and Shanghai SensiChip Tech & Infor Company for part of the data analysis. We specially thank Professor Roberto Tuberosa and Ms Kay Tuberosa for comments and revisions on an earlier version of the manuscript.

Author notes

†

These authors contributed equally to this work.