The whole genome protein sequences in Arabidopsis and soybean were downloaded from TAIR database (¹ TAIR release 10) and SoyBase database (² version wm82.a2. V1; Schmutz et al., 2010), respectively, while the whole genome protein sequences in cowpea and kidney bean were download from Legumeinfo database (³ Lonardi et al., 2019; version gnm1.ann1).

Transcriptome datasets were downloaded from Phytozome database⁴ for soybean flower and embryo, Legumeinfo database (see text footnote 3) for cowpea and kidney bean flower and embryo, NCBI⁵ for soybean seed, and NCBI⁶ for cowpea seed.

Seven hundred fifty-six mature miRNA sequences in soybean were downloaded from miRBase database⁷, while 656 mature miRNA sequences in cowpea were downloaded from Martins et al. (2020).

The OrthoFinder was used to identify homologous gene families (Emms and Kelly, 2015) when we input the whole genome protein sequence for Arabidopsis, soybean, and cowpea, and adopted the default parameters. Based on the homologous gene families, the number of unique, single-copy, multi-copy, and unclustered gene families and protein numbers in soybean and cowpea were counted.

The phylogenetic tree was constructed by Neighbor–Joining (NJ) approach using the MEGA7 software⁸ and protein sequence alignment was carried out by the ClustalW method. The number of replicates in Bootstrap Test was set as 1,000, and the other parameters were set as default values. The final result tree file was visualized and beautified using the ITOL online tool⁹.

The online tools GSDS2.0 (¹⁰ Gene Structure Display Server), pfam¹¹, Wolf PSORT¹², and MEME¹³ were used to graphically display gene structure, search the target protein sequences of conservative domain structure, predict the subcellular localization of target protein, and analyze conservative motif, respectively.

The collinearity among Arabidopsis thaliana, soybean, cowpea, and the collinearity within soybean species, were analyzed using TBtools software (Chen et al., 2020). The sequences of target proteins in Arabidopsis were submitted to the STRING database¹⁴ to search for interaction proteins with experimental evidence, active interaction sources in soybean were set up as “all,” and combined score >0.4 was regarded as the cut-off point of significance. In cowpea, there is no interacted protein database available. Thus, interaction protein network in mung bean (Vigna radiata), which is the closest species of cowpea, was used to represent the one in cowpea. Here, comparative genomics analysis between cowpea and mung bean was conducted on Phytozome and OrthoFinder. The protein interaction network was beautified by the Cytoscape software. Gene function annotation was conducted on SoyBase and Phytozome for soybean and TAIR for Arabidopsis.

The expression patterns of genes in flower, leaf, stem, root, pod, and seed in soybean, cowpea, and kidney bean were analyzed via standardized expression levels, log₂(y + 1), where y was real expression level. Relative expression levels were defined as the ratio of the expression level of one gene to average expression level of all the genes in the species (deleting the genes with expression level less than 1.0) (Zhang et al., 2018; Zhou et al., 2019; Cheng et al., 2021) and used to compare the differences of expression levels of genes in the flowers of soybean and cowpea. The TBtools and SigmaPlot software packages were used to draw HeatMap and relative expression levels, respectively.

The online tool psRNATarget¹⁵, a plant small RNA target analysis server, was used to predict micoRNAs, implemented by Schema V2 (2017 release) where expected value was set as 4 and other parameters were set as their default values. The number of predicted miRNA for each gene was counted and then plotted by SigmaPlot.

In this study, we focused on two types of patterns in gene expansion: tandem and segmental duplications. The above collinearity results were used to determine duplicated gene pairs, and these gene pairs were compared with the gene pairs downloaded from PlantDGD¹⁶ to validate the predicted gene pairs. The formula T = Ks/2λ was used to calculate the date of occurrence of repeated events, where λ is equal to 6.1 × 10^–9 (Lynch and Conery, 2000).

