High-Throughput Sequencing Reveals H2O2 Stress-Associated MicroRNAs and a Potential Regulatory Network in Brachypodium distachyon Seedlings

Oxidative stress in plants can be triggered by many environmental stress factors, such as drought and salinity. Brachypodium distachyon is a model organism for the study of biofuel plants and crops, such as wheat. Although recent studies have found many oxidative stress response-related proteins, the mechanism of microRNA (miRNA)-mediated oxidative stress response is still unclear. Using next generation high-throughput sequencing technology, the small RNAs were sequenced from the model plant B. distachyon 21 (Bd21) under H2O2 stress and normal growth conditions. In total, 144 known B. distachyon miRNAs and 221 potential new miRNAs were identified. Further analysis of potential new miRNAs suggested that 36 could be clustered into known miRNA families, while the remaining 185 were identified as B. distachyon-specific new miRNAs. Differential analysis of miRNAs from the normal and H2O2 stress libraries identified 31 known and 30 new H2O2 stress responsive miRNAs. The expression patterns of seven representative miRNAs were verified by reverse transcription quantitative polymerase chain reaction (RT-qPCR) analysis, which produced results consistent with those of the deep sequencing method. Moreover, we also performed RT-qPCR analysis to verify the expression levels of 13 target genes and the cleavage site of 5 target genes by known or novel miRNAs were validated experimentally by 5′ RACE. Additionally, a miRNA-mediated gene regulatory network for H2O2 stress response was constructed. Our study identifies a set of H2O2-responsive miRNAs and their target genes and reveals the mechanism of oxidative stress response and defense at the post-transcriptional regulatory level.

Oxidative stress in plants can be triggered by many environmental stress factors, such as drought and salinity. Brachypodium distachyon is a model organism for the study of biofuel plants and crops, such as wheat. Although recent studies have found many oxidative stress response-related proteins, the mechanism of microRNA (miRNA)-mediated oxidative stress response is still unclear. Using next generation high-throughput sequencing technology, the small RNAs were sequenced from the model plant B. distachyon 21 (Bd21) under H 2 O 2 stress and normal growth conditions. In total, 144 known B. distachyon miRNAs and 221 potential new miRNAs were identified. Further analysis of potential new miRNAs suggested that 36 could be clustered into known miRNA families, while the remaining 185 were identified as B. distachyon-specific new miRNAs. Differential analysis of miRNAs from the normal and H 2 O 2 stress libraries identified 31 known and 30 new H 2 O 2 stress responsive miRNAs. The expression patterns of seven representative miRNAs were verified by reverse transcription quantitative polymerase chain reaction (RT-qPCR) analysis, which produced results consistent with those of the deep sequencing method. Moreover, we also performed RT-qPCR analysis to verify the expression levels of 13 target genes and the cleavage site of 5 target genes by known or novel miRNAs were validated experimentally by 5 ′ RACE. Additionally, a miRNA-mediated gene regulatory network for H 2 O 2 stress response was constructed. Our study identifies a set of H 2 O 2 -responsive miRNAs and their target genes and reveals the mechanism of oxidative stress response and defense at the post-transcriptional regulatory level.

