Small RNA and Degradome Analyses Uncovering the Extensive Effects of miRNAs and Targets in Early Developing Grains of Common Wheat

Qiaoyan Chen Henan Institute of Science and Technology Lina Xu Henan Institute of Science and Technology Yuanyuan Guan Henan Institute of Science and Technology Zeeshan Ali Buttar Henan Institute of Science and Technology Gan Li Henan Institute of Science and Technology Guang Zhang Henan Institute of Science and Technology Chunxiao Li Henan Institute of Science and Technology Jinyong Huang Zhengzhou University Mingjiu Liu Henan Institute of Science and Technology Wenhui Wei (  whwei88@hotmail.com ) Henan Institute of Science and Technology

the detected miRNAs were up-regulated from 5 DAP to 20 DAP, and miR164 and miR160 increased abundantly during grain development [27]. Under the condition of over-expression, miR397 gave rise to an increase of grain yield by enlarging grain size and promoting panicle branching in rice [28].
MiRNAs in wheat were identi ed from seeding [29], ag leaf and developing seed [27,30], callus [31], germ extract, leaves, roots, spikes and multiple tissues [32]. The identi cation methods used mainly include high-throughput sequencing, in silico prediction and wet-lab validation [33]. To date, it has been found that wheat miRNAs take part in the regulations of grain development [34], stripe rust resistance [30], powdery mildew resistance and heat stress response [35], as well as cold stress response [16], etc.
Degradome sequencing is an effective means in studies of degraded mRNAs. This technology could identify the target mRNAs of miRNAs on a large scale when it is combined with high-throughput sequencing and bioinformatics analysis to search for miRNAguided cleavage sites in mRNAs [36]. In the same way, combined with small RNA sequencing, some novel miRNAs with low abundance and their targets could be identi ed successively [37,38]. In wheat, some miRNAs and their targets have been effectively studied with small RNA and degradome sequencing [39,40].
The early stage of grain development is an important period for determining grain yield and quality characteristics for wheat harvest. However, the regulation of grain development is unclear during the early stage. In the present study, wheat grains at 7 DAP and 14 DAP, which represent two transitional early developmental phases of wheat grains, were isolated. High-throughput small RNA and degradome sequencing were performed to explore the miRNAs and their target genes possibly participating in the regulation of early grain development.

Results
Phenotypic analysisof developing grain The development of wheat grain was evaluated by monitoring the pattern of increased grain weight and volume. Here, the growth of grain weight and volume was relatively slow at the early stage (before 11 DAP), increased sharply across 11 to 14 DAP and continued to rise until about 26 DAP (Fig. 1A). The appearances of developing grains at some representative time points are shown in Fig. 1B. Based on the pattern, we focused mainly on the transitional phases at 7 DAP and 14 DAP, the key periods for early grain development. To reveal the relevance between these changes and miRNAs abundance during early grain development, we then used the small RNA and degradome libraries from 7 DAP and 14 DAP grains to compare the expression of miRNAs and analyze the possible functions of their targets.
Deep-sequencing of sRNAs in early developing grains To uncover the roles of miRNAs during wheat grains development at early stage after pollination, the grains of cultivar "Bainong 4199" at 7 DAP and 14 DAP were used to construct two small RNA libraries and then they were sequenced. Firstly, 22,900,708 raw reads from 7DAP and 20,279,914 raw reads from 14 DAP were produced, respectively. Then, the low-quality reads were ltered, adapter sequences and sequences of less than 18 nt and more than 30 nt were removed. Finally, 14,020,684 and 10,379,004 clean sRNAs were obtained from 7 DAP and 14 DAP libraries, respectively ( Table 1). The length of clean sRNAs ranged from 18 to 30 nt in the two libraries. Majority of the clean sRNAs were 21-24 nt in length. It is found that the most common length is 21 nt in the 7 DAP library and 24 nt in the 14 DAP library (Fig. 2).