To identify orthologous genes (OGs) in soybean and cowpea, all the genes in Arabidopsis, soybean, and cowpea were clustered using the OrthoFinder software. As a consequence, 113,233 protein-coding genes from the three species were clustered into 21,582 OGs (Supplementary Table 1), with each OG representing a gene family. Among these gene families, 1,541 (7.14%) were identified as soybean-specific gene families, and only 451 (2.09%) were identified as cowpea-specific gene families. The two proportions (9.23%) were very low, indicating the very high similarity in the evolutionary process between soybean and cowpea. Meanwhile, 549 (2.54%) OGs with single-copy soybean and multi-copy cowpea genes was significantly lower than 10,196 (47.24%) OGs with single-copy cowpea and multi-copy soybean genes. This means that soybean underwent a unique allotetraploidy event.

Among all the above OGs, five GmPIN1 genes in soybean and three VuPIN1 genes in cowpea were identified. The genes GmPIN1a, GmPIN1b, GmPIN1c, VuPIN1a, VuPIN1b, and A. thaliana PIN1 were found to be in one gene family, and GmPIN1b and VuPIN1b had the highest relative expression levels in the flowers of soybean and cowpea (Table 1) and the closest genetic distance to A. thaliana PIN1 gene (Figure 1A). All the PIN1 genes were predicted by the Wolf PSORT software to be localized in the plasma membrane.

In the conserved motif analysis via the MEME software, almost all the PIN1 proteins were structurally relatively conserved (Figure 1B), except for GmPIN1c and GmPIN1e, which may be the repeat proteins with the lack of the sequences at the 3′ end of GmPIN1b and at the 5′ end of GmPIN1d, respectively. In the comparison of relative expression levels, GmPIN1a and GmPIN1b had significantly lower relative expression levels than VuPIN1b (Figure 1C). In the collinearity analysis in Arabidopsis, soybean, and cowpea via the TBtools software (Chen et al., 2020), GmPIN1a was collinear with A. thaliana PIN1 gene; GmPIN1d and GmPIN1e were collinear with VuPIN1c; and GmPIN1a, GmPIN1b, and GmPIN1c were collinear with VuPIN1a and VuPIN1b. The results were consistent with those in the evolutionary tree (Figure 1D).

Ten proteins in A. thaliana were found in STRING database to be interacted with A. thaliana PIN1 based on experimental evidence. These proteins were used to mine homologous genes in soybean and cowpea. If one protein had multiple homologous genes in soybean or cowpea, the gene with the highest relative expression level was selected. In the large gene family of TOPP4; meanwhile, further phylogenetic analysis was conducted to determine the closest homologous relationship (Supplementary Figure 1A). The relative expression levels of homologous genes selected in soybean and cowpea were compared to find the differences of these homologous genes in soybean and cowpea. As a result, MPK6, PID, D6PK, UNH, and RCN1 had significant differences of relative expression levels in soybean and cowpea (Figure 1F), while ABCB1, ABCB19, DL1, PGP1, and TOPP4 had similar relative expression levels in soybean and cowpea (Figure 1E and Supplementary Table 1). Here we focused on the first five proteins.

First, we discussed their functions. PID and D6PK have been confirmed to positively regulate auxin transport through phosphorylation at the PIN1 serine S1 (S231), S2 (S252), S4 (S271), and S3 (S290), which plays an important role in ovular formation (Benjamins et al., 2001; Lee and Cho, 2006; Zourelidou et al., 2014). The S1 to S4 phosphate sites in A. thaliana were also found in the five soybean and cowpea PIN1 proteins, which are in the same OG as A. thaliana PIN1 (Supplementary Figure 2). RCN1 is found to inhibit the phosphorylation of PIN1 protein to antagonize the polar transport of auxin (Rashotte et al., 2001; Zhou et al., 2004; Michniewicz et al., 2007; Dai et al., 2012), while MPK6-mediated phosphorylation of PIN1 leads to the loss of the plasma membrane localization of PIN1, affecting auxin polar transport (Dory et al., 2018). Up to now, only one phosphate site of MPK6, S337, has been confirmed (Jia et al., 2016), and this site was found to be conserved in soybean and cowpea PIN1 proteins (Supplementary Figure 2). UNH is important in reducing PIN1 expression level in marginal cells, possibly through the localization of PIN1 into vacuoles for degradation (Pahari et al., 2014).