BACKGROUND
Many environmental stress factors, including high light, UV irradiation, heat, salinity, drought, and cold, can cause plant cells to produce reactive oxygen species (ROS), leading to acceleration of lipid peroxidation and leaf senescence (Mittler et al., 2004;Upadhyaya et al., 2007). Hydrogen peroxide (H 2 O 2 ) is a kind of ROS mainly produced by the following parts of the cell: the mitochondrion, peroxisome, chloroplast, apoplast, and plasma membrane (Apel and Hirt, 2004). Superoxide radicals (O − 2 ), another kind of ROS, can also be rapidly dismutated to H 2 O 2 spontaneously or catalyzed by superoxide dismutase (SOD). In contrast to other ROS, such as O − 2 and hydroxyl radicals (OH − ), H 2 O 2 can easily pass through membranes (Foyer et al., 1997;Uchida et al., 2002;de Azevedo Neto et al., 2005;Wahid et al., 2007) and is relatively stable, so it is suitable for its roles as an important component of cell signaling cascades (Mittler, 2002;Vranová et al., 2002) and an indispensable second messenger in biotic and abiotic stress responses (Pastori and Foyer, 2002). H 2 O 2 stress can affect fluctuation of the Ca 2+ concentration in plants, thereby inducing the production of appropriate amounts of antioxidants (Rentel and Knight, 2004). Global transcript profiling under H 2 O 2 stress in tobacco revealed that redox homeostasis associated proteins are upregulated while proteins related to normal growth and development are downregulated (Vandenabeele et al., 2003). H 2 O 2 can also trigger acclimation and cross-tolerance phenomena Pastori and Foyer, 2002). Overall, at low concentrations H 2 O 2 may play a role as a signaling molecule, whereas at high concentrations H 2 O 2 will cause programmed cell death (Quan et al., 2008).
MicroRNAs (miRNAs) are small, non-coding RNAs that have been demonstrated to be involved in many responses to and defenses against various biotic and abiotic stresses in plants (Bartel, 2004). With the development of high-throughput sequencing technology and bioinformatics, many plant miRNAs have been identified. In animals, at least 60% of protein-coding genes can be regulated by miRNAs (Friedman et al., 2009), but known target genes of miRNAs in plants are far fewer (∼1% of the protein-coding genes) (Addo-Quaye et al., 2008;Li et al., 2010). Even so, the regulatory role of miRNAs in plants cannot be underestimated, because most known target genes are transcription factors (TFs) (Jones- Rhoades and Bartel, 2004;Jones-Rhoades et al., 2006). Under environmental stresses, plants up-or down-regulate certain miRNAs or synthesize new miRNAs to respond to or defend against stresses (Khraiwesh et al., 2012). Abiotic stresses, such as high light, UV, heat, heavy metals, drought, and salinity, can elevate ROS levels (Mittler et al., 2004). To date, only a few studies have focused on oxidative stress-triggered miRNA expression changes. miR398, whose target genes are Cu-Zn superoxide dismutases, is a well-studied miRNA related to the response to oxidative stress triggered by high light (Sunkar et al., 2006). Iyer et al. (2012) identified 22 ozone-induced oxidative stress miRNA families using a plant miRNA array in Arabidopsis; most of them were also reported as UV-B responsive miRNAs (Zhou et al., 2007). Jia et al. (2009) identified 24 UV-B responsive miRNAs (13 upregulated and 11 downregulated) in Populus tremula through a miRNA filter array. Some of these upregulated miRNAs (miR156, miR160, miR165/166, miR167, miR398, and miR168) were also reported in UV-B-stressed Arabidopsis (Zhou et al., 2007). Similarly, six miRNAs were identified as UV-B-responsive miRNAs in wheat (Wang et al., 2013), in which miR159, miR167a, and miR171 are upregulated and miR156, miR164, miR395 are downregulated. As a ROS, H 2 O 2 can also act as a secondary messenger during stress response and defense signal transduction. However, only a few H 2 O 2 -responsive miRNAs (miR169, miR397, miR1425, miR408-5p, miR827, miR528, and miR319a.2) have been identified in plants (Li et al., 2011). Thus, the miRNAs related to oxidative stress caused by abiotic stressors are far from being completely elucidated in plants.
Brachypodium distachyon, as a model plant for crops such as wheat and barley, has been sequenced (Vogel et al., 2010). Several studies under abiotic stress have been performed at different levels, including the transcriptome (Priest et al., 2014), proteome, and phosphoproteome (Lv et al., 2014b). In addition, recent miRNA studies using high-throughput sequencing have identified many stress responsive miRNAs under several kinds of abiotic stress, such as cold stress (Zhang et al., 2009), dehydration stress (Budak and Akpinar, 2011), and drought stress (Bertolini et al., 2013). In particular, Jeong et al. (2013) sequenced 17 small RNA libraries that represented different tissues and stressors and identified many previously unreported and B. distachyonspecific miRNAs. However, the identified miRNAs are far from sufficient for B. distachyon, especially for oxidative stress regulated miRNAs.
Thus, in this study, we identified miRNAs and their potential target genes related to H 2 O 2 stress using high-throughput sequencing, reverse transcription quantitative polymerase chain reaction (RT-qPCR) and 5 ′ RACE, combined with bioinformatics methods. The differentially expressed miRNAs observed between B. distachyon seedlings grown under control and H 2 O 2 -treated conditions, as well as the miRNA-directed regulatory network, provide new insights that will inform the genetic improvement of stress tolerance in plants.

Plant Materials
Seedlings of B. distachyon 21 (Bd21) were grown in a growth chamber at 25/20 • C (16 h day/8 h night) and 70% relative humidity, as reported previously (Lv et al., 2014a). For H 2 O 2 treatment, seedlings at the three leaves stage were treated with 20 mM H 2 O 2 for 6 h in plastic containers and collected at 2, 4, or 6 h based on our previous study (Bian et al., 2015). Untreated seedlings were used as a control. All of the H 2 O 2 treated and untreated samples had three biological replicates and at least 100 seedlings were used in each replicate. All samples were snapfrozen in liquid nitrogen and then stored at −80 • C until RNA extraction.

