Accepted Manuscript Title: Construction of regulatory networks mediated by small RNAs responsive to abiotic stresses in rice (Oryza sativa) Author: Jingping Qin Xiaoxia Ma Zhonghai Tang Yijun Meng PII: DOI: Reference:

S1476-9271(15)30012-8 http://dx.doi.org/doi:10.1016/j.compbiolchem.2015.05.006 CBAC 6431

To appear in:

Computational Biology and Chemistry

Received date: Revised date: Accepted date:

20-10-2014 16-4-2015 22-5-2015

Please cite this article as: Qin, Jingping, Ma, Xiaoxia, Tang, Zhonghai, Meng, Yijun, Construction of regulatory networks mediated by small RNAs responsive to abiotic stresses in rice (Oryza sativa).Computational Biology and Chemistry http://dx.doi.org/10.1016/j.compbiolchem.2015.05.006 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Construction of regulatory networks mediated by small RNAs responsive to abiotic stresses in rice (Oryza sativa) Jingping Qin1,&, Xiaoxia Ma2,&, Zhonghai Tang1, *, Yijun Meng2, * 1 College of Bioscience and Biotechnology, Hunan Agricultural University, Changsha, Hunan 410128, China 2 College of Life and Environmental Sciences, Hangzhou Normal University, Hangzhou 310036, P. R. China

*

Corresponding authors:

Zhonghai Tang College of Bioscience and Biotechnology, Hunan Agricultural University, Changsha, Hunan 410128, China E-mail: [email protected]

Yijun Meng Xuelin Street 16#, Xiasha, Hangzhou 310036, P. R. China Tel: +86-571-28865198 E-mail: [email protected]

&These authors contributed equally to this work.

Graphical abstract fx1

Highlights 

Thousands of sRNAs responsive to abiotic stresses were identified in the rice.



The sRNAs enriched in Argonaute 1 were extracted for target identification.



Twelve sRNA—target lists were obtained for network construction. 1



Within certain subnetworks, some targets were supported by microarray data.



Literature mining indicated that some targets were involved in stress response.

Abstract Plants have evolved exquisite molecular mechanisms to adapt to diverse abiotic stresses. MicroRNAs play an important role in stress response in plants. However, whether the other small RNAs (sRNAs) possess stress-related roles remains elusive. In this study, thousands of sRNAs responsive to cold, drought and salt stresses were identified in rice seedlings and panicles by using high-throughput sequencing data. These sRNAs were classified into 12 categories, including “Panicle_Cold_Down”, “Panicle_Cold_Up”, “Panicle_Drought_Down”, “Panicle_Drought_Up”, “Panicle_Salt_Down”, “Panicle_Salt_Up”, “Seedling_Cold_Down”, “Seedling_Cold_Up”, “Seedling_Drought_Down”, “Seedling_Drought_Up”, “Seedling_Salt_Down” and “Seedling_Salt_Up”. The stress-responsive sRNAs enriched in Argonaute 1 were extracted for target prediction and degradome sequencing data-based validation, which enabled network construction. Within certain subnetworks, some target genes were further supported by microarray data. Literature mining indicated that certain targets were potentially involved in stress response. These results demonstrate that the established networks are biologically meaningful. We discovered that in some cases, one sRNA sequence could be assigned to two or more categories. Moreover, within certain target-centered subnetworks, one transcript was regulated by several stress-responsive sRNAs assigned to different categories. It implies that these subnetworks are potentially implicated in stress signal crosstalk. Together, our results could advance the current understanding of the biological role of plant sRNAs in stress signaling. 2

Key words: abiotic stress; Argonaute 1 (AGO1); degradome; high-throughput sequencing (HTS); network; rice small RNA (sRNA); signal interaction

Introduction Plants under abiotic stresses such as heat, cold, drought and salt will suffer from physiological and metabolic damages harmful to growth and development. As compensation, the plants have evolved exquisite regulatory mechanisms to response to specific stress signals, thus improving their tolerance (Mirouze and Paszkowski, 2011; Seki et al, 2007; Shinozaki et al, 2003; Xiong et al, 2002). These stress-related molecular mechanisms have been extensively studied in plants these years (Hirayama and Shinozaki, 2010; Zhu, 2002). MicroRNAs (miRNAs), one of the small RNA (sRNA) species of ~21 nt in length, were demonstrated to play an important role in stress signal transduction in plants (Sunkar et al, 2007). For examples, the nitrate-responsive miR393—AFB3 regulatory module is involved in controlling root system architecture of Arabidopsis (Arabidopsis thaliana) (Vidal et al, 2010). miR395 is inducible by sulphur starvation and is important for sulfate assimilation in Arabidopsis (Jones-Rhoades and Bartel, 2004; Kawashima et al, 2011; Kawashima et al, 2009). miR399 is involved in phosphate signaling in both Arabidopsis and rice (Bari et al, 2006; Zhou et al, 2008). With the widespread application of high-throughput sequencing (HTS) in sRNA biology, an unexpectedly huge universe of plant endogenous sRNAs, in addition to the miRNAs, has been uncovered. However, whether these largely unclassified sRNAs have stress-related biological roles remain unclear. To date, only a few cases have been reported. An impressive example was provided by Borsani et al. (2005). The nat-siRNAs [short interfering RNAs derived from natural antisense transcript (NAT) pairs] produced from a cis-NAT pair SRO5-P5CDH are involved in salinity stress response of Arabidopsis 3

(Borsani et al, 2005). A recent study on rice sRNAs generated tens of sRNA HTS data sets including those prepared from seedlings and panicles under specific abiotic stress treatments (drought, salt and cold) (Jeong et al, 2011). In that study, Jeong and his colleagues mainly focused on the rice miRNAs for their stress responsiveness and target genes. However, we recognized that these sRNA HTS data sets were invaluable resources to evaluate the potential roles of the other unclassified sRNA species in abiotic stress response and tolerance of rice. In this consideration, we set out to do a comprehensive search for the stress-responsive sRNAs in the rice seedlings and the panicles. As a result, thousands of stress-responsive sRNAs were identified, which were classified into 12 categories, “Panicle_Cold_Down” (the sRNAs down-regulated under the cold treatment in panicles), “Panicle_Cold_Up”, “Panicle_Drought_Down”, “Panicle_Drought_Up”, “Panicle_Salt_Down”, “Panicle_Salt_Up”, “Seedling_Cold_Down”, “Seedling_Cold_Up”, “Seedling_Drought_Down”, “Seedling_Drought_Up”, “Seedling_Salt_Down” and “Seedling_Salt_Up”. For target identification, only the stress-responsive sRNAs enriched in Argonaute 1 (AGO1)-associated silencing complexes were included. Based on the target lists supported by degradome sequencing data, regulatory networks mediated by specific stress-responsive sRNAs were constructed. Deeper investigation on certain subnetworks based on microarray data and literature mining enabled us to gain further supportive evidences for the potential involvement of the target genes in stress signaling. On the other hand, it indicates that the networks are biologically meaningful. We also discovered that for some sRNA sequences, they could be assigned to two or more of the 12 categories mentioned above. Moreover, within certain target-centered subnetworks, one transcript was simultaneously regulated by several stress-responsive sRNAs assigned to different categories. This kind of regulation sometimes exists in the panicles or the seedlings specifically. Based on these observations, we proposed that these subnetworks were potentially implicated in stress signal interactions, some of which might occur with organ- or developmental stage-specific patterns. Taken together, our results could advance the current 4

understanding of the biological functions of the sRNAs, in addition to the miRNAs, in stress response, stress tolerance and signal interactions in plants.

Results and Discussion Identification of the stress-responsive sRNAs in rice The sRNA HTS data sets were classified into two groups, i.e. the one prepared from rice seedlings and the one from panicles (see Table S1 for detailed information). Specifically, the “seedling” group includes four “control” data sets (GSM309691, SC1I, SC2I and SC3D), three data sets prepared under “drought” treatment (GSM309692, SD1I and SD2D), three data sets prepared under “salt” treatment (GSM309693, SNa1I and SNa2D), and two data sets prepared under “cold” treatment (SCd1I and SCd2D). For the “panicle” group, it contains three “control” data sets (PC1I, PC1C and PC2C), three data sets prepared under “drought” treatment (PD1I, PD1I and PD2C), three data sets prepared under “salt” treatment (PNa1I, PNa1C and PNa2C), and two data sets prepared under “cold” treatment (PCd1I and PCd2D) (Figure 1). Similar rules were applied to extract sRNAs responsive to a specific stress treatment in seedlings or panicles. For example, to identify sRNAs up-regulated under the drought treatment in the seedlings: (1) the sRNA should be present in at least one of the “drought” data sets of the seedlings (GSM309692, SD1I or SD2D); (2) its normalized accumulation level should be 3 RPM (reads per million; the normalization step is similar to that for the degradome sequencing data; see Materials and Methods for detail) or higher, and should be three times or more than the levels of the same sequence in all of the “control” data sets of the seedlings (GSM309691, SC1I, SC2I and SC3D). To identify sRNAs down-regulated under the drought treatment in the seedlings: (1) the sRNA should be present in at least one of the “control” data sets of the seedlings (GSM309691, SC1I, SC2I and SC3D); (2) its normalized level should be 3 RPM or higher, and should be three times or more than the levels of the same sequence in all of the “drought” data sets of the seedlings (Figure 1). As a result, in the seedlings, 84,394 sRNAs (including 42 miRNAs) and 40,993 5

