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₂O₂) 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 (O2−), another kind of ROS, can also be rapidly dismutated to H₂O₂ spontaneously or catalyzed by superoxide dismutase (SOD). In contrast to other ROS, such as O2− and hydroxyl radicals (OH⁻), H₂O₂ 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; Neill S. et al., 2002; Neill S. J. et al., 2002; Vranová et al., 2002) and an indispensable second messenger in biotic and abiotic stress responses (Pastori and Foyer, 2002). H₂O₂ stress can affect fluctuation of the Ca²⁺ concentration in plants, thereby inducing the production of appropriate amounts of antioxidants (Rentel and Knight, 2004). Global transcript profiling under H₂O₂ 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₂O₂ can also trigger acclimation and cross-tolerance phenomena (Neill S. et al., 2002; Pastori and Foyer, 2002). Overall, at low concentrations H₂O₂ may play a role as a signaling molecule, whereas at high concentrations H₂O₂ 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₂O₂ can also act as a secondary messenger during stress response and defense signal transduction. However, only a few H₂O₂-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. distachyon-specific 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₂O₂ 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₂O₂-treated conditions, as well as the miRNA-directed regulatory network, provide new insights that will inform the genetic improvement of stress tolerance in plants.

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₂O₂ stress. Two sRNA libraries were constructed using total RNAs isolated from control seedlings (CS) and H₂O₂-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).

To examine the H₂O₂-responsive miRNAs, two sRNA libraries were constructed from a mixture of 12-day old B. distachyon seedlings treated with 20 mM H₂O₂ and a control sample, and were then subjected to next-generation deep sequencing. We identified a set of known and novel H₂O₂-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₂O₂ stress at the post-transcriptional 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₂O₂ 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₂O₂ 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 antagonism is the downregulation of Bradi2g55497 (Figure 4), so novel_mir_222 may play a main role under H₂O₂ 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 DNA-binding 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₂O₂ stress (log2 fold change = 2.88). bdi-miR528 is also upregulated in Bd21 under drought stress (Bertolini et al., 2013) and in rice under H₂O₂ stress (Li et al., 2010). 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₂O₂ 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₂O₂ 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₂O₂ 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 regulated by miR159 (Achard et al., 2004; Alonso-Peral et al., 2010). In our study, bdi-miR159a/b (target genes Bradi2g53010 and Bradi1g36540, two MYB genes) was upregulated under H₂O₂ stress, which is consistent with the results in rice and B. distachyon under drought stress (Zhou et al., 2010; Bertolini et al., 2013). Auxin response factors (ARFs) are a class of auxin-responsive TFs mediated and regulated by miR160 (Mallory et al., 2005). miR160 regulates the expression of ARF genes by combining with the complementary sites of a non-coding region of ARF genes (Rhoades et al., 2002). In our study, gene ARF22 was targeted by four downregulated members of the bdi-MIR160 family (Table 1). Interestingly, gene ARF5 was targeted by the upregulated bdi-miR397a, instead of bdi-miR160. In addition, the expression of some miRNAs can also be regulated by TFs through specific binding to the promoter region of the miRNA. For example, miR160 can be transcriptionally regulated by proteins ARF6 and ARF17 in Arabidopsis (Gutierrez et al., 2009). Thus, there are fine tuning mechanisms in plants to modulate gene expression through miRNA-TF-mediated feedback loops (Meng et al., 2011; Iyer et al., 2012). bdi-miR169 was significantly downregulated under H₂O₂ stress in our study and a previous study in rice (Li et al., 2010); its target gene encodes a nuclear TF Y subunit A-2 (NF-YA-2) or NF-YA-4, which belongs to the CCAAT binding TF family. In Arabidopsis, the upregulation of NF-YA gene depends on the regulation of miR169 under drought stress (Li et al., 2008). Thön et al. (2010) found that CCAAT-binding TFs are involved in the response to oxidative stress. Therefore, miR169 may play a crucial role in the process of resistance to oxidative stress in plants. The network analysis showed that NF-YA could interact with V-type proton ATPase subunit d (encoded by Bradi2g42100), transcription initiation factor TFIID subunit 12 (encoded by Bradi2g55497), 2-oxoglutarate-dependent dioxygenase DAO (encoded by Bradi2g24715), and TF MYB (encoded by Bradi2g53010 and Bradi1g36540). One of the target genes of the upregulated bdi-miR408-5p encodes TF TCP15, whose activity is inhibited by oxidative stress (Viola et al., 2013). Overexpression of miR408 in Arabidopsis can improve its tolerance to oxidative stress (Ma et al., 2015).

Among the differentially expressed novel miRNAs, novel_mir_160, a potential new member of the bdiMIR171 gene family, was downregulated (log2 fold change = −2.12) under H₂O₂ stress. miR171 is also significantly downregulated under drought stress in rice and potato (Zhou et al., 2010; Hwang et al., 2011). Although we did not identify the target genes of this potential novel bdiMIR171 gene family member, 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₂O₂ 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 ubiquitin-protein 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.

In this study, we sequenced and analyzed the sRNA of model plant Bd21 seedlings under H₂O₂ 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₂O₂-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₂O₂ 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.

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.