High throughput sequencing of small RNAs reveals dynamic microRNAs expression of lipid metabolism during Camellia oleifera and C. meiocarpa seed natural drying

Background Camellia species are ancient oilseed plants with a history of cultivation over two thousand years. Prior to oil extraction, natural seed drying is often practiced, a process affecting fatty acid quality and quantity. MicroRNAs (miRNA) of lipid metabolism associated with camellia seed natural drying are unexplored. To obtain insight into the function of miRNAs in lipid metabolism during natural drying, Illumina sequencing of C. oleifera and C. meiocarpa small-RNA was conducted. Results A total of 274 candidate miRNAs were identified and 3733 target unigenes were annotated by performing a BLASTX. Through integrated GO and KEGG function annotation, 23 miRNA regulating 131 target genes were identified as lipid metabolism, regulating fatty acid biosynthesis, accumulation and catabolism. We observed one, two, and four miRNAs of lipid metabolism which were specially expressed in C. Meiocarpa, C. oleifera, and the two species collectively, respectively. At 30% moisture contents, C. meiocarpa and C. oleifer produced nine and eight significant differentially expressed miRNAs, respectively, with high fatty acid synthesis and accumulation activities. Across the two species, 12 significant differentially expressed miRNAs were identified at the 50% moisture content. Conclusions Sequencing of small-RNA revealed the presence of 23 miRNAs regulating lipid metabolism in camellia seed during natural drying and permitted comparative miRNA profiles between C. Meiocarpa and C. oleifera. Furthermore, this study successfully identified the best drying environment at which the quantity and quality of lipid in camellia seed are at its maximum. Electronic supplementary material The online version of this article (doi:10.1186/s12864-017-3923-z) contains supplementary material, which is available to authorized users.


Background
C. oleifera and C. meiocarpa belong to the camellia family, Theaceae, and are known for their high quality oilseed that is dubbed as the Eastern olive oil [1,2]. Some of the camellia species are ancient oilseed plants with history of cultivation for over two thousand years [3]. Camellia oil is known for its edible and medicinal uses with an oleic acid content reaching above 80%, a high content of monounsaturated lipid, and the lowest level of saturated fats [4].
Camellia oil aids in cholesterol reduction, resistance to stress, oxidation reduction, reduced inflammation, and improved human immunity [5]. The drying management of camellia seed prior to oil extraction is a fundamental factor affecting its fatty acid quality and quantity [6][7][8].
MicroRNAs (miRNAs) are small regulatory molecules that have been shown to be involved in a wide range of biological pathways by modulating expression of specific mRNAs [9]. Sequencing small RNA libraries demonstrated the role of miRNAs in lipid metabolism in plants (safflower [10]; Camelina sativa [11]; soybean [12]; wheat [13]; Jatropha [14]; Arabidopsis [15]), animals [16,17], and insects [18]). Different oil crops were differentially expressed with different seed oil content and composition. There were 28 miRNAs regulated lipid metabolism from soybean [12], 30 miRNAs from Camelina sativa [11], and 13 miRNAs differentially expressed from two safflower genotypes that have difference to regulate oleic acid accumulation [10]. So there may be different miRNA expression between C. Meiocarpa and C. oleifera.
In order to explore the potential role of miRNAs in lipid metabolism during C. Meiocarpa and C. oleifera natural drying, miRNA expression profiles of seed samples at different moisture content (10,20,30,40, and 50%) were investigated using high throughput next generation small RNA sequencing technology, so the differentially expressed miRNAs are unraveled to help understand their involvement in lipid metabolism.

Plant material
In 2012, mature fruits of C. meiocarpa and C. oleifera were collected from the four cardinal directions of 3 superior trees' crowns per species. The trees are growing at the Minhou Tongkou State Forest Farm (26°05′ N, 119°17′ E), Fujian Province, China. The collected fruits were placed in a ventilated room until they naturally cracked and seeds were extracted by manual shell cutting. The seed moisture content at the time of extraction was closed to 50%. Seeds were naturally dried and their moisture content was determined daily. Over time, seed samples were collected at moisture content of 50,40,30,20, and 10% and were sequentially identified as S01 to S05 and S06 and S10 for C. meiocarpa and C. oleifera, respectively, and sampled three times, then were flash frozen in liquid nitrogen and stored at −80°C until RNA extraction.