sRNAs (including 11 miRNAs) were identified to be down- and up-regulated under the drought treatment, respectively. Under the salt treatment, 83,513 sRNAs (including 49 miRNAs) and 76,082 sRNAs (including five miRNAs) were down- and up-regulated, respectively. Under the cold treatment, 86,434 sRNAs (including 30 miRNAs) and 4,207 sRNAs (including seven miRNAs) were repressed and activated, respectively. In the panicles, 20,948 sRNAs (including 26 miRNAs) and 7,611 sRNAs (including 17 miRNAs) were down- and up-regulated under the drought treatment, respectively. Under the salt treatment, 26,889 sRNAs (including 55 miRNAs) and 13,522 sRNAs (including 19 miRNAs) were down- and up-regulated, respectively. Under the cold treatment, 30,682 sRNAs (including 63 miRNAs) and 21,014 sRNAs (including 26 miRNAs) were down- and up-regulated, respectively (Figure 1 and Data S1). The sequence variants of the miRNAs termed isomiRs were discovered in recent years (Guo et al, 2012; Li et al, 2012; Llorens et al, 2013). In animals, some of the isomiRs have been demonstrated to have biological roles (Cloonan et al, 2011). Here, based on the above sRNA lists, the stress-responsive accumulation patterns of the rice miRNAs along with their 3’ isomiRs (1 to 5 nt shorter or longer than the miRNAs, and could be perfectly mapped to the corresponding miRNA precursors) were analyzed together. The result showed that for certain miRNAs, their isomiRs displayed similar accumulation patterns (Table S2). For instances, both osa-miR1318-5p and osa-isomiR1318-5p (+1) (1 nt longer than its miRNA at 3’ end) were down-regulated under the cold treatment in panicles. osa-miR1320-3p and osa-isomiR1320-3p (-1) (1 nt shorter than its miRNA at 3’ end) were both down-regulated under the drought and the salt treatments in seedlings. However, contrary accumulation patterns were also observed between specific isomiRs and their miRNAs (Table S2). For examples, in seedlings, osa-miR437 was down-regulated under the drought treatment while osa-isomiR437 (+2) (2 nt longer at 3’ end) was up-regulated. More interesting patterns were observed for miR159. In seedlings, osa-miR159b and osa-isomiR159b (+1) were down-regulated under cold treatment while osa-isomiR159b (-2) and osa-isomiR159b (-3) were up-regulated. Taken 6

together, the stress-responsive accumulation patterns between the miRNAs and their isomiRs are quite complex. Some of the plant sRNAs, such as the miRNAs, regulate their targets through cleavages (Chen, 2009; Jones-Rhoades et al, 2006; Voinnet, 2009). To perform cleavages, the sRNAs should incorporate into the RNA-induced silencing complexes (RISCs) associated with specific AGO proteins. It was demonstrated that the AGO1 protein possesses RNA slicer activity, and could selectively recruit miRNAs and certain short interfering RNAs in plants (Baumberger and Baulcombe, 2005; Vaucheret, 2008). Based on this notion, we set out to identify the AGO1-enriched sRNAs from the stress-responsive sRNAs. The following rules were adopted: (1) the sRNAs should be present in at least one of the AGO1-associated sRNA HTS data sets (GSM455962, GSM455963 or GSM455964); (2) its normalized abundance should be 3 RPM or higher, and should be three times or more than the level of the same sequence in the control set (GSM455965). As a result, in seedlings, 128 sRNAs (including osa-miR171b/c-3p/d-3p/e-3p/f-3p and osa-miR810b.1), and 122 sRNAs (including osa-miR396a-5p/b-5p and osa-miR166i-3p/j-3p) were identified to be down- and up-regulated under the drought treatment, respectively. Under the salt treatment, 158 sRNAs (including osa-miR319a-3p.2-3p/b, osa-miR394 and osa-miR171b/c-3p/d-3p/e-3p/f-3p), and 247 sRNAs were down- and up-regulated, respectively. Under the cold treatment, 192 sRNAs (including osa-miR319a-3p.2-3p/b), and 58 sRNAs were repressed and activated, respectively. In the panicles, 63 sRNAs (including osa-miR394, osa-miR1433, osa-miR171e-5p and osa-miR171b/c-3p/d-3p/e-3p/f-3p) and 60 sRNAs (including osa-miR396h/i/g, osa-miR396a-5p/b-5p, osa-miR396c-5p and osa-miR396d/e-5p) were down- and up-regulated under the drought treatment, respectively. Under the salt treatment, 107 sRNAs (including 13 miRNAs) and 79 sRNAs (including osa-miR399d) were downand up-regulated, respectively. Under the cold treatment, 135 sRNAs (including nine miRNAs) and 88 sRNAs (including osa-miR169b/c and osa-miR399j) were downand up-regulated, respectively (Figure 1 and Data S2). The sequence characteristics of the stress-responsive sRNAs and of those enriched 7

in the AGO1 complexes were analyzed. For the stress-responsive miRNAs, 5’ Us (uridines) occupy a large portion (~40% on average) except for the miRNAs up-regulated under the salt treatment in the seedlings (Figure 2A). It is quite different for the other stress-responsive sRNAs. The sRNAs starting with 5’ As (adenines) and the ones starting with 5’ Gs (guanines) together occupy ~60% of the stress-responsive sRNAs excluding miRNAs (Figure 2B). Besides, the 21-nt miRNAs occupy more than 40% of the stress-responsive miRNAs (Figure 2C). However, for the other stress-responsive sRNAs, the 24-nt ones occupy a dominant portion, especially for those identified in the panicles (~60%) (Figure 2D). After AGO1 enrichment analysis, the sequence characteristics were extensively altered. For most of the stress-responsive miRNAs, they start with 5’ Us, and are 20 or 21 nt in length (Figure S1A and Figure S1C). For the other stress-responsive sRNAs, ~55% start with 5’ Us, and ~96% are 19 to 21 nt in length (Figure S1B and Figure S1D). A previous study reported that the sRNAs associated with the AGO1 silencing complexes favored 5’ Us and were predominantly 21 nt in length (Mi et al, 2008). From this point of view, our AGO1 enrichment filtering tended to selectively extract 5’ U-started, 19 to 21 nt ones from the stress-responsive sRNAs.

Identification of the targets of the AGO1-enriched, stress-responsive sRNAs, and network construction We speculated that the AGO1-enriched, stress-responsive sRNAs possessed great potential of performing target cleavages. In this regard, target prediction was performed for these sRNAs by using miRU algorithm (Dai and Zhao, 2011; Zhang, 2005) with default parameters. As a result, under drought treatment, 2,542 transcripts were predicted to be targeted by 128 sRNAs down-regulated in seedlings, 3,726 transcripts were predicted to be targeted by 118 sRNAs up-regulated in seedlings, 1,409 transcripts were predicted to be targeted by 62 sRNAs down-regulated in panicles, and 1,288 transcripts were predicted to be targeted by 59 sRNAs up-regulated in panicles. Under salt treatment, 3,386 transcripts were predicted to be 8

