Integration Analysis of Small RNA and Degradome Sequencing Reveals MicroRNAs Responsive to Dickeya zeae in Resistant Rice

Rice foot rot disease caused by the pathogen Dickeya zeae (formerly known as Erwinia chrysanthemi pv. zeae), is a newly emerging damaging bacterial disease in China and the southeast of Asia, resulting in the loss of yield and grain quality. However, the genetic resistance mechanisms mediated by miRNAs to D. zeae are unclear in rice. In the present study, 652 miRNAs including osa-miR396f predicted to be involved in multiple defense responses to D. zeae were identified with RNA sequencing. A total of 79 differentially expressed miRNAs were detected under the criterion of normalized reads ≥10, including 51 known and 28 novel miRNAs. Degradome sequencing identified 799 targets predicted to be cleaved by 168 identified miRNAs. Among them, 29 differentially expressed miRNA and target pairs including miRNA396f-OsGRFs were identified by co-expression analysis. Overexpression of the osa-miR396f precursor in a susceptible rice variety showed enhanced resistance to D. zeae, coupled with significant accumulation of transcripts of osa-miR396f and reduction of its target the Growth-Regulating Factors (OsGRFs). Taken together, these findings suggest that miRNA and targets including miR396f-OsGRFs have a role in resisting the infections by bacteria D. zeae.

Most recently, miRNAs were demonstrated to be a critical component of sophisticated plant defense system [1][2][3]. Plants are known to use a multifaceted defense system to prevent infections by pathogenic microbes through a co-evolutionary mechanism with pathogens. Plant immunity is involved in two tiers of defenses where the first tier is a pathogen-associated molecular pattern (PAMP) trigged immunity (PTI) [4,5]. The second tier of defense is often involved in pathogen effector-triggered immunity (ETI) mediated mostly by R proteins with leucine rich repeat and nucleotide binding site

D. zeae-Responsive miRNAs in Resistant Rice
To identify miRNAs in resistant rice that respond to D. zeae, we analyzed the differential expressions of miRNAs in the four libraries with the normalized reads from deep sequencing. A total of 79 differentially expressed miRNAs was found in all four libraries including 51 known miRNAs and 28 novel miRNAs under the criterion of normalized reads ≥10 ( Figure 2A; Table S3). Furthermore, 34 miRNAs were up-regulated and 32 miRNAs were down-regulated at six hours post inoculation (hpi). At 12 hpi, 36 miRNAs were up-regulated and 32 miRNAs were down-regulated. At 48 hpi, 26 miRNAs were up-regulated and 39 miRNAs were down-regulated ( Figure 2B; Table S3). These results suggest that certain miRNAs were expressed at some time points after D. zeae infection in an incompatible interaction suggesting that they may play an important role during pathogen infection.
highest frequencies were identified ( Figure 1; Table S2). The miRNAs mapped to the rice genome but did not match to known pre-miRNAs as novel miRNAs. 250 novel miRNAs were predicted in four libraries ( Figure 1; Table S2). The majority of novel miRNAs was found to be slowly accumulated after inoculation with D. zeae. However, osa-MIR5083-p5 was the most abundant one with reads of 135.04 (K0), 484.00 (K1), 713.11 (K2) and 233.71 (K3) in four libraries, respectively, suggesting its unique role for responding to D. zeae infection in a resistant rice variety. CountsMIRb, the counts of miRNAs from miRBase. Expression level, low indicates <10, middle indicates >10 but less than average, high indicates over average.