Then, we compared their relative expression levels. As a result, PID and D6PK, which are conducive to the phosphorylation of PIN1, had lower relative expression levels in soybean than in cowpea, while RCN1 and UNH, which affect plasma membrane localization, and MPK6, which affects the phosphorylation of PIN1, had higher relative expression levels in soybean than in cowpea. We speculated that an important reason for the SNPP difference between soybean and cowpea lies in the difference of relative expression level of PIN1. Owing to the difference of relative expression levels of genes encoding the interaction proteins of PIN1, PIN1 protein in soybean was less located on plasma membrane, and there was lower phosphorylation level in soybean. Thus, lower auxin transport efficiency results in less auxin maximum zone and fewer ovules and SNPP.

More importantly, the above method was used to predict the proteins that interact with PIN1, CKX3, and CKX5 homologies in soybean and cowpea. The results are shown in Figure 2A and Supplementary Table 7. In the proteins interacted with PIN1 homologies, two soybean proteins, GmABCB19 (Glyma.13G063700) and GmPID (Glyma.13G220100), were consistent with those in Arabidopsis, while six cowpea proteins, VuABCB1 (Vigun01g162000), VuABCB2 (Vigun07g072700), VuABCB19 (Vigun04g051400), VuMPK6 (Vigun03g181200), VuD6PK (Vigun06g148700), and VuPID (Vigun06g179800), were consistent with those in Arabidopsis. Based on the functions and annotations of these interacted genes, they commonly focus on auxin transport and serine phosphorylation. Interestingly, Glyma.02G186700 and Glyma.10G106900 with high expression levels in soybean ovary may be involved in the hydrolysis of serine, which is the active site of PIN1, because their homology AT2G41530 (AtSFGH) had been proven to encode a serine hydrolase in Arabidopsis (Cummins et al., 2006). CKX3 and CKX5 homologies in the two legumes commonly focus on embryo development and substance metabolism (Figure 2B). GmSLD5 (Glyma.08g194400, and Glyma.07g011200) and VuSLD5 (Vigun10g174500) are the homologies of AT5G49010 (SLD5), whose mutant can cause defective embryo development in Arabidopsis (Meinke, 2020). Another pair of genes, GmACX4 (Glyma.18g202800) and VuACX4 (Vigun10g081200), are the homologies of ACX4 (AT3G51840), which encodes a short-chain acyl-CoA oxidase. ACX4 is essential at early embryo development stages, and the acx3acx4 double mutants abort during the first embryo development phase (Rylott et al., 2003). VuACX4 (1.98) had much higher relative expression level than GmACX4 (0.51). This may lead to better embryo development and more seeds in cowpea than in soybean. Thus, we infer that the hindrance of ovular formation via serine hydrolysis in soybean and the promotion of embryo development via VuACX4 in cowpea may be responsible for the SNPP difference in the two legumes.

Using the database PlantDGD (Qiao et al., 2019), we obtained 13 pairs of duplicate genes in soybean, which were consistent with our collinearity analysis results. Among these duplicate genes, nine pairs were normal and four pairs were abnormal, namely, three pairs of GmCKX3a and GmCKX3b on chromosomes 9, 13, and 17 owing to the mismatch in collinearity analysis, and one pair of GmCKX6-2 and GmCKX1-2 (Supplementary Table 6). Here, we replaced GmCKX3b with its adjacent GmCKX3a and calculated the Ks values of these gene pairs to estimate their replication times (Supplementary Table 5). As a result, the duplications for six, six, and one pairs of GmCKX genes occurred, respectively, between 8 and 18 mya, between 70 and 74 mya, and between 133 and 169 mya. This indicates that most GmCKX duplications occurred at 10–15 mya [a soybean-specific whole genome duplication (WGD) event] and 59 mya (a legume-specific WGD event), and individual GmCKX duplications occurred at approximately 150 mya.

Using the database PlantDGD (Qiao et al., 2019), we checked the fragment duplications of CKX gene family in kidney bean and Arabidopsis. As a result, a pair of PvCKX1 and PvCKX6 was obtained, while three pairs of possible fragment replicators CKX2 and CKX3, CKX2 and CKX4, and CKX3 and CKX4 in Arabidopsis were observed. Thus, we speculated that CKX2 and CKX4 were a segmental duplication event that occurred in other species after legume differentiation, while one copy of CKX6 gene was a segmental duplication of CKX1 in legumes. The CKX6 copy from CKX1 is different from the other CKX6 copies in expression pattern. The former was not expressed in flowers, while the latter was expressed in flowers.