targeted by 157 sRNAs down-regulated in seedlings, 4,946 transcripts were predicted to be targeted by 245 sRNAs up-regulated in seedlings, 1,830 transcripts were predicted to be targeted by 105 sRNAs down-regulated in panicles, and 1,160 transcripts were predicted to be targeted by 79 sRNAs up-regulated in panicles. Under the cold treatment, 4,352 transcripts were predicted to be targeted by 191 sRNAs down-regulated in seedlings, 1,913 transcripts were predicted to be targeted by 56 sRNAs up-regulated in seedlings, 2,437 transcripts were predicted to be targeted by 133 sRNAs down-regulated in panicles, and 1,391 transcripts were predicted to be targeted by 86 sRNAs up-regulated in panicles. Degradome sequencing is a high-throughput method for sRNA—target pair validation (German et al, 2009; German et al, 2008). According to the analytical workflow (Meng et al, 2011), four degradome sequencing data sets (Table S1) were included for target validation. Considering the previous observations that some evident cleavage signals resided out of the canonical slicing sites (10th to 11th nt of the regulatory sRNAs) (Addo-Quaye et al, 2008; Addo-Quaye et al, 2009; Li et al, 2010; Llave et al, 2002), the binding sites with prominent cleavage signals resided within 8th to 12th nt of the regulatory sRNAs were retained. Among the validated sRNA—target pairs, we observed that many binding sites were supported by evident degradome signatures resided within the canonical cleavage sites (10th to 11th nt of the regulatory sRNAs) (Figure 3 and Figure S2). Based on the target lists, the sRNA-mediated regulatory networks were constructed by using Cytoscape (Shannon et al, 2003) (Figure S3). According to the stress responsiveness of the sRNAs, the networks were classified into 12 groups, i.e. “Panicle_Cold_Down” [a total of 158 sRNA--target interactions exist in this network which involves 105 target transcripts and 29 sRNAs (including 14 miRNAs) down-regulated under cold treatment in panicles], “Panicle_Cold_Up” [203 interactions exist in this network which involves 70 targets and 21 sRNAs (including osa-miR169b and osa-miR169c) up-regulated under cold treatment in panicles], “Panicle_Drought_Down” [47 interactions involving 46 targets and 12 sRNAs (including osa-miR171b, osa-miR171c-3p, osa-miR171d-3p, osa-miR171e-3p, 9

osa-miR171f-3p and osa-miR394) down-regulated under drought treatment in panicles], “Panicle_Drought_Up” [191 interactions involving 60 targets and 23 sRNAs (including osa-miR396a-5p, osa-miR396b-5p, osa-miR396c-5p, osa-miR396d, osa-miR396e-5p, osa-miR396g, osa-miR396h and osa-miR396i) up-regulated under drought treatment in panicles], “Panicle_Salt_Down” [146 interactions involving 75 targets and 32 sRNAs (including 19 miRNAs) down-regulated under salt treatment in panicles], “Panicle_Salt_Up” [65 interactions involving 25 targets and 11 sRNAs up-regulated under salt treatment in panicles], “Seedling_Cold_Down” [110 interactions involving 81 targets and 20 sRNAs (including osa-miR319a-3p.2-3p and osa-miR319b) down-regulated under cold treatment in seedlings], “Seedling_Cold_Up” [39 interactions involving 20 targets and nine sRNAs up-regulated under cold treatment in seedlings], “Seedling_Drought_Down” [92 interactions involving 57 targets and 17 sRNAs (including osa-miR171b, osa-miR171c-3p, osa-miR171d-3p, osa-miR171e-3p and osa-miR171f-3p) down-regulated under drought treatment in seedlings], “Seedling_Drought_Up” [60 interactions involving 59 targets and 11 sRNAs (including osa-miR166i-3p, osa-miR166j-3p, osa-miR396a-5p and osa-miR396b-5p) up-regulated under drought treatment in seedlings], “Seedling_Salt_Down” [160 interactions involving 129 targets and 27 sRNAs (including osa-miR171b, osa-miR171c-3p, osa-miR171d-3p, osa-miR171e-3p, osa-miR171f-3p, osa-miR319a-3p.2-3p, osa-miR319b and osa-miR394) down-regulated under salt treatment in seedlings] and “Seedling_Salt_Up” [102 interactions involving 89 targets and 23 sRNAs up-regulated under salt treatment in seedlings]. The result indicates that in addition to the miRNAs, some of the other unclassified sRNAs could guide the AGO1 complexes to perform target cleavages under specific stress treatments in rice. Besides, these degradome sequencing data-supported networks could serve as a basis for in-depth studies on the sRNA-mediated stress signaling in rice.

Target expression- and literature mining-based analysis 10

Considering the repressive role of the AGO1-enriched sRNAs in gene expression, we deduced that the target genes should display contrary stress-responsive expression patterns compared to the corresponding sRNA regulators. To support this notion, two sets of the rice gene microarray data (OS10 and OS37) were obtained from PLEXdb (Plant Expression Database) (Dash et al, 2012). OS10, a gift from the previous report by Jain et al. (2007), contains the gene expression data of rice seedlings under drought, salt and cold treatments. OS37, provided by a recent study (Wang et al, 2011), contains the gene expression data under drought treatment in panicles. Notably, for both data sets, there are three replicates for each stress treatment. In this expression-based analysis, the mean value of the three replicates was calculated for each target gene. Several pieces of supportive evidences were obtained for certain sRNA—target pairs (Figure 4 and Data S3). Specifically, within the network mediated by the sRNAs up-regulated under drought treatment in panicles (Figure 4A), the expression level of LOC_Os02g53690 (regulated by miR396 and Panicle_Drought_Up_sRNA5602) in the panicles was 51.54 under drought treatment versus 91.25 under control. LOC_Os12g41860 regulated by Panicle_Drought_Up_sRNA5361, Panicle_Drought_Up_sRNA5364, Panicle_Drought_Up_sRNA7099 and Panicle_Drought_Up_sRNA7393 was expressed at 176.02 under drought treatment in panicles, while expressed at 467.96 under control. LOC_Os10g33960, also regulated by the above four sRNAs, was expressed at 288.64 under drought treatment in panicles versus 628.12 under control. Similar supportive evidences were found for the target genes LOC_Os01g48060, LOC_Os01g54990, LOC_Os03g51330, LOC_Os03g51970, LOC_Os07g32170, LOC_Os12g29980 and LOC_Os12g41950, which were expressed much lower under drought treatment than control in panicles (Data S3). For the network mediated by the sRNAs down-regulated under the drought treatment in the panicles (Figure 4B), LOC_Os05g03574 and LOC_Os05g51140 both regulated by Panicle_Drought_Down_sRNA17639 were expressed at 80.20 and 2258.63 under drought treatment versus 48.57 and 1460.65 under control, respectively. Similar supportive evidences were found for LOC_Os01g69940 and 11

LOC_Os08g19114 (Data S3). In seedlings, within the network mediated by the sRNAs repressed by drought stress (Figure 4C), LOC_Os05g48590 targeted by Seedling_Drought_Down_sRNA51230 was expressed at 308.67 under drought treatment in seedlings versus 102.80 under control (Data S3). For LOC_Os02g17940 regulated by Seedling_Drought_Up_sRNA38390 (Figure 4D), it was expressed at 23.70 under drought and at 59.76 under control (Data S3). Within the network mediated by the sRNAs down-regulated under salt treatment in seedlings (Figure 4E), LOC_Os02g39910 and LOC_Os06g44620, both targeted by Seedling_Salt_Down_sRNA46587, were expressed at 105.13 and 1235.80 under salt treatment while at 55.55 and 805.51 under control, respectively (Data S3). Within the network mediated by the sRNAs up-regulated under salt treatment in seedlings (Figure 4G), LOC_Os02g17940 targeted by Seedling_Drought_Up_sRNA38390 and Seedling_Salt_Up_sRNA72907 was expressed at 30.94 under salt treatment while at 59.76 under control (Data S3). Within the network mediated by the sRNAs down-regulated under cold treatment in seedlings (Figure 4F), LOC_Os01g59660, simultaneously targeted by Seedling_Cold_Down_sRNA37423, Seedling_Cold_Down_sRNA56879 and Seedling_Cold_Down_sRNA57515, was expressed at 173.70 under cold treatment versus 71.58 under control (Data S3). LOC_Os04g46860 regulated by Seedling_Cold_Down_sRNA54170, was expressed at 307.84 under cold treatment in seedlings versus 173.53 under control (Data S3). More microarray data-based supportive evidences could be found in Data S3. Summarily, contrary expression patterns under specific stress treatments were observed between certain target genes and the regulatory sRNAs, which was consistent with the repressive role of the AGO1-enriched sRNAs in gene expression. We intended to find some stress-related evidences for the target genes based on the functional annotations and literature mining. Fortunately, some positive hints were obtained. Within the network mediated by the sRNAs down-regulated under cold treatment in seedlings (Figure 4F), LOC_Os01g59660 targeted by Seedling_Cold_Down_sRNA37423, Seedling_Cold_Down_sRNA56879 and Seedling_Cold_Down_sRNA57515 encodes an MYB family transcription factor 12

