We retrieved the genome assemblies and protein annotations for 52 Fabaceae species from multiple public databases, including NCBI, plantGIR (Liu et al., 2024) and LGRPv2 (Yu et al., 2025). Detailed information was provided in Supplementary Table S1. To establish a robust species evolutionary framework, a consensus topology (branching structure) was obtained from the TimeTree 5 web server (Kumar et al., 2022) and LGRPv2 (Yu et al., 2025). Documented WGD and whole genome triplications (WGT) events were manually annotated onto the species tree based on published literature (Jiao et al., 2012; Yu et al., 2025; Zhang et al., 2020; Zhao et al., 2021b).

The HMM profile of the conserved SRF-TF domain (PF00319) was used to search against the proteomes of all 52 Fabaceae species. To provide a broader evolutionary context, we also searched the proteomes of selected outgroup species: Polygala tenuifolia, Vitis vinifera, Amborella trichopoda, Nymphaea colorata, Arabidopsis thaliana, and Oryza sativa (Supplementary Table S1). All searches were conducted using HMMER (v3.4) with a cutoff E-value cutoff of 1×10^-5 (Potter et al., 2018). To ensure high accuracy, all candidate sequences were manually curated. Sequences shorter than 100 amino acids were removed, and the presence of the MADS domain in the remaining candidates was verified using the NCBI Conserved Domain Database (NCBI-CDD) web server.

To elucidate the evolutionary relationships within the MADS-box family, we performed a multi-step phylogenetic analysis. The full-length amino acid sequences of all identified MADS-box proteins were aligned using MAFFT (v7.505) (Katoh and Standley, 2013) and then trimmed to remove poorly aligned regions using trimAl (v1.4) (Capella-Gutiérrez et al., 2009). A maximum likelihood (ML) phylogenetic tree was constructed from the trimmed alignment using FastTree (v2.1.11) with the JTT (Jones-Taylor-Thornton) protein substitution model. Local branch support was assessed using the Shimodaira-Hasegawa-like (SH-like) method, based on 1,000 resamples, as implemented in FastTree (Price et al., 2009). Based on the topology of the phylogenetic tree and established domain architectures from foundational studies (Kaufmann et al., 2005; Parenicová et al., 2003). Tree visualization was done with Evolview (Subramanian et al., 2019).

Gene structure statistics, including gene length, amino acid length, intron length, intron/exon number were retrieved from genomic annotation GFF3 files for each species. To investigate the conservation and divergence of protein architecture, we analyzed the conserved motifs within each MADS-box subfamily using the MEME suite (v5.5.8) (Bailey et al., 2015). The analysis was set to identify a maximum of 5 motifs and only those that appear in over 50% of all sequences were selected for visualization.

The protein sequences were aligned using Diamond v2.0.5.143 in blastp mode with a E-value cutoff 0.001 and max-target-seqs 5 (Buchfink et al., 2021). Synteny analysis was conducted for each species using MCScanX with the following parameters: a minimum of colinear genes were required to define a block (-s 5), and the maximum gene gap allowed between collinear genes was set to (-m 25) (Wang et al., 2012). Based on the synteny information and genomic locations, duplicated gene pairs were classified into five types: WGD, tandem, proximal, dispersed, and transposed, using the DupGen-finder (v1.0.0) (Qiao et al., 2019). The rates of synonymous (Ks) and non-synonymous (Ka) substitutions were calculated using KaKs_Calculator (v2.0). To mitigate the impact of noise and artifacts on evolutionary rate estimates, only gene pairs with Ks values < 2 and P-Value < 0.05 were included in subsequent analyses. The ratio Ka/Ks was then used to infer the selection pressure acting on these duplicated genes.

To assess the conservation and dispensability of MADS-box subfamilies across the Fabaceae, we conducted a pangenome analysis. For each subfamily, protein sequences from all 52 species were clustered into orthologous groups (OGs) using CD-HIT (v4.8.1) a 40% sequence identity threshold (Fu et al., 2012). The resulting clusters were categorized into four frequency-based groups according to established pangenome nomenclature (Tong et al., 2025): (i) core (present in all 52 Fabaceae species), (ii) soft-core (present in >90% of species), (iii) shell (present in 10%–90% of species), (iv) cloud (present in <10% of species).