To confirm whether the above SNPP relationship and the difference of expression levels of PIN1 and CKX gene families exists in kidney bean, we compared their relative expression levels in soybean and kidney bean. As a result, Phvul.004G150600 had much higher relative expression level in kidney bean flowers than GmPIN1a and GmPIN1b in soybean flowers, while GmCKX gene family had higher relative expression levels in flowers than PvCKX gene family (Supplementary Figure 5). The results were consistent with those in soybean and cowpea. The differences of relative expression levels of PIN1 and CKX gene families in soybean and cowpea, along with the differences in soybean and kidney bean, may be an important reason for the SNPP difference.

Candidate SNPP genes in this and previous (Schwarz et al., 2020; Fang et al., 2021) studies and known gene Ln were used to construct interaction networks in soybean and cowpea. All the results are shown in Figure 5 and Supplementary Table 7. In soybean network, low SNPP may be due to two reasons. First, four-seed-pod-related gene Glyma.10G002200, in Fang et al. (2021), was interacted with GmCBP-1 and GmCBP-2, while its homology Vigun07g002900 in cowpea was interacted with VuCBP-1, VuCP1, and VuCAM4. The five interacted genes in soybean and cowpea belong to the calmodulin and calcium-binding protein gene (CBP) family. In previous studies, calmodulin and calcium-binding proteins in plants not only directly affected SNPP (Midhat et al., 2018), but were also critical for the biosynthesis of brassinosteroid (Du and Poovaiah, 2005), which plays an important role in determining the number of ovules and seeds via positive regulator BZR1 (Huang et al., 2013). However, VuCBP-1 (5.67), VuCP1 (4.30), and VuCAM4 (4.94) had much higher relative expression levels than GmCBP-1 (0.38) and GmCBP-2 (0.90), which are inhibited by gma-miR4405. Although Vigun03g412600 was inhibited by 016048_minus, its binding degree was relatively low (Supplementary Table 8). Thus, we speculated that the high expression may increase SNPP in cowpea.

Then, SNPP-related genes, GmCKX5-1 and GmCKX5-2, were interacted with soybean genes GmIPT-1 (Glyma.11G188100), GmIPT-2 (Glyma.13G271500), and GmIPT-3 (Glyma.12G086300). Their homologies in Arabidopsis all belong to Isopentenyltransferase (IPT) gene family, which plays an important role in cytokinin biosynthesis (Miyawaki et al., 2004, 2006), and overexpressing IPT can increase the levels of endogenous cytokinins (Décima Oneto et al., 2016). Their low expression levels in soybean may lead to a decrease in cytokinin synthesis, which further affects ovule numbers in soybean. The negative regulation was not found in cowpea.

In addition, a known SNPP gene Ln was interacted with Glyma.03G194700 (GmTFL1-1) and Glyma.19G194300 (GmTFL1-2), for which their Arabidopsis homology AT5G03840 (TFL1) determines seed size, and loss-of-function mutants exhibit a large seed phenotype (Zhang et al., 2020). Low expression levels of GmTFL1-1 (0.10) and GmTFL1-2 (0.01) may lead to large seed in soybean. The interacted network was not found in cowpea.

In this study, all the known yield-related soybean genes in Zhang et al. (2021), along with the above known and candidate SNPP genes, were used to construct a comprehensive network for the two legumes. The results are shown in Figure 6 and Supplementary Table 9.

In the network, some negative regulation was found in soybean. First, known seed size/weight gene Glyma.07G081700 was found to interact with Glyma.09G029400, Glyma.01G199800, and Glyma.15G135200 of cytochrome P450 subfamily 88A (CYP88A) gene family, which catalyzes the conversion of KA (ent-kaurenoic acid) to GA12 (the precursor of all gibberellins) and catalyzes three steps of gibberellin biosynthesis pathway in Arabidopsis (Helliwell et al., 2001; Regnault et al., 2014). It should be noted that GA was found to negatively regulate the number of ovules in Fuentes et al. (2012) and Carrera et al. (2012). Then, Glyma.07G081700 was found to interact with Glyma.U002000 that encodes gibberellin oxidase, whose expression decreases cytokinin activity (Jasinski et al., 2005). Finally, known seed-size-gene Glyma.05g019200 was found to interact with Glyma.11G067700 and Glyma.02G057500, which encode 22α hydroxylase that is an inhibitor of Brassinosteroid (BR) biosynthesis in Arabidopsis (Fujiyama et al., 2019). The negative regulation may result in lower SNPP in soybean than in cowpea.