Total RNA Isolation
Total RNA was extracted from the frozen seedlings with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. Prior to nucleic precipitation, two extra chloroform washes were performed. A 1% agarose gel stained with ethidium bromide was run to determine the preliminary integrity of the RNA. All RNA samples were quantified and examined using an ND 1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) for contamination with either protein (A260/A280 ratios) or reagent (A260/A230 ratios). The RNA integrity number (RIN) was >8, as determined with a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

Construction of Small RNA (sRNA) Libraries and Deep Sequencing
For sRNA library construction and deep sequencing, RNA samples were prepared as follows: equal quantities (10 mg) of total RNA isolated from Bd21 seedlings treated with 20 mM H 2 O 2 for 2, 4, and 6 h were mixed together to construct the TS library, and 30 mg of total RNA prepared from the control sample (without H 2 O 2 treatment) were used to construct the CS library. Then, total RNA was separated by 15% TBE-urea denaturing polyacrylamide gel electrophoresis (PAGE), and RNA molecules in the range of 18-30 nt were enriched and ligated with proprietary adapters to the 5 ′ and 3 ′ termini by T4 RNA ligase. The samples were used as templates for cDNA synthesis by Super-Script II Reverse Transcriptase (Invitrogen) and the resulting cDNA was amplified to produce sequencing libraries. The final quality of the cDNA library was ensured by examining its size, purity, and concentration with a 2100 Bioanalyzer. The sequencing was performed by the Beijing Genomics Institute (BGI, Shenzhen, China). The two libraries were run on the Illumina HiSeq TM 2000 platform side by side.

Differential Expression Analysis
The relative miRNA expression levels of the two libraries were compared and the differentially expressed miRNAs were screened based on a previously established method (Audic and Claverie, 1997). The frequency of miRNAs in the two libraries was normalized to one million by the total number of miRNAs in each library (transcripts per million (TPM) normalized expression = initial miRNA count * 1,000,000/total count of clean reads). Following normalization, if the miRNA gene expression in both libraries was zero, then it was revised to 0.01; if the miRNA gene expression in both libraries was less than 1, owing to its too low expression, it was excluded from further differential expression analysis. Fold change = log2(the normalized H 2 O 2 treatment reads/the normalized control reads). The P-value was calculated as described previously (Wu et al., 2014).

