Comparative de novo Transcriptome Analysis of Two Cultivars With Contrasting Content of Oil and Fatty Acids During Kernel Development in Torreya grandis

Vegetable oil is an indispensable nutritional resource for human health and mainly characterized by the composition and content of fatty acids (FAs). As a commercial species of gymnosperm, Torreya grandis produces oil-rich nuts with high unsaturated fatty acids content in the mature kernels. In this study, two cultivars, T. grandis ‘Xifei’ and T. grandis ‘Dielsii,’ with distinct oil content were employed to compare the profiles of FAs accumulation during kernel development. The accumulation rate of oil content was significantly different between ‘Xifei’ and ‘Dielsii.’ Besides, the final oil content of ‘Xifei’ (52.87%) was significantly higher than that of ‘Dielsii’ (41.62%) at maturity. The significant differences in main FAs were observed at almost each kernel development stages between the two cultivars. C16:0, C18:1, and C20:3 FA exhibited different accumulation patterns between cultivars. The content and the initiation of accumulation of C20:3 FA were different between the two cultivars. To explore the molecular mechanism associated with different content of oil and FAs between two cultivars, de novo transcriptome of kernels was compared between ‘Xifei’ (high oil) and ‘Dielsii’ (low oil) at three stages of oil accumulation, respectively. Totally 142,213 unigenes were assembled and 16,379 unigenes with a length of over 1,000 nt were successfully annotated, including 139 unigenes related to FA biosynthesis, elongation, and metabolism. Compared with ‘Dielsii,’ totally 1,476, 2,140, and 1,145 differentially expressed genes (DEGs) were upregulated in ‘Xifei’ at the stage of the initiative, the rapid rise, and the stationary oil accumulation, respectively; the number of downregulated DEGs reached 913, 1,245, and 904, respectively. Relative expressions of 11 DEGs involved in FAs biosynthesis and metabolism were confirmed by RT-qPCR. Abundant differentially expressed transcription factors and pathway DEGs were correlated to oil and FAs according to Pearson’s correlation analysis between transcriptome and metabolites (oil and FAs), suggesting their contributions to the differential oil and FAs between the two cultivars during kernel development of T. grandis. To conclude, our findings can provide novel insights into the developmental differences in metabolites and de novo transcriptome correlated to lipid accumulation and FA synthesis of kernels between cultivars with contrasting oil deposits and demystify the regulatory mechanism of high oil accumulation in T. grandis.

Vegetable oil is an indispensable nutritional resource for human health and mainly characterized by the composition and content of fatty acids (FAs). As a commercial species of gymnosperm, Torreya grandis produces oil-rich nuts with high unsaturated fatty acids content in the mature kernels. In this study, two cultivars, T. grandis 'Xifei' and T. grandis 'Dielsii,' with distinct oil content were employed to compare the profiles of FAs accumulation during kernel development. The accumulation rate of oil content was significantly different between 'Xifei' and 'Dielsii.' Besides, the final oil content of 'Xifei' (52.87%) was significantly higher than that of 'Dielsii' (41.62%) at maturity. The significant differences in main FAs were observed at almost each kernel development stages between the two cultivars. C 16:0 , C 18:1 , and C 20:3 FA exhibited different accumulation patterns between cultivars. The content and the initiation of accumulation of C 20:3 FA were different between the two cultivars. To explore the molecular mechanism associated with different content of oil and FAs between two cultivars, de novo transcriptome of kernels was compared between 'Xifei' (high oil) and 'Dielsii' (low oil) at three stages of oil accumulation, respectively. Totally 142,213 unigenes were assembled and 16,379 unigenes with a length of over 1,000 nt were successfully annotated, including 139 unigenes related to FA biosynthesis, elongation, and metabolism. Compared with 'Dielsii,' totally 1,476, 2,140, and 1,145 differentially expressed genes (DEGs) were upregulated in 'Xifei' at the stage of the initiative, the rapid rise, and the stationary oil accumulation, respectively; the number of downregulated DEGs reached 913, 1,245, and 904, respectively. Relative expressions of 11 DEGs involved in FAs biosynthesis and metabolism were confirmed by RT-qPCR. Abundant differentially expressed transcription factors and pathway DEGs were correlated to oil and FAs according to Pearson's correlation analysis between transcriptome and metabolites (oil and FAs), suggesting their contributions to the differential oil and FAs between the two cultivars during kernel development of T. grandis. To conclude, our

