The protein-coding sequences of all the genes in Arabidopsis, soybean, peanut, pea, and chickpea were downloaded from https://www.arabidopsis.org/index.jsp (TAIR10), https://www.soybase.org/ (Wm82.a2.v1), https://legumeinfo.org/data/public/Arachis_hypogaea (arahy.Tifrunner.CCJH), https://legumeinfo.org/data/public/Pisum_sativum (Cameor.gnm1. ann1.7SZR), and https://phytozome.jgi.doe.gov/pz/portal.html (Carietinum_492_v1.0), respectively. If there were multiple transcripts for a gene, the longest one was selected.

Transcriptome datasets in soybean, chickpea, peanut, and pea were downloaded from https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE42871 (GSE42871; Jones and Vodkin, 2013), https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE79719 and https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE79720 (GSE79719 & GSE79720; Garg et al., 2017), https://peanutbase.org/gene_expression/ atlas (Clevenger et al., 2016 and Alves-Carvalho et al., 2015), respectively. The soybean datasets included seven developmental stages: whole seed 4 days after fertilization (DAF) (R2); whole seed 12–14 DAF (R3); whole seed 22–24 DAF (R4); whole seed 5–6 mg in weight (R5); cotyledons 5–6 mg in weight, seed coats 5–6 mg in weight, cotyledons 100–200 mg in weight (R6); cotyledons 400–500 mg in weight (R7); and dry whole seed (R8). The chickpea datasets included S1–S7 seed development stages and leaf. The stages S1–S3 were the embryo development stages (embryogenesis), stages S4 and S5 were early and mid-maturation stages, respectively (grain filling), and stages S6 and S7 were late mature stages (seed desiccation). Based on the definition of soybean vegetative and reproductive growth (Ritchie et al., 1982), whole seed 12–14 DAF (R3), whole seed 22–24 DAF (R4), whole seed 5–6 mg in weight (R5), cotyledons 400–500 mg in weight (R7), and dry whole seed (R8) in soybean were almost consistent with the stages S2, S3, S4, S5, and S7 in chickpea, respectively. These stages were uniformly defined as the stages t1, t2, t3, t4, and t5 in this study, respectively. Thus, the five stages were selected to conduct expressional trend analysis. The peanut datasets included Pattee 5 seed, Pattee 6 seed, Pattee 7 seed, Pattee 8 seed, and Pattee 10 seed. The pea datasets included Stem_BC_LN, Pods_C_LN, seeds_5dai, seeds_5dai_mut, Seeds_12dap, Seeds_12dap_mut, RootSys_A_HN, RootSys_A_LN, Leaf_B_LN, and Flower_B_LN. The gene expression levels [reads per kilobase per million mapped (RPKM) reads] were normalized and quantified by the DESeq package in Bioconductor. The package DESeq provides methods to test for differential expression by estimating variance-mean dependence in count data from high-throughput sequencing assays and test for differential expression based on a model.

Zhang et al. (2018) described that the relative expression level was the ratio of the expression level of each gene to the average expression level of all genes in this species.

In this study, all the non-redundant whole-genome protein sequences in Arabidopsis, soybean, peanut, pea, and chickpea were clustered by the software package OrthoFinder 2.0 (Emms and Kelly, 2019) to identify orthologous gene families. All the genes of soybean and chickpea were analyzed by MCscan (Tang et al., 2008) to identify homologous collinearity gene fragments.

