A B. napus double haploid line (F117) with stable oil content over 3 years was used in this study. The plants were grown under natural conditions in the experimental field of the Chongqing Engineering Research Center for Rapeseed, Southwest University in Beibei, Chongqing, China (106.40°E, 29.80°N) from October 2014 to May 2015. Developing seeds from different F117 plants were collected in the middle of a light cycle at 14, 21, and 28 DAF and immediately frozen in liquid nitrogen and stored at −80°C until use. Total RNA was isolated using TRIzol H (Invitrogen, USA) according to the manufacturer's instructions, and the RNA quality was evaluated by electrophoresis on a 1% agarose gel (Han et al., 2014). Total RNA (>10 μg) was sent to Beijing Genome Institute (BGI, Shenzhen, China) for sRNA library construction and Solexa sequencing using standard protocols with the Illumina HiSeq 2000 platform.

Small RNA libraries were constructed and sequenced for the three stages (14, 21, and 28 DAF); all raw sequences were filtered with the SOAPnuke software (http://soap.genomics.org.cn/; Li et al., 2009). Low-quality reads, reads smaller than 18 nt, adaptor sequences, and contamination by adaptor–adaptor ligation were removed according to the software's default settings. The raw sequences were categorized to unique reads and annotated using the Rfam database (http://www.sanger.ac.uk/software/Rfam) and the GenBank non-coding RNA database (http://www.ncbi.nlm.nih.gov/). Small RNAs were then aligned to miRNA precursors of rapeseed in miRBase 21.0 (Kozomara and Griffiths-Jones, 2014), and the expression of known miRNAs was assessed.

To identify novel miRNAs, the software Mireap (http://sourceforge.net/projects/mireap/) developed by BGI was used to predict the unannotated small RNA reads mapping to the B. napus genome. A small RNA was regarded as a novel miRNA candidate if it met certain criteria described previously (Wang et al., 2011; Ding et al., 2012). Potential targets for the miRNAs were predicted using the psRobot software with default parameters (Wu et al., 2012). A previously defined scoring system was used to evaluate all predicted target genes, and genes with a score less than 3.0 were considered miRNA targets (Srivastava et al., 2014).

To better understand miRNA target functions and classifications as well as the metabolic regulatory networks associated with B. napus miRNAs and their targets, all target genes were mapped to Gene Ontology (GO) terms (http://www.geneontology.org/), and the number of genes for each term was calculated. To identify significantly enriched GO terms, a hypergeometric test was utilized to compare the target gene candidates with the reference gene background to determine the P-value (Sha et al., 2014). GO terms with a P-value less than the threshold of 0.05 were considered to be significantly enriched. GO annotation results were plotted using WEGO (http://wego.genomics.org.cn/cgi-bin/wego/index.pl). Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg) was used to analyse metabolic pathway assignments. The test and threshold values for estimating significantly enriched metabolic pathways and signal transduction pathways were the same as those used for the GO analysis (Geng et al., 2015).

Quantitative real-time PCR (qRT-PCR) for miRNAs and their targets was performed using a CFX96 Real-time System (BIO-RAD, USA). Total RNA from each sample was extracted as described above. Briefly, 1 μg of RNA from each sample was used to generate single-stranded miRNA cDNA by reverse transcription with miRcute miRNA First-Strand cDNA synthesis Kit (TIANGEN, Beijing, China) and miRNA-specific primers provided with the kit. Next, the expression levels of miRNAs involved in fatty acid biosynthesis were analyzed in three seed developmental stages using qRT-PCR and miRNA-specific primers with a CFX96 Real-time System (BIO-RAD, USA) and SYBR® Premix (TIANGEN, Beijing, China). U6 snRNA was used as the reference gene in qRT-PCR.

Predicted target genes were validated by quantitative RT-PCR using specific primers designed with the software Primer Premier 5.0 (PREMIER Biosoft Int., Palo Alto, CA, USA). qRT-PCR was performed with a CFX96 Real-time System (BIO-RAD, USA) using SYBR® Premix (TIANGEN, Beijing, China). Actin7 was used as an endogenous control. All samples were subjected to three technical replicates.

Deep sequencing of three small RNA libraries from developing B. napus seeds produced 12,225,750, 11,419,839, and 11,427,691 raw sequence reads. After the removal of low-quality reads and 3′ adapter, 5′ adapter, corrupted adapter (reads < 10 nt or >30 nt long) and other contaminating sequences, 12,120,056 (99.49%), 11,334,399 (99.63%), and 11,335,373 (99.57%) clean reads were obtained from the 14, 21, and 28 DAF libraries, respectively (Table S1). After the further removal of unannotated small RNAs and non-coding RNAs, such as tRNAs, rRNAs, siRNAs, snRNAs, snoRNAs and other non-coding RNAs, 393,346 (3.25%), 746,228 (6.59%), and 1,083,239 (9.56%) miRNA sequences were identified in the three libraries (14, 21, and 28 DAF, respectively; Table 1). The meaningful feature of the size profile permitted the miRNAs to be distinguished from other small RNAs. The miRNA length distribution (18–28 nt) of the original reads revealed that those 20–24 nt in length were the most abundant (Figure 1).

By mapping unique sRNA sequences to miRBase 21.0 with a maximum of two mismatches, a total of 85 unique sequences belonging to 30 known miRNA families were identified in the three libraries. Among the known miRNA families, seven, six, and seven members were found from the miR156, miR166, and miR171 families, respectively. As the most abundant in the 21 DAF library, 10 members of the miR169 family were identified. In addition, only one member was found for 10 miRNA families (Figure 2A; Table S2).

The number of reads differed drastically among the 30 known miRNA families. Extraordinarily high expression levels of a few known miRNA families, such as miR156, miR166, and miR167, were identified in all three libraries. MiR156 was the most abundant, with 343,028 (14 DAF), 1,572,529 (21 DAF), and 2,959,044 (28 DAF) reads accounting for 27.4, 59.7, and 74.1% of all known miRNA reads, respectively (Figure 2B). Several miRNA families, including miR164, miR168, and miR390, exhibited moderate abundance. In contrast, a few known miRNA families, such as miR161, miR393, miR2111, miR395, miR396, miR399, miR6035, and miR6036, showed relatively lower expression levels and were represented by < 50 reads in the three libraries. Among these miRNAs, 69 miRNAs were expressed at all three developmental stages, with only 1, 5, and 3 co-expressed at 14 and 21 DAF, 14 and 28 DAF, 21 and 28 DAF, respectively. For example, bna-miR2111c is stage-specifically expressed only at 14 DAF, bna-miR169i, bna-miR169j and bna-miR169l only at 21 DAF, and bna-miR395a, bna-miR395b and bna-miR395c only expressed at 28 DAF. Sixty-nine miRNAs were expressed at all three developmental stages, some of which demonstrated little variation throughout seed development, suggesting that they perhaps fulfill housekeeping functions.

To predict novel miRNAs, Mireap was used with strict criteria (Li et al., 2013) that include the characteristic hairpin structures of miRNA precursors, Dicer cleavage sites, and minimum free energy. In total, 1610 novel miRNAs were predicted from the three libraries; the lengths of the novel miRNAs ranged from 20 to 24 nt, with 24 nt being the most common in all three libraries (Table S3; Figure S1). More than half of the novel predicted miRNAs begin with a 5′ uridine, and these miRNAs accounted for more than 80% of 20 and 21 nt small RNAs. Compared with known miRNA families, the abundance of novel miRNAs was very low, and the majority were present in less than 50 reads. Nonetheless, these miRNAs comprised 83.26% (378/454), 91.01% (769/845), and 92.03% (924/1004) of the 14 DAF, 21 ADF and 28 DAF libraries, respectively. The most abundant novel miRNA was novel_mir_146, which was sequenced in 10,168, 13,227, and 16,128 reads of the 14, 21, and 28 DAF libraries, respectively (Table S2). Unlike known miRNAs, different types of novel miRNAs were expressed in the three independent libraries: 106 miRNAs were expressed at all three developmental stages; 40, 32, and 177 were co-expressed at 14 and 21 DAF, 14 and 28 DAF, 21 and 28 DAF, respectively; 200, 446, and 609 were stage specifically expressed at 14, 21, and 28 DAF, respectively (Figure 3); 106 were expressed at all three developmental stages, some of which showed little variation throughout seed development, suggesting housekeeping functions for these miRNAs.

The identification of miRNA target genes using bioinformatic approaches is essential for understanding the regulatory function of miRNAs (Zhao et al., 2013). In this study, the software psRobot was used to predict miRNA targets of known and novel miRNAs; using a cut-off threshold of 3.0, 2582, and 10,032 putative targets were found, respectively (Table S4). Conversely, no target genes were predicted for the remaining 180 novel miRNAs. Among the known miRNAs identified in our analysis, bna-miR156d, bna-miR156e, and bna-miR156f have 134 putative target genes with different functions, indicating that these three miRNAs are involved in regulating the expression of multiple genes in B. napus.

Using the criteria of an absolute fold change value ≥1.0 and a P-value ≤ 0.05, 702, and 509 miRNAs showed significantly different expression between the 14 and 21 DAF as well as 21 and 28 DAF libraries, respectively. Comparing the 14 and 21 DAF libraries indicated 5013 significantly altered genes, with 3283 up-regulated and 1730 down-regulated, and 1873 genes were detected between the 21 and 28 DAF libraries, with 1170 up-regulated and 703 down-regulated (Tables S2, S4). GO analysis was used to classify the functions of the target genes of the miRNAs differentially expressed during seed development based on the three main categories: biological process, cellular component, and molecular function). For 14 vs. 21 DAF, the target genes are involved in 12 different molecular functions, 22 biological processes, and 10 cellular components, and for 21 vs. 28 DAF, 23 biological processes, 11 molecular functions, and 10 cellular components were identified (Figure 4). Many biological processes were found to be involved, including cellular process (GO: 0009987), biological regulation (GO: 0065007), and metabolic process (GO: 0008152). Figure 4 shows up-regulated target genes specifically enriched during different seed developmental stages, involving amino acid biosynthesis, pigment accumulation, embryonic development and others.

To illuminate the relationship between miRNA and putative gene function, we constructed miRNA regulatory networks for fatty acid biosynthesis, pigment accumulation, embryonic development, sugar conversion, amino acid metabolism, plant hormones, and signaling pathways during seed development (Figure S2). Following these analysis, KEGG pathway analysis identified significant enrichment of 85 pathways with predominant enrichment of seed development-related pathways (Table S5).

To examine the molecular mechanism of fatty acid biosynthesis during seed development, we investigated target genes related to the fatty acid biosynthesis pathway of differentially expressed miRNAs. As shown in Figure 5, de novo synthesis of fatty acids utilizes acetyl-CoA as a substrate and malonyl-ACP as an elongator, and 27 targets, which encode 10 catalytic enzymes, are involved in plastid acetyl-CoA conversion to fatty acids. The formation of malonyl-CoA from acetyl-CoA and bicarbonate by acetyl-CoA carboxylase (ACC) has long been considered a key regulatory step of fatty acid biosynthesis (Turnham and Northcote, 1983; Harwood, 1996), and the miRNA target BnaA06g06030D encodes the carboxyltransferase alpha subunit of acetyl-CoA carboxylase (α-CT, EC: 6.4.1.2) and the biotin carboxylase subunit (BC, EC: 6.3.4.14). The malonyl-CoA produced by plastidial ACC constitutes the carbon donor for each cycle of the fatty acid biosynthesis pathway (Hannapel and Ohlrogge, 1988; Bonaventure and Ohlrogge, 2002). Malonyl-thioester undergoes a series of condensation reactions with acetyl-CoA percycle, steps that are catalyzed by 3-ketoacyl-ACP synthase of type III (KASIII, EC: 2.3.1.179), 3-ketoacyl-ACP reductase (KAR, EC: 1.1.1.100), 3-hydroxyacyl-ACP dehydratase (HAD, EC: 4.2.1.59), and enoyl-ACP reductase (ER, EC: 1.3.1.9), to produce a saturated fatty acid with two additional carbons. Among the four steps, 14 targets were found for three steps (EC: 2.3.1.179, EC: 1.1.1.100 and EC: 1.3.1.9); four targets (BnaA03g37760D, BnaA06g13360D, BnaC01g32050D, and BnaC05g14920D) encode KASIII, eight targets encode KAR, and two targets (BnaA03g38220D and BnaA07g04370D) encode enoyl-ACP reductase. However, none of the targets of differentially expressed miRNAs encode 3-hydroxyacyl-ACP dehydratase. After 7 cycles, the saturated 16-carbon acyl-ACP can either be hydrolysed by FATB acyl-ACP thioesterase (EC: 3.1.2.14; BnaA06g04900D and BnaAnng26510D) or further elongated by KASII to generate 18:0-ACP, which is then desaturated to 18:1-ACP by stearoyl-acyl-carrier-protein desaturase (EC: 1.14.19.2; BnaA05g33500D and BnaC05g48250D) and hydrolysed by FATA thioesterase (EC: 3.1.2.14; BnaA06g04900D and BnaAnng26510D) (Bates et al., 2013).

We further explored differentially expressed miRNAs associated with fatty acid biosynthesis, fatty acid desaturation, and fatty acid elongation pathways during seed development. As shown in Figure 6, 24 miRNAs (Table 2), which regulate 10 catalytic enzyme-encoding genes, are involved in the plastid fatty acid biosynthesis pathway. Among these catalytic enzymes, expression of KASII and KASIII is regulated by known the miRNA bna-miR159. The gene encoding KAR is regulated by known the miRNAs bna-miR156b, bna-miR156c, bna-miR156g, and bna-miR6029. Novel miRNAs are found at all 10 steps of the fatty acid biosynthesis pathway, among which KAR and fatty acyl-ACP thioesterases B (FATB) with five miRNAs have the greatest number of regulating miRNAs. In contrast, acetyl Co-enzyme, carboxylase biotin carboxylase subunit (BC), the carboxyltransferase alpha subunit of acetyl-CoA carboxylase (α-CT), malonyl CoA-acyl carrier protein transacylase (MCMT), KASI, and enoyl-ACP reductase (ER) are regulated by a single miRNAs. According to our sequencing results, the expression levels of bna-miR159, novel_mir_19, novel_mir_555, novel_mir_702, and novel_mir_2163 increased significantly, possibly indicating positive regulatory roles for these miRNAs. However, the expression levels of other miRNAs were negatively correlated with the content and composition of fatty acids during the middle and late seed developmental stages, indicating negative regulatory roles. In addition, 16 miRNAs regulating 5 catalytic enzyme-encoding genes are involved in the fatty acid desaturation and fatty acid elongation pathways (Figure S3); steadily up-regulated at 21 DAF and 28 DAF, bna-miR395d, bna-miR395e, and bna-miR395f, which regulate 3-ketoacyl-CoA synthase (KCS), may be involved in fatty acid elongation.

Compared with the size of other miRNA precursors (typically 40–200 nt) reported in a previous study (Xie et al., 2007), the novel miRNA precursors in our work are more diverse in structure but smaller in size (Table 2; Figure 7). The length of miRNA precursors involved in fatty acid biosynthesis varied from 76 to 101 nt, with an average of 89.9 ± 7.1, and approximately 89.5% of these precursors are 80–100 nt in length. The differences in size of the identified miRNAs within different families suggest that they may carry out unique functions in regulating miRNA biogenesis or gene expression (Zhang et al., 2006). The more diverse in structure include novel_mir_702 and novel_mir_2163, which are simultaneously located at the 5′ and 3′ ends of miRNA precursors; in contrast, novel_mir_1407, novel_mir_1706, novel_mir_173, novel_mir_1758, and novel_mir_555 are located at the 5′ end and the others at the 3′ end of miRNA precursors (Figure 7).

In addition to the miRNAs described above, which may regulate functional genes directly involved in KEGG pathway fatty acid biosynthesis, certain other miRNAs are indirectly involved in fatty acid metabolism by regulating a large number of transcription factors. To understand the regulatory mechanism of miRNAs indirectly involved in fatty acid metabolism, we constructed a miRNA-mediated gene regulatory network for 31 miRNAs and their 11 targets (Figure 8). We analyzed the connection distribution of the network and found that SPL9 and ZFP have the highest number of connections (8 and 6, respectively); ZFP is co-regulated by the bna-miR2111 family, miR172 family and novel_mir_1758, and SPL9 is mainly regulated by the miR156 family. Interestingly, three members of the miR156 family (bna-miR156b, bna-miR156c, and bna-miR156g) directly participate in the regulation of fatty acid biosynthesis; this was also found for transcription factor PEX, which is regulated by bna-miR159 and four novel miRNAs. In addition, the miR172 family regulates targets AP2 and TOE2, and novel_mir_1758 participates in the regulation of GL2 and a mitochondrial substrate carrier family protein. As a key connection, novel_mir_104 regulates five targets. Co-regulated targets of different novel miRNAs can also be observed. Of 11 targets, 10 are regulated by novel miRNAs; more novel miRNAs are involved in the miRNA-mediated gene regulatory network of fatty acid metabolism.

To confirm the sequencing results and examine the dynamic expression patterns of the miRNAs involved in fatty acid biosynthesis at different stages of seed development (14, 21, and 28 DAF) in B. napus, the expression patterns of five known and 11 novel miRNAs and their corresponding predicted targets were validated by qRT-PCR (Tables S6, S7; Figure S4). As expected, the qRT-PCR data of miRNAs showed a high degree of agreement with the expression profiles obtained by sRNA sequencing of the 14, 21, and 28 DAF libraries, and the qRT-PCR data of the corresponding target gene was just the opposite (Table S2; Figure 9). Regarding known miRNAs, transcripts of bna-miR159 and bna-miR395d were up-regulated, steadily increasing at 21 and 28 DAF; this was particularly obvious for bna-miR159, which exhibited expression changes greater than 2.2-fold between 14 DAF and 28 DAF (p < 0.01). By contrast, the gene encoding KASII, FAB1 (BnaA06g13360D), a target of bna-miR159, was down-regulated, gradually declined at 21 and 28 DAF, and had greater than 2-fold expression changes between 14 and 28 DAF. Moreover, bna-miR6029 first sharply declined at 21 DAF and then increased at 28 DAF. Novel_mir_1407, novel_mir_173, novel_mir_1706 and novel_mir_104 were down-regulated, remaining at an extremely low expression level at the three stages. The other miRNAs were sharply increased at 21 DAF and then sharply declined at 28 DAF, except that novel_mir_555 was slightly increased at 28 DAF (p < 0.001). On the contrary, BnaA01g09630D first sharply increased at 21 DAF (p < 0.01) and then declined at 28 DAF. BnaA03g13780D, BnaA05g33500D, BnaA06g06030D and BnaA03g37760D were up-regulated, remaining at an extremely high expression level at the three stages; among these, BnaA03g37760D was up-regulated at 28 DAF to 5 times the level at 14 DAF, BnaA06g06030D and BnaA03g13780D were in the same situation.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.