D. zeae-Responsive miRNAs in Resistant Rice
To identify miRNAs in resistant rice that respond to D. zeae, we analyzed the differential expressions of miRNAs in the four libraries with the normalized reads from deep sequencing. A total of 79 differentially expressed miRNAs was found in all four libraries including 51 known miRNAs and 28 novel miRNAs under the criterion of normalized reads ≥10 ( Figure 2A; Table S3). Furthermore, 34 miRNAs were up-regulated and 32 miRNAs were down-regulated at six hours post inoculation (hpi). At 12 hpi, 36 miRNAs were up-regulated and 32 miRNAs were down-regulated. At 48 hpi, 26 miRNAs were up-regulated and 39 miRNAs were down-regulated ( Figure 2B; Table  S3). These results suggest that certain miRNAs were expressed at some time points after D. zeae infection in an incompatible interaction suggesting that they may play an important role during pathogen infection.   K3, respectively ( Figure 3). When compared with K0, 165 unique miRNAs were identified at 6 hpi, ( Figure 3A), 163 unique miRNAs were identified at 12 hpi, 94 unique miRNAs were detected at 48 hpi, respectively ( Figure 3B,C). When compared among four libraries, 30 (K0), 44 (K1), 41 (K2) and 18 (K3) miRNAs were detected ( Figure 3D). The presences of different miRNAs at different time points after inoculation suggest miRNAs have a role to D. zeae in an incompatible interaction.

Target Prediction of miRNAs by Degradome Sequencing and Analysis
To identify the relevant targets of miRNAs, high-throughput degradome sequencing was used and the data were analyzed with the CleavLand 3 pipeline. Degradome sequencing analysis showed that 10931903 and 8679632 raw reads were detected from the control (T0 h) and infected (T6-48 h) rice samples, respectively (Table S4). After removing the adaptor reads and mapping the rice transcriptome, 35890 and 32146 covered transcripts were identified in the control (T0 h) and infected (T6-48) rice samples, respectively (Table S4). Further, a total of 799 targets were predicted to be cleaved by the 168 identified miRNAs in an incompatible interaction (Table S5). Many transcript targets mediated by miRNAs were closely related to the signal transduction involved in plant immunity, such as PR5 (LOC_Os03g14030.1), leucine-rich repeat (LRR) kinase (LOC_Os12g10740.1), calcium-dependent kinase (LOC_Os07g22640.), WRKY transcription factor (LOC_Os04g51560.1), AUX/IAA transcriptional regulator (LOC_Os12g40900.1), ethylene-responsive binding protein (LOC_Os07g42510.1), growth-regulating factor 5 (LOC_Os02g53690.1) and etc ( Table 2). These findings demonstrate that miRNAs-target pairs may be involved in resistance responses through multiple defense pathways in rice.

Functional Classification of Predicted Targets
A gene ontology (GO) analysis was used to classify the biochemical functions of predicted targets ( Figure 4). A total of 858 transcript targets were predicted as protein kinase, 852 as transferase, 171 as signal transducers, and 81 as protein receptors, suggesting their potential involvements in signal recognition and transduction in plant immunity (Figure 4). Similarly, a group of targets was predicted to be involved in multiple molecular function and all of them were predicted to be located in diverse cellular components of cells, suggesting that predicted target genes may also play a role in biological processes including hormone stimulus, signal transducer, sequence-specific DNA binding, secondary metabolic process, cell growth, and multicellular organismal development ( Figure 4). of targets was predicted to be involved in multiple molecular function and all of them were predicted to be located in diverse cellular components of cells, suggesting that predicted target genes may also play a role in biological processes including hormone stimulus, signal transducer, sequence-specific DNA binding, secondary metabolic process, cell growth, and multicellular organismal development (Figure 4).

Co-expression Analysis and Validation of miRNA-Targets by qRT-PCR
We analyzed the miRNAs obtained by miRNA sequencing, and verified the data with qRT-PCR and degradome sequencing. Firstly, 29 (36.71%) of 79 differentially expressed miRNAs detected by miRNA-seq were detected from the degradome sequencing (Table S6). A total of 82 genes resulting in 132 transcripts were predicted to be targets of these 29 miRNAs (Table S6). Moreover, 9 differentially expressed miRNAs based on the number of normalized reads were selected from 29 miRNAs (Table 3). Interestingly, we found that among 9 miRNAs, 3 miRNAs were induced and 6 were reduced (Table S2). Moreover, the secondary structures of 9 miRNAs precursor were predicted by RNAfold respectively ( Figure 5A). Despite the difference in the fold changes, the expression patterns of 9 miRNAs analyzed by qRT-PCR were consistent with that of miRNA-seq ( Figure 5B).

