Boron Stress Responsive MicroRNAs and Their Targets in Barley

Boron stress is an environmental factor affecting plant development and production. Recently, microRNAs (miRNAs) have been found to be involved in several plant processes such as growth regulation and stress responses. In this study, miRNAs associated with boron stress were identified and characterized in barley. miRNA profiles were also comparatively analyzed between root and leave samples. A total of 31 known and 3 new miRNAs were identified in barley; 25 of them were found to respond to boron treatment. Several miRNAs were expressed in a tissue specific manner; for example, miR156d, miR171a, miR397, and miR444a were only detected in leaves. Additionally, a total of 934 barley transcripts were found to be specifically targeted and degraded by miRNAs. In silico analysis of miRNA target genes demonstrated that many miRNA targets are conserved transcription factors such as Squamosa promoter-binding protein, Auxin response factor (ARF), and the MYB transcription factor family. A majority of these targets were responsible for plant growth and response to environmental changes. We also propose that some of the miRNAs in barley such as miRNA408 might play critical roles against boron exposure. In conclusion, barley may use several pathways and cellular processes targeted by miRNAs to cope with boron stress.


Introduction
MicroRNAs (miRNAs) are a class of single strand, endogenous, non-coding small RNA molecules, which post-transcriptionally regulate gene expression in many organisms by targeting mRNAs for cleavage or translation suppression [1,2,3]. Increasing evidence demonstrates that miRNAs play an important role in many biological and metabolic processesincluding regulation of plant growth, development and response to biotic and abiotic stresses via interactions with their specific target mRNAs [4,5,6,7,8]. Boron (B) is an essential element for plants, and its deficiency generally causes growth defects mainly in young and growing parts of the plants, while excessive levels of B are toxic to plants [9,10]. A number of physiological processes are shown to be altered by B exposure. Deterioration of cell wall biosynthesis, metabolic deterioration by binding to the ribose moieties of ATP, NADH and NADPH, and inhibition of cell division and elongation are the most distinct symptoms of B toxicity [11,12,13]. However, plants also evolve mechanims to cope with the presence of excessive amounts of metal ion. Although several studies have been performed on small RNAs and metal stressors such as mercury, cadmium, and aluminum [14,15,16], no studies have been reported on boron stress. Barley (Hordeum vulgare L.) is one of the most important grain crops grown and cultivated worldwidely [17]. Additionally, it is a well-studied model plant for triticacea research in terms of genetics, genomics, and breeding [18,19]. Although miRNAs in barley were identified in previous studies [19,20,21,22], compared with the number of identified miRNAs in other grain crops such as rice and maize, the number of known miRNAs in barley is still very insufficient. Initially, conventional approaches were extensively used for miRNA identification and contributed considerably to the miRNA exploration [8,23].
The purpose of this study is to identify tissue specific expression of miRNAs and their potential targets in barley exposed to high levels of boron. To achieve this goal, we identified miRNAs from the entire transcriptome RNA-seq data, which included more than 208 million reads generated from control and B-exposed roots and leaves of B-tolerant barley seedlings. Some of the identified barley miRNAs were validated in leaf and root tissues by quantitative RT-PCR. Additionally, 'degradome sequencing' approach was also employed for miRNA target identification in barley.

Plant Materials and Boron Treatment
Barley (Hordeum vulgare L. cultivar Sahara) seeds were sterilized and placed into Petri dishes for germination at room temperature. Then, four-day-old seedlings were transferred into liquid culture flasks including nutrient solutions. The treatments were repeated at least three times with triple biological replicates. For toxicity experimets, toxic (1000 mM) and nontoxic (50 mM) concentrations of B were added to different flasks. Germinated seedlings were exposed to B-toxic or B-nontoxic conditions for approximately 24 hours.

