An integrated omics analysis reveals molecular mechanisms that are associated with differences in seed oil content between Glycine max and Brassica napus

Zhang et al. BMC Plant Biology (2018) 18:328<br /> https://doi.org/10.1186/s12870-018-1542-8<br /> <br /> <br /> <br /> <br /> RESEARCH ARTICLE Open Access<br /> <br /> An integrated omics analysis reveals<br /> molecular mechanisms that are associated<br /> with differences in seed oil content<br /> between Glycine max and Brassica napus<br /> Zhibin Zhang1,2, Jim M. Dunwell3 and Yuan-Ming Zhang1*<br /> <br /> <br /> Abstract<br /> Background: Rapeseed (Brassica napus L.) and soybean (Glycine max L.) seeds are rich in both protein and oil, which<br /> are major sources of biofuels and nutrition. Although the difference in seed oil content between soybean (~ 20%) and<br /> rapeseed (~ 40%) exists, little is known about its underlying molecular mechanism.<br /> Results: An integrated omics analysis was performed in soybean, rapeseed, Arabidopsis (Arabidopsis thaliana L. Heynh),<br /> and sesame (Sesamum indicum L.), based on Arabidopsis acyl-lipid metabolism- and carbon metabolism-related genes.<br /> As a result, candidate genes and their transcription factors and microRNAs, along with phylogenetic analysis and<br /> co-expression network analysis of the PEPC gene family, were found to be largely associated with the difference<br /> between the two species. First, three soybean genes (Glyma.13G148600, Glyma.13G207900 and Glyma.12G122900)<br /> co-expressed with GmPEPC1 are specifically enriched during seed storage protein accumulation stages, while the<br /> expression of BnPEPC1 is putatively inhibited by bna-miR169, and two genes BnSTKA and BnCKII are co-expressed<br /> with BnPEPC1 and are specifically associated with plant circadian rhythm, which are related to seed oil biosynthesis. Then,<br /> in de novo fatty acid synthesis there are rapeseed-specific genes encoding subunits β-CT (BnaC05g37990D) and BCCP1<br /> (BnaA03g06000D) of heterogeneous ACCase, which could interfere with synthesis rate, and β-CT is positively regulated by<br /> four transcription factors (BnaA01g37250D, BnaA02g26190D, BnaC01g01040D and BnaC07g21470D). In triglyceride synthesis,<br /> GmLPAAT2 is putatively inhibited by three miRNAs (gma-miR171, gma-miR1516 and gma-miR5775). Finally, in rapeseed<br /> there was evidence for the expansion of gene families, CALO, OBO and STERO, related to lipid storage, and<br /> the contraction of gene families, LOX, LAH and HSI2, related to oil degradation.<br /> Conclusions: The molecular mechanisms associated with differences in seed oil content provide the basis for<br /> future breeding efforts to improve seed oil content.<br /> Keywords: Glycine max, Brassica napus, Acyl-lipid biosynthesis, Transcription factor, miRNA, Gene network<br /> <br /> <br /> Background is almost equal to the protein (~ 20%) and oil (~ 40%) con-<br /> Seed storage lipids not only provide food for human dietary tents in rapeseed [4]. As we know, most of the raw material<br /> consumption, but are also increasingly used as renewable required for seed oil and protein biosynthesis in rapeseed<br /> sources for biofuels [1, 2]. In oil crops, such as Arabidopsis, and soybean are derived from carbohydrate degradation<br /> soybean, rapeseed and sesame, seed oil content varies from [5]. And it should be noted that substrate competition be-<br /> 20 to 60%. Interestingly, the total seed storage reserves in tween seed oil and protein synthesis exists in oilseed crops<br /> soybean seed, consisting of ~ 20% oil and ~ 40% protein [3], [6, 7]. This is because phosphoenolpyruvate (PEP), a carbon<br /> compound derived from glycolysis, is not only used to<br /> synthesize acetyl-Coenzyme A (acetyl-CoA), which serves<br /> * Correspondence: soyzhang@mail.hzau.edu.cn<br /> 1<br /> Crop Information Center, College of Plant Science and Technology,<br /> as a substrate in the first step of de novo fatty acid synthe-<br /> Huazhong Agricultural University, Wuhan 430070, China sis, but is also required for the synthesis of oxaloacetate<br /> Full list of author information is available at the end of the article<br /> <br /> © The Author(s). 2018 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0<br /> International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and<br /> reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to<br /> the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver<br /> (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.<br /> Zhang et al. BMC Plant Biology (2018) 18:328 Page 2 of 15<br /> <br /> <br /> <br /> <br /> (OAA), which serves as a substrate in amino acid synthesis. can significantly increase soybean seed oil content [45].<br /> Thus, carbon metabolism is related to oil synthesis, and The down-regulation of CaFAD2 and CaFAE1 in crambe<br /> boosting the carbon flow to lipid synthesis can significantly with the FAD2-FAE1 RNAi vector led to a significant in-<br /> increase seed oil content [8]. crease in the seed oil to 80% compared to 13% for the wild<br /> In the past several decades, more than 700 acyl-lipid type [46]. Seed-specific simultaneous overexpression of<br /> metabolism-related genes and several hundred genes BnGPDH, BnGPAT and ScLPAAT genes in transgenic<br /> participating in carbohydrate metabolism have been iden- rapeseed may further enhance the desirable oil content<br /> tified in Arabidopsis thaliana [9, 10]. Among these genes, relative to single-gene overexpression [47]. Moreover, Yu<br /> more than 280 have been confirmed in A. thaliana et al. [48] developed a complete analysis platform of func-<br /> mutants as associated with acyl-lipid metabolism (http:// tional annotation for the soybean genes involved in<br /> aralip.plantbiology.msu.edu) [11]. Meanwhile, many genes acyl-lipid metabolism, and this makes the study of acyl<br /> have been experimentally validated to be closely related to lipid metabolism more efficient and accurate. However,<br /> seed oil content. For example, phosphoenolpyruvate none of the above studies were conducted at the whole<br /> carboxylase (PEPC) in cotton [12], acetyl-CoA carboxylase genome level.<br /> (ACCase) in rapeseed [13] and potato [14] participate in With the rapid development of sequencing technol-<br /> de novo fatty acid biosynthesis; fatty acylthioesterase B ogy, more and more plant genomes have been se-<br /> (GmFatB) in soybean [15] and patatin-related phospholip- quenced, and this accelerates the progress of research<br /> ase As (pPLAs) in Arabidopsis [16] are involved in fatty on acyl-lipid metabolism [49–53]. Troncoso-Ponce et<br /> acid elongation; glycerol-3-phosphate dehydrogenase al. [10] showed that the expression stoichiometry of<br /> (GPDH) in rapeseed [17], glycerol-3-phosphate acyltrans- most key lipid-related genes was relatively conserved<br /> ferase (GPAT) in Arabidopsis [18], 2-lysophosphatidic acid during seed development in Ricinus communis, Brassica<br /> acyltransferase (LPAAT) in Arabidopsis [19] and cotton napus, Euonymus alatus and Tropaeolum majus. Wang<br /> [20], acyl-CoA: diacylglycerol acyltransferase (DGAT) in et al. [54] analyzed the expression differences of key<br /> Arabidopsis [21, 22], maize [23], and rapeseed [24] are re- ovule-specific genes between non-fibrous Raymond’s<br /> lated to TAG synthesis; oleosins (OLE1) in Arabidopsis cotton and the upland cotton with fiber. Zhang et al.<br /> [25] participates in lipid droplet assembly and storage. In [55] dissected the molecular mechanisms of differences<br /> addition, some transcription factors (TFs) have been in seed oil content between four high-oil dicotyledons<br /> found to be associated with seed oil content, i.e., WRIN- and three low-oil grass plants. However, little is known<br /> KLED1 (WRI1) [26], LEAFY COTYLEDON1 (LEC1) [27, about the molecular mechanism for the difference of<br /> 28], LEAFY COTYLEDON2 (LEC2) [29], FUSCA3 seed oil contents between soybean and rapeseed.<br /> (FUS3) [30], GmDof4 and GmDof11 [31], GmbZIP123 To understand the molecular mechanisms of the differ-<br /> [32], GmMYB73 [33], GmDREBL [34], GmNFYA [35], ence of seed oil content between rapeseed and soybean,<br /> GmZF351 [36], and ABSCISIC ACID INSENSITIVE3 an integrated omics analysis was performed in Arabidop-<br /> (ABI3) [37, 38]. However, all the above studies involved sis, rapeseed, soybean and sesame based on Arabidopsis<br /> only a single lipid-related gene or transcription factor. acyl-lipid metabolism- and carbon metabolism-related<br /> Seed oil content is typically a quantitative trait regulated genes. The integrated omics analysis included gene copy<br /> by multiple genes. As we know, these genes have been number variation, expression pattern, microRNA/tran-<br /> identified in the form of quantitative trait loci in soybean scription factor, phylogenetic and co-expressional network<br /> and rapeseed in the past decades [39–43]. analyses. Thus, candidate genes, transcription factors, and<br /> In reality, acyl-lipid metabolism is a complex biological microRNAs that may be responsible for the difference<br /> process that includes at least conversion of sucrose to were identified. These results provide a novel explanation<br /> pyruvate, plastidial de novo fatty acid (FA) synthesis, for differences at the whole genome level and the basis for<br /> endoplasmic triacylglycerol (TAG) biosynthesis, and future breeding efforts to improve seed oil content in oil-<br /> oil-body assembly. It is therefore important to determine seed crops.<br /> whether specific combination of multiple genes from<br /> multiple metabolic pathways can increase seed oil con- Results<br /> tent more effectively as compared with the manipulation Identification of candidate genes related to lipid<br /> of an individual gene. For example, it was found that biosynthesis<br /> Arabidopsis seed-specific overexpression of WRI1 and To identify orthologous genes in soybean and rape-<br /> DGAT1 combined with suppression of SDP1 leads to seed, we used OrthoMCL to cluster putative OGs of<br /> higher seed oil content than the manipulation of each genes across Arabidopsis, soybean, rapeseed and<br /> gene individually [44]. Additionally, the simultaneous sesame. As a consequence, 172,626 (81.83%) protein-<br /> overexpression of GmFabG (Glyma.12G092900), GmACP coding genes from the four species were clustered into<br /> (Glyma.09G060900) and GmFAD8 (Glyma.03G056700) 27,236 OGs (Additional file 1: Table S1), with each<br /> Zhang et al. BMC Plant Biology (2018) 18:328 Page 3 of 15<br /> <br /> <br /> <br /> <br /> group representing a gene family. Among these gene and PKp involved in carbon metabolism from sucrose to<br /> families, 11,314 (41.54%) were defined as rapeseed- pyruvate; PDK1, ACCase, KASII, HAD, KAR, FATA,<br /> specific paralogous gene clusters without soybean SAD and FAD2 in de novo fatty acid biosynthesis; and<br /> genes, and only 3391 (12.45%) were soybean-specific PAP, PDCT participating in TAG synthesis, as well as<br /> families (Additional file 1: Table S2). In addition, 1687 oil-body proteins OBO, CALO, STERO, and oil degrad-<br /> OGs were identified to have one copy of a soybean ation genes LOX, LAH, HSI2 and DSEL (Fig. 1).<br /> gene and multiple copies for rapeseed genes, and 446 To further determine the functions of the above 484<br /> OGs with one copy for the rapeseed gene and multiple genes, KEGG enrichment analysis was conducted using<br /> copies for soybean genes (Additional file 1: Table S2). KOBAS 2.0 [59]. As a result, the top 10 KEGG pathways<br /> To date, more than 700 Arabidopsis acyl-lipid metab- for soybean and rapeseed candidate genes were obtained<br /> olism genes have been detected, of which 135 are dir- (Table 1). It was found that the two crops had eight<br /> ectly involved in the processes of de novo fatty acid KEGG pathways in common, namely pyruvate metabol-<br /> synthesis, triglyceride biosynthesis and lipid droplet for- ism, carbon metabolism, biosynthesis of secondary me-<br /> mation. Using the method as described by Troncoso-P tabolites, glycolysis, purine metabolism, biosynthesis of<br /> once et al. [10], we also extracted 238 Arabidopsis amino acids, carbon fixation in photosynthetic organ-<br /> carbon-metabolism genes. Subsequently, searches using the isms, and glycerophospholipid metabolism. In addition,<br /> 373 genes as queries were performed to obtain lipid fatty acid biosynthesis in soybean is similar to fatty acid<br /> biosynthesis-related homologous gene families (Additional metabolism in rapeseed. The results are consistent with<br /> file 1: Table S3). As a result, 230 putative OGs related to those in Troncoso-Ponce et al. [10] and Ohlrogge and<br /> lipid synthesis were identified, which include 781 soybean Browse [60], and thus ensure the reliability of candidate<br /> genes and 1267 rapeseed genes (Additional file 1: Table S4). genes in the next analysis.<br /> <br /> Copy number variation and expression clustering of Expression profiles of candidate genes responsible for the<br /> candidate genes related to lipid biosynthesis difference of seed oil content between rapeseed and<br /> In order to eliminate the effect of species ploidy, the soybean<br /> relative copy number of a gene was used to measure the The transcriptomic datasets from the seed developmental<br /> difference of copy number of homologous genes be- stages in soybean (GSE42871) and rapeseed (GSE77637)<br /> tween species. Of the above 230 OGs related to lipid downloaded from the GEO (Gene Expression Omnibus)<br /> synthesis, 44 were found to have differences in relative database were used to validate the above candidate genes.<br /> copy number of a gene between soybean and rapeseed To compare the expressional profiles of each candidate<br /> (Additional file 1: Table S5). gene in soybean and rapeseed, relative expression content<br /> To cluster and visualize the expressional patterns of for each gene was adopted in this study; this is defined as<br /> all 2048 rapeseed and soybean genes in the 230 OGs, the ratio of the expression of each gene to average expres-<br /> we exploited STEM software [56] to analyze the expres- sion of all the genes in the species. As a result, a majority<br /> sion data of four seed development stages in rapeseed of candidate genes in rapeseed, except for phosphoenol-<br /> (GSE77637, [57]) and soybean (GSE42871, [58]). In this pyruvate carboxylase (PEPC), Ribulose-1,5-bisphosphate<br /> study, R3, R4, R7 and R8 stages in soybean and 2, 4, 6, carboxylase/oxygenase small subunit (RBCS1A), lipoxy-<br /> and 8 weeks after pollination (WAP) in rapeseed were genase (LOX) and steroleosin (STERO), had higher rela-<br /> defined as t1’, t2’, t3’ and t4’ in soybean and t1, t2, t3 tive expression than those in soybean, especially for PK<br /> and t4 in rapeseed, respectively. Results showed that and ACCase (Additional file 2: Figure S3, Additional file 1:<br /> 2048 genes were grouped into 20 clusters, including Table S6).<br /> three up-regulation patterns (cluster13, cluster16 and More importantly, we noted some interesting phe-<br /> cluster18) and one down-regulation profile (cluster3) dur- nomena. First, GmPEPC had higher relative expression<br /> ing stages with rapid accumulation of seed oil (Additional at the early and middle seed development stages than<br /> file 2: Figure S1). Additionally, there were 23, 202, 86 and BnPEPC (Fig. 3), indicating that PEP may be more likely<br /> 22 genes in cluster3, cluster13, cluster16 and cluster18, to be used to synthesize protein in soybean seed,<br /> respectively (P-value < 0.05) (Additional file 2: Figure S2). because PEPC, a member of carboxyl lyase family, cata-<br /> Based on the above two results, 192 soybean and 292 lyzes phosphoenolpyruvate (PEP) to produce oxaloacetic<br /> rapeseed candidate genes were inferred to be related to acid (OAA) for amino acid biosynthesis. Then, the<br /> the differences of seed oil content between rapeseed and relative expression contents of rapeseed genes encoding<br /> soybean. According to the sequence homology with four subunits (α-CT, β-CT, BC and BCCP) of heteroge-<br /> Arabidopsis genes, 484 genes were found to putatively neous acetyl-CoA carboxylase (ACCase), catalyzing the<br /> encode a series of core enzymes. For example, GRF2 and first and committed reaction of de novo fatty acid bio-<br /> RBCS1A during photosynthesis; PGK, ApS1, SUC, PEPC synthesis in plastids, were higher than those of soybean<br /> Zhang et al. BMC Plant Biology (2018) 18:328 Page 4 of 15<br /> <br /> <br /> <br /> <br /> Fig. 1 Candidate genes for the difference of seed oil content between soybean and rapeseed. a Relative copy number variation analysis of candidate<br /> genes contributed to the difference of seed oil contents between soybean and rapeseed. The red genes in X-axis indicate more relative gene copies<br /> in rapeseed over in soybean, and the opposite situation is expressed in black. b PK in plastid is composed of Alpha (α) and Beta (β) subunits. ACCase<br /> contains homogeneous structure ACC2 and heterogeneous ACCase complex including α-CT, β-CT, BC, BCCP1 and BCCP2 subunits<br /> <br /> <br /> genes in the oil rapid accumulation stages. Especially, zinc finger TFs, namely BnaA01g37250D, BnaA02g2619<br /> BCCP1 and β-CT were not expressed during soybean 0D, BnaC01g01040D and BnaC07g21470D (Table 2;<br /> seed development (Additional file 2: Figure S4). Finally, Additional file 2: Figure S5) [62].<br /> GmPEPC1 (Glyma.06G277500) and GmPEPC3 (Glyma. Meanwhile, 639 soybean and 92 rapeseed mature<br /> 06G229900) had higher relative expression than PKp microRNAs were downloaded from miRBase (version<br /> and ACCase in soybean. Conversely, PKp-β (BnaC02g4 21). Their target genes were predicted using psRNATar-<br /> 4850D), PKp-α (BnaA01g24280D) and ACCase had get, a plant small RNA target analysis server. As a result,<br /> higher relative expression than PEPC in rapeseed devel- 4411 soybean and 1780 rapeseed miRNA-Target gene<br /> opment stages (Additional file 2: Figure S4). Therefore, pairs were obtained. Among these genes, 116 and 61<br /> we deduced that PEPC, PKp and ACCase are most likely were associated with lipid synthesis in soybean and rape-<br /> to be the key genes that regulate the distribution of car- seed, respectively. Note that bna-miR169 inhibits the ex-<br /> bon sources in soybean and rapeseed seeds. pression of the BnPEPC gene; this may facilitate a<br /> greater carbon flow to de novo fatty acid synthesis, and<br /> Transcription factors and microRNAs regulatory network the expression of GmLPAAT2 gene is putatively inhib-<br /> analysis of the candidate genes ited by gma-miR171, gma-miR1516 and gma-miR5775<br /> To clarify the differences of regulatory networks of key (Table 2; Additional file 2: Figure S5).<br /> candidate genes in soybean and rapeseed, we identified<br /> transcription factors (TFs) related to lipid biosynthesis Evolutionary analysis of PEPC gene family<br /> in seed development. TFs and their target genes were Phylogenetic analysis and conserved motifs analysis of<br /> downloaded from PlantTFDB v3.0 [61]. Together with PEPC gene family<br /> the results mentioned above, it was found that the Three plant-type PEPC genes (PTPCs) (AtPEPC1,<br /> rapeseed-specific gene (BnaA10g13960D), encoding AtPEPC2 and AtPEPC3) and one bacterial-type PEPC<br /> β-CT subunit of ACCase, is positively regulated by four genes (BTPCs) (AtPEPC4) exist in A. thaliana [63]. To<br /> Zhang et al. BMC Plant Biology (2018) 18:328 Page 5 of 15<br /> <br /> <br /> <br /> <br /> Table 1 KEGG pathway enrichment analysis for candidate genes related to the differences of seed oil content between rapeseed<br /> and soybean<br /> KEGG pathways Soybean Rapeseed<br /> ID Input Background P-value Corrected ID Input Background P-value Corrected<br /> number number P-value number number P-value<br /> Pyruvate metabolism ath00620 43 85 4.92E-80 3.65E-77 ath00620 80 85 1.03E-133 1.64E-132<br /> Carbon metabolism ath01200 41 262 6.29E-58 2.34E-55 ath01200 75 262 2.99E-94 2.39E-93<br /> Biosynthesis of secondary ath01110 42 1076 1.04E-35 3.52E-34 ath01110 85 1076 3.90E-65 1.25E-64<br /> metabolites<br /> Glycolysis ath00010 21 117 3.37E-30 4.24E-29 ath00010 38 117 1.89E-48 5.03E-48<br /> Purine metabolism ath00230 21 158 1.06E-27 9.51E-27 ath00230 38 158 3.43E-44 7.84E-44<br /> Biosynthesis of amino acids ath01230 21 255 1.18E-23 7.77E-23 ath01230 38 255 3.31E-37 5.29E-37<br /> Carbon fixation in photosynthetic ath00710 12 69 1.97E-17 9.34E-17 ath00710 17 69 2.36E-20 2.90E-20<br /> organisms<br /> Glycerophospholipid metabolism ath00564 15 86 1.43E-21 8.06E-21 ath00564 30 86 4.48E-39 7.97E-39<br /> Fatty acid biosynthesis ath00061 17 41 6.76E-30 8.37E-29<br /> Fatty acid metabolism ath01212 46 67 8.95E-71 3.58E-70<br /> Glycerolipid metabolism ath00561 11 53 8.03e-17 3.53e-16<br /> Biotin metabolism ath00780 8 14 1.32E-12 1.40E-12<br /> <br /> <br /> <br /> investigate the evolution of the PEPC gene family in soy- 2361.163 if the evolution rates changed across different<br /> bean and rapeseed, the full amino acid sequences branches (Additional file 3: Table S7; Additional file 2:<br /> encoded by 33 PEPC genes in soybean, rapeseed and Figure S6). Clearly, there was a significant difference be-<br /> Arabidopsis were used to construct a phylogenetic tree tween the above two models (Additional file 3: Table S7,<br /> using a neighbor-joining method. As a result, the PEPC P-value = 4.278e-06). We also found that the ω value for<br /> genes were grouped into two distinct families with four the bacterial-type PEPC soybean sub-branch (ω = 0.495)<br /> subfamilies, which are consistent with those in A. thali- was significantly higher than the ω0 value (= 0.340) and<br /> ana (Fig. 2). In Fig. 2, GmPEPC1 is close to GmPEPC3 the ω values for the other sub-branches were signifi-<br /> in its evolutionary relationship, and their expression pat- cantly lower than the ω0 value (= 0.340). This indicates<br /> terns during seed development in Additional file 2: Fig- that the BTPC genes have experienced positive selection<br /> ure S4 are complementary. Meanwhile, the conserved and the PTPC genes have experienced purifying selec-<br /> motifs of PEPC genes were further analyzed using tion. Furthermore, the BSM model was used to identify<br /> MEME [64]. Results showed that the motif structure of positively selected sites. As a result, we observed signifi-<br /> GmPEPC3 was more conserved than that of BnPEPC3, cant positive selection for soybean bacterial-type PEPC<br /> and the motif structures of BnPEPC1 and BnPEPC2 were genes and the amino acid sites of positive selection, as-<br /> also relatively more conserved than that of BnPEPC3 (Fig. sociated with GmPEPC4, were 56F and 61 V (Additional<br /> 2). Note that there are distinct differences of BTPCs be- file 3: Table S8).<br /> tween soybean and rapeseed, indicating the existence of<br /> its extensive functional differentiation. Differential analysis of PEPC1 gene co-expression<br /> networks between soybean and rapeseed<br /> Evolutionary rate and positive selection analysis of PEPC The co-expression network of one gene is frequently<br /> genes constructed by Pearson’s correlation coefficient [65, 66].<br /> To determine whether the genes of the PEPC gene fam- In the present study this method was used to construct<br /> ily are under different evolutionary constraints in soy- the co-expression networks of the PEPC1 gene in soy-<br /> bean and rapeseed, the ω (Ka/Ks) values for the above bean and rapeseed. The differences between the two net-<br /> genes were calculated using the branch model (BM) and works were also used to identify extra candidate genes.<br /> the branch-site model (BSM) of the Codeml program in As a result, 121 soybean and 133 rapeseed genes were<br /> PAML. As a result, the evolution rate ω0 (Ka/Ks) was es- co-expressed with GmPEPC1 and BnPEPC1, respect-<br /> timated to be 0.340 and log-likelihood was − 2343.655 if ively. Among these co-expressed genes, 17 were ortholo-<br /> the evolution rates at all branches were the same; six gous. The other genes were used to conduct KEGG<br /> evolution rates (ω) were estimated to be 0.078, 0.054, pathway enrichment analysis. In the top 10 KEGG path-<br /> 0.106, 0.093, 0.495 and 0.175 and log-likelihood was − ways for soybean or rapeseed, the soybean-specific<br />  Zhang et al. BMC Plant Biology<br /> <br /> <br /> <br /> <br /> Table 2 The microRNAs regulating rapeseed PEPC gene (BnaA03g33640D) and soybean LPAAT2 gene (Glyma.03G139700)<br /> miRNA_Acc. Target_Acc. Expectation UPE miRNA_start miRNA_end Target_start Target_end miRNA_aligned_fragment Target_aligned_fragment Inhibition<br /> bna-miR169a BnaA03g33640D 4 21.988 1 20 204 223 CAGCCAAGGAUGACUUGCCG AGGCAAGCCAAACUUGGCUG Translation<br /> bna-miR169b BnaA03g33640D 4 21.988 1 20 204 223 CAGCCAAGGAUGACUUGCCG AGGCAAGCCAAACUUGGCUG Translation<br /> (2018) 18:328<br /> <br /> <br /> <br /> <br /> bna-miR169c BnaA03g33640D 3.5 21.988 1 20 204 223 UAGCCAAGGAUGACUUGCCU AGGCAAGCCAAACUUGGCUG Translation<br /> bna-miR169d BnaA03g33640D 3.5 21.988 1 20 204 223 UAGCCAAGGAUGACUUGCCU AGGCAAGCCAAACUUGGCUG Translation<br /> bna-miR169e BnaA03g33640D 3.5 21.988 1 20 204 223 UAGCCAAGGAUGACUUGCCU AGGCAAGCCAAACUUGGCUG Translation<br /> bna-miR169f BnaA03g33640D 3.5 21.988 1 20 204 223 UAGCCAAGGAUGACUUGCCU AGGCAAGCCAAACUUGGCUG Translation<br /> bna-miR169g BnaA03g33640D 3.5 21.988 1 22 202 223 UAGCCAAGGAUGACUUGCCUGC GAAGGCAAGCCAAACUUGGCUG Translation<br /> bna-miR169h BnaA03g33640D 3.5 21.988 1 22 202 223 UAGCCAAGGAUGACUUGCCUGC GAAGGCAAGCCAAACUUGGCUG Translation<br /> bna-miR169i BnaA03g33640D 3.5 21.988 1 22 202 223 UAGCCAAGGAUGACUUGCCUGC GAAGGCAAGCCAAACUUGGCUG Translation<br /> bna-miR169j BnaA03g33640D 3.5 21.988 1 22 202 223 UAGCCAAGGAUGACUUGCCUGC GAAGGCAAGCCAAACUUGGCUG Translation<br /> bna-miR169k BnaA03g33640D 3.5 21.988 1 22 202 223 UAGCCAAGGAUGACUUGCCUGC GAAGGCAAGCCAAACUUGGCUG Translation<br /> bna-miR169l BnaA03g33640D 3.5 21.988 1 22 202 223 UAGCCAAGGAUGACUUGCCUGC GAAGGCAAGCCAAACUUGGCUG Translation<br /> bna-miR169n BnaA03g33640D 4 21.988 1 20 204 223 CAGCCAAGGAUGACUUGCCG AGGCAAGCCAAACUUGGCUG Translation<br /> gma-miR171b-3p Glyma.03G139700 3 24.213 1 20 168 186 CGAGCCGAAUCAAUAUCACU AGUGAUAUUGAUU-GGCUUG Cleavage<br /> gma-miR1516a-5p Glyma.03G139700 3 24.131 1 23 283 305 CAAGUUAUAAGCUCUUUUGAGAG CUCUCAAAAGCACUUAUGGCUUG Cleavage<br /> gma-miR1516b Glyma.03G139700 1 14.793 1 21 264 284 AGCUUCUCUACAGAAAAUAUA UAUAUUUUCAGUAGAGAAGCU Cleavage<br /> gma-miR5775 Glyma.03G139700 3 23.458 1 21 279 299 AUAAGCUCUUUUGAGAGCUUC GAAGCUCUCAAAAGCACUUAU Cleavage<br /> gma-miR171b-3p Glyma.19G142500 3 22.553 1 20 210 228 CGAGCCGAAUCAAUAUCACU AGUGAUAUUGAUU-GGCUUG Cleavage<br /> Page 6 of 15<br /> Zhang et al. BMC Plant Biology (2018) 18:328 Page 7 of 15<br /> <br /> <br /> <br /> <br /> Fig. 2 Phylogenetic tree and gene motif analysis of PEPC gene family. a Phylogenetic tree of PEPC gene family is constructed from the complete<br /> alignment of 33 PEPC protein sequences in Arabidopsis, soybean, and rapeseed using the neighbor-joining method with 1000 bootstrap replicates<br /> with the MEGA 7.0 software program. The bootstrap scores are indicated on the nodes, and the 4 PEPC branches, all of which are based on<br /> Arabidopsis PEPC orthologous genes, are indicated in four color boxes. The representative gene of each branch is shown followed by an additional<br /> abbreviation. b Conserved domains analysis of PEPCs. The domains of soybean genes in AtPEPC3 branch are relatively more conservative compared<br /> with rapeseed genes. And the domains of rapeseed genes in AtPEPC1 and AtPEPC2 branches are relatively conservative. However, the differences of<br /> gene domain in the bacterial AtPEPC4 branch are significant between soybean and rapeseed<br /> <br /> <br /> biological process is “the synthesis of valine, leucine and content between the two species (Fig. 4). The reasons<br /> isoleucine”, involving Glyma.13G148600, Glyma.13G207 are as follows.<br /> 900 and Glyma.12G122900, and rapeseed-specific bio- First, GmPEPC1 has higher relative expression at the<br /> logical process is “plant circadian rhythms”, involving early and middle stages of seed development than<br /> genes BnSTKA (BnaC08g48660D, BnaA09g42220D, Bna BnPEPC1, and bna-miR169 putatively inhibits the ex-<br /> A01g21040D, BnaC08g34660D, BnaC01g42660D) and pression of BnPEPC. Although the bacterial-type PEPC<br /> BnCKII (BnaC08g30500D, BnaC02g33100D, BnaC04g0 (BTPC) gene in A. thaliana can inhibit the expression<br /> 5080D, BnaA02g24960D) (Fig. 3). of the plant-type PEPC (PTPC) gene [63, 70–73], BTPC<br /> genes in soybean have experienced positive selection<br /> Discussion (Additional file 3: Tables S7 and S8) and it is possible<br /> PEPC, along with its miRNA and co-expressed genes, to lose the function of inhibiting the expression of<br /> which can affect the flow of carbon sources in seeds, may PTPC gene (GmPEPC1). In addition, Xu et al. [12]<br /> contribute to the difference of seed oil content between increased the accumulation of cotton seed oil by the<br /> soybean and rapeseed down-regulation of GhPEPC1 via RNA interference in<br /> Seed oil content is almost negatively correlated to seed transgenic cotton plants. These studies provide<br /> protein content in soybean and rapeseed [67–69]. As we evidence for greater carbon flow to amino acid metab-<br /> know, PEP is used to synthesize acetyl-CoA under the olism in soybean seed and to de novo fatty acid synthe-<br /> catalysis of pyruvate kinase (PK) and acetyl coenzyme A sis in rapeseed seed. This may partly explain why there<br /> carboxylase (ACCase) so that the PEP enters into the are high seed protein content in soybean and high seed<br /> fatty acid synthesis pathway. Additionally, PEP is also oil content in rapeseed.<br /> used to synthesize oxaloacetate (OAA) under the cataly- Secondly, the expression of LPAAT2-encoding gene,<br /> sis of phosphoenolpyruvate carboxylase (PEPC) so that involved in triacylglycerol synthesis in soybean seed, is<br /> the PEP enters into the amino acid synthesis pathway. putatively inhibited by three miRNAs (gma-miR171,<br /> Results in this study showed that PEPC genes, together gma-miR1516 and gma-miR5775) based on the results<br /> with their miRNA and co-expressed genes, may increase of bioinformatics analysis (Table 2).<br /> the flow of carbon to the biosynthesis of amino acids in Finally, gene co-expression network analysis helps us<br /> soybean seed and to the de novo fatty acid synthesis in to understand different biological pathways in soybean<br /> rapeseed seed, resulting in the difference of seed oil and rapeseed. KEGG enrichment analyses of genes<br /> Zhang et al. BMC Plant Biology (2018) 18:328 Page 8 of 15<br /> <br /> <br /> <br /> <br /> Fig. 3 Comparison of PEPC1 gene co-expression networks in rapeseed (a) and soybean (b). 17 genes in light pink are orthologous genes in rapeseed<br /> and soybean. More than 100 genes in red were enriched in the same processes based on KEGG pathway enrichment analysis. The blue<br /> nodes represent rapeseed genes enriched-specific in “plant circadian rhythms” (a) and soybean genes enriched-specific in “Valine, leucine<br /> and isoleucine biosynthesis” (b), respectively<br /> <br /> <br /> <br /> <br /> Fig. 4 Molecular mechanisms for the difference of seed oil content between soybean and rapeseed. Candidate genes contributed to the differences<br /> of seed oil content between soybean and rapeseed obtained in the study were marked with red color. GRF2 and RBCS1A: photosynthesis; PGK, ApS1,<br /> SUC, PEPC and PKp: carbon metabolism from sucrose to pyruvate; PDK1, ACCase, KASII, HAD, KAR, FATA, SAD and FAD2: in de novo fatty<br /> acid biosynthesis; PAP and PDCT: TAG synthesis; OBO, CALO and STERO: oil-body protein genes; LOX, LAH, HSI2 and DSEL: oil degradation<br /> genes. Among these candidate genes, BCCP1 (BnaA03g06000D) and β-CT (BnaC05g37990D) in heterogeneous acetyl-CoA carboxylase (ACCase) are<br /> rapeseed-specific genes, and β-CT is positively regulated by four transcription factors (BnaA01g37250D, BnaA02g26190D, BnaC01g01040D and<br /> BnaC07g21470D). The gene expression of PEPC1 in rapeseed is putatively inhibited by bna-miR169, while LPAAT in soybean putatively inhibited<br /> by gma-miR171, gma-miR1516 and gma-miR5775 in triglyceride synthesis. The pink genes are speculated related specifically to high seed oil<br /> content of rapeseeds, and the purple speculated specifically related to high seed protein content in soybean, which were both identified by<br /> PEPC co-expression network analysis. Soybean genes participated in Branched-Chain Amino Acid (BCAA) synthesis may contribute to seed high protein content<br /> by adjusting the flow of PEP and downstream protein biosynthesis. Rapeseed genes BnSTKA and BnCKII are likely to promote the triglyceride<br /> synthesis by phosphorylating circadian TFs cca1/lhy and thus increase the seed oil content of rapeseed<br /> Zhang et al. BMC Plant Biology (2018) 18:328 Page 9 of 15<br /> <br /> <br /> <br /> <br /> co-expressed with GmPEPC1 showed that three soybean predictably, which is consistent with the results of Jin et<br /> genes (Glyma.13G148600, Glyma.13G207900 and al. [62]. Similarly, Li et al. [36] demonstrated that overex-<br /> Glyma.12G122900) were enrich-specific in the “leucine, pression of GmZF351, a gene encoding a tandem CCCH<br /> isoleucine and valine” synthesis pathway. Leucine, iso- zinc finger protein, can activate lipid biosynthesis genes<br /> leucine and valine are the three major branched-chain and increase seed oil accumulation in soybean. Moreover,<br /> amino acids for protein synthesis. The content of Li et al. [84] found that transfer DNA insertional alleles<br /> branched-chain amino acids in seeds is positively corre- that completely eliminate the accumulation of BCCP2<br /> lated with the protein content in general, which can have no perceptible effect on fatty acid accumulation,<br /> effectively maintain the accumulation of storage proteins while reducing the BCCP1 accumulation can dramatic-<br /> in seeds [74]. Therefore, it is speculated that the expres- ally decreases fatty acid accumulation in Arabidopsis<br /> sion of such soybean genes may be beneficial for the thaliana. This partly supports that rapeseed-specific<br /> accumulation of storage proteins in seeds. Meanwhile, BCCP1 gene may associate with high seed oil content of<br /> KEGG enrichment analysis of genes co-expressed with rapeseed. It should also note that RNA levels don’t always<br /> BnPEPC1 revealed that nine rapeseed genes, encoding equate to protein and/or lipid metabolite levels in plants<br /> serine/threonine-protein kinase (BnSKTA; BnaC08g48 [85].<br /> 660D, BnaA09g42220D, BnaA01g21040D, BnaC08g3466<br /> 0D and BnaC01g42660D) and Casein kinase II subunit The expansion of gene families associated with lipid<br /> beta (BnCKII; BnaC08g30500D, BnaC02g33100D, BnaC storage and the contraction of gene families related to<br /> 04g05080D and BnaA02g24960D), were specifically lipid degradation may contribute to high seed oil content<br /> enriched in the “plant circadian rhythm” category, which in rapeseed<br /> can regulate seed oil metabolism and hormone signaling Seed triglyceride is mainly stored in lipid droplets, and<br /> pathway [75, 76]. Lipid metabolism is subject to diurnal the size of the lipid droplets and the spatial distribution<br /> regulation at the early stages of seed development in of their assembly proteins affect seed oil content [86,<br /> Arabidopsis [77]; diurnal differential expression of genes 87]. In this study, it was found that the relative copy<br /> encoding JcDof1, a dof TF of Jatropha curcas in numbers of genes encoding STERO, CALO and OLEs in<br /> response to light signal, β-hydroxy-3-methylglutaryl- rapeseed are significantly higher than those in soybean<br /> CoA reductase and Cyp7A1, regulates seed oil synthesis (Fig. 1, Additional file 1: Table S6), and such genes in<br /> and accumulation [78, 79]. CIRCADIAN CLOCK ASSO- rapeseed are obviously up-regulated during stages of<br /> CIATED1 (CCA1) and LATE ELONGATED HYPO- rapid lipid accumulation (t2~t3) (Additional file 2:<br /> COTYL (LHY) from the core clock system can affect Figure S3, Additional file 1: Table S6). On the other<br /> the reserve mobilization of storage lipid [77], but this hand, the gene families LOX, LAH and HSI2, related to<br /> process is affected by the phosphorylation of protein lipid degradation, have contracted in rapeseed. In other<br /> kinase (CK2) [80]. Therefore, the phosphorylation of words, the relative copy numbers of these genes are<br /> genes BnSTKA and BnCKII may promote the storage of much smaller than those in soybean (31/4 < 41/2, 10/4 <<br /> seed oil. 24/2 and 8/4 < 8/2, (gene absolute copy numbers) /<br /> (species polyploidy)). This relationship was also found<br /> Rapeseed-specific genes encoding β-CT and BCCP1 between soybean and sesame. Specifically, this latter spe-<br /> subunits of acetyl-CoA carboxylase and transcription cies shows contraction of gene families (LOX, LAH and<br /> factors may be associated with higher seed oil content FAR1) related to lipid degradation and expansion of<br /> in rapeseed gene families (LTP1 and SUT) related to lipid storage<br /> β-CT and BCCP are important components for hetero- [51]. Therefore, it was speculated that the contraction of<br /> geneous acetyl-CoA carboxylase (ACCase) [81, 82]. gene families related to lipid degradation and the expan-<br /> Overexpression of ACCase subunit genes can signifi- sion of gene families related to lipid storage may be an<br /> cantly increase fatty acid content in oil crop seed [13, important reason for the higher seed oil content in rape-<br /> 14, 83]. In this study, soybean β-CT and BCCP1 was not seed than in soybean.<br /> expressed, while BCCP2, BC and α-CT showed high ex- In order to further ascertain whether the degradation<br /> pression during seed development stages (Additional file of seed storage materials in oilseed crop is specialized in<br /> 2: Figure S4). This is consistent with the results of Zhang an evolutionary sense, we investigated gene families<br /> et al. [55]. Meanwhile, genes BCCP1, BCCP2, BC, α-CT related to protein degradation in soybean, rapeseed and<br /> and β-CT showed high expression during rapeseed seed sesame seeds. As we know, the degradation of protein in<br /> development stages. Especially, the β-CT subunit gene plant cells is mainly mediated by the ubiquitin prote-<br /> (BnaA10g13960D) was positively regulated by four zinc asome, lysosomal and caspase pathways. Among the<br /> finger C2H2 transcription factors (BnaA01g37250D, three pathways, the ubiquitination proteasome pathway<br /> BnaA02g26190D, BnaC01g01040D and BnaC07g21470D) is the main pathway of storage protein degradation in<br /> Zhang et al. BMC Plant Biology (2018) 18:328 Page 10 of 15<br /> <br /> <br /> <br /> <br /> oilseed crop seeds [88], and it mainly involves ubiquitin- content. Finally, the expansion of gene families related<br /> activating enzyme (E1), ubiquitin-conjugating enzyme to lipid storage, and the contraction of gene families re-<br /> (E2) and ubiquitin ligases (E3) [89]. In this study, we lated to oil degradation may play important roles on the<br /> found that the relative copy numbers of genes encoding difference in seed oil content.<br /> E1, E2 and E3 in soybean, rapeseed and sesame were 4/<br /> 2 = 8/4 > 2/2, 44/2 < 105/4 > 23/2 and 33/2 > 52/4 > 10/2,<br /> Methods<br /> which are not consistent with the protein contents in<br /> Data sources<br /> soybean, rapeseed and sesame seeds (~ 40%, ~ 20%, and<br /> Sequences were collected using the similar method de-<br /> ~ 17%). This indicates that the phenomenon, which has<br /> scribed by Tatusov et al. [91]. Protein-coding transcripts<br /> a bias to consume protein or oil mainly to power the life<br /> of Arabidopsis (TAIR release 10, https://www.arabidopsi-<br /> activities during seed development, does not occur in<br /> s.org/), rapeseed (release 4.1, http://www.genoscope.<br /> the evolution of oil crops.<br /> cns.fr/brassicanapus/), soybean (release Wm82.a2.v1,<br /> https://www.soybase.org/) and sesame (release 1.0, http:<br /> More evidence for candidate genes that are associated<br /> //ocri-genomics.org/Sinbase/) were downloaded, respect-<br /> with the seed oil content difference between soybean<br /> ively. If a gene has multiple transcripts, the longest was<br /> and rapeseed<br /> selected.<br /> Many candidate genes predicted in this study could be<br /> The transcriptome data of soybean (G. max Williams<br /> responsible for the difference of seed oil content be-<br /> 82) [58] (https://www.ncbi.nlm.nih.gov/geo/query/<br /> tween soybean and rapeseed have been experimentally<br /> acc.cgi?acc=GSE42871) and rapeseed (B. napus Darmor-<br /> confirmed to be related to seed oil content. In addition<br /> bzh) [57] (https://www.ncbi.nlm.nih.gov/geo/query/<br /> to those mentioned above, with the up-regulated expres-<br /> acc.cgi?acc=GSE77637. ) were downloaded from Gene<br /> sion for the genes BnGRF2 and BnRBCS1A and the<br /> Expression Omnibus (GEO). Rapeseed transcriptome<br /> down-regulated expression of gene BnPDK1, Hu et al.<br /> data included four seed developmental stages: 2, 4, 6,<br /> [69] cultivated a rapeseed line YN171 with a super high<br /> and 8 weeks after pollination (WAP), and soybean tran-<br /> seed oil content of 64.8% (Fig. 1). Similarly, with the<br /> scriptome data included seven seed developmental<br /> down-regulated expression of gene GmFAD2–1 by RNA<br /> stages: whole seed 4 days after fertilization (DAF), whole<br /> interference, seed oleic acid content in soybean in-<br /> seed 12–14 DAF, whole seed 22–24 DAF, whole seed 5–<br /> creased to 94.58% and the linoleic acid content de-<br /> 6 mg in weight, cotyledons 100–200 mg in weight, coty-<br /> creased to < 3% (Fig. 1) [15].<br /> ledon 400–500 mg in weight, and dry whole seed. Of<br /> Differentially expressed genes, associated with seed oil<br /> which, whole seed 12–14 DAF, whole seed 22–24 DAF,<br /> content and identified among cultivars with different<br /> cotyledon 400–500 mg in weight, and dry whole seed in<br /> seed oil content, also provide relevant evidence. Among<br /> soybean are almost respectively equal to 2, 4, 6, and 8<br /> the 33 differentially expressed genes identified in rape-<br /> DAF based on the definition of soybean vegetative and<br /> seed by Xu et al. [90], PDAT (Additional file 1: Table S4)<br /> reproductive growth [92], which are consistent with 2, 4,<br /> and OBO (Fig. 1) are also found in the present study.<br /> 6, and 8 WAP in rapeseed. Thus, we selected the four<br /> Among the 28 core enzymes involved in lipid synthesis<br /> stages of soybean and rapeseed seed development men-<br /> in soybean [55], 8 were also found in the present study<br /> tioned above for subsequent analysis. The genes expres-<br /> (Additional file 3: Table S9).<br /> sion level (RPKM: reads per kilobase per million<br /> In this study, we are focusing on the difference of total<br /> mapped reads) were normalized and quantified by the<br /> seed oil content between soybean and rapeseed. How-<br /> DESeq package in Bioconductor [93].<br /> ever, their other differences exist as well, i.e., seed oil<br /> composition, grown climatic environment, nitrogen<br /> fixation, and species characteristics, which likely affect Delimitation of orthologous genes<br /> the conclusion in this study. Identification of orthologous groups (OGs) in Arabidop-<br /> sis, soybean, rapeseed and sesame was conducted using<br /> Conclusion OrthoMCL software with default parameters [94]. Based<br /> In this study, we identified candidate genes and their on all-against-all BLASTP of non-redundant protein se-<br /> transcription factors and microRNAs to explain the quences, clusters were obtained according to reciprocal<br /> difference in seed oil content between soybean and rape- best similarity pairs between and within species, using<br /> seed. First, PEPC, along with its microRNAs and co-ex- OrthoMCL software implemented by the Markov clus-<br /> pression genes, affect the carbon source flow in seeds, tering algorithm (MCL; http://micans.org/mcl/) [95]. To<br /> which may lead to differences in seeds oil content. Then, obtain more accurate results, two other known methods,<br /> BCCP1 and β-CT and its transcription factors that are namely Proteinortho [96] and Inparanoid 8 [97], were<br /> characteristic of rapeseed may result in high seeds oil also used to determinate OGs of soybean and rapeseed.<br /> Zhang et al. BMC Plant Biology (2018) 18:328 Page 11 of 15<br /> <br /> <br /> <br /> <br /> Prediction of candidate genes related to carbon tools/meme) with the following parameters: - the width<br /> metabolism and lipid biosynthesis of a motif was between 6aa and 50aa, and the number of<br /> Acyl-lipid biosynthesis process is mainly involved in fatty motifs was no more than 20.<br /> acid synthesis and elongation from pyruvate, TAG syn-<br /> thesis, and oil-body storage. In Arabidopsis, 135 acyl-<br /> lipid biosynthesis-related genes were downloaded from Selective pressure and positive selection analyses<br /> ARALIP (http://aralip.plantbiology.msu.edu/) [9], and The amino acid sequences were aligned using MUSCLE<br /> 238 carbon metabolism-related genes were also obtained [99], alignment gaps were manually deleted, and then<br /> using the method of Troncoso-Ponce et al. [10] with a used for following calculations. The ratio (ω value) of<br /> slight modification. Such genes were used as query, nonsynonymous substitution rate (Ka) to synonymous<br /> along with OGs, to identify carbon metabolism- and substitution rate (Ks) of homologous gene pairs was<br /> lipid biosynthesis-related OGs. To further identify candi- computed with the maximum likelihood method of the<br /> date genes for the differences of seed oil content, gene branch model in Codeml from the PAML package<br /> expression pattern clustering and interspecific relative (version 4.9) [104].<br /> copy numbers variation analysis were carried out. Gene To test for the variation of the ω ratio among differ-<br /> expression clustering analysis was performed using Short ent branches in gene trees, a branch-specific model was<br /> Time-series Expression Miner (STEM, http:// used and conducted in Codeml. The branch-specific<br /> www.cs.cmu.edu/~jernst/stem/) [56] with the following model allows the ω ratio to vary among branches in the<br /> parameters: Log normalize a time series vector of gene phylogeny (model = 2, NSsites = 0), and it could be used<br /> expression values (v0, v1, v2, ..., vn) to (0, log2(v1/v0), to test whether there are different ω values on particu-<br /> log2(v2/v0), ⋯, log2(vn/v0)), Minimum Absolute Expres- lar lineages [105]; thus, this model can be compared<br /> sion Change 2, −p 0.05. with the one-ratio model (model = 0, NSsites = 0) that<br /> assumes a constant ω value across all branches using<br /> the likelihood ratio test (LRT). The datasets used in Ka/<br /> KEGG pathway enrichment analysis Ks ratio estimation were further used in the next posi-<br /> Kyoto Encyclopedia of Gene and Genome (KEGG) path- tive selection analysis of the branch-site model (BSM)<br /> way enrichment analysis was performed using the online using the Bayes empirical Bayes method described by<br /> tool KOBAS (version 2.0; http://kobas.cbi.pku.edu.cn/ Yang et al. [104].<br /> index.php) [59]. The P-values for each KEGG biological<br /> process was calculated by Fisher’s exact test [59]. To<br /> control the false discovery rate (FDR ≤ 0.05), the Benja- Transcription factor (TF)- and microRNA-targets analysis<br /> mini-Hochberg method was used to conduct multiple Soybean and rapeseed microRNAs were downloaded<br /> testing correction [98]. In addition, the small term cutoff from miRBase (release 21, http://www.mirbase.org/)<br /> value was set at 5. [106]. psRNATarget (http://plantgrn.noble.org/psRNA-<br /> Target) [107] was used to identify miRNA targets with<br /> default parameters except for the Expectation (e) and<br /> Phylogenetic analysis and motif analysis Max UPE, which were set at 3 and 25, respectively. The<br /> The full-length amino acid sequence alignments were transcription factors (TFs) and TF-target pairs were<br /> performed using MUSCLE [99] with default parameters downloaded directly from PlantTFDB 3.0 (http://<br /> and then phylogenetic tree reconstruction was con- planttfdb.cbi.pku.edu.cn/) [61]. To ascertain whether<br /> ducted with both Neighbor Joining (NJ) and Maximum miRNAs controls target-genes expression in seed devel-<br /> Likelihood (ML) approaches in MEGA 7.0 [100]. In the opment stages, bioinformatic analysis software miRDB<br /> NJ method, parameter setups were as follows: - model: (http://mirdb.org/miRDB/) [108] was used to prelimin-<br /> poisson correction; Bootstrap: 1000 replicates; and gap/ arily verify whether there is a putative binding site for<br /> missing data: pairwise deletion. To ensure the accurate- miRNAs in the 3′-UTR of target-genes mRNA.<br /> ness of ML tree, which is constructed to eliminate the<br /> long-branch attraction (LBA) caused by distant species,<br /> we also used maximum likelihood approaches with Gene co-expression network analysis<br /> PhyML v3.0 [101], and estimated the best-fitting models Pearson’s correlation coefficients (r) were calculated<br /> with the jModeltest software [102]. The phylogenetic using the ‘cor’ function of R package. The gene<br /> tree was displayed, annotated and managed using iTOL expression data (Reads Per Kilobase per Million<br /> (https://itol.embl.de/) [103]. Conserved functional motifs mapped reads: RPKM) was used to calculate the<br /> were identified using the program Multiple Em for Motif correlation coefficients between genes. The criteria for<br /> Elicitation [64] (MEME v4.11.2, http://meme-suite.org/ determining co-expressional genes were set at r ≥ 0.9 or<br /> Zhang et al. BMC Plant Biology (2018) 18:328 Page 12 of 15<br /> <br /> <br /> <br /> <br /> r ≤ − 0.9 and P-values ≤ 0.05 [65]. Graphical visualization Oleosin; OGs: Orthologous groups; PAP: Phosphatidate phosphatase;<br /> of the gene co-expression network was performed PDCT: Phosphatidylcholine:diacylglycerol cholinephosphotransferase;<br /> PDK1: Pyruvate dehydrogenase kinase isozyme 1; PEP: Phosphoenolpyruvate;<br /> using Cytoscape 3.4.0 (http://www.cytoscape.org/) [109]. PEPC: Phosphoenolpyruvate carboxylase; PGK: Phosphoglycerate kinase;<br /> Genes co-expressed with the target gene, meeting the fil- PK: Pyruvate kinase; PKp: Plastidial pyruvate kinase; PYR: Pyruvic acid;<br /> ter criteria, were further used to conduct KEGG pathway RBCS1A: Ribulose bisphosphate carboxylase small chain 1A; SAD: Stearyl-<br /> ACPdesaturase; STERO: Steroleosin; SUC: Sucrose transport protein;<br /> enrichment analysis using KOBAS 2.0 [59]. WAP: Weeks after pollination<br /> <br /> Acknowledgements<br /> Additional files Not applicable.<br /> <br /> Additional file 1: Table S1. 27,236 OGs of all the protein-coding genes Funding<br /> in Arabidopsis thaliana, Glycine max, Brassica napus and Sesamum indicum. This study was supported by the National Natural Science Foundation of<br /> Table S2. Comparisons of sequence similarity-based protein families China (31571268, 31871242), Huazhong Agricultural University Scientific &<br /> between Glycine max (gma) and Brassica napus (bna). Table S3. List Technological Self-innovation Foundation (Program No. 2014RC020), and<br /> of selected genes related to carbohydrate metabolism and lipid biosynthesis State Key Laboratory of Cotton Biology Open Fund (CB2017B01). Each of<br /> in Arabidopsis thaliana. Table S4. 230 candidate orthologous groups related the funding bodies granted the funds based on a research proposal. They<br /> to oil synthesis of seed for soybean, rapeseed and Arabidopsis. Table S5. had no influence over the experimental design, data analysis or interpretation,<br /> Candidate orthologous groups (OGs) for the difference of seed oil content or writing the manuscript.<br /> between rapeseed and soybean. Note: Genes with red color were negatively<br /> correlated with seed oil content and were obtained by cluster analysis of Availability of data and materials<br /> gene expression. The 44 OGs with black color differed in gene relative copy The datasets supporting the conclusions of this article are included within<br /> number between soybean and rapeseed. Table S6. Expressional contents the article and its additional files.<br /> and relative copy number of candidate genes associated with the difference<br /> of seed oil content in rapeseed and soybean. Note: The stages t1, t2, t3 and Authors’ contributions<br /> t4 were defined as R3, R4, R7 and R8 in soybean, and 2, 4, 6, and 8 weeks YZ conceived and supervised the study. ZZ carried out the experimental<br /> after pollination (WAP) in rapeseed, respectively. (XLSX 37341 kb) works and analyzed the data. ZZ, YZ and JMD wrote and approved the<br /> manuscript.<br /> Additional file 2: Figure S1. The expression patterns for the genes of<br /> 230 gene families related to oil biosynthesis. The expression clustering<br /> Ethics approval and consent to participate<br /> analysis of 2048 soybean and rapeseed genes in the 230 gene families<br /> Not applicable.<br /> was performed using Short Time-series Expression Miner (STEM, http://<br /> www.cs.cmu.edu/~jernst/stem/) [56]. Here, t1 represents the seed oil ini-<br /> Consent for publication<br /> tial synthesis stage; t2 to t3 represent the rapid accumulation period of<br /> Not applicable.<br /> seed oil biosynthesis; t4 represents the gradual decline stage after the<br /> seed oil accumulation content reaches the peak. In the end, all 2048<br /> Competing interests<br /> genes were clustered into 20 clusters. Figure S2. The expression profiles<br /> The authors declare that they have no competing interests.<br /> (A-D) of candidate genes related to oil biosynthesis. One down-regulated<br /> trend (profile 3) (A) and three up-regulated trends from t2 to t3 stages of<br /> seed oil biosynthesis (profile 13, 16 and 18, respectively) (B, C, D). Figure Publisher’s Note<br /> S3. Comparison of the expression patterns of the candidate genes be- Springer Nature remains neutral with regard to jurisdictional claims in published<br /> tween rapeseed and soybean. Note: t1-t4 and t1’-t4’ represent four seed maps and institutional affiliations.<br /> development stages in rapeseed and soybean, respectively. PKp-α and<br /> PKp-β denote Alpha (α) and Beta (β) subunits of PK in plastid, respect- Author details<br /> ively. ACCase contains homogeneous structure ACC2 and heterogeneous 1<br /> Crop Information Center, College of Plant Science and Technology,<br /> ACCase complex, which are composed of α-CT, β-CT, BC and BCCP. Fig- Huazhong Agricultural University, Wuhan 430070, China. 2Zhengzhou<br /> ure S4. Comparison of the expression patterns of genes encoding en- Research Base, State Key Laboratory of Cotton Biology, Zhengzhou<br /> zymes PEPC, PK and ACCase. t1’, t2’, t3’ and t4’ represent R3, R4, R7 and University, Zhengzhou 450000, China. 3School of Agriculture, Policy and<br /> R8 at soybean seed development stages, and t1, t2, t3 and t4 represent 2, Development, University of Reading, Reading RG6 6AS, UK.<br /> 4, 6 and 8 weeks after pollination (WAP) at rapeseed seed development<br /> stages, respectively. Figure S5. Transcriptional regulation of key candi- Received: 19 January 2018 Accepted: 20 November 2018<br /> date genes for the difference of seed oil content between rapeseed and<br /> soybean. Figure S6. Evolutionary rate of each branch of PEPC gene fam-<br /> ily. ω0 = 0.340 represents the evolutionary rate when the evolutionary References<br /> rate of each branch is assumed to be the same. (PDF 1291 kb) 1. Durrett TP, Benning C, Ohlrogge J. Plant triacylglycerols as feedstocks for<br /> Additional file 3: Table S7. LRT results for selective pressure branch the production of biofuels. Plant J. 2008;54(4):593–607.<br /> model (Model 0