Co-expression Analysis and Validation of miRNA-Targets by qRT-PCR
We analyzed the miRNAs obtained by miRNA sequencing, and verified the data with qRT-PCR and degradome sequencing. Firstly, 29 (36.71%) of 79 differentially expressed miRNAs detected by miRNA-seq were detected from the degradome sequencing (Table S6). A total of 82 genes resulting in 132 transcripts were predicted to be targets of these 29 miRNAs (Table S6). Moreover, 9 differentially expressed miRNAs based on the number of normalized reads were selected from 29 miRNAs (Table 3). Interestingly, we found that among 9 miRNAs, 3 miRNAs were induced and 6 were reduced (Table S2). Moreover, the secondary structures of 9 miRNAs precursor were predicted by RNAfold respectively (Figure 5A). Despite the difference in the fold changes, the expression patterns of 9 miRNAs analyzed by qRT-PCR were consistent with that of miRNA-seq ( Figure 5B). Moreover, the degradome sequencing data revealed an opposite expression pattern of the cognate targets of these 9 miRNAs (Tables S2 and S6) suggesting that these miRNAs may regulate target gene expression at different time points post inoculation in an incompatible interaction. Reduced expressions of three target genes, LOC_Os06g02560.1, LOC_Os03g51970.1 and LOC_Os03g47140.1) predicted to be regulated by osa-miR396f were verified by qRT-PCR, suggesting that the osa-miR396f-targets pair may play a role for resistance to D. zeae in rice ( Figure 6). Osa-miR5072 _L-4 LOC_Os05g07050.1 Pre-mRNA processing splicing factor 45.74 0 down 5043 expression at different time points post inoculation in an incompatible interaction. Red ssions of three target genes, LOC_Os06g02560.1, LOC_Os03g51970.1 and LOC_Os03g471 cted to be regulated by osa-miR396f were verified by qRT-PCR, suggesting that iR396f-targets pair may play a role for resistance to D. zeae in rice ( Figure 6).

Overexpression of Osa-miR396f Precursor Enhanced Rice Resistance to D. zeae
To validate if osa-miR396f was involved in defense response to D. zeae, the precursor sequences of osa-miR396f were isolated from the resistant cultivar Nanjing 40 and were expressed under the control of the constitutive CMV35S promoter and transformed into susceptible japonica rice cultivar Nipponbare. Seven positive osa-miR396f precursor-overexpressing rice plants were inoculated with D. zeae. The transgenic rice plants with the osa-miR396f precursor showed increased resistance to D. zeae in comparison with that of Nipponbare ( Figure 7A). Consistently, the transcripts of osa-miR396f precursor in all seven independent transgenic lines were higher than that in wild-type plants ( Figure  7B). As predicted, transcripts of all three targets OsGRFs of osa-miR396f were significantly reduced in these transgenics revealed by qRT-PCR ( Figure 7C). These results suggest that miRNA396f is involved in resisting D. zeae infection in rice.

Overexpression of Osa-miR396f Precursor Enhanced Rice Resistance to D. zeae
To validate if osa-miR396f was involved in defense response to D. zeae, the precursor sequences of osa-miR396f were isolated from the resistant cultivar Nanjing 40 and were expressed under the control of the constitutive CMV35S promoter and transformed into susceptible japonica rice cultivar Nipponbare. Seven positive osa-miR396f precursor-overexpressing rice plants were inoculated with D. zeae. The transgenic rice plants with the osa-miR396f precursor showed increased resistance to D. zeae in comparison with that of Nipponbare ( Figure 7A). Consistently, the transcripts of osa-miR396f precursor in all seven independent transgenic lines were higher than that in wild-type plants ( Figure 7B). As predicted, transcripts of all three targets OsGRFs of osa-miR396f were significantly reduced in these transgenics revealed by qRT-PCR ( Figure 7C). These results suggest that miRNA396f is involved in resisting D. zeae infection in rice.