RNA Isolation, cDNA Library Construction and Sequencing for Transcriptome Analysis
Total RNAs were extracted from barley root and leaf tissues using the TRIZOL Reagent (Invitrogen) according to the manufacturer's instructions. The extractions were performed separately for each sample with three independent biological replicates and same amount of total RNA was subsequently pooled based on their concentration. The quality and quantity of purified RNAs were assessed with a Nanodrop 2000c spectrophotometer (Nanodrop Technologies, USA) and the presence of ribosomal RNA bands was determined by Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). All RNA samples were stored at 280uC until further processing. The cDNA library construction and Illumina (Solexa) based-transcriptome sequencing experi- ments were conducted by the BGI (Beijing Genomics Institute, Hong Kong). In brief, for each library, the polyadenylated RNA (mRNAs) was isolated from 20 mg of each RNA pool using oligo(dT) 25 magnetic beads (Invitrogen) according to the manufacturer's protocol. Following purification step, the isolated mRNAs were fragmented into small pieces using fragmentation buffer. These mRNA fragments resulted were used as templetes for first strand cDNA synthesis by reverse transcriptase and random hexamer primers. After completion of first-strand synthesis, the second-strand cDNA was synthesized by using DNA polymerase I, RNaseH, dNTPs and buffer. Subsequently, newly synthesized short fragments were purified with a QIAquick PCR extraction kit and further subjected to end reparation and adding poly(A) using EB buffer, T4 DNA polymerase, the Klenow fragment, and T4 polynucleotide as well as Klenow 39 to 59 exopolymerase and then, the sequencing adapters were connected with those short fragments. After size selection and purification through agarose gel electrophoresis, the selected fragments were enriched by PCR amplification with appropriate primers and eventually, the library sequencing experiments were performed using Illumina HiSeq TM 2000 instrument.

De novo Assembly of Transcriptome and Data Processing
Firstly, the image data obtained from sequencing platform were converted by base calling into sequence data which is commonly known as raw reads and typically stored in the fastaq file format. In order to acquire high-quality clear reads, all raw sequence reads were filtered to remove adaptors, the reads containing unknown nucleotides larger than 5% and low quality reads. Then, remaining clear reads were subjected to de novo transcriptome assembly using Trinity assembler with default settings [24]. Briefly, Inchworm, one of software module of Trinity, assembles reads with definite length of overlap in order to generate longer fragments which are termed contigs. Then, the reads are are realigned to the newly formed contig so as to construct scaffolds which are basically derived from the contigs from the same transcript. By using paired-end information, these reads (pair-end clean reads) are mapped back to the resultant scaffolds to fill the intra-scaffolds gaps. Subsequently, the sequences generated by assembly of scaffolds are defined as unigenes which cannot be extended on either end by further assembly process. After obtaining non-redundant unigenes as long as possible using sequence clustering software, some of those unigenes were termed by singletons which are are not part of any contigs.