INTRODUCTION
Vegetable oil is indispensable and consumed in the human diet that requires high-quality edible lipids mainly characterized by the composition and proportion of saturated fatty acids (SFAs) and unsaturated fatty acids (UFAs). In most cases, fatty acids (FAs) have evolved during the de novo synthesis, elongation, and desaturation process. Most nature oil from plants frequently contains polyunsaturated fatty acids (PUFAs), which contain a long chain of at least 18 carbon atoms and 1 methylene (Lisa et al., 2007) or polymethylene interrupted double bonds (UPIFAs). Sciadonic acid (SA) is a cis5, 11, 14-eicosatrienoic acid ( 5-UPIFAs) with the first unsaturation site at the fifth carbon atom that differs from the structure of other PUFAs (Endo et al., 2009;Xie et al., 2016). Due to the beneficial function for health such as regulating the blood pressure (Sugano et al., 1994;Asset et al., 1999;Endo et al., 2006), anti-inflammatory (Chuang et al., 2009;Chen et al., 2015), and prophylactic impacts on angiosclerosis and obesity (Meng et al., 2020), the plant oils with a certain amount of 5-UPIFAs such as SA are essential for dietary modification of human to stem some chronic illnesses.
Torreya grandis is a gymnosperm species of six members in Torreya genus belonging to the Taxaceae family with wide distribution in South China (Shi et al., 2018). T. grandis produces nut-seeds with approximately 50% of dry weight (DW) highquality oils. Due to dioecious flowers, T. grandis possesses abundant variation of seedlings, from which a considerable number of landraces and cultivars were developed and identified with distinct SSR fingerprints (Chen et al., 2020a,b;Meng et al., 2020;Zhang et al., 2021). The commercial value of the rare conifer tree derives mainly from the oil-rich nuts with high UFA content in the mature kernels (Ni et al., 2015;Chen and Jin, 2019a). Great variation in seeds phenotype and qualities of mature nuts exist extensively between cultivars in T. grandis (Ni and Shi, 2014;Chen and Jin, 2019b). The oil content of T. grandis kernels was reported to range from 51.52 to 55.64% depending on the certain cultivar applied in the test (Ni and Shi, 2014). The content of FAs in mature kernels was significantly different between landraces or cultivars in T. grandis (Shi et al., 2018). Oleic (C 18:1 ) and linoleic acids (C 18:2 ) were dominant FAs in mature kernels based on 10 landraces of T. grandis, ranging from 17.62 to 35.11% and 39.77 to 46.06% of total FAs, respectively. The beneficial SA ranged from 9.18 to 18.15% of the total FAs in tested T. grandis oil. The difference in the oil content and FA composition influences the flavor of nuts in T. grandis.
Underlying the final difference of oil deposition is a complex transcriptional network. Genes that encode enzymes related to de novo lipid synthesis pathways, including acetyl-CoA carboxylase, ketoacyl-ACP synthases (KAS), stearoyl-ACP desaturase, acyl-ACP thioesterases (FATA and FATB), and acyl-CoA synthesis (LACS) are characterized and cloned in oil crops such as Arabidopsis, Brassica napus, and Zea mays. In addition, acyltransferases like glycerol-3-phosphate acyltransferase (GPAT), lysophosphatidic acid acyltransferase (LPAT), diacylglycerol acyltransferase (DGAT), or phospholipid: diacylglycerol acyltransferase are responsible for assembly of triacylglycerol (TAG). The oil content could be raised in LPAT overexpression transgenic Arabidopsis and B. napus (Rao and Hildebrand, 2009). Research showed that DGAT was dominantly responsible for TAG accumulation (Shiu-Cheung and Weselake, 2006;Meng et al., 2009), oil, and oleic acid content (Zheng et al., 2008). Although these pathway genes contributed directly to the distinct biosynthesis of lipids between species and cultivars, little was known about the profiles of regulation of their expression.
The further identification of the transcription factors (TFs) such as LEAFY COTYLEDON1 (LEC1), WRINKLED1 (WRI1), FUSCA3 (FUS3), APETALA2 (AP2)/ETHYLENE RESPONSE FACTOR family, and ABSCISIC ACID3 (ABI3) promote a deeper understanding of oil synthesis during seed development. These TFs were well known in seed storage reserve of nutrients such as lipids, protein, and starch, exhibiting distinct timing patterns, and differentiation of responsibilities that trigger the oil pathway genes (Ruuska et al., 2002). WRI1 was reported to activate oil accumulation by facilitating the carbon flux by glycolysis (Focks and Benning, 1998;Cernac and Benning, 2004). A homeobox TF GLABRA2 was discovered to prevent the accumulation of seed oil in Arabidopsis (Shen et al., 2006). ABI3, AtbZIP10, and AtbZIP25 were found for the activation of transcription of storage protein (Lara et al., 2003). In an attempt at developmental control of storage product accumulation, LEC1 and LEC2 were found expressed early in transgenic Arabidopsis and subsequently induced FUS3, ABI3, and ABI5 expression (Wang et al., 2007), indicating an upstream role of regulation on FUS3 in oil deposition aided by related increase of FUS3 transcripts and transcripts of enzyme genes involved in oil syntheses during seed development.
Recently, a few studies reported differences in FA composition between landraces and amounts of differentially expressed genes (DEGs) between fruits and vegetable samples (leaf, root, and stem; Wu et al., 2018;Ding et al., 2020), indicating a final profile that coordinates oil deposition at mature kernels in T. grandis. Although profiles of oil and FAs are clear at mature kernels, little information is available on dynamic patterns and molecular mechanisms leading to a final difference of oil and FAs with kernel development between cultivars in T. grandis. Torreya grandis 'Xifei' is a representative cultivar with high oil content and is applied extensively in main producing areas for long periods (He et al., 2016), while T. grandis 'Dielsii' is well known for rootstock in the propagation of T. grandis. A significant difference in oil accumulation was observed in mature kernels between the two cultivars in practice. Due to the large and high heterozygous genome of gymnosperm, a complete genome sequence of T. grandis has not been reported up to now. In this study, integrative strategies of lipid measurement and de novo transcriptome were employed to compare and understand the dynamic variation of oil and FA biosynthesis and underlying molecular mechanism during kernel development between T. grandis 'Xifei' and T. grandis 'Dielsii.'

