Abstract

   Tea oil camellia (Camellia oleifera), an important woody oil tree, is a
   source of seed oil of high nutritional and medicinal value that is
   widely planted in southern China. However, there is no report on the
   identification of the miRNAs involved in lipid metabolism and seed
   development in the high- and low-oil cultivars of tea oil camellia.
   Thus, we explored the roles of miRNAs in the key periods of oil
   formation and accumulation in the seeds of tea oil camellia and
   identified miRNA–mRNA regulatory modules involved in lipid metabolism
   and seed development. Sixteen small RNA libraries for four development
   stages of seed oil biosynthesis in high- and low-oil cultivars were
   constructed. A total of 196 miRNAs, including 156 known miRNAs from 35
   families, and 40 novel miRNAs were identified, and 55 significantly
   differentially expressed miRNAs were found, which included 34
   upregulated miRNAs, and 21 downregulated miRNAs. An integrated analysis
   of the miRNA and mRNA transcriptome sequence data revealed that 10
   miRNA–mRNA regulatory modules were related to lipid metabolism; for
   example, the regulatory modules of ath-miR858b–MYB82/MYB3/MYB44
   repressed seed oil biosynthesis, and a regulation module of
   csi-miR166e-5p–S-ACP-DES6 was involved in the formation and
   accumulation of oleic acid. A total of 23 miRNA–mRNA regulatory modules
   were involved in the regulation of the seed size, such as the
   regulatory module of hpe-miR162a_L-2–ARF19, involved in early seed
   development. A total of 12 miRNA–mRNA regulatory modules regulating
   growth and development were identified, such as the regulatory modules
   of han-miR156a_L+1–SPL4/SBP2, promoting early seed development. The
   expression changes of six miRNAs and their target genes were validated
   using quantitative real-time PCR, and the targeting relationship of the
   cpa-miR393_R-1–AFB2 regulatory module was verified by luciferase
   assays. These data provide important theoretical values and a
   scientific basis for the genetic improvement of new cultivars of tea
   oil camellia in the future.

   Keywords: Camellia oleifera, small RNA sequencing, miRNA–mRNA
   regulatory modules, lipid metabolism, seed development