Discussion
Recently, rice bacterial foot rot disease caused by D. zeae has the potential to become one of the most important bacterial diseases, resulting in the loss of yield and grain quality in China, in the southeast of Asia and worldwide. In the present study we identified and characterized a number of miRNAs-target pairs involved in defense responses in rice by miRNA, degradome sequencing and qRT-PCR.
The fast evolved high-throughput sequencing technology allows rapid identification of miRNAs and targets in plants that are involved in plant immunity [3,31]. Accumulated studies suggest that miRNAs played an important role for regulating the biological process and stress responses in plants [32,33]. It is known that plant miRNAs expressed differently during the infection process of pathogens [7]. In the present study, we found 652 miRNAs including 79 significantly differentially expressed miRNAs at different time points post inoculation in an incompatible interaction (Figure 2A; Table S3). For instance, two miRNAs miR827n-5p and miR164 we identified were also identified in a resistant rice variety to M. grisea, as negative and positive regulator respectively [21]. Moreover, the expression levels of seven miRNAs, miR156a, 159b, 166e-3p, 394, 396c-3p, 812 and 827 were decreased in rice in resistance response to Xoo [20]. We showed a similar responsive mode of these miRNAs to fungal and bacterial pathogens infection demonstrating that miRNAs were indeed involved in resistance response to D. zeae in rice. We also validated the expression profiles of nine miRNAs (miR2118, 393, 396, 166, 171, 156, 535, 159 and 5072) related to defense responses by qRT-PCR ( Figure 5). The relevant targets of these nine miRNAs were also identified and verified by degradome sequencing and qRT-PCR (Table S6; Figure 6). These 9 miRNAs were previously shown to respond to the pathogens infection in an