To investigate the functional divergence of MADS-box genes, particularly among duplicated pairs, we analyzed a large-scale RNA-seq dataset for soybean (Glycine max). We retrieved the expression profiles of 3,638 RNA-seq samples of soybean derived from various tissues and abiotic and biotic stresses from the soybean RNA-seq Database (https://plantrnadb.com/soybean/) (Yu et al., 2022) (Supplementary Table S2). FPKM values were extracted for all identified soybean MADS-box genes. These values were used to analyze tissue-specific expression patterns and to compare the expression divergence between duplicated gene pairs derived from WGD and SSD events.

To identify organ-specific expression modules, hierarchical clustering was performed using the pheatmap package in R. Prior to clustering, FPKM values were log2-transformed and standardized using z-scores to highlight expression trends rather than absolute abundance. A distance matrix was computed using the Euclidean method, and clustering was conducted using the Ward’s minimum variance method (Ward.D2). The resulting dendrogram was cut to define four distinct expression modules.

To validate the expression patterns, soybean seeds were germinated on moist paper for 5 days and then transferred into a light chamber under a 16h:8h, light: dark photoperiod with 60% relative humidity at 25 °C. The seedlings were grown in pots containing a 2:1 mixture of forest soil: vermiculite for 2 weeks. Select plants with uniform growth and divide them into four groups. For drought treatment, the seedlings roots were rinsed with water to remove adherent soil without damaging root hairs. Then excess surface water was blotted using filter papers for the induction of rapid drought for 0 (CK), 1, 3, 8 and 12 hours. For NaCl treatment, the seedlings were treated separately with 1/2 Hoagland nutrient solution, 150 mM NaCl solution and 200 mM NaCl solution on different periods (0, 1, 3, 8 and 12 hours), respectively. Leaf samples are collected into test tubes, rapidly frozen with liquid nitrogen, and subsequently used for RNA extraction. Total RNA was extracted using the Ultrapure RNA Kit (CWBIO, China), which includes TRIzon Reagent as a key component, and cDNA was synthetized with M-MLV reverse transcriptase kit (Takara). Quantitative real-time polymerase chain reaction (qRT-PCR) was performed with SYBR Premix Ex (TaKaRa) on the Roche Light Cycler 480 system (Roche, Germany) with PCR kit (Roche, Germany). The soybean endogenous gene TUBULIN (Glyma.05G157300) was used as an internal control, and the relative expression levels of examined genes were calculated using the 2^-ΔΔCt method (Livak and Schmittgen, 2001). Three biological replications were performed in each test. The primers are listed in Supplementary Table S3.

To systematically map the MADS-box gene family across the Fabaceae, we conducted genome-wide search in 52 species of legumes that represent three major subfamilies: 38 species from the Papilionoideae subfamily, 11 from the Caesalpinioideae, and 3 from the Cercidoideae (Supplementary Table S4). A total of 4,872 high-confidence MADS-box genes were identified, comprising 2,032 Type I and 2,840 Type II (MIKC) members (Supplementary Table S4). To elucidate the evolutionary relationships more clearly, we included several evolutionarily significant outgroup species in our analysis, encompassing, including Polygala tenuifolia, Vitis vinifera, Amborella trichopoda, Nymphaea colorata, Arabidopsis thaliana, Oryza sativa, and conducted identification for these species as well (Supplementary Tables S1, S4). The resulting phylogeny clearly partitions the family into 16 subfamilies (Figure 1). The Type I lineage is the largest branch, comprising the Mα, Mβ, and Mγ clades, with Mβ and Mγ showing a closer relationship. The Type II lineage consists of 13 distinct clades, with clear clustering observed among core floral development genes. Among them, AG/STK, SOC1, FLC, AGL6, SEP, and AP1/FUL are clustered together, suggesting a closer phylogenetic relationship. Additionally, the AP3/PI, BS, and AGL12 are also clustered together, indicating a closer relationship. The other clades in Type II are relatively independent, indicating that they may have undergone independent evolution.

Beyond their phylogenetic placement, a profound structural dichotomy between Type I and Type II genes became immediately apparent, providing a physical basis for their divergent evolutionary fates.

We observed significant variation in gene architecture across both lineages and subfamilies (Figure 2). A comprehensive structural survey of the MADS-box gene family among 52 legume species uncovered evolution patterns among different MADS-box subfamilies in three subfamilies (Papilionoideae, Caesalpinioideae and Cercidoideae) (Figures 2A–E). Gene length exhibits significant variation across different lineages, with Caesalpinioideae possessing longer gene (Figure 2A) and intron lengths (Figure 2C). This suggests that differences in intron length among lineages may contribute to variations in gene length. In contrast, amino acid sequence length shows no significant differences among the three lineages (Figure 2B), highlighting the high conservation of coding regions within the Fabaceae. Regarding intron number (Figure 2D) and CDS number (Figure 2E), Cercidoideae consistently have slightly higher values than Papilionoideae, underscoring the divergent gene structures present in different species lineages.

Type II genes are structurally complex, with an average length of 8,325 bp and typically containing 6 to 8 introns; some subfamilies, such as AGL17 (12,673 bp), SOC1 (13,148 bp), and FLC (12,916 bp), are particularly long due to extensive introns (Figure 2F). In stark contrast, Type I genes are structurally streamlined, with averaging over 9 introns per gene. Their average gene length is only 1,251 bp, and a remarkable 68.6% of genes are intronless and 91.04% have one or fewer introns. This pattern holds across Type I subfamilies, with 75.38% of Mα genes, 65.58% of Mβ genes, and 61.20% of Mγ genes lacking introns entirely (Figure 2F).

Motif analysis of the MADS domain further reinforces this dichotomy (Figure 3). Key MADS domains SRF-like domain (CDD:238166) were identified in Type-I (Figure 3A) and MEF2-like domain (CDD:238165) were identified in Type-II (Figure 3B). The MADS domain, typically ranging from 71 to 83 amino acids in length, is highly conserved and merits a detailed analysis of its sequence. All Type I proteins contain a highly conserved SRF-like domain (CDD:238166) characterized by two consensus motifs: “MGRKKIELKKISNDSARKVTFSKRKKGLFKK” and “ASELSTLCGVEACAIVFSPGD”. All Type II proteins feature a MEF2-like domain (CDD:238165), also defined by two motifs, but their primary consensus sequences “MGRGKIEJKRIENKTNRQVTFSKRRNGLLKKAYELSVLCDAEVALIIFSS” and “TGKLYEYASSSMMEK” are more highly conserved than its Type I counterpart and shows distinct positional variations among subfamilies. For example, in the AG/STK subfamily, it starts at the 21st amino acid from the N-terminus, while in the FLC subfamily, it begins at the 11th amino acid. In the AGL6 subfamily, there are 10 amino acids with higher divergence between the two motifs. It suggests fine-scale structural evolution even within this conserved domain.

The total MADS-box gene count varies dramatically across the Fabaceae, from 60 in Cercis canadensis to 177 in Melilotus albus, reflecting a dynamic evolutionary history of gene gain and loss (Figure 4). Gene number expansion is often linked to polyploidy; the recent WGD in Glycine genus, which contrasts with other species. For example, the Type-I subfamily in soybean comprises 70 MADS-box genes, with the Mα, Mβ, and Mγ subfamilies containing 38, 10, and 22 genes, respectively. The Type-II subfamily has a total of 96 genes, with the number of genes in different subfamilies ranging from 2 to 12. Furthermore, we observed specific expansions in some subfamilies in some species. Specifically, the number of genes in the Type-I Mα and Mγ subfamilies in Glycine, Medicago, Melilotus, and Lupinus is significantly higher than in other species. Similarly, the number of genes in the SVP and AGL17 subfamilies of Type-II in Prosopis exceeds that in other species.

Conversely, gene loss has also sculpted this family. TM8 genes were lost in all Fabaceae species surveyed but present in other species, such as grape and Amborella. FLC orthologues were mainly detected in Glycine, Dalbergia (Figure 4) and BS, as a sister group of AP3/PI, also lost in 6 species (Acacia crassicarpa, Bauhinia championii, Biancaea sappan, Delonix regia, Senna siamea, Vicia faba, and Vigna umbellata), suggesting recurrent, lineage-specific gene loss events.

To unravel the mechanisms behind this dynamic history, we traced the duplication origins of all genes (Figures 5, 6). This revealed a striking dualistic evolutionary model. We categorized the MADS-box genes into six duplication types, including WGD, tandem, dispersed, transposed, proximal, and singleton (Supplementary Figure S1). Among these, WGD was the predominant type in legumes, accounting for 42.2% of the cases. This was followed by dispersed duplication (16.2%) and tandem duplication (16.1%), which exhibited comparable proportions.

To understand the molecular basis for this dualistic pattern, we contrasted the selection pressures on 912 unique WGD-derived gene pairs and unique 635 SSD-derived gene pairs. Analysis of evolutionary rates revealed that WGD pairs exhibited significantly higher synonymous (Ks) and non-synonymous (Ka) substitution rates (Figure 7A). The median Ks values of WGD pairs were substantially higher than those of SSD pairs (Wilcox test, P < 0.001), indicating that the former diverged earlier.

To reconstruct the fine-scale temporal history of these lineages, we further analyzed the distribution of Ks across subfamilies (Figure 7B). This analysis revealed a clear temporal stratification. The conserved Type II subfamilies (e.g., SEP, AP1/FUL, AG/STK) exhibit a bimodal Ks distribution. In addition to recent duplications, they display a distinct concentration of gene pairs peaks. This suggests that the expansion of the core developmental machinery was a direct legacy of this ancient polyploidy event.

Consistent with these distinct evolutionary trajectories, the analysis of selection pressure showed a completely opposite trend. The Ka/Ks ratios of WGD pairs were significantly lower than those of SSD pairs (Wilcox test, P < 0.001) (Figure 7A), suggesting the MADS-box genes from WGDs are functionally constrained by strong purifying selection, while dynamic genes from SSDs experience more relaxed constraints, affording them greater freedom to accumulate mutations and explore new functions.

The functional consequences of these divergent evolutionary paths are evident at the pangenome level (Figure 8). We clustered all genes into 202 orthologous gene groups (OGGs), comprising 7 core clusters (1,891 genes), 5 soft-core clusters (946 genes), 53 shell clusters (1,666 genes), 137 cloud clusters (369 genes) (Figures 8A, B). We further investigated the relationship between gene duplication types and the conservation of MADS-box genes in the legume pangenome (Figures 8C, D). We found that WGD-derived genes were enriched in the core gene set (57.80%), with their prevalence decreasing toward the cloud category (18.70%). Conversely, SSD-derived genes exhibited an opposite trend: they were most abundant in the cloud (81.30%) and least frequent in the core (42.20%). These results suggest that WGD is a major contributor to highly conserved genes, while SSD-derived genes are more often variable or accessory. In addition, our investigation of evolutionary conservation at the subfamily level revealed that the “ABCDE model” genes—AG/STK, AGL12, AGL17, AP1/FUL, SEP, and SOC1—were predominantly classified as core or soft-core genes, The near-universal presence across the 52 species highlights their critical and indispensable functional roles. For example, the high composition proportion of core genes observed in the AGL12 (100% core genes), SEP (97.664%), AGL17 (95.05%), SOC1 (94.75%), AP1/FUL (93.79%), and AG/STK (81.39%). Conversely, subfamilies like Mβ, AGL6, AGL15 and FLC subfamilies showed lower conservation, predominantly comprised of shell genes. These results reflect their rapid turnover and variable presence across lineages, and mirror our dualistic model.

A comprehensive global analysis of 3,638 soybean transcriptomes, sourced from diverse tissues under various conditions, has unveiled four distinct organ-specific modules with higher expression levels (Figures 9A, B). Module 1, comprised of 31 genes predominantly from the SEP and AG/STK families, displayed expression highly specific to reproductive tissues (seed coat, flower, endosperm) and but absent in leaves and roots, likely regulating seed maturation and floral organogenesis. We also found a leaf-specific module of seven genes (Module 2), likely involved in photosynthesis and flowering time control. A distinct root-specific module of 12 genes (Module 3), enriched in AGL12 and AGL17-like genes, points to specialized roles in nutrient uptake and root development. Module 4 comprises 13 genes, predominantly from the SOC1 and SVP families, that are widely expressed in both roots and leaves. Our investigation into the cross-module distribution of gene duplication states revealed that the different copies of duplicated genes were systematically partitioned into different modules except for Module 1 (Figure 9C). This pattern implies that expression differentiation in root or leaf has occurred for these duplicated genes.

Furthermore, we focused on the expression divergence between the two main types of gene duplication: tandem duplication and whole-genome duplication (WGD). We focused on tandem duplication as the primary SSD representative because both WGD and tandem duplication are mechanistically distinct and algorithmically clear-cut duplication modes (i.e., syntenic blocks vs. local arrays). In contrast, other SSD categories like dispersed represent heterogeneous, catch-all classifications, making them less suitable for a robust comparison. Tandem duplicated pairs generally showed low expression across most tissues but were specifically upregulated in the endosperm (Figure 9D; Supplementary Table S5). In contrast, WGD-derived pairs exhibited higher and more coordinated expression across multiple tissues, particularly in flowers, seed coats, and pods, supporting their role in conserved developmental programs (Figure 9E). We also assessed the distribution of “high-expression, large-divergence” gene pairs—defined as those with expression levels > 10 and fold-change > 5—across different tissues. The highest number of such pairs was identified in the endosperm (17 pairs), followed by the nodule (11 pairs), and the root and seed coat (10 pairs each) (Supplementary Table S6).

In addition, we investigated the evolutionary and expression divergence within soybean transcriptomes in response to stress treatments. Under abiotic stress conditions, WGD-derived gene pairs sustained elevated overall expression levels, yet exhibited marked divergence among paralogs, particularly in response to ozone, salt, and drought stress (Figure 9F). This pattern of asymmetrical expression was even more pronounced under biotic stress. Notably, WGD-derived gene pairs demonstrated exceptionally high expression in response to infections by Sclerotinia sclerotiorum and Rotylenchulus reniformis (Figure 9G). Under biotic stress conditions, the average fold change of 3.24 for WGD pairs higher than the 2.76-fold change observed under abiotic stress conditions (Supplementary Tables S7, S8). This indicates that not only are asymmetric expression patterns prevalent under biotic stress, but they are also more pronounced compared to abiotic stress.

We focused on six MADS-box genes responding to abiotic stress treatment (Figure 10). The Glyma.08G05400, which is broadly expressed across multiple tissues within module 4, and the Glyma.19G034600, which shows leaf-specific high expression in module 2. Two pairs of WGD-derived genes displaying marked expression divergence under abiotic stress—Glyma.05G227200 vs. Glyma.08G033900 and Glyma.01G020500 vs. Glyma.09G201700—were selected for transcriptome validation. Expression analysis revealed that Glyma.08G05400 was strongly induced under salt stress, with a more than 110-fold increase in leaves after 12 hours relative to the control. Under drought stress, this gene was significantly activated at 3 hours and peaked at 8 hours, exhibiting over 27-fold upregulation. The Glyma.19G034600 was notably induced at 3 hours under drought stress and reached a maximum at 8 hours, with expression levels exceeding 150-fold that of the control. Among the duplicated gene pairs, both Glyma.05G227200 and Glyma.08G033900 were upregulated under drought stress, but their induction levels differed substantially: Glyma.05G227200 showed over 66-fold induction, whereas Glyma.08G033900 reached a maximum of only 25-fold. A similar divergence was observed under salt stress, where Glyma.05G227200 was induced more than 9.13-fold compared to only 3.49-fold for Glyma.08G033900. We further analyzed the duplicated pair Glyma.01G020500 and Glyma.09G201700. Under salt stress, both genes were comparably induced, with expression increases exceeding two-fold. Under drought stress, however, Glyma.01G020500 was upregulated over 15-fold, while Glyma.09G201700 was downregulated. Collectively, these results clearly demonstrate expression divergence between duplicated gene copies under stress conditions and provide experimental evidence for functional differentiation following whole-genome duplication.

Our comprehensive analysis of the MADS-box gene family across 52 Fabaceae species reveals a profound dualistic evolutionary pattern that provides a compelling solution to the post-polyploidy paradox of stability versus innovation. To anchor this model in a temporal context, our Ks distribution analysislocalized the primary burst of WGD-derived gene retention to approximately 59 Mya. Notably, this timing places the shared legume WGD event in the immediate aftermath of the K-Pg mass extinction (Koenen et al., 2021; Vanneste et al., 2014). This temporal coincidence supports the hypothesis that polyploidy served as a critical survival strategy, providing the ancestral legume lineage with a surplus of genetic raw material to rapidly colonize ecological niches left vacant by the extinction event.

Crucially, however, the explosive radiation of legumes following this K-Pg boundary WGD was not driven by a monolithic expansion of its master developmental regulators. Instead, we propose a two-speed evolutionary architecture: a conserved core of Type II (MIKC) genes, safeguarded by strong purifying selection after being duplicated by WGD, provided developmental stability. Simultaneously, a dynamic periphery, primarily composed of Type I genes, rapidly evolved through SSDs under relaxed selection, offering a continuous source of genetic novelty. This division of evolutionary roles allowed the ancestral legume genome to reconstruct order from the genomic plasticity of WGD, ensuring the robustness of essential programs while fostering the adaptive potential that fueled its subsequent success.

Our finding that core MIKC-type subfamilies expanded primarily via WGD provides strong support for the gene dosage balance hypothesis (Birchler and Veitia, 2014) (Birchler and Veitia, 2010; Conant et al., 2014; Shi et al., 2020). MADS-box proteins rarely act alone; they function by assembling into stoichiometrically precise multiprotein complexes, such as the heterotetramers of the floral quartet model, to regulate downstream gene networks (Sheng et al., 2019; Ye et al., 2021).

According to the dosage balance hypothesis, the random, small-scale duplication of a single component would disrupt this delicate stoichiometry, leading to non-functional complexes and deleterious phenotypes, and would thus be purged by strong purifying selection (Manzoor et al., 2024; Qiao et al., 2019; Wilson and Liberles, 2023). In contrast, a WGD event duplicates the entire network proportionally, preserving the balance and providing a viable pathway for the coordinated expansion of these interconnected regulatory modules (Clark and Donoghue, 2018; Soltis et al., 2015).

Our results, which show both the preferential retention of these genes post-WGD and the intense purifying selection acting upon them (low Ka/Ks ratio), empirically validate this theory. This demonstrates how a fundamental biochemical constraint—the need for stoichiometric balance—acts as a powerful evolutionary filter, preserving the integrity of core developmental machinery through deep time.

In stark contrast to the conserved core, the dynamic periphery driven by SSDs acts as a crucible for adaptive evolution (Defoort et al., 2019; Ezoe et al., 2021; Glover et al., 2015; Magadum et al., 2013). Our results show that Type I subfamilies, along with specific Type II lineages have undergone rapid, lineage-specific expansions via tandem and other SSDs. This mode of evolution, coupled with the relaxed purifying selection we observed, creates a “birth-and-death” scenario where new gene copies are constantly generated, providing raw material for neofunctionalization or subfunctionalization (Dort et al., 2024; Zhao et al., 2015). Crucially, our analysis suggests that this relaxed selection facilitates distinct functional trajectories depending on the lineage. For the enigmatic Type I subfamily, rapid evolution appears channeled primarily into reproductive specialization. Consistent with their preferential expression in seeds and gametophytes, the high sequence divergence of Type I genes may be driven by the parental conflict hypothesis or the need to establish rapid reproductive barriers, rather than abiotic stress tolerance in general (Mora-Garcia and Goodrich, 2000). In contrast, the functional shift towards environmental resilience is most prominent in SSD-expanded Type II lineages.

A compelling case for this environmental adaptation is the massive tandem expansion of the SVP subfamily in the arid-adapted genus Prosopis. As SVP homologs are known integrators of environmental signals that control flowering time, we hypothesize that this expanded repertoire of SVP genes may have allowed Prosopis to evolve a more sophisticated and resilient regulatory network to fine-tune its reproductive strategy in a harsh, unpredictable desert environment (Quesada-Traver et al., 2022; Shu et al., 2013). This dynamic nature is further reflected in the physical gene architecture itself. The structural analysis revealed that the Papilionoideae lineage possesses significantly more streamlined genes (i.e., fewer and shorter introns) compared to the more ancestral, intron-rich structures seen in Caesalpinioideae. This derived trait of genomic compaction, particularly prevalent in the Type I periphery (Bemer et al., 2010; De Bodt et al., 2003), would lower the energetic and temporal costs of transcription. Such streamlining acts as a powerful evolutionary facilitator for the rapid “birth-and-death” cycles and adaptive exploration that characterize these SSD-driven “accessory” genes. This provides a powerful example of the evolutionary chain linking genome structural variation, gene family expansion, and species-specific ecological adaptation.

Our findings reshape the understanding of WGD’s role, casting it not as a deterministic endpoint but as a stochastic “starting point” that triggers divergent, contingent evolutionary trajectories. The stark contrast between Bauhinia, which largely retained its WGD duplicates, and its close relative Cercis, which lost them and relied on dispersed duplications for expansion, vividly illustrates this post-polyploidy divergence. This discrepancy likely reflects fundamental differences in genomic stability and evolutionary strategy. Cercis is often characterized as possessing a highly stable genome, closely resembling the ancestral legume karyotype (Hyun-oh et al., 2024). This suggests that Cercis may operate under a conservative evolutionary regime where gene dosage balance is strictly enforced, leading to the rapid fractionation of WGD duplicates to restore a diploid-like state.

Conversely, Bauhinia (and the Papilionoideae lineages) appears to have leveraged genomic plasticity, retaining a larger repertoire of WGD-derived regulators. This retention may have provided the necessary genetic modularity to evolve complex traits suited for diverse tropical environments. Thus, the fate of WGD genes is not uniform but is sculpted by lineage-specific constraints—balancing the immediate cost of genomic instability against the long-term benefit of adaptive potential.

Crucially, this phenomenon of contingent evolution is visible even at finer scales. Within the genus Prosopis, two closely related species employed entirely different SSDs strategies (tandem vs. transposed/dispersed) to shape their MADS-box repertoires. This demonstrates that the final architecture of a gene family is not predetermined by WGD alone. Rather, it is the product of a complex interplay between lineage-specific selective pressures and the particular modes of SSDs that predominate, confirming that the period of massive gene loss and rearrangement following a WGD is a critical window for evolutionary innovation (Ren et al., 2018; Wilson and Liberles, 2023).

The functional data from our Glycine max case study provide a compelling, present-day validation of this evolutionary narrative. The expression patterns of duplicated genes directly reflect their divergent evolutionary origins. WGD-derived pairs within core subfamilies exhibit high, stable, and coordinated expression in key developmental organs, confirming their role as the conserved machinery (Carretero-Paulet and Fares, 2012; Zhao et al., 2020; Zhu et al., 2013). In contrast, the rampant asymmetric expression observed in many duplicated pairs, especially under abiotic and biotic stress, is a clear signature of functional divergence (Gu et al., 2004; Ha et al., 2007; Wei et al., 2025). The divergent or opposing stress responses of duplicated pairs (e.g., Glyma.05G227200/Glyma.08G033900 and Glyma.01G020500/Glyma.09G201700) suggest their post-duplication subfunctionalization or neofunctionalization (Qian and Zhang, 2014; Zhao et al., 2020). This extensive expression divergence generates a reservoir of genes with novel expression patterns, providing the functional basis for the adaptive plasticity that has allowed legumes to thrive in diverse environments (De Smet et al., 2017; Ebadi et al., 2023; Ha et al., 2007; Kou et al., 2022).

It is important to note that the “Dualistic Model” proposed here represents a dominant evolutionary trend rather than a rigid dichotomy. We observed instances where specific lineages defy general rules under unique selective pressures. A prime example is the SVP subfamily (Type II) in the arid-adapted genus Prosopis. Although Type II genes are typically constrained to the conserved core driven by WGD, the SVP lineage in Prosopis has “crossed the boundary” into the dynamic periphery, undergoing SSD. This suggests that the structural constraints preventing Type II SSD expansion are not absolute. When the adaptive value of diversifying a specific regulator outweighs the cost of genomic instability, evolution can override these constraints.

While our study provides a comprehensive evolutionary framework, these exceptions highlight that our current insights represent a starting point rather than a conclusion. Our conclusions on the adaptive significance of specific gene expansions are primarily based on genomic correlation and expression data, while our deep transcriptomic insights are confined to Glycine max. Furthermore, our focus on a single gene family necessarily simplifies what was undoubtedly a genome-wide phenomenon. Therefore, the most critical next step is to move from evolutionary inference to direct functional validation. CRISPR/Cas9-based functional genomics will be essential for testing the key adaptive hypotheses proposed here, such as dissecting the role of the massively expanded SVP subfamily in the drought tolerance of Prosopis and characterizing the highly stress-responsive, asymmetrically expressed gene pairs from soybean to confirm their roles in climate resilience.

Beyond direct validation, our findings open avenues to address fundamental questions in legume biology and evolution. A major challenge is to unravel the enigma of the rapidly evolving Type I MADS-box genes. An even broader frontier lies in understanding how the MADS-box regulatory network co-evolved with other key legume innovations, such as symbiotic nitrogen fixation, by exploring how these master regulators of root development were potentially co-opted to enable the formation of a novel organ—the nodule.