1. Introduction

   Tea oil camellia (Camellia oleifera Abel.) is an important woody oil
   tree in southern China and has a cultivation history of over two
   thousand years. Together with oil palm, olive, and coconut, tea oil
   camellia is well known as one of the four major woody oil trees in the
   world [[38]1,[39]2]. Tea oil camellia seeds have a high oil content of
   up to 40%; are rich in unsaturated fatty acids (the oleic acid content
   accounts for over 75% of all the fatty acids in the oil of the seeds);
   exhibit an abundance of bioactive compounds, such as vitamins,
   camelliaside, and tea polyphenols; and are an important source of
   high-grade edible vegetable oil [[40]3]. The oil extracted from tea oil
   camellia seeds is not easily oxidized, as it has high chemical
   stability and can be preserved for a long time; it can significantly
   contribute to the prevention of cardiovascular diseases and cancer and
   reduce cholesterol and blood lipids [[41]4,[42]5]. However, the seed
   yield of tea oil camellia is not stable, and the seed oil content
   varies greatly among different germplasms, which seriously restricts
   the sustainable development of the tea oil camellia industry. One
   important strategy is the breeding of tea oil camellia cultivars with
   high seed oil contents and high and stable seed yields.

   MicroRNAs (miRNAs) are a class of endogenous noncoding small RNAs with
   lengths of approximately 20–25 nucleotides, and they play a crucial
   role in mediating the posttranscriptional regulation of gene expression
   by targeting mRNAs for degradation or inhibiting their translation in
   eukaryotes. MiRNAs are derived from primary miRNA transcripts
   (pri-miRNAs) that contain a stem–loop secondary structure; the
   pri-miRNA is processed in the nucleus by DCL1, a Dicer-like protein, to
   create an miRNA–miRNA* duplex, where miRNA* is a passenger strand
   complementary to the miRNA. The duplex is then separated by helicase,
   and the mature miRNA is incorporated into an ARGONAUTE 1 (AGO1) protein
   to form an RNA-induced silencing complex (RISC) [[43]6,[44]7]. Target
   genes that contain a sequence with almost complete complementarity to
   the miRNA are cleaved by the RISC at a specific site opposite to the
   10th or 11th nucleotide in the miRNA [[45]8]. Studies have confirmed
   the key roles of miRNAs in various biological and metabolic processes
   in plants, such as fatty acid biosynthesis [[46]9,[47]10], lipid
   metabolism [[48]11,[49]12], growth and development [[50]13], responses
   to stresses [[51]14], and cellular proliferation and differentiation
   [[52]15]. As of 3 December 2021, a total of 38,589 miRNA sequences have
   been deposited in the database, miRBase [[53]16,[54]17], and many
   conserved and novel miRNAs have been identified in soybeans [[55]18],
   Arabidopsis [[56]19], Brassica napus [[57]20], rice [[58]21], peonies
   [[59]22], and maize [[60]23]. In walnuts, jre-miRn105, jre-miRn434,
   jre-miR477d, and jre-miR156a.2 are key miRNAs that regulate FA
   synthesis, and jre-miRn411 and jre-miR399a.1 are closely related to oil
   accumulation [[61]10]. MiRNAs control and influence a variety of
   physiological processes by regulating target genes. For example, the
   KAS and KAR genes, targeted by miR159 and miR156, respectively, are
   important for lipid biosynthesis [[62]24]. The FAD2 gene, targeted by
   miR159b and miR5026, regulates and influences FA biosynthesis
   [[63]18,[64]25]. Enhancing microRNA167A expression in the seeds
   decreased the alpha-linolenic-acid content and increased the seed size
   in Camelina sativa [[65]26]. MiR172 plays a highly important role in
   seed development by modulating its target ARF2 gene [[66]11].

   Transcription factors play critical roles in regulating lipid
   biosynthesis [[67]27], such as WRINKLED1 and LEAFY COTYLEDON, MYB, SPL,
   ARF, and AP2 [[68]28], which have been identified as targets of the
   major miRNAs [[69]29]. On the basis of the high-throughput small RNA
   and degradome analyses of soybean seeds 15 days after flowering, 55
   annotated and 26 novel miRNAs, which targeted 145 genes, were
   identified [[70]18]; a total of 82% of the targeted genes were
   transcription factors, including the ARF, MYB, TCP, GRF, and NAC
   families [[71]18]. Known transcription-factor-encoding genes involved
   in seed size/weight determination, including SPL, GRF, MYB, ARF,
   HAIKU1, SHB1, KLUH/CYP78A5, and E2Fb, along with novel genes, were
   found to be targeted by the predicted miRNAs in chickpeas [[72]30].
   MiRNA–target expression profiles evidenced that some miRNAs could tune
   distinct seed developmental stages by targeting the HD-ZIP, ARF, SPL,
   and NF-Y transcription-factor families in Phaseolus vulgaris [[73]31].
   MiR160 negatively regulates ARFs that significantly affect seed
   development in A. thaliana [[74]32]. SPL10 and SPL11 are targeted by
   miR156 and are involved in the early morphogenesis of Arabidopsis
   embryos [[75]33]. MYBs targeted by miR159 are involved in seed size,
   and a double mutation of miR159 (miR159ab) increased MYB33 and MYB65
   expression, which promoted the formation of small seeds [[76]34]. The
   transcription factor, WRI1, targeted by hrh-miRn215 in sea buckthorn
   seeds, involved in oil biosynthesis and fatty acid formation and
   accumulation [[77]35], ethylene-responsive transcription factor 3
   (ERF3), MADS domain protein, AGAMOUS-LIKE 61 (AGL61), and WRKY
   transcription factor 41 (WRKY41) were targeted by miR171k-5p_1,
   miR7760-p3_1, and Xso-miRn80 in the seeds of Xanthoceras sorbifolium
   [[78]11].

   Currently, miRNA–mRNA regulatory modules controlling multiple
   biological and metabolic processes have been revealed via the
   functional analysis of woody oil plants, such as hazel [[79]36], oil
   palm [[80]37], olive [[81]38], peony [[82]22], walnut [[83]10], and
   yellowhorn [[84]11]. In developing sea buckthorn seeds, 19 (4 known and
   15 novel) and 22 (6 known and 16 novel) miRNAs were found to be
   involved in lipid biosynthesis and seed size, respectively [[85]39],
   and an integrated analysis of the mRNA and miRNA transcriptomes and the
   qRT-PCR revealed key miRNAs and their targets (miR164d-ARF2,
   novelmiRNA-108-ACC, novelmiRNA-23-GPD1, novelmiRNA-58-DGAT1, and
   novelmiRNA-191-DGAT2) that are potentially involved in the seed size
   and lipid biosynthesis in the seeds of sea buckthorn. In developing
   seeds of yellowhorn, 141 differentially expressed miRNAs (120 known and
   21 novel miRNAs) were identified, along with several miRNA–target
   regulatory modules involved in the lipid biosynthesis of yellowhorn
   seeds, such as miR319p_1-FAD2–2, miR5647-p3_1-DGAT1, and
   miR7760-p5_1-MED15A [[86]11]. However, there are no reports on the
   roles of miRNAs in the different developmental stages of the seeds of
   high- and low-oil cultivars of tea oil camellia.

   To understand the potential roles of miRNAs and their target genes in
   the developing seeds of tea oil camellia, high-throughput Illumina
   sequencing technology was used to reveal the global profile of small
   RNAs in the developing seeds of high- and low-oil tea oil camellia
   cultivars. Then, we identified some novel and known miRNAs and their
   target genes by the integrated analysis of the mRNA and miRNA
   transcriptomes. Finally, we further analyzed differentially expressed
   miRNAs and the GO and KEGG enrichments of these miRNAs target genes. In
   addition, reverse transcription–quantitative PCR (qRT-qPCR) was used to
   validate the sequencing results for a selected group of miRNAs and the
   expression profiles of their targets, and the targeting relationship of
   the miRNA–mRNA regulatory module was verified by luciferase assays. Our
   datasets and results will promote the development of tea oil camellia
   miRNAs in public resource databases, and will provide a scientific
   resource and basis for further investigations of the identified miRNAs
   of tea oil camellia and their regulatory roles in the genetic
   improvement of high-quality and high-yielding tea oil camellia in the
   future.

