Fatty acid accumulation analysis based on transcriptome sequencing of developmental herbaceous peony ‘Hangshao’ ( Paeonia lactiflora ‘Hangshao’) seeds

DOI: https://doi.org/10.21203/rs.2.11287/v1

Abstract

Background: Paeonia lactiflora ‘Hangshao’ is widely cultivated in China because its root can be used to produce raw materials for traditional Chinese medicine ‘Radix Paeoniae Alba’. Due to the presence of abundant unsaturated fatty acids in its seed, it also can be regarded as a new oil plant. However, the process of the biosynthesis of unsaturated fatty acid in herbaceous peony ‘Hangshao’ remained largely unknown. Therefore, transcriptome analysis is helpful to better understand the molecular mechanisms. Results: Five main fatty acids, stearic acid, palmitic acid, oleic acid, linoleic acid and α-linolenic acid, were detected, and their absolute contents increased first and then decreased during seed development. A total of 150,156 Unigenes were obtained by transcriptome sequencing, with an average length of 1,030 bp. There were 1,550 Unigenes annotated in the seven functional databases including NR, NT, GO, KOG, KEGG, SwissProt and InterPro. Based on KEGG database, 1,766 Unigenes were annotated in the lipid metabolic pathways, among which 103, 74 and 70 Unigenes are annotated into fatty acid biosynthesis pathway, fatty acid elongation pathway and unsaturated fatty acid synthesis pathway; respectively. A total of 1480 DEGs were detected. Among them, 83 DEGs were enriched in the fatty acid metabolism pathway, including 12 DEGs involved in the fatty acid biosynthesis and 1 DEG involved in fatty acid elongation. Furthermore, qRT-PCR was used to analyze the expression patterns of nine fatty acid biosynthetic related genes including FBCP, BC, FabD, FabF, FATB, KCR, FAD2, FAD3 and FAD7, and it showed that they all highest expressed at 45 DAF. Conclusions: This study provides the first comprehensive genomic resources characterizing herbaceous peony seeds gene expression at the transcriptional level. These data lay the foundation for elucidating the molecular mechanism of the lipid biosynthesis and fatty acid accumulation for herbaceous peony 'Hangshao'.

Abstract

Background: Paeonia lactiflora ‘Hangshao’ is widely cultivated in China because its root can be used to produce raw materials for traditional Chinese medicine ‘Radix Paeoniae Alba’. Due to the presence of abundant unsaturated fatty acids in its seed, it also can be regarded as a new oil plant. However, the process of the biosynthesis of unsaturated fatty acid in herbaceous peony ‘Hangshao’ remained largely unknown. Therefore, transcriptome analysis is helpful to better understand the molecular mechanisms.

Results: Five main fatty acids, stearic acid, palmitic acid, oleic acid, linoleic acid and α-linolenic acid, were detected, and their absolute contents increased first and then decreased during seed development. A total of 150,156 Unigenes were obtained by transcriptome sequencing, with an average length of 1,030 bp. There were 1,550 Unigenes annotated in the seven functional databases including NR, NT, GO, KOG, KEGG, SwissProt and InterPro. Based on KEGG database, 1,766 Unigenes were annotated in the lipid metabolic pathways, among which 103, 74 and 70 Unigenes are annotated into fatty acid biosynthesis pathway, fatty acid elongation pathway and unsaturated fatty acid synthesis pathway; respectively. A total of 1480 DEGs were detected. Among them, 83 DEGs were enriched in the fatty acid metabolism pathway, including 12 DEGs involved in the fatty acid biosynthesis and 1 DEG involved in fatty acid elongation. Furthermore, qRT-PCR was used to analyze the expression patterns of nine fatty acid biosynthetic related genes including FBCP, BC, FabD, FabF, FATB, KCR, FAD2, FAD3 and FAD7, and it showed that they all highest expressed at 45 DAF.

Conclusions: This study provides the first comprehensive genomic resources characterizing herbaceous peony seeds gene expression at the transcriptional level. These data lay the foundation for elucidating the molecular mechanism of the lipid biosynthesis and fatty acid accumulation for herbaceous peony 'Hangshao'.

Keywords

Paeonia lactiflora ‘Hangshao’, fatty acid biosynthesis, transcriptome, differentially expressed gene, qRT-PCR

Background

As a Chinese traditional famous flower, herbaceous peony (Paeonia lactiflora Pall.) is cultivated more than 4000 years. Owing to the beauty of the color, shape and fragrance, it was regarded as ‘the Prime-Minister of flowers’ and was called as ‘two superb of flowers’ with tree peony (Paeonia suffruticosa Andr.) which is regarded as ‘the King of flowers’. Herbaceous peony was cultivated as a medicinal and ornamental plant in China. So far there are more than 500 cultivars distributed in different regions of China [1]. The dried roots of two herbaceous peony species are the main raw materials for traditional Chinese medicine ‘Radix Paeoniae Alba’ and ‘Radix Paeoniae Rubra’ [2, 3], which have the functions of antithrombus [4], antitumor [5], anti-post-traumatic stress disorder [6], alleviating inflammation [7], removing pathogenic heat from the blood, treating blood stasis and relieving pain [8] and so on. As an important medicinal plant, herbaceous peony ‘Hangshao’, which root was the traditional Chinese medicine ‘Radix Paeoniae Alba’, is mainly cultivated in Zhejiang province, and widely cultivated in China [9]. And it is a single-lobe and has higher seed-setting rate. However, the tremendous seeds of herbaceous peony ‘Hangshao’ are mainly wasted every year except that a handful of seeds is used for propagation [10]. Previous studies showed that the seed oil content of ‘Hangshao’ was about 25% with more than 90% unsaturated fatty acids (UFAs) [11]. Therefore, herbaceous peony ‘Hangshao’ has the potential to be cultivated as a multi-functional plant not only as ornamental, medicinal plant, but also as an edible oil resource.

Herbaceous peony ‘Hangshao’ seed oil predominantly contains five fatty acids which are oleic acid (C18:1Δ9, OA), linoleic acid (C18:2Δ9,12, LA), α-linolenic acid (C18:3Δ9,12,15, LA), palmitic acid (C16:0, PA) and stearic acid (C18:0, PA) [10-11]. OA, LA and ALA are classified to UFA which included monounsaturated fatty acid (MUFA) and polyunsaturated fatty acid (PUFA). OA which can improve diastolic heart function is belonged to MUFA [12]. LA and ALA are belonged to PUFA which can reduce thrombosis, decrease cardiovascular and resist cancer [13]. According to the location of the last double carbon bond relative to the terminal methyl end of the molecule, PUFA can be divided into omega-3 (n-3) series and omega-6 (n-6) series. ALA, an important omega-3 PUFA, is the necessary fatty acid which cannot be synthesized in the human body and must be obtained from dietary intake and is the precursor of docosahexaenoic acid (DHA) and eicosapentaenoic acid (EPA), which are essential to the growth and development of the brain and retina [14]. As well, LA, an important omega-6 PUFA, is the necessary fatty acid which cannot be produced by human body and is the precursor to the omega-6 series of fatty acids such as γ-Linolenic acid (C18:3Δ6,9,12, GLA) and Arachidonic acid (C18:4Δ5,8,11,14, AA) [15]. Nevertheless, the high intake of omega-6 PUFA with low intake of omega-3 PUFA might promote the pathogenesis of many diseases, including cardiovascular disease, cancer, and inflammatory and autoimmune diseases [16]. Therefore, the intake of the ratio of omega-6/omega-3 was lower than 5:1 suggested by Food and Agriculture Organization (FAO) of the Untied Nations and the World Health Organzation (WHO) [17].

The previous researches showed that the ratio of omega-6/omega-3 is about 1.5:1 in herbaceous peony ‘Hangshao’ mature seed oil. However, the ratio of omega-6/omega-3 is increased during the seed development because that the relative content of ALA was decrease gradually and the relative content of LA is increased slightly [10, 11]. Nevertheless, few research of the fatty acid accumulation at a molecular level for herbaceous peony seed oil had been reported. Fatty acid biosynthesis in plant is a very complicated biological pathway regulated by various enzymes. RNA sequencing (RNA-Seq) has provided convenience that the key genes of fatty acid biosynthesis and unsaturated fatty acid biosynthesis will be rapidly found at the molecular level. De novo transcriptome sequencing had been applied in identifying and profiling the different expression gene of some oil plants such as soybean [18], peanut [19], rape [20], sunflower [21], olive [22], perilla [23], oil palm [24], tree peony [17], tea-oil camellia [25], Eucommia ulmoides [26], sea buckthorn [27], Plukenetia volubilis [28] and so on. However, herbaceous peony research was mainly concentrated on the thermo-tolerant related differently expressed genes [29], codon usage patterns [30] and anthocyanin biosynthetic genes [31] based on transcriptom sequencing. The fatty acid metabolism of herbaceous peony remained unknown. Meanwhile, the genetic control of the unsaturated fatty acid accumulation of herbaceous peony was unexplored. In this study, the absolute content of five main fatty acids for herbaceous peony ‘Hangshao’ developing seed was measured by gas chromatography-mass spectrometry (GC-MS), and the first herbaceous peony seed transcriptome is generated to uncover genes related to fatty acid biosynthesis using high-throughput Illumina sequencing technology. In total, 150156 unigenes were obtained from nine samples of herbaceous peony ‘Hangshao’ seeds transcriptome within three developing stages (30, 60 and 90 days after flower, DAF) and most of genes in fatty acid biosynthesis were identified. Nine differentially expressed genes were selected for quantitative real-time PCR (qRT-PCR) analyses at five developmental stages (30, 45, 60, 75 and 90 DAF) which provided a large-scale transcriptome annotation of herbaceous peony ‘Hangshao’ and valuable gene resources related to fatty acid biosynthesis in seeds of herbaceous peony ‘Hangshao’.

