Interspecific Variation in the Unsaturation Level of Seed Oils Were Associated With the Expression Pattern Shifts of Duplicated Desaturase Genes and the Potential Role of Other Regulatory Genes

Seed oils are of great economic importance both for human consumption and industrial applications. The nutritional quality and industrial value of seed oils are mostly determined by their fatty acid profiles, especially the relative proportions of unsaturated fatty acids. Tree peony seed oils have recently been recognized as novel edible oils enriched in α-linolenic acid (ALA). However, congeneric species, such as Paeonia ostii and P. ludlowii, showed marked variation in the relative proportions of different unsaturated fatty acids. By comparing the dynamics of fatty acid accumulation and the time-course gene expression patterns between P. ostii and P. ludlowii, we identified genes that were differentially expressed between two species in developing seeds, and showed congruent patterns of variation between expression levels and phenotypes. In addition to the well-known desaturase and acyltransferase genes associated with fatty acid desaturation, among them were some genes that were conservatively co-expressed with the desaturation pathway genes across phylogenetically distant ALA-rich species, including Camelina sativa and Perilla frutescens. Go enrichment analysis revealed that these genes were mainly involved in transcriptional regulation, protein post-translational modification and hormone biosynthesis and response, suggesting that the fatty acid synthesis and desaturation pathway might be subject to multiple levels of regulation.