2. Materials and Methods

2.1. Study Site and Sampling

   The sample trees were from two cultivars of tea oil camellia, the
   cultivar, ‘M3’, with a high seed oil content, and the cultivar, ‘M8’,
   with a low seed oil content (voucher No. 1303 and No. 1308,
   respectively, identified by Chengjiang Ruan and deposited at Dalian
   Minzu University) growing in the tea oil camellia orchard in Yuping
   Dong Autonomics County, Guizhou Province, China. It is located at
   27°28′–27°31′ north latitude, and 108°34′–109°09′ east longitude. The
   orchard had a mean annual rainfall of 1174.1 mm, a mean annual
   temperature of 16.4 °C, and a mean annual sunshine rate of 1227.8.
   These two cultivars were selected from the tea oil camellia orchard
   because they had very similar phenotypic characteristics and genetic
   backgrounds, but not seed oil content.

   For tea oil camellia, the key period of seed oil biosynthesis and
   accumulation is from June to September [[87]40]. The oil contents in
   the seeds of these two cultivars showed a gradually upward trend from
   early June to early July; then, they began to rise rapidly from early
   July to early September and, subsequently, remained stable. The fruits
   were collected at random from the different plants of both cultivars on
   2 June, 4 July, 5 August, and 3 September 2016, which were the same
   dates as those for the published mRNA-seq data [[88]41], and the four
   sampling periods were referred to as T1~T4, for each cultivar in the
   data analysis; thus, we recorded them as M3T1, M3T2, M3T3, and M3T4,
   and as M8T1, M8T2, M8T3, and M8T4 for the samples of the two cultivars
   of the four periods, respectively. The seeds obtained by stripping the
   fruit shells were immediately wrapped in tin foil and placed into
   liquid nitrogen. The frozen seeds were transported to the laboratory
   and stored at −80 °C for the following experiments. The collection of
   all the samples fully complied with local and national legislation.