Meanwhile, some positive regulation was found in cowpea. First, the homology Vigun01g142900 of known plant-architecture gene Glyma.19G164600 was found to be interacted with Vigun07g265600 and Vigun05g280100, while the homology Vigun03g358000 of known seed size gene Glyma.17G112800 was found to be interacted with Vigun05g297000. Vigun07g265600 (VuIAA34), Vigun05g297000 (VuTPL), and Vigun05g280100 (VuIAA18) focused on auxin synthesis and transport (Liscum and Reed, 2002). Then, the homology Vigun01g173000 of known plant-architecture gene Glyma.19G194300 was found to be interacted with Vigun05g004800 and Vigun09g178000 in 14-3-3 gene family, which is involved in PIN auxin carrier, auxin transport-related development, and brassinosteroid signal transduction (Gampala et al., 2007; Keicher et al., 2017). In this study, Vigun05g004800 and Vigun09g178000 were found to have high expression levels at ovule developmental stages. Third, the homology Vigun03g411300 of known seed-weight gene Glyma.17g036300 was found to be interacted with cell wall invertase genes Vigun01g175600 and Vigun07g062400, in which their homology AT2G36190 (CWIN4) may regulate ovule formation by modulating downstream auxin signaling and MADS-box transcription factors in Arabidopsis (Liao et al., 2020). Finally, the homology Vigun06g114600 of known seed size gene Glyma.17G112800 was found to be interacted with Vigun06g114600 and Vigun07g156800, in which their homology AT5G28640 (AN3) can regulate seed embryo development together with AT1G55600 (MINI3), and its mutant line had lower seed/silique number, silique length, and seed/silique weight than wild-type plants (Meng et al., 2016). The positive regulation results in higher SNPP in cowpea than in soybean.

Known seed size or weight genes Glyma.01G061100 (GmCYP78A70), Glyma.02G119600 (GmCYP78A57), Glyma.19G240800 (GmCYP78A72), and Glyma.05G019200 (GmCYP78A10) belong to the CYP78A gene family, which is found to be associated with seed size in Arabidopsis (Adamski et al., 2009; Fang et al., 2012). In cowpea, there were only two homologous copies Vigun02g047800 and Vigun03g343100 of the four soybean genes. We found that the first three soybean genes had high expression levels at middle and later seed development stages, while all the two cowpea copies had low expression levels (Supplementary Table 10). This may explain why soybean seed is larger than cowpea seed.

In summary, species-specific traits in crops may be derived from species-specific gene networks.

In this study, we observed four interesting results. The results are showed in Figure 7. First, PIN1 had lower expression level in soybean flowers than in cowpea flowers. Then, among the proteins, PID, D6PK, RCN1, UNH, and MPK6 in Arabidopsis, that interacted with PIN1 and their homologies were differentially expressed between soybean and cowpea, lower PID and D6PK and higher RCN1 expression levels in soybean resulted in lower phosphorylation level in soybean flowers, while high UNH and MPK6 expression levels decreased plasma membrane localization level in soybean flowers, as compared with those in cowpea flowers. Thus, we speculate that lower PIN1 expression level and lower phosphorylation and plasma membrane localization levels derived from the above five differential expression proteins interacted with PIN1 make auxin transport efficiency lower in soybean flowers than in cowpea flowers, which forms a lower auxin maximum zone (Benková et al., 2003; Ceccato et al., 2013). Next, some differential genes in the interaction networks on PIN1 and CKX gene families were found. SFGH, which hydrolyzes serine, was found in soybean rather than in cowpea, and ACX4 had much higher relative expression level in cowpea than in soybean. These may lead to better embryo development and more seeds in cowpea. Finally, in the interaction networks on yield-related genes, higher expression levels of the CBP genes in cowpea, as compared with those in soybean, increase SNPP; low expression levels of three soybean-specific IPT genes may inhibit CK synthesis. The above results may lead to lower number of ovules and ultimately lower SNPP in soybean than in cowpea.