INTRODUCTION
Plant storage oils are an essential part of the human diet, providing a number of nutrients the body needs including essential fats and vitamins (Damude and Kinney, 2008). Plant oils are also utilized as a major source of calories in our food supply and as feedstocks for non-food uses such as soaps and polymers (Dyer et al., 2008). Due to the concern on the depletion of fossil fuels and the environmental problems caused by the use of fossil fuels, plant oils have recently attracted considerable attention as cleaner, renewable, and sustainable sources for biofuel production (Durrett et al., 2008).
Plant seed oils consist principally of triacylglycerols (TAG's) and contain both saturated and unsaturated (mono-and polyunsaturated) fatty acids, with the content and fatty acid (FA) composition varying considerably across plants (Voelker and Kinney, 2001). Variations in the ratios of saturated to unsaturated fatty acids (SFA: UFA), and monounsaturated to polyunsaturated fatty acids (MUFA: PUFA) have been widely detected in different plants (Zhang et al., 2015;Ohlrogge et al., 2018). The level of fatty acid unsaturation of oil can directly affect its value and utilization. Polyunsaturated fatty acids, such as linoleic acid (LA), and α-linolenic acid (ALA), are considered to be critical nutrients for human health (Damude and Kinney, 2008), but they are undesirable while being used as biofuel and chemical feedstocks because of their instability in processing, storage and use (Kodali, 2002). The requirements on different optimal fatty acid compositions for different uses have directly led to the growing interest in identifying genetic factors responsible for the variation in seed oil composition to better understand how the relative proportions of fatty acids are regulated in the seed oils derived from different plants.
A considerable amount of research has been conducted to identify genes potentially associated with the differential accumulation of unsaturated fatty acids in seed oils for optimizing the proportions of fatty acids through biotechnology for different industrial and nutrition goals (Rao et al., 2008;Haslam et al., 2013;Chen et al., 2014). Previous studies, however, have mostly focused on the genes encoding for three desaturases, i.e., stearoyl-ACP desaturase (SAD), fatty acid desaturases 2 (FAD2), and fatty acid desaturases 3 (FAD3) (Thambugala et al., 2013;Rajwade et al., 2014Rajwade et al., , 2016. These three enzymes catalyze stepwise desaturation of fatty acids during oil biosynthesis and determine the relative proportions of three major unsaturated fatty acids, oleic acid (OA), LA and ALA, in seeds. Although a great deal of evidence has shown the conserved roles of SAD, FAD2, and FAD3 across plant species, as well as the positive relations between high levels of fatty acid unsaturation and the upregulation of SAD, FAD2, and FAD3 genes in developing seeds (Ohlrogge and Browse, 1995), it remains elusive how these genes are differentially regulated in different plants so as to collectively control the ratios of SFA: UFA and MUFA: PUFA. Duplicate desaturase genes with altered expression patterns during seed development have also been detected in many plants (Cao et al., 2013;Wang and Liu, 2014;Yurchenko et al., 2014;Han et al., 2017). It is unclear, however, whether the differential transcriptional activities of duplicates are responsible for the differential accumulation of PUFAs in different plants and to what degree the duplicates have diverged within and between species. The question still remains: which genes other than desaturase genes are temporally expressed corresponding to the dynamic synthesis of unsaturated fatty acids in developing seeds, regulating the transcription and catalytic activity of desaturation enzymes and potentially contributing to the variation in unsaturated fatty acids accumulation in different plants. The newly developed methods of genome wide association studies (GWAS) and genome-wide transcription analysis allow for a more efficient discovery of additional genes involved in fatty acid synthesis and desaturation (Li et al., 2013;Branham et al., 2015;Menard et al., 2017;Zhao et al., 2019).
As a traditional ornamental and medicinal plant in China, tree peony has attracted wide attention in the recent past for its relatively high oil content (27%) in seeds and high ratio of unsaturated fatty acids (over 90%) . The oil of tree peony has now been recognized as novel edible oil enriched in omega-3 polyunsaturated fatty acid (Han et al., 2016). Comparative studies of oil content and fatty acid composition have revealed significant variations in the relative proportions of three major unsaturated fatty acids among different Paeonia (section Moutan DC) species. The polyunsaturated fatty acid ALA seems to be the dominant component in the oils of the species from the subsection Vaginatae, while the monounsaturated fatty acid OA is predominant in the species of the subsection Delavayanae (Yu et al., 2016). Genes associated with lipid accumulation and polyunsaturated fatty acid synthesis have been cloned and characterized in some tree peony species (Sun et al., 2018;Xiu et al., 2018;Yin et al., 2018). Gene expression patterns associated with differential accumulation of unsaturated fatty acids have also been examined between P. rockii (subsection Vaginatae) and P. lutea (subsection Delavayanae) (Zhang Q. et al., 2018;Zhang et al., 2019). However, their results showed that the predominant unsaturated fatty acid was ALA instead of OA in the oil of P. lutea, while multiple studies have shown that OA was predominated over LA and ALA in P. ludlowii, another species of the subsection Delavayanae (Han et al., 2016;Yu et al., 2016;Xiu et al., 2018). Differences in the ALA, LA, and OA content proportions are in fact indicative of the degree of unsaturation of different seed oils. Comparing the dynamics of fatty acid accumulation and gene expression between congeneric species with similar oil content but varied fatty acid composition can provide valuable clues for identifying genes potentially involved in the regulation of fatty acid desaturation.
In this study, we characterized the fatty acid profiles and the time-course gene expression patterns of P. ostii and P. ludlowii. Based on a comprehensively comparative analysis of homologous genes in sequence divergence, differential expression patterns and their co-expression relationships associated with phenotypic variation, we identified genes encoding desaturases and acyltransferases that were differentially expressed between two species in developing seeds. Network-based comparative transcriptome analysis of the data from the species of Paeonia and two other ALA-rich plants, Camelina sativa and Perilla frutescens, revealed genes potentially involved in regulating ALA biosynthesis and accumulation in plants.

Plant Materials
P. ostii was grown in Shanghai Chenshan Botanical Garden, Shanghai (31 • 4 52 N, 121 • 10 14 E). P. ludlowii was collected from Linzhi,Tibet (29 • 20 21.9 N,94 • 22 35.8 E). Seeds were harvested every 10 days after flowering (DAF), spanning a total of 120 days. Three samples were collected for each species at each developmental stage. Samples used for FA analysis were dried to a constant weight at 60 • C, while those used for RNA extraction were soaked in RNAlater solution and stored at −40 • C prior to RNA extraction. Young stems and mature leaves were also collected and used for comparative analysis.

Total RNA Extraction, cDNA Library Construction, and High Throughput Sequencing
Based on the dynamic pattern of oil accumulation, five stages of seed development (0, 40, 50, 60, and 100 DAF) were selected for transcriptome sampling. Total RNA's were extracted with TIANGEN RNA Prep Plant kit (Tiangen Biotech Co., Ltd.), following the manufacturer's protocol. The quality and quantity of obtained RNA's were determined with Nanodrop 2000C spectrophotometry (Thermo scientific) and 2100 Bioanalyzer (Agilent Technologies). RNA samples from five developmental stages of three replicates were pooled for each species and used for de novo transcriptome sequencing and assembly, to obtain comprehensive data on expressed gene sequences. cDNA library construction and Illumina sequencing were conducted at Beijing Genomics Institute (BGI) under the Illumina manufacturer's instructions (Illumina). In brief, mRNA was isolated from total RNA using Oligo (dT), and subsequently treated with fragmentation buffer. The short-fragmented mRNA was then utilized as template to generate first-strand cDNA. Followed by second strand cDNA synthesis, end reparation, adaptor connection and PCR amplification, the cDNA libraries were sequenced on an Illumina HiSeq 4000. To establish a temporal map of gene expression for each species, RNA's from the five developmental stages of three replicates were also individually used for cDNA library construction and sequencing. A total of fifteen independent cDNA libraries were generated for each species. The sequence raw data from this study have been submitted to the NCBI Sequence Read Archive (SRA) 1 under the BioProject ID PRJNA602603.

De novo Assembly, Functional Annotation, and Expression Quantification
After quality filtering, de novo transcriptome assembly was performed using Trinity (Grabherr et al., 2011). The software TGICL (Pertea et al., 2003) was then used to cluster the assembled transcripts to generate a single non-redundant unigene set for each species.
The BLASTx program was used to align the assembled unigenes against the NCBI non-redundant protein (Nr) database 2 , the Swiss-Prot protein sequence database 3 and the Arabidopsis protein database at the Arabidopsis Information Resource (TAIR 4 ) with an E-value threshold of 10 −10 . Functional annotation by Gene Ontology (GO 5 ) terms were obtained using the Blast2GO program (Conesa et al., 2005). In addition, we obtained metabolic pathway annotations for each hit by searching against the Kyoto Encyclopedia of Genes and Genomes (KEGG 6 ) pathway database (Ogata et al., 1999) using the BLASTx program with an E-value cut off of 10 −10 .
The clean reads generated from 15 independent cDNA libraries were separately aligned to the de novo assembled transcriptome for each species, respectively. The program RSEM (Li and Dewey, 2011) was used for estimating the expression levels of each unigene, normalizing into Fragments Per Kilobase per Million reads mapped (FPKM). The NOISeq R package was used to screen genes differentially expressed at different stages of seed development (Tarazona et al., 2011). Hierarchical clustering of genes was performed using the hclust function in R v.3.5.3 (Ihaka and Gentleman, 1996). The expression stoichiometries of oil accumulation-related genes were calculated with the methods used in previous studies (Troncoso-Ponce et al., 2011). The coefficient of variation (CV) was calculated to evaluate the expression consistence among growth stages and species.
To validate the reliability of the gene expression data obtained by RNA-Seq, the expression levels of genes of interest were measured by real-time quantitative reverse transcription polymerase chain reaction (RT-qPCR) with the LightCycler 96 system (Roche) using SYBR Green PCR Master Mix (TAKARA). Three independent biological replicates and three technical replicates were performed for each reaction. Based on previous studies (Li S.S. et al., 2015;Liu et al., 2015) and the coefficient of variation of expression level, a total of 5 candidate reference genes (ubiquitin, PUF1639, RPS9, Unigene15510_F/Unigene12802_D, and CL5272.Contig1_F/Unigene12410_D) were evaluated for the normalization of RT-qPCR. Three widely-used programs, GeNorm3.5 (Vandesompele et al., 2002), BestKeeper (Pfaffl et al., 2004), and NormFinder (Andersen et al., 2004), were used to rank the stability of the 5 selected reference genes. CL5272.Contig1_F/Unigene12410_D were consistently ranked as the most stably expressed in different tissues, at different developmental phases, and between two species, and were thus used as the reference genes. Gene-specific primers (Supplementary Table 1) were designed using Primer Premier 5.0 (Lalitha, 2000) based on assembled sequences. Relative expression levels of target genes were calculated using the 2 − Ct method (Livak and Schmittgen, 2001).

Full-Length cDNA Cloning, Sequencing, and Phylogenetic Analysis
To explore the divergence and phylogenetic relationships of the genes of interest, the full-length coding sequences of interested genes were cloned and sequenced. RT-PCR was conducted to clone and validate the sequences of genes with complete open reading frames (ORFs) assembled. Primers (Supplementary Table 1) were designed based on the predicted coding sequences and alignments against the GenBank database. The first-strand cDNA's were synthesized using the PrimeScript RT (Perfect Real Time) Kit (Takara Bio). Purified PCR products were cloned into pMD TM 18-T vector (Takara Bio) and sequenced. For partial coding unigenes, the RACE (Rapid Amplification of cDNA Ends) method was implemented using the SMARTer R RACE Amplification Kit (Takara Bio). Genespecific primers for the 5 -and 3 -RACE reactions were designed based on conserved regions identified from the multiple sequence alignments of different genes. Following cloning and sequencing the overlapping 5 and 3 -RACE fragments, pairs of "end-to-end" primer combinations were designed and used to amplify the full-length cDNA. The amplified products were then purified and cloned into pMD TM 18-T vector for sequencing. Sequence identities were confirmed by BLAST searches against the GenBank database. The gene sequences obtained in this study have been submitted to the GenBank database with the accession codes: MH748779-MH748792. Phylogenetic analyses were carried out using the neighborjoining method implemented in MEGA (Tamura et al., 2007).

Construction and Comparison of Gene Co-expression Networks
To identify genes potentially associated with the varied accumulation of unsaturated fatty acids in different species, network-based comparative analyses were carried out based on time course gene expression data. TransDecoder 7 was used to predict the coding sequences (CDS) and the peptide sequences from the assembled transcripts. Cd-hit was used to cluster contigs sharing highly identical CDS (with identity of over 90% and coverage of over 80%) (Li and Godzik, 2006). Expression values of the members included in a cluster were summed up, and the longest transcript was retained as the representative sequence. The Reciprocal Best Hits (RBH) procedure was applied to identify one to one orthologous pairs between P. ostii and P. ludlowii, which allowed for a detailed cross-species comparison (Movahedi et al., 2012). The NOISeq R package was used to screen interspecifically differentially expressed genes (Tarazona et al., 2011). Gene co-expression networks were constructed using the WGCNA R software package (Langfelder and Horvath, 2008). The matrices of pairwise Pearson correlation coefficients between all pairs of transcripts were generated and transformed into adjacency matrices using the formula: adjacency value = |(1 + correlation)/2| β . β represents the soft-thresholding power for the correlation matrix, which gives greater weight to strongest correlations while maintaining the gene-gene connectivity. The resulting adjacency matrices was then converted to topological overlap (TO) matrices via TOM similarity algorithm, and genes were hierarchically clustered based on TO similarity. The dynamic tree-cutting algorithm was used to cut the hierarchal clustering dendrogram, and modules were defined as branches from the tree cutting. The modules containing genes encoding core fatty acid synthesis enzymes and showing similar expression patterns were selected and merged to generate the oil biosynthesis module. The robustness of the merged module was validated using 5,000 permutation tests (P value < 10 −3 ; Zhan et al., 2015). The hyper-geometric test was performed for validating the selected gene set associated with oil biosynthesis (P value < 10 −6 ). The program Cytoscape was used to visualize the network (Shannon et al., 2003).
To confirm the potential roles of predicted genes in promoting polyunsaturated fatty acid synthesis and accumulation, the expression data of seed development of two other ALA-rich plants, C. sativa and P. frutescens, were integrated into analysis (Kagale et al., 2016;Kim et al., 2016). The PhyloMCL algorithm was used to sort the predicted peptide sequences of different species into groups of putative orthologs on the basis of sequence similarity with default parameters (Zhou et al., 2020). Gene coexpression networks were constructed and analyzed using the aforementioned methods.

Accumulation Dynamics of FA's in Developing Seeds
Paeonia ostii and P. ludlowii exhibited a similar trend of change in oil content, though the total oil contents in mature seeds were a bit different between species (Supplementary Figure 1). At the early stage of seed development (0∼40 DAF), the oil content remained constant and low in both species. A critical transition occurred at around 40 DAF, with a sharp increase in oil content that peaked at 90 DAF and dropped slightly afterward ( Figure 1A). Developmental stage-specific changes were detected not only in oil contents but also in oil compositions. Fatty acid profiling revealed five predominant components of the tree peony seed oil, i.e., stearic acid (SA), palmitic acid (PA), OA, LA, and ALA. The profiles of FA composition did not show significant difference between species at the early stages of seed development. Interspecific changes in FA composition emerged along with the massive accumulation of oil starting from 40 DAF to seed maturation. The major difference lay in the relative proportions of mono-and polyunsaturated fatty acids in the oils of different species. ALA was the most abundant fatty acid in the oil of P. ostii followed by LA and OA, whereas in the oil of P. ludlowii, OA was the predominant fatty acid ( Figure 1A). Compared with that of P. ludlowii, the oil of P. ostii had a higher degree of fatty acid unsaturation. Interestingly, this pattern was reversed in stem and leaf tissues. Comparing the fatty acid profiles of mature seeds with those derived from stem and leaf tissues of P. ostii and P. ludlowii revealed that P. ludlowii contained a remarkably higher proportion of PUFA (LA and ALA) than P. ostii in stems and leaves ( Figure 1B). The results also showed that long chain SFA's including arachidic acid (C20:0) and behenic acid (C22:0) were relatively abundant in stem tissues, leading to a higher ratio of SFA:UFA in stems. By contrast, the high concentration of OA in mature seeds increased the ratios of MUFA:PUFA in seed oils, especially in P. ludlowii. Calculating the product-substrate ratio revealed a relatively lower omega-6 desaturation efficiency in P. ludlowii than in P. ostii, but not for omega-3 desaturation reaction (Supplementary Figure 2).

De novo Transcriptome Assembly and Time-Course Gene Expression Profiling
De novo transcriptome assembly using 73.71 million and 73.63 million clean reads generated 54,655 and 57,665 unigenes in P. ostii and P. ludlowii, respectively, with an average length of 911 and 868 nt, N50 length of 1,447 and 1,328 nt, and overall annotation rate of 63.59 and 59.31% (Supplementary Table 2). The analogous distribution patterns of categorized annotations suggested that the expressed genes of two species were involved in similar functional categories and metabolic pathways (Supplementary Figure 3). The results of BLAST searches against the Arabidopsis Lipid Gene Database 8 revealed 819 and 759 unigenes in P. ostii and P. ludlowii respectively, to be associated with lipid metabolism and signaling.
A total of 30 independent cDNA libraries were generated and sequenced for evaluating temporal gene expression patterns during seed development in P. ostii and P. ludlowii. Over 20 million clean reads were generated for each library. The clean reads were then aligned respectively, to P. ostii and P. ludlowii reference transcriptomes assembled. More than 90% of the reads were successfully aligned to the reference transcriptomes, with the least mapped sample showing a ratio of 86.11% (Supplementary Table 3).
A total of 8,089 and 14,218 unigenes identified in P. ostii and P. ludlowii, respectively, were differentially expressed at different seed development stages. Hierarchical clustering of differentially expressed genes (DEGs) generated ten clusters of DEG with similar expression dynamics in each species (Supplementary Figure 4). The unigenes included in each cluster were listed in Supplementary Tables 4, 5. Of them, the unigenes included in the cluster1, 4, and 5 of P. ostii and those from the cluster1, 4, and 7 of P. ludlowii showed bell-shaped expression patterns in agreement with the seed oil accumulation dynamics in both species. Functional enrichment analyses indicated that these 8 http://aralip.plantbiology.msu.edu/ genes were majorly involved in starch and sucrose metabolism, ALA metabolism, LA metabolism, biosynthesis of secondary metabolites, etc., (Supplementary Table 6).

Cross-Species Comparison of Expression Stoichiometry of Oil Biosynthesis Pathway Genes
The key biological processes involved in oil synthesis have been well studied in plants: FA's are exclusively synthesized in plastids and then the activated FA's are exported to endoplasmic reticulum (ER) for modification and TAG assembly (Bates et al., 2013). Genes encoding the core enzymes involved in FA synthesis, modification and subsequent incorporation into TAG's have been identified, which are conserved among plants (Troncoso-Ponce et al., 2011). It is thus possible to detect the potential rate-limiting steps and enzymes associated with interspecific variation in oil components by calculating and comparing the transcriptional abundance of genes encoding for the core enzymes based on seed-specific genome-wide expression data. The results showed that the genes involved in de novo FA synthesis in plastids displayed similar stoichiometries between species (Figure 2), and maintained strong correlations throughout seed development with r 2 values ranging from 0.74 to 0.96, with an average inter-species expression stoichiometry CV value of 0.37. In contrast to the globally conserved expression stoichiometry of genes involved in de novo FA synthesis, genes associated with FA modification and TAG assembly underwent divergence in transcriptional abundance. Of them, the genes encoding FAD2, FAD3, DGAT (acyl-CoA: diacylglycerol acyltransferase), and PDAT (phospholipid: diacylglycerol acyltransferase) exhibited distinct stoichiometries in different species, with the CV values of stoichiometry of 0.49, 0.69, 0.65, and 0.41, respectively. The correlations of expression stoichiometries of various genes involved in TAG assembly were not statistically significant during oil-accumulation stages ( Figure 2C).

Sequence and Expression Validation of Distinctly Expressed Desaturase and Acyltransferase Genes
Based on assembled transcripts, the complete coding sequences of FAD2 and FAD3 were cloned and sequenced in both P. ostii and P. ludlowii. The expression levels of all desaturase and terminal acyltransferase genes throughout seed development were validated by RT-qPCR. The results revealed the existence of distinct transcripts for all of the genes examined. Sequence comparison and phylogenetic analysis showed that P. ostii and P. ludlowii had at least two FAD2 and three FAD3 paralogous genes, designated respectively, as Po/PlFAD2-1, Po/PlFAD2-2, Po/PlFAD3-1, Po/PlFAD3-2, and Po/PlFAD3-3. The tree peony FAD2 paralogous genes formed two separate groups with those from other plants in the phylogenetic tree, suggesting that they were derived from a duplication event happened in the common ancestor of different plants (Supplementary Figure 5). Instead, the FAD3 paralogs pairs from P. ostii and P. ludlowii clustered together in the phylogenetic tree, forming an independent clade with three paralogous branches (Supplementary Figure 6A). Two lineage-specific duplication events might have taken place subsequently prior to the divergence of tree peony species to produce the three FAD3 paralogous genes. All tree peony FAD3 genes contained three histidine rich motifs. The FAD3-1 branch was distinct from the other two branches for lacking of the isoleucine to valine substitution in the third histidine rich motif (Supplementary Figure 6B).
A significant correlation was found between the results of RT-qPCR and RNA-Seq, as evaluated by linear regression analysis (r 2 = 0.91; p < 0.0001; Supplementary Figure 7). Differential expressions of paralogous genes were detected in both FAD2 and FAD3 genes. As shown in Figure 3A, FAD2-2 was highly expressed in leaves and young seeds, and was down-regulated at seed maturing stages. In contrast, FAD2-1 was exclusively expressed in seeds. It was up regulated during seed maturation, and showed significantly higher expression levels in P. ostii than in P. ludlowii. Similarly, FAD3-1 was mainly expressed in leaves and young seeds, while FAD3-2 and FAD3-3 were mainly expressed in maturing seeds. The expressions of FAD3-2 and FAD3-3 were evidently increased in the oil accumulating stages. The relative expression level of FAD3-3 was much higher than FAD3-2 in both species, and the total level of expression of FAD3-2 and FAD3-3 in P. ostii was significantly higher than that in P. ludlowii.
DGAT's were a group of enzymes that play critical roles in TAG synthesis and accumulation. Previous studies have shown that DGAT's were encoded by multiple genes that evolved independently but became functionally convergent during evolution (Liu et al., 2012). In this study, three analogous DGAT transcripts were identified in each species, Po/PlDGAT1, Po/PlDGAT2, and Po/PlDGAT3, which were homologous to Arabidopsis At2g19450, At3g51520, and At1g48300, respectively. Among them, DGAT3 was mainly expressed during the oilaccumulating development stage, and showed a much higher FIGURE 2 | Cross-species comparison of expression stoichiometry of genes involved in fatty acid synthesis (FAS) and triacylglycerol (TAG) assembly. (A) Relative transcript abundance of genes encoding the core enzymes (indicated in bold) involved in lipid synthesis, with the grid showing the ratios of transcripts of each gene to the sum of FAS and TAG assembly transcripts at different developmental stages (0, 40, 50, 60, and 100 DAF from left to right) in P. ostii (upper) and P. ludlowii (lower), respectively; (B) Correlations of the expression stoichiometries of the 14 enzymes involved in de novo fatty acid synthesis at 50 DAF between P. ostii and P. ludlowii; (C) Correlations of the expression stoichiometries of the 10 enzymes involved in TAG assembly at 50 DAF between two species. Abbreviations are as follows: PDHC, plastidial pyruvate dehydrogenase complex; ACC, acetyl-CoA carboxylase; CoA, coenzyme A; ACP, acyl carrier protein; MCMT, malonyl-CoA ACP transferase; KAS, ketoacyl-ACP synthase; KAR, 3-ketoacyl-ACP reductase; HAD, 3-hydroxyacyl-ACP dehydratase; ENR, 2-enoyl-ACP reductase; FATA, acyl-ACP thioesterase A; FATB, acyl-ACP thioesterase B; SAD, stearoyl-ACP desaturases; LACS, long-chain acyl-CoA synthetase; FAD2, fatty acid desaturase 2; FAD3, fatty acid desaturase 3; GPAT9, glycerol 3-phosphate acyltransferase 9; LPAAT, lysophosphatidic acid acyltransferase; LPCAT, lysophosphatidylcholine acyltransferase; DGAT, acyl-CoA: diacylglycerol acyltransferase; PDAT, phospholipid: diacylglycerol acyltransferase; CPT, CDP-choline: DAG cholinephosphotransferase; PDCT, phosphatidylcholine:diacylglycerol cholinephosphotransferase; G3P, glycerol-3-phosphate; LPA, lysophosphatidic acid; PA, phosphatidic acid; LPC, lysophosphatidylcholine; PC, phosphatidylcholine; DAG, diacylglycerol; TAG, triacylglycerol. expression level in P. ostii than in P. ludlowii. DGAT1 and DGAT2 were expressed in both leaves and seeds with a lower expression level compared to DGAT3. As shown in Figure 3B, the expression levels of DGAT1 and DGAT2 in P. ostii were lower than those in P. ludlowii during the early stages of seed development, but the situation reversed at the later stages. PDAT has been proposed to be involved in TAG biosynthesis via an acyl-CoA independent pathway, enhancing oil accumulation (Bates et al., 2013). The PDAT transcripts homologous to Arabidopsis PDAT1 (At5g13640), PDAT2 (At3g44830), and the PDAT-related gene (At4g19860) respectively, were identified in both P. ostii and P. ludlowii. The expression level of PoPDAT2 in mature seeds was much higher than that of PlPDAT2. PDAT1 was clearly down-regulated with progressing seed development, whereas the PDAT-related gene was continuously up-regulated in P. ostii but not in P. ludlowii.

Network-Based Identification of Genes Potentially Associated With PUFA Biosynthesis and Accumulation
The gene co-expression network was constructed using WGCNA. A total of 23 co-expression modules were generated for P. ostii, consisting of 21,842 genes. Highly correlated modules were merged to generate the oil biosynthesis module, which contained 3,745 genes. Among the 3,745 co-expressed genes, 687 genes were differentially expressed between P. ostii and P. ludlowii, suggesting their probable involvement in the regulation of FA desaturation and promoted ALA accumulation in P. ostii.
To validate the potential regulatory roles of genes co-expressed with the oil biosynthesis pathway genes in P. ostii, the gene co-expression networks of two ALA-rich species C. sativa and P. frutescens were also constructed based on the gene expression profiles during seed development. The co-expression networks comprised 21,294 genes in C. sativa, and 24,754 in P. frutescens. The oil biosynthesis modules of C. sativa and P. frutescens consisted of 4,454 and 2,329 genes, respectively. Comparing the co-expressed gene sets with that of P. ostii revealed 265 orthogroups that were shared by three ALA-rich species. A total of 362 genes co-expressed in P. ostii were also co-expressed in other two species. Of them, 91 were differentially expressed between P. ostii and P. ludlowii (Supplementary Table 7). Go enrichment analysis revealed that these genes were majorly involved in transcriptional regulation (such as genes encoding the oil store regulating genes WRI1, the basic helix-loop-helix proteins bHLH113, and the ethylene responsive factor ERF72), protein modification (such as the protein kinase genes and the ubiquitin-conjugating enzyme gene UBC5), signal transduction (such as JAZ1 involved in the jasmonic acid mediated signaling pathway), and metabolic processes (Supplementary Table 8).
A few genes were found to be conservatively co-expressed with FAD2 and FAD3 across three ALA-rich species (Figure 4 and Supplementary Figure 8). Apart from some genes previously known to regulate desaturase genes, such as FUSCA3 (FUS3) and LEAFY COTYLEDON1-like (L1L), these genes included genes encoding the flavin-containing monooxygenase family protein (YUC10), the auxin-responsive GH3 family protein (DFL1), gibberellin-regulated family protein (GASA6), and genes encoding the proteins from the major facilitator superfamily. Few genes were detected to be conservatively co-expressed with DGAT3 and PDAT2 in three ALA-high species.

DISCUSSION
Given the indispensability of PUFA for human health and the distinct properties required for industrial uses, it is of great necessity to decipher key regulators and possible mechanisms underlying variations in seed oil fatty acid composition (Dyer et al., 2008). P. ostii and P. ludlowii are congeneric species but have diverged in FA accumulation patterns, with ALA and OA being the most abundant fatty acids in P. ostii and P. ludlowii seed oils, respectively. The higher percentage of polyunsaturated fatty acid indicated a greater degree of fatty acid unsaturation in the P. ostii seed oil. However, within stems and leaves, the levels of LA and ALA in P. ostii were lower than in P. ludlowii. It was also found that the abundance of OA remained very low in stems and leaves, which could hardly be detected in both species. These results not only suggested that the differences in seed oil composition between P. ostii and P. ludlowii were unlikely to be plastic changes induced by environmental variation since previous studies have shown increased proportions of unsaturated fatty acids in plants at lower temperatures/higher altitudes (Linder, 2000;Zhang et al., 2015;Guerin et al., 2020), but also indicated that the reactions of fatty acid desaturation were controlled by distinct regulatory mechanisms in different organs.
The relative proportions of unsaturated FA's in seed oils were largely determined by the stepwise desaturation processes in the oil biosynthesis pathway. The expression levels of individual FAD's and their specificity and efficiency in catalyzing fatty acid desaturation play a key role in controlling the component ratios and the unsaturation level of seed oils. The FAD2 and FAD3 genes of tree peony have experienced gene duplication events during evolution. The varied expression patterns of the FAD2 and FAD3 duplicates in different organs clearly showed subfunctionaliztion of FAD2 and FAD3 paralogs following duplication. The significantly higher expressions of FAD2-1, FAD3-2, and FAD3-3 in P. ostii during seed maturing period were consistent with the high proportion of PUFA's in the P. ostii seed oil. Instead, the downregulation of FAD2-1 in the maturing seeds of P. ludlowii probably induced reduction in omega-6 desaturation efficiency, decreasing the conversion rate of OA to LA and consequently leading to the higher accumulation of OA in P. ludlowii.
The last step in TAG biosynthesis was catalyzed by two types of enzymes: DGAT and PDAT. DGAT catalyzed the acyl-CoAdependent biosynthesis of TAG, which has been proposed to be a bottleneck for oil accumulation in some plant species (Cahoon et al., 2007). Different types of DGAT's with similar catalytic properties and no sequence homology have been identified (Liu et al., 2012;Bates et al., 2013). Gene expression profiling revealed that the tree peony DGAT genes homologous to Arabidopsis DGAT1, DGAT2 and DGAT3, respectively, were differentially expressed at different stages of seed development and between two peony species. The expression profiles of different genes shared a common feature that the expression levels of each gene in P. ostii were significantly higher than in P. ludlowii from mid to late stages of seed development . DGAT3 showed the most significant difference between species. Previous studies have shown that DGAT3 was distinct from DGAT1 and DGAT2 for its cytoplasmic localization due to lack of transmembrane domains. Although the soluble DGAT3 has been reported to participate in the cytosolic pathway of TAG synthesis, its specific role in TAG synthesis remains unclear. As demonstrated in peanut, soybean and Arabidopsis (Saha et al., 2006;Hernandez et al., 2012;Turchetto-Zolet et al., 2016), the putative DGAT3 homologs in peony species were ubiquitously expressed, with higher expression levels in maturing seeds. The main difference lay in the most abundant accumulation of the DGAT3-related transcripts occurring at the mid but not the late stages of seed development. The correlations between the expression levels of DGAT3 and the distinct fatty acid profiles observed in different peony species were yet to be clarified. Unlike DGAT's, PDAT catalyzed the acyl-CoA-independent formation of TAG in plants. Previous studies on ALA-rich species, such as Linum usitatissimum and Perilla frutescens, have suggested the activity of PDAT to preferentially transfer ALA into TAG (Pan et al., 2013;Kim et al., 2016). Thus, compared to DGAT, PDAT might contribute more to TAG synthesis in seeds that were high in polyunsaturated fatty acids. The increased expression of peony PDAT2 and the PDAT-related gene in P. ostii seemed to support this issue, though functional and expression divergence of PDAT paralogs had been detected in this study and in other plants (Pan et al., 2013(Pan et al., , 2015Zang et al., 2019).
The congruent patterns of variation at the phenotypic and molecular levels suggested that the relative transcriptional abundances of the enzymes involved in FA desaturation and TAG assembly were correlated with the differential accumulation of unsaturated FA's in different peony species. The higher expression of FAD2-1, FAD3-2, and FAD3-3 possibly promoted sequential FA desaturation in P. ostii, pushing the reaction toward the formation of PUFA. The functional relationships between the raised expression of DGAT2, DGAT3, and PDAT2 and the preferential accumulation of ALA in P. ostii remain to be elucidated. Previous studies have shown the ALA-preference of a few DGAT and PDAT enzymes in other plants (Pan et al., 2013;Kim et al., 2016;Xu et al., 2018). Experimental evaluation of the selectivity of DGAT's and PDAT's in peony species will be helpful to clarify whether the enhanced expression of DGAT2, DGAT3, and PDAT2 have contributed to the enrichment of ALA in P. ostii seed oil by preferentially transferring ALA into TAG. It was evident that the alterations in the relative proportions of the component fatty acids in seed oils and the developmental stage-specific changes in fatty-acid concentrations were not dependent solely on the activity of one single gene encoding enzyme, but on the synergistic effects of all pathway genes. Not only was the mRNA level of each gene encoding enzyme under the transcriptional control during seed development, but there could be some regulators participated in the spatialtemporal co-regulation of different genes and the coordination of upstream and downstream reaction modules. Deciphering the components that coordinately regulated the expression of multiple lipid biosynthetic genes and enzyme activity was essential for dissecting the regulatory mechanisms underlying variations in FA content and composition in seed oils.
Genes with coordinated expression were usually assumed to work cooperatively in related pathways, sharing similar function and regulation, especially when the co-expression relationships were conserved across species (Su et al., 2011;Medema and Osbourn, 2016). Gene co-expression analysis has been increasingly used to prioritize genes that might underlie a phenotype. In this study, we screened co-expressed gene pairs conserved across evolutionarily distant species and identified genes that were conservatively co-expressed with genes known to be involved in FA synthesis and desaturation. The results clearly showed that some of the co-expressed genes were well known as key regulators of seed oil generation, such as the transcription factor WRI1, which was identified to govern the flux of carbon through glycolysis and fatty acid synthesis by regulating genes associated with fatty acid biosynthesis in plastids (Cernac and Benning, 2004;Baud et al., 2007). The higher expression level of WRI1 of P. ostii than that of P. ludlowii was also in accordance with the slightly higher oil content of P. ostii.
In addition to WRI1, a dozen other transcription factors of different families were also tightly co-expressed with the oil biosynthesis pathway genes, including the basic helix-loop-helix proteins (bHLH113), the ethylene responsive factors (ERF72), the MYB transcription factors (MYB61, MYB86) and the NAC transcription factors (NAC100), the B3 protein family (HSL1). Most of them or their homologs have been found to be associated with diverse aspects of lipid metabolism, though very few of them have been shown to modify fatty acid proportions. For instance, bHLH113 has been shown to be up-regulated at different stages of mesocarp development to produce higher oil yield in oil palm and to be co-expressed with major lipid genes in Hiptage benghalensis (Wong et al., 2017;Tian et al., 2019). The role of hormones in regulating seed development and storage compound accumulation has been widely recognized (Yeap et al., 2017;Grimberg et al., 2018). Mutations in hormonerelated transcription factors could also induce changes in the fatty acid profiles of seed oil, showing varied SFA: UFA (Hughes et al., 2008) and MUFA:PUFA ratios (Geilen et al., 2017). The functional relationships of the NAC and MYB transcription factors identified in this study to the high-level accumulation of ALA in seed oils remained largely unknown, without empirical evidences showing how they were connected. Nevertheless, such connections were supported by the genomewide association study of A. thaliana for identifying determinants of variation in seed oil composition, which detected a group of MYB genes (MYB10, MYB64, MYB67, and MYB103) linked to the SNP's associated with variation in fatty acid proportions (Branham et al., 2015).
Some genes involved in protein post-translational modification, such as protein phosphorylation (Protein kinase superfamily protein) and ubiquitination (UBC5), were also tightly co-expressed with fatty acid biosynthesis genes. Previous studies have shown that a serine/threonine/tyrosine protein kinase participated in the regulation of the content and composition of the Arabidopsis seed oil by mediating phosphorylation of oil body proteins (Ramachandiran et al., 2018). It has also been revealed that the seed-expressed casein kinase I enhanced the bHLH-mediated transactivation of the FAD2 gene promoter via phosphorylation of the bHLH transcription factor (Kim et al., 2010), while the ubiquitin proteasome pathway modulated FAD3 protein amounts in response to temperature (O'Quin et al., 2010). The results of this study were thus consistent with previous studies that suggested that post-translational regulatory mechanisms played an important role in modulating FAD2 and FAD3 protein abundance and stability (Martz et al., 2006;Bourassa, 2008).
Network-based comparative transcriptome analysis also revealed genes with diverse functions to be conservatively coexpressed with FAD2 and FAD3 across different species. Some of them (FUS3 and L1L) were previously known to participate in the regulation of FAD3 via (bZIP67:L1L:NF-YC2) transcriptional complex, whose induction was triggered by FUS3 (Mendes et al., 2013). Though bZIP67 was only found to be co-expressed with FAD2 and FAD3 in C. sativa, another group A bZIP member, AREB3, which was closed related to bZIP67 was found to be co-expressed with FAD2 and FAD3 in P. frutescens and P. ostii, and was differentially expressed between P. ostii and P. ludlowii, previous research has also identified that AREB3 was able to enhance FAD3 expression in the presence of L1L and NF-YC2 (Mendes et al., 2013). In P. ostii, NF-YC2 gene was also found to be co-expressed with FAD3 and FAD2 while in P. frutescens and C. sativa, some other NF-YC family members were found to exhibit similar expression patterns with FAD2 and FAD3, suggesting the potential role of other NF-YC family members in regulating desaturase genes. Apart from them, other co-expressed genes have not yet been reported to date to play roles in regulating fatty acid proportions in plants. Of them, YUC10 (Jia et al., 2014), DFL1 (Nemhauser et al., 2006) and GASA6 (Qu et al., 2016) were associated with hormone biosynthesis or response. The function of the major facilitator superfamily proteins remained unclear in plants, though it has been shown that the major facilitator superfamily domaincontaining protein was an omega-3 fatty acid transporter in human (Zhang W. et al., 2018;Zhou et al., 2019). However, their relationships with the variation of fatty acid composition were not only demonstrated by the conserved co-expression patterns in different plants shown in this study, but were also supported by early genome-wide association studies (Menard et al., 2017;Zhao et al., 2019). The roles of these genes in regulating fatty acid composition remained to be functionally characterized in future research.

CONCLUSION
Fatty acid synthesis and desaturation involved a coordinated series of enzymatic reactions. Lineage-specific duplication and subsequent expression alterations of desaturase genes obviously contributed to the phenotypic divergence in the relative proportions of unsaturated fatty acids in the seed oils of peony species. However, the variation in fatty acid proportions was not entirely dependent on desaturase genes. A few genes associated with transcriptional and translational control (such as NAC100, MYB61, and UBC5) and hormonal regulation (such as YUC10, DFL1, and GASA6) seemed to be also involved to regulate the desaturase activity and coordinate various enzymatic steps. These genes were highly correlated with phenotypic variation, and were tightly co-expressed with desaturase genes, providing new targets for future functional characterization and bioengineering to extend our understanding of the regulatory mechanisms underpinning variations in fatty acid composition and unsaturation levels, and to help improve the composition of nutritional and industrial oils.

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
MW performed the wet lab work and data analysis and drafted the manuscript. GL, CZ, and JJ helped to set up experiments. ZX helped the material collection. YW, WZ, and ZS participated in the data analysis. LG, YH, and JY conceived the idea, participated in the design of the study, and finalized the manuscript. All of the authors read and approved the final manuscript.

FUNDING
This study was supported by the Chenshan Key Scientific Research Projects of Shanghai Municipal Administration of Landscaping and City Appearance (F132419) and the grant from Science and Technology Commission of Shanghai Municipality (14DZ2260400).