Results

Fatty acid composition and their absolute content change

Five main fatty acids, PA, SA, OA, LA and ALA, were detected through GC-MS analysis of the fatty acid components of the seeds from herbaceous peony 'Hangshao' (Fig. 1A). Through quantitative analysis of the contents of five fatty acids (Fig. 1B), it was found that the content of ALA, LA and PA showed a trend of first rising and then falling during seed development, and the maximum value appeared in the S4 period. While the content of OA and SA increased all the time during seed development, and the maximum value appeared in S5 period. However, the absolute contents of the five fatty acids in the five periods showed that ALA content was the highest and SA content was the lowest, and the contents were successively ALA > LA > OA > PA > SA. The total contents of five fatty acids from S1 to S5 seeds were 44.9 mg·g-1, 121.70 mg·g-1, 211.9 mg·g-1, 261.1 mg·g-1 and 247.5 mg·g-1, respectively. Meanwhile, the content of ALA was exceeded 100 mg·g-1 in both S4 and S5 stages.

Sequencing reads filtering and de novo assembly

Nine cDNA libraries containing three biological repeats were created from three stages of developing herbaceous peony ‘Hangshao’ seeds (the early, middle and late stages of seed development, that is 30, 60, and 90 days after flower respectively) collected from the same strain and sequenced using Illumina high-throughput sequencing platform in order to explore the molecular mechanism of fatty acid biosynthesis and accumulation for herbaceous peony ‘Hangshao’. The sequencing reads which containing low-quality, adaptor-polluted and high content of unknown base (N) reads, should be processed to be removed before downstream analyses. The ratio of the amount of clean reads for nine samples was more than 82% (Additional file 1). Due to without reference genome of herbaceous peony, clean reads had to be assembled after sequencing to get reference sequence for subsequent analysis. After reads filtering, we used Trinity to perform de novo assembly with clean reads [32]. Then, Tgicl were used twice on cluster transcripts to remove abundance and get the final unigenes list for downstream analyses [33]. After assembling all samples together and filtering the abundance, we got 150,156 unigenes, the total length, average length, N50, and GC content of unigenes were 154,674,160 bp, 1,030 bp, 1,825 bp and 39.97% respectively (Table 1). The size of unigenes’ sequences mainly ranged from 300 to 3000 nt, and unigene number gradually decreased without obvious separation as the size of sequences increased. These results mean that a best continuity and high quality of RNA sequencing was conducted (Additional file 2). Of the unigenes, 58,645 (39.1%) were short sequences between 300 nt to 400 nt, whereas only 8927 (6.9%) were longer than 3,000 nt. These results showed that the short sequences advantaged over other sequences of the unigenes of herbaceous peony. All clean reads were uploaded to National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database under accession number SRP148668.

Functional annotation of non-redundant unigenes

After assembly, unigenes functional annotation were performed with seven functional database (NR, NT, GO, KOG, KEGG, SwissProt and InterPro). A total of 67,096 (44.68%), 45,488 (30.29%), 42,533 (28.33%), 49,630 (33.05%), 51,612 (34.37%), 48,907 (32.57%) and 29,413 (19.59%) unigenes had the most significant BLAST matches with known proteins in the NR, NT, GO, KOG, KEGG, SwissProt and InterPro database, respectively. Of the unigenes, 73,399 (48.88%) unigenes were annotated by any of the seven functional databases, whereas 15,005 (9.99%) unigenes were annotated by all the seven functional databases (Table 2).

The similarity analysis between herbaceous peony ‘Hangshao’ seed unigenes and NT Nuleic acid database and NR protein database of NCBI were conducted using Blastn [34]. Based on the NR database function annotation result, the ratio of different species on unigene annotataion was calculated and the distributions map was drew (Fig 2A). The unigenes of herbaceous peony ‘Hangshao’ seed were matched to those of four species such as Vitis vinifera (32.79%), Nelumbo nucifera (5.84%), Theobroma cacao (5.45%) and Jatropha curcas (3.66%).

Based on the 25 functional groups in KOG database, the unigenes were annotated and the unigene distribution was calculated (Fig 2B). Among them, unigenes mapped to ‘General function prediction only’ category was the most with the number of 17,454 (25.89%), followed by ‘Signal transduction mechanisms’ category with the number of 6,039 (8.96%), and ‘Posttranslational modification, protein turnover, chaperones’ category with the number of 5,185 (7.69%), whereas the unigenes mapped to ‘Cell motility’ category (88, 0.13%) were the least. The number of unigenes mapped to ‘Lipid transport and metabolism’ category was 1,868 (2.77%).

The unigenes which aligned to NR database had been annotated to GO databases with Blast2GO [35]. According to GO functional annotation, 29,413 unigenes were obtained and categorized into 55 functional groups in three domains including biological process, cellular component and molecular functions (Fig 2C). The biological process category was further classified into 24 groups, metabolic process (15,254 ) and cellular process (14,489) were the predominant groups, followed by single-organism process (8,816), biological regulation (3,594), localization (3,412) and regulation of biological process (3,063). The cellular component category was further classified into 17 groups, cell (11,981) and cell part (11,901) were the largest two groups, followed by membrane (9,845), organelle (8,516) and membrane part (7,363). The molecular function category was classified into 14 groups, the predominant groups were catalytic activity (14,405) and binding (13,721), followed by transporter activity (1,877) and structural molecule activity (872). However, there were less than 10 unigenes in the GO terms of ‘cell killing’ (9), ‘locomotion’ (7), ‘biological adhension’ (5), ‘metallochaperone activity’ (5), ‘protein tag’ (5) and ‘biological phase’ (2).

Pathway-based analysis can provide convenience for further understanding the biological functions and interactions of genes. The 49,630 unigenes were annotated into six KEGG categories (including cellular processes, environmental information processing, genetic information processing, human diseases, metabolism and organismal systems) and 21 sub-categories (Fig 2D). Among the six categories, ‘metabolism’ had the largest number of unigenes (30,888), followed by ‘genetic information processing’ (12,739 unigenes), ‘environmental information processing’ (3,668 unigenes), ‘cellular processes’ (2,403 unigenes), ‘organismal systems’ (1,831 unigenes) and ‘human diseases’ (139 unigenes). Among the 21 sub-categories, ‘global and overview maps’ was the maximum group with 12,390 unigenes, followed by ‘carbohydrate metabolism’ (4,905 unigenes) and ‘translation’ (4,321 unigenes), and the smallest group containing only 17 unigenes was ‘drug resisistance: antimicrobial’. Whereas, the number of the unigenes for ‘lipid metabolism’ was 1,766.

Unigenes related to fatty acid biosynthesis

Based on the analysis of the unigenes in the KEGG pathway and functional annotation, the critical enzymes related to lipid metabolism pathways had been classified and listed in Table 3. According to the previous studies and references, the map of fatty acid biosynthesis pathway was constructed on the basis of these identified enzymes, including de novo synthesis of fatty acid, fatty acid elongation and biosynthesis of unsaturated fatty acid pathway (Fig 3). A total of 247 unigenes related to fatty acid biosynthesis were obtained, including 103 unigenes related to de novo synthesis of fatty acid, 74 unigenes related to fatty acids elongation and 70 unigenes related to biosynthesis of undesaturated fatty acid.