Identi cation of known miRNAs
To identify known miRNAs, we compared the unannotated sRNA sequences that matched perfectly to reference genome against miRBase 22.0 (http://www.mirbase.org) based on perfect match criterion. A total of 89 known miRNAs from wheat were detected in miRBase/Triticum aestivum in the two sRNA libraries, of which 46 were expressed in 7 DAP library and 87 were expressed in 14 DAP library (Additional le 2). Of the known miRNAs from other plant species, their secondary structure and miRNAs* were predicted based on RNAfold (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi), and these miRNAs were predicted based on MIREAP (http://sourceforge.net/projects/mireap/) by analyzing DCL1 cleavage sites and the minimum free energy. A total of 16 known miRNAs from other plant species were identi ed in the present two sRNA libraries, with 13 expressed in 7 DAP library and 16 expressed in 14 DAP library (Table 3; Additional le 2). The secondary structures of these miRNAs are shown in Additional le 3.

Predicted novel miRNAs
To predict novel miRNAs, the unannotated sRNAs were analyzed on the basis of pre-miRNA sequences that can form a canonical stem-loop hairpin structure [41]. A total of 79 novel miRNAs were predicted, of which 32 were expressed in 7 DAP library and 78 were expressed in 14 DAP library ( Table 4). The minimum free energy of pre-miRNAs ranged from -187.00 kcal mol −1 to -37. More novel miRNAs were expressed in 14 DAP than in 7 DAP, and 31 novel miRNAs were commonly expressed in two libraries, with 1 special in 7 DAP library and 47 specials in 14 DAP library. Majority of novel miRNAs presented a relatively low expression level. The most abundant novel miRNAs were novel-m0064_5p, novel-m0492_5p and novel-m0661_5p (Additional le 4).

Degradome sequencing and data summary
After RNA was extracted and poly (A+) RNA was puri ed from the developing wheat grains at 7 DAP and 14 DAP, a degradome library was constructed and subsequently subjected to degradome sequencing. A total of 10,620,205 clean reads and 4,612,965 unique reads were obtained (Table 5). A total of 2,776,485 (60.19%) unique reads were matched to the reference genome. Using BLASTN search against GenBank and Rfam databases, the structural RNAs (rRNAs, tRNAs, snRNAs, and snoRNAs) were removed, and the remaining reads were used to detect candidate targets of miRNAs.

Differentially expressed miRNAs between 7 DAP and 14 DAP grains
To identify miRNAs associated with early wheat grain development, the differentially expressed miRNAs were analyzed based on the statistical method reported by Audic and Claverie [42]. A total of 19 known and 20 novel miRNAs were differentially expressed (P < 0.05) between the 7 DAP and 14 DAP grains, and 18 known miRNAs and 13 novel miRNAs were up-regulated in 14 DAP grains (Additional le 7).
To verify the miRNA expression levels and the deep-sequencing results, eight known and two novel miRNAs were randomly selected for quantitative real-time polymerase chain reaction (qRT-PCR). The results showed that the expression patterns of these miRNAs were consistent with those from deep-sequencing, which indicated that our sRNA sequencing was credible (Fig. 6).

Functions of miRNA targets
To understand the functions of miRNAs in the regulatory network of early wheat grain development, miRNA target genes were analyzed across small RNA and degradome libraries. A total of 266 targets for 40 known wheat miRNAs, 152 targets for 13 other known plant miRNAs and 258 targets for 25 novel miRNAs were predicted (Additional le 8).
Functional annotations for the target genes of differentially expressed miRNA were performed by BLAST analysis, and it was found that the targets included those genes encoding mitogen-activated protein kinase, transcription factor GAMYB, WRKY transcription factor 33, F-box/kelch-repeat protein, transcription factor PCF6, NAC domain-containing protein, ethylene-responsive transcription factor, SPX domain-containing protein, vegetative cell wall protein gp1, extensin precursor, cytokinin dehydrogenase 5, calciumdependent protein kinase and squamosa promoter-binding-like protein (SPL), etc. (Additional le 9). GO analysis of the differentially expressed miRNAs between 7 DAP and 14 DAP grains was applied to classify their target genes according to their cellular component, molecular function, and biological processes. A total of 10 molecular functions were identi ed, of which binding, catalytic activity, and nucleic acid transcription factor activity were the three most frequent functions. For biological processes, 20 categories were identi ed, of which cellular processes, metabolic processes and single organism processes were the three most frequent ones. For the cellular component, 11 categories were identi ed, of which cell part, cell and organelle were the three most frequent processes (Additional le 9, Fig. 7).
Furthermore, gene-speci c qRT-PCR was performed to detected the expression levels of potential targets. The 8 target genes of 8 miRNAs ( ve known and three novel) were selected randomly. The expression analysis indicated that the target genes were negatively correlated with the expression of their corresponding miRNAs (Fig. 8). These results suggest that miRNAs are involved in the developmental processes of wheat grains by regulating their target genes expression.

Discussion
It is known that miRNAs are involved in diverse developmental processes. Previous studies have showed that miRNAs were associated with development and stress response in wheat. Here, we investigated the miRNAs and their targets to nd the association between miRNAs and early grain development of wheat, which is one of the most important crops for human foods.
In this study, 105 known and 79 novel miRNAs were identi ed from two libraries of grains at 7 DAP and 14DAP. Expression analysis showed that 19 known and 20 novel miRNAs were differentially expressed, and 18 known and 13 novel miRNAs were up-regulated at 14 DAP. During rice grain development, most miRNAs showed increased expression from 5 to 10 DAP [26]. A few researches have been reported in wheat. Most of the differentially expressed miRNAs were up-regulated with grain and leaf development; only little parts were down-regulated [27,40]. These results suggest that most miRNAs were highly expressed during organ development in crops.
Of the differentially expressed known miRNAs, miR156 was up-regulated at 14 DAP compared to 7 DAP in the present study. Pioneering research showed that miR156 directly repressed the expression of SPL transcription factor genes that play an important role in plant growth and development [43]. In rice, OsSPL16 controls grain size, shape and quality [44], and OsSPL13 positively regulates the cell size of grain hull and then controls grain size [45]. TaSPL16, a target of miR156 in wheat, is highly expressed in developing young panicles, and lowly expressed in developing grains [46], because the expression level of miR156 was up-regulated with the development of wheat grain. Therefore, miR156 may well play an important role in regulating wheat grain development.
It has been reported that miR164 targets NAC (NAM, ATAF, and CUC) transcription factor genes that play major roles in the proper formation and separation of plant organs [47], auxin signaling and defense [48]. Auxin is crucial during the grain/seed development, including pattern formation, cell division and cell expansion [49]. In wheat grain, NAC genes regulating senescence have the functions to improve protein, zinc and iron contents [50]. In this study, miR164 is up-regulated from 7 DAP to 14 DAP grains. It is also shown that miR164 increased in abundance from 5 DAP to 20 DAP [27]. Furthermore, NAC transcription factors were predicted to be involved in regulating the timing of organ formation (GO:0048504), negative regulation of cell division, regulation of cell size, response to temperature stimulus and light stimulus (Additional le 9). It may be that miR164 takes part in the regulation of early wheat grain development.
MiR319 (bdi-miR319b-3p, tae-miR319) targeted TCP (TEOSINTE BRANCHED 1, CYCLOIDEA and PCF) transcription factor genes in Arabidopsis thaliana. And the research showed that TCP 3 regulated the activities of miR164 during the differentiation of leaves in Arabidopsis [51]. In our study, the miR319 was down-regulated at 14 DAP, perhaps with the target up-regulated. The target in turn regulated miR164 expression. And we predict that the targets of miR319 are involved in negative regulation of cell proliferation and differentiation, positive regulation of ovule development (GO:0048481), mucilage biosynthesis involved in seed coat development, response to auxin and abscisic acid, response to gibberellin, jasmonic acid-mediated signaling pathway and ethylene-activated signaling pathway (Additional le 9). A putative miRNA regulatory network in wheat grains showed that miR156, miR164, miR319 and miR396 play a role in cell proliferation [40]. In the present study, however, miR396 showed no difference during the early development of wheat grains.
MiR167 targets ARF (AUXIN RESPONSE FACTOR) transcription factor genes, which play critical roles in regulating plant growth and development [52], ARF8 is a negative regulator of fruit initiation in Arabidopsis [53]. In this study, miR167 is up-regulated from 7 DAP to 14 DAP. In previous research, it was up-regulated from 5 DAP to 15 DAP [1], 7DAP to 14 DAP and down-regulated from 14 DAP to 28 DAP in wheat grain [40]. These results showed that MiR167 is important in wheat grain development. Research showed that ARF members are targeted by miR164, whereas many TFs such as TCP members are the regulators of miR167 [54], indicating mutual regulation between TFs and miRNAs.
Known miRNA such as tae-miR7757-5p targets the genes encoding WRKY DNA-binding domain and NB-ARC domain-containing proteins involved in jasmonic acid-mediated signaling pathway, regulation of plant-type hypersensitive response, response to hydrogen peroxide, and regulation of multi-organism process. MiR9662 (tae-miR9662a-3p, tae-miR9662b-3p) targets the genes related to vegetative cell wall protein, protein kinase activity (GO:0004672), protein phosphorylation (GO:0006468), plant-type cell wall organization (GO:0009664), multicellular organismal process (GO:0032501), and regulation of cellular process (GO:0050794) (Additional le 9). This shows the extensive effects of miRNAs in early wheat grain development and might participate in regulating wheat grain development and metabolism. Previous reports indicated that tae-miR1127b was important for early grain development and played a crucial role in the uptake of amino acids into the endosperm [27].
MiR156, miR164, miR319 and miR167 were differentially expressed in developing grains, as found in previous studies. However, some miRNAs, such as miR159 and miR396 targeting the growth-regulating factors (GRFs), were expressed at the same level in 7 DAP and 14 DAP grains. In previous research, miR396 was under-expressed and miR159 was over-expressed from 7 DAP to 14 DAP. They were predicted to target the mRNA of a GAMYB-like transcription factor gene involved in gibberellin signaling and plant anther development [40]. We guess that some miRNAs presented adaptive expression levels under different genetic backgrounds and environmental conditions. The novel-miR0018 targets WRKY transcription factor genes. Previous research postulated that WRKY23 participated in the regulation of plant stem cell speci cation via auxin-dependent or auxin-independent signaling pathway [55]. WRKY23 could regulate auxin distribution patterns through controlling avonol biosynthesis during Arabidopsis root development [56]. Perhaps novel-miR0018 plays important roles during the development of early wheat grain.
The novel-miR0036 was signi cantly down-regulated at 14 DAP compared to 7 DAP (Additional le 7), and it was predicted to target DHHC palmitoyl transferase gene. DHHC proteins regulated cell function and in uenced cell physiology and pathophysiology [59], indicating that novel-miR0036 played an important role during wheat grain development. And other novel miRNAs target the genes that are involved in lipid transport and metabolism (novel miR0005), signal transduction mechanisms (novel miR0079), posttranslational modi cation (novel miR0078), chromatin structure and carbohydrate transport (novel miR0034, novel miR0070), ribosomal structure and biogenesis (novel miR0075), replication, recombination and repair (novel miR0031) (Additional le 9). Perhaps these novel miRNAs widely participate in the regulations of early wheat grain development and metabolism.

Conclusion
This study suggests that quite a few known and novel miRNAs and their targets play extensive roles during the early grain development of common wheat. Understanding of miRNA-mediated regulatory network involved in wheat grain development will help us to elucidate the molecular mechanisms underlying wheat grain development and carry out ingenious molecular improvements in wheat breeding.

Plant materials and sample preparation
The experiments were performed using the 7DAP and 14DAP grains of common wheat cultivar "Bainong 4199" (Triticum aestivum L.) bred and friendly provided by the Centre for Wheat Breeding, Henan Institute of Science and Technology. "Bainong 4199" is a high yield cultivar with plump grain and is widely planted in China. Seeds of "Bainong 4199" were rstly soaked in the tap water for 24 h, then disposed in 4°C dark chamber for one month. After vernalization, they were planted in the greenhouse maintaining 75% relative humidity, 26/20 °C day/night temperature, 12-h light/dark photoperiod, and 10,000 lux light intensity. Day of pollination was recorded when half of the plants reached the owering stage. For volume and fresh weight measurements, immature grains were collected for every 3-5 days, starting from 5 DAP to 34 DAP. And the grains were collected from the middle four rows of spikes at 7 and 14 DAP for deep sequencing and miRNAs analysis. Three biological replicates were made, and each consists of 100 grains. All samples were snap-frozen in liquid nitrogen and stored at −80°C for subsequent experiments.

Small RNA and Degradome Sequencing
Total RNAs were isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions.
Quali ed RNA samples must meet the conditions of OD260/280 = 1.8-2.2 and RNA integrity number > 8.0. RNA samples in three biological replicates were equally mixed to get respective RNA pool of 7 DAP and 14 DAP grains. These two RNA pools were used to construct small libraries by TruSeq Small RNA Library Prep Kit (Illumina, USA), according to the manufacturer's protocols. In brief, 700-μg RNA sample of each RNA pool was rstly fractionated on a 15% denaturing polyacrylamide gel, then the sRNAs with 18-30 nt were recovered. The 5' and 3' RNA adapters were ligated to the 5' and 3' ends of these sRNAs using T4 RNA ligase (Takara, Dalian, China). Puri ed ligation products were converted to cDNAs, which were used to construct the cDNA tag libraries by RT-PCR ampli cation. The size, purity and concentration of cDNA tag libraries were detected using an Agilent2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The cDNA tag libraries were sequenced using a HiSeq™ 2000 Sequencing System (Illumina, San Diego, CA, USA), according to the manufacturer's instructions.
The samples for miRNA sequencing were used to construct the degradome libraries. Degradome libraries were constructed by ligating polyA-enriched RNAs to the custom RNA adapter containing a 3' Mme I site. This is followed by reverse transcription, second-strand synthesis, Mme I digestion, ligation of 3' dsDNA adapter, gel-puri cation and PCR ampli cation. Ampli ed degradome tag libraries were then sequenced using a Solexa/Illumina genome analyzer [60].

MiRNAs identi cation and expression analysis
The reads from cDNA tag libraries sequencing were processed by Phred and Crossmatch (http://www.phrap.org/phredphrapconsed.html) [61]. Clean reads were obtained after low-quality reads and adapter sequences were removed. The 18-30 nt clean reads (tags) were matched to the sequences of Rfam [62], GenBank, and RepBase [63] databases to distinguish rRNA, scRNA, snoRNA, snRNA and tRNA from clean reads. After those reads having more than 90% sequence similarity to above RNAs were removed, the remaining sequences were matched to wheat expressed sequence tag (EST) database (http://www.ncbi.nlm.nih.gov/nucest/?term=wheat) or wheat genome sequences (http://mips.helmholtzmuenchen.de/plant/wheat/uk454survey/index.jsp) [64]. Ultimately, the non-coding tags were used to identify the candidate miRNAs.
Known miRNAs and novel miRNAs were identi ed from above non-coding tags referring to the methods reported by Chu et al. [31], with some modi cations. In brief, the non-coding tags were rstly matched to the sequences of miRBase 22.0, the perfectly matched tags were the known miRNAs, and the others were used to predict novel miRNAs. The predicted miRNAs should have the potential forming a hairpin secondary structure when analyzed using RNAfold (http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi) [65].
For eliminating the disturbing of gene length and sequencing deep to the count of miRNA actual expression level, the frequency of each read count was rstly normalized according to the formula TPMi = (Ni/Li) × 1,000,000/sum (Ni/Li +……..+ Nm/Lm); Ni is the number of reads mapped to gene i, and Li is the length sum of the exons of gene i. The fold-change of miRNA expression between the 14 DAP and 7DAP libraries was calculated as log 2 (14DAP/7DAP). If the P-value was ≤ 0.01 and the normalized sequence count had a fold change of > 2 or < 0.5, the given miRNA was recognized to be differentially expressed.

MiRNA target annotation
Prediction and annotation of miRNA targets were performed as the methods reported by Chu et al. [31]. Identi ed miRNAs were aligned against wheat sequence databases. MiRNA targets were predicted according to the TargetFinder (https://github.com/carringtonlab/TargetFinder?). BLASTN hits with less than four mismatches were chosen as candidate targets; for the non-target predicted miRNA, the psRNATarget software (http://plantgrn.noble.org/psRNATarget/) (version 12) was used to predict the targets in wheat transcripts with prediction score cutoff value = 3.0, length for complementarity scoring = 20, and target accessibility = 25.
For the degradome sequencing data, 20-21 nt sequences of high quality were collected for subsequent analysis. The unique reads that perfectly matched wheat expressed sequence tag (EST) database from NCBI or the contig sequences from WGS assembly were retained. Approximately 15 nt upstream and downstream of 5' wheat EST sequences, mapped by degradome reads, were extracted to generate 31 nt target signatures as 't-signatures' [66]. The CleaveL and pipeline were used to identify miRNA targets [60]. Alignments with scores up to four where G:U pairs scored 0.5 and no mismatches were found at the site between the 10th and 11th nucleotides of the corresponding miRNAs were considered potential targets. Consistent mRNAs obtained from both methods were chosen as miRNA targets. To understand their functions, the putative target genes of the miRNAs were BLASTX and subjected to GO analysis.

MiRNA and Target validation by real-time PCR
QRT-PCR was performed to determine the validities of RNA Seq and target analysis. A total of 10 miRNAs and 8 target genes were selected randomly. RNA was extracted from three independent biological samples of 7DAP and 14DAP grains, individually, and used for transcription (RT) reactions. The One Step Primer Script miRNA cDNA Synthesis Kit (Takara) and PrimeScript® RT reagent Kit with gDNA Eraser (Perfect Real Time; Takara) were used for the RT reactions of miRNAs and targets following the manufacturer's instructions, respectively. RT-PCR was performed using a Bio-Rad IQ5 Real-Time PCR Detection System (BIO-RAD, Hercules, CA, USA) with SYBR® Premix Ex Taq II™ (Takara). The total volume of each reaction is 25 μL [2.0 μL of diluted product, 2.0 μL of primers, 12.5 μL of SYBR® Premix Ex Taq™ (Perfect Real Time; Takara), and 8.5 μL of nuclease-free water]. All reactions were performed rstly at 95 °C for 30 s, then followed by 40 cycles ( 95 °C for 5 s, 61 °C for 30 s, and 72 °C for 30 s). All reactions were performed three times.
Wheat U6 [GenBank: X63066] snRNA and the actin gene (GenBank: AB181991) were used as the endogenous control for miRNAs and target genes. The relative abundance of each miRNA was calculated by a comparative C T method (ΔΔC T ) using the formula 2 -ΔΔC T [67]. The miRNAs and target samples in the 0 DAP library with the C T value were selected as the calibrator, and the expression level was set as 1.0. The relative expression levels of the same miRNAs and targets were normalized through comparison. All primers used in this study are presented in Additional le 10.

Consent to publish
Not applicable.

Availability of data and materials
The datasets used and/or analyzed during the current study are available from the Additional les or corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests.