Email updates

Keep up to date with the latest news and content from BMC Plant Biology and BioMed Central.

Open Access Research article

Sequence mining and transcript profiling to explore differentially expressed genes associated with lipid biosynthesis during soybean seed development

Huan Chen12, Fa-Wei Wang1, Yuan-Yuan Dong1, Nan Wang1, Ye-Peng Sun1, Xiao-Yan Li1, Liang Liu2, Xiu-Duo Fan2, Hai-Long Yin1, Yuan-Yuan Jing1, Xin-Yue Zhang1, Yu-Lin Li1, Guang Chen2* and Hai-Yan Li12*

Author Affiliations

1 Ministry of Education Engineering Research Center of Bioreactor and Pharmaceutical Development, Jilin Agricultural University, Changchun, Jilin, 130118, China

2 College of Life Sciences, Jilin Agricultural University, Changchun, Jilin, 130118, China

For all author emails, please log on.

BMC Plant Biology 2012, 12:122  doi:10.1186/1471-2229-12-122

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2229/12/122


Received:11 December 2011
Accepted:27 June 2012
Published:31 July 2012

© 2012 Chen et al.; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Soybean (Glycine max L.) is one of the most important oil crops in the world. It is desirable to increase oil yields from soybean, and so this has been a major goal of oilseed engineering. However, it is still uncertain how many genes and which genes are involved in lipid biosynthesis.

Results

Here, we evaluated changes in gene expression over the course of seed development using Illumina (formerly Solexa) RNA-sequencing. Tissues at 15 days after flowering (DAF) served as the control, and a total of 11592, 16594, and 16255 differentially expressed unigenes were identified at 35, 55, and 65 DAF, respectively. Gene Ontology analyses detected 113 co-expressed unigenes associated with lipid biosynthesis. Of these, 15 showed significant changes in expression levels (log2fold values ≥ 1) during seed development. Pathway analysis revealed 24 co-expressed transcripts involved in lipid biosynthesis and fatty acid biosynthesis pathways. We selected 12 differentially expressed genes and analyzed their expressions using qRT-PCR. The results were consistent with those obtained from Solexa sequencing.

Conclusion

These results provide a comprehensive molecular biology background for research on soybean seed development, particularly with respect to the process of oil accumulation. All of the genes identified in our research have significance for breeding soybeans with increased oil contents.

Keywords:
Gene expression; Lipid; RNA-sequencing; Soybean; Unigene

Background