Some critical enzymes in the fatty acid biosynthesis of herbaceous peony ‘Hangshao’ seed were indentified. Firstly, acetyl-CoA was catalyzed to malonyl-CoA by the enzyme acetyl-CoA carboxylase (ACCase, EC: 6.4.1.2) which consisted of four subunits: biotin carboxylase (BC), carboxyl transferase subunit alpha (α-CT), carboxyl transferase subunit beta (β-CT) and biotin carboxyl carrier protein (BCCP). We identified 11, 8, 6 and 17 unigenes for BC, α-CT, β-CT and BCCP, respectively. Then, malonyl-CoA was catalyzed to malonyl-ACP which was the primary substrate for a subsequent cycle of condensation reactions by malonyl-CoA-acly carrier protein transaclase (FabD, EC: 2.3.1.39). Only 1 unigene was identified to FabD. Subsequently, malonyl-ACP combined with acetyl-CoA were condensed to β-ketobutyryl-ACP by 3-oxoacyl-ACP synthase III (FabH, EC: 2.3.1.180). Only 1 unigene for FabH was identified. Then, β-ketobutyryl-ACP was generated to β-hedroxybutyryl-ACP on reduction by 3- oxoacyl-ACP reductase (FabG, EC: 1.1.1.100), and 26 unigenes were indentified to FabG. Next, the β-hedroxybutyryl-ACP was dehydrated by 3-hydroxyacyl-ACP dehydratase (FabZ, EC: 4.2.1.59) and was generated to trans2-butenoyl-ACP. 5 unigenes were identified to FabZ. Then, trans2-butenoyl-ACP was reduced to butyryl-ACP (C4:0-ACP) by enoyl-ACP reductase I (FabI, EC:1.3.1.9), and 5 unigenes were indentified to FabI. Subsequently, after 6 cycles for a series of condensation, reduction, dehydration and reduction reaction by 3-oxoacyl-ACP synthase II (FabF, EC:2.3.1.179), FabG, FabZ and FabI, C4:0-ACP was generated to palmitoyl-ACP (C16:0-ACP). Meanwhile, after 7 cycles for those reactions, C4:0-ACP was generated to stearoyl-ACP (C18:0-ACP). 18 unigenes were indentified to FabF. Under the role of fatty acyl-ACP thioesterase B (FATB, EC:3.1.2.14 3.1.2.21) and fatty acyl-ACP thioesterase A (FATA, EC:3.1.2.14), C16:0-ACP and C18:0-ACP were catalyzed to palmitic acid (C16:0) and stearic acid (C18:0) respectively. 4 unigenes were indentified to FATB, however, only 1 unigene was indentified to FATA.

Fatty acid elongation were conducted from long-chain acyl-CoA (C≥16). Fatty acid (C≥16) which catalyzed by thiolase from fatty acyl-ACP, was catalyzed to long-chain acyl-CoA (C≥16) by long-chain acyl-CoA syntetase (fadD, EC:6.2.1.3) and 30 unigenes were identified to fadD. After a series of complicated biochemical reactions, long-chain acyl-CoA (C≥16) was catalyzed to long-chain acyl-CoA (C+2) by 3-ketoacyl-CoA synthase (KCS, EC:2.3.1.199), very-long-chain 3-oxoacyl-CoA reductase (KCR, EC:1.1.1.330), very-long-chain (3R)-3-hydroxyacyl-CoA dehydratase (HCD, EC:4.2.1.134) and very-long-chain enoyl-CoA reductase (ECR, EC:1.3.1.93). 23, 7, 8 and 6 unigenes were indentified to KCS, KCR, HCD and ECR respectively.

Biosynthesis of unsaturated fatty acid is conducted at two pathway including plastid pathway and endoplasmic reticulum pathway, and was started from C18:0-ACP. Firstly, with the dehydrogenation by stearoyl-ACP desaturase (SAD, EC:1.14.19.1 ), C18:0-ACP was synthesized to oleoyl-ACP (C18:1-ACP) in plastid. And 37 unigenes were identified to SAD. Then, some C18:1-ACP was dehydrogenized to linoleoyl-ACP (C18:2-ACP) by FAD6 in plastid. 2 unigenes were indentified to FAD6. Next, C18:2-ACP was dehydrogenized by FAD7 and FAD8 (EC:1.14.19.35) to form alpha-linolenony-ACP (C18:3-ACP) in plastid. 8 and 2 unigenes were encoded to FAD7 and FAD8. However, some C18:1-ACP was catalyzed to C18:1 by FATA, which was catalyzed to C18:1-CoA by fadD, and C18:1-CoA was transferred into the endoplasmic reticulum. Then, C18:1-CoA and lysophosphatidylcholine (LPC) were catalyzed to C18:1-PC by lysophosphatidylcholine acyltransferase (LPCAT). Next, C18:1-PC was dehydrogenized to C18:2-PC by FAD2, and C18:2-PC was dehydrogenized to C18:3-PC by FAD3. 16 and 5 unigenes were indentified to FAD2 and FAD3.

Analysis of differentially expressed genes for developmental herbaceous peony ‘Hangshao’ seed

In order to deeply understand the patterns of the seed development and excavate the key genes of the fatty acid biosynthesis, differentially expressed genes (DEGs) by comparing FPKM value of unigenes between the different developmental seeds of herbaceous peony ‘Hangshao’ were analyzed. Based on the gene expression level, DEGs were identified between different groups (Fig 4A). Seeds at 30 DAF were served as the control, a total of 4,635 and 18,291 differentially expressed unigenes were indentified at 60 and 90 DAF, respectively. And seeds at 60 DAF were served as the control, a total of 12,304 differentially expressed unigenes were also identified at 90 DAF. There were 2,410, 7,385 and 4,210 up-regulated unigenes and 2,225, 10,906 and 8,094 down-regulated unigenes in 30 vs 60 DAF, 30 vs 90 DAF and 60 vs 90 DAF, respectively. Two hierarchical cluster analysis were conducted based on the RMKM value of these differentially expressed unigenes for three group of 30 vs 60 DAF, 60 vs 90 DAF and 30 vs 90 DAF (Fig 4B). The DEGs intersection of three difference group (Inter) and the DEGs union of three difference group were clustered into three groups. DEGs in the same cluster meant that they have identical or similar expression patterns during the seed developmental stage.

The DEGs at different developmental stages were integrated and mapped into the GO database (Additional file 3 and Additional file 4). Overall, the number of DEGs annotated to the same GO term between the three different group was increased with the developmental stages. The DEGs of the group for 30 vs 90 DAF were annotated to 23 biological processes, 17 cellular components and 14 molecular functions. In contrast, the DEGs of the group for 60 vs 90 DAF were annotate to 21 biological processes (except for ‘biological adhesion’ and ‘locomotion’), 17 cellular components and 13 molecular functions (except for ‘metallochaperone activity’). And the DEGs of the group for 30 vs 60 DAF were annotate to 21 biological processes (except for ‘biological adhesion’ and ‘rhythmic process’), 14 cellular components (except for ‘nucleoid’, ‘virion’ and ‘virion part’) and 12 molecular functions (except for ‘metallochaperone activity’ and ‘protein tag’).

DEGs at the three different group of 30 vs 60 DAF, 60 vs 90 DAF and 30 vs 90 DAF were annotated into five KEGG categories (including cellular processes, environmental information processing, genetic information processing, metabolism and organismal systems) and 19 KEGG sub-categories (Additional file 5 and Additional file 6). The number of DEGs annotated to KEGG pathway was 1,698, 4,321 and 6,035 for the three group of 30 vs 60 DAF, 60 vs 90 DAF and 30 vs 90 DAF, respectively. Similarly, with the maturity of seeds, the number of DEGs annotated to the 19 sub-categories was increased. For example, the number of DEGs annotated to ‘Lipid metabolism’ was 191, 424 and 544 for the group of 30 vs 60 DAF, 60 vs 90 DAF and 30 vs 90 DAF, respectively.

In addition, a total of 1,480 DEGs were obtained from the intersection analysis of the three groups of 30 vs 60 DAF, 60 vs 90 DAF and 30 vs 90 DAF (Fig 5A). And the total of 1480 DEGs was annotated to 17 biological processes, 14 cellular components and 12 molecular functions in the GO database (Fig 5B). In biological process, the largest group was ‘metabolic process’ which was annotated by DEGs, followed by ‘cellular process’ and ‘single-organism process’. And in cellular component, the largest group was ‘membrane’ which annotated by DEGs, followed by ‘cell’ and ‘cell part’. In addition, in molecular function, the two largest group were ‘catalytic activity’ and ‘binding’. Moreover, the total of 1480 DEGs was annotated to 19 KEGG pathways (Fig 5C). The largest group was ‘global and overview maps’ which was annotated by DEGS, followed by ‘carboydrate metabolism’ and ‘lipid metabolism’.

Analysis of differentially expressed genes for Fatty acid biosynthesis