Samples Collection
Nut seeds with different developmental stages were collected from three individual adult trees of 'Xifei' and 'Dielsii, ' respectively. Sampling was performed on the 408, 426,438,448,458,468,475,485,498,510, and 521 day after pollination (DAP). Nut seeds were harvested and immediately preserved at -80 • C before RNA preparation, and the corresponding samples for lipid extraction were stored at -20 • C in darkness until needed.

Lipid Extraction and Determine
Frozen-dried powders were individually incubated in petroleum ether at 50 • C for 8 h to extract the lipids of kernels harvested with development of the two cultivars. The oil content was shown as a percent of sample DW. The fatty acid composition of the lipid sample was measured by gas chromatography (Waters H-class UPLC, Agilent Technologies, United States) as described in Wolff et al. (1999a) based on a 37-kinds of fatty acid methyl ester standard mixture. The relative content of FA composition was individually determined by the retention time of the FAmethyl ester standard.

Transcriptome Assembly and Functional Annotation
High-quality RNAs derived from different developmental kernels was prepared and then sequenced on the Illumina HiSeq2500 platform to generate raw reads with reads length of PE125. Following the removal of reads with the adaptor, reads with more than 10% poly-N, and low-quality reads [>50% with Phred quality score (Q) _5 bases] from the raw sequencing data according to in-house Perl scripts (Biomarker Biological Information Technology Co., Ltd., Beijing, China), high-quality clean data were acquired and assembled by Trinity software, in which min_kmer_cov was set 2 and the other parameters are set to default (Grabherr et al., 2011). Briefly, clean reads were assembled into contigs by the de Bruijn graph using k-mer 25 after trials in varying k-mer sizes. Then, mapping of the reads to contigs containing paired-end reads was made to examine contigs of the same transcript and calculate the distance of these contigs. Eventually, the contigs were connected to obtain assembled transcripts unable to be extended on both end. The longest transcripts were considered the unigenes of gene functional annotation.
Gene function annotation was made on the basis of public protein databases: Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences); KOG (clusters of orthologous groups of proteins); and Swiss-Prot (a protein sequence database with manual annotation and review) with NCBI blast 2.2.28 + with e-value 1.0e-5; Pfam (Protein family) via HMMER 3.0 package using e-value 0.01as well as KO (KEGG Ortholog database) with KEGG Automatic Annotation Server using e-value 1.0e-10. According to the annotation for Nr and Pfam databases, the assignment of unigenes was made in Gene Ontology (GO) with Blast2GO 2.5 (Conesa and Gotz, 2008) using e-value 1.0e-6. GO classification was distributed with the Web Gene Ontology Annotation Plot, which included biological process, molecular function, and cellular component.
The calculation method of the fragments per kilobase per million fragments (FPKM) was employed to determine the transcript abundance of each gene (Mortazavi et al., 2008). The DEGs between compared groups were evaluated with the threshold of FC (fold change) >2 and FDR (false discovery rate) <0.01 (Benjamini et al., 2001).