Soybean [Glycine max (L.) Merrill], one of the most important oil crops in the world that accounts for less than a third (27.7%) of the total global vegetable oil production. At present, the main producers of soybean are the United States, Brazil, Argentina, India, and China. Although world soybean oil production has steadily increased, rising from 37.73 million metric tons in 2007 to 43.65 million metric tons in 2012 (http://www.fas.usda.gov/psdonline/psdreport.aspx?hidReportRetrievalName=BVS&hidReportRetrievalID=708&hidReportRetrievalTemplateID=8 webcite), the production of soybean oil is still insufficient to meet the consumer demand, and the low oil content in soybean seeds is the main factor restricting yield. Therefore, an increase in soybean oil yield is desirable and has been one of the most major goals of oilseed engineering.

To increase the production of soybean oil, it is not feasible to extend the planting area of the crop, because of increased competition for land by the rapidly rising population. Therefore, it is a better strategy to increase the seed oil content than to increase the planting area [1]. In recent years, molecular genetics approaches have been used to modify seed oil content. Over-expression of a diacylglycerol acyltransferase (AtDGAT) cDNA in wild-type A. thaliana enhanced oil deposition and average seed weight [2]. The research of Wang et al. indicated that the oil content of soybean seeds could be increased by upregulation of two soybean Dof-type transcription factor (GmDof) genes, which are associated with fatty acid biosynthesis [3]. The soybean genome also has been examined with regard to the relationships between gene expression profiles and gene function. Although the biochemical pathways that produce different storage components are well characterized, there is still no integrated model describing differentially expressed genes (DEGs) involved in soybean lipid biosynthesis. With the development of high-throughput technologies, including the newly developed Solexa/Illumina RNA-sequencing and DEG high-throughput deep sequencing approaches, new genes have been discovered and specific transcripts analyzed. In Severin’s research, RNA Seq-Atlas provided a record of high-resolution gene expression in a set of various tissues. They also found dramatic tissue-specific gene expression of both the most highly-expressed genes and the genes specific to legumes in seed development and nodule tissues [4]. Furthermore, these technologies are useful for estimating overall gene expressions at different development stages and/or in different tissues, such as on rice transcriptome and disease [5-7].

In this study, we used Illumina Solexa sequencing to investigate gene expression in soybean seeds at different developmental stages (15, 35, 55, and 65 days after flowering; DAF), and then compared transcript reads with the most recent release of the G. max genome sequence (assembly Glyma1.01). Illumina Solexa de novo sequencing technology identified a total of 11592, 16594, and 16255 differentially expressed unigenes between the 15 and 35 DAF seeds, the 15 and 55 DAF seeds, and the 15 and 65 DAF seeds, respectively (Additional file 1: Tables S1, Additional file 2: Table S2, Additional file 3: Table: S3). Of these, 9905 unigenes existed in all three of these contrast groups and represented co-expressed unigenes. Furthermore, 124 candidate unigenes were screened from the three contrasting cDNA libraries and characterized as those responsible for lipid synthesis. Our researches provide the whole picture of DEG transcription patterns and expression levels during soybean seed development. The elucidation of DEG transcription patterns at specific stages of seed development also lays the foundation for understanding the molecular mechanisms underlying oil production. This research provides direction for controlling genes over-expression to increase soybean oil content.

Additional file 1. Table S1. Differentially expressed genes between 35 DAF and 15 DAF (FDR ≤0.001 and |log2ratio| ≥1).

Format: XLS Size: 2.4MB Download file

This file can be viewed with: Microsoft Excel ViewerOpen Data

Additional file 2. Table S2. Differentially expressed genes between 55 DAF and 15 DAF (FDR ≤0.001 and |log2ratio| ≥1).

Format: XLS Size: 3.4MB Download file

This file can be viewed with: Microsoft Excel ViewerOpen Data

Additional file 3. Table S3. Differentially expressed genes between 65 DAF and 15 DAF (FDR ≤0.001 and |log2ratio| ≥1).

Format: XLS Size: 3.4MB Download file

This file can be viewed with: Microsoft Excel ViewerOpen Data

Methods

Plant culture and collection

The soybean cultivar Jiyu-72 was used as the experimental material. The plants were grown in a green house with a 15-h light (200 μEm-2 s-1, 25°C)/9-h dark photoperiod (23°C), with relative humidity controlled at 75%. The developmental processes of soybean seeds from flowering to seed maturity were observed from July to October 2010. Pods were harvested at 15 DAF (immature stage), and then every 5 days until 70 DAF (pods containing mature seeds). After removing the seed coat, the seeds were used for oil extraction or frozen at −80°C until mRNA extraction and sequencing.

Measurement of oil content

To extract the oil (or lipid), seeds harvested at 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, and 70 DAF were oven-dried to constant weight at 85°C. The dry samples were ground to a fine powder in a disintegrator, and the powder was transferred into 10-mL glass tubes for oil extraction. Oil was extracted using ligarine to determine total lipids (TL) gravimetrically [8]. Oil extraction was performed using a SER148 3/6 Extraction apparatus (VELP Scientifica, Italy). Extractions were carried out using triplicate samples for each stage, and mean values were determined. Errors are shown as standard deviations.

Total RNA extraction, library construction, and deep sequencing

RNA was isolated from the seeds harvested at 15, 35, 55, and 65 DAF, respectively, using a TIANGEN RNA Prep Pure Plant kit (Tiangen Biotech Co. Ltd., Beijing, China) following the manufacturer’s protocol, followed by a chloroform extraction to improve RNA purity. The yield and quality of total RNA samples were determined using a ND-1000 NanoDrop spectrophotometer. For RNA library construction and deep sequencing, equal quantities of RNA samples from seeds at the four developmental stages were pooled. Approximately 6 μg RNA representing each group was submitted to Solexa (now Illumina Inc.) for sequencing.

Differential expression (DE) detection

The raw data were filtered to remove adaptor reads, low quality reads, and reads of copy number = 1, yielding a dataset consisting of clean reads. For the tissue-specific analyses, raw digital gene expression data were normalized using a variation of the reads/Kb/million (RPKM) method [9,10]. The RPKM method corrects for biases in total gene exon size and normalizes for the total short read sequences obtained in each tissue library. Genes with data P-value <0.005, false discovery rate (FDR) ≤0.001, and estimated absolute log2-fold change ≥1 in sequence counts across libraries were considered to be significantly differentially expressed. Differentially expressed genes that were co-expressed during these four stages were then subjected to cluster analysis using the R program (a language and environment for statistical computing and graphics). Functional categorization of stress-regulated genes was performed using BGI WEGO (http://wego.genomics.org.cn/cgi-bin/wego/index.pl webcite).

Quantitative real-time PCR (qRT-PCR) analysis

To verify the data obtained by Solexa RNA-seq, qRT-PCR was performed on 12 genes with log2ratios ranging from 2 to 11 (Additional file 4: Table S4). The RNA samples used for the qRT-PCR assays were the same as those used for the DEG experiments. qRT-PCR was performed on a Mx3000P instrument (Stratagene), with SYBR-Green detection (Quant qRT-PCR Kit, Tiangen Biotechnology), according to the manufacturer's instructions. The soybean tubulin gene was used as an internal control (forward primer: 5’-ATGTTCCGTGGCAAGATGAG-3’, reverse primer: 5’-CATTGTTGGGAATCCACTC-3’) [11]. Primers were designed using the Primer Premier 5.0 and Oligo 6 programs to amplify products approximately 150-bp long. Thermal cycle conditions were as follows: 2 min at 95°C followed by 40 cycles of 15 s at 95°C, 15 s at 56-57°C, and 15 s at 72°C. Each cDNA was analyzed three times, after which the average threshold cycle (Ct) was calculated per sample. Standard curves were established for all investigated genes using a series of amplicon dilutions. The relative expression levels were calculated as 2- (ΔCt of 35, 55, 65 DAF -ΔCt of 15 DAF).

Additional file 4. Table S4. Details of primers used for qRT-PCR.

Format: DOC Size: 60KB Download file

This file can be viewed with: Microsoft Word ViewerOpen Data

Results

Accumulation of oil at different stages of seed development

To explore lipid accumulation during development of soybean seeds, we quantified lipid content in soybean seeds harvested from 15 to 70 DAF. As shown in Figure 1, the oil content began to increase from 15 DAF, then markedly increased from 15 to 20 DAF, and then remained stable from 20 to 25 DAF before again increasing rapidly from 25 to 35 DAF. From 35 to 55 DAF, the oil content showed a steady increase, peaking at 60 DAF before dropping at 65 DAF. These results indicated that for the soybean cultivar JiYu-72, 15 DAF marked the beginning of oil accumulation, 35 DAF was the key stage for the rapid increase in the oil content. 55 DAF was the stage at which oil content was high and stable, and 65 DAF was the stage at which oil content decreased and nearly unchanged. Therefore, in order to explore the differentially expressed genes associated with lipid biosynthesis during soybean seed development, the following sequencing and qRT-PCR analyses were performed using samples from these four different developmental stages.

thumbnailFigure 1. Changes in lipid content during seed development. Lipid content was determined every 5 days. One star stands for the discrepancy is significant and double stars stand for the discrepancy is extremely significant.

Screening of differentially expressed genes (DEGs) from massive datasets

Initially, Solexa RNA-sequencing revealed a total of 11592, 16594, 16255 differentially expressed unigenes between 35 DAF and 15 DAF seeds, 55 DAF and 15 DAF seeds, and 65 DAF and 15 DAF seeds, respectively (Table 1). Among them, 9905 unigenes were present in all three of these contrast groups, and 1687, 6689, and 6350 were differentially expressed only in the 35 vs 15 DAF, 55 vs 15 DAF, and 65 vs 15 DAF groups, respectively. From the large datasets, we found that down-regulated genes were more abundant than up-regulated genes during soybean seed development. The normalized data are shown as scatter plots of log10-transformed transcription counts for the four different samples (Additional file 5: Figure S1).

Additional file 5. Figure S1. Scatter plot of differentially expressed genes in seeds harvested at 35 DAF, 55 DAF, and 65 DAF, compared with 15 DAF. (a). Scatter plot of differentially expressed genes between 35 DAF and 15 DAF. (b). Scatter plot of differentially expressed genes between 55 DAF and 15 DAF. (c). Scatter plot of differentially expressed genes between 65 DAF and 15 DAF. TPM = Transcript per million (normalized expression level of genes).

Format: TIFF Size: 17KB Download fileOpen Data

Table 1. Differentially expressed genes between seeds of 35 DAF and 15 DAF, 55 DAF and 15 DAF, and 65 DAF and 15 DAF

There were 680, 1004, and 850 up-regulated unigenes and 10912, 15590, and 15405 down-regulated unigenes in the three respective groups (Table 1). Among the up-regulated unigenes, 352 were co-expressed genes showing increases of at least 2-fold (log2ratio ≥ 1), and 269, 258, and 65 of the up-regulated genes were specific to each of the three respective groups. Among the down-regulated unigenes, 9476 were co-expressed genes showing decreases of at least 2-fold (log2ratio ≤ 1), and 672, 1278, and 769 of the down-regulated genes were specific to each of the three respective groups (Figure 2).

thumbnailFigure 2. Venn diagram showing overlaps among up-regulated genes during soybean seed development. (a). All the unigenes that showed increases of at least 2-fold (log2ratio ≥ 1) in the 3 contrasts, (b). All the unigenes that showed decreases of at least 2-fold (log2ratio ≤ 1) in the 3 contrasts. Number in one circle denotes stage-specific genes, and number in two or more intersecting circles denotes overlapped genes.

Differentially expressed genes provide clues about the molecular events related to seed development. Further investigation of highly expressed genes may be warranted to determine what functional roles the highly expressed sequences may play in soybean seeds. To understand the relationship between the time at which unigenes are co-expressed and their biological significance, we carried out Gene Ontology enrichment analysis.

Gene ontology category analysis and DEGs in each category

As a useful tool for gene functional annotation, WEGO (Web Gene Ontology Annotation Plot) has been widely used in many important studies, including the rice genome project and the silkworm genome project [12-14]. It has become one of the most commonly used tools for downstream gene annotation analysis, especially in comparative genomics studies. In this research, WEGO analysis assigned the DEGs (FDR ≤ 0.001 and |log2ratio| ≥ 1) to three functional categories; Cellular Component, Molecular Function, or Biological Process. A total of 11963 unigenes had at least one GO functional category. In the three contrast groups, all differentially expressed genes with a Gene Ontology annotation were further classified into subsets. There were 11 subsets within the Cellular Component category, 10 subsets within the Molecular Function category, and 23 subsets within the Biological Process category (Figure 3). Thus, the most abundant unigenes were related to binding and catalytic activity in the Molecular Function category, cell, cell parts and organelles in the Cellular Component category, and cellular and metabolic in the Biological Process category. The results showed that these types of genes were highly enriched in our soybean transcriptomes. In terms of the categories, the three contrast groups showed similar patterns.

thumbnailFigure 3. Functional categorization of significantly differentially expressed genes during seed development. Functional categorization was performed using BGI WEGO (Web Gene Ontology Annotation Plot). (a). The functional categorization of differentially expressed genes between 35 DAF and 15 DAF. (b). The functional categorization of differentially expressed genes between 55 DAF and 15 DAF. (c). The functional categorization of differentially expressed genes between 65 DAF and 15 DAF.

Because we are interested in lipid biosynthesis, we focused on genes associated with this process that were differentially expressed during seed development. From the Component Ontology data, eight differentially co-expressed lipid particle unigenes were identified. Among them, seven unigenes were up-regulated, and one unigene was down-regulated during seed development. From the Function Ontology data, unigenes associated with lipid biosynthesis could be further classified into three categories: those with lipid kinase activity (9 differentially co-expressed down-regulated unigenes); those with fatty acid ligase activity (6 differentially co-expressed down-regulated unigenes); and those with fatty acid synthase activity (18 differentially co-expressed unigenes, 4 of which were up-regulated and 14 down-regulated). From the Process Ontology data, unigenes associated with lipid biosynthesis were in two categories; lipid biosynthetic process (78 differentially co-expressed unigenes, 5 of which were up-regulated, and 73 down-regulated), and fatty acid biosynthetic process (49 differentially co-expressed unigenes, 4 of which were up-regulated, and 45 down-regulated). Interesting, the 49 co-expressed unigenes detected in the fatty acid biosynthetic process were all included in unigenes detected in the lipid biosynthetic process. Among the 113 unigenes analyzed by Gene Ontology, 15 co-expressed unigenes that showed high levels of expression (log2fold values ≥ 1) could be responsible for lipid biosynthesis (Table 2).

Table 2. DEGs identity, expression ratios, and Gene Ontology annotations

Annotation of lipid-relevant genes in fatty acid pathway

Genes usually interact with each other to carry out certain biological functions. Knowledge of the expressions of multiple genes and their regulation during oil biosynthesis is required to further understand the regulatory mechanisms controlling oil metabolism. Pathway-based analysis helps to clarify the biological functions of genes, and identifies significantly enriched metabolic pathways or signal transduction pathways associated with the DEGs compared with the whole genome background. For our research, 124 biological pathways, including the starch and sucrose metabolism pathway, the citrate cycle pathway, the fatty acid biosynthesis pathway, and many others were identified by KEGG pathway analysis of unigenes. A total of 4836, 6796, and 6758 DEGs with pathway annotations were identified in the three respective contrast groups. From those pathways, we selected the fatty acid biosynthesis pathway (Figure 4) for deep analysis. In the three contrast groups, 24 DEGs associated with the fatty acid pathway were co-expressed (Table 3). Interestingly, 13 of the identified unigenes were the same as those identified from the GO analysis. Among those 24 DGEs, only four ungenes (BT098182, AK245394, BT093926, BT093412) were up-regulated, which are encoding 3-oxoacyl-[acyl-carrier protein] reductases (FabG) and were annotated as having fatty acid synthase activity. The other 20 unigenes were down-regulated, which may be the negatively controlled genes in the fatty acid biosynthesis pathway. Figure 4 shows the locations of the differentially co-expressed unigenes in the fatty acid pathway. Those shown in red and green represent genes that were up-regulated and down-regulated, respectively, during seed development.

thumbnailFigure 4. Fatty acid pathway in soybean. Processes in pink rectangles are shown in abbreviated form. Red or green rectangles indicate up-regulated or down-regulated genes, respectively. FASN, Fas: fatty acid synthase; FAS1, 2: fatty acid synthase subunit beta, alpha; FabA: 3-hydroxydecanoyl-[acyl-carrier-protein] dehydratase; FabB, FabF, FabH: 3-oxoacyl-[acyl-carrier-protein] synthase I, II, III; FabD: [acyl-carrier-protein] S-malonyltransferase; FabG: 3-oxoacyl-[acyl-carrier protein] reductase; FabI, FabK, FabL: enoyl-[acyl-carrier protein] reductase I, II, III; FabZ: 3R-hydroxymyristoyl ACP dehydrase; MdcH: malonate decarboxylase epsilon subunit; FatA, FatB: fatty acyl-ACP thioesterase A, B.

Table 3. DEGs involved in fatty acid pathway

Clustering of candidate genes

A total of 124 unigenes, including those detected by GO and pathway analyses, were selected for clustering analysis. All of these unigenes showed different expression profiles. The genes were clustered according to the similarity of their expression patterns. A subset of the data is shown in Figure 5. The heatmap of the RPKM normalized log2-transformed transcription count was generated using R programs. Up-regulated unigenes are shown in red, while down regulated ones are shown in green. This analysis allowed us to define common qualitative patterns in gene expression changes over time. Our results suggested that only a few genes were expressed differentially in the three contrast groups, and that the majority of the transcriptome had approximately similar expression levels. In this manner, similar patterns were clustered together in the same manner as a taxonomic tree. Therefore, we can speculate on the roles of unigenes of unknown function by comparison to the annotated genes with known functions in the same cluster.

thumbnailFigure 5. Expression changes and cluster analysis of genes (124 in total) that were differentially expressed between 35 DAF and 15 DAF seeds (control), 55 DAF and 15 DAF seeds, and 65 DAF and 15 DAF seeds. Cluster analysis of genes was performed using R program. Rows represent differentially expressed genes, columns represent different contrast groups. Red, green, and black boxes represent genes showing increased, decreased, or equal expression levels, respectively. Color saturation reflects magnitude of log2 expression ratio for each gene.

Confirmation of differentially expressed genes by qRT-PCR

To confirm the expression patterns determined by Solexa RNA-sequencing analysis, we used qRT-PCR analyses to analyze expressions of 12 candidate genes (Additional file 4: Table S4). The soybean tubulin gene was used as an internal control. Although the Solexa log2-fold values showed slight variations compared with the corresponding values from the qRT-PCR analyses, the expression data from Solexa RNA-seq analysis was largely consistent with those obtained by qRT-PCR (Additional file 6: Figure S2).

Additional file 6. Figure S2. Confirmation of differential gene expression by qRT-PCR.

Format: TIFF Size: 3.8MB Download fileOpen Data

Discussion

Lipid synthesis depends on the correct spatial and temporal activity of many gene products. These genes execute their function in three stages: fatty acid synthesis in the plastid, triacylglycerol (TAG) synthesis in the endoplasmic reticulum (ER), and assembly into an oil body. Significant improvements in oil accumulation must be accompanied by changes in activity of the genes involved in fatty acid biosynthesis in developing seeds. Identification of these genes and their regulatory pathway would provide not only new genetic information for understanding soybean seed development, but also for controlling gene expression in developing seeds to alter oil accumulation.

This study has provided a new data set identifying the expression of DEGs during oil accumulation in developing seeds. Massive parallel sequencing identified 11592, 16594 and 16255 differentially expressed unigenes from three contrasting libraries covering four stages of seed development. We examined gene expression levels in detail, and found significant differences among the four growth stages. Among the differentially expressed genes, more were down-regulated than up-regulated, and some showed a differential expression pattern in all three contrast groups, indicating that there were overlaps at the transcriptional level. The fact that there were many down-regulated genes indicates that there are more negatively regulated genes than positively regulated ones with functions in the fatty acid pathway. However, it does not mean that lower expression of these genes leads to lower oil contents, because of the complexity of the lipid biosynthetic pathway. To identify the genes associated with lipid biosynthesis, we determined gene expression patterns at specific stages of seed development, and conducted a GO analysis. To explore the genes with unknown functions, the expression patterns of 124 unigenes were analyzed by hierarchical clustering according to similarities in expression profiles across all conditions.

In the lipid biosynthetic pathway, acetyl CoA carboxylase (ACCase) plays an important regulatory role. ACCase catalyzes the carboxylation of acetyl-CoA to malonyl-CoA [15], which is then transacylated by malonyl-CoA-acyl carrier protein transacylase (FabD or MCAT,[EC 2.3.1.39]) to the acyl carrier protein (ACP), forming malonyl ACP. The latter adds a 2-carbon acetyl unit to the nascent or growing fatty acyl chain. The basic reaction of fatty acid synthesis is to combine molecules composed of two carbon units into longer chains to form fatty acids. The initial condensation of 2-carbon units is catalyzed by β-ketoacyl-ACP synthase III (KASIII, EC2.3.1.180) which uses acetyl-CoA and malonyl-ACP as substrates. After the two reductions and dehydration reactions, a 4-carbon fatty acid, butyrate, is produced. This is poorly condensed by KASIII but is a good substrate for KASI, which elongates 4-carbon chains to 14-carbon chains.

In this study, we studied in detail the genes with key roles in lipid biosynthesis. Very frequently, the same enzymatic function is redundantly encoded by several unigenes. This may be the result of different proteins referenced with the same EC number or it may represent different transcripts encoding specific enzyme subunits. This situation was significant in the present study. Most plants have two forms of ACCase, the homomeric form in the cytosol, composed of a single large polypeptide catalyzing the individual carboxylation reactions [16], and the heteromeric form in plastids, composed of four subunits; biotin carboxylase (BC) [17], biotin carboxyl carrier protein (BCCP) [18], α-carboxyltransferase (α-CT) [19] and β-carboxyltransferase (β-CT) [20]. In the present study, seven co-expressed unigenes (GenBank accession nos: BT094733, AK245399, U40979, AF165159, CX703216, BM188175, AW830581) encoding ACCases were identified as DEGs in the three contrast groups. All of them were down-regulated (log2ratio ≤ 1) in the three contrast groups except one gene (AK245356) encoding an ACCase was up-regulated (log2ratio = 1.4) at 35 DAF. Three co-expressed unigenes (GenBank accession nos: BM188175, CX703216, AW830581, EC6.3.4.14) that are associated with the transformation of BCCP to carboxybiotin-carboxyl-carrier protein were down-regulated during seed development.

The next enzyme in fatty acid pathway is FabD, which catalyzes the transfer of malonyl-CoA to the holo acyl carrier protein (ACP), generating malonyl-ACP [21]. In the present study, the co-expressed unigene (AW34878) encoding FabD was down-regulated during seed development.

KASIII (FabH) catalyzes the subsequent condensation and transacylation of acetyl-CoA with malonyl-ACP and has a universal role in fatty acid biosynthesis. Transgenic B. napus seeds overexpressing KASIII driven by napin also contained lower oil levels compared to what was found in the wild-type [22]. In the present study, the co-expressed unigene (AF260565) encoding KASIII was significantly down-regulated in the fatty acid pathway with a log2ratio of −1.10, -2.05, and −2.09 in the three contrast groups. So compared to 15 DAF, the differentially expressed levels of the gene encoding KASIII at 35DAF, 55DAF, 65DAF negatively correlated to fatty acid synthesis during the seed development, which is consistent with previous research [23]. The same down-regulated patterns were also observed for FabF (AF243183, CX702997, AK244605, AK285347), FabZ (BT098415), FabI (BI970856), FatB (AK244101) and oleoyl-[acyl-carrier-protein] hydrolase (EC 3.1.2.14, AK244101).

In most plant tissues, Acyl-ACP thioesterase is the major determinant of chain length and level of saturated fatty acids [24]. It plays an important role by influencing the fatty acid composition of the produced oil and then mainly the ratio of 16 C to 18 C fatty acids and the level of saturated fatty acid [25]. Two distinct but related thioesterase gene classes exist in higher plants: FatA is an acyl-specific thioesterase, with specificity for 18:1> > 18:0> > 16:0 fatty acids [26], which is considered an essential “housekeeping” enzyme in all plant cells; FatB is a thioesterase, which shows specificity for 16:0 > 18:1 > 18:0 fatty acids [27]. In this study, Solexa sequencing analysis showed that the expression of gene (CX706542) encoding FatA was unchanged at 35 DAF, but down-regulated at 55 DAF and 65 DAF with log2ratios of −1.7 and −1.5. While differentially co-expressed unigene (AK244101) encoding FatB showed constant decreased expressed levels with log2ratios of −1.4, -2.4, -2.7, respectively. Therefore, the expression level variations of CX706542 and AK244101 would influence the 16 and 18 carbon fatty acids synthesis, which needs to be confirmed by experiment.

Soybean seeds contain three lipoxygenase isozymes; lipoxygenases 1, 2, and 3. Lipoxygenases (linoleate: oxygen oxidoreductase; EC 1.13.11.12) catalyze the oxidation of unsaturated fatty acids to produce conjugated unsaturated fatty acid hydroperoxides, which are converted to volatile compounds associated with undesirable flavors [28]. Eliminating this enzyme from seeds could lead to better quality soybean protein and oil products. Three co-expressed lipoxygenase genes (X67304, U50081, D13949) identified in this study were among those with the highest expression levels. To our knowledge, this is the first report of a close correlation between lipoxygenase expression and fatty acid accumulation.

The P450 family is a large and diverse group of isozymes that mediate a diverse array of oxidative reactions. The activities of most of these enzymes are associated with biosynthetic processes such as phenylpropanoid, terpenoid, and fatty acid biosyntheses. Ten alkane-inducible P450 genes from Candida tropicalis (ATCC20336), which were responsible for omega-hydroxylation of n-alkanes and fatty acids, were cloned [29]. In their research, these enzymes were believed to be at least in part limiting in the conversion of fatty acid to diacids, but their relative oxidative activity toward other fatty acids was not known. The two unigenes (AF022464, DQ340246) encoding soybean P450 monooxygenase were identified in this research. Both of them showed down-regulated expressions during seed development, which indicated that they are negatively correlated to the fatty acid accumulation. For the other differentially co-expressed unigenes shown in Table 2, there is little information available about their relationship with lipid accumulation. Their log2ratios indicate that they have significant functions in regulating soybean seed development. However, our results cannot validate a direct relationship between these unigenes and oil accumulation. This topic requires further research.

Because we are interested in oil production in soybeans, we selected the fatty acid biosynthesis pathway for deep analysis from among the 124 biochemical pathways identified by Solexa sequencing. Twenty four co-expressed unigenes in the fatty acid pathway, including 13 that overlapped with the unigenes identified in the GO analysis, showed significant correlations with fatty acid accumulation (Table 3).

The up-regulated genes that were significantly correlated with fatty acid accumulation included ACCase and lipoxygenase. The down-regulated genes that were significantly correlated with fatty acid accumulation included FabD, KAS III, FabF, FabZ, FabI, FatB, FatA, P450, and oleoyl-[acyl-carrier-protein] hydrolase genes. FabG genes were both up- and down-regulated during seed development.

The analysis of the genes involved in the fatty acid biosynthetic pathway provides a basis to identify key regulatory processes controlling oil accumulation in soybean. However, biosynthetic pathways involve the cooperation of multiple genes. It is difficult to increase seed oil content by overexpressing a single gene. The large-scale characterization of unigenes described in this study shows comprehensive correlations between DEGs and fatty acid accumulation in soybean.

Conclusion

In this study, 11592, 16594, and 16255 differentially expressed unigenes were identified at 35, 55, and 65 days after flowering of soybean, respectively. Gene Ontology analyses detected 113 co-expressed unigenes associated with lipid biosynthesis. Of these, 15 showed significant changes in expression levels (log2fold values ≥ 1) during seed development. Pathway analysis revealed 24 co-expressed transcripts involved in lipid biosynthesis and fatty acid biosynthesis pathways. These results provide a comprehensive molecular biology background for research on soybean seed development, particularly with respect to the process of oil accumulation. All of the genes identified in our research have significance for breeding soybeans with increased oil contents.

Abbreviations

ACCase: acetyl CoA carboxylase; ACP: acyl carrier protein; BC: biotin carboxylase; BCCP: biotin carboxyl carrier protein; DAF: days after flowering; DE: differential expression; DEGs: differentially expressed genes; ER: endoplasmic reticulum; FabA: 3-hydroxydecanoyl-[acyl-carrier-protein] dehydratase; FabB: FabF, FabH or KASIII: 3-oxoacyl-[acyl-carrier-protein] synthase I, II, III; FabD or MCAT: [acyl-carrier-protein] S-malonyltransferase; FabG: 3-oxoacyl-[acyl-carrier protein] reductase; FabI: FabK, FabL: enoyl-[acyl-carrier protein] reductase I, II, III; FabZ: 3R-hydroxymyristoyl ACP dehydrase; FASN: Fas: fatty acid synthase; FAS1: 2: fatty acid synthase subunit beta, alpha; FatA: FatB: fatty acyl-ACP thioesterase A, B; FDR: false discovery rate; MdcH: malonate decarboxylase epsilon subunit; RPKM: reads/Kb/million; TAG: triacylglycerol; TL: total lipids.

Authors’ contributions

HC, FWW, YYD, NW, and YPS performed data analysis. HC, HYL wrote the manuscript. HYL and GC conceived the study. HC, HLY, YYJ, XYZ and YLL prepared the samples. HC, XYL, LL and XDF participated in the qRT-PCR verification. All the authors approved the final manuscript.

Acknowledgments

This work was supported by grants from the Program for the Special Program for Research of Transgenic Plants (Grant No. 2011ZX08010-002), the National Natural Science Foundation of China (Grant No. 30971804), the National High Technology Research and Development Program of China (863 Program) (Grant No. 2011AA100606), the New Century Excellent Talents in University (Grant No. NCET-08-0693), the Excellent innovation Team Project of Jilin Province, China (Grant No.20111815); the Science and Technology Development Project of Changchun City, Jilin, China (Grant No 2009024), and the State Key Laboratory of Crop Biology at Shandong Agricultural University, China (Grant No. 2010KF02).

References

  1. Zhou SQ, Zhang DJ, Luan HH, Yu FY, Xin XJ, Hu GH: Primary study on protein and lipid accumulation in high oil content soybean varieties.

    Chin. J. Oil Crop Sci 2006, 28(2):214-216. OpenURL

  2. Jako C, Kumar A, Wei YD, Zou JT, Barton DL, Giblin EM, Covello PS, Taylor DC: Seed-Specific Over-Expression of an Arabidopsis cDNA Encoding a Diacylglycerol Acyltransferase Enhances Seed Oil Content and Seed Weight1.

    Plant Physiol 2001, 126(2):861-874. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  3. Wang HW, Zhang B, Hao YJ, Huang J, Tian AG, Liao Y, Zhang JS, Chen SY: The soybean Dof-type transcription factor genes, GmDof4 and GmDof11, enhance lipid content in the seeds of transgenic Arabidopsis plants.

    The Plant Journal 2007, 52(4):716-729. PubMed Abstract | Publisher Full Text OpenURL

  4. Severin AJ, Woody JL, Bolon YT, Joseph B, Diers BW, Farmer AD, Muehlbauer GJ, Nelson RT, Grant D, Specht JE, Graham MA, Cannon SB, May GD, Vance CP, Shoemaker RC: RNA-Seq Atlas of Glycine max: A guide to the soybean transcriptome.

    BMC Plant Biol 2010, 10:160. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  5. Levin JZ, Berger MF, Adiconis X, Rogov P, Melnikov A, Fennell T, Nusbaum C, Garraway LA, Gnirke A: Targeted next-generation sequencing of a cancer transcriptome enhances detection of sequence variants and novel fusion transcripts.

    Genome Biol 2009, 10(10):R115. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  6. Zhang GJ, Guo GW, Hu XD, Zhang Y, Li QY, Li RQ, Zhuang RH, Lu ZK, He ZQ, Fang XD, Chen L, Tian W, Tao Y, Kristiansen K, Zhang XQ, Li SG, Yang HM, Wang J, Wang J: Deep RNA sequencing at single base-pair resolution reveals high complexity of the rice transcriptome.

    Genome Res 2010, 20(5):646-654. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  7. Jones SI, Gonzalez DO, Vodkin LO: Flux of transcript patterns during soybean seed development.

    BMC Genomics 2010, 11:136. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  8. Dong XL, Bai PL, Wang JM, Ruan CJ: Comparative study on determination of seed oil content of energy plants by using NMR and Soxhlet Extraction.

    Renewable Energy Resources 2011, 29(3):21-24. OpenURL

  9. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNA-Seq.

    Nat Methods 2008, 5(7):621-628. PubMed Abstract | Publisher Full Text OpenURL

  10. Nagalakshmi U, Wang Z, Waern K, Shou C, Raha D, Gerstein M, Snyder M: The Transcriptional Landscape of the Yeast Genome Defined by RNA Sequencing.

    Science 2008, 320(5881):1344-1349. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  11. Hu RB, Fan CM, Li HY, Zhang QZ, Fu YF: Evaluation of putative reference genes for gene expression normalization in soybean by quantitative real-time RT-PCR.

    BMC Molecular Biology 2009, 10:93. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  12. Yu J, Hu SN, Wang J, Wong GKS, Li SG, Liu B, Deng YJ, Dai L, Zhou Y, Zhang XQ, Cao ML, Liu J, Sun JD, Tang JB, Chen YJ, Huang XB, Lin W, Ye C, Tong W, Cong LJ, Geng JN, Han YJ, Li L, Li W, Hu GQ, Huang XG, Li WJ, Li J, Liu ZW, Li L, et al.: A Draft Sequence of the Rice Genome (Oryza sativa L. ssp. indica).

    Science 2002, 296(5565):79-92. PubMed Abstract | Publisher Full Text OpenURL

  13. Yu J, Wang J, Lin W, Li SG, Li H, Zhou J, Ni PX, Dong W, Hu SN, Zeng CQ, Zhang JG, Zhang Y, Li RQ, Xu ZY, Li ST, Li XR, Zheng HK, Cong LJ, Lin L, Yin JN, Geng JN, Li GY, Shi JP, Liu J, Lv H, Li J, Wang J, Deng YJ, Ran LH, Shi XL, et al.: The Genomes of Oryza sativa: A History of Duplications.

    PLoS Biol 2005, 3(2):e38. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  14. Xia QY, Zhou ZY, Lu C, Cheng DJ, Dai FY, Li B, Zhao P, Zha XF, Cheng TC, Chai CL, Pan GQ, Xu JS, Liu C, Lin Y, Qian JF, Hou Y, Wu ZL, Li GR, Pan MH, Li CF, Shen YH, Lan XQ, Yuan LW, Li T, Xu HF, Yang GW, Wan YJ, Zhu Y, Yu MD, Shen WD, et al.: A Draft Sequence for the Genome of the Domesticated Silkworm (Bombyx mori).

    Science 2004, 306(5703):1937-1940. PubMed Abstract | Publisher Full Text OpenURL

  15. Sasaki Y, Nagano Y: Plant acetyl-CoA carboxylase: structure, biosynthesis, regulation, and gene manipulation for plant breeding.

    Biosci Biotechnol Biochem 2004, 68(6):1175-1184. PubMed Abstract | Publisher Full Text OpenURL

  16. Roesler K, Shintani D, Savage L, Boddupalli S, Ohlrogge J: Targeting of the Arabidopsis homomeric acetyl-coenzyme A carboxylase to plastids of rapeseeds.

    Plant Physiol 1997, 113(1):75-81. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  17. Shintani D, Roesler K, Shorrosh B, Savage L, Ohlrogge J: Antisense expression and overexpression of biotin carboxylase in tobacco leaves.

    Plant Physiol 1997, 114(3):881-886. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  18. Thelen JJ, Mekhedov S, Ohlrogge JB: Brassicaceae express multiple isoforms of biotin carboxyl carrier protein in a tissue-specific manner.

    Plant Physiol 2001, 125(4):2016-2028. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  19. Elborough KM, Swinhoe R, Winz R, Kroon JT, Farnsworth L, Fawcett T, Martinez-Rivas JM, Slabas AR: Isolation of cDNAs from Brassica napus encoding the biotin-binding and transcarboxylase domains of acetyl-CoA carboxylase: assignment of the domain structure in a full-length Arabidopsis thaliana genomic clone.

    Biochem J 1994, 301(Pt 2):599-605. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  20. Roesler KR, Shorrosh BS, Ohlrogge JB: Structure and expression of an Arabidopsis acetyl-coenzyme a carboxylase gene.

    Plant Physiol 1994, 105(2):611-617. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  21. Simon JW, Slabas AR: cDNA cloning of Brassica napus malonyl-CoA: ACP transacylase (MCAT) (fabD) and complementation of an E. coli MCAT mutant.

    FEBS Lett 1998, 435(2–3):204-206. PubMed Abstract | Publisher Full Text OpenURL

  22. 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(2):1103-1114. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  23. Jaworski JG, Clough RC, Barnum SR: A cerulenin insensitive short chain 3-ketoacyl-acyl carrier protein synthase in Spinacia oleracea leaves.

    Plant Physiol 1989, 90(1):41-44. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  24. Salas JJ, Ohlrogge JB: Characterization of substrate specificity of plant FatA and FatB acyl-ACP thioesterases.

    Arch Biochem Biophys 2002, 403(1):25-34. PubMed Abstract | Publisher Full Text OpenURL

  25. Li-Beisson Y, Shorrosh B, Beisson F, Andersson MX, Arondel V, Bates PD, Baud S, Bird D, DeBono A, Durrett TP, Franke RB, Graham IA, Katayama K, Kelly AA, Larson T, Markham JE, Miquel M, Molina I, Nishida I, Rowland O, Samuels L, Schmid KM, Wada H, Welti R, Xu C, Zallot R, Ohlrogge J: Acyl-Lipid metabolism.

    The Arabidopsis Books 2010, 8:e0133. OpenURL

  26. Mandal MNA, Santha IM, Lodha ML, Mehta SL: Cloning of acyl-acyl carrier protein (ACP) thioesterase gene from Brassica juncea.

    Biochem Soc Trans 2000, 28(6):967-969. PubMed Abstract | Publisher Full Text OpenURL

  27. Dörmann P, Voelker TA, Ohlrogge JB: Accumulation of palmitate in Arabidopsis mediated by the acyl-acyl carrier protein thioesterase FATB1.

    Plant Physiol 2000, 123(2):637-644. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  28. Start WG, Ma Y, Polacco JC, Hildebrand DF, Freyer GA, Altschuler M: Two soybean seed lipoxygenase nulls accumulate reduced levels of lipoxygenase transcripts.

    Plant Mol Biol 1986, 7(1):11-23. Publisher Full Text OpenURL

  29. Craft DL, Madduri KM, Eshoo M, Wilson CR: Identification and Characterization of the CYP52 Family of Candida tropicalis ATCC 20336, Important for the Conversion of Fatty Acids and Alkanes to α, ω-Dicarboxylic Acid.

    App Environ Micro 2003, 9(10):5983-5991. OpenURL