In order to search the key genes in the fatty acid synthesis of developmental herbaceous peony ‘Hangshao’ seed, DEGs in three KEGG pathway of fatty acid synthesis including de novo synthesis of fatty acid, fatty acid elongation and biosynthesis of unsaturated fatty acid were analyzed (Additional file 7). In the group of 30 vs 60 DAF, 19 DEGs were indentified in the KEGG pathway of de novo synthesis of fatty acid which included 2 up-regulated DEGs and 17 down-regulated DEGs. The 2 up-regulated DEGs were indentified to FabG, and the 17 down-regulated DEGs were indentified to BC(1), α-CT(2), BCCP(4), FabD(1), FabZ(4), FabI(3) and FabF(2), respectively. In addition, 6 DEGs were indentified in the KEGG pathway of fatty acid elongation, in which 1 up-regulated DEG was indentified to KCR and 5 down-regulated DEGs were indentified to KCS. However, in the KEGG pathway of biosynthesis of unsaturated fatty acid, no DEGs were indentified. In the group of 60 vs 90 DAF, 29 DEGs were indentified in the KEGG pathway of de novo synthesis of fatty acid which included 4 up-regulated DEGs and 25 down-regulated DEGs. The all 4 up-regulated DEGs were indentified to FabG, and the 25 down-regulated DEGs were identified to BC(3), α-CT(1), BCCP(3), FabD(1), FabG(7), FabZ(3), FabI(2), FabF(2), FATB(2) and FATA(1), respectively. In addition, 9 DEGs were indentified in the KEGG pathway of fatty acid elongation, in which 1 up-regulated DEGs was indentified to fadD(1) and 8 down-regulated DEGs were indentified to KCS(2), KCR(3), HCD(2) and ECR(1), respectively. Moreover, 22 DEGs were indentified in the KEGG pathway of biosynthesis of unsaturated fatty acid, in which 2 up-regulated DEG was indentified to SAD and 20 down-regulated DEGs were indentified to SAD(5), FAD2(8), FAD3(1), FAD7(4) and FAD8(2). In the group of 30 vs 90 DAF, 39 DEGs were indentified in the KEGG pathway of de novo synthesis of fatty acid which included 6 up-regulated DEGs and 33 down-regulated DEGs. The 6 up-regulated DEGs were indentified to BCCP(1) and FabG(5), and the 33 down-regulated DEGs were identified to BC(3), α-CT(3), BCCP(5), FabD(1), FabH(1), FabG(7), FabZ(4), FabI(3), FabF(4), FATB(1) and FATA(1), respectively. In addition, 24 DEGs were indentified in the KEGG pathway of fatty acid elongation, in which 4 up-regulated DEGs were indentified to fadD(2) and TER(2), and 20 down-regulated DEGs were indentified to fadD(1), KCS(12), KCR(4), HCD(2) and ECR(1), respectively. Moreover, 23 DEGs were indentified in the KEGG pathway of biosynthesis of unsaturated fatty acid, which were all down-regulated DEGs, and they were indentified to SAD(5), FAD2(9), FAD3(2), FAD7(5) and FAD8(2). Furthermore, 13 DEGs were identified in the all three groups of 30 vs 60 DAF, 60 vs 90 DAF and 30 vs 90 DAF which were belonged to BC(1), α-CT(1), BCCP(2), FabD(1), FabG(1), FabZ(3), FabI(2), FabF(1) and KCR(1).

Experimental validation and analysis of nine DEGs involved in fatty acid biosynthesis

To validate these DEGs, the relative expression levels of 9 randomly selected DEGs were analyzed by qRT-PCR. Among them, there are 5 DEGs related to de novo synthsis of fatty acid, namely BCCP, BC, FabD, FabF and FATB, and 1DEGs was involved in fatty acid elongation (KCR).There were 3 DEGs involved in the biosynthesis of unsaturated fatty acids, including FAD2, FAD3 and FAD7 (Fig 6). Overall, the expression patterns of 9 DEGs verified by qrt-pcr were highly consistent with the results of transcriptome sequencing, indicating that the results of DEGs analysis were reliable.

Discussion

Lipid, which was produced by fatty acid and glycerinum as substrate, was stored in the seed in the form of triacylglycerol (TAG). The accumulation of fatty acids is closely related to the accumulation of lipids, that is, the increase of fatty acid content will accelerate the accumulation of lipids [36]. The proportion of fatty acids, especially the proportion of unsaturated fatty acids, is one of the evaluation criteria for oil quality [37]. However, the content of LA and ALA constituted different proportions of essential fatty acids have different effects on human physiological functions [38, 39]. At present, many studies have shown that the intake ratio of LA/ALA has the most significant physiological function at 1:1 [40]. In this study, five main fatty acids were detected by GC-MS in the seed of herbaceous peony 'Hangshao', that is PA, SA, OA, LA and ALA, and which was consistent with the study on the component of other peony [11, 17]. The absolute contents of the five fatty acid components during the seed development were all the highest for ALA and the lowest for SA. The specific manifestations were ALA>LA>OA>PA>SA. The ratio of LA/ALA was lower during the developmental seeds of herbaceous peony 'Hangshao', maintained between 0.55 and 0.71, close to the optimal ratio of 1:1.

At present, the studies on fatty acid of herbaceous peony are limited to the extraction of seed oil and its component, the related gene cloning and expression of ALA synthesis [11], and few research of the fatty acid biosynthesis for herbaceous peony seed had been reported. With the development of sequencing technology, species without reference genome can use RNA-seq to directly sequence the total mRNA of samples. Transcriptome sequencing is also sensitive to changes of gene expression, which has important reference significance for species with complex genomic backgrounds such as herbaceous peony [17]. In this study, Illumina HiSeqTM2000 sequencing platform was used to carry out transcriptional sequencing on the seeds of herbaceous peony 'Hanghao' at different development stages, and 100.07 Gb data were obtained. After assembling and removing the redundancy, 150,156 unigenes were obtained, and 67,096 unigenes were annotated in the NR database. In order to obtain comprehensive transcriptome information of herbaceous peony 'Hangshao' seeds, the following measures were adopted in this study. Firstly, the seeds were collected from the same plants at early (30 DAF), middle (60 DAF) and mature (90 DAF) stages, and cDNA libraries were constructed from the total RNA of the seeds in 3 follicles at each stages, and transcriptome was sequenced in parallel. Secondly, for unigenes annotation, seven databases (NR, NT, SwissProt, KOG, KEGG, GO and InterPro) were used to get the most comprehensive annotation information possible.

According to the annotation information of unigene from transcriptome database, two important pathway related to fatty acid biosynthesis were found in the developmental seed of herbaceous peony 'Hangshao'. One pathway of fatty acid biosynthesis is in the plasmid, and the other is in the endoplasmic reticulum. In terms of fatty acid biosynthesis in the plasmids, the genes catalyzing the generation of malonyl-CoA including 11 BC, 8 α-CT, 6 β-CT and 17 BCCP were detected in this study. FabD, which catalyzes the formation of malonyl-ACP, has only one unigene. And the genes catalyzing the generation of butyryl-ACP (C4:0-ACP) including 1 FabH, 26 FabG, 5 FabZ and 5 FabI were detected. In addition of the above FabG, FabZ and FabI, the genes catalyzing the generation of palmityl-ACP (C16:0-ACP) and stearoyl-ACP (C18:0-ACP) also include 18 FabF. Furthermore, the fatty acid dehydrogenase gene including 37 SAD, 2 FAD6, 8 FAD7 and 2 FAD8. In terms of fatty acid biosynthesis in the endoplasmic reticulum, the genes catalyzing the generation of long-chain fatty acid including 3 KCS, 7 KCR, 8 HCD and 6 ECR were detected. And the fatty acid dehydrogenase gene including 16 FAD2 and 5 FAD3 were also detected in the endoplasmic reticulum. Fatty acid dehydrogenase FAD2/FAD6 catalyzes the dehydrogenation of OA at Δ12 position to form LA [41], and catalyze LA to form ALA by Δ15 fatty acid dehydrogenase including FAD3, FAD7 and FAD8 [42]. In this process, FAD genes play a key role in the formation of ALA [43]. After the synthesis of fatty acids, the TAG was assembled in the endoplasmic reticulum through Kennedy pathway, and finally formed into oil body by budding for storage [44].