Target Gene Prediction of miRNAs and Functional Analysis
Target genes were predicted using the MIREAP program developed by the BGI, combined with psRNATarget online software (http://www.plantgrn.org/psRNATarget/) (Dai and Zhao, 2011), and obeying the rules described in Allen et al. (2005) and Schwab et al. (2005). The criteria for using MIREAP and psRNATarget followed a previous study (Bertolini et al., 2013). Only the shared predictions of the two softwares were considered as the final target genes. The biological processes, molecular functions, and cellular components of the target genes were examined using the agriGO online tool (Du et al., 2010) to perform Gene Ontology (GO) annotation and GO enrichment analysis. The statistical test method was set as Fisher and the multi-test adjustment method was set as Bonferroni. The threshold of significance was defined as p < 0.01 and the false discovery rate (FDR) as <0.01. The dataset containing protein sequences of B. distachyon genome was set as the background dataset.

H 2 O 2 -Responsive miRNA-Mediated Network Analysis
The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database of physical and functional interactions (Szklarczyk et al., 2011) was used to analyse the protein-protein interactions (PPI) of the proteins encoded by target genes of differentially expressed miRNAs. The miRNA-regulated PPI network was displayed by the Cytoscape (version 3.1.1) software (Shannon et al., 2003).

Expression Validation of miRNAs and Their Targets
We verified the patterns of expression of seven conserved B.
distachyon miRNAs (miR159a-3p, miR159b-3p.1, miR160a/b/c/d-5p, miR169b, miR169d, miR397a, and miR528-5p). A miRcute miRNA First-strand cDNA Synthesis Kit (TIANGEN) was used for the RT reactions. The thermocycling program was adjusted to run for 60 min at 37 • C, 5 s at 85 • C, and then 4 • C forever. For each miRNA, three biological replicates were used. The miR168 and 5.8S genes served as the endogenous controls (Bertolini et al., 2013;Jeong et al., 2013). All primers are listed in Table S1. RT-qPCR was conducted on a CFX96 Real-Time PCR Detection System (Bio-Rad). Each reaction included 2 µL of product from the diluted RT reactions, 1.0 µL of each primer (forward and reverse), 12.5 µL of SYBR R Premix Ex Taq TM (Perfect Real Time; TaKaRa), and 8.5 µL of nuclease-free water. The reactions were incubated in a 96-well plate at 95 • C for 30 s, followed by 40 cycles of 95 • C for 5 s, 60 • C for 30 s, and 72 • C for 10 s. All reactions were run in triplicate for each sample. All data were analyzed using the CFX Manager software (Bio-Rad). We also selected 13 target genes to validate their expression profiles in the CS and TS libraries via RT-qPCR following the method described by Lv et al. (2014b) according to the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines. The Actin and SamDC genes served as the endogenous controls (Lv et al., 2014b). All primers are listed in Table S1. Statistical analysis was performed using the SPSS 17.0 software. Statistical differences amongst the two libraries were assessed using the independent two-sample t-test. P < 0.05 were considered statistically significant.
Modified RNA Ligase-Mediated (RLM) 5 ′ Race for the Mapping of mRNA Cleavage Sites To identify cleavage sites in the target mRNAs, a modified RLM-5 ′ -RACE was performed using a FirstChoice RLM-RACE Kit (Ambion, Austin, TX, USA). All the steps followed the manufacturer's instructions, except that the calf intestinal phosphatase treatment was omitted to maintain the cleaved transcripts. Nested PCR amplifications were performed using the general sense primers and gene specific nested antisense primers that were listed in Table S1. The amplification products were gel purified, cloned, and sequenced, and 10 independent clones were sequenced.

RESULTS
In this work, the Illumina Solexa sequencing platform was used to investigate the genome-wide identification and expression profiles of miRNAs in B. distachyon under H 2 O 2 stress. Two sRNA libraries were constructed using total RNAs isolated from control seedlings (CS) and H 2 O 2 -treated seedlings (TS). sRNA sequencing yielded a total of 18,220,106 and 19,373,978 high-quality raw sequence reads from the CS and TS libraries, respectively. The raw reads of the two libraries were uploaded to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA; accession numbers: SRX1542485 and SRX1542460). After removing low quality reads, adapters, poly-A sequences, and short RNA reads smaller than 18 nucleotides (nt), 17,811,109 (97.76%) and 17,708,762 (91.40%) clean reads representing 3,830,474 and 3,548,088 unique sRNAs were obtained from the CS and TS libraries, respectively (Table S2). Among the unique sequences, 1,963,866 (51.27%) and 1,689,996 (47.63%) generated from the CS and TS libraries, respectively, were mapped to the B. distachyon genome using SOAP2. To reveal the sequence distribution of the sRNAs, all clean reads were queried against the B. distachyon genome database at Phytozome (http://www.phytozome.net/), Rfam (http://rfam.sanger.ac.uk/), and miRBase v20.0 (http://www. mirbase.org/), and classified into seven annotation categories: non-coding RNAs (tRNA, rRNA, snRNA, and snoRNA), miRNA, exon-sense, exon-antisense, intron-sense, intron-antisense, and unknown sRNAs ( Table S2). The length distribution of the total sRNA reads revealed that the majority of reads from each library were 20-25 nt in length, of which, 24-nt reads were the most abundant, followed by 21-nt reads ( Figure 1A). Compared to the CS library, the TS library contained more 21-nt sRNAs and fewer 24-nt sRNAs. Of the unique sRNAs, 24-nt sRNAs accounted for 47.91 and 43.78% in the CS and TS libraries, respectively, while 21-nt sRNAs accounted for 5.10 and 5.33% ( Figure 1B).

Identification of Known and Novel miRNAs in B. distachyon
To identify known miRNAs from the CS and TS libraries, sRNA sequences generated from each library were independently aligned with currently known and experimentally validated mature B. distachyon miRNAs deposited in miRBase v20.0. Finally, a total of 144 known miRNAs were identified from the CS (133) and TS (138) libraries (Table S3); among them, 127 (88.19%) were detected in both libraries. For the total reads of each known miRNA in our study, 54.86% (79 miRNAs) had more than 100 reads, 29.17% (42) had more than 1000 reads, 12.5% (18) had more than 10,000 reads, and only 3.47% (5) had more than 100,000 reads. Among them, bdi-miR168-5p possessed the highest expression level in each library (556,771 reads in the CS library and 571,620 reads in the TS library), followed by bdi-miR156c, bdi-miR156b-5p, bdi-miR156d-5p, and bdi-miR528-5p.
After identifying the known miRNAs, the remaining sequences from the two libraries, which were classified as "unannotated" (excluding known miRNAs and other noncoding RNAs classified by Rfam), were used to discover novel and potential B. distachyon-specific miRNA candidates. To accomplish this, these small RNA sequences were aligned with the B. distachyon genome to identify genomic regions potentially harboring pre-miRNA sequences, whose hairpin-like structures are widely used to distinguish miRNAs from other small noncoding RNAs. The minimum free energy (MFE) of the secondary structures was also considered as a criterion for the prediction of potential pre-miRNAs. After aligning these unannotated sequences to the genome, we obtained a total of 221 novel miRNA candidates from the CS (156) and TS (132) libraries (Table S4). Among the 221 novel miRNA candidates, 36 had homologs in other plant species or their pre-miRNA sequences belonged to known B. distachyon pre-miRNAs (Table S4), while the remaining 185 novel miRNA candidates were B. distachyonspecific (Table S4). In agreement with previously reports, the cytosine and uracil nucleotides were dominant in the first position of the 5 ′ end for the majority of these newly determined putative novel miRNAs ( Figure S1A). In detail, first positions included 11,629 cytosine nucleotides (67.52%) and 4279 uracil nucleotides (24.85%) in the CS library and 16,115 cytosine nucleotides (78.05%) and 3455 uracil nucleotides (16.73%) in the TS library. The first nucleotide bias analysis showed that cytosine was the most frequently used first nucleotide for novel 21-nt miRNAs and uracil was most frequent for novel 20-, 22-, and 23-nt miRNAs ( Figure S1B). Our sequence analyses of the two libraries showed that the putative pre-miRNAs of each library greatly varied in length from 69 to 361 nt in the CS library and from 65 to 372 nt in the TS library. These pre-miRNA sequences were applied to predict the characteristic stem-loop secondary structure of pre-miRNA and their locations were also determined in the genomic loci (Table S4). We also calculated the minimum folding free energies of putative miRNA precursors for each library, which ranged from −18.1 to −180.64 kcal/mol with an average of −62.07 kcal/mol for the CS library and from −19.9 to −192.4 kcal/mol with an average of −60.46 kcal/mol for the TS library (Table S4). In contrast with the known miRNAs, most of the predicted novel miRNAs were expressed at very low levels. Only 9.50% (21) of the 221 novel miRNAs had more than 100 reads, and only two (novel_mir_59 and novel_mir_54) had more than 1000 reads (Table S4).

Target Gene Prediction and GO Annotation Analysis
A total of 352 putative known miRNA target transcripts (corresponding to 284 target genes) and 554 putative novel miRNA target transcripts (corresponding to 460 target genes) were obtained ( Table S5). The number of targets for each known miRNA and novel miRNA varied, ranging from 1 to 54 and 1 to 61, respectively, and the percent of novel miRNAs with more than 10 predicted target transcripts was 33.08%, while that number was 7.34% for known miRNAs. For comprehensive annotation, all putative target genes were analyzed by GO terms with the aid of the Blast2GO program with default parameters. Genes with a known function were categorized by biological process, molecular function, and cellular component according to the ontological definitions of the GO terms ( Figure S2). For biological process, genes were mainly in the single-organism process (10.12%), response to stimulus (9.97%), localization (6.49%), multicellular organismal process (6.17%), and developmental process (6.01%) categories. For molecular function categories, nucleic acid binding transcription factor activity (4.11%), transporter activity (3.01%), and structural molecule activity (1.42%) were highlighted. For cellular component, they were mainly localized in the organelle (26.90%), membrane-bounded organelle (25.00%), and membrane (13.61%).

Screening of H 2 O 2 -Responsive miRNAs
In this study, 31 known and 30 novel miRNAs were observed with a more than two-fold change in response to H 2 O 2 treatment in B. distachyon seedlings (Figure 2, Tables 1, 2). As reported in a previous study (Li et al., 2011), a series of known H 2 O 2 -responsive miRNAs, including miR159, miR160, miR169, miR397, and miR528, were also identified in our study. Further analysis revealed that 10 of the 31 known miRNAs were downregulated in the TS library compared to the CS library, whereas 21 were upregulated (Figure 2A). The 10 downregulated known miRNAs were composed of three miRNA families, including five members of miR169, four members of miR160, and miR7770. Among them, miR7770 displayed a dramatic (log2 fold change = −7.92) decrease. For the 21 upregulated known miRNAs, miR395 was the major family containing 14 members. bdi-miR159b-3p.1 showed the highest upregulation (log2 fold change = 16.21).
Among the 30 differential novel miRNAs, 17 downregulated and 13 upregulated miRNAs were found in the TS library compared to the CS library ( Figure 2B). All but three novel miRNAs (novel_mir_4, novel_mir_148, and novel_mir_161) showed dramatic changes (log2 fold change >4 or < −4), in contrast to known differential miRNAs, of which only four exhibited dramatic changes.
We performed GO enrichment analysis for the target genes of these H 2 O 2 -responsive miRNAs and found that the proteins encoded by these target genes were mainly involved in the categories of multicellular organismal development

A miRNA-Mediated Regulatory Network for H 2 O 2 Stress Response
Based on the miRNA target prediction and protein-protein interaction (PPI) analysis, the H 2 O 2 -responsive miRNAs and the proteins encoded by their target genes were used to construct the miRNA-mediated regulatory network for H 2 O 2 stress response (Figure 4). Generally, the more proteins that a protein interacts with, the more important are that protein's functions in the network. Two MYB TFs (encoded by genes Bradi2g53010 and Bradi1g36540) were centered in the network and each could interact with 24 proteins. Thus, the two MYB TFs may play critical roles during H 2 O 2 -triggered oxidative stress response.
Some target genes were regulated by two or more miRNAs that belonged to the same miRNA family or different families, and the proteins encoded by these genes always interacted with more proteins. For example, Bradi2g52840 (target gene of the bdi-miR395 family) encodes disease resistance protein RGA4, which interacts with seven proteins encoded by the target genes of H 2 O 2 stress-responsive miRNAs, including the two MYB TFs mentioned above (Figure 4). Gene Bradi2g55497 is regulated by three novel miRNAs (novel_mir_98, novel_mir_120, and novel_mir_222) and encodes transcription initiation factor TFIID subunit 12, which can interact with four proteins, including three members of the NFYA TF family (encoded by      genes Bradi4g01380, Bradi1g11800, and Bradi3g57320) and a transcriptional regulator SAGA-associated factor 29 homolog (encoded by gene Bradi4g08140).

RT-qPCR Validation of B. distachyon miRNAs and Target Genes
We applied RT-qPCR analysis for further experimental verification of the presence of several conserved miRNAs and compared the expression patterns of these miRNAs with deep sequencing results. Analysis of seven H 2 O 2 -responsive miRNAs (miR159a-3p, miR159b-3p.1, miR160a/b/c/d-5p, miR169b, miR169d, miR397a, and miR528-5p) by RT-qPCR ( Figure 5) showed that all of the relative expression profiles exhibited the same trends as their deep sequencing results, although there were some extent differences between the results obtained from deep sequencing and the RT-qPCR experiment. In detail, miR528-5p, miR397a, miR159a-3p, and miR159b-3p.1 were up-regulated under H 2 O 2 stress, while the expression of miR160a/b/c/d-5p, miR169d, and miR169b decreased during H 2 O 2 treatment.
RT-qPCR was also used for detection and quantification of the predicted targets of 13 H 2 O 2 -responsive miRNAs (Figure 6). Our results revealed that for most of the miRNAs, there was a negative correlation between the miRNA level and the levels of their target genes (Figure 6), with the exception of Bradi1g36540. Bradi1g36540 and Bradi2g53010 were both regulated by the upregulated bdi-miR159a-3p and bdi-miR159b-3p.1. Bradi2g53010 was downregulated, but Bradi1g36540 showed no significant change under H 2 O 2 stress. Bradi2g55497 was mediated by three novel miRNAs (Figure 4). Two of the miRNAs (novel_mir_98 and novel_mir_120) were downregulated, while another (novel_mir_222) was upregulated, but the expression level of Bradi2g55497 was unexpectedly downregulated. This phenomenon was also observed in the profiles of Bradi4g10380 (the target gene of novel_mir_98, novel_mir_120, and novel_mir_222) and Bradi1g14810 (the target gene of novel_mir_98 and novel_mir_222).

Verification of miRNA-Guided Cleavage of Target mRNAs in B. distachyon
To verify the miRNA-guided target cleavage, RLM-5 ′ -RACE experiment was performed to detect cleavage product of 3 known (bdi-miR160b-5p, bdi-miR159b-3p.1, and bdi-miR397a) and two novel bdi-miRNAs (novel_mir_152 and novel_mir_222). As shown in Figure 7, all five of the B. distachyon miRNAs guided the target cleavage, often at the tenth nucleotide (Figure 7).

DISCUSSION
To examine the H 2 O 2 -responsive miRNAs, two sRNA libraries were constructed from a mixture of 12-day old B. distachyon seedlings treated with 20 mM H 2 O 2 and a control sample, and were then subjected to next-generation deep sequencing. We identified a set of known and novel H 2 O 2 -responsive miRNAs with notable expression pattern changes and we also provide a miRNA-mediated oxidative stress response PPI network. These results provide useful information for elucidating the response and defense mechanisms for H 2 O 2 stress at the posttranscriptional level in plants.
Given that the expression of target genes is negatively regulated by miRNAs, the expression patterns of target genes always show an inverse correlation with those of miRNAs. However, our miRNA-target network showed that the regulatory mechanism of miRNA may be more complex. Some miRNAs can each regulate several target genes with the same function or different functions, and some genes can be regulated by several miRNAs. Among the known differential miRNAs in our study, some regulated two or more target genes with the same function, such as bdi-miR159a and bdi-miR169d/k, whereas more miRNAs, such as bdi-miR395, bdi-miR408, and bdi-miR528, regulated two or more target genes with different functions (Figure 4 and Table 1). Among the novel differential miRNAs, novel_mir_152 was downregulated (log2 fold change = −8.13) under H 2 O 2 stress and mediated four target genes encoding a cell differentiation protein RCD1 homolog associated with RNA degradation and cell differentiation. The 20 potential genes regulated by novel_mir_120 are involved in various functions. Gene Bradi2g55497, which encodes transcription initiation factor TFIID subunit 12, was regulated by three novel miRNAs (novel_mir_98, novel_mir_120, and novel_mir_222). Two of them (novel_mir_98 and novel_mir_120) were downregulated under H 2 O 2 stress, while novel_mir_222 was upregulated. Compared with novel_mir_98 and novel_mir_120, novel_mir_222 may have the opposite effect on the expression of Bradi2g55497. Furthermore, RT-qPCR analysis demonstrated that the final result of the FIGURE 4 | A miRNA-mediated H 2 O 2 -responsive regulatory network. Triangles represent the differentially expressed miRNAs and the color gradient shows the fold change (red: upregulation; green: downregulation). Blue circles represent the target genes of differentially expressed miRNAs; yellow circles represent target genes that encode transcription factors.
antagonism is the downregulation of Bradi2g55497 (Figure 4), so novel_mir_222 may play a main role under H 2 O 2 stress. Four of the target genes of downregulated novel_mir_120 and upregulated novel_mir_222 are the same, which seems to imply that they are a pair of miRNAs with opposing effects. Additionally, novel_mir_91 and novel_mir_197 also exhibited the same pattern. Previous studies also found similar phenomenon. For example, genes SPL and AP2, encoded DNAbinding transcription factors, can be regulated by miR156 and miR172, which showed opposite expression patterns (Ding et al., 2014).
The upregulation of a miRNA always results in the aggravated degradation of its target gene in plants, and vice versa. These target genes of upregulated miRNAs may be associated with stress response, whereas the target genes of downregulated miRNAs are always involved in stress resistance. miR528, a monocot-specific miRNA, was upregulated under H 2 O 2 stress (log2 fold change = 2.88). bdi-miR528 is also upregulated FIGURE 5 | RT-qPCR validation of the miRNAs. Blue bar: relative gene expression level in the control library. Red bar: relative gene expression level in the H 2 O 2 stress library. The data are derived from three biological repeats and represent mean ± standard deviation (n = 3). "*" and "**" indicate significant difference at P < 0.05 and 0.01 level, respectively.
in Bd21 under drought stress (Bertolini et al., 2013) and in rice under H 2 O 2 stress . Bradi3g19170, the target gene of bdi-miR528, codes the E3 ubiquitin-protein ligase XBOS35, which may be involved in ubiquitination for proteasomal degradation and the ethylene (ET) signaling pathway during abiotic stress response (Sobeih et al., 2004;Prasad et al., 2010). Thus, miR528 may be a vital miRNA involved in water and oxidative stress response in monocots. As the largest MIR gene family in B. distachyon, the MIR395 gene family includes 15 members. In our study, 13 members of the MIR395 gene family were overexpressed under H 2 O 2 stress and their target genes were Bradi1g09030 (encoded ATP-sulfurylase 1), Bradi1g24110 (encoded Ribosome-recycling factor), and Bradi2g52840 (encoded disease resistance protein RGA4), which may be involved in oxidative stress response. Similar to our results, all of the identified MIR395 gene family members are also upregulated under drought stress in B. distachyon (Bertolini et al., 2013). All of the target genes of bdi-miR397a, except Bradi3g03407, encode laccase, which plays an important role during the formation of lignin in the cell wall (Mayer and Staples, 2002); a previous study also confirmed that miR397 is a negative regulator of laccase genes (Lu et al., 2013). bdi-miR827-3p was upregulated under H 2 O 2 stress, and a previous study (Jeong et al., 2013) found that this miRNA was also upregulated under phosphate starvation conditions in shoots. In contrast to our results, a previous study showed that miR397 and miR827 were downregulated in rice by H 2 O 2 treatment, which indicates that the miRNA response to the same abiotic stress occurs in a genotype/species-dependent manner (Zhang, 2015).
In plants, miRNAs participate in the response to abiotic stress by mediating key components of complex gene networks (Zhang, 2015). Many of the target genes of miRNAs identified in our study are TFs, which is consistent with previous studies (Jones- Rhoades et al., 2006;Baev et al., 2011;Zhang, 2015). Thus, a miRNA can indirectly regulate the expression of downstream target genes through regulation of the expression of TFs. Previous studies revealed that the mechanism of the MYB TF involved in gibberellic acid (GA) signal transduction was FIGURE 6 | RT-qPCR analysis of the miRNA target genes. The Actin and SamDC genes were served as the endogenous controls. Blue bar: relative gene expression level in the control library. Red bar: relative gene expression level in the H 2 O 2 stress library. The data are derived from three biological repeats and represent mean ± standard deviation (n = 3 ). "*" and "**" indicate significant difference at P < 0.05 and 0.01 level, respectively. a previous study in potato showed that GRAS family TF, involved in development and stress responses, such as drought stress (Ma et al., 2010), is the putative target gene of miR171 (Hwang et al., 2011). The pre-miRNA of novel_mir_162 was the same as the precursor of bdi-miR395c, which may be a novel member of the bdiMIR395 family. Similar to other known bdi-miR395 family members, this miRNA was also upregulated ( Table 2). In addition, novel_mir_40 was identified as a novel member of the bdiMIR169 family, which was also downregulated (log2 fold change = −10.28), as were the five known differentially expressed bdiMIR169 family members ( Table 1). In contrast to typical members of the bdiMIR169 family, whose target genes always encode NF-YA TFs, the target gene of novel_mir_40 is Bradi1g25517, which encodes glucan endo-1,3-beta-glucosidase 3.
H 2 O 2 can trigger various phytohormone signaling pathways involved in abiotic and biotic stress responses and there are complex crosstalks among different phytohormone signaling pathways (Harrison, 2012;Saxena et al., 2016). Recent study found out some immune hormone marker genes involved in salicylic acid (SA), jasmonic acid (JA) and ET signaling pathways in B. distachyon (Kouzai et al., 2016). In our study, we also found some genes related with phytohormone signaling pathways were the target of miRNAs. For example, the two MYB genes involved in GA pathway were regulated by bdi-miR159a/b, ARF22, and ARF5 genes involved in auxin pathway were the target of bdi-miR160 and bdi-miR397a, respectively, and the E3 ubiquitinprotein ligase XBOS35 and ethylene-responsive transcription factor RAP2-3 genes related with ET signaling pathways were mediated by bdi-miR528 and novel_mir_120/novel_mir_222, respectively.

CONCLUSION
In this study, we sequenced and analyzed the sRNA of model plant Bd21 seedlings under H 2 O 2 stress and normal growth conditions using large-scale sequencing technology. Finally, we identified 39 known and 221 novel miRNAs, of which 31 known miRNAs and 30 novel miRNAs were involved in H 2 O 2 -stress response and resistance. Moreover, RT-qPCR analysis of several representative miRNAs and their target genes and cleavage site analysis through 5 ′ RACE validated our sequencing and bioinformatic results. The PPI network mediated by miRNAs revealed the regulation mechanism of signal transduction and oxidative stress resistance under H 2 O 2 stress. Further analysis of the differentially expressed miRNAs and their target genes will help us understand the mechanism of oxidative stress response and tolerance in plants.

AUTHOR CONTRIBUTIONS
DL carried out all experiments, data analysis and wrote the manuscript. SZ contributed the RLM-5 ′ -RACE experiment. GZ, YB, and GC performed the RNA extraction and RT-qPCR, CH and ZY conducted GO annotation and enrichment analyses. DL and YY conceived the study, participated in the design and coordination, and in interpretation of the dataset. All authors read and approved the final manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 01567  Table S1 | All the RT-qPCR and RACE primers used in this study.