Oil content analysis
To obtain the camellia seed oil content, the Soxhlet method was used by a fatty acid instrument (SZF-06A, Shanghai Hongji Instrument Equipment company) [28]. Camellia seeds at different moisture content were cut into thin slices, dried by silica gel and liquid nitrogen, milled into a powder by a pulverizer (FW100, Tianjin Taisite Instrument company). About 2 g of powder were weighed out (M1), packed in a folded filter paper bag and bound with a skim cotton thread, placed into the Soxhlet extraction thimble with 15 mL of petroleum ether for 1 h. Soxhlet extraction bottle were weight (M2) into which extraction thimble was placed, extracted for 5 h using 70 ml petroleum ether (75°C). After extraction, the solvent was evaporated in Soxhlet extraction bottle and dried at 75°C drying oven. Then Soxhlet extraction bottle was weighed again (M3).
All experiments were carried out at least in triplicate and data were analyzed using the SPSS statistics 17.0 software.

RNA isolation
Total RNA was extracted from the seed samples of the two camellia species using RNA kit (RNA simply total RNA kit, Tiangen, Beijing, China) according to the manufacturer's instructions. The purity and quality of the RNA were determined by assessing the absorbance ratio OD260/280 using NanoDrop ND1000 Spectrophotometer (NanoDrop, Wilmington, DE). The RNA was quantified with a Qubit 2.0 Fluoremeter (Invitrogen Corporation, Carlsbad, CA, USA) in accordance with the manufacturer's instructions. The integrity of the RNA samples was verified using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA).