At present, 17 lipid metabolism pathway and a series of related protein have been annotated in the KEGG database, including 33, 23 and 29 related genes in the pathway for de novo synthesis of fatty acid, fatty acid elongation and biosynthesis of unsaturated fatty acid, respectively [45]. In the de novo synthesis of fatty acids, acetyl CoA is condensed with CO2 to form malonyl CoA, and then malonyl CoA is catalyzed by FabD to form malonyl-ACP. Then, malonyl-ACP is reacted with acetyl CoA to form butyryl-ACP (C4:0-ACP) through condensation, reduction, dehydration, reduction by a series of enzymes. In this study, the unigenes of 1 BC (down), 1 α-CT (down), 2 BCCP (down), 1 FabD (down), 1 FabG (up), 3 FabZ (down) and 1 FabI (down) were also differentially expressed in the two stages of 30 d vs 60 d and 60 d vs 90 d. Butyryl-ACP and malonyl-ACP undergo 6 or 7 cycles of condensation, reduction, dehydration and reduction by a series of enzymes to generate palmityl-ACP (C16:0-ACP) and stearoyl-ACP (C18:0-ACP), respectively. 1 FabF (down) was also differetially expressed in the two stages of 30 d vs 60 d and 60 d vs 90 d except for FabG, FabZ and FabI mentioned above. Fatty acid elongation was carried out in the endoplasmic reticulum, and the catalytic enzymes were KCS, KCR, HCD and ECR. In the two stage of 30 d vs 60 d and 60 d vs 90 d, no differentially expressed genes were also found. However, differentially expressed genes including 4 KCS (down) and 1 KCR (up) were found in the stage of 30 d vs 60 d, and 1 KCS (down), 3 KCR (down), 2 HCD (down) and 2 ECR (down) were found in the stage of 60 d vs 90 d. The synthesis pathway of unsaturated fatty acid was partially carried out in plastids, C18:0-ACP was desaturated to generate C18:1-ACP by SAD, and then was desaturated to generate C18:2-ACP and C18:3-ACP by FAD6 and FAD7/FAD8, respectively. In addition, the synthesis pathway of unsaturated fatty acid was partially in endoplasmic reticulum. C18:1-CoA was catalyzed by LPCAT to generate C18:1-PC, then was desaturated by FAD2 and FAD3 to generate C18:2-PC and C18:3-PC. No differentially expressed genes were also found in the two stages of 30 d vs 60 d and 60 d vs 90 d. Furthermore, no differential expression of SAD, FAD2, FAD3, FAD6, FAD7 and FAD8 were found in the stage of 30 d vs 60 d. However, 22 DEGs including 7 SAD (2 up and 5 down), 8 FAD2 (down), 1 FAD3(down), 4 FAD7 (down) and 2 FAD8 (down) were found in the stage of 60 d vs 90 d. The above DEGs expression showed a general downward trend, which may be due to the fact that the peak of lipid accumulation was 75DAF, while the maximum rate of unsaturated fatty acid accumulation was 45 DAF (Fig 1). Although the accumulation amount of fatty acid increased slowly, the accumulation rate of fatty acid decreased gradually. According to transcriptome, the expression of desaturase (SAD, FAD2, FAD3, FAD6, FAD7 and FAD8) showed no significant change at the stage of 30 d vs 60 d, possibly due to the maximum accumulation was happened at 45 DAF (Fig 6). However, with the decrease of fatty acid accumulation rate, the expression level of FAD gene family was gradually reduced from 60 d to 90 d. In the genes of BC, α-CT, BCCP, FabD, FabH and FabG that catalyzed the synthesis of malonyl-CoA in the plastid, the expression of FabG (CL2376.Contig2_All) showed a gradually increasing trend during seed development, however, the expression of other genes showed a decreasing trend, which was consistent with the gradually decreasing rate of fatty acid synthesis.

In the transcriptome of herbaceous peony 'Hangshao' seed, the unigene annotated to FabD, FabH and FATA was only 1, respectively. FabD is malonyl CoA-acyl carrier protein transacylase (MACT), which gene encodes malonyl CoA-ACP transallase that binds to malonyl-CoA to form a stable malonyl-MACT intermediate. Then, the acyl protein binds to the malonyl-MACT intermediate, and the malonyl group is transferred to the mercaptan at the end of the acyl carrier protein by hydrophobic interaction, so as to synthesize malonyl-ACP [46]. Malonyl ACP is a carbon chain donor that extends the carbon chain in the de novo synthesis of fatty acids. Under the catalysis of the complex enzyme system of FAS, it eventually forms C16:0-ACP and C18:0-ACP through multiple cycles. So, FabD is a key enzyme in the fatty acid synthesis pathway, and its expression level is closely related to the lipid accumulation. At present, the overexpression of FadD gene in rape and other plants has a good effect on promoting the accumulation of C16 and C18 fatty acid [47]. FabD belongs to DEGs in the transcriptome of herbaceous peony 'Hangshao' seed, and shows a downward trend in the seed development, which is consistent with the change trend of fatty acid accumulation rate. FabH is ketoacyl-ACP synthase III (KASIII), a key enzyme that initiates de novo synthesis of fatty acid. It acts on the condensation in the first cycle of fatty acid synthesis, which catalyzes condensation of malonyl ACP with acetyl CoA to produce β-ketobutyryl-ACP. Then after reduction, dehydration and reduction, β-ketobutyryl-ACP was catalyzed to butyryl-ACP (C4:0-ACP) [48]. At present, FabH gene has been cloned from plants such as tung tree [49], sunflower [50] and jatropha [51]. Condensation reaction by FabH was conducted that acetyl CoA receptor substrate release CoA and bind C2 unit offered by malonyl-ACP after decarboxylation to form β-ketobutyryl-ACP. In addition, some studies showed that FabH enzymes can use other receptor substrate to generate different final products of fatty acid synthesis, which indicate that FabH enzyme is one of the main factors determining the carbon chain structure of fatty acid [52, 53]. FabH, as the initator of de novo synthesis of fatty acid, plays a second role in detemining the length of the carbon chain of fatty acid [54]. Moreover, some studies believe that the synthesis rate of fatty acid is affected by FabH enzymatic activity, indicating that it may be one of the rate-limiting enzymes in the fatty acid synthesis pathway [55]. The FPKM value of FabH gene in the transcriptome of herbaceous peony 'Hangshao' seeds decreased with the development of seeds, but had no significant difference, which was consistent with the accumulation rate of fatty acid in herbaceous peony 'Hangshao' seeds. FATA is acyl-ACP thioesterase A, a type of acyl-ACP thioesterase. De novo synthesis of fatty acid was terminated by acyl-ACP thioesterase (FAT), which hydrolyzes the thiol-lipid bond located between acyl and ACP to form fatty acid. Then, fatty acid are transported from the plastids and then are esterified to form acly-CoA in the plasmid membrane [56]. The specificity of FAT largely determines the chain length and unsaturated degree of vegetable fatty acid in glycerolipid and TAG [57]. FAT is a soluble enzyme encoded by nuclear genes, which can be classified into FATA and FATB according to sequence alignment. FATA and FATB have different substrate specificity. According to all the tested FATA homologous body, they have similar substrate specificity among the different species. And FATA has the highest catalytic activity of C18:1Δ9-ACP. The FPKM value of FATA gene in the transcriptome of herbaceous peony 'Hangshao' seeds also decreased with the development of seeds, but has no obvious different in the stage of 30 d vs 60 d. However, FATA was belonged to DEGs at the stage of 30 d vs 60 d, which FPKM value has obvious different. The tendency of FATA expression was consistent with the total accumulation rate of C18:1Δ9, and its later product C18:2Δ9,12 and C18:3Δ9,12,15 in the seed of herbaceous peony 'Hangshao'.

Conclusions

GC-MS analysis revealed that ‘Hangshao’ seeds have five main fatty acids, including palmitic acid, stearic acid, oleic acid, linoleic acid and α-linolenic acid, respectively. The order of α-linolenic acid, linoleic acid, oleic acid, palmitic acid and stearic acid was arranged in the order of highest to lowest content. The total content of the five fatty acids increased first and then decreased, and reached the maximum at 75 DAF. A total of 150,156 Unigenes were obtained by transcriptome sequencing, with an average length of 1,030 bp. There were 1,550 Unigenes annotated in the seven functional databases including NR, NT, GO, KOG, KEGG, SwissProt and InterPro. Based on KEGG database, 1,766 Unigenes were annotated in the lipid metabolic pathways, among which 103, 74 and 70 Unigenes are annotated into fatty acid biosynthesis pathway, fatty acid elongation pathway and unsaturated fatty acid synthesis pathway; respectively. By pairwise compared the transcription of ‘Hangshao’ seeds of three development stages, a total of 2,410 (30 DAF vs 60 DAF), 4,210 (60 DAF vs 90 DAF) and 7,385 (30 DAF vs 90 DAF) up-regulated DEGs were obtained, respectively, while a total of 2,225 (30 DAF vs 60 DAF), 8,094 (60 DAF vs 90 DAF) and 10,906 (30 DAF vs 90 DAF) down-regulated DEGs were obtained, respectively. Combined with the comparison, a total of 1480 DEGs were detected. Among them, 83 DEGs were enriched in the fatty acid metabolism pathway, including 12 DEGs involved in the fatty acid biosynthesis and 1 DEG involved in fatty acid elongation. Furthermore, qRT-PCR was used to analyze the expression patterns of nine fatty acid biosynthetic related genes including FBCP, BC, fabD, fabF, FATB, KCR, FAD2, FAD3 and FAD7, and it showed that they all highest expressed at 45 DAF.

Methods

Plant materials

Seeds of herbaceous peony 'Hangshao' were collected in 2017 at peony germplasm resource garden of Yangzhou University, China (32°39′N, 119°42′E). It had been introduced to the garden and grown in the same eco-environmental and cultivation conditions for eight years. Due to the limited fruit and low oil extraction rate at early stage, the follicle were hand-collected in the 30 days after flowering (DAF) for the first time and at intervals of 15 days until full maturity (including 30 DAF, 45 DAF, 60 DAF, 75 DAF and 90 DAF) for fatty acid absolute content analysis by GC-MS (Fig 7). Seeds used for tanscriptome sequencing were collected from the same strain at developing stages for 30, 60 and 90 DAF. The collected samples were immediately frozen in liquid nitrogen and stored at -80℃ cryogenic refrigerator until further use.