based on the annotation from TIGR rice (the rice genome annotation project established by the institute for genome research) (Yuan et al, 2003). Involvement of the rice MYB genes in cold tolerance has been demonstrated by transgenic experiments in both Arabidopsis and rice (Dai et al, 2007; Park et al, 2010; Vannini et al, 2004; Yang et al, 2012). Both LOC_Os03g29760 and LOC_Os12g42400, regulated by Seedling_Cold_Down_sRNA26378 and Seedling_Cold_Down_sRNA51152, encode the subunits of nuclear transcription factor Y, which are homologous to the subunit A of nuclear factor Y (NF-YA) in Arabidopsis. Based on a recent report, overexpression of NF-YA confers tolerance to diverse abiotic stresses including cold in Arabidopsis (Leyva-Gonzalez et al, 2012). More interestingly, the Aux/IAA family is down-regulated under the cold treatment in Arabidopsis (Hannah et al, 2005), whereas LOC_Os05g48590 encoding OsIAA19 was induced by the cold stress signal in rice based on the microarray data (Data S3) and a recent report (Jain and Khurana, 2009). Whether different members of the Aux/IAA family possess opposite roles in cold stress response needs further characterization. Based on a previous study, many auxin-related genes [including ARF (Auxin Response Factor) and Aux/IAA genes] were responsive to drought and salt stress signals in rice (Jain et al, 2009). This correlates well with the existence of certain ARF genes (LOC_Os01g48060, LOC_Os01g54990, LOC_Os12g41950) and Aux/IAA genes (LOC_Os05g48590) within the networks mediated by the sRNAs responsive to drought and salt stresses (Figure 4A, 4C and 4E). Within the networks mediated by the sRNAs responsive to drought stress (Figure 4A to 4D), several literature-based evidences were obtained. LOC_Os01g69940 regulated by osa-miR394 and Panicle_Drought_Down_sRNA4816 encodes OsFBX32, an F-box domain containing protein. According to the annotation from TIGR rice, the homologous gene of LOC_Os01g69940 in Arabidopsis is AT1G27340. Based on a recent report, overexpression of gma-miR394a of soybean increased drought tolerance of the transgenic Arabidopsis. Notably, AT1G27340 was the target of gma-miR394a (Ni et al, 2012). Consistently, overexpression of LOC_Os02g44990, another F-box gene, reduced the drought tolerance in transgenic rice (Yan et al, 2011). LOC_Os05g34700, 13

targeted by Seedling_Drought_Up_sRNA38390, encodes a GDSL-like lipase. Interestingly, overexpression of CaGLIP1 (a pepper GDLS-type lipase gene) in Arabidopsis enhanced the drought tolerance of the transgenic plants (Hong et al, 2008). LOC_Os03g51330, encoding a GRAS family transcription factor domain containing protein, is targeted by Panicle_Drought_Up_sRNA1813. Based on the recent reports, overexpression of BnLAS (a GRAS family gene of Brassica napus) or PeSCL7 (a poplar GRAS gene) could increase drought tolerance of the transgenic plants of Arabidopsis (Ma et al, 2010; Yang et al, 2011). Three genes (LOC_Os02g53690, LOC_Os03g51970 and LOC_Os12g29980, all of which were targeted by osa-miR396 and Panicle_Drought_Up_sRNA5602) encoding growth regulating factor (GRF) proteins were included in the drought-related networks (Figure 4A and 4D). In Arabidopsis, overexpression of miR396 targeting GRF genes increased drought tolerance of the transgenic plants (Liu et al, 2009). Both LOC_Os03g48970 (regulated by Seedling_Drought_Down_sRNA47087) and LOC_Os07g41720 (regulated by Seedling_Drought_Down_sRNA24360 and Seedling_Drought_Down_sRNA47087) encode the nuclear transcription factor Y subunits. In both Arabidopsis and maize, the nuclear factor Y family genes were demonstrated to be implicated in drought resistance (Leyva-Gonzalez et al, 2012; Li et al, 2008; Nelson et al, 2007). LOC_Os08g32500, targeted by Seedling_Drought_Up_sRNA31040, encodes a nucleobase-ascorbate transporter. In a recent study, the involvement of a rice nucleobase-ascorbate transporter (encoded by LOC_Os09g15170) in drought signal transduction was reported (Moumeni et al, 2011). LOC_Os07g32170 encoding OsSPL13 (a SBP-box family member) was regulated by five sRNAs up-regulated under the drought treatment in the rice panicles (Figure 4A). Potential role of SPL genes in drought tolerance was proposed in Medicago truncatula (Trindade et al, 2010). LOC_Os10g33960 (regulated by Panicle_Drought_Up_sRNA5361, Panicle_Drought_Up_sRNA5364, Panicle_Drought_Up_sRNA7099 and Panicle_Drought_Up_sRNA7393) and LOC_Os12g41860 (regulated by osa-miR166i-3p/j-3p, Panicle_Drought_Up_sRNA5361, Panicle_Drought_Up_sRNA5364, 14

Panicle_Drought_Up_sRNA7099 and Panicle_Drought_Up_sRNA7393) encode START domain containing proteins. In Arabidopsis, overexpression of an HD-START protein could improve drought tolerance (Yu et al, 2008). Within the networks mediated by the sRNAs responsive to salt stress in the seedlings (Figure 4E and 4G), LOC_Os01g04860 encoding a DUF647 domain containing protein was previously demonstrated to be highly responsive to salt stress (Pandit et al., 2011). LOC_Os05g34700 encoding a GDSL-like lipase was regulated by Seedling_Salt_Up_sRNA72907. According to a previous study, overexpression of a GDSL-motif lipase of Arabidopsis could increase salt tolerance in yeast and the transgenic plants (Naranjo et al, 2006). LOC_Os06g06560 encoding a starch synthase was regulated by Seedling_Salt_Down_sRNA53510. A previous study by Chen et al. (2008) showed that the starch content was decreased in the rice seedling leaves upon the NaCl-stressed treatment. And, this process was mediated by granule-bound starch synthase. LOC_Os09g34250 regulated by Seedling_Salt_Down_sRNA53510 encodes a UDP-glucosyltransferase. In Arabidopsis, the UDP-Glucosyltransferase UGT74E2 could improve water stress tolerance through perturbation of indole-3-butyric acid homeostasis (Tognetti et al, 2010). Summarily, all these evidences indicate that the networks established in this study are biologically meaningful. In recent years, many stress-responsive sRNAs were successfully identified from rice and other plants. To gain better understanding of the common biological roles of plant sRNAs in stress signaling, a comprehensive comparison was performed between the rice stress-responsive sRNAs discovered in this study and the previously identified stress-responsive plant sRNAs. In rice, it has been reported that miR166 and miR396 were down-regulated under drought stress (Ding et al., ; Zhou et al.), while miR169 and miR393 were up-regulated under water deficit (Contreras-Cubas et al., ; Ding et al., ; Zhou et al., ; Zhao et al., 2007). Consistently, in our study, we discovered that most of the miR166 family members (including osa-miR166a-5p, osa-miR166b-5p, osa-miR166d-5p, osa-miR166e-5p, osa-miR166g-5p and osa-miR166n-5p) and osa-miR396c-3p were greatly repressed by drought treatment in either seedlings or panicles. Moreover, according to the results of this study, osa-miR169a, 15

osa-miR169f.2, osa-miR169n, osa-miR169o and osa-miR393b-3p were identified to be up-regulated by drought stress. However, contradictory cases were also discovered in our analysis. For examples, miR168 was reported to be repressed by drought stress (Ding et al., ; Zhou et al.). But, we found that only one member of rice miR168 family (osa-miR168a-5p) was responsive to drought, and this miRNA was up-regulated by drought stress in panicles. Interestingly, in Arabidopsis, miR168 was demonstrated to be induced by drought stress (Ding et al., ; Liu et al., 2008). For cetain miRNA families, it has not been reported that these miRNAs were responsive to abiotic stresses in rice, but they were validated to be stress-responsive in other plant species according to previous studies. For instances, miR171 of Arabidopsis is induced by drought stress (Contreras-Cubas et al., ; Liu et al., 2008). In our study, we found that many miR171 family members (osa-miR171b, osa-miR171c-3p, osa-miR171d-3p, osa-miR171e-5p, osa-miR171e-3p, osa-miR171f-5p and osa-miR171f-3p) were also drought-responsive in rice, but they were repressed by water deficit. In common bean (Phaseolus vulgaris), miR2118 is up-regulated under water deficit (Contreras-Cubas et al., ; Arenas-Huertero et al., 2009). Notably, we discovered that osa-miR2118f, osa-miR2118h, osa-miR2118j, osa-miR2118k and osa-miR2118m were also inducible by drought stress in rice panicles. In Brassica juncea, the accumulation level of miR172 is repressed under salt or drought stress (Bhardwaj et al.). Interestingly, in rice seedlings, we found that miR172a, miR172b and miR172d-3p were also down-regulated under both stresses. In tomato (Solanum habrochaites), miR159 and miR390 were reported to be repressed by chilling stress, while miR399 was up-regulated by cold stress (Cao et al.). Consistently, in rice panicles, we identified that osa-miR159a.1, osa-miR159b and osa-miR390-5p were down-regulated by low-temperature treatment, while osa-miR399j was induced by chilling stress. Furthermore, we noticed that except for the miRNAs, all of the other rice stress-responsive sRNAs identified in our study had not been reported in any previous report. Thus, for the future experimental studies, it is interesting to see whether these stress-responsive sRNAs are conserved among different plant species or they are just rice-specific, and the biological roles of these sRNAs need to be analyzed in-depth. 16