In a previous study, one locus for soybean four-seed pod was located by Yang et al. (2013) on chromosome 13. Around this locus, there was one candidate gene GmMIF2 (Glyma.13G063100). Here we found that soybean gene GmMIF2 had much higher relative expression level in flower than cowpea gene VuMIF2 (Vigun04g052300). More importantly, their homologous gene MIF2 (AT3G28917) in Arabidopsis was found to hinder its ovule development (Bollier et al., 2018). Thus, MIF2 may be a potential gene for SNPP difference between soybean and cowpea.

In addition, we explained the reasons for low expression level of PIN1 gene in soybean. In this study we observed another two interesting phenomena. First, CKX genes and their positive regulator ROCK1 had higher relative expression levels in soybean flower than in cowpea flower. Then, fewer miRNAs were predicted to be targeted with soybean CKX genes than with cowpea CKX genes. These two results lead to higher cytokinin degradation level in soybean than in cowpea, which decreases cytokinin level in soybean (Supplementary Table 11) and affecting the expression level of PIN1 in soybean (Bencivenga et al., 2012; Zuñiga-Mayo et al., 2018).

In previous studies, three aspects of low auxin transport efficiency in Arabidopsis have been observed. First, PIN1 expression level affects auxin transport efficiency. Okada et al. (1991) compared the polar transport efficiency of auxin in the inflorescence tissue between pin1 mutants and their wild type, indicating the significant decreases of polar transport efficiency between pin1-1 (↓ 86%) and pin1-2 (↓ 93%) mutants and their wild type. Then, plasma membrane localization level of PIN1 affects auxin transport efficiency. Wisniewska et al. (2006) modified PIN1 polarity and examined auxin translocation direction in Arabidopsis thaliana. As a result, PIN1 polarity determines a primary direction in auxin transport of meristematic tissues. Finally, some kinases have been found to positively regulate auxin transport through PIN1 phosphorylation in Arabidopsis, such as PID and D6PK (Benjamins et al., 2001; Lee and Cho, 2006; Zourelidou et al., 2014). In addition, GmPIN1 is expressed and polarly localized in nodule primordium cells, and controls nodule formation by transporting auxin to form an auxin maximal zone in soybean (Gao et al., 2021). As we know, nodule primordium and ovule primordium are meristematic tissues. Thus, these findings provide evidence for the role of GmPIN1 in the formation of soybean ovule primordium. In this study, we found that the relative expression level of PIN1 gene was much lower in soybean flowers than in cowpea flowers, and the differences of relative expression levels of the above five proteins interacting with PIN1 might cause lower phosphorylation and plasma membrane localization levels in soybean than in cowpea. These results may cause lower auxin transport efficiency in soybean than in cowpea.

During ovule primordium formation, Benková et al. (2003) found that auxin accumulated in large amounts at the apex of ovule primordium, and this auxin maximum zone is a necessary condition for ovule formation. Bencivenga et al. (2012) showed that lower auxin transport efficiency significantly reduced ovule number per pistil from wild-type Col-0 (48 ovules) to weak mutant pin1-5 (9 ovules). Carlson and Lersten (2004) and Yang et al. (2017) observed that ovule number difference could cause SNPP difference. In summary, auxin transport efficiency difference in soybean and cowpea may cause lower formation efficiency of auxin maximum zone in soybean, resulting in lower ovule number and SNPP in soybean.

