Integrative Metabolic and Transcriptomic Profiling in Camellia oleifera and Camellia meiocarpa Uncover Potential Mechanisms That Govern Triacylglycerol Degradation during Seed Desiccation

Camellia seed oil is a top-end quality of cooking oil in China. The oil quality and quantity are formed during seed maturation and desiccation. So far, it remains largely unresolved whether lipid degradation occurs and contributes to Camellia oil traits. In this study, three different Camellia germplasms, C. oleifera cv. Min 43 (M43), C. meiocarpa var. Qingguo (QG), and C. meiocarpa cv Hongguo (HG) were selected, their seed oil contents and compositions were quantified across different stages of seed desiccation. We found that at the late stage of desiccation, M43 and QG lost a significant portion of seed oil, while such an event was not observed in HG. To explore the molecular bases for the oil loss In M43, the transcriptomic profiling of M43 and HG was performed at the early and the late seed desiccation, respectively, and differentially expressed genes (DEGs) from the lipid metabolic pathway were identified and analyzed. Our data demonstrated that different Camellia species have diverse mechanisms to regulate seed oil accumulation and degradation, and that triacylglycerol-to-terpenoid conversion could account for the oil loss in M43 during late seed desiccation.


Introduction
Camellia oleifera Abel., a member of the Theaceae family, is widely distributed in the subtropical mountain areas of the Yangtze River basin and South China [1,2]. As evergreen broadleaf shrubs or small trees, Camellia is tolerant to drought and barren land [3]. C. oleifera is mainly used for oil production from its seed kernel, which contains~50% of oilseed dry mass [4], and oleic acid account for~80% of its seed oil. Thus, Camellia seed oil is reputed as the "Oriental olive oil" [5]. In addition, the Camellia oil is rich in bioactive compounds, including polyphenol, squalene, and vitamin E, etc. [6]. Thus, Camellia oil shows potent activities of antioxidation, antiaging, anti-inflammation, anticancer, and antibacterial. The long-term consumption of Camellia oil could lower the levels of serum triglycerides and cholesterol, enhance immunity, stabilize blood sugar levels, reduce cardiovascular disease, and prevent hypertension [5].
There are a lot of interests to understand the factors affecting Camellia oil synthesis. Recently, transcriptomic technology has been applied to understand lipid synthesis and response to environmental stresses in Camellia oleifera [3,[7][8][9][10][11]. Zeng et al. (2014) reported that the mRNA levels of fructose-1,6-bisphosphate aldolase (CoFBA) and stearoyl-ACP desaturases (CoSAD) are correlated with oil content, whereas fatty acid desaturase 2 (Co-FAD2) gene expression levels are correlated with fatty acid composition [12]. Feng et al. (2017) analyzed transcriptomic changes during postharvest seed drying [10]. The authors found that natural drying improved the quality and quantity of Camellia seed oil; in addition, putative transcripts were identified that are potentially involved in lipid metabolism during postharvest drying. Lin et al. (2018) performed the transcriptomic profiling of Camellia seeds at four developmental stages, and found that the distribution of the DEG numbers was associated with seed oil accumulation patterns, the authors suggest that the regulation of CoSAD and CoFAD2 is critical to boost oleic acid level at the late stages of seed development [13]. Wu et al. (2019) reported that the coordinated higher expression levels of the upstream gene (HAD, EAR and KASI) are directly associated with increased levels of palmitic acid (C16:0), and continuous higher expression levels of the SAD gene could accelerate oleic acid synthesis and accumulation, while coordinated lower expression levels of the downstream genes (FAD2, FAD3, FAD7, FAD8 and FAE1) are correlated with a decreased oleic acid conversion [11]. Gong et al. (2020) applied iso-seq to obtain a fulllength transcriptome from Camellia oleifera seeds [14]. The authors found that transcript variants could be involved in seed oil biosynthesis. Lin et al. (2022) performed a genomewide association study, and found that during Camellia domestication, sugar-dependent triacylglycerol lipase 1 (CoSDP1), β-ketoacyl-acyl carrier protein synthase III (KAS III), and CoSAD play important roles in enhancing the yield and quality of Camellia seed oil [15].
For oil seeds such as B. napus, Arabidopsis, Crambe abyssinica, and Nicotiana tabacum, lipid contents were reported to decrease during seed maturation [16][17][18]. In B. napus, the enzyme activities and abundance of β-oxidation, glyoxylate cycle, and phosphoenolpyruvate carboxykinase were increased during seed desiccation compared with the oil accumulation stage, seed lipid degradation was not associated with net gluconeogenesis activity, and partial lipids were released as CO 2 [19]. These data suggest that B. napus seed oils were actively degraded or converted into other metabolites other than sugars. Zhou et al. (2013) and Liu (2020) reported that harvest time affects seed oil content in C. oleifera, and that it remains unclear whether seed oil degradation accounts for it [20,21].
In China, C. oleifera is the dominant cultivated oil-tree species [22]. Until now, limited attention has been paid to other Camellia oil-tree species, including C. meiocarpa and C. chekiangoleosa [10,23]. In this study, two varieties of C. meiocarpa, were included, they are: C. meiocarpa var. Hongguo (HG) and C. meiocarpa var. Qingguo (QG). The obvious phenotypic difference between these two germplasms is the timing of fruit color changes: on September, HG turns its green fruit color to red when seeds start to mature, while QG keeps the green color of its fruits at same stage (Supplementary Figure S1). To investigate whether there are species-or variety-specific effects on Camellia seed oil accumulation or degradation, the seeds from HG, QG, and C. oleifera cv. Min 43 (M43) were harvested weekly from late September to mid-November, which span seed maturation and desiccation. Seed oil contents and compositions were then measured and compared. We found that at the late stage of desiccation, M43 and QG lost a significant portion of seed oil, while such oil loss was not observed in HG. To elucidate the underlying molecular bases, transcriptomic profiling from M43 and HG was performed at the early (26 September, 47 weeks after anthesis) and the late (1 November, 52 weeks after anthesis) desiccation stages, respectively, and differentially expressed genes (DEGs) in the lipid metabolic pathway were identified and analyzed. Our data provide novel insights about the potential mechanisms responsible for seed oil degradation in M43 during preharvest seed desiccation stages.
The seed oil contents among these three germplasms showed diverse accumulation patterns during seed maturation and desiccation: at the 46th week after anthesis, both QG and HG had similar seed oil contents; in the following week, both showed large differences in oil accumulation rates where the oil content of HG increased rapidly, followed by a much gentle increase. In contrast, the seed oil content from QG showed a slower but steady increase until the 53rd week, followed by a sharp decrease. From the 53rd to the 54th week, QG lost~25% of its seed oil reserve. The seed oil content from M43 showed biphasic on the 49th and the 52nd week, respectively. From the 52nd to 54th week, M43 lost~30% of its seed oil reserve. In contrast, at the same period, the seed oil content of HG even slightly increased. Despite their diverse oil accumulation patterns, these three cultivars showed similar seed oil contents on the 54th week after anthesis ( Figure 1).
The seed oil contents among these three germplasms showed diverse accum patterns during seed maturation and desiccation: at the 46th week after anthesis, b and HG had similar seed oil contents; in the following week, both showed large ences in oil accumulation rates where the oil content of HG increased rapidly, fo by a much gentle increase. In contrast, the seed oil content from QG showed a slow steady increase until the 53rd week, followed by a sharp decrease. From the 53rd 54th week, QG lost ~25% of its seed oil reserve. The seed oil content from M43 s biphasic on the 49th and the 52nd week, respectively. From the 52nd to 54th wee lost ~30% of its seed oil reserve. In contrast, at the same period, the seed oil conten even slightly increased. Despite their diverse oil accumulation patterns, these thre vars showed similar seed oil contents on the 54th week after anthesis ( Figure 1).

The Accumulation Patterns of Individual Fatty Acid Species during Seed Desiccation
The contents of individual fatty acid species were quantified during various st seed desiccation; the data were presented in Figure 2 and Supplementary Table oleic acid, the dominant fatty acid, showed diverse changing patterns among thes germplasms ( Figure 2). In QG, the oleic acid showed an almost linear increase fr 46th to the 53rd week after anthesis, then followed by a 26.4% decrease on the 54t Figure 1. The total seed oil contents during seed desiccation. Data were expressed as mean ± SD (n = 3).

The Accumulation Patterns of Individual Fatty Acid Species during Seed Desiccation
The contents of individual fatty acid species were quantified during various stages of seed desiccation; the data were presented in Figure 2 and Supplementary Table S1. The oleic acid, the dominant fatty acid, showed diverse changing patterns among these three germplasms ( Figure 2). In QG, the oleic acid showed an almost linear increase from the 46th to the 53rd week after anthesis, then followed by a 26.4% decrease on the 54th week after anthesis ( Figure 2, left upper panel). In HG, the oleic acid showed a rapid increase from the 46th to the 47th week after anthesis, then remained almost stagnant until the 49th week. At the 50th week, the oleic acid (18:1 ∆9 ) slightly decreased, followed by a constant increase until the 54th week. Interestingly, on the 50th week, the drop in 18:1 ∆9 is accompanied by the concurrent increase in linoleic acid, suggesting that the decline in oleic acid could be due to its accelerated conversion into linoleic acid at this time point (Figure 2, middle upper panel). In M43, the oleic acid showed two peaks, with the first small peak presented on the 49th week, the second large peak presented on the 52nd week, then followed by a 26.7% decrease on the 54th week after anthesis. The linoleic acid content did not show coordinated changes with oleic acid as we saw in HG ( Figure 2, right upper panel). The linoleic acid (18:2) is the major polyunsaturated fatty acid, and accounts for 3.9% to 4.4% of seed dry mass; palmitic acid (16:0) is the dominant saturated fatty acid in Camellia seed oil, and accounts for 3.3% to 4.0% of seed dry mass. During seed desiccation, the linoleic acid and palmitic acid contents kept fairly constant ( Figure 2, upper panel). after anthesis ( Figure 2, left upper panel). In HG, the oleic acid showed a rapid increase from the 46th to the 47th week after anthesis, then remained almost stagnant until the 49th week. At the 50th week, the oleic acid (18:1 Δ9 ) slightly decreased, followed by a constant increase until the 54th week. Interestingly, on the 50th week, the drop in 18:1 Δ9 is accompanied by the concurrent increase in linoleic acid, suggesting that the decline in oleic acid could be due to its accelerated conversion into linoleic acid at this time point (Figure 2, middle upper panel). In M43, the oleic acid showed two peaks, with the first small peak presented on the 49th week, the second large peak presented on the 52nd week, then followed by a 26.7% decrease on the 54th week after anthesis. The linoleic acid content did not show coordinated changes with oleic acid as we saw in HG ( Figure 2, right upper panel). The linoleic acid (18:2) is the major polyunsaturated fatty acid, and accounts for 3.9% to 4.4% of seed dry mass; palmitic acid (16:0) is the dominant saturated fatty acid in Camellia seed oil, and accounts for 3.3% to 4.0% of seed dry mass. During seed desiccation, the linoleic acid and palmitic acid contents kept fairly constant ( Figure 2, upper panel).   The oil contents for steric acid (18:0), cis-11-octadecenoic acid (18:1 ∆11 ), linolenic acid (18:3), and arachidic acid (20:0) were in the range of 0.1-1.4% of seed dry mass. Beginning on the 50th week, 18:0 showed a general upward trend, while 18:1 ∆11 showed a downward trend from these three germplasms. Beginning on the 52nd week, the 18:0 contents from QG and M43 started to decrease, while its content in HG kept rising. The contents for 18:3 and 20:0 only showed small variations during seed desiccation (Figure 2 middle panel). Cis-11-hexadecenoic acid (16:1 ∆11 ) and paullinic acid (20:1) accounted for less than 0.1% of seed dry mass ( Figure 2, lower panel). Their variations would not significantly affect the total seed oil contents.
Next, we expressed the seed oil content data as mol percent (Mol%), then examined fatty acid compositional changes during seed desiccation. The proportion of oleic acid from QG and M43 showed general increasing trends with seed desiccation. Meanwhile, the relative proportion of palmitic acid and linoleic acid gradually decreased ( Figure 3). The oleic acid proportion from HG showed two peaks, with the first peak appearing on the 48th week after anthesis, followed by a decrease until the 50th week, then increase again until the 54th week after anthesis. Accordingly, the linoleic acid proportion showed an opposite changing trend to that of oleic acid ( Figure 3, middle panel). The other six fatty acids were minor components, and only showed small variations during seed desiccation ( Figure 3). on the 50th week, 18:0 showed a general upward trend, while 18:1 Δ11 showed a downw trend from these three germplasms. Beginning on the 52nd week, the 18:0 contents f QG and M43 started to decrease, while its content in HG kept rising. The contents for and 20:0 only showed small variations during seed desiccation ( Figure 2 middle pan Cis-11-hexadecenoic acid (16:1 Δ11 ) and paullinic acid (20:1) accounted for less than 0.1% seed dry mass ( Figure 2, lower panel). Their variations would not significantly affect total seed oil contents. Next, we expressed the seed oil content data as mol percent (Mol%), then exami fatty acid compositional changes during seed desiccation. The proportion of oleic from QG and M43 showed general increasing trends with seed desiccation. Meanwh the relative proportion of palmitic acid and linoleic acid gradually decreased (Figur The oleic acid proportion from HG showed two peaks, with the first peak appearing the 48th week after anthesis, followed by a decrease until the 50th week, then incre again until the 54th week after anthesis. Accordingly, the linoleic acid proportion show an opposite changing trend to that of oleic acid ( Figure 3, middle panel). The other fatty acids were minor components, and only showed small variations during seed de cation ( Figure 3). The unsaturated fatty acid (USFA)-to-saturated fatty acid (SFA) ratios showed a g erally increasing trend during desiccation among the three germplasms ( Figure 4, up panel). The monounsaturated fatty acid (MUFA)-to-polyunsaturated fatty acid (PUFA tios from HG showed large variations during seed desiccation, with the first peak app ing on the 48th week after anthesis, then reached to nadir on the 50th week. After that The unsaturated fatty acid (USFA)-to-saturated fatty acid (SFA) ratios showed a generally increasing trend during desiccation among the three germplasms ( Figure 4, upper panel). The monounsaturated fatty acid (MUFA)-to-polyunsaturated fatty acid (PUFA) ratios from HG showed large variations during seed desiccation, with the first peak appearing on the 48th week after anthesis, then reached to nadir on the 50th week. After that, the ratio kept increasing until the 54th week. In contrast, QG and M43 showed a stable increase until the 53rd week ( Figure 4, middle panel). The ratios of 18C to 16C fatty acid shared similar changing trends as the USFA-to-SFA ratios ( Figure 4, lower panel), largely due to the fact that oleic acid (18:1 ∆9 ) and palmitic acid (16:0) are the dominant USFA and SFA, respectively.

The Correlations among Total Seed Oil Contents and Individual Fatty Acid Species
Fatty acids are the major components in Camellia seed oil (triacylglycerol, TAG). To better understand the relationships between oil contents and fatty acid species, below, we briefly describe the major steps and enzymes for fatty acid synthesis ( Figure 5a). In plastid, beta-ketoacyl-ACP synthase I (KAS I) uses acetyl-CoA (Ac-CoA) as substrate to synthesize fatty acids with a chain length less than 16C, then KAS II converts 16C fatty acid (16:0) into 18C fatty acid (18:0) [24,25]. 18:0 is either exported into cytoplasm or desaturated to 18:1 by plastidial steryl-ACP desaturase (SAD) before exportation [26]. After exportation, partial 18C fatty acid is converted to 20C fatty acids by the endoplasmic reticulum (ER)associated fatty acid elongase complex (FAE) [27]. Thus, almost all the 18C and 20C fatty acids are synthesized through KAS II except 18:1 Δ11 and 20:1 Δ13 , which are synthesized by
Among individual fatty acid species, both positive and negative correlations were identified. In all three germplasms, 18:1 Δ9 is positively correlated with 18:0, 18:3, and 20:1 (Figure 5b-d); the positive correlation between 18:2 and 18:3 was only observed in QG (Figure 5b). In both varieties of C. meiocarpa, 18:1 Δ11 was negatively correlated with 18:0 ( Figure 5b,c), while such correlation was not found in M43 (Figure 5d).  The correlations among total seed oil contents and individual fatty acid contents were analyzed (Supplementary Figure S4). In all three germplasms, total seed oil contents were significantly and positively correlated with the contents of 18:0, 18:1 ∆9 , and 18:3, but not correlated with the contents of 16:1 ∆11 , 20:0, and 18:2 (Figure 5b-d). In both of the C. meiocarpa germplasms (QG and HG), 20:1 contents were also positively correlated with total seed oil contents (Figure 5b,c), while no such correlation was observed in C. oleifera cv. Min 43 (Figure 5d). In QG and M43, total seed oil contents were also positively correlated with 16:0 contents (Figure 5b

The Synchronization of Enzyme Activity Fluctuation during Seed Desiccation
As described in Figure 5a, all the 18C and 20C fatty acids (except 18:1 ∆11 and 20:1 ∆13 ) are synthesized through KAS II. When the absolute contents of 18C and 20C fatty acids (except 18:1 ∆11 and 20:1 ∆13 ) are added up at each sampling week, the difference relative to the previous week reflect the average activity changes of KAS II during that week in planta. Similarly, the palmitic acid content changes relative to the previous week represent FatB activity changes. Since all the 18:1 ∆9 , 18:2, and 18:3 are synthesized through SAD, and 18:2 and 18:3 are synthesized through FAD2, the total seed content changes of 18:1 ∆9 , 18:2, and 18:3 relative to the previous week reflect SAD activity changes. The total seed content changes of 18:2 and 18:3 relative to the previous week reflect FAD2 activity changes, while the 18:3 content changes relative to the previous week reflect FAD3 activity changes. The total content changes of cis-11-octadecenoic acid (18:1 ∆11 ) and paullinic acid (20:1 ∆13 ) relative to the previous week reflect FAE activity changes in planta (Figure 5a). Based on the absolute seed fatty acid contents as presented in Figure 2 and Supplementary Table S1, the weekly average enzyme activity fluctuations for KAS II, SAD, FatB, FAD2, FAD3, and FAE were plotted and presented in Figure 6. In QG and M43, all these enzymes showed highly synchronized fluctuation patterns ( Figure 6, left and middle panel). In HG, these synchronized fluctuation patterns were sustained until the 50th week after anthesis, after that, FAD2 showed discordant change trends with KAS II and SAD ( Figure 6, right panel). Considering that QG and M43 showed rapid oil accumulation, whereas the oil accumulation in HG was much slower (Figure 1), we reasoned that these concerted enzyme activity fluctuations in fatty acid synthesis could maximize seed oil biosynthetic efficiency.

The Synchronization of Enzyme Activity Fluctuation during Seed Desiccation
As described in Figure 5a, all the 18C and 20C fatty acids (except 18:1 Δ11 and 20:1 Δ13 ) are synthesized through KAS II. When the absolute contents of 18C and 20C fatty acids (except 18:1 Δ11 and 20:1 Δ13 ) are added up at each sampling week, the difference relative to the previous week reflect the average activity changes of KAS II during that week in planta. Similarly, the palmitic acid content changes relative to the previous week represent FatB activity changes. Since all the 18:1 Δ9 , 18:2, and 18:3 are synthesized through SAD, and 18:2 and 18:3 are synthesized through FAD2, the total seed content changes of 18:1 Δ9 , 18:2, and 18:3 relative to the previous week reflect SAD activity changes. The total seed content changes of 18:2 and 18:3 relative to the previous week reflect FAD2 activity changes, while the 18:3 content changes relative to the previous week reflect FAD3 activity changes. The total content changes of cis-11-octadecenoic acid (18:1 Δ11 ) and paullinic acid (20:1 Δ13 ) relative to the previous week reflect FAE activity changes in planta (Figure 5a). Based on the absolute seed fatty acid contents as presented in Figure 2 and Supplementary Table S1, the weekly average enzyme activity fluctuations for KAS II, SAD, FatB, FAD2, FAD3, and FAE were plotted and presented in Figure 6. In QG and M43, all these enzymes showed highly synchronized fluctuation patterns ( Figure 6, left and middle panel). In HG, these synchronized fluctuation patterns were sustained until the 50th week after anthesis, after that, FAD2 showed discordant change trends with KAS II and SAD ( Figure 6, right panel). Considering that QG and M43 showed rapid oil accumulation, whereas the oil accumulation in HG was much slower (Figure 1), we reasoned that these concerted enzyme activity fluctuations in fatty acid synthesis could maximize seed oil biosynthetic efficiency.  cantly higher than that of HG, then decreased dramatically, while such a reduction was not observed in HG (Figure 1). To explore the potential mechanisms governing seed oil loss in M43 at the late phase of seed desiccation, seeds were harvested on the 47th and the 52nd week after anthesis for transcriptomic analysis.
The four libraries of HG contained 38.99 G clean reads, and bioinformatics analysis identified 46,620 unigenes (Supplementary Table S2 Tables S6 and S7). To identify homologous genes between M43 and HG, their protein sequences were blasted to each other by using reciprocal best BLAST hit (RBH) method, and~30,000 homologous genes were identified between M43 and HG ( Supplementary Tables S8 and S9). Compared to the 47th week after anthesis, on the 52nd week, 5940 and 7644 unigenes were differentially expressed from HG and M43, respectively (Supplementary Tables S10 and S11). GO term enrichment analysis demonstrated that in HG, the downregulated unigenes are enriched in monooxygenase activity, oxidoreductase activity, and toxin activity (Supplementary Figure S5a), whereas the upregulated unigenes are enriched in response to chitin (Supplementary Figure S5b). In M43, the downregulated unigenes are enriched in photosystem I stabilization (Supplementary Figure S6a), while the upregulated unigenes are enriched in chitin binding, chitinase activity and monooxygenase activity, defense response, response to wounding, response to chitin, camalexin biosynthetic process, and cellular response to sulfur starvation (Supplementary Figure S6b).
The genes involved in lipid pathway were further screened. In total, 382 homologous unigenes for lipid metabolism were functionally annotated (Supplementary Table S12). Within lipid pathway, the DEGs between the 47th and the 52nd week samples were screened. In total, 75 and 107 DEGs were identified from HG and M43, respectively (Supplementary Table S13). Comparing to the 47th week after anthesis, in HG, 29 out of the 75 genes were upregulated and 46 genes downregulated; in M43, 60 of the 107 genes were upregulated and 47 genes downregulated (Supplementary Table S13).
To identify the potential genes involved in the TAG degradation of M43, we reasoned that those putative genes would exhibit relatively higher expression levels in M43 compared with that of HG on the 52nd week after anthesis. Thus, we calculated the gene expression fold changes of M43 relative to HG by dividing their respective expression ratios. By using >2 folds' changes as cutoff, 58 candidate genes were identified (Table 1). Based on respective expression patterns within individual germplasms, these genes can be divided into five groups. Group I includes six genes, which were upregulated from M43 and HG during seed desiccation, and M43 showed a much higher upregulation compared with that of HG. Group II includes five genes, which were downregulated in M43 and HG during seed desiccation, and M43 showed a smaller reduction compared with HG. Group III includes six genes, which were upregulated in M43 but downregulated in HG. One gene (TRINITY_DN21501_c0_g1_i1_1) is involved in steroid biosynthesis. Group IV includes 15 genes, which kept constant in M43 but downregulated in HG. Within this group, four genes are involved in fatty acid or glycan degradation, they are: TRIN-ITY_DN23188_c1_g1_i1_1, TRINITY_DN28894_c0_g2_i1_1, TRINITY_DN29729_c0_g2_i1_1 and TRINITY_DN26102_c2_g4_i1_1; and two genes are involved in steroid biosynthesis, including TRINITY_DN23273_c1_g1_i1_1 and TRINITY_DN27207_c0_g1_i3_1 (Table 1). Finally, group V includes 26 genes, which were upregulated in M43, but kept constant in HG. Five of them are involved in steroid biosynthesis, including: TRINITY_DN17413_c0_g1_i1_2, TRIN-ITY_DN18365_c0_g2_i1_1, TRINITY_DN26549_c2_g1_i1_1, TRINITY_DN28994_c1_g1_i1_1, and TRINITY_DN4285_c0_g1_i1_2 (Table 1).

qRT-PCR Validated RNA-seq Results
To confirm the RNA-seq results, at least one DEG was selected from each group in Table 1. A total of eight DEGs (Supplementary Table S14) were selected and subjected to qRT-PCR analysis, then the results were compared with the RNA-seq data. We found that both datasets showed generally similar changing trends (Figure 7).

R PEER REVIEW
14 of 22

qRT-PCR Validated RNA-seq Results
To confirm the RNA-seq results, at least one DEG was selected from each group in Table 1. A total of eight DEGs (Supplementary Table S14) were selected and subjected to qRT-PCR analysis, then the results were compared with the RNA-seq data. We found that both datasets showed generally similar changing trends (Figure 7).

C. meiocarpa and C. oleifera Accumulate High Levels of Oleic Acid
Both C. meiocarpa and C. oleifera accumulated 65% to 72% oleic acid (Figure 3), it is interesting to understand why Camellia seed kernel preferably accumulates oleic acid to such high levels. In hickory seeds, a high level of SAD with a low level of FAD2 expression was associated with high levels of oleic acid accumulation [34]. In olive, the expression of FAD2 is repressed by a newly evolved small RNA, and contributes to its high-level oleic acid content [35]. Zeng et al. (2014) showed that in Camellia seeds, fructose-1,6-bisphosphate aldolase (CoFBA) and CoSAD mRNA levels were well-correlated with oil content, whereas CoFAD2 gene expression levels were correlated with fatty acid composition [12]. The authors propose that CoFBA and CoSAD are two important factors to determine Camellia seed oil content, while CoFAD2 gene expression levels were correlated with fatty acid composition. In this study, we found that during seed maturation and desiccation, the coordinated changes in enzyme activities of KAS II, SAD, FatB, FAD2, FAD3, and FAE were associated with seed oil accumulation rates (Figures 1 and 6). Correlation analysis demonstrated that seed oil contents were significantly and positively correlated with 18:0, 18:1 Δ9 , and 18:3 (Figure 5b-d). Notably, the 18C fatty acids were eight-folds higher than that of 16C fatty acids (Figure 4, lower panel), suggesting that KAS II, a key enzyme to CoGAPDH was used as reference gene for qRT-PCR, the relative expression level of H-2 sample was set as 1.

C. meiocarpa and C. oleifera Accumulate High Levels of Oleic Acid
Both C. meiocarpa and C. oleifera accumulated 65% to 72% oleic acid (Figure 3), it is interesting to understand why Camellia seed kernel preferably accumulates oleic acid to such high levels. In hickory seeds, a high level of SAD with a low level of FAD2 expression was associated with high levels of oleic acid accumulation [34]. In olive, the expression of FAD2 is repressed by a newly evolved small RNA, and contributes to its high-level oleic acid content [35]. Zeng et al. (2014) showed that in Camellia seeds, fructose-1,6-bisphosphate aldolase (CoFBA) and CoSAD mRNA levels were well-correlated with oil content, whereas CoFAD2 gene expression levels were correlated with fatty acid composition [12]. The authors propose that CoFBA and CoSAD are two important factors to determine Camellia seed oil content, while CoFAD2 gene expression levels were correlated with fatty acid composition. In this study, we found that during seed maturation and desiccation, the coordinated changes in enzyme activities of KAS II, SAD, FatB, FAD2, FAD3, and FAE were associated with seed oil accumulation rates (Figures 1 and 6). Correlation analysis demonstrated that seed oil contents were significantly and positively correlated with 18:0, 18:1 ∆9 , and 18:3 (Figure 5b-d). Notably, the 18C fatty acids were eight-folds higher than that of 16C fatty acids (Figure 4, lower panel), suggesting that KAS II, a key enzyme to convert 16C to 18C fatty acids, would play critical roles for high oleic acid accumulation in both C. oleifera and C. meiocarpa. Besides KAS II, FatA, SAD, and FAD2 are located on the critical positions of the fatty acid synthesis pathway; thus, they are expected to play important roles for high oleic acid accumulation. Consistent with this notion, Gong et al. (2020) reported that SAD and FAD2 showed high levels of expression during seed maturation [14]. In Arabidopsis, several plastidial stearoyl-ACP desaturase (SAD) paralogs, including FAB2, AAD1, AAD5, and AAD6, work redundantly to catalyze 18:0-ACP into 18:1-ACP during seed maturation [36][37][38]. This could also be the case in Camellia seeds, considering that 18:3 only accounts for a minor fraction of seed oil (Figure 3), thereby suggesting that FAD3 enzymes are either not active or their genes are expressed at very low levels, thus restricting 18:1 ∆9 from being converted into 18:2 and contributing to high oleic acid accumulation. FAD6, FAD7, and FAD8 are plastidial galactolipid desaturases, which sequentially desaturate 18:1 (16:1) to 18:2 (16:2) to 18:3 (16:3) [39][40][41]. Camellia are 18:3 plants, their leaf galactolipids are dominated by 18:3 fatty acids [42,43], suggesting that FAD6, FAD7, and FAD8 play important roles for galactolipid desaturation in Camellia leaf tissue. However, Camellia seed oil only contains trace amounts of 18:3, while 16:2 or 16:3 fatty acids are barely detectable (Figure 2 and Figure S3). These data suggest that FAD6, FAD7, and FAD8 do not make significant contributions to Camellia seed oil synthesis.

The Exocarp Photosynthesis May Contribute to Camellia Seed Oil Accumulation
During seed maturation and desiccation, HG exhibited distinct seed oil accumulation patterns that differentiate itself from QG and M43 (Figure 1). In view of its roles in the trophic connection between seeds and the mother plant, green exocarp could possess photosynthetic capacity and contribute to seed oil accumulation. In chloroembryos such as oilseed rape, seed photosynthesis plays a role in the accumulation of storage lipids [44,45]. The Camellia embryo is embedded in fruit peels, including exocarp, mesocarp, and endocarp. With seed maturation, the endocarp lignifies into a dark seed coat. The seed coat in combination with mesocarp would prevent light penetration into the embryo (Supplementary Figure S2). Thus, unlike chloroembryos, the photosynthesis from Camellia embryo would not make any significant contributions to seed oil accumulation. However, the fruit exocarp contains chlorophyll (Supplementary Figure S2) and stomata [46], and has the potential to assimilate atmospheric CO 2 by photosynthesis to fuel Camellia seed oil accumulation. In accordance with this notion, it was reported that tomato fruit photosynthesis is restricted to the green phases of development, and contributes to net sugar accumulation and growth of the fruits [47,48]. QG and M43 showed green exocarp during seed maturation and desiccation; in contrast, the exocarp of HG turned red in color at the early phase of seed desiccation, which may reduce its photosynthetic capacity of the exocarp. In accordance with the fruit color changes, the seed oil accumulation from HG was much slower compared with QG or M43 (Figure 1). These observations suggest that the exocarp photosynthesis could be an important factor that affects seed oil accumulation in Camellia seeds.

The Harvest Time of Camellia Seeds Is Important to Obtain Higher Oil Yield
In this study, we found a pronounced seed oil loss from QG and M43 at the late phase of seed desiccation (Figure 1), and most individual fatty acid contents also reduced concurrently ( Figure 2). In contrast, the total oil contents from HG was not reduced, instead, 18:1 ∆9 content continually increased at the expense of 16:0 and 18:2 reductions (Figure 2, upper middle panel). Although QG and M43 accumulated much higher seed oils than that of HG between the 51st and the 53rd week after anthesis, these three germplasms contained similar levels of seed oil (~43% of seed dry mass) when their seeds were harvested on the 54th week ( Figure 1). Thus, the seed oil accumulation from individual germplasms in the Theaceae family is diverse and highly dynamic, and harvesting the seeds at the right time is therefore important to obtain high oil yield.

The Potential Pathways Governing Triacylglycerol Degradation at the Late Phase of Seed Desiccation
Previous studies demonstrated that triacylglycerol synthesis and degradation are two concurrent processes during oil seed filling [19,49]. During seed desiccation, the Camellia seed endocarp is lignified into a seed coat, which could prevent nutrients from freely diffusing from the mother plant to the seed kernel. However, some biosynthetic processes in the seed could still persist as other oil seeds. For example, in B. napus, the enzymes of β-oxidation, the glyoxylate cycle, and phosphoenolpyruvate carboxykinase were present in embryos during oil accumulation, and increased in activities and abundance with seed maturation and desiccation [10,17,50]. In this study, we performed transcriptomic profiling at the early and the late phases of seed desiccation from HG and M43, and identified four genes in fatty acid or glycan degradation pathway (Table 1). These data suggest that the transcriptional upregulation of fatty acid degradation genes could play an important role in governing seed oil degradation in M43. In Arabidopsis and rape, triacylglycerol lipases (SDP1 and SFAR) have been demonstrated to be involved in TAG degradation during seed maturation [51,52]. Surprisingly, these lipase homologs were not present from the list in Table 1. We speculate that there are two potential explanations: (1) our selection criteria for DEGs could be too harsh such that some candidate genes could be missed out and/or (2) Arabidopsis or rape seeds are type I seeds, which are small and tolerant to desiccation. In contrast, Camellia seeds are type II seeds, which are large and only partially tolerant to desiccation [53,54]. It remains an open question whether these two types of oil seeds could leverage different lipase for TAG degradation during seed development. Chia et al. (2005) demonstrated that seed oil degradation in B. napus could be fueled into seed protein synthesis [19]. In M43, eight genes in terpenoid (steroid/sesquiterpenoid/ triterpenoid) synthesis pathway were upregulated (Table 1). This raises the possibility that seed oil degradation in M43 could be channeled into terpenoid synthesis. It has been reported that during Camellia seed maturation, seed saponin continuously increased to more than 20% of its seed dry mass [20,55,56]. Our data provide potential molecular evidences that TAG degradation could be associated with terpenoid accumulation at the late phase of seed desiccation.

The Physiological Implications of Camellia TAG Degradation at the Late Phase of Seed Desiccation
As we discussed above, Camellia seeds are only partially tolerant to desiccation [53,54]; thus, mature Camellia seeds have a relatively short time window to germinate under natural drought conditions. We speculate that the oil degradation at the late phase of seed desiccation could make Camellia seeds better prepared for seed germination and seedling establishment and thus may confer evolutionary fitness for Camellia species. This notion is in accordance with a previous report that ripe tea seeds germinated quickly and showed a high rate of germination, while the unripe seeds germinated slowly and unevenly, and showed relatively low germination rates [57]. Besides facilitating seed germination, the saponins synthesized at the late phase of seed desiccation could have functions in the deterrence of animal predators or inhibition of pathogen infection, thereby enhancing seed survival and spreading.

Plant Materials and Growth Conditions
Camellia oleifera cv.  Figure S1). C. oleifera and C. meiocarpa fruits show similar developmental timelines: both blossom and fertilize in late October or early November, the fertilized fruits stay in dormancy during the winter, then start expansion and dry mass accumulation the following May. In September, the average seed weight reaches a constant, the water content continues to decline and is accompanied with dry mass increase, and the seed enters the desiccation stage. At early to mid-November, the fruits are harvested for oil production. In this study, the fruits were harvested weekly from 20 September to 14 November 2015, which correspond to the 46th to the 54th week after anthesis. At each sampling week, the harvested fruits were selected based on the following two criteria: (1) the fruits were positioned on the same direction of the tree so that they would likely receive similar amounts of sunlight and (2) the fruits were uniform in size and color appearance within individual germplasms. Thirty fruits were harvested at 10:00 AM from each germplasm, the mesocarp and seed coat were removed and the kernels from the same germplasm were mixed and randomly divided into six parts, with each part containing 20 seeds as one replicate. Three parts were frozen in liquid nitrogen, then stored in −80 • C freezer for total RNA extraction; the other three parts were cut into flakes with a blade, then freeze-dried (Labcoon, KS, USA) for oil content measurement.

Seed Oil Content Analysis
Kernel powder (~10 mg) was methylated in 1.0 mL of methanol containing 5% sulfuric acid (v/v); 75 µg margaric acid was added as internal standard. After heating at 95 • C in a dry bath for 90 min, the fatty acid methyl esters (FAMEs) were extracted into hexane, concentrated under gentle nitrogen stream, then injected into gas chromatography mass spectrometry in a split ratio of 20:1. The FAMEs were separated in RT-2560 column (0.25 mm × 30 m × 0.25 µm, RESTEK, Bellefonte, PA, USA) and detected by flame ionization detector (FID) and quadrupole mass analyzer (MS) (MDGC/GCMS-2010, Shimadzu, Kyoto, Japan), respectively. Helium flow rates for MS and FID were 1.0 mL min −1 . The temperature for the injection port, ion source, and detector was 240 • C, 200 • C, and 220 • C, respectively. The flow rates of hydrogen, nitrogen, and zero air were 40, 30, and 400 mL min −1 , respectively. After 5 min of solvent delay at 150 • C, oven temperature was ramped at 2 • C min −1 to 200 • C, held for 5 min, then ramped to 150 • C for next sample injection. The fatty acids were identified from MS data by Quest NIST14 chemical library, and quantified from FID data by normalizing to the peak area of internal standard. The kernel oil content was calculated by dividing the total FAMEs to kernel dry mass weight; the data were also transformed into mole percent to represent oil fatty acid composition. At each sampling point, triplicates were used for fatty acid analysis.

Total RNA Extraction
Total RNA was isolated from kernel by a modified CTAB method [58]. Kernel was first ground into powder in the presence of liquid nitrogen, 1.5 g of powder was transferred into 50 mL centrifuge tube, and then 10 mL of CTAB solution and 450 µL of β-mercaptoethanol were added, mixed, heated in 65 • C water bath for 20 min, and then centrifuged at 12,000× g for 10 min (4 • C). The supernatant was transferred into a prechilled tube containing equal volumes of chloroform/isopropanol (24:1, v/v). The mixture was incubated in ice-water bath for 10 min, then centrifuged at 12,000× g for 15 min (4 • C). The above extraction steps were repeated, the supernatant was transferred into a clean tube, and half of a volume of prechilled 8 M LiCl solution and 1% β-mercaptoethanol were added, mixed, then stored in −20 • C freezer overnight. The next day, after centrifugation, the pellet was washed twice with prechilled 75% ethanol, then air dried. The dried pellet was dissolved in 300 µL DEPC-treated water and stored in −80 • C freezer. The quality and concentration of the total RNA was analyzed by 1% agarose gel electrophoresis (Supplementary Figure S7) and spectrophotometer (NanoDrop 2000, Thermo Scientific, Pleasanton, CA, USA), respectively.

Library Construction and Transcriptome Sequencing
Based on their oil accumulation patterns, seeds harvested on the 47th and the 52nd week after anthesis from HG and M43 were selected for transcriptomic profiling, at each sampling point, two replicates were used. For library construction, mRNA was purified from 6.0 µg of total RNA by oligo-dT-bound magnetic beads, the sequencing libraries were constructed by using NEBNext Ultra TM RNA Library Prep Kit (Illumina, NEB, San Diego, CA, USA) and following the manufacturer's instructions. The libraries were sequenced on an Illumina Hiseq X platform, 150 bp paired-end reads were generated.

Bioinformatics Analyses
Raw reads were processed using Trimmomatic [59]; the reads containing poly-N and the low-quality reads were removed. After removing adaptor and low-quality sequences, the clean reads were assembled into expressed sequence tag clusters (contigs), then de novo assembled into transcripts by using Trinity in paired-end method [60].
The function of the unigenes was annotated by alignment with a database (NCBI non-redundant (NR), KOG, GO, SwissProt, eggNOG, and KEGG) by Blastx, the threshold E-value was set at 10 −5 [61,62]. The proteins with the highest hits to the unigenes were used to assign functional annotations.
Total reads per kilobases per million reads (FPKM) [63] and read count values of each unigene were calculated using bowtie 2 [64] and eXpress [65]. DEGs were identified using the DESeq [66] functions estimate size factors and nbinom test. p value < 0.05 and fold change > 2 or fold change < 0.5 were set as the threshold for significantly differential expression. The comparative transcriptome analysis was performed using a workflow modified from a previous study [67]. Briefly, the protein sequence of M43 and HG were blast to each other to identify homologous genes by using reciprocal best BLAST hit (RBH) method. Gene ontology (GO) classification was performed by the mapping relation between SwissProt and GO term. The unigenes were mapped to the Kyoto Encyclopedia of Genes and Genomes database (KEGG) to annotate their potential metabolic pathways.

Quantitative RT-PCR Analysis
qRT-PCR was applied to validate RNA-seq results. The names, accession numbers, and primer sequences of the selected eight genes were provided in Supplementary Table S14. Two micrograms of total RNA were reverse-transcribed by M-uLV reverse transcriptase to obtain cDNA, then qPCR was performed in a 20 µL volume by using a SYBR Premix Es Taq kit (with Tli RNase H) (Mei5 Biotechnology, Beijing, China). The PCR program was: initial denaturation was at 95 • C for 30 s, followed by 40 cycles of 95 • C for 15 s and 55 • C for 15 s; the PCR was finished by 1 cycle of 95 • C for 60 s, 55 • C for 30 s, and 72 • C for 30 s. CoGAPDH was used as the reference gene. Three biological replicates and three technical replicates were performed. Fold difference was calculated using 2 −∆∆Ct method [68].

Correlation Analysis
Multiple regression analysis was performed by SPSS (V17.0; SPSS, IBM, Armonk, NY, USA). The significance of correlations between different parameters were determined by bivariate correlations based on Pearson's correlation (two-tailed).

Conclusions
During seed desiccation, the oil accumulation in C. oleifera and C. meiocarpa showed diverse trends; the fatty acid compositions were continuously modified with the increased ratios of unsaturated fatty acids to saturated fatty acids and monounsaturated fatty acids to polyunsaturated fatty acids. At the late phase of seed desiccation, C. oleifera cv. Min 43 and C. meiocarpa var. Qingguo lost more than 25% of their seed oil reserves, while such an event was not observed in C. meiocarpa var. Hongguo. Transcriptomic profiling between C. oleifera cv. Min 43 and C. meiocarpa var. Hongguo suggests that partial oils in M43 were degraded and could be fed into terpenoid synthesis at the late phase of seed desiccation.

Supplementary Materials:
The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/plants12142591/s1, Figure S1: Photographs of Camellia tree with fruits on 26 September 2015; Figure S2: Photographs of Camellia fruits and seeds harvested on 26 September 2015; Figure S3: GC-MS chromatography of fatty acid methyl esters (FAMEs) derived from Camellia seeds; Figure S4: Correlation matrix (Pearson's R 2 values) between different fatty acid components of Camellia oil. White to blue color represents negative correlation while white to red color represents positive correlation. Values labeled with asterisk (*) are significant at p < 0.05; Figure S5: GO term enrichment analysis of downregulated DEGs (a) and upregulated DEGs (b) from C. meiocarpa var. Hongguo. Figure S6: GO term enrichment analysis of downregulated DEGs (a) and upregulated DEGs (b) from C. oleifera cv. Min 43. Figure S7: Agarose electrophoresis of total RNAs isolated from seed kernels of C. meiocarpa var. Hongguo and C. oleifera cv. Min 43 on the 47th week (26 September) and the 52nd week (1 November) after anthesis; Table S1: Fatty acid content changes during seed desiccation in C. meiocarpa var. Qingguo, C. meiocarpa var. Hongguo, and C. oleifera cv. Min 43. Table S2: Unigene list identified from C. meiocarpa var. Hongguo; Table S3: CDS sequences identified from C. meiocarpa var. Hongguo by blast; Table S4: CDS sequences identified from C. meiocarpa var. Hongguo by EST scan; Table S5: Unigene list identified from C. oleifera cv. Min 43; Table S6: CDS sequences identified from C. oleifera cv. Min 43 by blast; Table S7: CDS sequences identified from C. oleifera cv. Min 43 by EST scan; Table S8: RBH gene list from C. meiocarpa var. Hongguo; Table S9: RBH gene list from C. oleifera cv. Min 43; Table S10: Identified DEGs list from C. meiocarpa var. Hongguo; Table S11: Identified DEGs list from C. oleifera cv. Min 43; Table S12: The list of unigenes in Camellia seed lipid metabolic pathway; Table S13: Differentially expressed lipid genes between C. meiocarpa var. Hongguo and C. oleifera cv. Min 43; Table S14: Selected genes and their primer sequences for qRT-PCR analysis.