Dry matter accumulation and appearance of maturing grain

To evaluate the developmental process for wheat grains, we monitored the pattern by which dry matter was accumulated in those tissues. Here, seed weights increased sharply between 12 and 18 days after pollination (DAP) and continued to rise until 30 DAP (Figure 1A). Based on the pattern we observed, we continued our focus mainly on activities between 15 and 25 DAP, the key period for grain-filling, as well as on grains sampled at both 5 DAP and 30 DAP so that we could compare the expression profiles of miRNAs near the beginning and end of that stage of formation. The size and appearance of developing caryopses at those time points are shown in Figure 1B.

Deep-sequencing of sRNAs in developing grains

To investigate the enrichment of miRNAs, we generated four libraries from developing wheat grains sampled on four dates. After sequencing via the Illumina Genome Analyzer, we removed low-quality reads and corrupted adapter sequences (reads <18 nt or >30 nt long). In all, we obtained 13,525,513 reads (3,492,987 unique) representing 5 DAP, 14,460,398 reads (6,015,357 unique) for 15 DAP, 15,310,251 reads (9,243,757 unique) for 25 DAP, and 13,282,907 reads (6,200,456 unique) for 30 DAP. Our assessment of the size distributions demonstrated that approximately 98% of the detected sRNAs were 18 to 25 nt long, and that 24-nt sRNAs were prevailing in all stages whereas those that were 21 nt long were less abundant (Figure 2A). Further analysis revealed that the distribution of redundant sRNAs in various size classes was quite similar among the libraries from 5, 15, and 25 DAP grains. We found it interesting that the abundance of 24-nt redundant sRNAs was significantly decreased in 30 DAP grains while the proportion of those that were 18 nt to 22 nt were obviously increased (Figure 2A). The distribution patterns for unique sRNAs were similar among all four libraries (Figure 2B). The patterns of distribution between redundant and unique sRNAs indicated that the abundance of sRNAs fluctuated during the period of seed development. Therefore, we evaluated their relative abundance from 18 to 25 nt based on the ratio of redundant/unique and found that expression was generally altered over time (Figure 2C). Transcript levels for most sRNA sizes were obviously decreased in the 25 DAP samples. This was especially true for 24-nt sRNAs, which were gradually repressed. The canonical heterochromatic siRNAs are 24 nt long and these results suggested that siRNA-mediated gene regulation might be involved in the control of grain formation.