Potential stress signal interactions within the networks We noticed that certain target genes supported by three pieces of evidences (i.e. degradome sequencing data, microarray data and literatures) were simultaneously regulated by several sRNAs responsive to different stresses. For example, LOC_Os05g34700 was regulated by Seedling_Drought_Up_sRNA38390 and Seedling_Salt_Up_sRNA72907. LOC_Os05g48590 was targeted by Seedling_Cold_Down_sRNA55719, Seedling_Drought_Down_sRNA51230 and Seedling_Salt_Down_sRNA51251. This interesting observation points to the potential interactions among different stress signals in rice. Based on this consideration, we first searched for the sRNA sequences that could be assigned to two or more of the 12 categories (“Panicle_Cold_Down”, “Panicle_Cold_Up”, “Panicle_Drought_Down”, “Panicle_Drought_Up”, “Panicle_Salt_Down”, “Panicle_Salt_Up”, “Seedling_Cold_Down”, “Seedling_Cold_Up”, “Seedling_Drought_Down”, “Seedling_Drought_Up”, “Seedling_Salt_Down” and “Seedling_Salt_Up”). All the stress-responsive sRNAs identified above (Data S1) were analyzed. As a result, thousands of sRNA sequences were identified to be responsive to multiple stress signals in the seedlings and/or the panicles (Data S4). More detailed statistical result is shown in Figure 5, and only the sRNA sequences simultaneously assigned to four or more categories were included. Specifically, 386 sRNA sequences could be simultaneously assigned to six categories, and they were consistently down-regulated under different stress treatments in both seedlings and panicles. Similarly, 164 sRNAs were assigned to “Panicle_Cold_Down”, “Panicle_Salt_Down”, “Seedling_Cold_Down”, “Seedling_Drought_Down” and “Seedling_Salt_Down”, 112 sRNAs were assigned to “Panicle_Cold_Down”, “Panicle_Drought_Down”, “Seedling_Cold_Down”, “Seedling_Drought_Down” and “Seedling_Salt_Down”, 197 sRNAs were assigned to “Panicle_Cold_Down”, “Seedling_Cold_Down”, “Seedling_Drought_Down” and “Seedling_Salt_Down”, and 101 sRNAs were assigned to “Panicle_Cold_Down”, “Panicle_Salt_Down”, “Seedling_Cold_Down” 17

and “Seedling_Salt_Down”. These results indicate that the stress responsiveness of certain sRNAs is not significantly altered between the vegetative seedlings and the reproductive panicles. Interestingly, these sRNAs were all repressed under the stress treatments. However, whether these sRNAs play an essential role in maintaining consistent stress responsiveness during developmental stage transition remain to be clarified. Contrary response patterns of a sRNA sequence between the seedlings and the panicles were also discovered. There are 160 sRNA sequences simultaneously assigned to “Panicle_Cold_Up”, “Panicle_Salt_Up”, “Seedling_Cold_Down”, “Seedling_Drought_Down” and “Seedling_Salt_Down”, 410 sRNA sequences assigned to “Panicle_Cold_Up”, “Seedling_Cold_Down”, “Seedling_Drought_Down” and “Seedling_Salt_Down”, and 520 sRNA sequences including osa-miR435 and osa-miR1425-3p assigned to “Panicle_Salt_Up”, “Seedling_Cold_Down”, “Seedling_Drought_Down” and “Seedling_Salt_Down”. It indicates that distinct response patterns under specific stress treatment exist between the vegetative seedlings and the reproductive panicles, which might be mediated by these novel stress-responsive sRNAs. Next, a comprehensive search was performed to extract the subnetworks simultaneously mediated by different categories of the stress-responsive sRNAs from the established networks (Figure S3). As a result, 192 transcription forms belonging to 133 target genes were regulated by the stress-responsive sRNAs assigned to two or more categories (Table S3). As discovered above, in many cases, one sRNA sequence was able to response to multiple stress signals. It includes eight miRNA sequences osa-miR159a.1/b, osa-miR160f-5p, osa-miR171b/c-3p/d-3p/e-3p/f-3p, osa-miR319a-3p.2-3p/b, osa-miR394, osa-miR396a-5p/b-5p, osa-miR396c-5p and osa-miR396d/e-5p (Table S4). To gain deeper insights, certain featured, target gene-centered subnetworks potentially implicated in stress signal interactions were extracted (Figure 6). Potential signal crosstalk was observed for the LOC_Os01g11430.1 (encodes an expressed protein)-centered (Figure 6A), the LOC_Os02g44360.1 (a scarecrow transcription factor family protein)-centered (Figure 6D), the LOC_Os06g40330.1 (a MYB family transcription factor)-centered 18

(Figure 6E), and the LOC_Os02g45570.1 (a growth regulating factor protein)-centered (Figure 6H) subnetworks. Each of these transcripts was targeted by the cold-responsive, the drought-responsive and the salt-responsive sRNAs simultaneously. We also noticed that for certain subnetworks, potential stress signal interactions appeared to be seedling-specific or panicle-specific (Table S3 and Figure 6). For examples, within the LOC_Os01g69940.1 (OsFBX32, an F-box domain containing protein)-centered (Figure 6B), the LOC_Os02g41800.1 (an auxin response factor)-centered (Figure 6F), the LOC_Os04g43910.1 (an auxin response factor)-centered (Figure 6I), and the LOC_Os10g39970.1 (a harpin-induced protein 1 domain containing protein)-centered (Figure 6J) subnetworks, signal interactions mediated by the stress-responsive sRNAs specifically exist in the panicles. Whether stress signal crosstalk could be mediated by these sRNAs, and whether certain signal interactions indeed display panicle- or seedling-specific patterns need further validations.