2.2. Small RNA Library Construction and Sequencing

   The total RNA of each sample was extracted by using TRIzol reagent
   (Invitrogen, Carlsbad, CA, USA), according to the manufacturer’s
   instructions. The quality and quantity of the total RNA of each sample
   were analyzed using the Bioanalyzer 2100 (Agilent, Santa Clara, CA,
   USA) and RNA 6000 Nano LabChip Kit (Agilent, Santa Clara, CA, USA),
   ensuring that the RIN (RNA integrity number value) was >7.0. The small
   RNA library for each sample was constructed with the TruSeq Small RNA
   Sample Prep Kit (Illumina, San Diego, CA, USA), in accordance with the
   manufacturer’s protocol. Then, we performed single-end sequencing (36
   bp or 50 bp) on an Illumina HiSeq2500 at LC-BIO (Hangzhou, China),
   following the vendor’s recommended protocol.

2.3. Identification and Analysis of miRNA Sequencing Data

   To obtain accurate and reliable miRNA sequencing data, the raw read
   quality was evaluated using Illumina FastQC to obtain Q30 data. Then,
   the raw reads were processed to obtain valid reads with an in-house
   program, ACGT101-miR (LC Sciences, Houston, TX, USA) [[89]42], to
   filter out adapter-dimer reads, junk reads, low-quality reads,
   low-complexity reads, reads for common RNA families (rRNA, tRNA, snRNA,
   and snoRNA), and repeats. Then, the unique sequences, with lengths of
   18–25 nt, were selected to map to the precursors of 66 specific species
   ([90]Table S1) in miRBase 21.0 ([91]http://www.mirbase.org/, accessed
   on 3 August 2017)), through a BLAST search, to identify known miRNAs
   and novel 3p- and 5p-derived miRNAs [[92]43]. Because miRNAs are
   conserved, related species of tea oil camellia were selected as
   specific species. In the alignment analysis, the seed length was set to
   16, and the length variation at both the 3′ and 5′ ends, and one
   mismatch inside the sequence, were allowed. The unique sequences
   mapping to the hairpin arms of the mature miRNAs of specific species
   were identified as known miRNAs, and the unique sequences mapping to
   the other arms of known specific species’ precursor hairpins opposite
   to the annotated mature miRNA-containing arms were considered to be
   novel 5p- or 3p-derived miRNA candidates [[93]44]. The remaining
   sequences were mapped to other selected species precursors (with the
   exclusion of specific species) in miRBase 21.0 by a BLAST search, and
   the mapped pre-miRNAs were further BLAST-searched against the specific
   species genomes to determine their genomic locations. We defined the
   above two as known miRNAs. The unmapped sequences were BLAST-searched
   against specific genomes, and the hairpin RNA structures containing
   sequences were predicted from the flanking 120 nt sequences using the
   RNAfold software ([94]http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi,
   accessed on 13 October 2017)). The criteria for the secondary structure
   prediction were: (1) The number of nucleotides in one bulge in the stem
   (≤12); (2) The number of base pairs in the stem region of the predicted
   hairpin (≥16); (3) The cutoff of free energy (kCal/mol ≤−15); (4) The
   length of the hairpin (up and down stems + terminal loop ≥50); (5) The
   length of the hairpin loop (≤200); (6) The number of nucleotides in one
   bulge in the mature region (≤4); (7) The number of bias errors in one
   bulge in the mature region (≤2); (8) The number of biased bulges in the
   mature region (≤2); (9) The number of errors in the mature region (≤4);
   (10) The number of base pairs in the mature region of the predicted
   hairpin (≥12); and (11) The percent of mature sequence in the stem
   (≥80). After alignment analysis, the miRNAs that met the above criteria
   were considered novel miRNAs.