Small RNA library construction and sequencing
Equal amount of total RNA from the three samples at each moisture content of the two camellia species, was mixed to construct a transcriptome library using an Illumina TruSeq RNA Sample PrepKit following the manufacturer's instructions. Small RNAs of 18-30 nt in length were separated and purified by denatureing polyacrylamide gel electrophoresis. After dephosphorylation and ligation of a pair of Solexa adaptors to their 5′ and 3′ ends, the products were reverse-transcribed and amplified by RT-PCR and gel purification. After the library was constructed, the Qubit 2.0 Fluorometer (Invitrogen Corporation, Carlsbad, CA, USA) were used to calculate the molar concentration and confirm the insert size. The cDNA libraries were sequenced using the Illumina HiSeq2500 Genome analyzer (Illumina Inc., San Diego, CA, USA).
The unannotated sequences were then analyzed by miRDeep-P software package to predict miRNAs, which was developed by modifying miRDeep2 [32]. All mature sequences, star sequences and precursor sequences cored by miRDeep2 were retained and regarded as putative miR-NAs and used for further analysis to identify known and novel miRNAs [33,34]. Known miRNAs were annotated by identifying their homologous imRNAs in miRBase database (http://www.mirbase.org/) using the following criteria: 1) seed region, nucleotides 2-7 must be identical; and 2) the remainder of the sequence alignment must contain no more than two mismatches [35,36]. The putative miRNAs produced by miRDeep2 analysis were regarded as conservative miRNAs when it met the above criteria. Those miR-NAs produced by miRDeep2 analysis that did not meet the above criteria were considered as novel miRNAs. In order to predict novel miRNAs with high confidence, only those with a miRDeep-P score higher than ≥0 were retained as true novel miRNAs [34,37].

Screening of differentially expressed miRNAs
Differentially expressed miRNAs were identified using the TPM [38] and IDEG6 software [39]. TPM (Tags Per Million reads) is a standardized method for calculating miRNA expression levels. TPM values were calculated using the following equation: TPM = number of mapped miRNA reads/number of clean sample reads × 10 6 . In order to calculate the levels of differential expressed miRNAs, normally the value was set as 0.01 by default when the sequencing read is 0 (no reads) [40]. We calibrated miRNA expression levels using multiple hypothesis tests with a false discovery rate (FDR) ≤0.01, performed generalized chi-square tests for differential miRNA expression using the IDEG6 software (http://telethon.bio.unipd.it/bioinfo/ IDEG6/), and screened the miRNAs for those with Pvalues less than 0.01 or TPM ratios between samples that were greater than 1 (fold change ≥1) or FDR ≤0.01. The miRNAs that met these criteria were identified as being differentially expressed [41].

miRNA target prediction
The putative target sites of miRNAs were identified by aligning mature miRNA sequences with a draft genome sequence using TargetFinder, (http://carringtonlab.org/resources/targetfinder). miRNA targets were computationally predicted essentially as described [38,42,43]. Briefly, potential targets from FASTA searches were scored using a position-dependent, mispair penalty system. Penalties were assessed for mismatches, bulges, and gaps (+1 per position) and G:U pairs (+0.5 per position). Penalties were doubled if the mismatch, bulge, gap, or G:U pair occurred at positions 2 to 13 relative to the 59 end of the miRNAs. Only one single-nucleotide bulge or singlenucleotide gap was allowed, and targets with penalty scores of four or less were considered to be putative miRNA targets [42,44].
Functional annotation of predicted target genes was assigned using Nr (non-redundant protein database, NCBI), Nt (non-redundant nucleotide database, NCBI), Swiss-Prot, GO (gene ontology, http://www.geneontology.org/) and COG (clusters of orthologous groups) databases. BLASTX was employed to identify related sequences in the protein databases based on an Evalue of less than 10 −5 [45]. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis were performed with package GO stats (http://www.geneontology.org/) of P value <0.05 was set as the cut-off to select out significantly enriched terms [46].

Quantitative real-time PCR analysis of miRNAs
qRT-PCR was used to validate the miRNAs identified using deep sequencing and to analyze their expression patterns. Total RNA of C. meiocarpa and C. oleifera seeds samples, was extracted using TRUzol Universal Reagent according to the manufacturer's protocol. They were then reverse-transcribed into cDNA using the microRNA cDNA First Strand cDNA Synthesis Kit (CWBIO, Beijing, China) according to the manufacturer's instructions. The cDNA was quantified by microRNA Real-Time PCR Assay Kit (CWBIO, Beijing, China) using a 20 μL reaction mixture, which consisted of 1 μL of diluted cDNA, 0.25 μM forward and reverse primer, and 10 μL of 2 × SYBR Green PCR Master Mix (CWBIO, Beijing, China). All reactions were performed under the following conditions: 95°C for 5 min, 40 cycles of 95°C for 15 S, 62°C for 45 S. Melting curve analysis was performed to verify specific amplification (from 72 to cycles at 60°C for 15 S). Each sample was processed in triplicate, and 5.8S rRNA was used as an internal control [47,48]. The equation ratio 2 -ΔΔCT was applied to calculate the relative expression level of miRNAs. The qRT-PCR primers are listed in file (Additional file 1: Table S1) and Ct value of 5.8S (Additional file 1: Table S2).

Results
Oil content of two camellia species during natural seed drying Camellia oil content were determined by Soxhlet extraction method during natural seed drying ( Table 1). The two species showed increased in oil content at 50 to 30% moisture contents followed by a slow decline at 30 to 10% moisture content, with 30% moisture content producing the highest oil content. Oil content accumulation was largest when the seed moisture content dropped from 50 to 40% and 40 to 30% for C. meiocarpa and C. oleifera, respectively. At 30% moisture content, C. meiocarpa and C. oleifera showed increase in their seed oil content of 8.80 and 10.23%, respectively, indicating that the effect of appropriate seed natural drying on oil accumulation can be great. While the percentages of oil content increase seem different between the two species, the relative increase amounted to 4.40%, indicating that natural drying can promote oil accumulation in both camellia specie.

Sequence analysis of small RNAs
To obtain a comprehensive profile of the sRNAs involved in natural drying, camellia seed from both species were subjected to Solexa high-throughput sequencing, with 5 libraries for each species corresponding to the sampled moisture contents. Average total of 30,641,435 and 31,012,304 reads were generated from C. meiocarpa and C. oleifera, respectively (Table 2). After filtering the low quality reads, such as 3′ insert null, poly(A), length < 18 nt or length > 30 nt, and other artifacts, the majority of the small RNAs were 21 to 24 nt in length. A total of 24,070,601 and 21,653,584 clean reads of 18-30 nt were obtained for C. meiocarpa and C. oleifera, respectively ( Table 2). GC content of clean reads were more than 47.75 and 48.99 and the Q30 (meaning 1 error in 1000 reads) of clean reads were more than 85.29 and 85.25% for C. meiocarpa and C. oleifera, respectively (Table 2). Through quality control, each sample of clean data were greater than 19.40 M, indicating the high sRNA quality (Tables 2, 3
The majority of small RNAs were 21 to 24 nt in length, producing similar length distributions in both species (Fig. 1). The 24 nt small RNAs were the most abundant representing 35.07 and 27.85% of small RNAs in C. meiocarpa and C. oleifera, respectively, the second was 21 nt representing 17.72 and 18.85%, third was 22 nt representing 16.35 and 17.97% (Additional file 1: Tables S3 and S4). The 40% moisture level (sample S02) for C. meiocarpa, produced the highest 24-nt RNA peak while 10% moisture level (sample S05) produced the highest 21-nt and 22-nt RNA peaks among the studied moisture content levels (Additional file 1: Table S3), conversely, 50% moisture (sample S06) produced the highest 24-nt, 21-nt, and 22-nt RNA peaks in C. oleifera (Additional file 1: Table S4).

Identification of miRNAs during natural drying
We used the software miRDeep2 to map the retained sequence reads to identify candidate miRNAs. Across the five moisture content levels, C. meiocarpa produced successful 2,355,539 reads (11.51%) that were mapped to the plant genomes, of which the 10 and 50% moisture content levels (S05 and S01) produced the most (2,724,598) and least (1,902,986) mapped reads, respectively. Similarly, on average, C. oleifera produced 2,396,805 (14.37%) successful mapped reads, of which 50 and 20% moisture content levels (S06 and S09) produced most (2,803,519) and least (2,121,984) mapped reads. In total, 2,288,508 (12.80%) mapped reads were successfully annotated (Table 5).
A total of 274 candidate miRNAs, 248 and 246 unique sequences were assigned to C. meiocarpa and C. oleifera, respectively (Additional file 1: Table S5). Out of the identified candidate miRNA sequences, 110 were identified to 64 families, 57 families belonging to each of C. meiocarpa (99 out of 248 candidate miRNA) and C. oleifera (98 out of 246 candidate miRNA) (Additional file 1: Tables S5, S6 and S7). The diversity of miRNA families could be determined from their members' number. As shown, MIR482 families were the largest with 10 members, followed by MIR159 and MIR535 with 5 members, and MIR160, MIR169_2 and MIR5272 with 4 members (Additional file 1: Table S6). Most of the conserved miRNA families had one member (65.63%) (Additional file 1: Table  S6). The miRNAs sequences ranged in length from 18 to 25 nt, with a peak of 24 nt (Additional file 1: Figure S1). The miRNA first nucleotide preference distributions are show in Fig. 2a. miRNAs of 24 nt tended to start with 5′-A while the 21 nt tended with 5′-U (Fig. 2a). Tall miRNAs tended to start with 5′-U and not 5′-G (Fig. 2b). During the seed natural drying process, the number of miRNAs across moisture content levels showed a decreasing trend which was followed by increase with a peak at 40% moisture content (S02) for C. meiocarpa (Table 6). C. oleifera, on the other hand, showed a trend of reduction, followed by increase, followed by reduction in the number of miR-NAs across the studied moisture content levels with a pronounced peak at 50% moisture content (S06) ( Table 6). Next we conducted sequence homology search of these candidate miRNAs against known mature miRNAs in miR-Base. Any miRNA that met the sequence homology criteria of Yang et al. [35] and Jain et al. [36] was considered a conserved miRNA gene. Through this analysis, we identified a total of 151 conserved miRNAs which belonging to 47 miRNA families in both camellia species (Additional file 1: Table S6). miRDeep2 score above 1.0 resulted in 98 pre-miRNAs (64.90%) of which 35 with a read count ranging between 10 and 100 (23.19%), 41 with read count above 100 (27.15%), and 75 with read count below 10 (49.67%) (Additional file 1: Table S6).
Those miRNA sequences, which met the threshold of miRDeep2 analysis but did not have any known homologous miRNA gene families in miRBase, were further analyzed to predict novel miRNAs in the two camellia species. The remaining miRNAs, which met the total score of >0, were considered to be true novel miRNAs. A total 123 novel mature miRNAs sequences were discovered and belonged to 36 miRNA families in both camellia species (Additional file 1: Table S7). miRDeep2 score above 1.0 had 87 pre-miRNAs (70.73%). The majority of pre-miRNA (69) read count ranged from 10 to 100 (56.10%), followed by 35 precursors in above 100 (28.45%) and 19 precursors below 10 (15.45%) (Additional file 1: Table S7).

Prediction and annotation of miRNAs target genes
To better understand the functions of the identified miRNAs, putative target genes were predicted using TargetFinder software [42]. A total of 6250 target genes were identified (   respectively (Table 7). A total of 2743 (73.48%) annotated target genes had the length of ≥1000 (Table 7).
To evaluate the potential functions of these miRNA target genes, we next applied gene ontology (GO) and KEGG pathway analyses to categorize the miRNA targets. The miRNA target genes were categorized according to biological process, cellular component and molecular function by Go analysis (Fig. 3). A total of 1860 target genes were categorized into 19 biological process. Based on molecular function, 1722 target genes were classified to 14 categories. A total of 1212 target genes were categorized as cellular components. Target genes related to lipid metabolism were found among 51 GO terms, in which 18 GO terms are related to biological process and 33 related to molecular function. There were 31 miRNAs involved in lipid metabolism and targeted 148 unigenes (Additional file 1: Table S8). The KEGG enrichment analysis revealed 12 pathways related to lipid metabolism, involved in 15 miRNA targeted 93 unigenes. There were 19 target genes in Glycolysis/Gluconeogenesis pathway, 12 in Fatty acid biosynthesis pathway, and 7 in biosynthesis of unsaturated fatty acids pathway (Additional file 1: Table S9). Integrated GO and KEGG function annotation, identified 23 miRNA regulating 131 target genes that were annotated as lipid metabolism (Additional file 1: Table S10). These miRNAs regulated the changes of seed oil content at different moisture content levels. There were high correlations between transcript abundance with added value of oil content, for example MIR482 (Additional file 1: Table  S11). Finally, miRNA of lipid metabolism only expressed were identified in C. Meiocarpa (Group1_Unigene_BMK. 30485_635795) and C. oleifera (Group1_Unigene_BMK. 37364_696840 and Group1_Unigene_BMK.38037_703962) ( Table 8).
The highest number of significantly differentially expressed miRNAs was observed between 40 and 50% (S06 vs. S07) moisture contents, with the highest up-and lowest downregulated numbers of 61 and 10, respectively (Table 9), with 8 significantly different miRNAs involved in lipid metabolism (Tables 11 and S13).
Comparing across the two species, the average number of miRNAs in C. oleifera seeds was higher than that of C. Meiocarpa (Table 9). Pairwise analysis of miRNA abundance between the two species for the same moisture level libraries indicated that there were 78, 51, 44, 58, and 61 significant differentially expressed miRNAs for the 50, 40, 30, 20, and 10% moisture contents, respectively. There were three differentially expressed miRNAs of lipid metabolism during the seed natural drying process of the studied two camellia species (Group1_Unigene_BMK.37987_703484, Group2_Unigene_BMK.25259_1025465, and Group2_Uni gene_BMK.25259_1025498) (Tables 12 and S14). The highest up-(69) and lowest down-regulated number (9) of significant differentially expressed miRNAs were detected for 50% moisture contents (S06 vs. S01) ( Table 9). This indicated that the greatest difference between the two species was observed at the 50% moisture content, with 12 significant differentially expressed miRNAs regulating lipid metabolism during seed natural drying (Tables 12 and S14).

Validation of the expression patterns of differentially expressed miRNAs related lipid metabolism by RT-qPCR
To validate the data obtained from the high-throughput sequencing, four significantly differentially expressed miRNAs  ne_BMK.63506_1315063, Group1_Unigene_BMK.37987_ 703484, and Group2_Unigene_BMK.38504_1137258) were predicted to target genes involved in lipid metabolism and their expression levels were quantified using stem-loop qRT-PCR (Fig. 4). The results were consistent with deep sequencing data (Table 8) and the majority of analyzed miRNAs showed moisture content-and species-specific expression. For C. meiocarpa, Group1_Unigene_BMK. 45675_802511 and Group2_Unigene_BMK.63506_1315063 peaked at 10 and 30% moisture content while the other two miRNAs (Group1_Unigene_BMK.37987_703484, and Group2_Unigene_BMK.38504_1137258) peaked at 20% moisture content. Additionally, all four miRNAs had different lowest point (Fig. 4). C. oleifera showed four miRNAs peaked at 50% moisture content but had different lowest point at 10 (Group1_Unigene_BMK.45675_802511 and Group2_Unigene_BMK.63506_1315063), 20 (Group2_ Unigene_BMK.38504_1137258), and 40% (Group1_Uni gene_BMK.37987_703484) moisture content (Fig. 4). These expression patterns indicate that lipid metabolism of the two camellia species were regulated by miRNA during the seed natural drying process.

miRNAs regulate lipid metabolism during camellia seed natural drying
Several studies have shown that miRNAs are differentially regulated in response to stress [56] and that this differential regulation varied among different plant species [57]. For example, distinct responses to drought stresses were reported for miRNAs in Arabidopsis [58], switchgrass [59], Populus [60] and Caragana intermedia [61]. Especially, drought stress in switchgrass [59] and Populus [60] were regulated by miRNAs related to lipid metabolism; however, the linkage between drought responses and lipid metabolism miRNAs changes is not established.
In the present study, C. Meiocarpa produced 9 significant differentially expressed miRNAs regulating lipid metabolism with only 4 at 20-30% moisture contents. These pre-miRNAs belonged to Group2_Unigene_BMK. 34335_1,093,229, Group1_Unigene_BMK.23434_588836, CL19777_Contig1_314088, and CL2440_Contig1_359627, with the former three showing more than 100 TPM (Table 10 and Additional file 1: Table S12) and targeting 3-ketoacyl-CoA synthase III, which catalyze the initial elongation step of fatty acid biosynthetic process [62] and glycerol-3-phosphate transporter, a precursor protein for phospholipid biosynthesis [63]. It is interesting to note that these three pre-miRNAs were down-regulated resulting in a reduction in the fatty acid biosynthesis, so seeds with 30% moisture content have high fatty acid synthesis and accumulation activities (Table 1).
Similarly, C. oleifera produced 8 significant differentially expressed miRNAs during seed natural drying, with upand down-regulated for the 40-50 and 30-40% moisture contents, respectively (Tables 11 and S13). The largest fold changes were observed for Group1_Unigene_BMK.23434_ 588836 and Group2_Unigene_BMK.34335_1,093,229, which target 3-ketoacyl-CoA synthase III (Additional file 1: Table S10 and Table S13). These pre-miRNAs regulate fatty acid biosynthesis with seeds at 40-50% moisture content showing low fatty acid biosynthesis activities while those at 30-40% moisture content exhibiting high activities (Table S13). This can be confirmed by the observed 1.06 and 9.17% increase in oil content at 40-50 and 30-40% moisture contents, and reaching the highest point at 30% moisture content (Table 1). So seeds with 30% moisture   content have also high fatty acid synthesis and accumulation activities.
Comparing the significant differentially expressed miR-NAs of C. meiocarpa (9) with those of C. oleifera (8), indicated that the two species share 6 miRNAs involved in lipid metabolism, with unique miRNAs belonging to each species (CL19455Contig1_54014 and Group2_ Unigene_BMK.38504_1137258 in C. oleifera, and CL2440 Contig1_359627, CL19777Contig1_314088 and Group1_ Unigene_BMK.45675_802511 in C. meiocarpa) (Tables 10   Table 12 Significant differentially expressed miRNAs of lipid metabolism during the two camellia species seed natural drying across the same moisture content level (C. oleifera vs. C.  Table S10). These indicate that oil content differences between the two camellia species are mainly due to their differential abilities miRNAs of fatty acid accumulation and catabolism. Collectively, the studied two camellia species produced 12 significant differentially expressed miRNAs to regulate lipid metabolism during seed natural drying (Tables 12 and S14). These pre-miRNAs indicated that C. meiocarpa has higher activities to regulate the lipid metabolism and this can be confirmed by its lower oil content as compared to C. oleifera (Table 1). It should be stated that all these 12 differentially expressed miRNAs were significantly expressed in the 50% moisture content stage (Tables 12  and S14). Seeds with 50% moisture content had only significant differentially expressed pre-miRNAs (CL19777Con tig1_314088, CL19455Contig1_54014, Group2_Unigene_ BMK.9543_1378570, CL9644Contig1_380257, Group2_ Unigene_BMK.38504_1137258 and Group2_Unigen e_BMK.38504_1137263) (Tables 8 and S14). CL19777 Contig1_314088 and CL19455Contig1_54014 target glycerol-3-phosphate transporter (Additional file 1: Table  S10). Group2_Unigene_BMK.9543_1378570 target acetyl-CoA-carboxylase (Additional file 1: Table S10), which catalyze the ATP-dependent carboxylation of acetyl-CoA to form malonyl-CoA, the rate limiting and first committed reaction in fatty acid synthesis [64]. So these three pre-miRNA control different fatty acid synthesis and accumulation in the two camellia species. For all significant differentially expressed miRNAs, the largest fold change was observed for Group2_Unigene_BMK.38504_1137263 (Fatty acid hydroxylase), followed by CL9644Contig1_ 380257 (Carboxylesterase), and Group1_Unigene_BMK. 23434_588836 (3-ketoacyl-CoA synthase III) and Group2_ Unigene_BMK.34335_1,093,229 (3-ketoacyl-CoA synthase III) (Additional file 1: Table S10 and Table S14). So the oil content of C. oleifera is higher than C. meiocarpa and this is attributable to four pre-miRNAs, of which Group2_ Unigene_BMK.38504_1137263 and CL9644Contig1_380257 regulating fatty acid catabolism and Group1_Unigene_ BMK.23434_588836 and Group2_Unigene_BMK.34335_ 1,093,229 regulating fatty acid synthesis.
To our knowledge, this work provides the first small RNA expression analysis of lipid metabolism in camellia seed during natural drying as well as the first comparative miRNA profiling analysis between C. Meiocarpa and C. oleifera that exhibit significantly different fatty acid profiles.

Additional file
Additional file 1: Figure S1. Length distributions of miRNAs in two camellia species. Table S1. qRT-PCR primer sequences. Table S2. Length distributions of small RNAs in C. meiocarpa during seed natural drying. Table S3. Length distributions of small RNAs in C. oleifera during seed natural drying. Table S4. miRNAs transcript abundance in the two camellia species during seed natural drying. Table S5. All conservative miRNA discovered in the two camellia species during seed natural drying. Table S6. All novel miRNA discovered in the two camellia species during seed natural drying. Table S7. GO terms related with the lipid metabolism in the two camellia species during seed natural drying. Table S8. KEGG enrichment related with the lipid metabolism in the two camellia species during seed natural drying. Table S9. miRNA of lipid metabolism targets and their putative functions. Table S10. Differentially expressed miRNAs of lipid metabolism during C. meiocarpa seed natural drying. Table S11. Differentially expressed miRNAs of lipid metabolism during C. oleifera seed natural drying. Table S12. Differentially expressed miRNAs of lipid metabolism between two camellia species during seed natural drying. (DOCX 168 kb)