For annotation, we used the sRNA datasets to query the non-coding RNAs deposited in the NCBI GenBank (http://www.ncbi.nlm.nih.gov), Rfam database [32], and miRBase (http://www.mirbase.org/). They were then classified into seven categories: miRNA, rRNA, siRNA, snRNA, snoRNA, tRNA, and those detected but without annotation (Additional file 1). The number of unique miRNAs increased gradually as the process of grain-filling continued (Figure 3A), implying that miRNA-mediated gene silencing is involved in developmental regulation. Moreover, the total read of miRNAs was relatively lower at 25 DAP (Figure 3B), leading to a lower ratio (redundant/unique) of miRNAs at that time point (Figure 3C). Small RNAs were more abundant in the 25 DAP library than in any other (Figure 3D), suggesting that miRNA expression was generally down-regulated at that later stage of maturation.

Conserved miRNAs differentially expressed in developing grains

To identify the conserved miRNAs in developing grains, we aligned the sRNA sequences with known mature miRNAs from plants in the miRBase. In all, 605 miRNAs representing 540 families were identified (Additional file 2). In general, we found that 424 conserved miRNAs (382 families) were differentially expressed (p-values <0.05) over time (Additional file 3). Their scatter plots illustrated the similarities during the stages of grain formation, as well as the broad relationships among miRNA expression profiles from 5 DAP through 30 DAP (Figure 4). However, some differences in profiles were apparent from this examination. When compared with values at 5 DAP, 231, 266, and 309 miRNAs were significantly up- or down-regulated (p <0.05) at 15 DAP, 25 DAP, and 30 DAP, respectively. It is likely that the number of up-regulated miRNAs was gradually increased (114, 147, and 211 at 15 DAP, 25 DAP, and 30 DAP, respectively) over time, demonstrating that mechanisms for miRNA-mediated repression are involved in this developmental regulation.

Conserved miRNAs associated with grain-filling

Among those 424 differentially expressed miRNAs, expression was lower for 284 miRNAs, i.e., those with fewer than 20 TPMs (transcripts per million). For our investigation, we selected 140 miRNAs (125 families) highly expressed (at least has a TPM ≥20 in the four stages detected) in developing grains (Additional file 3). To find common expression patterns, we performed hierarchical clustering based on fold-changes in expression from a base line of 5 DAP. In all, 86 miRNAs could be sorted into eight clusters (Figure 5, Additional file 4) that we considered to be potentially involved in the regulation of grain-filling.

The first four clusters (I-IV) of conserved miRNAs generally included those that were up-regulated in the process of grain-filling. Clusters I and II comprised 16 miRNAs with expression that was very low at 5 DAP, but which greatly increased at 15 DAP, and remained high for most miRNAs through 30 DAP. Considering that 15 DAP and 25 DAP are crucial stages during which a complex gene regulatory network is involved in endosperm development, we concluded that those significantly induced miRNA candidates in Clusters I and II might have important regulatory roles. Expression of the miRNAs in Cluster III was gradually increased over time. For example, miR167a, miR156a, and miR156c, were quite abundant, with expression levels at 5 DAP of 1825.73, 845.74, and 212.12 TPM, respectively (Additional file 4). Although they were up-regulated by only a few fold at 15 DAP, they could have a greater regulatory effect because of their more abundant basis. Cluster IV had 11 miRNAs that were only slightly up-regulated. Of these, eight were more highly expressed, with TPM values of >100 at 15 DAP and 25 DAP. For example, the TPMs for miR168a were 7509.36 at 5 DAP and 21776.72 at 15 DAP.

Four clusters of miRNAs showed repression, including six members that were rapidly down-regulated in Cluster V. For example, although miR473 was highly expressed at 5 DAP (TPM = 150.97), its level decreased to 0.13 TPM at 25 DAP and was undetectable in 30 DAP grains (Additional file 4). Clusters VI and VII contained 3 and 9 miRNAs, respectively, all of which were significantly repressed at 15 DAP and/or 25 DAP. Finally, Cluster VIII had 20 members, including miR2628 and miR1146, were gradually repressed over time.

Novel miRNAs identified in wheat grains

Novel miRNAs were predicted according to the characteristic hairpin structures of their precursors, which distinguishes them from other endogenous sRNAs [33, 34]. A total of 268 putative novel miRNAs (182 families) were identified (Additional file 5, Additional file 6). Among them, 24 showed at least a relatively higher expression (TPM ≥10) in developing grains (Table 1). We were interested to learn that Ta-miR034-3p and Ta-miR034-5p are homologous to complementary sequences in the same stem-loop (Figure 6). Moreover, Ta-miR007-5p and Ta-miR007-3p were also located in the same stem-loop generated from a wheat EST (GB#: CJ832040).

Further analysis revealed that these 24 highly expressed novel miRNAs could be sorted into Clusters I through VI (Figure 7, Table 1). Cluster I contained four miRNAs for which expression rose sharply at 15 DAP and was maintained at a high level throughout the process of grain-filling. For example, Ta-miR021-1-5p was strongly up-regulated over time, peaking at 30 DAP (TPM = 7480.74). Expression of the five miRNAs within Cluster II was obviously increased at 15, 25, and 30 DAP when compared with levels detected at 5 DAP.

In Cluster III, Ta-miR042-3p, Ta-miR107-2-3p, and Ta-miR106-5p were highly expressed at 15 DAP, but not at any other developmental stages. By contrast, Ta-miR051-3p and especially Ta-miR154-5p, both in Cluster IV, were significantly up-regulated at 25 and 30 DAP. In Cluster V, four novel miRNAs were down-regulated over time. Ta-miR068-5p and Ta-miR057-1-3p were gradually repressed in developing grains whereas Ta-miR007-5p and Ta-miR007-3p were highly expressed at 5 and 15 DAP, but undetectable at 25 and 30 DAP. This contrast in up- and down-regulation among these 18 novel miRNAs from Clusters I through V demonstrated their important regulatory roles.

Finally, the novel miRNAs included in Cluster VI were either expressed only at 30 DAP (near the completion of the filling stage) or else their transcript abundance fluctuated over time. Therefore, their activity did not seem to be associated with this type of developmental regulation. The exception to this was Ta-miR034-5p, also in Cluster VI, which showed strong expression at all sampled time points.

Validation of miRNAs in developing grains

We conducted quantitative real time PCR (qRT-PCR) to validate the expression patterns of miRNAs identified via high-throughput sequencing. Based on their patterns of development-regulated expression, as determined by deep-sequencing, we selected eight miRNAs for examination. Seven of them -- miR167a, miR397, miR156a, Ta-miR021-1-5p, Ta-miR004-1-5p, Ta-miR044-1-3p, and miR827a -- were predicted to be induced while miR1852 was expected to be repressed. In fact, miR167a, miR397, miR156a, Ta-miR004-1-5p, Ta-miR044-1-3p,and miR827a were significantly up-regulated by 15 DAP and 25 DAP, with expression then peaking at 30 DAP (Figure 8). Their expression profiles were quite similar to those determined by high-throughput sequencing (Figures 5, 7). Although Ta-miR021-1-5p was up- regulated over time while miR1852 was down-regulated, the degree to which their expression was altered was not as dramatic as had been demonstrated with high-throughput sequencing (Additional file 4, Table 1). Those responses may have been affected by the relative abundance of the other sRNAs. Nevertheless, our qRT-PCR results were generally consistent with the data obtained from our high-throughput sequencing, thereby indicating that it is possible to create a set of grain filling-associated miRNAs through deep-sequencing of wheat.

Targets of grain filling-associated miRNAs

The potential targets of 86 grain filling-associated conserved miRNAs (Additional file 4), as well as 24 highly expressed novel miRNAs (Table 1), were computationally predicted using the Triticum aestivum (wheat) unigene library, DFCI Gene Index (TAGI; version 12), and the psRNATarget program (http://plantgrn.noble.org/psRNATarget/) [35]. All of the targets predicted for those conserved and novel miRNAs are shown in Additional file 7 and Additional file 8. Functional annotations of those target genes were performed by BLAST analysis with NCBI and it is found that these grain filling-associated miRNAs potentially target to multiple wheat genes, which including transcription factors, proteins involved in the membrane transporting, ubiquitin-mediated proteolysis, carbohydrate metabolism, responding to stress, signal transduction and phytohormone signaling.

To validate the expression patterns of potential targets for grain filling-associated miRNAs, we selected six target genes (TC383723, TC384445, TC370322, TC400547, DR734203 and TC370450) of candidate miRNAs, which were significantly up-regulated with the development of grains, and gene specific qRT-PCR was performed. The expression analysis demonstrated that all detected target genes were strongly down-regulated during grain-filling (Figure 9), which were negatively correlated to the expression profiles of their corresponding miRNAs (Additional file 4). These results suggested that miRNAs play a crucial regulatory role during the developmental processes of grain-filling in wheat.

Plant materials

Seeds of wheat (Triticum aestivum L. cv. Yumai 18) imbibed water for 24 h and were then transferred to a cold chamber (4°C in the dark). After 4 weeks of vernalization, they were planted in pots and cultured in a naturally lit phytotron glasshouse (25°C/20°C day/night). Developing grains were harvested from the middle four rows of a head at 3 to 33 DAP. Each experiment comprised 100 grains, with three replicates. Samples were immediately frozen in liquid nitrogen and stored at −80°C.

Total RNA isolation

Total RNA was extracted 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 by ethidium bromide, was run to determine the preliminary integrity of the RNA. All RNA samples were quantified and examined with an ND 1000 spectrophotometer (Nanodrop) for contamination with either protein (A260 nm/A280 nm ratios) or reagent (A260 nm/A230 nm ratios). The RNA integrity number (RIN) was >8 for all samples, as determined with a 2100 Bioanalyzer (Agilent Technologies).

Construction of sRNA libraries and deep-sequencing

Total RNAs, extracted from developing wheat grains at 5, 15, 25, and 30 DAP, were purified by electrophoretic separations on a 15% TBE-Urea denaturing PAGE gel. The sRNA regions corresponding to 18- to 30-nucleotide bands were excised and recovered. Those sRNAs were then 5′ and 3′ RNA adapter-ligated by T4 RNA ligase. At each step, their lengths were verified before they were purified. The adapter-ligated sRNAs were transcribed into cDNA by Super-Script II Reverse Transcriptase (Invitrogen). We performed PCR-amplifications with primers that annealed to the ends of the adapters. Amplified cDNA constructs were also purified and recovered. The final quality of the cDNA library was ensured by examining its size, purity, and concentration with the 2100 Bioanalyzer. Four libraries, one for each sampling time point, were used for sequencing with the Illumina Genome Analyzer at the Beijing Genomics Institute, Shenzhen, China.

Analysis of sequencing data

Automated base calling of the raw sequence and vector removal were performed with PHRED and CROSS MATCH programs [51, 52]. Any low-quality reads as well as reads with adaptor contamination that were not ligated to any other sequences were discarded. The remaining, high-quality, sequences were trimmed of their adapter sequences; only sequences between 18 nt and 30 nt were analyzed further. The clean reads were mapped to all of the wheat EST sequences obtained from NCBI by SOAP2 (http://soap.genomics.org.cn/). They were aligned to known miRNA precursors (http://www.mirbase.org/) to determine how many identified miRNAs were present. The remaining sequences were used in a BLASTN search of the Rfam database (http://rfam.sanger.ac.uk/) to remove most of the non-siRNA and non-miRNA sequences. Putative origins for the remaining sequences were identified by a BLASTN search against the wheat EST database from NCBI. The protein-coding EST sequences were removed and the remaining non-coding candidate ESTs, with perfect matches to sRNA sequences, were used for fold-back via Mireap (http://sourceforge.net/projects/mireap) based on critical criteria described by Meyers et al. [33]. The secondary structure of precursors were predicted by the Mfold web server (http://mfold.rna.albany.edu/?q=mfold/RNA-Folding-Form) [53].

Expression analysis of miRNAs based on deep-sequencing data

Prediction and annotation of potential miRNA targets

Predictions of putative targets were performed by using the most abundant miRNA variants as queries against the TAGI database at psRNATarget (http://plantgrn.noble.org/psRNATarget/) [35]. The following default parameters were applied: maximum expectation value, 3; target accessibility – allowed maximum energy to unpair the target site (UPE), 25; flanking length around the target site for target accessibility analysis, 17 bp upstream and 13 bp downstream; and the range of central mismatch leading to translational inhibition, 9 to 11 nt. To acquire a better annotation of the potential targets, we performed BLAST searches for their sequences against the NCBI database.

Expressing validation of miRNAs and their targets

We verified the patterns of expression by five conserved wheat miRNAs -- miR167a, miR397, miR156a, miR1852, and miR827a – plus three novel miRNAs -- Ta-miR021-2-5p, Ta-miR004-1-5p, and Ta-miR044-1-3p. A One Step PrimeScript® miRNA cDNA Synthesis Kit (Perfect Real Time; TaKaRa) was used for the RT reactions. The temperature program was adjusted to run for 60 min at 37°C, 5 s at 85°C, and then 4°C. For each miRNA, three biological replicates were performed. qRT-PCR was conducted on a Bio-Rad IQ5 Real-Time PCR Detection System. 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® Premix Ex Taq™ (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 30 s. All reactions were run in three replicates for each sample. The actin gene (GB#: AB181991) served as the endogenous control. We also selected six target genes -- TC383723, TC384445, TC370322, TC400547, TC370450, and DR734203 -- for grain filling-induced miR1037, miR156a, miR172a, miR1855, and miR2093 in order to validate their expression profiles in developing grains via qRT-PCR. All primers are listed in additional file 9.