Abstract
MicroRNAs (miRNAs) and small interfering RNAs are important regulators of plant development and seed formation, yet their population and abundance in the oil crop Brassica napus are still not well understood, especially at different developmental stages and among cultivars with varied seed oil contents. Here, we systematically analyzed the small RNA expression profiles of Brassica napus seeds at early embryonic developmental stages in high-oil-content and low-oil-content B. napus cultivars, both cultured in two environments. A total of 50 conserved miRNAs and 9 new miRNAs were identified, together with some new miRNA targets. Expression analysis revealed some miRNAs with varied expression levels in different seed oil content cultivars or at different embryonic developmental stages. A large number of 23-nucleotide small RNAs with specific nucleotide composition preferences were also identified, which may present new classes of functional small RNAs.
Increasing numbers of small RNAs have been identified as regulatory RNAs that modulate gene expression at both the transcriptional and posttranscriptional levels. Currently, two major classes of endogenous small RNAs have been identified in plants, namely microRNAs (miRNAs) and small interfering RNAs (siRNAs; Brodersen and Voinnet, 2006; Jones-Rhoades et al., 2006; Mallory and Vaucheret, 2006; Chen, 2009; Voinnet, 2009). They are approximately 20 to 24 nucleotides long and have crucial regulatory functions in diverse biological processes. Plant miRNAs are produced from hairpin-shaped precursors and mostly 21 nucleotides long, whereas the majority of plant siRNAs are 24 nucleotides in length and generated from long double-stranded RNA duplexes or transcripts derived from inverted repeat regions (Chen, 2009; Voinnet, 2009). Endogenous plant siRNAs can be divided into several classes, including miRNA-induced transacting siRNAs (ta-siRNAs), heterochromatic siRNAs, natural antisense siRNAs, and millions of unclassified small RNA sequences (Choudhuri, 2009). Ta-siRNAs are in-phase-generated small RNAs whose production depends on the cleavage of non-protein-coding transcripts of ta-siRNA genes by miRNAs (Vazquez et al., 2004; Allen et al., 2005; Yoshikawa et al., 2005; Montgomery et al., 2008a, 2008b). Both miRNAs and ta-siRNAs can induce the cleavage of mRNAs with partially or fully complementary sequences to them. Heterochromatic siRNAs are usually produced from repeats and transposable elements (Xie et al., 2004; Lu et al., 2005; Kasschau et al., 2007; Zhang et al., 2007; Mosher et al., 2008) and inhibit gene expression at the transcriptional level by modulating DNA or histone methylation of homologous regions. Natural antisense siRNAs are generated from natural sense-antisense transcript pairs and modulate the expression of one gene in the pair (Borsani et al., 2005; Katiyar-Agarwal et al., 2006; Held et al., 2008; Zhou et al., 2009; Ron et al., 2010).
MiRNAs and siRNAs play key roles in the regulation of diverse processes related to plant development, such as seed germination (Liu et al., 2007; Reyes and Chua, 2007), organ formation (Mallory et al., 2004; Baker et al., 2005; Wang et al., 2005), developmental timing and patterning (Chuck et al., 2009; Poethig, 2009; Rubio-Somoza et al., 2009; Wollmann and Weigel, 2010), hormonal response (Achard et al., 2004; Mallory et al., 2005; Wang et al., 2005; Reyes and Chua, 2007), and stress resistance (Sunkar et al., 2007; Ruiz-Ferrer and Voinnet, 2009). Tightly regulated space- and time-specific expression patterns have been identified for some miRNAs and siRNAs. The recent application of massively parallel sequencing technology has significantly accelerated the discovery and functional study of small RNAs by producing millions of short reads in a single sequencing run (Moxon et al., 2008; Sunkar et al., 2008; Zhu et al., 2008; Hsieh et al., 2009; Ledandais-Brière et al., 2009; Li et al., 2010; Pantaleo et al., 2010). Nevertheless, except for several model species, the expression patterns and functions of most small RNAs are poorly understood in other plants.
Brassica napus is a major edible oil-producing crop. The oil accumulation of B. napus seeds usually starts around 3 weeks after flowering and peaks after another 3 weeks (Norton and Harris, 1975; Cummins et al., 1993; Eastmond and Rawsthorne, 2000). The seed oil content varies significantly among B. napus cultivars or under different growth environments. So far, there have been only a few reports on small RNA population and function analysis in B. napus (Wang et al., 2007; Xie et al., 2007; Buhtz et al., 2008; He et al., 2008; Pant et al., 2008, 2009; Huang et al., 2010; Wei et al., 2010). In addition, most reports on the oil production mechanisms in B. napus focused on the rapid oil accumulation stage; few studied the gene expression profiles of early embryos, which may have already determined the oil production abilities of B. napus seeds (Li et al., 2006; Sharma et al., 2008; Fu et al., 2009). To systematically identify miRNAs and other small RNAs that may be involved in the regulation of B. napus early embryonic development as well as its seed oil accumulation, we constructed small RNA libraries from the early siliques of high-seed-oil-content and low-seed-oil-content B. napus cultivars and profiled their small RNA expression using the Solexa sequencing technology (GenomeAnalyzer; Illumina). This work identified 50 conserved and 9 new miRNAs in B. napus, among which some miRNAs with developmentally regulated or cultivar-diverged expression patterns. We also validated several conserved targets and identified new targets for some miRNAs with expression differences in the high-seed-oil-content and low-seed-oil-content cultivars. Ta-siRNAs and a group of 23-nucleotide-long small RNAs with characteristic nucleotide composition were also identified and analyzed.
RESULTS
Deep Sequencing of B. napus Small RNAs
To investigate the small RNA expression profiles in B. napus cultivars with different oil contents as well as the impact of environmental factors on small RNA expression, a low-seed-oil-content cultivar (L cultivar; XZ088) and a high-seed-oil-content cultivar (H cultivar; XZ366) were planted in Nanjing, a major B. napus-growing riverside plain land in China, and Lhasa, an altiplanic area where the oil content of B. napus seeds usually increased about 4% as compared with those harvested in Nanjing, respectively (Supplemental Table S1). For each sample, small RNAs from developing siliques at 3, 7, 14, and 21 d after flowering (DAF) were collected and pooled together for sequencing using the Solexa sequencing technology. A total of 18,881,853 reads with lengths of 15 to 30 nucleotides were obtained from the four libraries after removing adaptors and low-quality reads, representing 1,994,899, 2,149,419, 1,232,566, and 1,811,074 nonredundant sequences from the four libraries, respectively. Since the B. napus genome sequence has not been released yet, we mapped these small RNAs to the available bacterial artificial chromosome (BAC) clone sequences of all Brassica family species. As a result, perfect matches were obtained for about 20% of the total small RNAs, representing 15% of nonredundant sequences (Supplemental Table S2).
Similar to the length distribution of small RNAs observed in other plants, the most abundant small RNAs with or without perfect genomic matches were 24 nucleotides long in all four libraries, and the second largest population was the 21-nucleotide RNAs (Supplemental Fig. S1). Unexpectedly, a relatively large number of 22- and 23-nucleotide small RNAs were obtained in our data sets, as compared with previous reports in Arabidopsis (Arabidopsis thaliana) and rice (Oryza sativa; Sunkar et al., 2008; Zhu et al., 2008; Hsieh et al., 2009).
Genomic location analysis showed that 13.2% of the total small RNAs with perfect genomic matches were generated from regions encoding reported B. napus miRNAs or sequences homologous to miRNA genes in other species, and 93.1% of these sequences were 21 nucleotides long (Supplemental Table S3). Small RNAs from unannotated genomic regions counted for 58.4% of the total mapped sequences and were dominated by 24-nucleotide small RNAs, followed by 22- and 23-nucleotide RNAs (Supplemental Table S3). rRNA-derived small RNAs also counted for 19.6% of the total population and were enriched in 21-nucleotide-long sequences (Supplemental Table S3).
Last Nucleotide Preference of 23-Nucleotide Small RNAs
Previous studies have shown that in plants, most miRNA sequences start with uridine (U), whereas the majority of 24-nucleotide siRNAs have adenosine (A) as their 5′ first nucleotide (Mi et al., 2008; Rubio-Somoza et al., 2009; Voinnet, 2009; Czech and Hannon, 2011). The same trends were also observed among our cloned B. napus small RNAs (Supplemental Fig. S2). Unexpectedly, a clear preference for adenosine at the last nucleotide was observed among the 23-nucleotide small RNAs, which counted for about 45.5% of total mapped 23-nucleotide small RNAs. Detailed analysis revealed that such last-nucleotide preference did not present in the sequences of other lengths, and the difference was statistically significant (Fig. 1A).
Nucleotide preferences of small RNAs. A, Percentage of adenosine at the last position of 21- to 24-nucleotide (nt) small RNAs. The asterisk indicates a difference of P < 0.01 (Student’s single-sample two-sided t test). B, Nucleotide preferences of the 24-nucleotide-independent 23-nucleotide small RNAs at each position. C, Percentage of adenosine at each position of the 23-nucleotide small RNAs of 162 Arabidopsis deep sequencing data sets deposited in the NCBI GEO database.
To investigate the cause of such nucleotide composition bias, we compared the genomic locations of the 23-nucleotide small RNAs with the 24-nucleotide ones in detail. Among the 519,119 23-nucleotide small RNAs with perfect genomic matches, 365,343 were substrings of 24-nucleotide small RNAs with identical 5′ or 3′ ends as the 24-nucleotide ones. Therefore, we termed this group of sequences the 24-nucleotide-dependent 23-nucleotide small RNAs. The rest of the 153,776 23-nucleotide small RNAs were 24-nucleotide independent, as they were not degradation products of the 24-nucleotide ones. Nucleotide distribution analysis revealed that the preference for adenosine at the last position was most significant among the 24-nucleotide-independent 23-nucleotide small RNAs, presented in about 51.6% of sequences (Fig. 1B), whereas the occurrence rate of adenosine at the last position of the 24-nucleotide-dependent 23-nucleotide small RNAs was about 40%, only slightly higher than that of sequences with other lengths (Supplemental Fig. S3).
To examine whether the last nucleotide preference of the 23-nucleotide small RNAs also presented among small RNAs from other species, we downloaded 162 Arabidopsis small RNA deep sequencing datasets from the National Center for Biotechnology Information Gene Expression Omnibus (NCBI GEO) database (Barrett et al., 2011) and analyzed their nucleotide compositions. Strong preference for adenosine was also observed among the Arabidopsis 23-nucleotide small RNAs, but not among sequences of other lengths (Fig. 1C; Supplemental Fig. S4).
Identification of Conserved and New miRNAs
We first analyzed the presence of miRNAs in the sequencing data, as miRNAs have been shown to play critical roles in many aspects of plant developmental regulation. Although hundreds of plant miRNAs have been identified, only 46 conserved miRNAs were reported in B. napus by previously computational or experimental studies, according to the miRBase release 16 record (Griffiths-Jones et al., 2008). Therefore, it would be essential to systematically identify both plant-conserved and species-specific miRNAs in B. napus. A total of 36 nonredundant reads with identical sequences to miRBase-recorded plant miRNAs and 14 additional sequences that were highly homologous to miRBase-recorded plant miRNAs had hairpin-shaped precursors and were selected as miRNAs (Supplemental Table S4). Among them, 37 had their complementary short RNAs (miRNA*s) sequenced in our datasets, and only 13 were recorded as B. napus miRNAs in the miRBase. In addition, we also predicted 9 new miRNA genes with well-formed hairpin-shaped secondary structures and cloned miRNA* sequences or confirmed targets (Supplemental Fig. S5). The distribution of small RNAs on these miRNA precursors was centered on the miRNA or miRNA* sequences, with the defined miRNAs or miRNA*s having the highest clone numbers (Supplemental Fig. S5).
Some miRNA*s Accumulated Higher Than Corresponding miRNAs
It has been widely accepted that after processing from precursors, the functional strand of the short double-stranded RNA duplex will be retained as miRNAs, whereas their complementary miRNA* sequences will be quickly degraded (Jones-Rhoades et al., 2006; Chen, 2009; Voinnet, 2009). Therefore, in most cases, the abundance of miRNAs is much higher than that of their corresponding miRNA*s. Unexpectedly, we detected that the abundance of the star sequences of two conserved miRNAs, miR160a* and miR171b*, was about seven times more than that of miRNA160a and miRNA171b (Fig. 2, A and B).
Expression of several miRNAs and their corresponding miRNA*s. A to C, Normalized sequence reads of miR160a (A), miR171b (B), and miR408 (C) and their miRNA*s in the four libraries. D, Northern blot hybridization of miR408 and miR408*.
The sequencing results showed that the expression abundance of miR408* was extremely higher than that of miR408 in all four samples (Fig. 2C). Moreover, the abundance of miR408* was higher in the H cultivar than in the L cultivar, as shown by both the sequencing results and northern-blot hybridization (Fig. 2, C and D).
Differentially Expressed miRNAs between the High-Oil-Content and Low-Oil-Content Cultivars
Among the identified miRNAs, 22 miRNA families had more than five reads and therefore were compared for their expression abundance in the two cultivars (Supplemental Tables S4 and S5; Supplemental Fig. S5). Four miRNA families, namely miR169, miR390, miR394, and miR6028, had consistently higher expression in the L cultivar when planted both in Nanjing and Lhasa (Fig. 3A; 1.5-fold change cutoff). To the contrary, miR408* and miR2111 were more abundant (1.5-fold change cutoff) in the H cultivar in both growth environments. Such expression differences were also confirmed by quantitative reverse transcription (RT)-PCR and northern-blot hybridization (Fig. 3, B and C).
Expression profile of miRNAs between the H and L cultivars. A, Comparison of miRNA abundance between the two cultivars. B and C, Quantitative RT-PCR (B) and northern-blot hybridization (C) analyses of selected differentially expressed miRNAs between the cultivars. Mixed samples containing all examined time points were used. Error bars indicate the sd of three replicates.
Most other miRNAs only had greater than 1.5-fold expression changes between the two cultivars in either Nanjing or Lhasa; thus, it is difficult to conclude whether the expression changes were due to genetic or environmental differences. For example, the sequencing results showed that the expression of miR156 was higher in the H cultivar than in the L cultivar in samples grown in Nanjing but had no difference between the two cultivars in samples collected in Lhasa. Such an expression difference was also confirmed by northern-blot hybridization (Fig. 3C). We also used quantitative RT-PCR to examine the expression of a new miRNA, miR6029. Consistent with the sequencing results, miR6029 had higher expression in the L cultivar than in the H cultivar in both growth areas, but only the expression change in the Nanjing sample was greater than 1.5-fold (Fig. 3B). These experiments confirmed the fidelity of the sequencing results and also indicated that miR156 and miR6029 may relate to B. napus oil production.
Expression Changes of miRNAs during Early Silique Development
As the small RNA sequencing data were obtained by pooling samples from four silique developmental stages together, we further examined the expression dynamics of several miRNAs at different time points. To exclude environmental effects, only samples collected from Nanjing were used in this analysis. Diversified expression patterns were observed for some miRNAs. For example, the expression of miR156 remained consistent in the L cultivar but increased significantly at 21 DAF in the H cultivar (Fig. 4A). Similarly, but reversed, the amount of miR167 increased gradually along with the examined time points and shot up significantly at 21 DAF in the L cultivar (Fig. 4B). The expression pattern of miR6029 in the L cultivar was very similar to that of miR167, but it remained consistently low in the H cultivar (Fig. 4B). On the other hand, the expression of miR390 was high at 3 DAF in both cultivars, decreased quickly at 7 DAF, and remained low until 21 DAF (Fig. 4A). The amount of both miR2111 and miR6028 increased consistently during development in the L cultivar, but miR6028 decreased gradually with time in the H cultivar, whereas miR2111 also increased in the H cultivar (Fig. 4).
Expression changes of miRNAs at different silique developmental stages. A, Quantitative RT-PCR analysis for the relative expression of miR156, miR390, miR2111, and miR6028 during early silique development. The expression level of each high-oil-content cultivar miRNA at 3 DAF is set to 1. B, Northern-blot hybridization of miR156, miR390, miR167, and miR6029 during early silique development. Error bars indicate the sd of three replicates.
Conserved and New Targets of B. napus miRNAs
We used previously reported plant miRNA target prediction criteria (Allen et al., 2005) to search for putative targets of the conserved and new B. napus miRNAs. A total of 346 putative targets were identified for all conserved miRNAs as well as all but two new miRNAs, namely miR6031 and miR6036 (Supplemental Table S6). The targets of most conserved miRNAs were homologous to those reported miRNA targets in Arabidopsis or rice. The newly identified targets with transcription factor functions are listed in Table I.
The newly identified miRNA targets with transcription factor functions
| miRNA | Target Gene | Target Site (Nucleotides to the 5′ End) | Target Gene Annotation | |
| Start | End | |||
| miR156 | TC80211 | 1,035 | 1,051 | MADS box transcription factor |
| TC93255 | 1,163 | 1,180 | ZIM-like 1 (ZML1) | |
| miR164 | TC67424 | 138 | 158 | Variant in methylation 1 (VIM1) |
| TC79889 | 679 | 659 | ZFHD transcription factor | |
| miR172 | TC104052 | 616 | 634 | MYB transcription factor |
| TC76001 | 84 | 103 | WRKY DNA-binding protein | |
| miR390 | TC75148 | 414 | 434 | K box transcription factor |
| miR395 | TC109054 | 36 | 16 | ASIL transcription factor |
| miR399 | TC82291 | 460 | 479 | Nuclear factor Y A3 (NF-YA3) transcription factor |
| miR860 | TC68388 | 903 | 884 | IAA transcription factor |
| miR6029 | TC73674 | 239 | 258 | VirE2-interacting protein 1 (VIP1) |
| miR6034 | TC101312 | 1,025 | 1,043 | Basic helix-loop-helix family protein (BIM1) |
| miRNA | Target Gene | Target Site (Nucleotides to the 5′ End) | Target Gene Annotation | |
| Start | End | |||
| miR156 | TC80211 | 1,035 | 1,051 | MADS box transcription factor |
| TC93255 | 1,163 | 1,180 | ZIM-like 1 (ZML1) | |
| miR164 | TC67424 | 138 | 158 | Variant in methylation 1 (VIM1) |
| TC79889 | 679 | 659 | ZFHD transcription factor | |
| miR172 | TC104052 | 616 | 634 | MYB transcription factor |
| TC76001 | 84 | 103 | WRKY DNA-binding protein | |
| miR390 | TC75148 | 414 | 434 | K box transcription factor |
| miR395 | TC109054 | 36 | 16 | ASIL transcription factor |
| miR399 | TC82291 | 460 | 479 | Nuclear factor Y A3 (NF-YA3) transcription factor |
| miR860 | TC68388 | 903 | 884 | IAA transcription factor |
| miR6029 | TC73674 | 239 | 258 | VirE2-interacting protein 1 (VIP1) |
| miR6034 | TC101312 | 1,025 | 1,043 | Basic helix-loop-helix family protein (BIM1) |
| miRNA | Target Gene | Target Site (Nucleotides to the 5′ End) | Target Gene Annotation | |
| Start | End | |||
| miR156 | TC80211 | 1,035 | 1,051 | MADS box transcription factor |
| TC93255 | 1,163 | 1,180 | ZIM-like 1 (ZML1) | |
| miR164 | TC67424 | 138 | 158 | Variant in methylation 1 (VIM1) |
| TC79889 | 679 | 659 | ZFHD transcription factor | |
| miR172 | TC104052 | 616 | 634 | MYB transcription factor |
| TC76001 | 84 | 103 | WRKY DNA-binding protein | |
| miR390 | TC75148 | 414 | 434 | K box transcription factor |
| miR395 | TC109054 | 36 | 16 | ASIL transcription factor |
| miR399 | TC82291 | 460 | 479 | Nuclear factor Y A3 (NF-YA3) transcription factor |
| miR860 | TC68388 | 903 | 884 | IAA transcription factor |
| miR6029 | TC73674 | 239 | 258 | VirE2-interacting protein 1 (VIP1) |
| miR6034 | TC101312 | 1,025 | 1,043 | Basic helix-loop-helix family protein (BIM1) |
| miRNA | Target Gene | Target Site (Nucleotides to the 5′ End) | Target Gene Annotation | |
| Start | End | |||
| miR156 | TC80211 | 1,035 | 1,051 | MADS box transcription factor |
| TC93255 | 1,163 | 1,180 | ZIM-like 1 (ZML1) | |
| miR164 | TC67424 | 138 | 158 | Variant in methylation 1 (VIM1) |
| TC79889 | 679 | 659 | ZFHD transcription factor | |
| miR172 | TC104052 | 616 | 634 | MYB transcription factor |
| TC76001 | 84 | 103 | WRKY DNA-binding protein | |
| miR390 | TC75148 | 414 | 434 | K box transcription factor |
| miR395 | TC109054 | 36 | 16 | ASIL transcription factor |
| miR399 | TC82291 | 460 | 479 | Nuclear factor Y A3 (NF-YA3) transcription factor |
| miR860 | TC68388 | 903 | 884 | IAA transcription factor |
| miR6029 | TC73674 | 239 | 258 | VirE2-interacting protein 1 (VIP1) |
| miR6034 | TC101312 | 1,025 | 1,043 | Basic helix-loop-helix family protein (BIM1) |
We next performed a 5′-RNA ligase-mediated (RLM)-RACE experiment to examine whether some miRNAs could indeed induce the cleavage of their predicted targets. The putative targets of miR156d were squamosa-promoter-binding protein-like (SPL) family genes. We identified several SPL genes containing a putative miR156d-binding site. Cleavage sites mediated by miR156d were detected in the transcripts of one SPL2 gene (TC82507) and two SPL10 genes (TC98126 and TC94184; Fig. 5A). Similarly, the miRNA::target relationships were also confirmed for miR167c and an ARF8-like gene (TC84482). We identified that transcripts of a basic-Leu zipper transcription factor family gene, a gene encoding VirE2 interacting protein1 (VIP1)-like protein, as the target of a new miRNA, miR6029 (Fig. 5A). It is worth noting that the cleavage of VIP1 transcript occurred at the ninth position of miR6029 instead of the commonly observed 10th position in all 21 RACE experiments. Similarly, higher cleavage frequency was also observed between the pairs of miR156d::TC98126 and miR156d::TC94184 at the ninth position.
Validation and expression of selected miRNA target genes. A, 5′-RLM-RACE analysis of the cleavage on target mRNAs by corresponding miRNAs. B, Quantitative RT-PCR analysis of miR156 target gene TC98126 at different developmental stages and miR6028 target gene CD814980 in different cultivars. Error bars indicate the sd of three replicates.
We also used quantitative RT-PCR to examine the abundance of several miRNAs and their targets. Contradictory to the expression profile of miR156 (Fig. 4), the abundance of its target, TC98126, decreased significantly at 21 DAF in the H cultivar but remained almost unchanged in the L cultivar (Fig. 5B). In another example, the abundance of CD814980, a putative target of a new miRNA, miR6028, was comparable in the H and L cultivars at 3 DAF (Fig. 5B) but decreased significantly in the L cultivar at 21 DAF, which was in accordance with the increased expression of miR6028 in the L cultivar at the same time (Fig. 4).
Identification of miRNA-Derived ta-siRNA Genes
MiRNA-derived ta-siRNA is another important class of regulatory small RNAs in plants, but they have not been identified in B. napus yet. We carried out a genome-wide systemic screen for ta-siRNA genes in our data sets using methods modified from previous reports (Vazquez et al., 2004; Allen et al., 2005; Yoshikawa et al., 2005; Montgomery et al., 2008a, 2008b). Phased small RNA clusters were identified from three genes with putative miRNA targeting sites, namely TC65909, S39172450, and a gene located in BAC (region 44,041–44,663, antisense strand; Supplemental Table S7). All three sequences had high homology to the Arabidopsis ta-siRNA gene 3 (TAS3; Supplemental Fig. S6) and therefore were considered as the B. napus TAS3 genes.
Similar to the observations in Arabidopsis, the B. napus TAS3 gene TC65909 also had two miR390-binding sites located at both the upstream and downstream boundaries of the ta-siRNA production region (Fig. 6A). Consecutively produced phased small RNAs covering both the sense and antisense strands of TC65909 were identified from the majority of putative ta-siRNA production loci, with TAS3 5′ D7(+) as the most abundant sense strand-produced ta-siRNA and TAS3 5′ D6(−) as the most abundant antisense strand-produced ta-siRNA. These sense and antisense strand-generated ta-siRNAs were mainly 21 nucleotides long and could form double-stranded duplexes with two-nucleotide overhangs at their 3′ end, consistent with the previously reported features of ta-siRNAs (Supplemental Fig. S6).
MiRNA-derived ta-siRNAs in B. napus. A, Expression abundance of ta-siRNAs generated from TAS3 gene TC65909. The locations of two miR390 target sites and TAS3 5′ D7(+) are labeled by short bars. B and C, Quantitative RT-PCR (B) and northern-blot hybridization (C) analyses of TAS3 5′ D7(+) during early silique development. D, Validation of cleavage on B. napus TC110025 (ARF3) mRNAs by TAS3 5′ D7(+) using 5′-RLM-RACE. Error bars indicate the sd of three replicates.
We next investigated the biogenesis and functions of these B. napus ta-siRNAs using TAS3 5′ D7(+), whose homolog is also the functional ta-siRNA in Arabidopsis. The abundance of TAS3 5′ D7(+) was highest at 3 DAF, then decreased significantly at 7 DAF, and remained low at the following time points (Fig. 6, B and C). Such an expression pattern was identical to that of miR390 (Fig. 4), suggesting that the expression of TAS3 5′ D7(+) was induced by miR390. Target prediction revealed that several genes of the Auxin Response Factor (ARF) family, namely ARF2 (TC81248), ARF3 (TC110025), and ARF4 (TC100380), were putative targets of TAS3 5′ D7(+). We were able to confirm the cleavage of ARF3 by TAS3 5′ D7(+) using 5′-RLM-RACE, with most identified cleavage occurring at the 10th position of the TAS3 5′ D7(+)-binding site (Fig. 6D).
DISCUSSION
MiRNAs have emerged as a new class of regulatory factors and attracted much attention during the past decade. As B. napus is an important crop for edible oil production, understanding the functions of miRNAs in regulating oil accumulation in B. napus seeds will be of great value for the cultivation of B. napus cultivars with increased oil content. In this work, we examined the expression profiles of miRNAs in high-oil-content and low-oil-content B. napus cultivars. Prior to this work, only 46 miRNAs were reported in B. napus, according to the miRBase release 16 record. Here, we enlarged the known population of B. napus miRNAs by identifying 37 unreported conserved miRNAs and 9 new miRNAs, as well as the star sequences for some miRNAs.
We identified putative target genes for 57 miRNAs using computational methods. The lack of predicted targets for other miRNAs is probably due to the unavailability of the B. napus genome sequence. The target genes of most conserved miRNAs have also been identified as miRNA targets in Arabidopsis, suggesting the functional conservation of these miRNAs. We confirmed that miRNA156 regulated SPL family genes and miR167 targeted transcripts of the ARF family in B. napus. We also identified and confirmed transcripts encoding VIP1 as the target of a new miRNA, miR6029. VIP1 is a member of basic-Leu zipper transcription factors and can regulate many genes by binding to their promoter regions (Pitzschke et al., 2009). It was not reported to be regulated by miRNA previously. Our 5′-RLM-RACE experiment obtained consistent miR6029 cleavage sites on VIP1 mRNAs in 21 individual sequencing reactions, demonstrated that VIP1 is an authentic target of miR6029. The pairing between miR6029 and the VIP1 transcript contained two mismatches and two insertions, further proving that gaps are also tolerable between plant miRNA and target pairs. In addition, all cleavage of VIP1 occurred at the ninth position of the miRNA sequence, instead of the commonly observed 10th position, probably due to the presence of an insertion at the 11th position of the pair.
It is well known that many miRNAs have time- and space-specific expression patterns (Jones-Rhoades et al., 2006; Mallory and Vaucheret, 2006; Chen, 2009; Chuck et al., 2009; Poethig, 2009; Rubio-Somoza et al., 2009; Wollmann and Weigel, 2010), yet the expression profiles of miRNAs in the early developmental stage of plant embryos are still poorly understood, especially in the oil production crop B. napus. By comparing the small RNA expression profiles of 3- to 21-DAF embryos using the massively parallel sequencing technology, we have identified some differentially expressed miRNAs between the high-oil-content and low-oil-content B. napus cultivars as well as some miRNAs with expression changes at different developmental stages within the same cultivar, and we have further confirmed their expression differences using northern-blot hybridization and quantitative RT-PCR experiments. The similar or reversed expression patterns of the selected miRNAs examined at individual time points suggested that they may have synergistic or antagonistic functions.
By comparing the high-oil-content and low-oil-content samples, we have identified and confirmed that miR156, miR167, miR390, miR2111, and two new miRNAs, miR6028 and miR6029, as well as the TAS3-produced ta-siRNAs were differentially expressed between the two cultivars in one or both growth environments. According to the varied expression levels of these small RNAs at different developmental stages and the functions of their targets, we propose that miR390 is involved in early embryo development regulation via the ta-siRNA-mediated auxin signaling pathways, whereas miR156, miR167, and miR6029 function around 21 DAF by regulating SPL, ARF, and VIP1 transcripts, which may affect the oil content of B. napus seeds in later stages through multiple pathways (Supplemental Fig. S7).
The plant genomes usually express large numbers of small RNAs, the majority of which were not classified yet. Unlike the small RNA profiles of most reported plant samples, a large number of 23-nucleotide small RNAs were identified in our data. A total of 29.6% of these 23-nucleotide small RNAs were not potential degradation products of 24-nucleotide small RNAs and exhibited a strong preference for adenosine at the 23-nucleotide position, indicating that they may be an uncharacterized class of functional small RNAs.
CONCLUSION
To the best of our knowledge, this work provides the first small RNA expression analysis in B. napus early embryos. It enlarged the population of known B. napus miRNAs and proved the presence of ta-siRNAs in B. napus. The comparison between high-oil-content and low-oil-content cultivars revealed some miRNAs that may be involved in the regulation of B. napus seed oil production. The detailed expression analysis of some miRNAs and their targets at different developmental time points provided clues for their varied roles in regulating the early development of B. napus seed, which may be conserved in other dicot plants as well. The 23-nucleotide small RNAs found in this work suggested that there may be other types of uncharacterized small RNAs remaining to be discovered.
MATERIALS AND METHODS
Silique Collection, RNA Isolation, and Small RNA Sequencing
Two Brassica napus cultivars (XZ088 and XZ366) were grown in the experimental fields in the plain land Nanjing and the altiplanic area Lhasa, respectively. Siliques of both cultivars were collected at 3, 7, 14, and 21 DAF separately. For each time point, the siliques of more than 20 plants were collected and pooled together to avoid the individual genotype differences. Total RNAs of each sample were isolated using Trizol reagent (Invitrogen). Equal amounts of total RNAs of the 3-, 7-, 14-, and 21-DAF samples of each cultivar were mixed before sending for sequencing. Small RNA sequences of XZ088 and XZ366 siliques collected in both growth environments were generated by Illumina’s Solexa sequencing technology.
Primary Analysis of the Deep Sequencing Data Sets
After trimming adapter sequences and removing sequences with low quality or length less than 15 nucleotides, the resulting small RNA sequences were compared with the Rfam database (Griffiths-Jones et al., 2005) and the NCBI nucleotide database to classify tRNA, rRNA, small nucleolar RNA, and protein-coding gene-derived small RNAs. The BAC clone sequences of B. napus, Brassica rapa, and Brassica oleracea collected in the NCBI nucleotide database by the BLASTN program (Altschul et al., 1990) were used to analyze the genomic origins of the remaining small RNAs. The BAC clone sequences of B. rapa and B. oleracea were included because the B. napus genome (AC genome) is a polyploid of the B. rapa genome (A genome) and the B. oleracea genome (C genome). All BAC clone sequences were first analyzed using the RepeatMasker software to identify transposons and repeat regions. Small RNA sequences were aligned to the BAC clones by the BLASTN program (Altschul et al., 1990), and only sequences with perfect genomic matches were retained for further analysis.
Identification of miRNAs and ta-siRNAs
Small RNAs of lengths between 20 and 24 nucleotides, with perfect matches to the nonrepeat-associated noncoding regions of BAC clones, as well as with no more than seven mapping loci on nonredundant BAC sequences, were used in miRNA prediction. According to the community-established criteria (Meyers et al., 2008), small RNA sequences satisfying all the following rules were selected as miRNAs: (1) intergenic location of miRNA candidates; (2) hairpin-shaped secondary structures of precursors; (3) stem region location of mature miRNAs on precursors; (4) fewer than seven gaps and mismatches between mature miRNAs and their pairing sequences; (5) the presence of miRNA* sequences or the confirmation of target cleavage by a small RNA. Putative precursor sequences of each small RNA were extracted by assuming that the small RNAs can be derived from either the 5′ or 3′ arm of the hairpin-shaped precursors, with a 20-nucleotide extension at either the 5′ or 3′ end of the small RNA loci and a 60- to 200-nucleotide extension (20-nucleotide increment each time) at the other end. The secondary structure of each putative small RNA precursor was evaluated by the mfold software (Markham and Zuker, 2008). Candidate miRNA precursor sequences were compared with each other using the BLASTN program to remove redundant sequences from different BAC clones. The conservation status of the identified miRNAs was analyzed by comparing their sequences with all plant miRNAs collected in miRBase (release 16).
In the process of ta-siRNA identification, all Brassica BACs, B. napus tentative consensus sequences, and UniGene sequences were screened by a 200-nucleotide width window, with a 10-nucleotide increment at each step. The number of phased small RNAs with lengths of 21 nucleotides within each window was calculated, and in-phase-produced small RNAs starting from putative miRNA target sites were identified as ta-siRNAs according to previously described criteria (Vazquez et al., 2004; Allen et al., 2005; Yoshikawa et al., 2005; Montgomery et al., 2008a, 2008b). Redundant genes were merged by the BLASTN program.
Target Gene Prediction for miRNAs and ta-siRNAs
The tentative consensus sequences from the Dana Farber Cancer Institute (DFCI) Gene Index Project of B. napus were downloaded and used for miRNA and ta-siRNA target prediction. MiRNAs and ta-siRNAs were aligned to these sequences using Smith-Waterman algorithm-based C programs (Smith et al., 1981). Putative miRNA- and siRNA-binding sites were identified using a scoring system modified from a previous report (Allen et al., 2005) as follows: (1) a penalty score of 1 for each mismatch or gap between the pairing of miRNAs and transcripts; (2) a penalty score of 0.5 for each G::U pair; (3) doubled penalty scores for mismatches, gaps, or G::U pairs located within the two to 13 positions of miRNAs. A total penalty score of no more than 5 was used as the cutoff to select miRNA and ta-siRNA targets.
Small RNA Northern-Blot Hybridization
Small RNA northern-blot hybridization was carried out following previously reported procedures (Liu et al., 2005) using 20 μg of total RNA in each experiment. The complementary sequences of small RNAs were labeled by [γ-32P]ATP and used as probes. The sequences of all probes used are listed in Supplemental Table S8.
Quantitative RT-PCR
The stem-loop RT-PCR method was used in the quantitative RT-PCR experiments of small RNAs following previously reported procedures (Chen et al., 2005; Varkonyi-Gasic et al., 2007). In the quantitative RT-PCR experiments of small RNAs, 3 μg of total RNAs treated with DNase I (New England Biolabs) was used in each reaction with U6 snRNA as the internal control. These reactions were performed using previously reported procedures (Chen et al., 2005; Varkonyi-Gasic et al., 2007). In the quantitative RT-PCR experiments of mRNAs, 1 μg of total RNAs treated with DNase I (New England Biolabs) was synthesized to cDNA by Moloney murine leukemia virus (New England Biolabs) using poly(dT) oligonucleotides with the mRNA of Bn-EF1 (for elongation factor 1-α; TC111637) as the reference gene. SYBR Green PCR Master Mix (Applied Biosystems) was used in all quantitative RT-PCR experiments. The relative fold expression changes of miRNAs and genes were both calculated using the 2 δ-δ Ct method (Livak and Schmittgen, 2001). Primers used in all quantitative RT-PCR experiments are listed in Supplemental Table S9.
5′-RLM-RACE
The RNA 5′-RLM-RACE experiments were carried out using the FirstChoice RLM-RACE Kit (Ambion) following previously reported procedures (Lu et al., 2008). All gene-specific primers used in the experiments are listed in Supplemental Table S10.
Sequence data from this article can be found in the GenBank/EMBL data libraries under accession number GSE34727.
Supplemental Data
The following materials are available in the online version of this article.
Supplemental Figure S1. Size distribution of small RNAs.
Supplemental Figure S2. Nucleotide preference at each position of small RNAs.
Supplemental Figure S3. Nucleotide preference at each position of the 24-nucleotide-dependent 23-nucleotide small RNAs.
Supplemental Figure S4. Percentage of adenosine at each position of 21-nucleotide, 22-nucleotide, and 24-nucleotide small RNAs of the 162 Arabidopsis small RNA deep sequencing data sets deposited in the NCBI GEO database.
Supplemental Figure S5. Secondary structures of the B. napus miRNA candidate precursors, and the expression and location of small RNAs mapped to these precursors.
Supplemental Figure S6. Size distribution of B. napus ta-siRNAs, and the sequence conservations of B. napus TAS3 genes.
Supplemental Figure S7. A proposed model of the regulation of B. napus early silique development via miRNAs, ta-siRNAs, and their target genes.
Supplemental Table S1. Seed oil contents of the two cultivars planted in Nanjing and Lhasa.
Supplemental Table S2. Summary of small RNA deep sequencing data.
Supplemental Table S3. The mapping loci distributions of small RNAs with different lengths.
Supplemental Table S4. Conserved B. napus miRNA candidates.
Supplemental Table S5. New B. napus miRNA candidates.
Supplemental Table S6. Predicted target genes of the B. napus miRNAs.
Supplemental Table S7. Ta-siRNA genes in B. napus.
Supplemental Table S8. Small RNA northern-blot hybridization probes used in this study.
Supplemental Table S9. Quantitative RT-PCR primers used in this study.
Supplemental Table S10. 5′-RLM-RACE primers used in this study.
ACKNOWLEDGMENTS
We thank Ms. Nima Zhuoma (Tibet Academy of Agricultural and Animal Husbandry Sciences) for help with plant culturing and sample collection in Lhasa.
LITERATURE CITED
Author notes
This work was supported by the National Natural Science Foundation of China (grant nos. 30725014 and 30921016 to X.-J.W. and grant no. 30800632 to M.W.) and the National High Technology Research and Development Program of China (grant no. 2009AA101105 to C.-K.Q.).
These authors contributed equally to the article.
These authors contributed equally to the article.
Corresponding author; e-mail xjwang@genetics.ac.cn.
The author responsible for distribution of materials integral to the findings presented in this article in accordance with the policy described in the Instructions for Authors (www.plantphysiol.org) is: Xiu-Jie Wang (xjwang@genetics.ac.cn).
The online version of this article contains Web-only data.
Open Access articles can be viewed online without a subscription.