Reverse Transcript Quantitative PCR
Total RNA was individually extracted from kernels of two cultivars sampled at 408, 438, 458, 475, 510, and 521 DAP with a modified method of RNAiso Plus (Takara Bio Inc, Japan), respectively. The cDNA template of each sample for reverse transcription quantitative PCR (qRT-PCR) was prepared by EasyScript One-Step gDNA removal and cDNA synthesis supermix (TRANSGEN). Pathway genes and TFs involved in FA metabolism were chosen for validation using qRT-PCR. Specific primers of qRT-PCR were obtained with Primer Quest 1 based on the assembly sequences in the transcriptome. The Actin gene was used as reference gene, and the primer sequences were described in the Supplementary Table 9. SYBR Premix Ex Taq TM II (Takara) was employed in qRT-PCR with 20 µl volumes on the Roche Light Cycler480 system (Roche Diagnostics, CA, United States). All samples were tested in triplicate. Each reaction was conducted as below: 95 • C for 30 s, followed by 40 cycles at 95 • C for 5 s, and 60 • C for 34 s.

Statistical Analysis
The determinations of metabolites were performed thrice and error bars in presented figures stand for SD. Statistical analysis was performed with SPSS 20.0. Data comparison was made with one-way analysis of variance (ANOVA), and differences among tested metabolites of two cultivars were evaluated with Duncan's multiple-range tests at p < 0.05. Pearson correlation coefficient (r) was calculated for correlation analysis on DEGs and metabolites.

