All relevant data are within the paper and its Supporting Information files.

Introduction {#sec001}
============

The angiosperms are the most diverse group of land plants with the number of species in the range of 250,000 to 400,000 \[[@pone.0154882.ref001]\]. They have dramatic differences not only in organ morphology (leaves, flowers, seeds, roots, and vascular tissues) but also in the chemical composition of the seed \[[@pone.0154882.ref002]\]. Interestingly, oil plants, such as peanut, sesame and soybean, are generally rich in oil of their seeds, while most cereals like rice, wheat and sorghum specifically accumulate starch and have a relatively low oil content in the seed. Since the formation of seed oil and starch is both dependent on the supply of photosynthetic carbon \[[@pone.0154882.ref003], [@pone.0154882.ref004]\], it is likely that there are divergent mechanisms in the evolution of these high- and low-oil plants to regulate the partitioning of carbon between oil and other storage products. Previous studies have demonstrated that increasing the carbon flow to lipid biosynthesis can significantly increase the seed oils \[[@pone.0154882.ref004]--[@pone.0154882.ref006]\]. However, the mechanism by which more carbohydrates flow to *de novo* fatty acid (FA) synthesis in high-oil content plants is unclear.

Seed oil is not only the major source of carbon and energy for germination and seedling growth but also provides humans with renewable sources of food, biofuel and industrial raw materials. Up to now acyl-lipid metabolism in Arabidopsis has been well studied and more than 600 genes have been predicted to encode the enzymes and regulatory factors associated with this process \[[@pone.0154882.ref007], [@pone.0154882.ref008]\]. Among these predicted genes, some have been shown to be associated with changes in seed oil accumulation \[[@pone.0154882.ref006]\]. The over-expression of many individual key enzymes altered seed oil content in various plants. For example, ACCase in potato \[[@pone.0154882.ref009]\], *Brassica napus* \[[@pone.0154882.ref010]\] and *Escherichia coli* \[[@pone.0154882.ref011]\]; glycerol-3-phosphate dehydrogenase (GPDH) in *B*. *napus* \[[@pone.0154882.ref012]\]; glycerol-3-phosphate acyltransferase (GPAT) in *Arabidopsis* \[[@pone.0154882.ref013]\]; 2-lysophosphatidic acid acyltransferase (LPAAT) in *Arabidopsis* \[[@pone.0154882.ref014]\]; and acyl-CoA: diacylglycerol acytransferase (DGAT) in *Arabidopsis* \[[@pone.0154882.ref015], [@pone.0154882.ref016]\], *Glycine max* \[[@pone.0154882.ref017]\], *B*. *napus* \[[@pone.0154882.ref018]\] and maize \[[@pone.0154882.ref019]\]. In addition to these key enzymes that participate in lipid synthesis, the expression of transcription factors (TFs) directly or indirectly regulating genes involved in carbohydrate and lipid metabolism can also evidently change seed oil content; such TFs include *WRINKLED1* (*WRI1*) \[[@pone.0154882.ref020]\], *LEAFY COTYLEDON1* (*LEC1*) \[[@pone.0154882.ref021], [@pone.0154882.ref022]\], *LEAFY COTYLEDON2* (*LEC2*) \[[@pone.0154882.ref023]\], *FUSCA3* (*FUS3*) \[[@pone.0154882.ref024]\], *GmbZIP123* \[[@pone.0154882.ref025]\], *GmMYB73* \[[@pone.0154882.ref026]\] and *ABSCISIC ACID INSENSITIVE3* (*ABI3*) \[[@pone.0154882.ref027], [@pone.0154882.ref028]\]. However, each of the above studies only focused on a single enzyme or TF involved in lipid metabolism. In reality, seed oil content is affected by multiple genes \[[@pone.0154882.ref029], [@pone.0154882.ref030]\] or multiple reactions \[[@pone.0154882.ref031]\]. For example, specific combination of expression of *WRI1*, *DGAT* and triacylglycerol lipase *SUGAR-DEPENDENT1* resulted in a higher percentage seed oil content than that obtained by manipulation of each gene individually \[[@pone.0154882.ref032]\]. More importantly, the seed oil content is influenced by multiple metabolic pathways, such as sucrose catabolism, glycolysis, pentose phosphate pathway, and related pathways which rely on the supply of carbon \[[@pone.0154882.ref033]--[@pone.0154882.ref038]\]. Therefore, it is necessary to consider the involvement of multiple genes and the interaction of multiple related pathways in order to understand the mechanism of high-oil content in high-oil plants. The recent sequencing of many plant genomes has provided an opportunity to investigate this mechanism using comparative genome and transcriptome analyses.

Analysis of genome-wide differential gene expression between developmental stages, and between sub-species (like 输入文字或网址, 即可翻译 cultivated and wild forms) could provide insights into biological pathways and molecular mechanisms that regulate seed development and nutrient accumulation. Seed development is also an important part of the reproductive (R) process in flowering plants. In soybean, this reproductive process is divided into eight stages from flowering (R1) to full maturity (R8) \[[@pone.0154882.ref039]\]. During the R4 to R7 stages, importantly, seeds grow rapidly, accumulating nutrients, lipids and storage proteins.

In this study, all the genes from three low-oil grasses (*Sorghum bicolor*, *Setaria italica*, and *Oryza sativa*) and four high-oil dicots (*Glycine max*, *Gossypium raimondii*, *Ricinus communis*, and *Arabidopsis thaliana*) were clustered in order to obtain a list of high-oil dicot-specific genes. A gene ontology (GO) enrichment analysis and a pathway level co-expression (PLC) network analysis were then conducted to identify genes or TFs that are likely to be associated with oil accumulation. Analyses of gene expression during the various stages of seed development in soybean and of RNA sequencing (RNA-seq) differential expression between high- and low-oil content soybean accessions were performed to identify differentially expressed genes (DEGs) in the core pathways associated with the deposition of seed oil. These results were used to further investigate how evolutionary divergence contributes to differences in seed oil content between the two kinds of plants and to discover new genes associated with the seed oil differences.

Results {#sec002}
=======

Identification and GO enrichment analysis of dicot-specific genes {#sec003}
-----------------------------------------------------------------

In this study, OrthoMCL was applied to construct potential orthologous groups (OGs) of proteins across four high-oil dicots (seed oil content: 20--50%) and three low-oil grasses (\<3%) ([S1 Table](#pone.0154882.s009){ref-type="supplementary-material"}), because it can group both orthologs and paralogs over multiple eukaryotic taxa by using a Markov Cluster algorithm (MCL) \[[@pone.0154882.ref040]\]. As a result, all the 252,443 genes from the above seven species were clustered into 29,095 OGs ([S1 Dataset](#pone.0154882.s001){ref-type="supplementary-material"}). Among these gene families, only 1,534 (5.27%) OGs appear to be specific to the high-oil dicot lineage and are defined as high-oil dicot-specific clusters since all the genes in these OGs come from all the four rosid species but not from any grasses ([S1 Dataset](#pone.0154882.s001){ref-type="supplementary-material"}). Note that OrthoMCL clusters proteins based on overall conservation but not on individual protein domains. Thus, high-oil dicot-specific OGs in this study contain families specific in high-oil dicots or families with much low similarities between dicots and grasses. Among the 1,534 high-oil dicot-specific clusters, 4,051, 2,758, 1,731 and 2,152 genes were found in *G*. *max*, *G*. *raimondii*, *R*. *communis*, and *A*. *thaliana*, respectively.

To understand the functions of these high-oil dicot-specific genes, a GO enrichment analysis was conducted for 4,051 soybean genes compared with all the annotated genes. 392 GO terms for biological processes, molecular functions and cellular components were identified; and these were distributed in 77 GO slim terms ([S2 Dataset](#pone.0154882.s002){ref-type="supplementary-material"}).

Among the 43 GO slims for biological processes, some were involved in metabolic pathways ([S2 Dataset](#pone.0154882.s002){ref-type="supplementary-material"}), such as biosynthetic process, carbohydrate metabolic process, catabolic process, generation of precursor metabolites and energy, lipid metabolic process, metabolic process, protein metabolic process and secondary metabolic process. In addition, slim GO:0006810 (transport) was involved in transport of many intermediates of carbohydrate degradation and glycolytic pathways, which includes triose phosphate transmembrane transport, phosphoenolpyruvate transport, acylglycerol transport, glucose-6-phosphate transport, phosphoglycerate transport, hexose phosphate transport, regulation of intracellular transport and triose phosphate transport.

Expression patterns of dicot-specific genes during seed development {#sec004}
-------------------------------------------------------------------

The transcriptomic data for soybean seeds at seven stages of development \[[@pone.0154882.ref041]\] were used to analyze the expression patterns of 4,051 dicot-specific genes. 3,155 (77.88%) genes were expressed in developing seeds across two biological replicates. Among these 3,155 genes, 3,150 (99.84%) were grouped into eight clusters based on Pearson's correlation coefficients, implemented by MCL \[[@pone.0154882.ref042]\]. All the eight clusters are shown in [Fig 1](#pone.0154882.g001){ref-type="fig"} and the genes in each cluster are listed in [S1 Dataset](#pone.0154882.s001){ref-type="supplementary-material"}.

![Average expression profile (*y* axis) for each cluster of dicot-specific genes at seven stages of soybean seed development (*x* axis).\
RPKM: reads per kilobase of gene model per million mapped reads; DAF: days after flowering. The dataset was downloaded from <http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42871>.](pone.0154882.g001){#pone.0154882.g001}

Identification of genes related to seed oil content in soybean {#sec005}
--------------------------------------------------------------

### Identification of high-oil dicot-specific genes directly participating in acyl-lipid metabolism {#sec006}

Compared with acyl-lipid metabolism genes in *Arabidopsis* \[[@pone.0154882.ref007], [@pone.0154882.ref008]\], 1,123 orthologous genes were identified in soybean ([S1 Dataset](#pone.0154882.s001){ref-type="supplementary-material"}). Among these 1,123 genes, 54 were high-oil dicot-specific genes and distributed in the above eight clusters ([S2 Table](#pone.0154882.s010){ref-type="supplementary-material"}). In particular, almost all the genes coding key enzymes of fatty acid synthesis like Biotin Carboxylase (BC), Biotin Carboxyl Carrier Protein 2 (BCCP2), Carboxyltransferase alpha subunit (α-CT), WRI1, and FUS3 were in cluster 3.

### Dicot-specific lipid-metabolism-related genes predicted from PLC analysis {#sec007}

Genes in the clusters 3 to 5 show highest expression levels at one of the stages from 5--6 mg to 400--500 mg, which tend to occur between R4 and R7. During these stages, lipids are rapidly accumulated in the seed \[[@pone.0154882.ref039]\]. To ensure prediction accuracy, these genes in clusters 3 to 5 were selected for further analysis. This is because that the genes in cluster 2 show highest expression at the dry whole seed stage and the ones in the remaining clusters are down-regulated during seed development. In the clusters 3 to 5, there were 828 genes ([S1 Dataset](#pone.0154882.s001){ref-type="supplementary-material"}), among which 23 genes have been determined to participate in the above lipid synthesis ([S2 Table](#pone.0154882.s010){ref-type="supplementary-material"}). GO enrichment analysis for the remaining 805 genes assigned 207 (25.71%) genes to the biological processes that may be associated with carbohydrate and lipid metabolic pathways ([S3 Dataset](#pone.0154882.s003){ref-type="supplementary-material"} and [S3 Table](#pone.0154882.s011){ref-type="supplementary-material"}), such as carbohydrate metabolic process, lipid metabolic process, transport, signal transduction, and catabolic process.

To further identify and prioritize novel candidate acyl-lipid metabolism members in these three clusters, we conducted a pathway-level coexpression (PLC) network analysis \[[@pone.0154882.ref043]--[@pone.0154882.ref046]\] between the above 207 genes and 1,123 genes in acyl-lipid metabolism. The results show that 93 dicot-specific genes are candidate genes involved in acyl-lipid metabolism pathways ([Table 1](#pone.0154882.t001){ref-type="table"}, [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"} and [S5 Dataset](#pone.0154882.s005){ref-type="supplementary-material"}).

10.1371/journal.pone.0154882.t001

###### Dicot-specific genes associated with acyl-lipid metabolism.

![](pone.0154882.t001){#pone.0154882.t001g}

  GO slims                                     Dicot-specific genes                                                                                                                                                                                                                                                                                                                                                                                                                                          Pathway annotation                       Pfam annotation                                                     Enzymes/Proteins                                            Reference
  -------------------------------------------- ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ---------------------------------------- ------------------------------------------------------------------- ----------------------------------------------------------- ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  GO:0005975; Carbohydrate metabolic process   *Glyma0165s00200*, *Glyma01g28520*, *Glyma03g08860*                                                                                                                                                                                                                                                                                                                                                                                                                                                    PF00686                                                             *CBM20-like* starch binding domian                          Marchler-Bauer et al. \[[@pone.0154882.ref067]\]; Southall et al. \[[@pone.0154882.ref070]\]; Rodriguez-Sanoja et al. \[[@pone.0154882.ref069]\]; Christiansen et al. \[[@pone.0154882.ref068]\]
                                               ***Glyma09g04330***, *Glyma15g15360*                                                                                                                                                                                                                                                                                                                                                                                                                          PWY-6902; chitin degradation II          PF00182                                                             chitinase                                                   
                                               *Glyma01g33440*, *Glyma03g03400*, *Glyma03g03460*, *Glyma06g47690*, *Glyma06g47690*, *Glyma12g00700*                                                                                                                                                                                                                                                                                                                                                          PWY-1081; homogalacturonan degradation   PF01095, PF04043                                                    pectinesterase                                              
                                               ***Glyma04g02660*, *Glyma06g02690***, *Glyma14g40400*, ***Glyma17g37750*, *Glyma17g37760***                                                                                                                                                                                                                                                                                                                                                                                                            PF02704                                                             Gibberellin-regulated family protein                        Aubert et al. \[[@pone.0154882.ref097]\]; Chen et al. \[[@pone.0154882.ref098]\]
                                               *Glyma01g37150*, ***Glyma11g08120***, *Glyma15g08960*, *Glyma05g36680*, *Glyma14g40110*, ***Glyma02g34870***, *Glyma03g36250*, ***Glyma10g10530***, ***Glyma19g38900***, *Glyma02g36150*, *Glyma10g08740*, *Glyma08g16880*, *Glyma09g01520*, *Glyma17g36200*, *Glyma14g08970*, *Glyma10g16100*, *Glyma07g05970*, *Glyma04g40630*, *Glyma06g14160*, *Glyma05g06410*, *Glyma19g07830*, *Glyma18g00560*, *Glyma02g47500*, ***Glyma14g01260***, *Glyma17g09260*                                            [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}                                                               
  GO:0006810; Transport                        *Glyma07g38830*, *Glyma13g27680*, ***Glyma15g11270***                                                                                                                                                                                                                                                                                                                                                                                                                                                  PF03151, PF00892                                                    *GPT2*                                                      Kammerer et al. \[[@pone.0154882.ref071]\]; Niewiadomski et al. \[[@pone.0154882.ref074]\]; Andriotis et al. \[[@pone.0154882.ref073]\]; Kunz et al. \[[@pone.0154882.ref075]\]; Bourgis et al. \[[@pone.0154882.ref077]\]
                                               *Glyma18g08740*                                                                                                                                                                                                                                                                                                                                                                                                                                                                                        PF03151                                                             triosephosphate translocator subfamily protein              
                                               *Glyma12g13070*, ***Glyma08g02510***, *Glyma03g05880*, *Glyma15g17310*, *Glyma06g40980*, *Glyma07g12460*, *Glyma09g06260*, *Glyma12g34020*, *Glyma18g14810*, *Glyma09g04870*, ***Glyma15g15990***, *Glyma04g06230*, *Glyma03g33150*, *Glyma17g09640*, ***Glyma05g37050***, *Glyma09g08410*, *Glyma12g33880*, *Glyma19g34820*, *Glyma03g33660*, *Glyma06g22240*, *Glyma14g33680*, *Glyma20g23120*, *Glyma17g09260*                                                                                      [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}                                                               
  GO:0006629; Lipid metabolic process          *Glyma02g45680*, *Glyma14g03130*                                                                                                                                                                                                                                                                                                                                                                                                                                                                       PF00067                                                             cytochrome P450, family 718                                 Pinot and Beisson \[[@pone.0154882.ref099]\]; Sun et al. \[[@pone.0154882.ref100]\]
                                               *Glyma18g45250*                                                                                                                                                                                                                                                                                                                                                                                                                                               PWY-2761; glyceollin biosynthesis I      PF01370                                                             pterocarpin synthase                                        Kim et al. \[[@pone.0154882.ref101]\]
                                               ***Glyma02g46220***, *Glyma14g02510*                                                                                                                                                                                                                                                                                                                                                                                                                                                                   [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}                                                               
  GO:0007165; Signal transduction              *Glyma09g07090*, ***Glyma15g18380***, *Glyma17g06290*                                                                                                                                                                                                                                                                                                                                                                                                                                                  PF00320                                                             GATA type zinc finger transcription factor family protein   Bi et al. \[[@pone.0154882.ref096]\]; Velmurugan et al. \[[@pone.0154882.ref095]\]
                                               *Glyma09g39570*, *Glyma10g12130*                                                                                                                                                                                                                                                                                                                                                                                                                              PWY-5035; gibberellin biosynthesis III   PF03171                                                             gibberellin 3β-dioxygenase                                  Chen et al. \[[@pone.0154882.ref098]\]
                                               *Glyma12g32910*, ***Glyma13g39340***, *Glyma02g47500*, *Glyma14g01260*, *Glyma13g29070*, *Glyma15g09980*                                                                                                                                                                                                                                                                                                                                                                                               [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}                                                               
  GO:0019538; Protein metabolic process        *Glyma05g23620*, *Glyma17g16690*, *Glyma01g03580*, *Glyma06g10110*, *Glyma08g39290*, *Glyma17g12200*, *Glyma18g19720*                                                                                                                                                                                                                                                                                                                                                                                  [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}                                                               

High-oil dicot-specific genes were predicted by PLC analysis. 17 genes with bold type were also differentially expressed in the high- and low-oil soybean seeds at the 0.01 significant level. All the information containing both annotation and P-value in differential expression analysis for each gene is given in [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}.

### Expression patterns of genes encoding core lipid synthetic enzymes {#sec008}

In soybean, the core pathways for accumulating seed oil operate through FA synthesis and the export of FAs from the plastid followed by triacylglycerol synthesis, and oil body formation, which included 156 genes ([S1 Dataset](#pone.0154882.s001){ref-type="supplementary-material"}). Among these genes, 113 (72.44%) had higher expression level than the average in at least one of the stages from 5--6 mg to 400--500 mg ([S4 Dataset](#pone.0154882.s004){ref-type="supplementary-material"}), indicating that most key genes for lipid synthesis were up-regulated during oil accumulation in *G*. *max*. This phenomenon was similar to that in *Arabidopsis* \[[@pone.0154882.ref047]\]. If the criterion used was twice average level, 40 (25.6%) genes were identified ([S4 Dataset](#pone.0154882.s004){ref-type="supplementary-material"}). These 40 genes were largely distributed in the initiation of FA synthesis, triacylglycerol synthesis, and oil-body formation ([Fig 2](#pone.0154882.g002){ref-type="fig"}).

![The distribution of 40 highly up-regulated genes and the respective DEGs in the framework of acyl-lipid metabolism.\
Enzymes with blue color are highly up-regulated and 40 genes coding them are listed. Among these 40 highly up-regulated genes, DEGs between the high- and low-oil soybean accessions are highlighted with red color. ABI, ABSCISIC ACID INSENSITIVE; ACP, acyl carrier protein; BCCP, biotin carboxyl carrier protein; DGAT, acyl-CoA: diacylglycerol acytransferase; DHAP, Dihydroxyacetone phosphate; DHLAT, dihydrolipoamide acetyltransferase; ENR, enoyl-ACP reductase; ER, endoplasmic reticulum; FAS, fatty acid synthase; FATA(B), fatty acylthioesteraseA(B); FUS, FUSCA; G3P, glycerol-3-phosphate; GPAT, glycerol-3-phosphate acyltransferase; GPDH, glycerol-3-phosphate dehydrogenase; HAD, hydroxyacyl-ACP dehydrase; KAR, ketoacyl-ACP reductase; KAS, ketoacyl-ACP synthase; LACS, long-chainacyl-CoA synthetase; LPD, dihydrolipoamide dehydrogenase; MCMT, malonyl-CoA:ACP malonyltransferase; OLE, oleosins; PDAT, phospholipid:diacylglycerolacyl transferase; PDH, pyruvatede hydrogenase; PP, phosphatidate phosphatase; SAD, stearoyl-ACP desaturase; STERO, steroleosin; TAG, triacylglycerol; WRI, WRINKLED.](pone.0154882.g002){#pone.0154882.g002}

Fourteen of the above 40 genes were observed in FA synthesis, and 10 of the 14 genes were involved in the initiation of FA synthesis. Among the 10 genes, five encoded two subunits of the plastidial pyruvate dehydrogenase complex (PDHC) that promotes pyruvate decarboxylation to acetyl-CoA and is the key enzyme linking carbohydrate metabolism to FA synthesis \[[@pone.0154882.ref004]\]; and the others encoded two subunits of ACCase. Of these 40 genes, two genes, those encoding ACP4 and LACS9, were involved in FA transportation. Eight of the above 40 genes were involved in TAG assembly; this group included 2, 3, 2, and 1 genes respectively coding phosphatidate phosphatase (PP), DGAT1, phospholipid: diacylglycerol acyltransferase (PDAT1), and PDAT2; these enzymes catalyze the consecutive steps after the second acylation of glycerol-3-phosphate \[[@pone.0154882.ref015], [@pone.0154882.ref019], [@pone.0154882.ref048]\]. Three and six of the above 40 genes, respectively encoding steroleosins (STEROs) and oleosins (OLEs), were observed in oil-body formation, and they were up-regulated, mostly by some hundredfold, during stages from 5--6 mg to 400--500 mg in seeds ([S4 Dataset](#pone.0154882.s004){ref-type="supplementary-material"}). The remaining genes in this group of 40 genes encoded the TFs WRI1, VAL1, FUS3, and ABI3. Kim et al. \[[@pone.0154882.ref049]\] divided Arabidopsis oleosin genes into three groups on the basis of their tissue-specific expression. In this study, all the genes in the S type (genes expressed only in maturing seed (siliques)) and SM type (genes expressed in both maturing seeds and florets (microspores)) were used to conduct a phylogenetic analysis ([Fig 3](#pone.0154882.g003){ref-type="fig"}). The results showed distinct differences of OLEs in patterns of gene duplication and loss between high-oil dicots and low-oil grasses. Genes coding oleosins were preferentially retained in high-oil dicots. In two angiosperm groups, there was only one copy in each grass species, while there were many duplicates in dicots. Remarkably, two high-oil dicot-specific groups were identified, which were resulted from gene loss in grasses after whole genome duplication of ancestor angiosperm. In dicot-specific group 1, *Glyma04g08220* and *Glyma06g08290* were up-regulated by a thousand-fold during stages from 5--6 mg to 400--500 mg ([S4 Dataset](#pone.0154882.s004){ref-type="supplementary-material"}) and may play an important role in soybean seed oil accumulation. In dicot-specific group 2, *At3G01570*, *At3G27660* and *At5G40420* in Arabidopsis were shown to be S type genes \[[@pone.0154882.ref049]\]. Meanwhile, *Glyma05g08880* and *Glyma19g00400* in this group were highly up-regulated in hundred-fold in the dry whole seed ([S4 Dataset](#pone.0154882.s004){ref-type="supplementary-material"}), indicating that this dicot-specific group may perform some specific functions in seed development. Furthermore, highly up-regulated and abundant OLEs probably play major roles in oil accumulation and/or oil body development, and are associated with high-oil content in seeds \[[@pone.0154882.ref050], [@pone.0154882.ref051]\].

![Phylogenetic tree of oil body protein OLE.\
The phylogenetic tree was constructed using Neighbor-Joining method. The numbers on the branches represent the bootstrap support, and square boxes indicate duplication events. S indicates genes expressed only in maturing seed and SM in both maturing seeds and florets. The highly up-regulated genes in [Fig 2](#pone.0154882.g002){ref-type="fig"} are highlighted in red color.](pone.0154882.g003){#pone.0154882.g003}

Interestingly, *OLEs*, *STEROs* and *ABI3* genes are all significantly over-expressed and co-expressed with each other ([S1 Fig](#pone.0154882.s006){ref-type="supplementary-material"} and [S4 Dataset](#pone.0154882.s004){ref-type="supplementary-material"}). Since *OLEs* were regulated by *ABI3* \[[@pone.0154882.ref028]\], *ABI3* might also up-regulate the expression of *STEROs*. To improve seed oil content in molecular breeding, it may be beneficial to increase the expression of *OLEs*, *STEROs* or *ABI3* genes.

### Identification of candidate TF genes related to seed oil accumulation {#sec009}

Several important transcription factors (TFs) have been found to participate in the regulation of seed oil accumulation. These TFs are distributed in the TF families like AP2, basic leucine zipper (bZIP), B3, NF-YB and Dof \[[@pone.0154882.ref021], [@pone.0154882.ref022], [@pone.0154882.ref024]--[@pone.0154882.ref027], [@pone.0154882.ref052], [@pone.0154882.ref053]\]. In soybean, there were 3,714 genes coding TFs \[[@pone.0154882.ref054]\], and 19 of these were found to be known lipid-related TFs ([S5 Dataset](#pone.0154882.s005){ref-type="supplementary-material"}). To identify additional potential TFs related to seed oil accumulation, we conducted a PLC network analysis between all the TFs in soybean and the above 156 genes coding core lipid synthetic enzymes. As a result, 327 TFs were found to be co-expressed with at least two genes coding core lipid synthetic enzymes (P-value\<1e-4) ([S5 Dataset](#pone.0154882.s005){ref-type="supplementary-material"}), and were classified into 42 gene families. Among the 327 TFs, 67 were distributed in high-oil dicot-specific clusters ([S5 Dataset](#pone.0154882.s005){ref-type="supplementary-material"}), and 70 had expression levels greater than twice the average level in at least one of stages from 5--6 mg to 400--500 mg. These 70 TFs were compared with the 69 seed-specific TFs described by Song et al. \[[@pone.0154882.ref025]\]. Seven genes were found to be seed-specific, respectively coding tandem CCCH zinc finger protein 4 (*GmTZF4*; *Glyma12g13300*, *Glyma06g44440*, *Glyma12g33320*), growth-regulating factor 5 (*GmGRF5*; *Glyma07g04290*), ABI5 (*Glyma10g08370*), ABI3 (*Glyma18g38490*) and a sequence-specific DNA binding TF (*Glyma03g34730*) ([Table 2](#pone.0154882.t002){ref-type="table"}). Since these seven seed-specific TFs were not only highly up-regulated during the rapid oil accumulation phase of seed development but also co-expressed with key enzymes of lipid synthesis genes, we deduced that they possibly play roles in regulating seed oil accumulation. It should be noted that ABI3 has been confirmed to be a regulator in oleosin gene expression \[[@pone.0154882.ref017], [@pone.0154882.ref028]\].

10.1371/journal.pone.0154882.t002

###### Pathway level co-expression analysis between seven seed-specific transcript factors and genes in the key pathways of lipid synthesis.

![](pone.0154882.t002){#pone.0154882.t002g}

  Transcript factor            Subfamily   P-value in differential expression analysis   Pathway-level co-expression analysis   Arabidopsis homolog of soybean transcript factor                                  
  ---------------------------- ----------- --------------------------------------------- -------------------------------------- -------------------------------------------------- --------- --------- ---------- -------------------------------
  *Glyma12g13300* (*GmTZF4*)   C3H         3.73e-34                                      3.10e-09                               *Glyma18g05840*                                    *VAL*     0.9728    5.87e-05   *At1g03790* (*AtTZF4*; *SOM*)
                                                                                                                                *Glyma04g08220*                                    *OLE*     0.9930    1.05e-06   
                                                                                                                                *Glyma06g08290*                                    *OLE*     0.9891    3.87e-06   
                                                                                                                                *Glyma09g07520*                                    *DGAT1*   0.9683    9.26e-05   
                                                                                                                                *Glyma08g01390*                                    *STERO*   0.9997    8.20e-11   
                                                                                                                                *Glyma07g32850*                                    *SAD*     0.9898    3.16e-06   
  *Glyma06g44440* (*GmTZF4*)   C3H         6.21e-31                                      1.29e-09                               *Glyma19g03530*                                    *BCCP2*   0.9757    4.22e-05   *At1g03790* (*AtTZF4*; *SOM*)
                                                                                                                                *Glyma04g08220*                                    *OLE*     0.9855    9.10e-06   
                                                                                                                                *Glyma09g07520*                                    *DGAT1*   0.9929    1.08e-06   
                                                                                                                                *Glyma11g19570*                                    *PDAT2*   0.9717    6.62e-05   
                                                                                                                                *Glyma08g01390*                                    *STERO*   0.9871    6.39e-06   
                                                                                                                                *Glyma07g32850*                                    *SAD*     0.9768    3.68e-05   
  *Glyma12g33320* (*GmTZF4*)   C3H         8.79e-37                                      1.14e-07                               *Glyma03g29600*                                    *LPAAT*   -0.9759   4.08e-05   *At1g03790* (*AtTZF4*; *SOM*)
                                                                                                                                *Glyma19g32420*                                    *LPAAT*   -0.9705   7.43e-05   
                                                                                                                                *Glyma18g38490*                                    *ABI3*    0.9904    2.66e-06   
  *Glyma07g04290* (*GmGRF5*)   GRF         2.20e-13                                      3.15e-08                               *Glyma09g07520*                                    *DGAT1*   0.9915    1.84e-06   *At3g13960* (*AtGRF5*)
                                                                                                                                *Glyma11g19570*                                    *PDAT2*   0.9912    2.06e-06   
  *Glyma03g34730*              Trihelix    1.83e-17                                      2.58e-15                               *Glyma13g06080*                                    *BCCP2*   0.9715    6.76e-05   
                                                                                                                                *Glyma19g03530*                                    *BCCP2*   0.9798    2.43e-05   
                                                                                                                                *Glyma09g07520*                                    *DGAT1*   0.9940    6.63e-07   
                                                                                                                                *Glyma11g19570*                                    *PDAT2*   0.9934    8.84e-07   
  *Glyma10g08370* (*ABI5*)     bZIP        4.66e-21                                      5.97e-09                               *Glyma03g29600*                                    *LPAAT*   -0.9737   5.30e-05   *At2g36270* (*ABI5*)
                                                                                                                                *Glyma19g32420*                                    *LPAAT*   -0.9831   1.43e-05   
                                                                                                                                *Glyma01g43780*                                    *STERO*   0.9751    4.49e-05   
                                                                                                                                *Glyma08g47240*                                    *ABI3*    0.9692    8.49e-05   
                                                                                                                                *Glyma18g38490*                                    *ABI3*    0.9927    1.18e-06   
  *Glyma18g38490* (*ABI3*)     B3          9.51e-15                                      1.26e-01                               *Glyma03g29600*                                    *LPAAT*   -0.9855   8.98e-06   *AT3G24650* (*ABI3*)
                                                                                                                                *Glyma19g32420*                                    *LPAAT*   -0.9787   2.85e-05   

Co-expression network analysis was conducted between seven seed-specific TFs and genes involved in the core lipid synthesis pathways at the 1e-04 level. Group H-L1: HanDou 5 (high-oil) and ZYD4364 (low-oil); Group H-L2: HanDou 5 and Y117249 (low-oil).

We investigated co-expression networks of the above 19 known lipid-related TFs with 156 genes coding core lipid synthetic enzymes. As a result, 14 lipid-related TFs were co-expressed with more than 10 lipid synthesis genes (P-value\<0.05) ([S5 Dataset](#pone.0154882.s005){ref-type="supplementary-material"}), such as two genes (*Glyma15g34770* and *Glyma08g24420*) coding WRI1 were found to be co-expressed with 19 lipid synthesis genes ([Table 3](#pone.0154882.t003){ref-type="table"}). Among the 19 lipid synthesis genes, genes coding BCCP2, BC, ENR and α-PDH have been verified to be regulated by *WRI1* in *Arabidopsis* \[[@pone.0154882.ref053], [@pone.0154882.ref020], [@pone.0154882.ref023]\], and 18 genes have an AW-box in their promoter ([S2 Fig](#pone.0154882.s007){ref-type="supplementary-material"}) that has been shown to be a direct target of WRI1 in *Arabidopsis* \[[@pone.0154882.ref053]\]. Note that *WRI1* has a high co-expression relationship with *Glyma06g11860* (*LACS9*; r = 0.9801, P-value = 2.32e-5), and *LACS9* is found to impact the biosynthesis of seed storage lipids in *Arabidopsis* \[[@pone.0154882.ref055]\], and is considered as the major LACS isoform involved in plastidial FA export for TAG formation \[[@pone.0154882.ref056]\]. We surmised that *WRI1* up-regulated *LACS9* and this regulation could lead to an increased export of FAs to the endoplasmic reticulum (ER), and a subsequent increase in the rate of FA synthesis and triacylglycerol synthesis.

10.1371/journal.pone.0154882.t003

###### Soybean genes co-expressed with *GmWRI1*.

![](pone.0154882.t003){#pone.0154882.t003g}

  Gene              Co-expression analysis with *Glyma08g24420*   Co-expression analysis with *Glyma15g34770*   Similar result in Arabidopsis                                                
  ----------------- --------------------------------------------- --------------------------------------------- ------------------------------- ---------- -------- ---------- ------------- --------------------------------------------------------------------------------------------------------------------------
  *Glyma13g06080*   *BCCP2*                                       4                                             0.8105                          1.53e-02   0.8491   8.24e-03   *At5g15530*   Baud et al. \[[@pone.0154882.ref020]\]; Maeo et al. \[[@pone.0154882.ref053]\]; Fukuda et al. \[[@pone.0154882.ref110]\]
  *Glyma18g50020*   *BCCP2*                                       4                                             0.8844                          3.93e-03   0.8771   4.66e-03   *At5g15530*   
  *Glyma19g03530*   *BCCP2*                                       3                                             0.7445                          3.36e-02   0.7842   2.17e-02   *At5g15530*   
  *Glyma05g36450*   *BC*                                          1                                             0.8222                          1.29e-02   0.7981   1.82e-02   *At5g35360*   Fukuda et al. \[[@pone.0154882.ref110]\]
  *Glyma08g03120*   *BC*                                          1                                             0.9482                          3.89e-04   0.9513   3.26e-04   *At5g35360*   
  *Glyma07g37050*   *LPD*                                         3                                             0.8636                          6.24e-03   0.8602   6.67e-03                 
  *Glyma15g15310*   *LPD*                                         1                                             0.8357                          1.04e-02   0.8177   1.38e-02                 
  *Glyma17g03560*   *LPD*                                         3                                             0.8222                          1.29e-02   0.8202   1.33e-02                 
  *Glyma08g45990*   *ENR*                                         2                                             0.8129                          1.48e-02   0.8077   1.59e-02   *At2g05990*   Baud et al. \[[@pone.0154882.ref023]\]
  *Glyma18g31780*   *ENR*                                         2                                             0.7291                          3.90e-02   0.7354   3.68e-02   *At2g05990*   
  *Glyma18g36130*   *FATA*                                        2                                             0.7311                          3.83e-02   0.7522   3.11e-02                 
  *Glyma07g05550*   *α-PDH*                                       2                                             0.9671                          1.03e-04   0.9646   1.28e-04   *At1g01090*   Baud et al. \[[@pone.0154882.ref023]\]
  *Glyma16g02090*   *α-PDH*                                       2                                             0.9464                          4.31e-04   0.9539   2.77e-04   *At1g01090*   
  *Glyma09g07520*   *DGAT1*                                       3                                             0.7525                          3.10e-02   0.7773   2.35e-02                 
  *Glyma06g11860*   *LACS9*                                       2                                             0.9801                          2.32e-05   0.9875   5.81e-06                 
  *Glyma18g42280*   *α-CT*                                        1                                             0.7890                          2.04e-02   0.7591   2.89e-02                 
  *Glyma18g42300*   *α-CT*                                        2                                             0.7972                          1.84e-02   0.7664   2.67e-02                 
  *Glyma11g19570*   *PDAT2*                                       0                                             0.8025                          1.71e-02   0.8144   1.45e-02                 
  *Glyma18g06500*   *MCMT*                                        2                                             0.8225                          1.29e-02   0.8069   1.61e-02                 

Two soybean genes (*Glyma08g24420* and *Glyma15g34770*) are homologous to *AtWRI1*, and all the genes co-expressed with the two soybean genes at the 0.05 level were listed in this table.

Interestingly, the expression of *WRI1* could increase seed oil content in both dicots and grasses, such as in *Arabidopsis* \[[@pone.0154882.ref020]\], oilseed rape \[[@pone.0154882.ref057]\], oil palm \[[@pone.0154882.ref058], [@pone.0154882.ref059]\] and maize \[[@pone.0154882.ref060]\]. However, distinct differences exist between dicots and grasses. These differences include low sequence similarity that results in dicot-specific and grass-specific clusters ([S1 Dataset](#pone.0154882.s001){ref-type="supplementary-material"}), different gene structures, and a low evolutionary rate in dicots (ω1 = 0.0861) as compared with that in grasses (ω0 = 0.2809) ([S3 Fig](#pone.0154882.s008){ref-type="supplementary-material"}). More importantly, *WRI1* is included in different regulatory networks in dicots and grasses \[[@pone.0154882.ref061]\]. In *Arabidopsis*, the expression of *WRI1* is up-regulated by *LEC1*, *LEC2* and *FUS3*, and *WRI1* is a direct target of *LEC2* \[[@pone.0154882.ref023]\], and possibly of *FUS3* \[[@pone.0154882.ref062]\]. In maize, however, no ortholog of *AtLEC2* was identified \[[@pone.0154882.ref060]\], and *WRI1* is able to regulate amino acid biosynthesis \[[@pone.0154882.ref063]\]. In addition, *LEC1* and *FUS3* in *Arabidopsis* have low sequence similarity compared with those in grasses.

Differential expression analysis between high- and low-oil soybean accessions {#sec010}
-----------------------------------------------------------------------------

A differential expression analysis to detect up-regulated genes in the high-oil materials may contribute to a better understanding of high-oil content mechanisms. In this study, RNA-seq differential expression analysis between high-oil cultivar Handou 5 (HD5; seed oil content: 22.3%) and two low oil wild soybeans ZYD4364 (11.9%) and Y117249 (12.5%) was conducted in the seeds at four developmental stages (15, 25, 35 and 55 days after flowering (DAF)). As a result, 8,356 DEGs between accessions HD5 and ZYD4364 (Group H-L1), and 5,551 DEGs between accessions HD5 and Y117249 (Group H-L2) were identified, and 3,997 common DEGs were observed (P\<0.01).

Among 1,123 acyl-lipid metabolism genes, 77 were differentially expressed in the above two groups (H-L1 and H-L2), and 28 DEGs encoded core lipid synthetic enzymes ([S2 Table](#pone.0154882.s002){ref-type="supplementary-material"} and [Fig 4B](#pone.0154882.g004){ref-type="fig"}). We found that 28 DEGs in core lipid synthesis pathways were up-regulated in all the three accessions at stages 25 and 35 DAF, during which seed oil is rapidly accumulated ([Fig 4B](#pone.0154882.g004){ref-type="fig"}), indicating that they play vital roles in seed oil accumulation. Interestingly, at early (15 DAF) and late (55 DAF) stages, the 28 DEGs still expressed at a relative high level in high-oil accession but at a relative low level in both low-oil accessions ([Fig 4B](#pone.0154882.g004){ref-type="fig"}). This phenomenon was also observed in seven seed-specific lipid-related TFs. Cluster analysis of expression patterns for the 1,123 acyl-lipid metabolism genes clearly distinguished the 15 and 55 DAF from the other two stages, and the high-oil accession from the low oil accessions ([Fig 4A](#pone.0154882.g004){ref-type="fig"}), indicating that divergent expression patterns of acyl-lipid metabolism genes in the early and late seed development stages played roles in determining different oil content \[[@pone.0154882.ref064]\]. The above 28 DEGs in core lipid synthesis pathways were compared with the above 40 highly up-regulated genes. As a result, 17 common genes were found; these are shown in red in [Fig 2](#pone.0154882.g002){ref-type="fig"}, and encode BCCP2, HAD, SAD, PP, DGAT1, PDAT1, STERO, OLE, and ABI3. These common genes were not only up-regulated during the rapid oil accumulation phase of seed development but also differentially expressed in high- and low-oil accessions, indicating the importance of these key enzymes in determining seed oil content. Clearly, there must be a relationship between high expression pattern and key proteins. However, there are some minor exceptions; for example, *PDAT1* and *PDAT2* in *Arabidopsis* have no effects on seed oil content and TAG synthesis, respectively, although the two genes are highly expressed during seed development \[[@pone.0154882.ref065], [@pone.0154882.ref066]\].

![Expression profiles of acyl-lipid metabolism genes in high- and low-oil soybean seeds.\
A: Hierarchical clusters of the soybean seed samples using expression levels of 1,123 genes in acyl-lipid metabolism. 15DAF (15 days after flowering), 25DAF, 35DAF and 55DAF are four stages of seed development. B: Expression patterns of 28 differentially expressed genes in core lipid synthesis pathways. Seven seed-specific TFs are also shown. HD5 (high-oil content in seed, 22.3%), ZYD4364 (low-oil content in seed, 11.9%) and Y117249 (low-oil content in seed, 12.5%) indicate the sample codes.](pone.0154882.g004){#pone.0154882.g004}

Differential expression analysis was also used to validate the candidate lipid-related genes predicted by PLC analysis. As a result, 17 of the 93 dicot-specific lipid-related genes ([Table 1](#pone.0154882.t001){ref-type="table"} and [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}), and 61 of the 327 lipid-related TFs were found to be differentially expressed in the above two groups (H-L1 and H-L2) ([S5 Table](#pone.0154882.s013){ref-type="supplementary-material"}). More importantly, all the seven seed-specific TFs were also differentially expressed in high- and low-oil soybean accessions, except *Glyma18g38490* coding ABI3, which was differentially expressed only between accessions HD5 and ZYD4364 ([Table 2](#pone.0154882.t002){ref-type="table"}).

Discussion {#sec011}
==========

Carbohydrate degradation and transport {#sec012}
--------------------------------------

The *de novo* synthesis of FAs in plants has a close connection with carbohydrate metabolism and transportation of its metabolic intermediates \[[@pone.0154882.ref006]\]. Interestingly, of the 93 dicot-specific genes putatively related to lipid metabolism, 42 (45.16%) and 27 (29.03%) genes were included in carbohydrate metabolic process and transport, respectively ([Table 1](#pone.0154882.t001){ref-type="table"} and [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}).

The 42 genes associated with carbohydrate metabolism were classified into three sub-groups. The first one had three genes (*Glyma0165s00200*, *Glyma01g28520*, *Glyma03g08860*), encoding proteins with a conserved domain PF00686 (starch binding domain; SBD) ([Table 1](#pone.0154882.t001){ref-type="table"} and [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}), which is also known as the family 20 carbohydrate-binding module (CBM20) and is found in great many starch degrading enzymes including alpha-amylase, beta-amylase, glucoamylase, and cyclodextrin glucanotransferase \[[@pone.0154882.ref067]\]. SBD could mediate the attachment between starch-active enzymes and starch granules and might disrupt the structure of the starch surface, thereby enhancing the amylolytic rate \[[@pone.0154882.ref068]--[@pone.0154882.ref070]\]. Accordingly, these three CBM20-like proteins possibly promote the hydrolysis of starch by combining with some amylolytic enzymes or other related glycosidases. Although all the three CBM20-like genes were not differentially expressed at the 0.01 probability level, *Glyma01g28520* was differentially expressed at the 0.05 level ([S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}).

The second sub-group includes the dicot-specific genes that participate in the degradation of other carbohydrates, such as homogalacturonan degradation (PWY-1081) and chitin degradation II (PWY-6902) ([Table 1](#pone.0154882.t001){ref-type="table"} and [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}). Note that *Glyma09g04330* in chitin degradation II was differentially expressed between high- and low-oil soybean accessions. These carbohydrates degradation genes might be involved in the increase in carbon sources during FA biosynthesis. The last sub-group comprised some dicot-specific genes with unknown functions ([Table 1](#pone.0154882.t001){ref-type="table"} and [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}).

Among the 27 dicot-specific genes in transport, four genes (*Glyma07g38830*, *Glyma13g27680*, *Glyma15g11270* and *Glyma18g08740*) belonged to *PF03151* (Triose-phosphate Transporter family). The first three genes were annotated as genes encoding the glucose 6-phosphate/phosphate translocator 2 (GPT2), which is located in the envelope of the plastids for the importation of Glc6P \[[@pone.0154882.ref071]\]. In non-green plastids, glucose 6-phosphate (Glc6P) can be used as a precursor for starch biosynthesis and FA synthesis and also as a substrate for the Oxidative Pentose Phosphate Pathway (OPPP) which supplies reducing power to drive FA synthesis \[[@pone.0154882.ref004], [@pone.0154882.ref072]\]. Although *AtGPT2* has no obvious effects on plant development \[[@pone.0154882.ref073], [@pone.0154882.ref074]\], its over-expression could increase the net import of Glc6P from cytosol to chloroplast and accelerate the accumulation of soluble sugars in *Arabidopsis* \[[@pone.0154882.ref075], [@pone.0154882.ref076]\]. Moreover, the expression of *GPT2* in oil palm increased notably during fruit ripening and was significantly higher than in date palm, a low oil content palm \[[@pone.0154882.ref077]\]. More importantly, *Glyma15g11270* coding GPT2 was also differentially expressed in high- and low-oil soybean accessions ([Table 1](#pone.0154882.t001){ref-type="table"} and [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}). In that case, the specific GPT2 in high-oil content dicots implies a strong funneling of carbon toward pyruvate in the plastid, a significant increase for FA synthesis and ultimately an increase in the seed oil content \[[@pone.0154882.ref077]\]. The remaining 23 transport genes had unclear functions ([Table 1](#pone.0154882.t001){ref-type="table"} and [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}).

FA synthesis and regulation of seed oil accumulation {#sec013}
----------------------------------------------------

### Functions of dicot-specific genes encoding key enzymes of FA synthesis {#sec014}

Among the 23 dicot-specific genes encoding key enzymes of FA synthesis in clusters 3 to 5 ([S2 Table](#pone.0154882.s010){ref-type="supplementary-material"}), 10 are involved in FA synthesis, i.e., *Glyma05g36450* and *Glyma08g03120* encode Biotin Carboxylase (*BC*); *Glyma13g06080*, *Glyma18g50020* and *Glyma19g03530* encode Biotin Carboxyl Carrier Protein 2 (*BCCP2*); *Glyma18g42280* and *Glyma18g42300* encode Carboxyltransferase alpha subunit (*α-CT*); *Glyma15g34770* and *Glyma08g24420* belong to the *AP2/EREBP* family and encode soybean *WRI1*; and *Glyma16g05480* encodes the FUSCA3 (FUS3). All these genes except the two genes encoding BC were highly up-regulated during the rapid oil accumulation phase of seed development. *Glyma19g03530* encoding BCCP2 was up-regulated in high-oil soybean accessions compared with the two low-oil accessions ([Fig 2](#pone.0154882.g002){ref-type="fig"}).

*BC*, *BCCP* and *α-CT* are the subunits of heteromeric acetyl-CoA carboxylase (ACCase), which is the rate limiting enzyme of FA synthesis \[[@pone.0154882.ref078]\]. In dicots, there are two forms of Acetyl-CoA Carboxylase (ACCase), a heteromeric form in the plastid and a homomeric form in the cytosol. In grasses, however, ACCase in both the plastid and the cytosol is the homomeric type \[[@pone.0154882.ref079]\]. Previous studies have shown that ACCase is able to control the rate of carbon flow in plant leaves and its expression level is related to oil content in seed \[[@pone.0154882.ref009],[@pone.0154882.ref011]\]. Specific heteromeric ACCase and its high expression level in high-oil dicots are probably key factors leading to the fact that oil content in high-oil dicots is significantly higher than that in grasses.

*WRI1* plays an important role in controlling the rate of carbon flow from carbohydrate metabolism to lipid synthesis, and is capable of affecting the seed oil accumulation by regulating a set of genes involved in lipid synthesis, glycolysis and photosynthesis \[[@pone.0154882.ref080], [@pone.0154882.ref057], [@pone.0154882.ref047]\]. Although *WRI1* could regulate seed oil accumulation in both grasses and high-oil dicots, different gene structures and divergent evolutionary rates between these two lineages were observed ([S3 Fig](#pone.0154882.s008){ref-type="supplementary-material"}). We deduced that different evolutionary mechanisms of *WRI1* in high-oil dicots and grasses might lead to different regulatory networks \[[@pone.0154882.ref061]\].

Although up-regulation of key enzymes in lipid synthesis like ACCase alone slightly increases oil content \[[@pone.0154882.ref011]\], its effect on increasing the oil content was much less than that by up-regulation of *WRI1* or other TFs \[[@pone.0154882.ref047], [@pone.0154882.ref059]\]. This partly indicates that oil accumulation is a complex biological process and increasing the expression of TFs may be an effective approach to significantly improve seed oil content.

### Functions for up-regulated genes of FA synthesis in seed {#sec015}

In the initial period of FA synthesis, apart from five dicot-specific genes coding ACCase, five genes coding PDHC were also found to be highly up-regulated during the rapid oil accumulation phase of seed development ([Fig 2](#pone.0154882.g002){ref-type="fig"}). In several species, the expression of *PDHC* has been reported to be associated with seed oil content \[[@pone.0154882.ref081]--[@pone.0154882.ref083]\]. In this study, *Glyma07g05550* encoding a subunit of PDHC was found to be a DEG between high- and low-oil soybean accessions ([Fig 2](#pone.0154882.g002){ref-type="fig"}). A similar phenomenon between high- and low-oil accessions in oat was also observed by Hayden et al. \[[@pone.0154882.ref084]\]. Therefore, we assumed that the up-regulated expression of *PDHC* and *ACCase* in soybean likely resulted in the increase of carbon flux to FA synthesis, and then to an increase of the efficiency of FA synthesis \[[@pone.0154882.ref084]\].

FA transportation, triacylglycerol synthesis and oil-body formation {#sec016}
-------------------------------------------------------------------

### FA transportation {#sec017}

Among the 54 dicot-specific genes involved in acyl-lipid metabolism, 12 genes encoded lipid transfer proteins (LTPs), including 5 DEGs between high- and low-oil soybean accessions ([S2 Table](#pone.0154882.s010){ref-type="supplementary-material"}). Apart from the above 12 high-oil dicot-specific LTPs, we also identified nine additional DEGs coding LTPs ([S2 Table](#pone.0154882.s010){ref-type="supplementary-material"}). Wang et al. \[[@pone.0154882.ref064]\] hypothesizes that an increase of the number of *LTP1* genes in sesame might enhance oil accumulation by strengthening the transport of FAs, acyl-CoAs, and other lipid molecules. On this basis, we proposed that abundant *LTP* genes in oil plants might possibly benefit oil accumulation. Among the 40 up-regulated genes during the stages from 5--6 mg to 400--500 mg ([Fig 2](#pone.0154882.g002){ref-type="fig"}), two genes *Glyma08g01180* and *Glyma06g11860* encoded ACP4 and LACS9, respectively. In addition, two genes (*Glyma13g11700* and *Glyma20g07280*) that belong to LACS family are up-regulated DEGs in high-oil accession ([Fig 4B](#pone.0154882.g004){ref-type="fig"}). Therefore, dicot-specific or up-regulation of *LTP*, *ACP4* and *LACS* in dicots could possibly increase the efficiency of plastidial fatty acid export for TAG synthesis and then consequently regulate seed oil content \[[@pone.0154882.ref055], [@pone.0154882.ref085], [@pone.0154882.ref086]\].

### Triacylglycerol synthesis and oil-body formation {#sec018}

In TAG assembly, eight genes were highly up-regulated. Among these eight, *Glyma04g04060* encoding PP, *Glyma09g07520* and *Glyma13g16560 en*coding DGAT1, and *Glyma13g16790* encoding PDAT1 were also differentially expressed in high- and low-oil soybean accessions ([Fig 2](#pone.0154882.g002){ref-type="fig"}). *DGAT1* has a principal role in TAG biosynthesis \[[@pone.0154882.ref087]\] and over-expression of *DGAT1* had been shown to enhance oil accumulation \[[@pone.0154882.ref016], [@pone.0154882.ref032], [@pone.0154882.ref088]\]. In *B*. *napus*, two domestication-related genes *BnaA01g32210D* and *BnaAnng30990D*, encoding PP and DGAT1, respectively \[[@pone.0154882.ref089]\], possibly played roles in increasinging seed oil content. In oil-body formation, three and six genes respectively encoding STEROs and OLEs were identified and were up-regulated by some hundredfold between the stages from 5--6 mg to 400--500 mg ([Fig 3](#pone.0154882.g003){ref-type="fig"} and [S4 Dataset](#pone.0154882.s004){ref-type="supplementary-material"}). Phylogenetic analysis showed more copies of genes coding OLEs in the dicots than in the grasses ([Fig 3](#pone.0154882.g003){ref-type="fig"}). Highly up-regulated and abundant OLEs possibly play vital roles in regulating seed oil content in dicots. Evidently, Siloto et al. \[[@pone.0154882.ref090]\] and Miquel et al. \[[@pone.0154882.ref051]\] showed that the size and spatial distribution of oil bodies affect the total lipid content and oil body proteins have specific functions in lipid accumulation.

We also found that 40 highly up-regulated genes and 28 DEGs in high- and low-oil soybean accessions were all enriched significantly in downstream section of the TAG biosynthesis pathway and the oil-body formation pathway (Figs [2](#pone.0154882.g002){ref-type="fig"} and [4B](#pone.0154882.g004){ref-type="fig"}). This phenomenon was also observed in sesame \[[@pone.0154882.ref064]\]. Therefore, we deduced that enzymes in the downstream section of the TAG biosynthesis pathway and the oil-body formation pathway played vital roles in the variation of seed oil content in soybean and sesame \[[@pone.0154882.ref064]\]. In other words, the efficient flow of fatty acids to formation of TAG and oil-body may ultimately influence seed oil content.

Signal transduction and other factors involved in lipid metabolism {#sec019}
------------------------------------------------------------------

Among the 54 dicot-specific genes involved in acyl-lipid metabolism ([S2 Table](#pone.0154882.s010){ref-type="supplementary-material"}), *Glyma07g01310*, *Glyma08g20710* and *Glyma15g02710* encode Phospholipase Dα (*PLDα*), which is able to hydrolyze phospholipids, producing signalling molecule phosphatidic acid \[[@pone.0154882.ref091]\]. A suppression of the expression of *PLDα* led to a significant decrease in triacylglycerol levels in *Arabidopsis* leaves \[[@pone.0154882.ref026]\] and could also slow the conversion of phosphatidylcholine to TAG in soybean seeds \[[@pone.0154882.ref092]\]. Two transgenic *B*. *napus* cultivars expressing an *Arabidopsis PLDα1* both demonstrated a 9% increase in seed total oil content \[[@pone.0154882.ref093]\]. In addition, *PLDα* is regulated by *GmMYB73*, which has been shown to be an important TF regulating lipid content \[[@pone.0154882.ref026]\]. Although these three dicot-specific *PLDα* are not DEGs in high- and low-oil soybean accessions, two of their homologous genes, *Glyma13g44170* and *Glyma08g22600*, are found to be differentially expressed genes ([S2 Table](#pone.0154882.s010){ref-type="supplementary-material"}). We suspect that up-regulated dicot-specific *PLDα* and the other two DEGs of *PLDα* would possibly have a significant effect on the seed oil content.

Among the 93 dicot-specific genes associated with acyl-lipid metabolism, 11 genes were in the category GO:0007165 (signal transduction) ([Table 1](#pone.0154882.t001){ref-type="table"} and [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}). Among these 11 genes, three (*Glyma09g07090*, *Glyma15g18380* and *Glyma17g06290*) encode GATA TFs, which play a role in light-mediated transcriptional regulation \[[@pone.0154882.ref094]\]. The three genes were homologous to the *Arabidopsis* gene *At5g56860* that is closely associated with lipid metabolism in green algae \[[@pone.0154882.ref095]\] and is capable of regulating carbon and nitrogen metabolism \[[@pone.0154882.ref096]\]. In this study, *Glyma15g18380* was validated to be a DEG in high- and low-oil soybean accessions. We deduced that these specific GATA TFs might affect the accumulation of seed oil in dicots. Two genes, *Glyma09g39570* and *Glyma10g12130*, are included in category PWY-5035 (gibberellin biosynthesis III) and might be related to the synthesis of the hormone gibberellin. Similarly, five dicot-specific genes, *Glyma04g02660*, *Glyma06g02690*, *Glyma14g40400*, *Glyma17g37750* and *Glyma17g37760* in GO:0005975 (carbohydrate metabolism) were annotated as PF02704 (Gibberellin regulated protein). The expression of genes encoding gibberellin regulated proteins is up-regulated by gibberellin \[[@pone.0154882.ref097]\]. It is known that the gibberellin signalling pathway is related to FA content \[[@pone.0154882.ref098]\]. More importantly, four genes (*Glyma04g02660*, *Glyma06g02690*, *Glyma17g37750* and *Glyma17g37760*) encoding gibberellin regulated proteins were also validated to be DEGs in high- and low-oil soybean accessions. Therefore, the above six gibberellin related genes may play roles in FA synthesis and thereby affect the oil content.

### Other factors involved in lipid metabolism {#sec020}

Apart from the above specific genes directly participating in acyl-lipid metabolism, we also found that five dicot-specific genes in soybean were putatively involved in the lipid metabolic process (GO:0006629) ([Table 1](#pone.0154882.t001){ref-type="table"} and [S4 Table](#pone.0154882.s012){ref-type="supplementary-material"}). Among these genes, *Glyma02g45680* and *Glyma14g03130* encoded cytochrome P450, which plays a role in the FA metabolism in plants \[[@pone.0154882.ref099], [@pone.0154882.ref100]\]); and *Glyma18g45250* encoded a protein participating in glyceollin biosynthesis, which was related to lipid peroxidation \[[@pone.0154882.ref101]\]. Since these genes were highly expressed and coordinated with genes participating in acyl lipid synthesis, we hypothesized that these genes and their associated pathways possibly play some roles in seed oil accumulation.

Conclusion {#sec021}
----------

Ninety-three dicot-specific genes, including 42 and 27 genes respectively in carbohydrate degradation and transport, were predicted to be candidate genes associated with acyl-lipid metabolism pathways. And seed-specific TF genes *GmGRF5*, *ABI5* and *GmTZF4* were also predicted to play roles in regulating seed oil accumulation. Furthermore, *ACCase*, *DGAT1*, *PP*, *OLEs* and *STEROs* were highly up-regulated not only in specific stages of seed development but also in high-oil accessions, which indicates that enzymes in initial fatty acid synthesis, downstream of TAG biosynthesis and oil-body formation, played vital roles in the variation of seed oil content. In particular, highly up-regulated and abundant *OLEs* possibly play vital roles in determining seed oil content. Most of the above key genes were further confirmed by differential gene expression analysis between high-oil cultivated and low-oil wild soybeans.

Materials and Methods {#sec022}
=====================

Genomic data and gene expression data {#sec023}
-------------------------------------

The genomic data of four high-oil dicot species (*Glycine max* (seed oil content: 20%), *Gossypium raimondii* (30%), *Ricinus communis* (50%), and *Arabidopsis thaliana* (35%)) and three low-oil grasses (*Sorghum bicolor* (3%), *Setaria italic* (1.7%), and *Oryza sativa* (3%)) were downloaded from Phytozome V9.1 (<http://www.phytozome.net/>) \[[@pone.0154882.ref102]\]. The longest encoded protein sequence was chosen for genes with multiple transcripts.

The transcriptome data with two biological replicates of *G*. *max* Williams 82 \[[@pone.0154882.ref041]\] were downloaded from the Gene Expression Omnibus database (<http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42871>). The data included seven stages of seed development: whole seed 4 DAF, whole seed 12--14 DAF, whole seed 22--24 DAF, whole seed 5--6 mg in weight, cotyledons 100--200 mg in weight, cotyledon 400--500 mg in weight, and dry whole seed.

GO annotation and GO enrichment analysis {#sec024}
----------------------------------------

The GO annotations of the *G*. *max* genes, including molecular function, molecular location and biological process, were conducted using the online tool Goanna (<http://agbase.msstate.edu/cgi-bin/tools/GOanna.cgi>) \[[@pone.0154882.ref103]\]. The GO enrichment analysis was performed using GOstats with a threshold P value of less than 0.01 \[[@pone.0154882.ref104]\]. The GO slims, which are a subset of GO terms for a high level summary of the ontology content, were summarized using GOSlimViewer (<http://www.agbase.msstate.edu/cgi-bin/tools/goslimviewer_select.pl>) \[[@pone.0154882.ref103]\].

OrthoMCL analysis and definition of lineage-specific clusters {#sec025}
-------------------------------------------------------------

Orthologous gene clusters were calculated from OrthoMCL comparisons of four dicots and three Grasses \[[@pone.0154882.ref040]\].

Based on all-against-all BLASTP comparisons of a set of protein sequences from genomes of interest, clusters of proteins were grouped according to reciprocal best similarity pairs between and within species, using OrthoMCL software implemented by the Markov Clustering algorithm (MCL; <http://micans.org/mcl/>) \[[@pone.0154882.ref105]\]. Here all the critical values were set as default values in the software. One OG was defined as dicot-specific if all the genes in the OG were coming from all the four rosid species but not from any grass species.

Prediction of genes involved in lipid metabolism of *Glycine max* {#sec026}
-----------------------------------------------------------------

More than 600 genes involved in acyl-lipid metabolism in Arabidopsis \[[@pone.0154882.ref007],[@pone.0154882.ref008]\] were downloaded from the website ARALIP (<http://aralip.plantbiology.msu.edu/>). Based on the OrthoMCL result, all the orthologous and paralogous genes participating in acyl lipid biosynthesis in soybean were identified.

Clustering by expression pattern and PLC network analysis {#sec027}
---------------------------------------------------------

Using the transcriptome data, gene models were clustered using BioLayout Express^3D^, implemented by the MCL \[[@pone.0154882.ref042]\]; and any two genes with an absolute value of the Pearson correlation coefficient greater than 0.7 were considered to be similar in their expression level \[[@pone.0154882.ref043], [@pone.0154882.ref106]\].

The PLC network analysis developed by Wei et al. \[[@pone.0154882.ref043]\] was used to identify new candidate pathway members in the lipid synthesis pathway. Genes co-expressed with more than two genes in the soybean acyl-lipid metabolism pathway were considered as candidate lipid synthesis pathway members. A relatively stringent correlation threshold was set at the probability value of 1e-4.

RNA-seq for transcriptome analysis of high- and low-oil soybean accessions {#sec028}
--------------------------------------------------------------------------

The materials used for RNA-seq to analyze lipid synthesis were three soybean accessions, one high-oil cultivar Handou 5 (HD5; seed oil content: 22.3%) and two low oil wild soybeans ZYD4364 (11.9%) and Y117249 (12.5%). Whole seeds at stages 15, 25, 35 and 55 DAF were harvested as samples. Total RNA of every sample was extracted from tissues using the TRIzol reagent (Invitrogen). The quality and the concentration of total RNA were quantified separately using Agilent 2200 TapeStation and ND-1000 Nanodrop. cDNA library were constructed using the same procedure described in Severin et al. \[[@pone.0154882.ref107]\] and sequenced using an Illumina HiSeq 2500 sequencing platform.

The raw reads were cleaned by removing reads with adapters and those of low quality. Clean reads were mapped to reference sequences using SOAPaligner/soap2 (<http://soap.genomics.org.cn/soapdenovo.html>). Mismatches no more than two bases were allowed in the alignment. The gene expression level was calculated by using RPKM method (Reads Per kb per Million reads) \[[@pone.0154882.ref108]\]. Fisher's exact test method in DEGseq package \[[@pone.0154882.ref109]\] was used to identify DEGs between high- and low-oil accessions, at the 0.01 significant level.

Supporting Information {#sec029}
======================

###### OGs of all the genes and genes in acyl-lipid metabolism of seven species.

Dicot-specific OGs shared with all the four dicots and 4051 dicots-specific genes in soybean are also shown.

(XLSX)

###### 

Click here for additional data file.

###### GO enrichment analysis in *Glycine max* for 4051 soybean genes in dicot-specific OGs.

GOslim summary of the enriched GO terms are also listed.

(XLSX)

###### 

Click here for additional data file.

###### GO enrichment analysis for 805 dicot-specific genes of clusters 3 to 5 in *Glycine max*.

GOslim summary of the enriched GO terms are also listed.

(XLSX)

###### 

Click here for additional data file.

###### Up-regulated genes of lipid synthesis core pathways in at least one of the seed developmental stages from 5--6 mg to 400--500 mg.

(XLSX)

###### 

Click here for additional data file.

###### Predicted dicot-specific genes and TFs that are associated with acyl-lipid metabolism and their co-expression network with acyl-lipid metabolism genes in soybean.

Co-expression network of 19 known lipid-related transcription factors with core enzymes in lipid synthesis (P\<0.05) are also shown.

(XLSX)

###### 

Click here for additional data file.

###### Co-expression network among soybean genes coding *ABI3*, *OLE* and *STERO*.

Two genes with a coordinated relationship were linked by regular (P-value \< 0.01) or bold (P-value \< 1e-04) lines.

(PDF)

###### 

Click here for additional data file.

###### AW-boxes, with arrowhead, in 5'-upstream sequences of *GmWRI1* coordinated genes involved in lipid synthesis pathways.

Schematic drawing of the sequence 2 kb upstream of the ATG start codon of *GmWRI1* co-expressed genes.

(PDF)

###### 

Click here for additional data file.

###### Phylogenetic analysis of *WRI1*.

A: Phylogenetic tree of *WRI1* constructed by MEGA 6.0 using Neighbor-Joining method and the bootstrap test was performed with 1,000 iterations. Square boxes indicate duplication events and numbers on the branches represent the bootstrap support. Genes structure visualizing positions of exons and introns are also shown; this was constructed by GSDS 2.0. B: Selection detection using branch model implemented by PAML.

(PDF)

###### 

Click here for additional data file.

###### Genomic information of the seven species used in this study.

(PDF)

###### 

Click here for additional data file.

###### Acyl-lipid metabolism genes that are high-oil dicot-specific or differentially expressed genes in high- and low-oil soybean accessions.

(PDF)

###### 

Click here for additional data file.

###### 207 genes in interested biological processes that were included in the 805 genes in clusters 3 to 5.

(PDF)

###### 

Click here for additional data file.

###### Annotation and differential expression analysis of 93 dicot-specific genes that were associated with oil accumulation.

(PDF)

###### 

Click here for additional data file.

###### 327 lipid-related TFs predicted by PLC analysis.

(PDF)

###### 

Click here for additional data file.

We thank Dr Guo Liang (Huazhong Agricultural University) for helpful comments on the revised version of this manuscript.

[^1]: **Competing Interests:**The authors have declared that no competing interests exist.

[^2]: Conceived and designed the experiments: YMZ. Performed the experiments: LZ SBW QGL JS YQH LZ YMZ. Analyzed the data: LZ JS YQH. Wrote the paper: LZ YMZ HQZ JMD.