Quantitative analysis of main fatty acid composition

Sample hydrolysis: 0.5g seeds were ground into granules and put into a flask. Then, 100 μL methyl nonadecanoate, 100 mg pyrogallic acid, several zeolites and 2 mL 95% ethanol were added and mixed. Subsequently, 10 mL hydrochloric acid solution were added and mixed well. Next, the flask was put into 70℃ water bath to hydrolyze for 40 min. After hydrolysis, the flask was cooled to room temperature.

Fat extraction: the hydrolyzed sample was added with 10 mL 95% ethanol and mixed. And the hydrolysate in the flask was transferred to the separation funnel. Then, the flask and the stopper were washed with 50 mL ethyl ether petroleum ether mixture. Next, the rinse solution was incorporated into the separation funnel, capped, shaken for 5 min, and set aside for 10 min. Subsequently, the extracting solution was collected into 250 mL flask. The hydrolysate was extracted 3 times as described in the previous steps. Lastly, separatory funnel was flushed by petroleum ether with ether mixture (volume ratio 1:1) and the flushing fluid were collected into a constant weight flash which would be boiled in water bath and dried in 100℃ for 2 h.

Fatty acid esterification: 2 mL of 2% sodium hydroxide methanol solution was added into the fat extract, and heated in 80℃ for 30 min. Then, 3 mL 14% boron trifluoride methanol solution was added and heated in 80℃ for 30 min. Next, the mixture was cooled to room temperature and put into a centrifugal tube and 1 mL n-hexane was added. And it was extracted by shock for 2 min, then left for 1 h, and then stratified. 100 μL supernatant was volumed to 1 mL with n-hexane and filtered by 0.45 μm organic phase filter to 2 mL chromatographic sample bottle.

Chromatographic and mass spectrometric conditions: the conditions was referred to the reference [10].

Calculation formula of fatty acid absolute content: W (mg•g-1) =C×V×N×K/ (1000×M), where W is the absolute content of fatty acid, C is the concentration of fatty acid methyl ester in the sample determination solution (mg• L-1), V is the constant volume (mL), N is the dilution factor, K is the conversion coefficient of fatty acid triglyceride into fatty acid, and M is the sample mass (g).

RNA isolation and cDNA library construction

Nine seed samples from three stage (30, 60 and 90 DAF, with three biological replicates) were snap-frozen in liquid nitrogen and ground to a fine powder. Total RNA samples were extracted through Total RNA Isolation System (Takara, Japan) according to the manufacturer’s instructions. The quality of the resulting RNA was verified by 2100 Bioanalyzer RNA Nano chip device (Agilent, Santa Clara, CA, USA). All RNA extractions were delivered an RNA integrity number(RIN) value of >8.0, and a 28S:18S ratio >1.5. After checking the A260/A280 nm ratios, and A260/A230 nm ratios by Nanodrop ND-430 1000 spectrophotometer, the nine total RNA samples were selected based on 28S/18S band intensity (1.5:1~2:1) and spectroscopic A260/A280 nm readings between 1.8~2.2, A260/A230 nm readings >2.0. Total RNA samples were treated with RNase-free DNaseⅠ(Takara, Japan) to degrade gDNA. After mixing with oligo (dT), the coated magnetic beads were used to concentrate the poly A mRNA.

The mRNA was degrade into ~200 nt fragments during incubation in a fragmentation buffer (The Beijing Genomics Institute). The first strand of cDNA was connected with random hexamer, the second strand was generated with DNA polymerase Ⅰ, buffer, dNTPs and RNase H. Then the double strand cDNA was purified with a QiaQuick PCR extraction kit. The end repair and addition of single nucleotide A were carried out under EB buffer. After sequencing adaptors were ligated into the fragments and agarose gel electrophore, suitable fragments were selected as templates for PCR.

RNA-Seq and identification of differentially expressed genes