Dynamic Changes in the Oil Content and Fatty Acids Composition Between T. grandis 'Xifei' and T. grandis 'Dielsii' During Kernel Development
Concerning 17 months for ripening duration in kernels, nutfruits derived from T. grandis 'Xifei' and T. grandis 'Dielsii' , respectively, were sampled from 408 to 521 DAP ( Figure 1A) to characterize the temporal variation of the oil content and FA composition during kernel development. Generally, the oil content was raised gradually with kernel development (Figure 1B and Supplementary Table 1). During the initiation phase, no difference in the oil content was observed and the curve line of the oil content almost overlapped from 408 to 438 DAP between the two cultivars. A wider gap occurred in the subsequent 10 days, and a significant difference in the oil content was observed at 448 DAP. Then, the oil content in 'Xifei' raised more rapidly and exceeded that in 'Dielsii.'     (2) c98678.graph_c0 When the fruits get mature at 521 DAP, the oil content is distinct between 'Xifei' and 'Dielsii, ' reaching 52.87 ± 0.19% and 41.62 ± 0.23%, respectively. To summarize, the oil content of 'Xifei' was raised more rapidly and higher than that of 'Dielsii' during kernel development, suggesting a differentially regulatory mechanism underlying the phenotype of oil between two cultivars of T. grandis. Meanwhile, FA composition and content of 'Xifei' and 'Dielsii' were measured and compared during the whole development of kernels (Figure 1C and Supplementary Table 2). Generally, palmitic acid (C 16:0 ) was the dominant SFA in both 'Xifei' and 'Dielsii.' C 16:0 FA of 'Xifei' fell down after a transitory rise before 438 DAP, while C 16:0 FA of 'Dielsii' kept plummeting during the whole development of kernels. C 16:0 FA of 'Xifei' was significantly less than that of 'Dielsii' at each stage of kernel development. C 18:0 FA of two cultivars was less than 4% during the whole development. UFAs recognized and studied by the current research included oleic acid (C 18:1 ), linoleic acid (C 18:2 9 , 12 and C 18:3 9 , 12 , 15 ), eicosadienoic acid (C 20:2 11 , 14 ), and SA (C 20:3 5 , 11 , 14 ). In general, dominant UFA was C 18:2 9 , 12 , C 18:1 , and C 20:3 5 , 11 , 14 in both cultivars, and total main UFAs were raised with kernel growing. Besides, the content of UFAs in 'Xifei' surpassed that in 'Dielsii' at almost all stages of kernel development and reached 88.27% significantly higher than that in 'Dielsii' (85.88%) at maturation. Among three main UFAs, C 18:2 accounted for the majority of UFA in both cultivars, and the significant differences in C 18:2 content were observed between cultivars at each developmental phase, which ranged from 16.69 ± 2.65% to 47.69 ± 0.7% and 9.10 ± 0.02% to 43.91 ± 0.15% in 'Dielsii' and 'Xifei' during kernel development, respectively. C 18:1 is the secondary UFA. The content of C 18:1 in 'Dielsii' was significantly exceeded that in 'Xifei' during the initiation phase. However, C 18:1 FA of 'Xifei' surpassed that of 'Dielsii' in the subsequent development of kernels. UFAs with 20 carbon atoms mainly include C 20:2 11 , 14 and C 20:3 5 , 11 , 14 in this study. C 20:2 11 , 14 FA contributed approximately 2% in 'Xifei' and 2-4% in 'Dielsii' to total FAs during the maturation of kernels. C 20:3 5 , 11 , 14 FA accounted for 8-14.16% and 2.17-21.76% of total FAs in 'Xifei' and 'Dielsii.' C 20:2 11 , 14 had been accumulated almost from the initiation phase, whereas C 20:3 5 , 11 , 14 was not detected until 475 and 485 DAP in 'Dielsii' and 'Xifei, ' respectively. When the kernels developed into maturation, the content of C 20:3 5 , 11 , 14 FA reached 15.76 ± 0.28% in 'Dielsii' and 10.56 ± 0.08% in 'Xifei.' In conclusion, a significant difference was obtained in C 16:0 , C 18:1 , and C 20:3 5 , 11 , 14 FA between two cultivars during kernel development.

De novo Assembly Sequencing and Annotations
To decipher the molecular mechanism responsible for the difference of oil and FAs between 'Xifei' and 'Dielsii' during kernel development, comparative transcriptome analysis was further performed between two cultivars at 438 DAP (initiation), 475 DAP (rapid rise), and 510 DAP (stationary) in accordance with distinct oil content in lipid measurement, respectively. Over 6.0 Gb of clean data were obtained and Q30 is more than 90% in each sample (Supplementary Table 3). The de novo transcriptome assembly was conducted to obtain 142,213 nonredundant unigenes and 22,976 unigenes with over 1 kb in length (Supplementary Figure 1). Annotation of unigenes was carried out based on the best-matched known sequences. The FIGURE 4 | RT-qPCR of 11 DEGs related to fatty acid biosynthesis and metabolism measured during development. Relative expression of DEGs was shown with means of triplicates. Samples were collected at 408, 438, 458, 475, 510, and 521 days after pollination (DAP). The dark square represents T. grandis 'Xifei.' White square represents T. grandis 'Dielsii.' The significant difference of the DEG expression between development stages were labelled with a to g according to Duncan's multiple comparison test (P < 0.05). best-matched species of the unigenes in T. grandis successively included Picea sitchensis (22%), Amborella trichopoda (7%), Nelumbo nucifera (6%), Vitis vinifera (5%), and Prunus persica (3%; Supplementary Figure 2A).
Functional classification was carried out to better understand the unigenes on biological function and interaction in development. A total of 6,485 unigenes were allocated to 119 pathways, and 139 unigenes were related to fatty acid biosynthesis, elongation, and metabolism pathways. A total of 17,783 unigenes were annotated into 52 subgroups in accordance with GO function terms including biological process, cellular component, and molecular function subgroups (Supplementary Table 4 and Supplementary Figure 2B). The dominating biological processes were metabolic, cellular, and single-organism processes. The top molecular function was catalytic activity and binding.

