Overview of the sequencing data from raw data to cleaned sequences
To comprehensively identify key small RNAs in the seeds of high- and low-oil cultivars of tea oil camellia at four different seed developmental stages, 16 small RNA libraries were sequenced on an Illumina HiSeq2500 platform, with two biological replicates for each seed developmental stage of each cultivar. A total of 13,572,633 and 14,183,198 raw reads were generated from the tea oil camellia cultivars ‘M3’ and ‘M8’, respectively. Further, the raw reads were analyzed with the LC sRNA analysis pipeline ACGT101-miR, and a total of 11,653,829 (85.67%) and 12,248,627 (86.23%) valid reads were obtained for ‘M3’ and ‘M8’, respectively (Table S3 and Table S4). Based on the analysis and statistics of the raw sequencing data, the analysis of the overall length distributions of the valid data after filtration (Fig. 1) showed that most of the valid data were 20-24 nt in length, which was consistent with the typical characteristics of dicer enzymatic cleavage.
Identification and analysis of known and novel miRNAs
To detect the miRNAs involved in the development of tea oil camellia seeds, the expressed miRNAs were classified into four groups (gp1, gp2a, gp2b and gp3) according to defined criteria (Table S5). The unique sequences were selected and mapped against miRBase 21.0 by BLAST search to identify the known miRNAs and novel miRNAs. A total of 196 expressed miRNAs were identified in the developing seeds of tea oil camellia, including seven known miRNAs belonging to the gp1 group, 149 conserved and known miRNAs belonging to the groups gp2a, gp2b and gp3 and 40 novel miRNAs belonging to gp4 (Table S5). Analyses of the number and length ratios of unique miRNA sequences of different lengths revealed that the lengths of miRNAs were mainly distributed within 20~24 nt (Table S6). The number of miRNAs with a length of 21 nt was frequent and accounted for 49.55% of all miRNAs, which is consistent with the definition of miRNAs (Table S8).
The results of the statistical analysis for the expression level of the detected miRNAs and the detection rates of miRNAs by evaluating redundant miRNAs showed that 196 unique miRNAs were identified, including 156 known miRNAs and 40 novel miRNAs (Table S7). Thirty-eight of all novel miRNAs had a medium expression level; in contrast, 38 known miRNAs had a high expression level (a high expression level indicates that the number of reads in the reported miRNAs was higher than the average copy number of the data set; a medium expression level indicates that the number of reads in the reported miRNAs was higher than 10 and less than the average copy number of the data set; and a low expression level indicates that the number of reads in the reported miRNAs was less than 10), and only 4 novel miRNAs had a low expression level, namely, PC-5p-236956_11, PC-3p-145806_25, PC-3p-253727_10 and PC-3p-150080_24. In addition, we detected the number of specifically expressed miRNAs and the number of specifically coexpressed miRNAs with a Venn diagram by comparing samples ‘M3’ (a) and ‘M8’ (b) at four different developmental stages (Fig. 2a, b). It was found that 54 and 63 coexpressed miRNAs were identified in the cultivars ‘M3’ and ‘M8’ at four different developmental stages, respectively. Mining these coexpressed miRNAs could provide a favorable foundation for the genetic improvement of tea oil camellia via breeding.
Identification of differentially expressed miRNAs
The differentially expressed miRNAs were identified by suitable differential test methods. Significantly differentially expressed miRNAs were defined with a significance threshold of P-value <0.05. There were 55 significantly differentially expressed miRNAs identified between the samples of one cultivar at two different developmental stages and between different samples of two cultivars at the same developmental stage. A total of 34 upregulated miRNAs (red) and 21 downregulated miRNAs (green) were identified among these groups (Fig. 3). Moreover, to identify the expression levels of significantly differentially expressed miRNAs at different seed developmental stages, the significantly differentially expressed miRNAs were compared between different groups (Table 1). Five miRNAs (mtr-MIR2586a-p3_1ss15TC, ath-MIR5645b-p5_2ss19TG21TC, ath-MIR5645b-p3_2ss19TG21TC, mtr-MIR2586a-p3_2ss4TG21CT, and hpe-miR166a_1ss5AC) were upregulated in ‘M3’ relative to ‘M8’ at the early seed developmental stage. At the second seed developmental stage, three miRNAs (csi-miR1515_1ss16AG, stu-miR8016_1ss7GT, and ppe-miR172a-5p_1ss21GA) were upregulated in ‘M3’ relative to ‘M8’. At the late seed developmental stage, four miRNAs (PC-3p-9080_426, mtr-MIR2586a-p3_1ss15TC, ghr-MIR7499-p3_2ss10GC17CG, and rco-miR167a_R+1) were downregulated in ‘M3’ compared with their expression in ‘M8’.
In addition, the number of common and specific differentially expressed miRNAs was visually shown with a Venn diagram comparison between different groups. The results showed that there were four significantly differentially coexpressed miRNAs (csi-miR1515_1ss16AG, ccl-miR167a, aly-miR167b-3p_1ss13AG, and ppe-miR172a-5p_1ss21GA) when the different stages in ‘M3’ were compared and M3T2 vs. M8T2 were compared and that there were six significantly differentially coexpressed miRNAs (csi-miR1515_1ss16AG, ccl-miR167a, aly-miR167b-3p_1ss13AG, ppe-miR172a-5p_1ss21GA, aly-miR167b-3p_1ss13AG, and ppe-miR172a-5p_1ss21GA) when ‘M8’ was separately compared with M3T2 vs. M8T2 (Fig. 4a). Moreover, six significantly differentially coexpressed miRNAs (aly-miR167b-3p_1ss13AG, ptc-miR172b-5p, ssl-miR828, csi-miR1515_1ss16AG, ath-miR858b and ppe-miR172a-5p_1ss21GA) were found in the comparison of the whole ‘M3’ and ‘M8’ groups, but surprisingly, only one differentially coexpressed miRNA (rco-miR167a_R+1) was found in ‘M3’ among all compared groups relative to M3T4 vs. M8T4 (Fig. 4b). Mining for these significantly differentially expressed and coexpressed miRNAs is helpful for exploring the expression differences and regulatory patterns between high- and low-oil cultivars of tea oil camellia at different developmental stages, which could provide a scientific basis for improved seed oil content and quality in tea oil camellia in the future. These results suggest that miRNA-mediated regulatory mechanisms may play an important role in various biological pathways during seed development in tea oil camellia.
Based on the analysis of the detected miRNAs, we further performed miRNA conservation analysis by mapping the miRNAs against selected species that have conserved evolutionary relationships and counted the frequency of occurrence of the same miRNAs in other species. The results showed that the conserved miRNA frequency was highest in the gma species, and a large amount of these conserved miRNAs also existed in the mdm, ptc and mes species (Table S9). In the family analysis, all known miRNAs were clustered into 35 miRNAs families, 24 of which (miR160, miR156, miR159, miR162, miR164, miR166, miR168, miR167, miR169, miR171, miR172, miR319, miR390, miR393, miR394, miR396, miR397, miR398, miR399, miR403, miR408, miR482, miR530, and miR535) contained many conserved miRNAs (Table S10), demonstrating that they may maintain a wide range of regulatory functions in seed development in tea oil camellia.
Prediction and identification of miRNA target genes
To further understand the functions of miRNAs in developing seeds of high- and low-oil cultivars of tea oil camellia, the target genes of the differentially expressed miRNAs were predicted using Target Finder software based on the miRNA sequences with mRNA transcriptome sequence data. A total of 17,166 target genes were predicted, and some genes were targeted by multiple miRNAs, resulting in a total of 33,418 predicted miRNA-target gene modules (Table S11). Furthermore, the predicted target genes of known and novel miRNAs were subjected to gene ontology (GO) analysis. The target genes were classified by the hypergeometric distribution algorithm based on the molecular function, biological process and cellular component categories by GO enrichment analysis (Fig. 5). The results showed that a total of 25,988 miRNA-target function modules were related to various biological processes, 21,132 of which were found to be involved in different cellular components, and 14,675 of which were identified to be involved in molecular function. Furthermore, the major biological functions of differentially expressed target genes were identified by the significant terms in the GO enrichment analysis. A total of 107 significant GO terms were identified, which were related to nucleus, sequence-specific DNA binding, regulation of transcription, DNA-template and electron carrier activity (Fig. 6a and Table S12). To further understand the biological functions of the target genes in the developing seeds of tea oil camellia, the target Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway genes were identified by KEGG pathway enrichment analysis. A total of 16,871 miRNA-target function modules were involved in 206 individual KEGG pathways, which were involved in linoleic acid metabolism, alpha-linolenic acid metabolism and glycosylated biosynthesis. Afterwards, the significant KEGG pathway enrichment analysis showed that 31 significant KEGG pathways (Fig. 6b and Table S13) were mainly involved in pyruvate metabolism, regulation of the actin cytoskeleton, fatty acid biosynthesis and citrate cycle.
Correlation analysis of differentially expressed miRNAs and their target genes involved in lipid metabolism
The seeds of the ‘M3’ cultivar of tea oil camellia have a high oil content and contain rich unsaturated fatty acids. The analyses of the differentially expressed miRNAs and their target genes involved in lipid metabolism revealed that multiple KEGG pathways were related to lipid metabolism during seed development by KEGG enrichment analysis (P<0.05) and that 1,869 miRNA-mRNA function modules were involved in 12 pathways related to lipid metabolism, including 1,723 predicted known miRNA-mRNA function modules and 146 novel miRNA-mRNA function modules, in which the same miRNA regulated different target genes. There were 188 target genes involved in the fatty acid biosynthesis pathway, 411 involved in the glycolysis/gluconeogenesis pathway, 118 involved in the biosynthesis of the unsaturated fatty acid pathway, and 159 related to the fatty acid metabolism pathway (Table S14). The significant KEGG enrichment analysis and function annotation were integrated, and four significant KEGG pathways related to lipid metabolism were identified, including pyruvate metabolism, fatty acid biosynthesis, glycolysis/gluconeogenesis, sphingolipid metabolism, and steroid biosynthesis. In addition, nine key miRNAs and their target genes related to the lipid metabolism pathway were identified, including one novel miRNA (PC-5p-21064_221) and eight known miRNAs that were involved in the seeds of high- and low-oil cultivars of tea oil camellia at four different stages. These genes are commonly involved in the seed oil biosynthesis pathway, such as ACC1 (comp59939_c1) targeted by tcc-miR162 and related to the pyruvate metabolism metabolic pathway (ko00620), KAS1 (comp67779_c0) and S-ACP-DES6 (comp67006_c0) targeted by csi-miR166e-5p, KAS3B (comp61049_c0) targeted by byrgl-miR5139_L-1, fabG (comp63911_c0) targeted by mtr-miR156h-3p_1ss8AC, Mcat (comp63026_c1) targeted by cpa-miR164d, FATB1 (comp67050_c0) targeted by PC-5p-21064_221, MOD1 (comp52017_c0) targeted by aly-miR393a-3p_1ss12TC, Δ9D (comp48800_c0) targeted by ath-miR172a and involved in the fatty acid biosynthesis metabolic pathway (ko00061), and GAPN (comp67185_c0) targeted by gma-MIR5368-p3_1ss21CT and involved in the glycolysis/gluconeogenesis metabolic pathway (ko00010) (Table 2).
To clarify the molecular mechanism of the above identified miRNAs and their target genes related to lipid metabolism in the high- and low-oil cultivars of tea oil camellia at four different stages, the expression patterns and expression differences of miRNAs and their target genes were analyzed by cluster analysis based on relevant differentially expressed miRNAs and their target gene expression levels at different seed developmental stages (Fig. 7). The results showed that the expression level of cpa-miR164d was the highest at the early seed developmental stage and was significantly higher in ‘M3’ than in ‘M8’ at the early seed developmental stage (June 2), and then it appeared to have a declining trend; in contrast, the expression level of the Mcat (comp63026_c1) gene targeted by cpa-miR164d was the lowest at the early seed developmental stage and then gradually increased in ‘M3’. These results indicated that there was a negative correlation between the expression level of cpa-miR164d and its target Mcat gene. The expression level of mtr-miR156h-3p_1ss8AC was significantly downregulated in ‘M3’ relative to ‘M8’ at the early seed developmental stage and then started to increase gradually at the late seed developmental stage; however, the fabG (comp63911_c0) gene targeted by mtr-miR156h-3p_1ss8AC was significantly upregulated. The expression level of fabG was highest in ‘M3’ at the early seed developmental stage, gradually decreased and reached the lowest expression level at the fourth developmental stage. The expression level of most miRNAs in the ‘M3’ cultivar decreased gradually, such as aly-miR393a-3p_1ss12TC, rgl-miR5139_L-1, PC-5p-21064_221, and ath-miR172a, from the second to third seed developmental stage; in contrast, the expression levels of their corresponding target genes, Mcat, KAS3B, FATB1 and Δ9D, appeared to gradually increase at these two stages.
The expression pattern of csi-miR166e-5p was similarly upregulated in the ‘M3’ and ‘M8’ cultivars; however, its target genes, KAS1 (comp67779_c0) and S-ACP-DES6 (comp67006_c0), appeared to have a trend of upregulation. The expression level of csi-miR166e-5p was highest at the fourth developmental stage, and there was a positive correlation between the expression levels of csi-miR166e-5p and its target genes. gma-MIR5368-p3_1ss21CT had a similar expression pattern in ‘M3’ and ‘M8’, and its expression level was low at the early seed developmental stage and then started to increase gradually. However, the GAPN (comp67185_c0) gene targeted by gma-MIR5368-p3_1ss21CT also appeared to have a similar expression pattern in ‘M3’ and ‘M8’, and its expression level was first downregulated at the early seed developmental stage and then upregulated, with a peak expression level at the late seed developmental stage. There was also a positively regulated relationship between the expression levels of gma-MIR5368-p3_1ss21CT and its target gene.
Based on these results, it could be speculated that these miRNAs and their target genes might be involved in lipid metabolism pathways. The analysis of their interactions and influence in different relevant lipid metabolism pathways in tea oil camellia will provide important information in the future.
Correlation analysis of differentially expressed miRNAs and their target genes involved in seed size
Seed size is one of the main factors affecting crop yield. Several miRNAs and their target genes belonging to MYB, CNR, and ARF transcription factors, which play important roles in controlling seed size, were identified (Table 3 and Fig. 8). The expression pattern of most miRNAs were slightly different between the high- and low-oil cultivars of tea oil camellia. For example, the expression levels of PC-5p-315953_8 were relatively low in ‘M8’ at the four seed developmental stages, but this low expression level appeared only at the first to third developmental stages in ‘M3’, as its expression increased sharply and showed a high expression level at the late seed developmental stage. The transcription factor CNR2 (comp59426_c0), a target of PC-5p-315953_8, was maintained at a low expression level in ‘M3’ at the late seed developmental stage, but its expression level was the highest at the third seed developmental stage; however, the expression level of CNR2 was the lowest in ‘M8’ at the third seed developmental stage.
The expression pattern of ath-miR858b appeared to have a gradually downregulated trend in ‘M3’ at different seed developmental stages, and its expression level was highest at the early seed developmental stage. The overall expression level of ath-miR858b was lower in the developing seeds of the ‘M8’ cultivar than in the developing seeds of the ‘M3’ cultivar; the lowest expression level of ath-miR858b appeared at the third developmental stage, and then it was gradually upregulated. MYB82 (comp278917_c0), a target transcription factor of ath-miR858b, showed a declining trend in the developing seeds of the ‘M3’ and ‘M8’ cultivars of tea oil camellia; its expression level was high at the early seed developmental stage and then declined gradually. The expression level of MYB3 (comp65157_c0), another target transcription factor of ath-miR858b, was high at the early seed developmental stage and then gradually declined; the expression level of MYB3 was lowest at the third seed developmental stage and was gradually upregulated from the third to the fourth developmental stage. At the third seed developmental stage, the expression level of MYB3 was lower in ‘M3’ than in ‘M8’, and the rate of increase in ‘M3’ was slower than that in ‘M8’ from the third to the fourth seed developmental stage. In addition, the expression level of MYB44 (comp31235_c0), another target transcription factor of ath-miR858b, showed a slight difference between the ‘M3’ and ‘M8’ cultivars. The expression level of MYB44 in ‘M3’ was highest at the early seed developmental stage and then decreased rapidly to the lowest value at the second seed developmental stage; after this, it gradually increased. In contrast, MYB44 consistently maintained a high expression level and did not change significantly in the developing seeds of ‘M8’.
The expression level of mtr-miR156h-3p_1ss8AC was lowest at the early seed developmental stage in the ‘M3’ and ‘M8’ cultivars; however, the expression level of its corresponding target transcription factor MYB44 (comp55270_c0) was highest in the ‘M3’ and ‘M8’ cultivars at the early seed developmental stage. The expression level of hpe-miR162a_L-2 in ‘M3’ was low during all seed developmental stages, but it was high in ‘M8’ at the early seed developmental stage and then gradually decreased to a low level at the late seed developmental stage. ARF19 (comp57222_c0), a target transcription factor of hpe-miR162a_L-2, had a similar expression pattern in the ‘M3’ and ‘M8’ cultivars from the early to third seed developmental stage, with high expression levels at the early seed developmental stage; after this, its expression gradually decreased but increased again at the third seed developmental stage. From the third to fourth seed developmental stage, there were obvious differences in the expression trends for ARF19 between the ‘M3’ and ‘M8’ cultivars. The expression level of ARF19 was downregulated in ‘M3’ but upregulated in ‘M8’.
mdm-miR167h_1ss22AT had a similar expression pattern in the ‘M3’ and ‘M8’ cultivars of tea oil camellia, with low expression levels at early seed developmental stages and high expression levels at late seed developmental stages. MED27 (comp135267_c0), a target transcription factor of mdm-miR167h_1ss22AT, had the highest expression level in ‘M3’ at the early seed developmental stage and the lowest expression level at the late seed developmental stage; however, the expression level of MED27 remained high in ‘M8’ throughout all seed developmental stages. The expression level of MED27 in ‘M8’ was lower than that in ‘M3’ at the early seed developmental stage, but it was higher than that in the ‘M3’ cultivar at other seed developmental stages.
miRNA-mediated regulation of multiple transcription factors is complex but coordinated with each transcription factor to enable their biological functions. The miRNAs that regulate the expression of target genes directly affect seed size; however, whether these predicted target genes can control seed size in tea oil camellia remains unknown.
Correlation analysis of differentially expressed miRNAs and their target genes involved in growth and development
According to the target prediction annotation information, 12 key miRNAs and their target genes related to growth and development were identified between the high- and low-oil cultivars of tea oil camellia during the different seed developmental stages. Based on the heat map analysis, the target genes targeted by relevant miRNAs were ubiquitously expressed in the developing seeds of tea oil camellia, and the different miRNA expression levels were significantly different in the ‘M3’ and ‘M8’ cultivars (Table 4 and Fig. 9). For example, the expression level of mtr-miR156h-3p_1ss8AC was significantly lower in ‘M3’ than in ‘M8’ at the early developmental stage and then increased slowly; however, the expression level of miR156h-3p_1ss8AC was consistently high in ‘M8’ throughout the whole seed developmental stage. TCP24 (comp65957_c0), a target of miR156h-3p_1ss8AC, had a high expression level in ‘M3’ at the early seed developmental stage and was significantly upregulated; however, it was consistently low in ‘M8’ at all seed developmental stages. The relationship between miR156h-3p_1ss8AC and its target, TCP24, was positively regulated.
The expression levels of PC-3p-9080_426 showed significant differences in the ‘M3’ and ‘M8’ cultivars. PC-3p-9080_426 was upregulated gradually in ‘M8’ from the early to third seed developmental stage but had the lowest expression level at the fourth seed developmental stage; in contrast, PC-3p-9080_426 consistently had a higher expression level in ‘M3’ than in ‘M8’ at the whole seed developmental stage, although it showed a slightly decreasing trend at the late seed developmental stage. The COL13 (comp63752_c0) gene targeted by PC-3p-9080_426 was downregulated in both cultivars of tea oil camellia at the late seed developmental stage, although its expression level was high in ‘M3’ at the early developmental stage and was then downregulated. The expression pattern of ERF115 (comp58165_c0), another target gene of PC-3p-9080_426, was similar to that of the COL13 gene in the ‘M3’ cultivar. However, in the ‘M8’ cultivar, the expression level of ERF115 was highest at the early seed developmental stage; after this, it decreased rapidly to a low expression level at the third seed developmental stage and then had a slightly increased expression level at the late seed developmental stage.
The expression pattern of han-miR156a_L+1 was similar in the ‘M3’ and ‘M8’ cultivars of tea oil camellia. Its expression level was significantly low at the early seed developmental stage and then increased sharply, reaching the highest level at the late seed developmental stage. The two genes SPL4 (comp60372_c0) and SBP2 (comp62492_c0) targeted by han-miR156a_L+1 had the highest expression levels at the early seed developmental stage in the two cultivars of tea oil camellia and the lowest expression levels at the late developmental stage. The expression pattern of rco-miR156e_L+1R-1 was similar to that of han-miR156a_L+1, but the expression patterns of its three target genes, ACS1 (comp66521_c0), SBP2 (comp62492_c0) and SBP1 (comp61590_c1), contrasted with that of the miRNA han-miR156a_L+1, showing a downward trend in developing seeds.
According to these results, it could be found that these miRNAs and their target genes have highly correlated relationships in expression levels during seed development, and the differences in the expression levels between the two cultivars or among the different seed developmental stages might exert a great influence on the growth and development of tea oil camellia seeds.
Correlation analysis of differentially expressed miRNAs and their target genes involved in resistance
According to the enrichment analysis of the target genes of differentially expressed miRNAs, a total of 11 functional modules of miRNAs and targeted genes related to stress resistance were identified, mainly involving two categories: abiotic stresses and biotic stresses (Table 5). In addition, two miRNA-mRNA functional modules, which may affect the yield and quality of tea oil camellia seed, were also identified. Their expression patterns and levels appeared to be obviously different between the ‘M3’ and ‘M8’ cultivars of tea oil camellia and among different seed developmental stages (Fig. 10).
The results of the heat map analysis showed that the majority of miRNAs involved in stress resistance displayed different expression profiles in two cultivars of tea oil camellia at four different seed developmental stages. The expression patterns of ath-miR858b and nta-MIR6149a-p5_2ss18CT21TA were different in the high- and low-oil cultivars. Their expression levels were consistently high in ‘M3’ throughout the four seed developmental stages; however, for the ‘M8’ cultivar, their expression levels were high at the early and second seed developmental stages, gradually decreased to the lowest level at the third seed developmental stage, and then slowly increased. In contrast, the expression level of MYB4 (comp31687_c0), a target gene of ath-miR858b, was lowest in ‘M8’ at the fourth seed developmental stage; the expression level of CPR30 (comp65951_c0), a target gene of nta-MIR6149a-p5_2ss18CT21TA, was lowest in ‘M3’ at the second seed developmental stage, but its highest expression level in ‘M8’ appeared at the third seed developmental stage relative to the other stages.
The expression level of mtr-miR156h-3p_1ss8AC was lowest in ‘M3’ at the early seed developmental stage and then increased gradually; however, the expression level of mtr-miR156h-3p_1ss8AC remained consistently high in ‘M8’ throughout the four different seed developmental stages. Interestingly, the two genes BAG7 (comp63682_c0) and NRT1.7 (comp66713_c0) targeted by mtr-miR156h-3p_1ss8AC showed low expression levels in ‘M3’ at the early seed developmental stage and then decreased gradually to the lowest expression level at the late seed developmental stage. Although the expression levels of these two target genes were also lowest in the ‘M8’ cultivar at the late seed developmental stage, they were higher than those in ‘M3’.
The expression levels of han-miR156a_L 1 in the ‘M3’ and ‘M8’ cultivars were lowest at the early seed developmental stage and then increased gradually. In contrast, the expression level of its target gene, SPL16 (comp56544_c0), was highest at the early seed developmental stage and then decreased. The expression level of ppe-miR172a-5p_1ss21GA remained high in ‘M3’ at all seed developmental stages, but its expression level in ‘M8’ was only high at the early seed developmental stage and then was significantly downregulated after the second seed developmental stage, with the lowest expression level at the fourth seed developmental stage. The expression of the PAB2 (comp65585_c0) gene targeted by ppe-miR172a-5p_1ss21GA was lowest in ‘M3’ at the late seed developmental stage, and its highest expression level in ‘M8’ appeared at the early seed developmental stage.
Validation of differentially expressed miRNAs and their target genes by qRT-PCR
To validate the high-throughput sequencing data, six functional modules of differentially expressed miRNAs and their target genes were selected, han-miR156a_L+1-SPL4, rco-miR156e_L+1R-1-SBP1, mtr-miR156h-3p_1ss8AC-TCP24, ppe-miR172a-5p_1ss21GA-PAB2, nta-MIR6149a-p5_2ss18CT21TA-CPR30 and PC-3p-9080_426-ERF115. The differences in their expression patterns and expression levels in the high- and low-oil cultivars of tea oil camellia at the four different seed developmental stages were identified by qRT-PCR analysis (Fig. 11 and Fig. 12). The results showed that six pairs of functional modules of miRNAs and their target genes showed differential expression, and the expression patterns of most miRNAs and their target genes were consistent with the results obtained by high-throughput small RNA sequencing data. All six selected miRNAs showed opposite expression trends relative to their target genes.
han-miR156a_L 1 was upregulated in high- and low-oil cultivars of tea oil camellia, and its expression level was lowest at the early seed developmental stage and highest at the late developmental stage (Fig. 11a). The expression level of han-miR156a_L 1 in the ‘M3’ cultivar was higher than that in the ‘M8’ cultivar. In contrast, the expression level of the target gene SPL4 of this miRNA was highest at the early seed developmental stage and then decreased sharply to the lowest expression level at the late seed developmental stage (Fig. 12a). The expression level of rco-miR156e_L 1R-1 was upregulated in the ‘M8’ cultivar; it gradually increased from the early to third seed developmental stage, peaking at the third seed developmental stage, and then decreased slightly (Fig. 11b). The expression level of rco-miR156e_L 1R-1 in the ‘M3’ cultivar was upregulated continuously and reached the highest value at the late seed developmental stage but was always lower than that in the ‘M8’ cultivar (Fig. 11a). Interestingly, the target gene, SBP1, of rco-miR156e_L 1R-1 was downregulated in the ‘M3’ and ‘M8’ cultivars of tea oil camellia, and its expression level was highest at the early seed developmental stage and lowest at the late developmental stage (Fig. 12b).
The expression level of mtr-miR156h-3p_1ss8AC was upregulated in the ‘M3’ cultivar and was significantly lower than that in the ‘M8’ cultivar at the early seed developmental stage; after this, it showed a rapidly increasing trend, with the highest expression level at the late seed developmental stage (Fig. 11c). In contrast, the expression level of mtr-miR156h-3p_1ss8AC in the ‘M8’ cultivar was relatively high at all seed developmental stages. Surprisingly, the expression level of the target gene TCP24 of mtr-miR156h-3p_1ss8AC in the ‘M3’ cultivar was also low at the early seed developmental stage (Fig. 12c) but was then rapidly upregulated; however, its expression level in the ‘M8’ cultivar remained low at all seed developmental stages (Fig. 12c). The expression level of miR156h-3p_1ss8AC was found to be positively correlated with the expression of its target gene in the high-oil ‘M3’ cultivar but negatively correlated with the expression of its target gene in the low-oil ‘M8’ cultivar.
The expression level of ppe-miR172a-5p_1ss21GA was highest in ‘M3’ at the second seed developmental stage (Fig. 11d), while the expression level of its target gene PAB2 was lower at the second seed developmental stage than at the early and third seed developmental stages. The expression pattern of ppe-miR172a-5p_1ss21GA was similar to that of its target gene in the ‘M8’ cultivar, and both were first upregulated and then downregulated (Fig. 11d and Fig. 12d). The expression level of ppe-miR172a-5p_1ss21GA was found to be negatively correlated with the expression of its target gene in high-oil ‘M3’ but positively correlated with the expression of its target gene in low-oil ‘M8’.
nta-MIR6149a-p5_2ss18CT21TA had the highest expression level in ‘M3’ at the third seed developmental stage, but its expression level was lowest in ‘M8’ at the third seed developmental stage and was significantly lower than that in ‘M3’ (Fig. 11e). Its corresponding target gene, CPR30, was low in the ‘M3’ cultivar at the fourth developmental stage. The expression of the target gene CPR30 of nta-MIR6149a-p5_2ss18CT21TA in ‘M3’ was higher than that in ‘M8’ at the early seed developmental stage but was lower in ‘M3’ than ‘M8’ at the other three seed developmental stages (Fig. 12e). The expression level of PC-3p-9080_426 peaked in ‘M3’ at the second seed developmental stage and in ‘M8’ at the third seed developmental stage (Fig. 11f). It was first upregulated and then downregulated in ‘M3’. The target gene ERF115 of nta-MIR6149a-p5_2ss18CT21TA showed the lowest expression level in ‘M8’ at the third seed developmental stage, and its expression level first increased and then decreased in low-oil ‘M3’ (Fig. 12f). The expression level of ERF115 in ‘M3’ was positively correlated with the expression level of PC-3p-9080_426.
The above results obtained by qRT-PCR validation are consistent with the sequencing results for the expression patterns and levels of miRNAs and their targeted genes (Fig. 11 and Fig. 12), indicating the accuracy and high efficiency of small RNA high-throughput sequencing results and reflecting the dependability of the qRT-PCR method used to identify the temporal and spatial expression characteristics and expression differences of related functional miRNAs and target genes in tea oil camellia at different seed developmental stages.
Relationship validation of the cpa-miR393_R-1-AFB2 functional module by luciferase activity assays
A target prediction analysis server psRNA target was used to assess the complementarity between cpa-miR393_R-1 and the target site of AFB2. The potential to target AFB2 3¢UTR of cpa-miR393_R-1 was predicted (Fig. 13a). Luciferase activity in 293 cells cotransfected with the cpa-miR393_R-1 recombinant expression vector and the expression vector containing the 3¢UTR of AFB2 fused with the reporter gene decreased by nearly 36.26% (p < 0.05) compared to that in the control group (Fig. 13b). These results indicate that AFB2 is one of the target genes of cpa-miR393_R-1.