The expressional levels of soybean- and chickpea-specific genes at the above seven different seed development stages and tissues were analyzed by R package WGCNA v1.68 to construct co-expressional networks. When calculating the scale-free topological overlap matrix (TOM), the parameters “unsigned” and “bicor” were used. The optimal soft thresholds in soybean and chickpea were calculated by the function “pickSoftThreshold,” and their thresholds were 10. The minimum module size was set to 30, and the “merge cut height” parameter of the module was set to 0.3 (Langfelder and Horvath, 2008). The KEGG enrichment analysis for the genes in the above co-expressional networks was conducted by KOBAS (http://kobas.cbi.pku.edu.cn/index.php, version 2.0; Xie et al., 2011).

The expressional levels of soybean and chickpea common genes at the t1–t5 stages were used to perform cluster analysis, as described in Zhang et al. (2018), using Short Time-series Expression Miner (STEM, http://www.cs.cmu.edu/~jernst/stem/) (Ernst and Bar-Joseph, 2006) with the following parameters: log normalize a time-series vector of gene expression values (v₀, v₁, v₂,…, v_n) to [0, log₂(v₁/v₀), log₂(v₂/v₀), …, log₂(v_n/v₀)], minimum absolute expression change 2, –p 0.05.

One hundred and thirty-five genes related to fatty acid synthesis, TAG synthesis, lipid droplet assembly, and storage were downloaded from ARALIP (http://aralip.plantbiology.msu.edu/; Li-Beisson et al., 2013; Zhang et al., 2018). Thirty-seven starch metabolism-related genes in Arabidopsis were mined from Schwarte et al. (2015). Using the method of Troncoso-Ponce et al. (2011), 238 genes related to carbon metabolism were obtained. The copy number difference analysis was carried out for the above 410 genes (Supplementary Table 4). Meanwhile, the 410 genes were compared with the ones in categories three and four of the STEM analysis to mine common genes related to carbon, lipid, and starch metabolism.

The microRNAs in soybean and chickpea were downloaded from miRbase (21st edition, http://www.mirbase.org/) and Jain et al. (2014), respectively. The online software psRNATarget (http://plantgrn.noble.org/psRNATarget) (Dai and Zhao, 2011) was used to identify miRNA targets with default parameters except for the Expectation (e), which was set at 3. The TFs were downloaded from PlantTFDB 3.0 (http://planttfdb.cbi.pku.edu.cn/) (Jin et al., 2017). The expressional levels of 88 gene families obtained from copy number analysis and the target TFs of microRNA were used to calculate the Pearson correlation coefficient.

The gene families of WRI1 in soybean, chickpea, Arabidopsis, peanut, and pea were used to construct a gene phylogenetic tree using the neighbor-joining (NJ) method, implemented by software MEGA7 (http://www.megasoftware.net/, Kumar et al., 2016); the bootstrap test was repeated 1,000 times, using Poisson distribution model and “pairwise deletion” option. The iTOL (http://itol.embl.de/) was used to visualize and post-beautify the final result tree file.

The online tool MEME (http://meme-suite.org/tools/meme, v4.11.2) was used to predict the conserved domains of key protein-coding genes (Bailey and Elkan, 1994). The Clustalw program in the MEGA7 software was used to perform a multiple sequence alignment of amino acid sequences of the WRI1 homologous gene family in Arabidopsis, soybean, chickpea, peanut, and pea with default parameters. The GeneDoc software (https://github.com/karlnicholas/GeneDoc, Nicholas et al., 1997) was used to compare and beautify the results.

The collinearity analysis tool MCscan was used to analyze the collinearity of soybean and chickpea genomes. As a result, there were 31,519 and 1,443 repeated gene pairs in the soybean and chickpea genomes, respectively. The latter was far fewer than the former, which may be due to the soybean-specific tetraploid (Wang et al., 2017). In soybean, 472, 305, and 152 blocks had more than 10, 20, and 50 collinear genes, which contained 29,101 (52.8%), 26,689 (48.4%), and 21,172 (39.5%) collinear genes, respectively. In chickpea, 45, 21, and 5 blocks had more than 10, 20, and 50 collinear genes, which contained 1,214 (4.3%), 830 (2.9%), and 380 (1.3%) collinear genes, respectively (Supplementary Table 2).

All the 223,609 protein-coding genes in soybean, chickpea, Arabidopsis, peanut, and pea were clustered by OrthoFinder to identify homologous gene families in soybean and chickpea. As a result, these genes were clustered into 26,575 orthologous groups (OGs) (Supplementary Table 1), and each OG represents a gene family. Among these OGs, there were 335 OGs with both soybean single-copy and chickpea multi-copy genes, 8,675 OGs with both soybean multi-copy and chickpea single-copy genes, 3,261 OGs with soybean and chickpea multi-copy genes, and 2,396 OGs with soybean and chickpea single-copy genes, which are regarded as common genes, while 18,657 (33.29%) were defined as soybean-specific paralogous gene clusters without chickpea genes and only 8,010 (28.33%) were chickpea-specific genes (Supplementary Table 3; Zhang et al., 2018).

The expressional levels of 18,657 soybean-specific genes at seven different seed development stages and tissues (Jones and Vodkin, 2013) were analyzed by R package WGCNA v1.68 to obtain different co-expression modules. As a result, seven co-expression modules with different expression patterns were identified. The genes in each module were used to conduct KEGG enrichment analysis. Results showed that four modules, blue (3,318 genes), black (900 genes), green (563 genes), and light green (250 genes), significantly enriched multiple lipid metabolism pathways (P < 0.05; Figure 1B). The genes in the blue module were highly expressed between R2 and R5 seed development stages, and significantly enriched not only the carbohydrate metabolism pathways of pyruvate metabolism (gmx00620), pentose phosphate pathway (gmx00030), carbon metabolism (gmx01200), and glycolysis (gmx00010), but also the glycerolipid metabolism pathway (gmx00561), glycerophospholipid metabolism (gmx00564), linoleic acid metabolism (gmx00591), and other metabolic pathways. Meanwhile, these genes enriched the photosynthesis-antenna protein (gmx00196) and carbon fixation (gmx00710) metabolic pathways in photosynthetic organisms. Therefore, the blue module is considered to be a co-expression module for oil accumulation during soybean seed storage material accumulation (Figure 1A).

The expressional levels of 8,010 chickpea-specific genes at seven different seed development stages and tissues (Garg et al., 2017) were also analyzed. As a result, nine co-expression modules with different expression patterns were identified. The genes in the black (573 genes), tan (190 genes), yellow (335 genes), and dark red (101 genes) modules were highly expressed at S3–S5 stages, and they enriched cell replication and synthesis-related metabolic pathways. This may be related to the rapid expansion of the number of seed cells (Figures 1C,D).

To obtain the expressional trends of soybean and chickpea common genes, the standardized expression levels of 37,387 soybean and 20,259 chickpea common genes at different seed development stages (Supplementary Figure 1) were used to conduct cluster analysis with the k-means algorithm via the STEM software. As a result, these common genes were clustered into six clusters. Among these clusters, the expression of all the genes in Cluster3–Cluster5 was upregulated during the t2–t4 stages, while the expression of all the genes in Cluster1 was downregulated during the t2–t4 stages.

All the genes in each cluster were used to conduct KEGG metabolic pathway enrichment analysis. As a result, all the soybean genes in Cluster3 and Cluster4 enriched pentose phosphate pathway (gmx00030), fructose and mannose metabolism (gmx00051), galactose metabolism (gmx00052), fatty acid biosynthesis (gmx00061), fatty acid degradation (gmx00071), starch and sucrose metabolism (gmx00500), pyruvate metabolism (gmx00620), glyoxylic acid and dicarboxylic acid metabolism (gmx00630), carbon fixation in photosynthetic organisms (gmx00710), fatty acid metabolism (gmx01212), and other metabolic pathways (Supplementary Figure 2A). Meanwhile, all the chickpea genes in Cluster4 mainly enriched fatty acid biosynthesis (cam00061), fatty acid metabolism (cam01212), glycolysis/gluconeogenesis (cam00010), citrate cycle (TCA cycle) (cam00020), carbon fixation in photosynthetic organisms (cam00710), glycerophospholipids metabolism (cam00564), starch and sucrose metabolism (cam00500), glyoxylate and dicarboxylate metabolism (cam00630), and other pathways (Supplementary Figure 2B). Clearly, the above enriched metabolic pathways were found to be related to seed lipid and starch synthesis, and the expression of these genes was upregulated during the rapid accumulation stages of these substances, indicating the existence of fatty acids and starch in soybean and chickpea.

Four hundred and ten carbon-, lipid-, and starch metabolism-related genes in Arabidopsis thaliana (Supplementary Table 4) were used as the seed to search for gene families in 26,575 OGs. As a result, 374 common gene families were mined. Among these gene families, six had a >2 copy number ratio of soybean to chickpea, while eight had a copy number ratio of <1. The former includes GRF2 during photosynthesis, CT-α and Acyl carrier protein (ACP) involved in the fatty acid synthesis, and CDS, PECT, and DGAT participating in TAG synthesis, while the latter includes STP7 involved in carbohydrate transport, PKp-β involved in glycolysis, Beta-hydroxyacyl-acyl carrier protein dehydratase (HAD) and β-Ketoacy-acyl carrier protein synthase II (KAS II) participating in the fatty acid synthesis, glycerol-3-phosphate acyltransferase (GPAT9), LPAAT3, and FAD2 participating in TAG synthesis, and Pyrophosphorylase (PPA) involved in starch metabolism (Supplementary Figure 3). Meanwhile, the above 410 genes were traversed and the genes in Cluster3 and Cluster4 were searched to obtain homologous gene families. As a result, 158 soybean genes and 121 chickpea genes were found to be homologous with the above 410 genes.

By merging the above 279 upregulated genes and the above 14 gene families with copy number difference, 184 soybean and 135 chickpea genes in 97 OGs were obtained (Supplementary Table 5), which are candidate oil- and starch biosynthesis-related common genes in soybean and chickpea, and some genes have been confirmed by molecular biological experiments (Table 1).

To compare the expressional profile difference of each gene in the above 97 OGs at t1–t5 stages in soybean and chickpea, the relative expression level was used in this study (Supplementary Table 6). In soybean, most starch synthesis-related genes were highly expressed at the early seed development stage, while most oil synthesis-related genes were highly expressed in the middle seed development stage, which is also the rapid oil accumulation stage. In chickpea, most starch synthesis-related genes were highly expressed at the middle seed development stage.

Based on the studies in Klaus et al. (2004), Tang et al. (2016), and Andriotis et al. (2010), common genes ACCase, AGPase, and BAM5 were important genes in seed oil and starch metabolism. Thus, we compared the expressional profile differences of these genes in soybean and chickpea. First, the relative expression levels of genes encoding subunits α-CT, BCCP, and BC of ACCase at rapid oil accumulation stages were higher in soybean than in the chickpea (Figures 2A–C). Because ACCase catalyzes the first step of the fatty acid synthesis, which is the formation of malonyl-CoA from acetyl-CoA (Klaus et al., 2004), ACCase may make more carbon metabolites flow into the fatty acid synthesis in soybean seeds. Then, the relative expression levels of genes encoding subunits ADP-glucose pyrophosphorylase (APS) and glucose-1-phosphate adenylyltransferase (APL) of AGPase at the t4 stage were much higher in chickpea than in soybean (Figure 2D).

Because AGPase is the first key regulatory enzyme and rate-limiting enzyme in starch synthesis (Tang et al., 2016), AGPase may play a significant role in chickpea starch synthesis. Finally, the relative expression level of BAM5 at the early seed development stage was much higher in soybean than in chickpea (Figures 2E,F). Because BAM5 catalyzes the production of maltose from linear glucans (Andriotis et al., 2010), it may play a major role in plant starch hydrolysis. This may result in low starch content in soybean seeds. Therefore, we infer that ACCase, AGPase, and BAM5 were more likely to be candidate key genes that affect the differences of seed oil and starch content in soybean and chickpea.

Using the starch biosynthesis genes in rice (17 genes; Ohdan et al., 2005; Tian et al., 2009) and Arabidopsis (13 genes; Schwarte et al., 2015) as bait, we performed an ortholog search in the soybean and chickpea genomes. In total, fewer starch biosynthesis genes were found in chickpea than in soybean (33 vs. 45; Supplementary Table 7), but there was no significant difference between these two genomes in the exact Fisher test (P = 1.00). The relative expression levels of these starch biosynthesis genes at seven seed development stages were used to compare their differences. We found that the average relative expression level in chickpea was significantly higher than that of soybean, e.g., Ca_03834, Ca_04774, Ca_07632, and Ca_00773 have a higher expressional level at stage t3 (Figure 3), although more starch biosynthesis genes were detected in the soybean genome (Supplementary Table 7). The above results are similar to those in Yang et al. (2015).

As described above, the soybean-specific genes in four of seven modules were highly expressed at the middle and late oil accumulation stages and significantly enriched in fatty acid, glycerolipid, and linoleic acid metabolisms. Among these genes, main lipid synthesis-related genes include Glyma.19G147400 (Lakhssassi et al., 2017), Glyma.07G207200 (Deng et al., 2020), Glyma.12G146700 (Carther et al., 2019), and Glyma.06G058500 (Yang et al., 2021), which have been confirmed by molecular biology experiments. In detail, Glyma.19G147400 is responsible for the conversion of oleic acid to linoleic acid, Glyma.07G207200 catalyzes the synthesis of monounsaturated oleic acid or palmitoleic acid in plastids, Glyma.12G146700 transforms the diacylglycerol into phosphatidic acid, and Glyma.06G058500 is involved in long-chain fatty acids synthesis. Meanwhile, around stable loci for soybean seed oil content in Yao et al. (2020), Glyma.05G226800, Glyma.05G227400, Glyma.08G102100, Glyma.10G063600, Glyma.11G148900, Glyma.11G190600, Glyma.13G112700, Glyma.13G148100, Glyma.14G150100, Glyma.15G051100, and Glyma.20G088000 were found to be soybean-specific genes in the above four modules. Thus, the above soybean-specific and highly expressed genes may result in high seed oil content.

Three soybean copies (Glyma.13G057400, Glyma.18G265300, and Glyma.19G028800) of common gene BCCP2 in soybean and chickpea were highly expressed during t3–t4 stages; in particular, Glyma.19G028800 was rapidly upregulated and reached a peak during t3–t4 stages. As described in Konishi and Sasaki (1994), BCCP2 is an important subunit of ACCase. Overexpression of ACCase can significantly increase fatty acid content in crop seed (Roesler et al., 1997; Klaus et al., 2004), while antisense expression of BCCP2 transcript in developing Arabidopsis seeds resulted in an average 38% reduction in BCCP2 protein, further leading to an average 9% reduction in fatty acid content in mature seeds (Thelen and Ohlrogge, 2002). It should be noted that Glyma.13G057400 and Glyma.19G028800 (BCCP2) in soybean had significantly higher expression than Ca_21112 and Ca_10464 (BCCP2) in chickpea. Thus, the high expression genes Glyma.13G057400 and Glyma.19G02880 in soybean may result in its high seed oil content.

In this study, WRI1 was found to be co-expressed with oil synthesis-related genes ACCase, KAS, KAR, ACP, and LACS in soybean and peanut. This may result in high seed oil content in soybean. The evidence is as follows. Some target genes of AtWRI1 are involved in Arabidopsis glycolysis and fatty acid synthesis, while their mRNAs are accumulated in a synergistic manner and regulated by WRI1 (Baud et al., 2009). Twenty-eight transcripts, such as KASIII, KAR, and LACS, were found to have increased expression levels in the transgenic plants of GmWRI1, which can specifically bind the 500 bp upstream sequences of ATG codon of the above genes and effectively activate their transcription (Chen et al., 2018).

As described in Scheidig et al. (2002), BAM degrades starch and amylopectin to form maltose. In the study of Andriotis et al. (2010), seed starch content was 1.5–2.5 times higher in Arabidopsis thaliana bam4 and bam1234 mutants than in wild-type embryos, and the starch still existed in mature embryos, while seed starch content in amy3, isa3, and phs1 mutants did not increase, indicating the necessity of BAM in embryo starch degradation. In this study, two BAM5 copies, Glyma.06G301500 and Glyma.12G102900, in soybean had more than 6 times expression levels at the early seed development stage than one chickpea BAM5 gene Ca_22584. This indicates that seed starch, produced at the early seed development stage, maybe degraded by Glyma.06G301500 and Glyma.12G102900. This results in low seed starch content in soybean. NAC36-, bZIP-, and DOF-homologous genes in chickpea have more co-expressed starch biosynthesis genes than those in soybean. This may result in higher seed starch content in chickpea than in soybean.

The above results are summarized in Figures 4, 6.

Although CaWRI1 in chickpea has a similar expression pattern as its homologous genes Glyma.08G227700 and Glyma.5G221600 in soybean, no chickpea oil synthesis-related genes were co-expressed with CaWRI1. The possible reason is 24 extra amino acids of the RAYD element of the first AP2 domain of CaWRI1 before the VYL motif. In previous studies, mutant seeds with T-DNA inserted in the 5th exon of AtWRI1 decreased seed oil content 72% in Arabidopsis than in its wild type (An et al., 2017), while wri1 mutant could not convert glucose and sucrose into precursors of fatty acid biosynthesis during seed development, resulting in the decrease of seed oil content (Focks and Benning, 1998). Thus, it is speculated that the structural variation of CaWRI1 leads to the failure to regulate glycolysis and lipid synthesis-related genes, which may be an important reason for low seed oil content in chickpea.

As described above, AGPase had more than 7 times relative expressional level in chickpea than in soybean. Its high expression is likely to catalyze more Glc-1-P to ADP glucose (Crevillén et al., 2005). As evidenced in rice (Tuncel et al., 2014) and potato (Hajirezaei et al., 2003), AGPase could effectively control starch flux. In molecular biology, the 100-grain weight of maize seeds with overexpressing AGPase was increased by 15%, and seed starch content was increased by 74% compared to its wild type (Li et al., 2011), while the weight of wheat seeds with overexpressing cytoplasmic AGPase large subunit gene was increased by 9.1%, and seed starch content was increased by 9.6% compared to its wild type (Kang et al., 2013). Thus, high expression of AGPase is more likely to result in higher seed starch content in chickpea.

The above results are summarized in Figure 6.

Yang et al. (2015) performed the comparative genomic and transcriptome analyses of seed starch and oil content between adzuki bean and soybean in the genome sequencing project of adzuki bean (Vigna angularis). In this project, they found no significant variation of starch biosynthesis genes in the two genomes, meanwhile, they checked the expression levels of five key genes that encode the proteins involved in the conversion from starch to oil synthesis. Although Yang et al. (2015) did not compare starch synthesis-related genes, key TFs, miRNAs, and their regulation networks, oil- and starch synthesis-related genes, TFs, miRNAs, and their regulatory networks in soybean and chickpea were investigated in this study to obtain the molecular mechanisms for the differences in seed oil and starch content.

As pointed out by Dragičevi et al. (2015), increasing seed oil content in chickpea can reduce seed phytic acid content and improve its digestibility. However, there is low seed oil content in the present chickpea cultivars. To address this issue, it is necessary to improve seed oil content in chickpea. Based on the above results, two genetic improvement methods are proposed.

First, deleting the 24 extra amino acids of the RAYD element in Ca_15711 (CaWRI1) may make the edited CaWRI1 co-express with oil synthesis-related genes to increase seed oil content in chickpea. In this study, we found that the relative expressional level of Ca_15711 (CaWRI1) at stage S3 in chickpea is higher than those of soybean genes Glyma.15G221600 and Glyma.08G227700, homologous to Ca_15711, at stage R4 (Figure 5B). Shen et al. (2010) had provided transgenic evidence. In detail, seed oil content in maize was significantly higher in transgenic ZmWRI1 endosperm (0.81%) than in null endosperm (0.47%), but seed starch content was reduced approximately 60%, although the expression of ZmWRI1 does not affect protein content in embryo and grain yield, indicating that ZmWRI1 may enhance oil biosynthesis by reducing carbon metabolites in starch biosynthesis.

Second, increasing the expression of miRNA159 and miRNA319 in chickpea may decrease the expression of MYB33, which downregulates starch synthesis-related enzymes AGPase, SS, and ISA, making more carbon metabolites flow into fatty acid synthesis. Hu et al. (2021) provided experimental evidence. In detail, Zma-miR159k-3p could negatively regulate the expression of ZmMYB138 and ZmMYB115, homologous to AtMYB33, in maize endosperm, while ZmMYB138 and ZmMYB115 may be positively correlated with amylopectin and amylose content in maize endosperm, respectively. Thus, Zma-miR159k-3p may indirectly regulate starch synthesis in maize endosperm by regulating the expression of ZmMYBs.