Differentially Expressed Genes Involved in Fatty Acid Biosynthesis and Metabolism During Kernel Development Between Two Cultivars
A total of 5,683 DEGs were obtained in comparison between cultivars (Supplementary Table 5). Compared to 'Dielsii, ' a total of 1,476, 2,140, and 1,145 DEGs were upregulated during the stage of initiation, rapid rise, and the stationary oil accumulation in 'Xifei, ' respectively. DEGs with significantly decreased expression sum up to 913, 1,246, and 904 at three stages in 'Xifei, ' respectively, (Figure 2). Interestingly, a total of 177 differentially expressed transcription factors (DETFs) were captured between cultivars, which contained 13 DETFs involved in three stages, 41 DETFs activated or suppressed at two stages, and 123 DETFs with functions only at one stage of kernel development.
By comparing development stages, DEGs associated with development were 5,933 and 5,406 in 'Xifei' (Supplementary Table 6) and 'Dielsii' (Supplementary Table 7), respectively. In 'Xifei, ' DEGs of upregulation and downregulation were up to 1,236 and 1,769 when the oil content increased rapidly. A total of 714 and 873 DEGs were increased and decreased in expression when oil was steadily stored, respectively. Similarly, in 'Dielsii, ' the up-and down-regulated DEGs reached 1,637 and 2,065 when oil accumulated rapidly, while 1,316 and 1,012 DEGs were up and down expressed during the stationary phase of oil accumulation. In addition, more DETFs were found in the comparison between stages of kernel development, which were captured by 213 DETFs in 'Xifei' and 206 DETFs in 'Dielsii.' KEGG pathway of DEGs was analyzed and the comparison between cultivars was focused on (Figure 3). Abundant DEGs between cultivars were dominantly involved in metabolism pathways (Supplementary Table 8). DEGs involved in fatty acid metabolism were enriched at both initiation ( Figure 3A) and rapid rise stage ( Figure 3B). DEGs involved in the biosynthesis of UFAs were enriched only at the rapid rise stage. More DEGs were found in the rapid rise phase of oil production between cultivars (Table 1). Besides, up-or down-streamed pathways of fatty acid biosynthesis, such as glycolysis, pyruvate metabolism, glycerolipid metabolism, and glycerophospholipid metabolism, were specifically enriched in the rapid rise stage when compared between cultivars.
Differentially expressed genes expression related to fatty acid biosynthesis and the metabolism was verified in kernel development (Figure 4 and Supplementary Table 9), which included 3 CYP family genes, FAD2-2, FAS1, KAS1, GPAT6, LACS4, FUS3, WRI1, and FATA. The expression patterns of DEGs were almost consistent with the transcriptome.