Cytokinin dehydrogenase (CKX) can specifically degrade cytokinin, although cytokinin can increase ovules number by promoting PIN1 expression (Bencivenga et al., 2012). Bartrina et al. (2011) compared single and double CKX gene mutations with their wild types in Arabidopsis thaliana. As a result, no significant change in the overall plant morphologies of single CKX gene mutants was observed, indicating the redundant role of CKX gene family. Multiple double mutations with CKX3-1 allele could form more flowers, especially for ckx3 ckx5 double mutant, which formed more ovules. This increased SNPP and led to 55% higher seed yield. Schwarz et al. (2020) obtained similar results in Brassica napus. In detail, compared with the wild type, the bnckx3 bnckx5 sixfold mutant increased the number of flowers, ovule number per pistil, and pod numbers on main stem, increasing seed yield by 20–32%. These results suggest that CKX gene family plays an important role in ovule formation. In this study, we identified 17 soybean and 10 cowpea CKX genes. Among these genes, most were not expressed in flowers, and multiple copy genes CKX3, CKX6, and CKX7 were much higher in soybean flowers than in cowpea flowers. These high expression CKX genes in soybean may lead to more cytokinin degradation and lower cytokinin content, decreasing PIN1 expression level, ovule number per pistil, and SNPP.

At present, there have been limited studies on CKX gene regulator, and only one regulator was reported by Niemann et al. (2015). In detail, the enhanced CKX activity in 35S:CKX1, 35S:CKX2, and 35S:CKX3 plants was reduced through rock1 introgression by 90.5, 64, and 100%, respectively. Meanwhile, rock1 mutant enhanced the activity of apical meristem and organ formation rate in Arabidopsis thaliana. Cytokinin content in inflorescence was increased and the numbers of flowers and pods on main stem were 50% higher in rock1 mutant than in its wild type, which was very similar to the phenotypic changes in the ckx3 ckx5 mutant. These results indicate that ROCK1 acts as a positive regulator of CKX protein. In this study, we found one ROCK1 homologous gene in soybean or cowpea, and ROCK1 had higher relative expression level in soybean flowers than in cowpea flowers. The higher ROCK1 gene expression level may increase the activity of soybean CKX protein. In addition, fewer miRNAs were predicted to target GmCKX genes. In other words, less miRNA regulation may be an important reason for higher GmCKX expression level in soybean. Thus, higher ROCK1 expression and fewer miRNAs enhanced cytokinin degradation by regulating GmCKX, so ovule number and SNPP in soybean decreased.

Molecular design breeding has been widely used in soybean with some success (Haun et al., 2014; Bao et al., 2019; Han et al., 2019; Le et al., 2020). As we know, cowpea has more SNPP and smaller seeds, and soybean has fewer SNPP and larger seeds. However, cowpea has significantly higher yield than soybean. We also notice that relatively few SNPP in wild and cultivated soybeans may hinder the increase of soybean yield. Therefore, it is possible and necessary to mine yield-related novel genes in cowpea to conduct molecular design breeding in soybean, especially in the current situation of very sharp contradiction between soybean supply and demand in China.

First, Gmckx3a quadruple mutant may be used to increase SNPP and yield in soybean. In previous studies, ckx3 ckx5 mutant in Arabidopsis and bnckx3 bnckx5 sixfold mutant in Brassica napus can increase their yields (Bartrina et al., 2011; Schwarz et al., 2020). In soybean, we found high expression of four CKX3a and no expression of two CKX5 and three CKX3 tandem-duplication-derived genes (CKX3b) in flowers (Figure 4). Thus, ckx3a quadruple mutant may be used to increase SNPP and soybean yield. In addition, GmCKX6-1 and GmCKX7-2 were found to be highly expressed in flowers in this study. Another possible way is to obtain gmckx3a gmckx6 (or gmckx7) mutant.

Second, we can identify elite alleles of GmCKX3a, GmCKX6 and GmCKX7 from existing four-seed pod cultivars, and transfer these elite alleles into current excellent cultivars via cross and backcross approaches. This method has been confirmed to be effective in Zhu et al. (2020).

Third, we may interfere with the expression of CYP88A genes (Glyma.09G029400, Glyma.01G199800, and Glyma.15G135200) and increase DELLA protein to decrease GA content in soybean. Carrera et al. (2012) reported that inhibiting the synthesis and function of GA via DELLA protein in tomato can increase ovule number.

Finally, we may over-express the CBP genes in soybean (Glyma.17G019400 and Glyma.07G259400) and transfer an excellent mutant gene GmBZL2(P216L) (GmBZL2*) into soybean (Zhang et al., 2016) to increase SNPP via enhancing the synthesis and signal transduction of BR.