Materials and Methods Data sources High-throughput sequencing data sets of rice sRNAs were retrieved from GEO (Gene Expression Omnibus; http://www.ncbi.nlm.nih.gov/geo/) (Barrett et al, 2009) and the Rice Small RNA Database 2 belonging to the Next-Gen Sequence Databases (http://mpss.udel.edu/rice_sRNA2/index.php?menu=../common/web/library_info.php ?SITE=rice_sRNA2&showAll=true) (Nakano et al, 2006). The AGO1-associated sRNA HTS data sets and the rice degradome sequencing data sets were downloaded from GEO. Please refer to Table S1 for the detailed information of these data sets. The microarray data sets of rice were downloaded from PLEXdb (Plant Expression Database; http://www.plexdb.org/plex.php?database=Rice) (Dash et al, 2012). Specifically, OS10 contains the expression data under specific stress treatments (control, drought, salt and cold) of rice seedlings. Each treatment has three replicates, i.e. OS10_H1, OS10_H2 and OS10_H3 for control, OS10_H4, OS10_H5 and 19

OS10_H6 for drought treatment, OS10_H7, OS10_H8 and OS10_H9 for salt treatment, and OS10_H10, OS10_H11 and OS10_H12 for cold treatment. OS37 contains the expression data under the drought treatment of rice panicles. OS37_H16, OS37_H17 and OS37_H18 are the three replicates of the control sets of panicles, and OS37_H34, OS37_H35 and OS37_H36 are the three replicates of panicles under drought treatment. For both OS10 and OS37, the mean expression value of the three replicates of a specific treatment was calculated for the Affymetrix 57k Rice GeneChip probe of a specific target gene to perform expression data-based analysis. The miRNA sequences were obtained from miRBase (release 21; http://www.mirbase.org/) (Griffiths-Jones et al, 2008). The transcripts of rice genes and the gene annotations were retrieved from the rice genome annotation project established by the institute for genome research (TIGR rice, release 6.1; http://rice.plantbiology.msu.edu/index.shtml) (Yuan et al, 2003).

Prediction and validation of the sRNA/miRNA targets Target prediction was performed by using miRU algorithm (Dai et al, 2011; Zhang, 2005) with default parameters. The degradome sequencing data were utilized to validate the predicted sRNA—target pairs. First, in order to allow cross-library comparison, the normalized read count (in RPM, reads per million) of a short sequence from a specific degradome library was calculated by dividing the raw count of this sequence by the total counts of the library, and then multiplied by 106. Then, two-step filtering was performed to extract the most likely sRNA—target pairs. During the first step, the predicted sRNA binding sites along with the 50-nt surrounding sequences at both ends were collected in order to reduce the BLAST time. For the BLAST, all the collected degradome data sets were utilized at the same time to do a comprehensive search. It was based on the scenario that a sRNA—target pair was considered to be the candidate once the cleavage signal(s) existed in any data set(s). The predicted targets met the following criteria were retained: (1) there must be perfectly matched degradome signatures with their 5’ ends resided within 8th to 12th nt 20

region from the 5’ ends of the regulatory sRNAs; and (2) for a specific position within the 8th to 12th nt region, which could be regarded as the potential cleavage site, there must be two or more distinct degradome signatures with accordant 5’ ends to support this position. These retained transcripts were subjected to a second BLAST, and the degradome signals along each transcript were obtained to provide a global view of the signal noise when compared to the signal intensity within a specific target binding site. Referring to our previous study (Meng et al, 2011), both global and local t-plots were drawn. Finally, exhaustive manual filtering was performed, and only the transcripts with cleavage signals easy to be recognized were extracted as the potential sRNA—target pairs.

Acknowledgements The authors would like to thank all the publicly available datasets and the scientists behind them.

Funding This work is financially supported by the National Natural Science Foundation of China [31100937], Zhejiang Provincial Natural Science Foundation of China [Y15C060021], and the Starting Grant funded by Hangzhou Normal University to Yijun Meng [2011QDL60].

References Addo-Quaye, C., Eshoo, T.W., Bartel, D.P., Axtell, M.J., 2008. Endogenous siRNA and miRNA targets identified by sequencing of the Arabidopsis degradome. Current Biol. 18, 758-762. Addo-Quaye, C., Snyder, J.A., Park, Y.B., Li, Y.F., Sunkar, R., Axtell, M.J., 2009. Sliced microRNA targets and precise loop-first processing of MIR319 hairpins revealed by analysis of the Physcomitrella patens degradome. RNA 15, 2112-2121. 21

Arenas-Huertero, C., Perez, B., Rabanal, F., Blanco-Melo, D., De la Rosa, C., Estrada-Navarrete, G., Sanchez, F., Covarrubias, A.A., Reyes, J.L., 2009. Conserved and novel miRNAs in the legume Phaseolus vulgaris in response to stress. Plant Mol. Biol. 70, 385-401. Bari, R., Datt Pant, B., Stitt, M., Scheible, W.R., 2006. PHO2, microRNA399, and PHR1 define a phosphate-signaling pathway in plants. Plant Physiol. 141, 988-999. Barrett, T., Troup, D.B., Wilhite, S.E., Ledoux, P., Rudnev, D., Evangelista, C., Kim, I.F., Soboleva, A., Tomashevsky, M., Marshall, K.A., et al. 2009. NCBI GEO: archive for high-throughput functional genomic data. Nucleic Acids Res. 37, D885-890. Baumberger, N., Baulcombe, D.C., 2005. Arabidopsis ARGONAUTE1 is an RNA Slicer that selectively recruits microRNAs and short interfering RNAs. Proc. Natl. Acad. Sci. U. S. A. 102, 11928-11933. Bhardwaj, A.R., Joshi, G., Pandey, R., Kukreja, B., Goel, S., Jagannath, A., Kumar, A., Katiyar-Agarwal, S., Agarwal, M., 2014. A genome-wide perspective of miRNAome in response to high temperature, salinity and drought stresses in Brassica juncea (Czern) L. PLoS One 9, e92456. Borsani, O., Zhu, J., Verslues, P.E., Sunkar, R., Zhu, J.K., 2005. Endogenous siRNAs derived from a pair of natural cis-antisense transcripts regulate salt tolerance in Arabidopsis. Cell 123, 1279-1291. Cao, X., Wu, Z., Jiang, F., Zhou, R., Yang, Z., 2014. Identification of chilling stress-responsive tomato microRNAs and their target genes by high-throughput sequencing and degradome analysis. BMC Genomics 15, 1130. Chen, X., 2009. Small RNAs and their roles in plant development. Annu. Rev. Cell Dev. Biol. 25, 21-44. Chen, H.J., Chen, J.Y., Wang, S.J., 2008. Molecular regulation of starch accumulation in rice seedling leaves in response to salt stress. Acta Physiologiae Plantarum 30, 135-142. Cloonan, N., Wani, S., Xu, Q., Gu, J., Lea, K., Heater, S., Barbacioru, C., Steptoe, 22

A.L., Martin, H.C., Nourbakhsh, E., et al. 2011. MicroRNAs and their isomiRs function cooperatively to target common biological pathways. Genome Biol. 12, R126. Contreras-Cubas, C., Palomar, M., Arteaga-Vazquez, M., Reyes, J.L., Covarrubias, A.A., 2012 Non-coding RNAs in the plant response to abiotic stress. Planta 236, 943-958. Dai, X., Xu, Y., Ma, Q., Xu, W., Wang, T., Xue, Y., Chong, K., 2007. Overexpression of an R1R2R3 MYB gene, OsMYB3R-2, increases tolerance to freezing, drought, and salt stress in transgenic Arabidopsis. Plant Physiol. 143, 1739-1751. Dai, X., Zhao, P.X., 2011. psRNATarget: a plant small RNA target analysis server. Nucleic Acids Res. 39, W155-159. Dash, S., Van Hemert, J., Hong, L., Wise, R.P., Dickerson, J.A., 2012. PLEXdb: gene expression resources for plants and plant pathogens. Nucleic Acids Res. 40, D1194-1201. Ding, Y., Tao, Y., Zhu, C., 2013. Emerging roles of microRNAs in the mediation of drought stress response in plants. J. Exp. Bot. 64, 3077-3086. German, M.A., Luo, S., Schroth, G., Meyers, B.C., Green, P.J., 2009. Construction of Parallel Analysis of RNA Ends (PARE) libraries for the study of cleaved miRNA targets and the RNA degradome. Nat. Protoc. 4, 356-362. German, M.A., Pillay, M., Jeong, D.H., Hetawal, A., Luo, S., Janardhanan, P., Kannan, V., Rymarquis, L.A., Nobuta, K., German, R., et al. 2008. Global identification of microRNA-target RNA pairs by parallel analysis of RNA ends. Nat. Biotechnol. 26, 941-946. Griffiths-Jones, S., Saini, H.K., van Dongen, S., Enright, A.J., 2008. miRBase: tools for microRNA genomics. Nucleic Acids Res. 36, D154-158. Guo, L., Li, H., Liang, T., Lu, J., Yang, Q., Ge, Q., Lu, Z., 2012. Consistent isomiR expression patterns and 3' addition events in miRNA gene clusters and families implicate functional and evolutionary relationships. Mol. Biol. Rep. 39, 6699-6706. Hannah, M.A., Heyer, A.G., Hincha, D.K., 2005. A global survey of gene regulation 23

during cold acclimation in Arabidopsis thaliana. PLoS Genetics 1, e26. Hirayama, T., Shinozaki, K., 2010. Research on plant abiotic stress responses in the post-genome era: past, present and future. Plant J. 61, 1041-1052. Hong, J.K., Choi, H.W., Hwang, I.S., Kim, D.S., Kim, N.H., Choi du, S., Kim, Y.J., Hwang, B.K., 2008. Function of a novel GDSL-type pepper lipase gene, CaGLIP1, in disease susceptibility and abiotic stress tolerance. Planta 227, 539-558. Jain, M., Khurana, J.P., 2009. Transcript profiling reveals diverse roles of auxin-responsive genes during reproductive development and abiotic stress in rice. FEBS J. 276, 3148-3162. Jain, M., Nijhawan, A., Arora, R., Agarwal, P., Ray, S., Sharma, P., Kapoor, S., Tyagi, A.K., Khurana, J.P., 2007. F-box proteins in rice. Genome-wide analysis, classification, temporal and spatial gene expression during panicle and seed development, and regulation by light and abiotic stress. Plant Physiol. 143, 1467-1483. Jeong, D.H., Park, S., Zhai, J., Gurazada, S.G., De Paoli, E., Meyers, B.C., Green, P.J., 2011. Massive analysis of rice small RNAs: mechanistic implications of regulated microRNAs and variants for differential target RNA cleavage. Plant Cell 23, 4185-4207. Jones-Rhoades, M.W., Bartel, D.P., 2004. Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol. Cell 14, 787-799. Jones-Rhoades, M.W., Bartel, D.P., Bartel, B., 2006. MicroRNAs and their regulatory roles in plants. Annu. Rev. Plant Biol. 57, 19-53. Kawashima, C.G., Matthewman, C.A., Huang, S., Lee, B.R., Yoshimoto, N., Koprivova, A., Rubio-Somoza, I., Todesco, M., Rathjen, T., Saito, K., et al. 2011. Interplay of SLIM1 and miR395 in the regulation of sulfate assimilation in Arabidopsis. Plant J. 66, 863-876. Kawashima, C.G., Yoshimoto, N., Maruyama-Nakashita, A., Tsuchiya, Y.N., Saito, K., Takahashi, H., Dalmay, T., 2009. Sulphur starvation induces the expression of microRNA-395 and one of its target genes but in different cell types. Plant J. 57, 24

313-321. Leyva-Gonzalez, M.A., Ibarra-Laclette, E., Cruz-Ramirez, A., Herrera-Estrella, L., 2012. Functional and transcriptome analysis reveals an acclimatization strategy for abiotic stress tolerance mediated by Arabidopsis NF-YA family members. PLoS One 7, e48138. Li, S.C., Liao, Y.L., Ho, M.R., Tsai, K.W., Lai, C.H., Lin, W.C., 2012. miRNA arm selection and isomiR distribution in gastric cancer. BMC Genomics 13 Suppl 1, S13. Li, W.X., Oono, Y., Zhu, J., He, X.J., Wu, J.M., Iida, K., Lu, X.Y., Cui, X., Jin, H., Zhu, J.K., 2008. The Arabidopsis NFYA5 transcription factor is regulated transcriptionally and posttranscriptionally to promote drought resistance. Plant Cell 20, 2238-2251 Li, Y.F., Zheng, Y., Addo-Quaye, C., Zhang, L., Saini, A., Jagadeeswaran, G., Axtell, M.J., Zhang, W., Sunkar, R., 2010. Transcriptome-wide identification of microRNA targets in rice. Plant J. 62, 742-759. Liu, D., Song, Y., Chen, Z., Yu, D., 2009. Ectopic expression of miR396 suppresses GRF target gene expression and alters leaf growth in Arabidopsis. Physiol. Plant. 136, 223-236. Liu, H.H., Tian, X., Li, Y.J., Wu, C.A., Zheng, C.C., 2008. Microarray-based analysis of stress-regulated microRNAs in Arabidopsis thaliana. RNA 14, 836-843. Llave, C., Xie, Z., Kasschau, K.D., Carrington, J.C., 2002. Cleavage of Scarecrow-like mRNA targets directed by a class of Arabidopsis miRNA. Science 297, 2053-2056. Llorens, F., Banez-Coronel, M., Pantano, L., Del Rio, J.A., Ferrer, I., Estivill, X., Marti, E., 2013. A highly expressed miR-101 isomiR is a functional silencing small RNA. BMC Genomics 14, 104. Ma, H.S., Liang, D., Shuai, P., Xia, X.L., Yin, W.L., 2010. The salt- and drought-inducible poplar GRAS protein SCL7 confers salt and drought tolerance in Arabidopsis thaliana. J. Exp. Bot. 61, 4011-4019. Meng, Y., Shao, C., Chen, M., 2011. Toward microRNA-mediated gene regulatory 25

networks in plants. Brief. Bioinform. 12, 645-659. Mi, S., Cai, T., Hu, Y., Chen, Y., Hodges, E., Ni, F., Wu, L., Li, S., Zhou, H., Long, C., et al. 2008. Sorting of small RNAs into Arabidopsis argonaute complexes is directed by the 5' terminal nucleotide. Cell 133, 116-127. Mirouze, M., Paszkowski, J., 2011. Epigenetic contribution to stress adaptation in plants. Curr. Opin. Plant Biol. 14, 267-274. Moumeni, A., Satoh, K., Kondoh, H., Asano, T., Hosaka, A., Venuprasad, R., Serraj, R., Kumar, A., Leung, H., Kikuchi, S., 2011. Comparative analysis of root transcriptome profiles of two pairs of drought-tolerant and susceptible rice near-isogenic lines under different drought stress. BMC Plant Biol. 11, 174. Nakano, M., Nobuta, K., Vemaraju, K., Tej, S.S., Skogen, J.W., Meyers, B.C., 2006. Plant MPSS databases: signature-based transcriptional resources for analyses of mRNA and small RNA. Nucleic Acids Res. 34, D731-735. Naranjo, M.A., Forment, J., Roldan, M., Serrano, R., Vicente, O., 2006. Overexpression of Arabidopsis thaliana LTL1, a salt-induced gene encoding a GDSL-motif lipase, increases salt tolerance in yeast and transgenic plants. Plant Cell Environ. 29, 1890-1900. Nelson, D.E., Repetti, P.P., Adams, T.R., Creelman, R.A., Wu, J., Warner, D.C., Anstrom, D.C., Bensen, R.J., Castiglioni, P.P., Donnarummo, M.G., et al. 2007. Plant nuclear factor Y (NF-Y) B subunits confer drought tolerance and lead to improved corn yields on water-limited acres. Proc. Natl. Acad. Sci. U. S. A. 104, 16450-16455. Ni, Z., Hu, Z., Jiang, Q., Zhang, H., 2012. Overexpression of gma-MIR394a confers tolerance to drought in transgenic Arabidopsis thaliana. Biochem. Biophys. Res. Commun. 427, 330-335. Pandit, A., Rai, V., Sharma, T.R., Sharma, P.C., Singh, N.K., 2011. Differentially expressed genes in sensitive and tolerant rice varieties in response to salt-stress. J. Plant Biochem. Biotechnol. 20, 149-154. Park, M.R., Yun, K.Y., Mohanty, B., Herath, V., Xu, F., Wijaya, E., Bajic, V.B., Yun, S.J., De Los Reyes, B.G., 2010. Supra-optimal expression of the cold-regulated 26

OsMyb4 transcription factor in transgenic rice changes the complexity of transcriptional network with major effects on stress tolerance and panicle development. Plant Cell Environ. 33, 2209-2230. Seki, M., Umezawa, T., Urano, K., Shinozaki, K., 2007. Regulatory metabolic networks in drought stress responses. Curr. Opin. Plant Biol. 10, 296-302. Shannon, P., Markiel, A., Ozier, O., Baliga, N.S., Wang, J.T., Ramage, D., Amin, N., Schwikowski, B., Ideker, T., 2003. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498-2504. Shinozaki, K., Yamaguchi-Shinozaki, K., Seki, M., 2003. Regulatory network of gene expression in the drought and cold stress responses. Curr. Opin. Plant Biol. 6, 410-417. Sunkar, R., Chinnusamy, V., Zhu, J., Zhu, J.K., 2007. Small RNAs as big players in plant abiotic stress responses and nutrient deprivation. Trends Plant Sci. 12, 301-309. Tognetti, V.B., Van Aken, O., Morreel, K., Vandenbroucke, K., van de Cotte, B., De Clercq, I., Chiwocha, S., Fenske, R., Prinsen, E., Boerjan, W., et al. 2010. Perturbation of indole-3-butyric acid homeostasis by the UDP-glucosyltransferase UGT74E2 modulates Arabidopsis architecture and water stress tolerance. Plant Cell 22, 2660-2679. Trindade, I., Capitao, C., Dalmay, T., Fevereiro, M.P., Santos, D.M., 2010. miR398 and miR408 are up-regulated in response to water deficit in Medicago truncatula. Planta 231, 705-716. Vannini, C., Locatelli, F., Bracale, M., Magnani, E., Marsoni, M., Osnato, M., Mattana, M., Baldoni, E., Coraggio, I., 2004. Overexpression of the rice Osmyb4 gene increases chilling and freezing tolerance of Arabidopsis thaliana plants. Plant J. 37, 115-127. Vaucheret, H., 2008. Plant ARGONAUTES. Trends Plant Sci. 13, 350-358. Vidal, E.A., Araus, V., Lu, C., Parry, G., Green, P.J., Coruzzi, G.M., Gutierrez, R.A., 2010. Nitrate-responsive miR393/AFB3 regulatory module controls root system 27

architecture in Arabidopsis thaliana. Proc. Natl. Acad. Sci. U. S. A. 107, 4477-4482. Voinnet, O., 2009. Origin, biogenesis, and activity of plant microRNAs. Cell 136, 669-687. Wang, D., Pan, Y., Zhao, X., Zhu, L., Fu, B., Li, Z., 2011. Genome-wide temporal-spatial gene expression profiling of drought responsiveness in rice. BMC Genomics 12, 149. Xiong, L., Schumaker, K.S., Zhu, J.K., 2002. Cell signaling during cold, drought, and salt stress. Plant Cell 14 Suppl, S165-183. Yan, Y.S., Chen, X.Y., Yang, K., Sun, Z.X., Fu, Y.P., Zhang, Y.M., Fang, R.X., 2011. Overexpression of an F-box protein gene reduces abiotic stress tolerance and promotes root growth in rice. Mol. Plant 4, 190-197. Yang, A., Dai, X., Zhang, W.H., 2012. A R2R3-type MYB gene, OsMYB2, is involved in salt, cold, and dehydration tolerance in rice. J. Exp. Bot. 63, 2541-2556. Yang, M., Yang, Q., Fu, T., Zhou, Y., 2011. Overexpression of the Brassica napus BnLAS gene in Arabidopsis affects plant development and increases drought tolerance. Plant Cell Rep. 30, 373-388. Yu, H., Chen, X., Hong, Y.Y., Wang, Y., Xu, P., Ke, S.D., Liu, H.Y., Zhu, J.K., Oliver, D.J., Xiang, C.B., 2008. Activated expression of an Arabidopsis HD-START protein confers drought tolerance with improved root system and reduced stomatal density. Plant Cell 20, 1134-1151. Yuan, Q., Ouyang, S., Liu, J., Suh, B., Cheung, F., Sultana, R., Lee, D., Quackenbush, J., Buell, C.R., 2003. The TIGR rice genome annotation resource: annotating the rice genome and creating resources for plant biologists. Nucleic Acids Res. 31, 229-233. Zhang, Y., 2005. miRU: an automated plant miRNA target prediction server. Nucleic Acids Res. 33, W701-704. Zhao, B., Liang, R., Ge, L., Li, W., Xiao, H., Lin, H., Ruan, K., Jin, Y., 2007. Identification of drought-induced microRNAs in rice. Biochem. Biophys. Res. 28

Commun. 354, 585-590. Zhou, J., Jiao, F., Wu, Z., Li, Y., Wang, X., He, X., Zhong, W., Wu, P., 2008. OsPHR2 is involved in phosphate-starvation signaling and excessive phosphate accumulation in shoots of plants. Plant Physiol. 146, 1673-1686. Zhou, L., Liu, Y., Liu, Z., Kong, D., Duan, M., Luo, L., 2010. Genome-wide identification and analysis of drought-responsive microRNAs in Oryza sativa. J. Exp. Bot. 61, 4157-68. Zhu, J.K., 2002. Salt and drought stress signal transduction in plants. Annu. Rev. Plant Biol. 53, 247-273.

29

Fig. 1. Workflow for the identification of the rice small RNAs (sRNAs) responsive to the abiotic stresses. The IDs of the sRNA high-throughput sequencing data sets utilized in this study along with the total numbers of the sRNA reads in each data set are shown in the top two tables. The total numbers of the identified stress-responsive sRNAs including microRNAs are shown in the middle table. The total numbers of the Argonaute 1 (AGO1)-enriched, stress-responsive sRNAs including microRNAs are shown in the bottom table. For the middle and the bottom tables, the numbers in the brackets indicate the microRNA sequence counts (please note that different members of a miRNA family may share the same sequence) identified in this study. The rules for the two-step filtering are summarized in the figure.

Fig. 2. Sequence characteristics of the small RNAs (sRNAs) and the microRNAs (miRNAs) responsive to the abiotic stresses. (A) 5’ terminal nucleotide compositions of the miRNAs responsive to the abiotic stresses. (B) 5’ terminal nucleotide compositions of the sRNAs responsive to the abiotic stresses. (C) Sequence length distribution of the miRNAs responsive to the abiotic stresses. (D) Sequence length distribution of the sRNAs responsive to the abiotic stresses. For (A) to (D), “P” represents “panicle” and “S” represents “seedling”. For examples, “P_Cold_Down” represents the sRNAs or the miRNAs down-regulated under the cold treatment, and “S_Drought_Up” represents the sRNAs or the miRNAs up-regulated under the drought treatment.

Fig. 3. Degradome sequencing data-based validation of the targets regulated by the rice small RNAs (sRNAs) responsive to the abiotic stresses. For all the figure panels, the IDs of the target transcripts and the regulatory sRNAs are listed on the top. The y axes measure the intensity (in RPM, reads per million) of the degradome signals, and the x axes indicate the positions of the cleavage signals on the target transcripts. The binding sites of the sRNA regulators on their target transcripts were denoted by gray horizontal lines, and the dominant cleavage signals were marked by gray arrowheads. The figure keys for the degradome signatures from four libraries are shown at the 30

bottom of the figure.

Fig. 4. Certain subnetworks mediated by the stress-responsive small RNAs (sRNAs) in rice. (A) The subnetworks mediated by the sRNAs up-regulated under the drought treatment in rice panicles. (B) The subnetworks mediated by the sRNAs down-regulated under the drought treatment in rice panicles. (C) The subnetworks mediated by the sRNAs down-regulated under the drought treatment in rice seedlings. (D) The subnetworks mediated by the sRNAs up-regulated under the drought treatment in rice seedlings. (E) The subnetworks mediated by the sRNAs down-regulated under the salt treatment in rice seedlings. (F) The subnetworks mediated by the sRNAs down-regulated under the cold treatment in rice seedlings. (G) The subnetworks mediated by the sRNAs up-regulated under the salt treatment in rice seedlings. Please note, the sRNA—target interactions are supported by degradome sequencing data. For all the subnetworks, the red nodes and the blue nodes represent the sRNAs up-regulated and down-regulated under the stress treatments, respectively. The gray nodes represent the targets. The target nodes highlighted in green or purple color represent the target genes with stress-responsive expression patterns supported by microarray data. And, the target nodes with red or blue borderlines indicate that the biological functions of the genes are potentially involved in stress response based on literature mining.

Fig. 5. Statistical result of the numbers of the small RNAs (sRNAs) responsive to multiple stress signals in rice. Twelve categories of the stress-responsive sRNAs are listed in the first column of the table. For example, “Panicle Cold_Down” represents the sRNAs down-regulated under the cold treatment in panicles, and “Seedling Salt_Up” represents the sRNAs up-regulated under the salt treatment in seedlings. The numbers of the sRNAs responsive to multiple stress signals which could be assigned to four or more of the 12 categories (assignments were indicated by colored boxes in each column; the blue color indicates “down-regulated” and the orange color indicates “up-regulated”) were counted. The numbers of the sRNA sequences 31

including miRNAs (in the brackets) are listed in the chart corresponding to each column of the table below.

Fig. 6. Target-centered subnetworks involving multiple stress signals potentially transmitted by specific small RNAs (sRNAs). (A) LOC_Os01g11430.1-centered subnetwork. (B) LOC_Os01g69940.1-centered subnetwork. (C) LOC_Os02g36924.1-centered subnetwork. (D) LOC_Os02g44360.1-centered subnetwork. (E) LOC_Os06g40330.1-centered subnetwork. (F) LOC_Os02g41800.1-centered subnetwork. (G) LOC_Os04g08415.1-centered subnetwork. (H) LOC_Os02g45570.1-centered subnetwork. (I) LOC_Os04g43910.1-centered subnetwork. (J) LOC_Os10g39970.1-centered subnetwork. For (A) to (J), the gray circular nodes represent the targets. The colored (pink for “up-regulated” and blue for “down-regulated”) circular, square and triangular nodes denote the sRNAs responsive to cold, drought and salt treatments respectively (see figure keys on the bottom right). Both the target IDs and the sRNA IDs were marked on the corresponding nodes. The stress-responsive sRNAs identified in seedlings and in panicles are displayed on the upper part and the lower part of each figure panel, separated by a gray horizontal line. For each panel, the nodes with the borderlines in the same color indicate that the sRNAs with different IDs share the same sequence. For all the subnetworks, the sRNA—target interactions were supported by the degradome sequencing data.

32

Supporting Information Data S1 Lists of stress-responsive (including “down-regulated” and “up-regulated”) small RNAs identified in the rice seedlings and the panicles. Data S2 Lists of stress-responsive small RNAs enriched in the Argonaute 1-associated silencing complexes of rice. Data S3 List of the target genes supported by the microarray data. Data S4 Statistical result of the small RNAs potentially involved in stress signal interactions. Figure S1 Sequence characteristics of the stress-responsive small RNAs enriched in the Argonaute 1-associated silencing complexes of rice. Figure S2 Results of the degradome sequencing data-based validation of the small RNA—target pairs. Figure S3 Networks mediated by the stress-responsive small RNAs in rice seedlings and panicles. Table S1 Detailed information of the high-throughput sequencing data sets used in this study. Table S2 Statistical result of the stress responsiveness of the microRNAs and their isomiRs in rice. Table S3 List of the target transcripts simultaneously regulated by several stress-responsive small RNAs assigned to two or more of the 12 categories in rice. Table S4 List of the rice microRNAs responsive to two or more stress signals within the networks established based on degradome sequencing data (Figure S3).

33

Figure 1

34

Figure 2

35

Figure 3

36

Figure 4

37

Figure 5

38

Figure 6

39