Integrative Analyses of Fatty Acid and Transcriptome
To further identify the related DEGs contributing to the oil and FAs differences, the analyses of the Pearson Correlation Coefficient (r) were carried out with the total available 5,683, 5,933, and 5,406 non-redundant DEGs identified from the comparison of cultivars and development in 'Xifei' and 'Dielsii, ' respectively, (Supplementary Tables 5-7). A total of 922 and 574 DEGs (p < 0.01) derived from the comparison of the two cultivars (cultivar-DEGs) were closely associated with the synthesis and metabolism of oil and FAs in 'Xifei' and 'Dielsii, ' respectively, (Supplementary Tables 10,11). The numbers of cultivar-DEGs correlated to different types of FA were extremely different between cultivars. Totally 24 and 21 correlated DETFs were identified in 'Xifei' and 'Dielsii, ' respectively, ( Table 2), and the expression profile of DETFs was individually correlated to different FAs according to the FPKM value of transcriptome (Figure 5 and Supplementary Table 12). Concerning distinct oil content between cultivars, WRKY45 (c99504.graph.c0, r = -0.99997, and p = 0.00482) and ERF017 (c99963.graph_c0, r = -0.99996, and p = 0.00599) were involved in oil content of kernels with negative coefficient in 'Xifei, ' exhibiting downregulated expression with kernel development in 'Xifei.' A total of 15 DETFs were correlated to syntheses of unsaturated FAs in 'Xifei, ' containing 5 involved in UFA, 3 DETFs in C 18:1 FA, 4 DETFs in C 18:3 FA, and 3 in C 20:3 FA. Among three DETFs of MYB type termed as GTL1, two exhibited positive correlation to C 18:1 FA (c99828.graph_c0, r = 1.00000, and p = 0.00154) and C 20:3 FA (c99348.graph_c0, r = 0.99998, and p = 0.00413), and one was negatively correlated to UFA (c99816.graph_c0, r = -0.99997, and p = 0.00525). FUS3 (c99613.graph_c0) was involved in C 20:3 FA syntheses in 'Xifei' with negative coefficient (r = -0.99993, p = 0.00766). Given that C 20:3 FA in 'Dielsii' was significantly different from that in 'Xifei, ' a total of 6 DETFs were concerned in C 20:3 FA syntheses with negative correlation, and only one DETF was positively correlated (ERF RAP2-3, r = 0.99988, and p = 0.009995) and significantly upregulated with the development of kernels in 'Dielsii.' To better understand the difference in regulation on the pathway of FA syntheses, DEGs were identified based on references and paralleled with the pathway of FA biosynthesis. The pathway of FA syntheses was initiated from malony-CoA and subsequently catalyzed by enzymes of FA syntheses. The profiles of DEGs were shown in the form of a heatmap (Figure 6 and Supplementary Table 13). Given that C 16:0 FA was differentially produced between two cultivars, the expression of FAS was upregulated at the initial development of kernels development and responsible for the nascent rise of C 16:0 FA in 'Xifei, ' whereas the declined expression of FAS accounted for C 16:0 FA in 'Dielsii.' When C 18:0 FA was converted into C 18:1 FA, FATA (c100682.graph_c0) was upregulated by kernel development in 'Xifei' and 'Dielsii, ' while FATA (c100688.graph_c0) was upregulated only at nascent syntheses of FA in 'Xifei, ' indicating a different rule accounting for the more accumulation of C18:1 between two cultivars. Given that a total of 10 DETFs were significantly associated with C 20:3 FA syntheses, the expression of these DETFs was different between cultivars during kernel development. Eight DETFs were upregulation at the early stage of kernel development in 'Xifei.' Only one DETF (c99475.graph_c0) was specifically upregulated at kernel maturation in 'Dielsii, ' whereas the expression of c99613.graph_c0 was upregulated at the early stage of 'Dielsii, ' suggesting the differential regulation rules on the deposition of C 20:3 FA between 'Xifei' and 'Dielsii.' Torreya grandis was one of the rare conifer trees in China, and was famous for the high content of oil in dry nuts. Although the final profiles of oil and FA have been reported in a few published reports (Niu et al., 2011;Ni and Shi, 2014;Shi et al., 2018), the rules of differential accumulation with kernel development remained unclear. In the current study, oil deposition was compared between two cultivars, 'Xifei' and 'Dielsii, ' suggesting a dynamic profile of oil and FA biosyntheses during kernel development. The oil content of the two cultivars was within the variant scale in previous reports (Shi et al., 2018). The oil content of 'Xifei' rose more rapidly with kernel growth and was significantly higher than that of 'Dielsii' at mature, suggesting a more efficient mechanism of regulation and biosynthesis of lipids in 'Xifei.' Interestingly, FA syntheses were significantly different between the two cultivars. The content of C 18:1 FA in 'Xifei' had been increased more rapidly than 'Dielsii' after 448 DAP responsible for the higher final deposition of C 18:1 FA in 'Xifei.' FIGURE 6 | Differential expression of unigenes associated with fatty acid biosynthesis with kernel development between T. grandis 'Xifei' (left) and T. grandis 'Dielsii' (right). Enzyme, unigene, and expression profiles are individually indicated at each primary step of biosyntheses. Colorful grids ranging from green to red in the row indicated the FPKM value of each unigene at the different development stages of kernels. Unigenes with significant correlation to C 20:3 FA were indicated with an asterisk (*) beside the sequence ID. C 20:3 FA (especially refers to SA) was reported to exist in almost all coniferophyte families (Wolff et al., 2001). Besides, T. grandis was a rare species that contain high levels of C 20:3 FA (Niu et al., 2011). Although different final deposition of C 20:3 FA has been reported between landraces, the developmental changes of C 20:3 FA were poorly understood. In this study, accumulation of C 20:3 FA was initiated at 468 DAP in 'Dielsii' and 485 DAP in 'Xifei, ' indicating a novel rule of biosyntheses different from other FAs and a delayed deposition in 'Xifei' accounting for the final lower content of C 20:3 FA. It was reported that the C 20:3 FA presence and its content varied between varieties (Ni and Shi, 2014;He et al., 2016;Wang et al., 2016;Shi et al., 2018), which are depending on certain cultivars of T. grandis and promising a potential correlation to the function of PUFA. As is well known for rootstock, 'Dielsii' has a lower content of oil and UFAs but higher content of C 20:3 FA, suggesting a smart strategy that enrichment of more beneficial polyUFAs such as C20:3 at the expanse of reduction in the oil content in non-edible cultivars of T. grandis.