The library was sequenced using an Illumina HiSeq located at the Beijing Genomics Institute (Shenzhen, China; http://www.genomics.cn/index). Datas were deposited in the US National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA, http://www.ncbi. nlm.nih.gov/Traces/sra; under accession number (SRP). Clean reads were obtained through data filtering after removing adaptor sequences, reads of unknown bases (N) and low quality reads. The clean reads were mapped onto the reference sequences using SOAP (2.21) software. A maximum of two mismatches was allowed in the alignment. The NCBI nonredundant protein (Nr) database (http://www.ncbi.nlm.nih.gov) and the Swiss-Prot protein database (http://www.expasy.ch/sprot) were used for blast search and annotation using an E-value cut off of 10-5. Functional annotation by gene ontology terms (GO, http://www.geneontology.org) was analyzed using the Blast2GO program [35]. The Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG, http://www.genome.jp/kegg), the major public pathway-related database [58] was also used to predict and classify possible functions. The unigene expression level was calculated using the fragments per kilobase of exon model per million mapped fragment (FPKM ) method [59] . The threshold for significantly differentially expressed genes (DEGs) was set at a fold change ≥ 2.0 and adjusted Pvalue ≤ 0.05, and DEseq2 method was used for DEGs detection [60]. DEG functions were explored through GO and KEGG pathway analysis and the terms which false discovery rate (FDR) no larger than 0.01 were defined as significant enrichment. This was performed to identify significantly enriched metabolic pathways.

qRT-PCR validation of differential transcription

Gene transcript levels were analyzed using qRT-PCR with a BIO-RAD CFX ConnectTM Optics Module (BioRad, USA). cDNA was synthesized from RNA using the PrimeScript® RT Reagent Kit With gDNA Eraser (TaKaRa, Japan). Actin (JN105299) was used as an internal control in P. lactiflora [61].All gene-specific primers in this study were shown in Additional file 8, and synthesized by Shanghai Sangon Biological Engineering Technology and Services Co., Ltd. (Shanghai, China). qRT-PCR was performed using SYBR® Premix Ex TaqTM (Perfect Real Time) (TaKaRa, Japan). The amplification was carried out under the following conditions: 55 ℃ for 2 min, followed by an initial denaturation step at 95℃ for 30 s, and 40 cycles at 95 ℃ for 5 s, 55 ℃ for 15 s, and 72℃ for 30 s. Relative expression levels of target genes were calculated by the 2-△△Ct comparative threshold cycle method [62]. The Ct values of the triplicate reactions were gathered using the Bio-Rad CFX Manager V1.6.541.1028 software (Bio-Rad, Hercules, CA, USA).

Statistical analysis

All experiments described here were repeated three times arranged in a completely randomized design. The results were analyzed for variance using the SAS/STAT statistical analysis package (version 6.12, SAS Institute, Cary, NC, USA). The data shown in the figures were means ± SDs and the different letters indicated significant differences (P < 0.05).

Abbreviations

DEG: Differentially expressed gene; UFA: Unsaturated fatty acid; OA: Oleic acid; LA: Linoleic acid; ALA: α-linolenic acid; PA: Palmitic acid; SA: Stearic acid; MUFA: Monounsaturated fatty acid; PUFA: Polyunsaturated fatty acid; DAF: Days after flower; qRT-PCR: Quantitative real-time PCR; GC-MS: Gas chromatography-mass spectrometry; GO: Gene ontology; KOG: Clusters of orthologous groups; KEGG: Kyoto encyclopedia of genes and genomes; BC: biotin carboxylase; BCCP: biotin carboxyl carrier protein; α-CT: carboxyl transferase subunit alpha; β-CT: carboxyl transferase subunit beta; FabD: malonyl-CoA-acly carrier protein transaclase; FabH: 3-oxoacyl-ACP synthase III; FabG : 3- oxoacyl-ACP reductase; FabZ : 3-hydroxyacyl-ACP dehydratase; FabI : enoyl-ACP reductase I; FabF : 3-oxoacyl-ACP synthase II; FATB : fatty acyl-ACP thioesterase B; FATA: fatty acyl-ACP thioesterase A; fadD: long-chain acyl-CoA syntetase; KCS: 3-ketoacyl-CoA synthase; KCR: very-long-chain 3-oxoacyl-CoA reductase; HCD: very-long-chain (3R)-3-hydroxyacyl-CoA dehydratase; ECR: very-long-chain enoyl-CoA reductase; SAD: stearoyl-ACP desaturase; FAD: fatty acid desaturase; LPC: lysophosphatidylcholine; LPCAT: lysophosphatidylcholine acyltransferase; FPKM: fragments per kilobase of exon model per million mapped fragment; TAG: triacylglycerol.

Declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

All data generated or analyzed during this study are included in this published article and its supplementary information files.

Competing interests

The authors declare that they have no competing interests.

Funding

Funding for this study came from the Three-Side Innovation Projects for Agricultural in Jiangsu Province (SXGC[2017]297).

Authors' contributions

JT designed and managed the experiments, and organized the manuscript. JSM carried out sequence data analysis and drafted the manuscript. JSM and YHT performed the experiments. YHT and JS contributed to manuscript revision. All authors have carefully read and approved the final manuscript.

Acknowledgements

Not Applicable

References

1. Qing KJ. Illustrtion of one hundred ornamental flowers bonsai-the herbaceous peony. Beijing: China Forestry Press;2004.(In Chinese)

2. Zhu MX, Li SN, You HD, Han B, Wang ZP, Hu YX, Liu YF. Simultaneous determination of paeoniflorin and albiflorin in Radix paeoniae Rubra by HPLC–DAD–ELSD. Acta Chromatogr. 2017;29(2):279-89.

3. Bi XX, Qu TG, Mu Y, Guan PP, Qu XD, Wang ZY, Huang XS. Anti-inflammatory effects, SAR, and action mechanism of monoterpenoids from Radix paeoniae Alba on LPS-stimulated RAW 264.7 cells. Molecules. 2017;22:715.

4. Xie PY, Cui LL, Shan Y, Kang WY. Antithrombotic effect and mechanism of radix paeoniae rubra. BioMed Res Internatinal, 2017;doi:10.1155/2017/9475074.

5. Lin MY, Chiang SY, Li YZ, Chen MF, Chen YS, Wu JY, Liu YW. Anti-tumor Effect of Radix paeoniae Rubra extract on mice bladder tumors using intravesical therapy. Oncol Lett. 2016;12(2):904-10.

6. Qiu ZK, He JL, Zeng J, Cheng JS, Nie H. Anti-PTSD-like effects of albiflorin extracted from Radix paeoniae Alba. J Ethnopharmacol. 2016;198:324-30.

7. Jiang DX, Chen YS, Hou XT, Xu JF, Mu X, Chen W. Influence of Paeonia lactiflora roots extract on cAMP-phosphodiesterase activity and related anti-inflammatory action. J Ethnopharmacol. 2011;137(1):914-20.

8. Zhang P, Zhang JJ, Su J, Qi XM, Wu YG, Shen JJ. Effect of total glucosides of paeony on the expression of nephrin in the Kidneys from diabetic rats. Am J Chinese Med. 2009;37(2):295-307.

9. Xu JZ, Sun YM, Yu XP, Wang Z. Study on method of Radix paeonia Alba producing and concocting integration processing. China J of Chinese Mater Med. 2014;39(13):2504-8.(In Chinese)

10. Meng JS, Jiang Y, Tao J. Fatty acid composition and PlFADs expression related to a-linolenic acid biosynthesis in herbaceous peony (Paeonia lactiflora Pall.). Acta Physiol Plant. 2017;39:222.

11. Ning CL, Jiang Y, Meng JS, Zhou CH, Tao J. Herbaceous peony seed oil: a rich source of unsaturated fatty acids and γ-tocopheral. Eur J Lipid Sci Tech. 2015;117(4):532-42.

12. Thandapilly SJ, Raj P, Louis XL, Perera D, Yamanagedara P, Zahradka P, Taylor CG, Netticadan T. Canola oil rich in oleic acid improves diastolic heart function in diet-induced obese rats. J Physiol Sci. 2017;67:425-30.

13. Patterson E, Wall R, Fitzgerald GF, Ross RP, Stanton C. Health implications of high dietary omega-6 polyunsaturated fatty acids. J Nutr Metab, 2012;doi:10.1155/2012/539426.

14. Kulkarni KP, Kim M, Song JT, Bilyeu KD, Lee JD. Genetic improvement of the fatty acid biosynthesis system to alter the ω-6/ω-3 ratio in the soybean seed. J Am Oil Chem Soc. 2017;94:1403-10.

15. Domenichiello AF, Kitson AP, Chen CT, Trépanier MO, Stavro PM, Bazinet RP. The effect of linoleic acid on the whole body synthesis rates of polyunsaturated fatty acids from α-linolenic acid and linoleic acid in free-living rats. J Nutr Biochem. 2016;30:167-76.

16. Simopoulos AP. The importance of the ratio of omega-6/mogea-3 essential fatty acids. Biomed Pharmacot. 2002;56:365-79.

17. Li SS, Wang LS, Shu QY, Wu J, Chen LG, Shao S, Yin DD. Fatty acid composition of developing tree peony (Paeonia section Moutan DC.) seeds and transcriptome analysis during seed development. BMC Genomics. 2015;16:208.

18. Goettel W, Ramirez M, Upchurch RG, An YC. Identification and characterization of large DNA deletions affecting oil quality traits in soybean seeds through transcriptome sequencing analysis. Theor Appl Genet. 2016;129:1577-93.

19. Yin DM, Wang Y, Zhang XG, Li HM, Lu X, Zhang JS, Zhang WK, Chen SY. De Novo assembly of the peanut (Arachis hypogaea L.) seed transcriptome revealed candidate unigenes for oil accumulation pathways. PLOS ONE. 2013;doi:org/10.1371/journal.pone.0073767.

20. Bancroft I, Morgan C, Fraser F, Higgins J, Wells R, Clissold L, Baker D, Long Y, Meng JL, Wang XW, Liu SY, Trick M.(2011) Dissecting the genome of the polyploidy crop oilseed rape by transciptome sequencing. Nat Biotechnol. 2011;29:762-6.

21. Livaja M, Wang Y, Wieckhorst S, Haseneyer G, Seidel M, Hahn V, Knapp SJ, Taudien S, Schön CC, Bauer E. BSTA: a targeted approach combines bulked segregant analysis with next-generation sequencing and de novo transcriptome assembly for SNP discovery in sunflower. BMC Genomics. 2013;14:628.

22. Bazakos C, Manioudaki ME, Therios Io, Voyiatzis D, Kafetzopouls D, Awada T, Kalaitzis P. Comparative transcriptome analysis of two olive cultivars in response to NaCl-stress. PLOS ONE. 2012;doi:10.1371/journal.pone.0042931.

23. Kim HUK, Lee KR, Shim D, Lee JH, Chen GQ, Hwang S. Transcriptome analysis and identification of genes associated with ω-3 fatty acid biosynthesis in Perilla frutescens (L.) var. frutescens. BMC Genomics. 2016;17:474.

24. Dussert S, Guerin C, Andersson M, Joët T, Tranbarger TJ, Pizot M, Sarah G, Omore A, Durand-Gasselin T, Morcillo F. Comparative transcriptome analysis of three oil palm fruit and seed tissues that differ in oil content and fatty acid composition. Plant Physiol. 2013;162:1337-58.

25. Xia EH, Jiang JJ, Huang H, Zhang LP, Zhang HB, Gao LZ. (2014) Transcriptome analysis of the oil-rich tea plant, Camellia oleifera, reveals candidate genes related to lipid metabolism. PLOS ONE. 2014;doi:10.1371/journal.pone.0104150.

26. Feng YZ, Wang L, Fu JM, Wuyun TN, Du HY, Tan XF, Zou F, Li FD. Transcriptome sequencing discovers genes related to fatty acid biosynthesis in the seeds of Eucommia ulmoides. Genes Genom. 2016;38:275-83.

27. Fatima T, Snyder CL, Schroeder WR, Cram D, Datla R, Wishart D, Weselake RJ, Krishna P. Fatty acid composition of developing sea buckthorn (Hippophae rhamnoides L.) berry and the transcriptome of the mature seed. PLOS ONE. 2012;doi:10.1371/journal.pone.0034099.

28. Fu QT, Niu LJ, Chen MS, Tao YB, Wang XL, He HY, Pan BZ, Xu ZF. De novo transcriptome assembly and comparative analysis between male and benzyladenine-induced female inflorescence buds of Plukenetia volubilis. J Plant Physiol. 2018;221:107-18.

29. Hao ZJ, Wei MR, Gong SJ, Zhao DQ, Tao Jun. Transcriptome and digital gene expression analysis of herbaceous peony (Paeonia lactiflora Pall.) to screen thermo-tolerant related differently expressed genes. Genes Genom. 2016;38:1201-15.

30. Wu YQ, Zhao DQ, Tao J. Analysis of codon usage patterns in herbaceous peony (Paeonia lactiflora Pall.) based on transcriptome data. Genes. 2015;6(4):1125-39.

31. Zhao DQ, Jiang Y, Ning CL, Meng JS, Lin SS, Ding W, Tao J. Transcriptome sequencing of a chimaera reveals coordinated expression of anthocyanin biosynthetic genes mediating yellow formation in herbaceous peony (Paeonia lactiflora Pall.). BMC Genomics. 2014;15:689.

32. Grahherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowhury R, Zeng QD, Chen ZH, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechno. 2011;29:644-52.

33. Pertea G, Huang XQ, Liang F, Antonescu V, Sultana R, Karamycheva S, Lee Y, White J, Cheung F, Parvizi B, Tsai J, Quackenbush J. TIGR gene indices clustering tools (TGICL): a software system for fast clustering of large EST datasets. Bioinformatics. 2013;19:651-2.

34. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403-10.

35. Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674-6.

36. Ohlrogge J, Browse J. Lipid biosynthesis. Plant Cell. 1995;7:957-70.

37. Liu ZG, Wang LM, Wang HL, Liu LJ. Recent advances in understanding the effects of polyunsaturated fatty acids on brain function. Food Sci. 2015;36(21):284-90. (In Chinese)

38. Bougnoux P, Maillard V, Chajès V. Omega-6/omega-3 polyunsaturated fatty acid ratio and breast cancer. World Rev Nutr Diet. 2005;94:158-65.

39. Russo GL. Dietary n-6 and n-3 polyunsaturated fatty acids: From biochemistry to clinical implications in cardiovascular prevention. Biochem Pharmacol. 2009;77:937-46.

40. Jin M, Lu Y, Pang TT, Zhu TT, Yuan Y, Sun P, Zhou F, Ding XY. Effects of dietary n-3 LC-PUFA/n-6 C18 PUFA ratio on growth, feed utilization, fatty acid composition and lipid metabolism related gene expression in black seabream, Acanthopagrus schlegelil. Aquaculture. 2019;500:521-31.

41. Okuley J, Lightner J, Feldmann K, Yadav N, Lark E, Browse J. Arabidopsis FAD2 gene encodes the enzyme that is essential for polyunsaturated lipid synthesis. Plant Cell. 1994;6:147-58.

42. Hitz WD, Carlson TJ, Booth JR, Kinney AJ, Stecca KL, Yadav NS. Cloning of a higher-plant plastid ω-6 fatty acid desaturase cDNA and its expression in a cyanobacterium. Plant Physiol. 1994;105:635-41.

43. Iba K, Gibson S, Nishiuchi T, Fuse T, Nishimura M, Arondel V, Hugly S, Somerville C. A gene encoding a chloroplast ω-3 fatty acid desaturase complements alterations in fatty acid desaturation and chloroplast copy number of the fad7 mutant of Arabidopsis thaliana. J Biol Chem. 1993;268:24099-105.

44. Huang AHC. Oil bodies and oleosins in seeds. Annu Rev Plant Physiol Plant Mol Biol. 1992;43:177-200.

45. Ding J. Transcript expression analysis involved in oil synthesis and accumulation of Hippophae rhamnoides pulp and seed. Haerbin: Northeast Forestry University; 2016.(In Chinese)

46. Natarajan S, Kim JK, Jung TK, Doan TTN, Ngo HPT, Hong MK, Kim S, Tan VP, Ahn S, Lee SH, Han Y, Ahn YJ, Kang LW. Crystal strcture of malonyl CoA-Acyl carrier protein transacylase from Xanthomanous oryzae pv. oryzae and its proposed binding with ACP. Mol Cells. 2012;33:19-25.

47. Simon JW, Slabas AR. cDNA cloning of Brassica napus malonyl-CoA:ACP transacylase (MCAT) (fab D) and complementation of an E. coli MCAT mutant. FEBS Lett. 1998;435: 204-6.

48. White SW, Zheng J, Zhang YM, Rock CO. The structural biology of type II fatty acid biosynthesis. Annu Rev Biochem. 2005;74:791-831.

49. González-Mellado D, Westtstein-Knowles P, Garcés R, Martínez-Force E. The role of β-ketoacyl-acyl carrier protein synthase III in the condensation steps of fatty acid biosynthesis in sunflower. Planta. 2010;231:1277-89.

49. Yan FF, Tan XF, Long HX. Cloning and sequence analysis of full-length cDNA encoding beta-ketoacyl-ACP synthase III from Vernicia fordii. Acta Agri Univ Jiangxiensis. 2013;35(4):775-81.(In Chinese)

51. Li J, Li MR, Wu PZ, Tian CE, Jiang HW, Wu GJ. Molecular cloning and expression analysis of a gene encoding a putative beta-ketoacyl-acyl carrier protein (ACP) synthase III (KAS III) from Jatropha curcas. Tree Physiol. 2008;28:921-7.

52. Choi KH, Heath RJ, Rock CO. β-ketoacyl-acyl carrier protein synthase III (FabH) is a determining factor in branched-chain fatty acid biosynthesis. J Bacteriol. 2000;182:365-70.

53. Abbadi A, Brummel M, Schütt BS, Slabaugh MB, Schuch R, Spener F. Reaction mechanism of recombinant 3-oxoacyl-(acyl-carrier-protein) synthase III from Cuphea wrightii embryo, a fatty acid synthase tyepe II condensing enzyme. Bioche J. 2000;345:153-60.

54. Brück FM. Brummel M, Schuch R, Spener F. In-vitro evidence for feed-back regulation of β-ketoacyl-acyl carrier protein synthase III in medium-chain fatty acid biosynthesis. Planta. 1996;198:271-8.

55. Dehesh K, Tai H, Edwards P, Byrne J, Jaworski JG. Overexpression of 3-ketoacyl-acyl-carrier protein synthase IIIs in plants reduces the rate of lipid synthesis. Plant Physiol. 2001;125:1103-14.

56. Ohlrogge JB, Jaworski JG. Regulation of fatty acid synthesis. Ann Rev Plant Bio. 1997;48:109-36.

57. Gibson S, Falcone DL, Browse J, Somerville C. Use of transgenic plants and mutants to study the regulation and function of lipid composition. Plant Cell Environ. 1994;17:627-37.

58. Minoru K, Michihiro A, Susumu G, Masahiro H, Mika H, Masumi I, Toshiaki K, Shuichi K, Shujiro O, Toshiaki T, Yoshihiro Y. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008;36:D480-4.

59. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.

60. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

61. Zhao DQ, Tao J, Han CX, Ge JT. Actin as an alternative internal control gene for gene expression analysis in herbaceous peony (Paeonia lactiflora Pall.). Afr J Agric Res. 2012;7:2153-2159.

62. Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative CT method. Nat Protoc. 2008;3:1101-8.

Figures, Tables and Additional Files

Fig 1

Fatty acid composition and content of seeds of Paeonia lactiflora 'Hangshao'

A: Fatty acid ion chromatography; B: Content of main fatty acid during developing seed

1: PA; 2: SA; 3: OA; 4:LA; 5: ALA; S1: 30 DAF; S2: 45 DAF; S3: 60 DAF; S4: 75 DAF; S5: 90 DAF

Fig 2

Functional distribution of Unigenes annotaed for seeds of Paeonia lactiflora 'Hangshao'

A: Distribution of NR annotated species

B: Functional distribution of GO annotation

C: Functional distribution of KOG annotation

D: Functional distribution of KEGG annotation

Fig 3

Herbaceous peony sequence associated with de novo fatty acid biosynthesis, fatty acid elongation, and biosynthesis of unsaturated fatty acid pathway.

Fig 4

Analysis of DEGs for seeds of Paeonia lactiflora 'Hangshao'

A: Gene analysis of DEGs in diffrent stage; B: Cluster analysis of DEGs in different stage

Fig 5

Analysis of DEGs obtained from intersection of three groups

A: Venn diagram between DEGs from three groups; B: Functional distribution of intersection DEGs annotaed GO; C: Functional distribution of intersection DEGs annotaed KEGG

Fig 6

qRT-PCR validations of expression levels of DEGs related to fatty acid biosynthesis for developmental seed of Paeonia lactiflora 'Hangshao'

The column chart showed the qRT-PCR validations in 5 periods; the broken line graph showed the transcriptome sequencing in 3 priods

S1: 30 DAF; S2: 45 DAF; S3: 60 DAF; S4: 75 DAF; S5: 90 DAF

Fig 7

Seeds of collecting period from Paeonia lactiflora 'Hangshao'

S1: 30 DAF; S2: 45 DAF; S3: 60 DAF; S4: 75 DAF; S5: 90 DAF

Table 1

Quality metrics of transcripts of seeds of Paeonia lactiflora 'Hangshao'

Table 2

Functional annotation of uingenes for seeds of Paeonia lactiflora 'Hangshao' in public functional database

Table 3

Number of Unigenes related to fatty acid biosynthesis in seeds of Paeonia lactiflora 'Hangshao'

Additional files

Additional file 1: Clean reads quality metircs

Additional file 2: Length distribution of unigenes for seeds of Paeonia lactiflora 'Hangshao'

Additional file 3: Number of DEGs for GO classification

Additional file 4: Functional distribution of DEGs annotaed GO for seeds of Paeonia lactiflora 'Hangshao'

Additional file 5: Number of DEGs for KEGG annotation

Additional file 6: DEGs for KEGG annotation

Additional file 7: Differentially expressed genes related to fatty acid biosynthesis in herbaceous peony 'Hangshao' seeds

Additional file 8: Gene-specific primers sequence for for qRT-PCR detection