Discussion
Recently, rice bacterial foot rot disease caused by D. zeae has the potential to become one of the most important bacterial diseases, resulting in the loss of yield and grain quality in China, in the southeast of Asia and worldwide. In the present study we identified and characterized a number of miRNAs-target pairs involved in defense responses in rice by miRNA, degradome sequencing and qRT-PCR.
The fast evolved high-throughput sequencing technology allows rapid identification of miRNAs and targets in plants that are involved in plant immunity [3,31]. Accumulated studies suggest that miRNAs played an important role for regulating the biological process and stress responses in plants [32,33]. It is known that plant miRNAs expressed differently during the infection process of pathogens [7]. In the present study, we found 652 miRNAs including 79 significantly differentially expressed miRNAs at different time points post inoculation in an incompatible interaction (Figure 2A; Table S3). For instance, two miRNAs miR827n-5p and miR164 we identified were also identified in a resistant rice variety to M. grisea, as negative and positive regulator respectively [21]. Moreover, the expression levels of seven miRNAs, miR156a, 159b, 166e-3p, 394, 396c-3p, 812 and 827 were decreased in rice in resistance response to Xoo [20]. We showed a similar responsive mode of these miRNAs to fungal and bacterial pathogens infection demonstrating that miRNAs were indeed involved in resistance response to D. zeae in rice. We also validated the expression profiles of nine miRNAs (miR2118, 393, 396, 166, 171, 156, 535, 159 and 5072) related to defense responses by qRT-PCR ( Figure 5). The relevant targets of these nine miRNAs were also identified and verified by degradome sequencing and qRT-PCR (Table S6; Figure 6). These 9 miRNAs were previously shown to respond to the pathogens infection in an incompatible interaction, probably resulting in the activation of multiple defense responses [34][35][36][37][38]. Our data suggest that these nine differentially expressed miRNAs may be involved in regulating the resistance to D. zeae in rice.
In the present study, we showed that overexpression of osa-miR396f precursor in rice enhanced resistance to D. zeae, correlated with significant reduction of transcripts of three predicted targets OsGRFs (Figure 7). Therefore, miR396 may modulate resistance to D. zeae by cleaving the transcripts of its targets OsGRFs. Rice miR396 family has eight members, osa-miR396a, 396b, 396c, 396d, 396e, 396f, 396g and 396h in the miRNA database (miRBase, http://www.mirbase.org/, 27 March 2018). Rice miR396-GRFs pair was previously shown to regulate plant growth and development by positively activating plant hormones indole-3-acetic acid (IAA), gibberellin (GA) and brassinolide biosynthesis pathways in rice [39][40][41][42]. The miR396 family members, osa-miR396e-3p and osa-miR396d/e-5p were induced by blast pathogen in a susceptible rice variety, and the expression level of osa-miR396c-5p was increased in a resistant rice variety [21]. Some members of miR396 were previously predicted to be a negative factor for the infection process of southern rice black-streaked dwarf virus (SRBSDV) in rice [38]. In contrast, overexpression of miR396a or miR396b significantly compromised the susceptibility to cyst Nematode in Arabidopsis [43]. The Growth-Regulating Factors (GRFs) were also the targets of miR396 family members and positively regulated the number and size of the rice grain [44]. The blocking of OsGRFs by miR396 generated a larger grain size, auxiliary branches and spikelets then enhanced grain yield through auxin in rice [41,45,46]. We speculate that miR396-OsGRFs pair may play a role for regulating the balance between yield and disease resistance via the auxin signaling pathway. However, the biological functions of miR396 on defense responses and rice grain development still need to be explored.
Plant miRNAs were predicted to be involved in the biological process and stress responses by inhibiting the transcripts of relevant targets [1,13,33]. In the present study, 132 miRNA responsive to the infection of D. zeae in an incompatible interaction were identified, and the predicted targets of 29 differentially expressed miRNAs were identified by degradome sequencing (Table S6). These identified targets might modulate rice resistance response to D. zeae. Most of the targets we identified were predicted to be involved in multiple defense response signal pathways such as OsGRFs, Leucine-rich repeat protein kinase, MYB domain protein, and auxin efflux carrier family protein (Tables 2 and 3;  Table S6). The genes, AtGRF1 and AtGRF3 were shown to be involved in multiple resistance-related signal pathways related with cell-wall modification, cytokinin biosynthesis and the accumulation of secondary metabolites in Arabidopsis [43,47]. Overexpression of miR396a or miR396b and miR396 targets grf1/grf2/grf3 triple mutants previously significantly compromised the susceptibility to cyst nematode in Arabidopsis [43]. The miR396-GRF6 interaction network was predicted to be involved in the development of inflorescence architecture by regulating the auxin biosynthesis and signaling pathways [41]. Moreover, the auxin-responsive GH3 family gene LOC_Os07g40290.1 as the target of osa-miR172d-5p_R-2 showed more cleavages in a resistant rice variety after inoculation with D. zeae, compared with that without inoculation (Table 2). Additionally, auxin as a negative regulator was shown to contribute to resistance to bacteria pathogen Xoo, which resulted in the change of the cell wall structure [45,46]. The Mybs1 gene encoding one MYB transcription factor was involved in a broad-spectrum resistance against M. oryzae in rice by increasing the accumulation of H 2 O 2 through inhibition of the expression level of bsr-d1 [48]. Taken together, these findings suggest that miRNA and targets are involved in multiple defense responses, including auxin and active oxygen signal pathways.
In summary, knowledge of miRNAs and targets opens a new avenue to understand the complexity of plant immunity. In the present study, a group of differentially expressed miRNAs and targets in an incompatible interaction at early interphases of host-pathogen interaction was identified by miRNA and degradome sequencing and qRT-PCR. These miRNAs and targets were predicted to be involved in multiple defense signal and diverse cellular pathways in rice. Through a transgenic approach, we demonstrated that one of these miRNAs, osa-miR396f enhanced the defense responses against D. zeae by inhibiting the transcripts of the three growth related transcription factor genes, OsGRFs. These results not only pave the road for controlling the newly emerging damaging foot rot disease but also promise a better understanding of plant innate immunity.

Plant Materials and Measurement of Resistance Reactions to D. zeae
The foot rot disease resistant Japonica rice variety Nanjing 40 was used to analyze the candidate miRNAs involved in disease resistance and the relevant targets by high-through sequencing technologies. The individual positive rice lines overexpressing the osa-miR396f precursor were chosen for analyzing disease reactions to D. zeae. All the plants used in this study were grown in a greenhouse in Jiangsu Academy of Agricultural Sciences, Nanjing, Jiangsu Province, China. Disease reactions to D. zeae were measured by the basal stem and root inoculation methods, respectively [30,49]. Disease reaction was classified as five groups by measuring the percentage of diseased area at 7-10 days after inoculation. No visible disease symptom in the stem and foot indicates immunity. The disease index less than or equal to five indicates highly resistance (HR). The disease index 5.1 to 12.4 indicates moderate resistance (MR). The disease index, 12.5 to 19.9 indicates moderate susceptibility (MS). The disease index greater than or equal to 20 indicates highly susceptibility (HS) [30,49].

Construction of Small RNA Library, Sequencing and Data Analysis
Total RNA was extracted from four mixed rice roots at 0 (K0), 6 (K1), 12 (K2) and 48 (K3) hpi with D. zeae using Trizol reagent (Invitrogen, Shanghai, China) following the manufacturer's procedure respectively. The rice roots sampled with two biological repeats were combined as K0, K1, K2 and K3 for total RNA extraction, small RNA library construction and sequencing. Approximately 1 mg of total RNAs was used to prepare the miRNA library according to the protocol of TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA, USA). The libraries were sequenced with the single-end sequencing (36 bp) on an Illumina Hiseq2500 at the LC-Bio Co., Ltd (Hangzhou, China) following the vendor's recommended protocol.
The raw reads were analyzed with the Illumina pipeline filter (Solexa 0.3), and the dataset was further processed with an in-house program, ACGT101-miR (LC Sciences, Houston, TX, USA) to remove adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA) and repeats. Subsequently, unique sequences with length in 18-25 nucleotides were mapped to specific species precursors in the miRBase database (http://www.mirbase.org/, 27 March 2018) using the BLAST algorithm to identify known miRNAs and novel 3p and 5p-derived miRNAs. The unique sequences mapping to specific species mature miRNAs in hairpin arms were identified as known miRNAs; the unique sequences mapping to the other arm of known specific species precursor hairpin opposite to the annotated mature miRNA-containing arm were considered to be novel 5p or 3p-derived miRNA candidates [50,51].

Identification and Function Analysis of Targets of miRNAs
The total RNA extracted from control T0h (0 hpi) and infected T6-48h (mixtures of 6, 12, 48 hpi) rice sample respectively were used for degradome sequencing on Illumina Hiseq 2500 at the LC-Bio Co., Ltd. (Hangzhou, China). The degradome reads were mapped to the rice transcriptome data. Publicly available software packages, CleaveL and 3.0 pipeline was used to predict and identify the potentially cleaved targets [50,52]. Then, the computational target prediction algorithms (Target Finder) were used to identify miRNA binding sites [53]. All of the predicted target genes were analyzed with NCBI BLASTX algorithm. Finally, the GO terms of these differentially expressed miRNA targets were also annotated by AgriGO program. The candidate miRNAs-related to defense responses and their targets with differential expression based on normalized deep-sequencing counts were analyzed by the Fisher exact test, Chi-squared test, Student's t test, and ANOVA based on the experiments design [54]. The level of significance threshold was set at 0.01 and 0.05 for each test.

Validation of Differentially Expressed miRNAs and Targets by qRT-PCR
To validate miRNAs and their target expressions determined by high-through sequencing technology, the root fragments at different time points after inoculation were used to extract total RNA using the Trizol reagent (TaKaRa, Dalian, China). The qRT-PCR analysis was performed with the Applied Biosystems 7500 Real Time PCR System and SYBR Premix Ex Taq TM (TaKaRa, Dalian, China) according to the manufacturer's instructions. The rice gene EF1-a and U6 was used as the internal reference gene to standardize RNA quantity for evaluating relative expression levels. For qRT-PCR assays, three independent biological samples were carried out with three technical replicates with a gene-specific primer (Table S7).

Vector Construction and Rice Transformation
The 176nt precursor sequence of osa-miR396f was isolated from Nanjing 40 by PCR amplification, and inserted into the pCAMBIA1301 binary vector driving by the constitutive CMV35s promoter for overexpression. The T-DNA recombinant plasmids were transformed into calli derived from the mature Nipponbare embryos by the Agrobacterium mediated transformation method [46].