Comprehensive Regulation of Differentially Expressed Genes Involved in Fatty Acid Biosynthesis
The supply of FAs is a factor restricting oil accumulation in embryo development (Bao and Ohlrogge, 1999). Based on the comparison of DEGs between two cultivars, the most enriched KEGG pathways were glycolysis and oxidative phosphorylation, implying a difference in supplies of substances of fatty acids biosynthesis or processing between 'Xifei' and 'Dielsii.' Acyl carrier protein (ACP) is crucial for elongation of acyl chain in de novo FA synthesis in plants. Acyl-ACP thioesterases are required to catalyze acyl-ACP hydrolysis to discharge free FA (Dussert et al., 2013). FatA and FatB are two classes of Acyl-ACP thioesterases (Jones et al., 1995). FatA enzymes preferentially hydrolyzed18:1-ACP whereas FatB usually catalyzed saturated acyl-ACP (Dussert et al., 2013). Temperate oil crops, such as Arabidopsis, are typically accepted to accumulate low saturated FAs in accordance with the low level of FatB transcripts (Troncoso-Ponce et al., 2011). T. grandis is a precious conifer tree with high lipid-rich nutfruits, the unsaturated FAs of T. grandis usually occupy over 85% approximately (Shi et al., 2018). In this study, two DEGs encoding FATA, c100688.graph_c0 in 'Xifei' and c100682.graph_c0 in 'Dielsii, ' differentially expressed with kernel development. FATA (c100682.graph_c0) was upregulated during kernel development, which was in contrast to the lower accumulation of C 18:1 FA in 'Dielsii.' FATA (c100688.graph_c0) in 'Xifei' was upregulated at initiation and then declined with the kernel growing, associating with the nascent rise of C 18:1 FA content in 'Xifei.' Interestingly, the expression of three DETFs correlated to C 18:1 FA was consistent with the presence of C 18:1 FA in 'Xifei' and 'Dielsii' (Figure 6). In this study, two DEGs encoding FATA were expressed differentially during kernel development, revealing a novel regulation mechanism of C 18:1 FA syntheses in T. grandis. C 20:3 FA (SA) was first found in Taxus cuspidate (Japanese yew; Itabashi and Takagi, 1982), mainly distributed in gymnosperm species with wide content variation between family, genus, or even the species (Aitzetmuller, 1995;Wolff et al., 1999a,b;Berger et al., 2002). Concerning the outstanding benefits to health, C 20:3 FA was widely studied in extraction, purification, characteristics, and pharmacodynamics (Tanaka et al., 2001;Berger et al., 2002). Subsequently, functional characterization of candidate genes encoding 5-desaturase in transgenic Arobidopsis was carried out (Sayanova et al., 2007). Recently, an attempt on the exploration of genes involved in the biosyntheses of C 20:3 FA was performed in T. grandis using mature seeds and vegetable organs (leaf, stem, and root; Wu et al., 2018). However, the regulation mechanism of C 20:3 FA biosyntheses during kernel development remains unclear. In this study, the DEGs between two cultivars with contrasting content of C 20:3 FA were explored by de novo transcriptome during the whole maturation of kernels. Combined with correlation analyses, abundant DEGs correlated to C 20:3 FA deposition were captured (Supplementary Tables 10,11), including particularly 10 DETFs with significantly different expressions between the two cultivars and responsible for different content of C 20:3 FA, promising a regulation mechanism in the process of C 20:3 FA biosyntheses.

ONE SENTENCE SUMMARY
Lipid synthesis and de novo transcriptome were compared between two cultivars during kernel development in Torreya grandis.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
CZ, HL, and MZ wrote the manuscript. CZ and MZ analyzed the data. HZ, WD, and CHZ performed the experiments. All authors contributed to the article and approved the submitted version.

FUNDING
This work was financially supported by the Basic Public Welfare Project of Zhejiang Province (LGN19C160005); Economy, Commerce and Technology Bureau of Jingning (2021A18); and Science Technology Department of Zhejiang Province (2021C02066-11).