Identification of miRNAs and Prediction of miRNA Precursors (Pre-miRNAs)
The RNA-seq generated more than 208 million clean reads (Kekec et al. Unpublished data) that were used to identify miRNAs and their targets. To identify potential miRNA precursors (pre-miRNAs), Blastn search was perfromed with an e-value of 1e-10 between unigene reads constructed from a total of four libraries for whole transcriptome of H. vulgare and previously known pre-miRNA sequences obtained from miRBase v19.0 (http://www.mirbase.org/) and plant miRNA database (PMRD, http://bioinformatics.cau.edu.cn/PMRD/), respectively. Briefly, we created our own blast-search database from a fasta sequence file including all transcripts of H. vulgare to be used with Blastn algoritm. As a query, all plant pre-miRNAs were used to search against the database generated via Blastall-2.2.22. The output file was manually investigated to extract the longest possible set of matches from the batch sequences. Then, the chosen potential pre-miRNA sequences obtained from the result of Blastall were subjected to the Zuker folding algorithm for in silico secondary The mapped read counts of each pre-miRNAs were normalized in terms of the length of pre-miRNA and total read numbers according to RPKM method (Reads Per kb per Million reads) [61]. doi:10.1371/journal.pone.0059543.t002 structure generation via the web-based computational software MFOLD 3.2 [25]. The default parameters of the software were adjusted for predicting secondary structure of the selected sequences and the minimal folding free energy index (MFEI) was calculated for each pre-miRNA sequence as described previously [26]. After identification of putative pre-miRNAs, we determined the localization of predicted mature miRNAs on the their own pre-miRNAs by mapping these mature miRNAs to the pre-miRNAs using BLAST search algorithm with default parameters. To consider whether these sequences are potential miRNA and pre-miRNA candidates, the following empirical criteria was adopted: (i) a pre-miRNA sequence can fold into an appropriate stem-loop hairpin secondary structure which contains 19-24 nucleotides mature miRNA in lengths, (ii) the mature miRNAs should take place in one arm of the hairpin structure, (iii) the minimum length of pre-miRNA sequences should be 60 nucleotides, (iv) mismatch between candidate mature miRNAs and query mature miRNAs are allowed to be up to 3 nucleotides, (v) it is strongly emphasized that MFEI and negative minimum folding energy (MFE) of potential miRNA precursors are relatively higher than those of other types of RNAs, and MFEI was considered to be one of the most important criterion for distinguishing miRNAs from other types of RNAs, e.g. pre-miRNA with approximately .0.67 has been identified as more likely to be a miRNA [27], (vi) there is no large loop or break in the miRNA:miRNA*, and (vii) the miRNA has less than six mismatches with the opposite miRNA sequence (miRNA*) on the opposite arm [28,29,30].

Computational Identification of miRNA Target Genes
Identification of miRNA-regulated gene targets is crucial for understanding miRNA functions. Therefore, the putative targets of H. vulgare miRNAs were identified by aligning the miRNAs with the high-quality unigene set obtained from the assembled transcripts and the singleton transcripts of barley de novo transcriptome libraries using the web-based psRNA Target Server (http://plantgrn.noble.org/psRNATarget/) with default parameters of user-submitted option. Alignment between Hvu-miRNA and its potential target(s) was evaluated by the parameters described by Zhang [31]. These computationally identified miRNA targets were further analyzed using BlastX searches with an e-value of 1e-10 against Hordeum EST sequences at NCBI database to identify putative gene homologs for confirmation.

Gene Ontology (GO) Analysis of Potential miRNA Targets
In order to better understand the functional roles of miRNAs in barley, Blast2Go (B2G) software v2.3.1 was used to assign gene ontology (GO) annotations of predicted target genes [32]. First, all miRNA target transcripts were arranged in a text file (Fasta format) as an input data and uploaded to the program for BlastX searches with an e-value of 1e-06. The BLAST results included: sequence length, gene name, e-value, similarity, Hit-length, Alignlength, GenBank and Uniprot accession number as well as Gene Ontology IDs belonging to each target sequences. Next, the output file (.dat file) obtained from BlastX analysis was used to retrieve GO terms associated with each blast hit and Gene Ontology annotations. Finally, all miRNA targets representing genes with known function were categorized by biological process, cellular component and molecular function according to the ontological definitions of the GO terms.

miRNA Target Cleavage Product (Degradome) Analysis
In order to characterize the miRNA cleaved target library (degradome) of H. vulgare, we evaluated a dataset derived from the output generated by CleaveLand (v2.0) software (a pipeline for using degradome data to find cleaved small RNAs). The miRNAdirected cleavage site in the miRNA:mRNA alignment is represented by red arrow (Table S1).  (Table S2). The miRNA stem-loop reverse transcription was carried out using 500 ng of total RNA samples of B-nontoxic (50 mM) and B-toxic (1000 mM) leaf and root samples (1 mL), 0.5 mL 10 mM dNTP mix, 1 mL stem-loop RT primer (1 mM) and 10.5 mL nuclease free water. Those components were also mixed separately for the different dilutions of total RNA stem-loop RT primer for cDNA synthesis and incubated for 5 min at 65uC, and then put into ice for 2 min. Thereafter, 4 mL first-strand buffer (56), 2 mL 1 M DTT, 0.1 mL RNAseOUT (40 units/mL), and 0.25 mL SuperScript III (200 units/mL) were added to each tube. The RT reactions were fulfilled as 30 min at 16uC followed by 60 cycles of 30uC for 30 s, 42uC for 30 s and 50uC for 1 s. The control tubes included all components without RT primer (no RT or -RT) and RNA template (no RNA or -RNA).

Quantitative Real-time PCR
To verify some of the predicted H. vulgare miRNAs experimentally, and to measure and compare the expression levels of the miRNAs in root and leaf tissues treated with boron, qRT-PCR was conducted using SYBR Green I Master Kit (Roche, Germany) on a LightCyclerH 480 Real-Time PCR System (Roche, Germany). For qRT-PCR analysis,10 mL 2X Master mix, 0.1 mL (100 pmol) forward and 0.1 mL (100 pmol) reverse primers, 7.8 mL nuclease-free water and 2 mL RT stem-loop cDNA products were used. Forward primers were specifically designed for each individual miRNA, but 59-GTGCAGGGTCCGAGGT-39 was used as the universal reverse primer [33] (Table S3). The qRT-PCR conditions were as follows: initial denaturation at 95uC for 10 min, followed by 41 cycles at 95uC for 10 s, 55uC for 20 s, and 72uC for 10 s. The melting curves were adjusted as 95uC for   5 s and 55uC for 1 min and then cooled to 40uC for 30 s. All reactions were repeated three times [30]. For each conditions, the qRT-PCR experiments were run as biological triplicates and expression levels were normalized according to pervious studies [8,19,28,29,33,34]. The relative fold change for each comparison was calculated by 2 ' -DCt after normalization [33,34]. Error bars were derived from the three experiments in triplicate and error bars represent standard deviation.

Validation of Barley miRNA Target mRNAs by qRT-PCR
To verify the expression levels of identified 11 barley miRNAs, the mature miRNAs were measured via qRT-PCR. Relative expression levels of predicted barley miRNAs were compared in root and leaf tissues treated with excess boron. The expression profile of these miRNA targets was also measured using qRT-PCR and their specific PCR primers were listed in the Table S3. The reverse transcription reaction was performed with Transcriptor High Fidelity cDNA Synthesis Kit (Roche, Germany) according to the manufacturer's protocol. The qRT-PCR experiment was carried out as reported previously [34,35]. Briefly, 2 mL cDNA was amplified with 0.1 mL specific primers in a total volume of 18 mL using SYBR Green I Master (Roche) with LightCyclerH 480 Real-Time PCR System. 18s rRNA (GenBank ID: AF147501) amplified with forward: GTGACGGGTGACGGAGAATT and reverse: GACACTAATGCGCCCGGTAT primers were used as a reference gene with triple replicates [36,37].

Identification of Boron Responsive miRNAs in Barley
According to sequence similarity to known plant miRNAs, 31 known and 3 new miRNAs were identified. Previously, miR157, miR165, and miR319 have been identified in other plant species, but so far they have been undetermined in barley. Identified miRNAs in barley were located on either arm of the predicted pre-miRNA sequences. Of the 34 identified H.vulgare miRNAs, 47% of mature sequences were located in the 59 arm of pre-miRNAs, while 53% were situated in the 39 arm ( Fig. 1; Fig. S1). The majority of these miRNAs were 21 nt long, followed by 22 nt, 20 nt and 23 nt, respectively (Table 1), which is consistent with miRNAs from other plant species [23,36]. In addition, our study showed that the average of MFEI was 0.86, which is higher than that of other types of RNA molecules such as tRNAs (0.64), rRNAs (0.59) and mRNAs (0.62-0.66) (Table 1) [29,37,38,39].

Boron Stress Induced Aberrant Expression of miRNAs in Barley
To identify the response of barley miRNAs to B treatment, we compared the expression profile of miRNAs between treated and untreated groups. The reads were normalized on the basis of transcripts per million obtained from high-throughput sequencing ( Table 2). Several conserved miRNAs (such as miR160 and miR171) and non-conserved miR5141 were found abundantly in both libraries, but many others were detected with only a few in both libraries or could not be found in either library. We also found that some miRNAs are only expressed in either root or leaf tissues. The miR156c and miR319a were highly expressed in root, whereas miR408 was only detected in leaf. In addition, some miRNAs such as miR156, miR169, miR172, and miR1121 were highly expressed in root but miR2004 was highly expressed in leaf. Expression of most miRNAs was significantly changed in a tissuespecific manner under boron stress whereas the remaining miRNAs were found to be responsive in both tissues. In root tissue responding to boron stress, miR165, miR2004, and miR5051 were up-regulated whereas miR444b and miR2024a were down-regulated. miR156, miR169c, miR171, miR171a, miR444a, miR444c, miR2023a were up-regulated while miR156d, miR397, miR408, miR1121, miR2014, miR5049, miR5141, miR5180, and miR5180a were down-regulated in leaf tissue upon boron stress. In addition, some miRNAs, such as miR172, miR399, miR2021, miR5053 and miR5066 were expressed in both root and leaf (Table 3).

Target Identification of miRNAs in Barley Using Degradome Analysis
A total of 934 genes targeted by miRNAs were identified in barley by CleaveLand (v2.0) ( Table 4). However, we could not identify the cleavage signature for some of the known miRNAs. The miRNA guided cleavage sites by degradome analysis are shown in Fig. 2 and Table S1. According to the results of blastn analysis of the identified miRNA targets, many of the targets were homologous to conserved target genes existing in other plants species; these targets included squamosa promoter-binding protein, auxin response factor (ARF), MYB transcription factor family, AP-2 Transcription Factors, NAC transcription factor (NAC), AGO1, and class III homeodomain-leucine zipper (HD-ZIP III) proteins. Most of these targets were found to be responsible for plant growth and response to environmental changes. For example, the target transcript of miR168 was ARGONAUTE1 protein (AGO1) family protein, which functions in plant development and in response to stress stimulus, such as NaCl and mannitol stress in rice. [40].

qRT-PCR Validation and Measurement of H. vulgare miRNA Levels and their Targets
Eleven identified barley miRNAs and their targets were further investigated using qRT-PCR. Both conserved barley miRNAs (miR156, miR159, miR164, miR166, miR168, miR171, miR395 and miR396) and non-conserved barley miRNAs (miR1120 and miR5048) were detected. The expression levels of barley miRNAs and their targets were comparatively shown in Fig. 2. The miR159, miR164, miR166, miR171, and miR414 were induced in leaf, but were inhibited in root tissues exposed to boron stress.   Although miR168 was induced, miR159, miR396, miR1120 and miR5048 were inhibited in both root and leaf upon excess boron exposure. The targets of miR159 and miR1120 were found to be up-regulated in both root and leaf upon boron stress, but miR395 and miR5048 target genes were down-regulated in root but remained at the same levels in leaf tissue upon boron stress. Additionally, miR171 target gene was down-regulated in leaf but up-regulated in root upon boron stress (Fig. 2).

Gene Ontology (GO) Analysis
According to the gene ontology analysis, the predicted targets were classified into three main categories: biological processes, cellular components, and molecular functions (Table 5). Of these, cellular and metabolic process in biological process, cell and organell part in cellular component, and binding and catalytic activity in molecular function were the most established categories.

Discussion
High-throughput sequencing technology has currently been successfully applied to identify miRNAs at whole genome scale in several plant species, including: soybean [41,42], peanut [43,44], barley [21,45], poplar [46], olive [47], Medicago [48], grapevine [49], rice [50], and cucumber [23]. However, almost all of the previous studies have been performed under normal growth conditions, few are associated with stress conditions. Li et al. [41] reported soybean miRNAs under three stress treatments (drought, salinity, and alkalinity) via high-throughput sequencing. Drought stress responsive miRNAs shows differential expression in response to heat stress in Populus euphratica and wheat [46,51]. In this study, we constructed RNA libraries from barley leaves and roots treated with boron stress compared to control conditions to idnetify boron stress-responsive miRNAs in barley using high-throughput sequencing.
Boron treatment affected the expression profiles of miRNAs in barley leaf and root tissues. The most striking oneswith 16-fold and 12-fold changes were miR408 and miR5180, respectively. The remaining changes in the expression of miRNAs ranged between 2-to 7-fold (Table 3).
Recently, miR408 was identified in barley, which targets Cubinding domain containing chemocyanin and blue copper protein [19]. In this study, we found that miR408 also potenially targets heterotrimeric G protein alpha (a) subunit and ATPase family gene 1 (AFG1). Heterotrimeric G proteins and ATPase gene family plays significant roles in signal transduction pathways in plants [52,53,54]. Fujisawa et al. [55] reported that suppression of a subunit gene expression causes abnormal morphology in rice. In response to water deficit, miR398 and miR408 were induced in Medicago truncatula [56]. In addition, expression of miR408 upon drought stress in barley was found to be induced in leaves, but unchanged in roots [19]. However, in Oryza sativa, miR408 expression was reported as 2.76-fold down-regulated 12 days after water withholding at tillering stage upon drought stress using microarray analysis [40]. In our study, expression of miR408 was down-regulated significantly (16-fold) upon excess boron treatment in barley leaves.
Previous studies reported the miRNA expression in a speciesspecific or tissue-specific manner [57,58]. miR168, miR319, miR396, and miR397 were induced by drought in Arabidopsis thaliana but were suppressed in Oryza sativa [59]. Additionally, the expression of miR399 was induced in shoots upon phosphate deficiency treatment, but it was accumulated in both shoots and roots [57]. In barley, miR166 was up-regulated in leaves, but was down-regulated in roots; miR171 level was induced in leaves, but it was not affected in roots [19]. In our study, miR169c, miR171, and miR399 were up-regulated in leaves whereas miR397, miR444b were down-regulated in roots after exposure to high B concentration. The miR172 was down-regulated 2-fold in roots but up-regulated 4-fold in leaves in response to boron stress. The miR169c and miR171 was determined to be 6-fold up-regulated and 2-fold up-regulated in leaves under boron stress, respectively. In Medicago truncatula, miR169 and miR172 were up-regulated but miR171 and miR390 were down-regulated upon mercury exposure [16]. Similarly, miR171 was down-regulated but miR172 was up-regulated by cadmium exposure in Brassica napus [15]. However, in response to Al 3+ treatment, miR171 was upregulated in Medicago truncatula [14].
Our study demonstrated that boron stress inhibited miR156a expression in barley leaves. However, we did not detect its expression in roots. In addition, the target of miR156a, SBP protein gene, was down-regulated in stressed leaves, but was unaltered in roots in response to boron stress (Fig. 2). This result is similar to the prvious report [19], whereas not affected in roots upon drought stress. Expression of miR156 has been investigated in many studies as down-regulated in Oryza sativa, Zea mays, Populus tremula, Populus trichocarpa in response to drought stress, salt stress, cold stress, mechanical stress, while up-regulated in Arabidopsis thaliana, Triticum aestivum, Nicotiana tabacum upon salt stress, heat stress, viral infection, respectively [59]. Our study indicates that miR156 was also boron stress responsive in leaves upon excess boron treatment.
For better understanding of the functions of miRNAs, gene ontology analysis for miRNA target transcripts was performed. Sixty genes targted by 34 miRNAs were found to be involved in 77 biological processes. These major processes are as follows: biological regulation, metabolic process, response to stimulus, cellular process, signaling, multicellular organismal process, reproduction, immune system process, developmental process, and localization. The most (24 out of 34) miRNAs participated in the cellular and metabolic processes, and the rest 12 miRNA families may be involved in other processes. For example, miR168 and miR2910 may have a role in plant reproduction, whereas miR160, miR2014 and miR2916 might be associated with signal transduction. Using gene ontology analysis, Mao et al. [23] reported that abscisic acid and salicylic acid stimulus might be regulated by miR159 and miR858 in cucumber. Furthermore, according to gene ontology analysis, 3 miRNAs (miR399, miR1122 and miR2014) were determined to be regulated in response to boron stress.
In conclusion, we identifed 32 known and 4 new barley miRNAs, as well as 934 target genes using recently developed degradome analysis. The majority of the identified miRNAs were significantly responsive to boron stress in barley. In particular, the signal transduction mechanism in leaves regulated by miR408 plays an important role in boron tolerance in barley consistent with previous reports [40,60]. Figure S1 The sequences, additional properties, and stem-loop secondary structure of pre-microRNAs of Hordeum vulgare (DOC)