2.4. Analysis of Differentially Expressed miRNAs

   The number of clean reads originating from each miRNA represents the
   relative expression abundance, or the level of the corresponding miRNA
   in small RNA deep sequencing. Normalized read counts were calculated
   for miRNA differential expression analysis according to the procedure
   in a previous study [[95]45], with minor modifications: (1) Find a set
   of sequences common among all the samples; (2) Construct a reference
   dataset. Each datum in the reference set is the copy number median
   value of a corresponding common sequence of all the samples; (3)
   Perform base-2 logarithm transformation on the copy numbers
   (log[2](copy#)) of all the samples and reference the dataset; (4)
   Calculate the (log[2](copy#)) difference (∆log[2](copy#)) between the
   individual sample and the reference dataset; (5) Form a subset of
   sequences by selecting |∆log[2](copy#)| < 2, which means less than a
   (2^2 =) four-fold change from the reference set; (6) Perform linear
   regressions between individual samples and the reference set on the
   subset sequences to derive linear equations, y = a[i]x + b[i], where
   a[i] and b[i] are the slope and intercept, respectively, of the derived
   line; x is the log[2](copy#) of the reference set; and y is the
   expected log[2](copy#) of sample i on a corresponding sequence; (7)
   Calculate the mid value, x[mid] = (max (x) − min (x))/2 of the
   reference set. Calculate the corresponding expected log[2](copy#) of
   sample i, y[i,mid] = a[i]x[mid] + b[i]. Let y[r,mid] = x[mid], and let
   ∆y[i] = y[r,mid] − y[i,mid], which is the logarithmic correction factor
   for sample i. We then derive the arithmetic correction factor, fi =
   2^Δyi, for sample i; (8) Correct the copy numbers for the individual
   samples by multiplying the corresponding arithmetic correction factor,
   fi, by the original copy numbers.

   A Student’s t-test was used to identify miRNAs differentially expressed
   between two samples on the basis of the normalized read counts. An
   analysis of variance (ANOVA) was used to identify miRNAs differentially
   expressed among samples at the four development stages in each line by
   using the mirnaTA software [[96]46]. The miRNAs with p-values < 0.05
   were considered differentially expressed miRNAs [[97]45]. The
   normalized read counts of some miRNAs were set to be 0.01 for further
   calculation if they had no reads in the library.

2.5. Prediction and Identification of miRNA–mRNA Regulatory Modules

   The target genes of differentially expressed miRNAs were identified by
   aligning mature miRNA sequences with published mRNA-seq sequences
   [[98]41] using TargetFinder
   ([99]https://github.com/carringtonlab/TargetFinder, accessed on 8 April
   2018) [[100]47], following the procedures described in previous studies
   [[101]48]. The TargetFinder algorithm is based on the miRNA and mRNA
   complementary-pairing principle. The alignment score value in the
   prediction result was used as the screening standard, which was the
   score of the prediction result. The comparison between the miRNA and
   mRNA was performed, and then the complementary pairing score was given
   for the comparison position: a total of 1 point was given for one
   mismatch or deletion, and 0.5 points were given for one G:U pairing;
   the penalty for a mismatched G:U pairing in the core region was
   doubled, and the core region began at the 2nd nt and continued to the
   13th nt of the miRNA (5′–3′). Finally, the results were obtained for
   scores less than or equal to 4. The significance of the alignment score
   was the degree of match between the target gene and the miRNA. The
   lower the score, the more complete the matching, and the more credible
   it was.

   The GO terms and KEGG pathways of genes targeted by these
   differentially expressed miRNAs were annotated. The GO function and
   KEGG information of this species were used to perform GO function
   annotation ([102]http://bioinfo.cau.edu.cn/agriGO/, accessed on 1
   September 2018) and KEGG signaling pathway annotation
   ([103]https://www.kegg.jp/kegg/pathway.html, accessed on 1 September
   2018) for the target genes of differentially expressed miRNAs. Fisher’s
   exact hypothesis test was performed on the GO and KEGG pathways of the
   targeted genes, and the enrichment analysis of each GO term and KEGG
   pathway was separately performed. According to the prediction results
   for the targeted genes of the differentially expressed miRNAs,
   significance differences in the GO and KEGG information were defined
   with a default threshold, p < 0.05.

2.6. Identification of the Expression of miRNAs and Their Target Genes by
qRT-PCR

   To validate the high-throughput sequencing data, six pairs of miRNA
   target genes, on the basis of the sequencing data, were selected to
   perform the qRT-PCR. The total RNA from each sample for the qRT-PCR
   analysis was extracted, as described above, for the sequencing
   experiments. First-strand cDNA was synthesized using the Mir-X miRNAs
   First-Strand Synthesis Kit (TaKaRa Biotech, Dalian, China), according
   to the manufacturer’s instructions. The U6 snRNA was selected as the
   reference gene, and the qRT-PCR for the miRNA experiments was performed
   in an ABI7500 real-time PCR instrument (Applied Biosystems, Foster
   City, CA, USA), using the Mir-X miRNA qRT-PCR SYBR Kit (TaKaRa Biotech,
   Dalian, China) and miRNA-specific primers ([104]Table S2). The qRT-PCR
   experiments for each sample were performed with three replicates.

   The predicted target genes of the miRNAs were selected for qRT-PCR
   analyses. First-strand cDNA was synthesized using the PrimeScript^TM RT
   Master Mix (TaKaRa Biotech, Dalian, China), according to the
   manufacturer’s instructions. All the specific primer pairs of the
   predicted target genes ([105]Table S3) were designed using the Primer
   Quest online software
   ([106]http://sg.idtdna.com/PrimerQuest/Home/Index, accessed on 10
   September 2018). Elongation factor 1-alpha (CoEF1α) was selected as a
   reference gene. The qRT-PCR experiments were performed using the SYBR
   Premix Ex Taq^TM II Kit (Tli RNaseH Plus; TaKaRa Biotech, Dalian,
   China). The qRT-PCR was conducted using the ABI7500 real-time PCR
   instrument (Applied Biosystems, Foster, CA, USA). The qRT-PCR analyses
   for each sample were performed in triplicate. The relative expression
   levels of each gene were calculated using the 2^−ΔΔCt method [[107]49].

2.7. Dual-Luciferase Reporter Assay

   Many studies have shown the AFB2 gene and miR393 to be involved in
   plant growth and development [[108]50,[109]51], but few have reported
   on the potential role of the miR393–AFB2 regulatory module in seeds.
   Thus, we investigated the targeting relationship between cpa-miR393_R-1
   and AFB2 in the seeds of tea oil camellia. Fragments from the 3′-UTR of
   AFB2 containing the predicted binding sequences for cpa-miR393_R-1 were
   amplified and subcloned into the pmirGLO luciferase promoter vector.
   The pCDNA3.1 plasmid was used as the template vector. The fragment
   containing the nucleotide sequences of the precursor of cpa-miR393_R-1
   was cloned into the vector to construct the recombinant vector
   expressing cpa-miR393_R-1 pCDNA3.1, as described previously. The
   pmirGLO vector containing the 3′-UTR of AFB2 was cotransfected with
   pCDNA3.1 or pCDNA3.1, containing pre-cpa-miR393_R-1 into HEK-293T cells
   [[110]52], using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA),
   according to the manufacturer’s protocol and a previous report
   [[111]53,[112]54,[113]55]. Forty-eight hours after treatment, the
   firefly and Renilla luciferase activities were measured using a
   luciferase reporter assay kit (BioVision, Inc., Milpitas, CA, USA).
   Renilla was used as the transfection control.

3. Results

3.1. Brief Outline for Sequencing Data, from Raw Data to Cleaned Sequences

   To comprehensively identify key small RNAs in the developing seeds of
   high- and low-oil cultivars of tea oil camellia, 16 small RNA libraries
   were sequenced on an Illumina HiSeq2500 platform for four different
   seed developmental stages, with two biological replicates for each seed
   developmental stage of each cultivar. To evaluate the consistency of
   the biological replicates, the Pearson correlation between the samples
   was calculated and is shown in [114]Figure S1. A high correlation
   coefficient (r) between every two replicates confirmed the reliability,
   as well as the operational stability, of the experimental results.
   Totals of 13,572,633 and 14,183,198 raw reads were obtained from the
   cultivars, ‘M3’ and ‘M8’, of tea oil camellia, respectively.
   Furthermore, the raw reads were analyzed with the LC sRNA analysis
   pipeline, ACGT101-miR, and totals of 11,653,829 (85.67%) and 12,248,627
   (86.23%) valid reads were obtained for ‘M3’ and ‘M8’ after removing the
   adaptor dimers, junk reads, low-complexity sequences, and sequences
   shorter than 18 nucleotides, or longer than 25 nucleotides,
   respectively ([115]Tables S4 and S5). The length distributions of the
   unique sRNAs for both cultivars at four developmental stages were then
   summarized. Most of the sRNA reads ranged from 21 to 24 nt in length,
   leading to similar length distributions in two tea oil camellia
   cultivars at the different developmental stages. sRNAs of 24 nt were
   the most abundant, accounting for 71.58% (M8T1) to 79.74% (M3T4) of the
   total ([116]Figure 1). In addition, 21-, 22-, 23-, and 25-nt sRNAs were
   common, which were more abundant than those of any other length,
   besides 24 nt.

Figure 1.

   [117]Figure 1
   [118]Open in a new tab

   Length distributions of unique miRNAs from the seeds of the ‘M3’ and
   ‘M8’ cultivars of tea oil camellia at four different developmental
   stages.

3.2. Identification of Known and Novel miRNAs in Developing Seeds of Tea Oil
Camelia

   To detect the miRNAs involved in the development of tea oil camellia
   seeds, the unique sequences were selected and mapped against miRBase
   21.0 by a BLAST search to identify the known miRNAs and novel miRNAs,
   and then a total of 196 expressed miRNAs were identified in the
   developing seeds of tea oil camellia ([119]Table S6). The expressed
   miRNAs were classified into five groups (gp1, gp2a, gp2b, gp3, and
   gp4), according to criteria defined in a previous study [[120]44], with
   minor modifications: seven known miRNAs belonged to gp1; 149 conserved
   and known miRNAs belonged to gp2a, gp2b, and gp3; and 40 novel miRNAs
   belonged to gp4 ([121]Table S6). Analyses of the number and length
   ratios of unique miRNA sequences of different lengths revealed that the
   lengths of the miRNAs were mainly distributed within 20~24 nt
   ([122]Table S7). MiRNAs with a length of 21 nt were frequent, and
   accounted for 49.55% of all the miRNAs, which is consistent with the
   definition of miRNAs ([123]Table S7). The results of the statistical
   analysis of the expression levels of the detected miRNAs, and the
   detection rates for miRNAs upon evaluating redundant miRNAs, showed
   that 196 unique miRNAs were identified, including 156 known miRNAs
   ([124]Table S8), and 40 novel miRNAs ([125]Table S9).

   The miRNA sequences were mapped to known miRNAs from 42 plant species,
   with the highest number of miRNAs mapped to known gma-miRNAs of Glycine
   max (17 (10.89%)), followed by ath-miRNAs of Arabidopsis thaliana (13
   (8.33%)), and aly-miRNAs of Arabidopsis lyrata (11 (7.05%)) ([126]Table
   S8). Among the identified known miRNA sequences, 156 were identified as
   belonging to 35 miRNA families, and the MIR171 family had the largest
   number of members (12), followed by the MIR167 (11), and the MIR156 and
   MIR5645 families, with 10 members each ([127]Table S8). The analysis of
   the distribution of the miRNA first-nucleotide preferences showed that