Developmental Changes in Barley Microrna Expression Profiles Coupled with Mirna Target Analysis

MicroRNAs are 19-to 24-nt-long single-stranded RNAs that are crucial regulators of gene expression which control plant development and response to environmental cues. We have analyzed microtranscriptomes of five barley developmental stages. Generally, during the barley development, miR168-3p and miR1432-5p levels increase while the 5'U-miR156-5p level decreases (with exception for the 2-week-old barley). We have identified two miR156-5p izomiRs (called 5'U-miR156-5p [20 nt] and 5'UU-miR156-5p [21 nt]), which were expressed differently during barley development. The 5' U-miR156-5p level decreased in 3-week-, 6-week-, and 68-day-old barley, when compared to the 1-week-old plants. Meanwhile , the 5' UU-miR156-5p level increased significantly in the 68-day-old barley plants. Moreover, only the 5' U-miR156 isomiR recognizes and guides unique transcription factor mRNAs from the Squamosa Promoter Binding Protein-Like (SPL) family. We identified many non-canonical microRNAs with changed expression levels during the barley development. Here, we present the profiles of microRNA expression characteristics for particular barley developmental stages. These analyses are accompanied by the experimental degradome analysis of miRNA targets.


INTRODUCTION
MicroRNAs are small molecules that act on complementary target mRNAs, leading to their degradation via the endonucleolytic activity of Argonaute (AGO).In Arabidopsis, there are ten AGO proteins (Garcia-Ruiz et al., 2015).AGO proteins that contain the DDH motif can cleave target mRNAs (Baumberger & Baulcombe, 2005;Carbonell et al., 2012).A microRNA bound with AGO creates RISC (RNA-Induced Silencing Complex).Hsp90 facilitates RISC assembly (Iki et al., 2010).The role of microRNA is predominantly related to guiding the AGO1 protein to the proper mRNA target.Canonical cleavage of target mRNA occurs at the position located opposite the 10-11th nucleotide (counting from the guiding miRNA 5' end).In Arabidopsis, more than half of the target mRNAs are transcription factors (TFs) (Fahlgren et al., 2007).The primary transcript of a micro-RNA is produced by RNA Pol II followed by the excision of the pre-miR stem-loop precursor by the DCL1 in the nucleus.Then, double-stranded microRNA/ microRNA* duplex is excised from the pre-miR precursor, exported to the cytoplasm, and bound by the AGO protein (Mallory & Vaucheret, 2010;Schirle & MacRae, 2012).There are 325 precursors and 427 mature Arabidopsis microRNAs, and only 69 and 71 barley micro-RNA precursors and mature microRNAs (respectively) are deposited in the miRBase (MiRBase, release 21: June 2014; (Kozomara & Griffiths-Jones, 2011)).The nucleotide sequence, miRNA length, 5' end terminal nucleotide of microRNAs, and miR/miR* duplex structure dictate the sorting and future action of a particular microRNA (Zhang et al., 2014).In Arabidopsis, a predominant fraction of microRNAs (76%) contains molecules that are 21 nt long, while the remnant microRNAs represent smaller populations: 11% -22 nt long; 9% -20 nt long; 3% -24 nt long; 1% -23 nt long; and only 0.23% -19 nt long.Additionally, in the case of 56% of 21 nt Arabidopsis microRNAs, the first nucleotide from the 5' end is occupied by uridine, which suggests miRNA association with AGO1 and its involvement in mRNA cleavage (unpublished data, Zhang et al., 2014).
Barley development is mostly described by using phenotypical features like leaf number, tiller numbers, and kernel stage (Zadoks et al., 1974).In our previous work, we described pri-microRNA expression profiles in five barley developmental stages (Kruszka et al., 2013).Moreover, in the mirEX 2.0 database (an integrated environment for the expression profiling of plant microRNAs), we included 140 barley pri-microRNAs and also mature microRNA expression profiles for these 5 barley developmental stages (Zielezinski et al., 2015).Although the pri-miRNA and miRNA expression profiles are deposited in the mirEX 2.0 database, they were never analyzed and compared in detail in the different barley developmental stages.Here, we present a global analysis of the barley microtranscriptome profiles identifying the mi-croRNAs specific for the various barley developmental stages.Moreover, we also identified miRNA-targeted mRNAs.In addition, we present a deeper microRNA isomiR analysis in the case of miR156-5p.
RNA isolation.A total RNA enriched in small RNAs was isolated using the previously described approach (Kruszka et al., 2013).The quality of the isolated RNA was verified using an Agilent 2100 Bioanalyzer and Nano Plant RNA assay (Agilent).
Small RNA library preparation.Four independent small RNA libraries were designed for each barley developmental stage.Three replicas were performed using an Illumina TruSeq Small RNA Library Preparation Kit (Illumina).The quantification of libraries was carried out using a Quant-iT PicoGreen dsDNA Assay kit (Molecular Probes) reagent and Tecan Infinity 200 Pro Spectrometer (Tecan).NGS was performed using a HiScan-SQ machine (Illumina).The fourth replicated library for each developmental stage was performed using a different protocol.In brief, 1 μg of total RNA was ligated to a 3' adenylated adapter using a truncated T4 RNA ligase2 (New England Biolabs).After 1 h of incubation at 22°C, an RNA 5' adapter was added with T4 RNA ligase (Ambion), ATP and incubated at 37°C for 1 hour.cDNA was synthesized using Superscript II Reverse Transcriptase (Invitrogen).PCR was performed using a Phusion Hot Start High Fidelity DNA polymerase (New England Biolabs).PCR was run at 98°C for 30 sec, followed by 15 cycles at 98°C for 15 sec, 62°C for 30 sec, and 72°C for 15 sec, followed by a 5 min incubation.Products were separated electrophoretically, and the desired miRNA fraction in size was excised and eluted.The quality was tested using an Agilent Bioanalyzer DNA1000 chip and quantified using an Qubit fluorometer (Invitrogen).Each library was diluted to a target concentration for cluster generation and then loaded onto a lane of an Illumina GAiiX flow cell.The fourth library and its sequencing was performed at the BC Cancer Agency, Vancouver, BC, Canada.
Dataset analysis.Four replications were used for each particular barley developmental stage.Deep-sequencing data analysis was carried out using the approach described previously (Alaba et al., 2015).Then the normalized counts for each sRNA were averaged from replicates and the fold change and p-value were calculated between conditions using the Student's t-test.In the case of replications where no reads were detected (0 reads), value of 1 was added to the average in order to enable fold change calculation.For further analysis of the expression pattern, we considered microRNAs that have abundance described by either log 2 (Fold Change) >1 (increase higher than 2 fold) or log 2 (FC) <-1 (decrease higher than 2 fold) with the p value <0.05 (Student's t-test).Barley pri-microRNA and microRNA expression profiles deposited in the mirEX 2.0 database are available at the following web page: http://www.combio.pl/mirex(Zielezinski et al., 2015).
Target prediction and verification based on the degradome data.Possible target identification was performed using the psRNATarget online tool, (http:// plantgrn.noble.org/psRNATarget)with default settings (Dai & Zhao, 2011).When possible, the predicted targets were verified with our data obtained from the degradome sequencing.
Barley stem-loop pre-miR structure identification.For construction of the stem-loop pre-microRNA precursors, the following sequences were used: Sequences presented above were found by using the IPK BLAST Server (IPK Gatersleben -http://webblast.ipk-gatersleben.de/barley/)and the ViroBLAST online tool (Deng et al., 2007).The obtained sequences were folded using Folder version 1.11 (RNAfold algorithm) software available at http://www.ncrnalab.dk/#rnafolder/rnafolder.php.The longest sequences which could form the stem-loop structures were used for pre-microRNA construction.The pre-microRNA structures with the lowest ΔG energy value were chosen.

RESULTS
Using the Illumina deep-sequencing approach, we analyzed the microtranscriptomes of five barley developmental stages with regard to phenotypic traits.For each stage, four small independent RNA libraries were constructed and sequenced.Normalized reads from a particular developmental stage were compared, and fold changes were calculated.Analyses were performed in three ways: (i) we identified the differences in microRNA expression profiles in 2-week-, 3-week-, 6-week-, and 68-dayold barley plants, when compared to the 1-week-old barley plants (Suppl.Table 1 at www.actabp.pl); (ii) we compared microRNA expression profiles between particular barley stages to find unique expression profiles for the selected developmental stage (Suppl.Table 2 at www.actabp.pl); Based on the statistically significant differentially expressed microRNAs presented in Suppl.Table 1, we analyzed 14 up-regulated and 11 down-regulated micro-RNAs in further steps.The same sets of microRNAs were analyzed for expression profiles in the barley developmental stages compared to 1-week old barley plants (Fig. 1) and other barley developmental stages compared between each other (Fig. 2).
(iii) we analyzed potential target mRNAs for the selected microRNAs based on the performed degradome analysis.Degradome data were obtained for 68-day-old barley plants (Table 1).Characterization of the identified target mRNAs was carried out by comparisons to known nucleotide/protein sequences (Table 2).

Barley microRNAs showing different expression profiles during development
MiR168-3p and miR1432-5p are more frequent, and 5'U-miR156-5p (20 nt) is less frequent in all of the older Developmental microRNA expression profiles barley stages (except for 5'U-miR156-5p in 2-week old barley), when compared to the 1-week stage (Fig. 1A-H, Suppl.Table 1).MiR168-3p was earlier considered as a passen-ger miRNA for miR168-5p.MiR168-3p has a higher abundance level at all stages when compared to the 1-week-old barley plants, with the highest 2.8 fold change in 2-week-Table 1. Barley microRNA targets.The targets were identified based on the degradome data obtained from PARE-Seq of 68-day-old barley.Accession no, Target id -based on Ensembl Plants accession number; length, length of the analyzed cDNA; cleavage, cleavage_site (deg_ score|deg_rank) -cleavage position, degradome score, degradome ranking; * cleavage position between nucleotides with match to 10th and 11th nucleotide of microRNA counting from the 5' end, + non canonical cleavage position; score -aln_score -alignment score between microRNA and target sequence (the lower the value, the better the score); energy, hybridization energy (kcal/mol).Only targets with the best score numbers were included in the old plants (Fig. 1A-D).The degradome-approved mRNA target for this miRNA encodes the regulator of chromosome condensation (RCC1) protein (Table 2).

MicroRNA expression level as a specific marker for a given developmental stage
Comparison of microRNA expression in 2-versus 1-week-old barley plants.
Comparison of microRNA expression in 3-versus 1-week-old barley plants.
Comparison of microRNA expression in 6-versus 1-week-old barley plants.

DISCUSSION
Deep sequencing of small RNAs is a powerful tool to discriminate between microRNA sequences and their abundance (Sobkowiak et al., 2012;Lukasik et al., 2013;Barciszewska-Pacak et al., 2015).Analysis of the changes in small RNA level between the particular developmental stages of barley plants revealed some interesting data, of which two observations seem to be the most profound: (i) the 20-nt 5'U-miR156-5p level decreases during barley development, while in spike (68-day-old plants), upregulation of the 21-nt 5'UU-miR156-5p is observed.Both isomiRs target the same mRNAs, but the shorter one additionally guides the cleavage of two unique mRNAs; (ii) the level of miR390-5p was specifically elevated only in the 2-week and 6-week-old plants (more than 2 fold).Some interesting observations concern the microRNA and its microRNA* levels.Previously, we identified miR159b-3p as an up-regulated microRNA during barley development (Kruszka et al., 2013).Here, we found that, apart from the miR159b-3p increase, the highest fold change is observed for its microRNA* -miR159b-5p.MiR168-3p is up-regulated at all stages as compared to the 1-week-old barley, while the 21 nt miR168-5p is significantly up-regulated only in the 68-day-old barley plants.Our previous work had shown that the level of miR168-5p was almost unchanged and the level of miR168-3p was even decreased when compared to the 1-week barley stage (Kruszka et al., 2013).However, the previous analysis was based on northern blot hybridization, which is less sensitive and often does not discriminate between microRNAs belonging to the same family.Thus, sRNA deep sequencing enables a much-deeper analysis of microRNAs derived from the same family or isomiRs.
It is known that the miR156-5p expression profile is the opposite of the miR172 expression pattern.In rice, miR156 targets 11 SPL genes (Xie et al., 2006;Xie et al., 2012).SPL TFs regulate miR172 expression (Spanudakis & Jackson, 2014).It is possible that the same regulation occurs in barley plants.In further studies, it will be interesting to establish which miR156-5p isomiR affects the miR172 expression in barley.As previously mentioned, 5'U-miR156-5p 20 nt is down-regulated, but its 21-nt isomiR with additional U residue at the 5' end is up-regulated in the 68-day-old barley when compared to 1-week-old plants.One additional nucleotide at the 5'end of the miR156-5p allows for the recognition of different SPL mRNAs.In the degradome dataset, we found indeed that different SPL mRNAs are targeted by these two isomiRs of the miR156-5p.Comparison with known Arabidopsis SPL TFs showed that SPL mRNA that are targets for miR156-5p belong to different SPL classes.In agreement with these observations, we found that 20-nt 5'U-miR156-5p recognizes SPL mRNAs similar to Arabidopsis SPL9 and SPL15 TF mRNAs (Fig. 4B).Thus, microRNA isomiRs derived from a single pri-miRNA can specifically regulate various SPL gene expression.
It will be interesting to find an exact relationship between the up-and down-regulated miR156-5p isomiRs and their SPL mRNA targets.The deep-sequencing data of sRNAs isolated from plants from several developmental stages confirm earlier observations published in our previous paper (Kruszka et al., 2013): Northern blot hybridization showed the presence of two microR-NAs156-5p that differed in size.The northern signals in the 68-day-old barley reflect up-regulated 21-nt-long 5'UU-miR156-5p.In the mirEX2.0database, there are expression profiles of three barley miR156 family members: miR156b, miR156c, and miR156d.These three microRNA156a/b/c are downregulated in the 6-weekand 68-day-old barley plants.This reflects the decreased expression of 5'U-miR156-5p.MiR156 expression decreased in shoots or whole rice plants; however, careful inspection of selected organs revealed that this microR-NA level is increased in older leaves (Xie et al., 2012).The authors did not discriminate between the isomiRs.
Here, we have shown that, in barley plants at the spike-formation stage (68-day-old plants), there are two miR156-5p izomiRs, which are differentially expressed.Compared to the 6-week-old plants, we observed a decreased level of miR390-5p in the 68-day-old barley plants.This microRNA regulates the production of Arabidopsis ta-siRNAs that, in turn, downregulate ARF2, ARF3, and ARF4 TF levels.It was shown that the decreased level of ARF2, ARF3, and ARF4 TFs promotes lateral root growth (Marin et al., 2010).It would be interesting to test the level of these transcription factors in the 68-day-old barley plants and during lateral root growth.One of the elevated microRNAs in older barley plants is the miR1432-5p that targets the mRNA-encoding protein containing 2-oxoglutarate (2OG) and Fe(II)dependent oxygenase domains.We blasted the barley MLOC_61497.1 (target of the miR1432-5p) to the NCBI protein database.The highest similarity to this barley sequence was found in the case of Naringenin, 2-oxoglutarate 3-dioxygenase from Triticum urartu (EMS57783.1,E value=0.0,identities=95%, 323/339).This enzyme catalyzes the 3-beta-hydroxylation of 2S-flavanones to 2R, 3R-dihydroflavonols, which are intermediates in the biosynthesis of flavonols, anthocyanidins, catechins, and proanthocyanidins.Altogether, we observed that during barley development, microRNA expression profiles change qualitatively and quantitatively.These changes influence the mRNA target levels and plant development.