Costus zingiberoides raw Illumina reads were downloaded from NCBI’s SRA repository under accession number SRR18516547 (Valderrama et al., 2022) and then subjected to genome assembly using MaSuRCA v.4.1.0 (Zimin et al., 2013) with default parameters. This approach allows creation of “super-reads”, which improve assembly of repetitive regions and low coverage from short-read data. Although contiguity remains limited, super-reads help preserve genetic regions of interest, making it a recommended approach for species without additional transcriptome evidence. Quality control was verified by FastQC v0.12.0 program (Wingett and Andrews, 2018) followed by adapter removal and quality filtering using Trimmomatic v0.39 program (Bolger et al., 2014) with the following parameters: ILLUMINACLIP: TruSeq3 adapter file (max 2 mismatches, min length 30 bp), SLIDINGWINDOW:5:30, LEADING:5, TRAILING:5, MINLEN:100, and AVGQUAL:30. Scaffolds shorter than 1000 bp were removed using SeqKit v0.13.2 (Shen et al., 2016). For C. bracteatus, we used the genome available on the CoGe platform (https://genomevolution.org/coge/) under project ID 63654 assembled by Valderrama et al. (2022) following a hybrid genome assembly approach with MaSuRCA v3.4.2 (Zimin et al., 2013) to utilize both the Illumina short reads and the error-corrected Nanopore long reads. Default parameters were used in the assembly, and the final assembly was polished using POLCA (Zimin and Salzberg, 2020) as implemented in MaSuRCA v3.4.2.

The AUGUSTUS v3.1 program (Stanke et al., 2004) was used for gene prediction using the BUSCO v5.5.0 pipeline (Waterhouse RM, et al., 2018). The completeness of the assembled genome was evaluated with a predefined set of universal orthologs available in the embryophyta.odb9 database from OrthoDB v9.1 (Zdobnov et al., 2017). To identify putative homologous sequences between both species, we employed a reference-guided reciprocal BLAST strategy. Predicted proteins from both Costus species were independently searched against the well-annotated Musa acuminata reference proteome (UniProt ID UP000012960) using DIAMOND v2.1.10.164 (Buchfink et al., 2014). Candidate homologs were then identified based on shared best hits to the same Musa proteins, allowing the identification of corresponding genes between the two Costus species. The analysis employed default parameters, with an e-value of ≤ 1 × 10^–5, an identity of ≥ 70%, coverage of ≥ 80%, and consideration of the highest bitscore among candidate sequences.

We retrieved sequences of twelve candidate genes previously reported as involved in floral development, including shape, color, size, and organ arrangement (Valderrama et al., 2022). Sequences were obtained primarily from GenBank, prioritizing those that met predefined reliability criteria established through BLAST search parameters (see Supplementary Information 2) (Altschul et al., 1990). The AGAMOUS sequence from Chamaecostus fusiformis was provided by the authors of the PRJNA738961 bioproject.

We created separate alignments for each candidate gene (AGAMOUS, PISTILATA, YABBY2, BIG BROTHER, LONGIFOLIA-like, AUXIN RESPONSE FACTOR 2-like, RADIALIS, AGAMOUS-LIKE6, CUP-SHAPED COTYLEDON-like, CHALCONE SYNTHASE-like, DIHIDROXIFLAVONOL 4-REDUTASE, and FLAVONOL SYNTHASE). Initially, the coding DNA sequences (CDS) were aligned using MAFFT v7.520 (Katoh and Standley, 2013) with the parameters –genafpair –maxiterate 1000. These CDS alignments were then translated into amino acid sequences to guide structural and comparative analyses. All alignments were subsequently manually refined using MEGA11 (Tamura et al., 2021) to ensure correct codon correspondence.

All maximum likelihood phylogenetic trees were inferred based on DNA alignment with IQ-TREE v2.0.7 (Minh et al., 2020) using the substitution model selected by ModelFinder (Kalyaanamoorthy et al., 2017). Branch support was assessed with 1000 bootstrap replicates. To place Costus zingiberoides and Costus bracteatus within a phylogenetic framework of land plants and to estimate their divergence time, we used conserved single-copy orthologs from the embryophyta_odb10 dataset available on BUSCO (Waterhouse RM, et al., 2018), which represents a curated collection of orthogroups across Embryophyta. Only genes without in-paralog copies were used. The Embryophyta species tree was inferred using ASTRAL (Mirarab et al., 2014), which employed individual gene trees generated by IQ-TREE v2.0.7 (Minh et al., 2020). Divergence times were estimated with the RelTime method implemented in MEGA11 (Tamura et al., 2021) following (Morris et al., 2018).

Scanning for episodic positive selection on selected genes was performed using the aBSREL method (Smith et al., 2015), as implemented in HyPhy v2.5.7 (Kosakovsky Pond et al., 2020). This method implements a branch-site model, identifying both amino acid sites and branches of the phylogeny with d_N/d_S rate ratio significantly greater than one by fitting model complexity to the dataset by determining the optimum number of rate categories for each branch using the Akaike Information Criterion.

The three-dimensional modeling of proteins was conducted through homology modeling using the Swiss-Model web server (https://swissmodel.expasy.org/; (Waterhouse A, et al., 2018). This platform performs homology-based modeling constrained to experimentally solved templates of closely related proteins, allowing greater control over conserved domains and ensuring model comparability among genes and species (Kiefer et al., 2009; Waterhouse A, et al., 2018). The analysis focused on Costus bracteatus and Costus zingiberoides, and included proteins encoded by genes showing evidence of positive selection, as well as homologous proteins from phylogenetically related species that did not exhibit positive selection, allowing a comparative assessment of structural divergence. For each protein, the best structural template was automatically selected based on the highest statistical scores provided by Swiss-Model (e.g., GMQE and QMEAN values), optimizing model reliability and accuracy (Kiefer et al., 2009).

Model quality was assessed using the Ramachandran plot (Ramakrishnan and Ramachandran, 1965) and QMEAN scores (Benkert et al., 2009) implemented in the Swiss-model webserver. The protein domains were predicted in silico using the PFAM database (Mistry et al., 2021) and structural alignment was carried out with UCSF Chimera v1.17.3 (Pettersen et al., 2004) using only the protein domains identified, as this allows us a more accurate comparison, while the N- and C-terminal regions can exhibit high variability and flexibility. Additionally, the Root Mean Square Deviation (RMSD) of alpha-carbon atoms (Cα) was calculated using UCSF Chimera’s “Match–align” tools to quantify structural similarity. RMSD-Cα provides an average measure of deviation between equivalent atoms in superimposed protein structures, serving as an indicator of conformational differences. Lower RMSD values suggest a higher structural similarity.

To evaluate the impact of amino acid replacements on protein structure under positive selection, an alanine scanning analysis was performed using the FoldX software (Guerois et al., 2002). This software estimates the free energy (ΔG) of proteins using a combination of empirical factors, such as electrostatic interactions, van der Waals forces, hydrophobicity, solvation effects, and hydrogen bonding. The mutational effect (ΔΔG_fold) is calculated as the difference between the free energy of the wild-type (ΔG_wt) and the mutant structure (ΔG_mut) using the formula:

A ΔΔG value less than 0 indicates a stabilizing mutation, while values greater than 0 suggest destabilization. Mutations were categorized into five stability-impact groups: highly stabilizing (ΔΔG < −1.84 kcal·mol^–1), slightly stabilizing (−1.84 to −0.46 kcal·mol^–1), neutral (−0.46 to +0.46 kcal·mol^–1), slightly destabilizing (+0.46 to +1.84 kcal·mol^–1), and highly destabilizing (ΔΔG > +1.84 kcal·mol^–1). This classification helps interpret the structural and functional consequences of mutations on protein stability.

Due to the lack of a publicly available genome assembly for C. zingiberoides, we performed a de novo assembly using Illumina short reads, which resulted in 76,483 contigs with an N50 of 16 kb while the number of scaffolds was 54,383 with an N50 of 20 kb. The total length of the genome was 706,297,564. The completeness assessment, including the percentage of completeness, the number of single-copy genes, and other relevant metrics, is provided in the supplementary materials (Supplementary Figure 1). The C. bracteatus genome assembled by Valderrama et al. (2022) after polishing resulted in 13,709 contigs with an N50 of 142,42 kb and BUSCO analysis showed 1,291 out of 1,375 complete BUSCOs with 1,163 single-copy and 128 duplicated, with an additional 33 fragmented BUSCOs and 51 missing BUSCOs.

Evolutionary relationships were inferred for twelve genes associated with floral development and morphological traits directly linked to pollination syndromes. These genes—LONGIFOLIA-like, FLS, PISTILLATA, RADIALIS, YABBY2, AGL6, CUC-like, AGAMOUS, AUXIN RESPONSE FACTOR 2-like, BIG BROTHER, CHALCONE SYNTHASE-like, and DIHIDROXIFLAVONOL 4-REDUTASE—were selected based on their regulatory roles in floral form and function. The dataset encompasses representatives from eight distinct monocot orders, and includes selected eudicot species from Arabidopsis thaliana, A. lyrata, Raphanus sativus, Brassica rapa, Camelina sativa, and Antirrhinum majus, which were used as outgroups to root individual gene trees and improve the inference of evolutionary directionality. The distribution of species per gene varied, with the number of sequences for each gene depending on the availability in GenBank (Supplementary Information 2).

We placed both Costus within a species-level phylogeny constructed from conserved single-copy orthologs retrieved from the embryophyta_odb10 dataset which includes well-characterized orthologous groups shared across land plants. In this framework, both Costus species were consistently recovered as sister taxa to Musa. Divergence time estimates calibrated using the split between Amborella and the remaining angiosperms (minimum age 125.0 Ma, maximum 247.2 Ma), resulting in a divergence time of approximately 5.7 million years ago for the two Costus species (Supplementary Figure 3). This phylogenetic framework supports subsequent analyses of lineage-specific evolutionary patterns among floral developmental genes.

Out of the twelve candidate genes analyzed (Supplementary Table 1), four exhibited evidence of episodic positive selection on at least one lineage, namely, the LNG-like (Zingiber officinale), AGL6 (Costus zingiberoides), CUC-like (Oryza sativa and Setaria viridis), and AGAMOUS (Hordeum vulgare, Zea mays, Lacandonia schismatica, Elaeis guineensis, Alpina oblongifolia and Lilium longiflorum) coding sequences (Figure 2). These were nonetheless included for downstream structural and mutational analyses in silico to enable comparative evaluation between Costus and other lineages under selection, thereby identifying both shared and unique molecular features. This approach allows a broader understanding of how different selective regimes may have shaped the diversification of floral traits and pollination strategies across monocots.

Three-dimensional protein structures were modeled with a targeted, hypothesis-driven strategy focusing on genes and lineages most relevant to pollination-related divergence. Structural models were generated for Costus bracteatus and C. zingiberoides, as well as for homologous proteins from selected species whose genes showed evidence of positive selection in the phylogenetic analyses. To provide a structural baseline, a limited set of phylogenetically related species lacking signatures of positive selection was also included.

Root Mean Square Deviation (RMSD) analyses were performed jointly among all modeled proteins corresponding to each gene, in addition, leave-one-out RMSD comparisons were performed to evaluate the relative structural divergence of the two Costus species to maximize biological interpretability by highlighting lineage-specific structural changes potentially associated with adaptive evolution, while avoiding unnecessary expansion of structurally conserved comparisons that would not directly inform the study’s central hypotheses. This strategy enabled a more refined assessment of lineage-specific structural changes associated with differential selective pressures.

The following genes exhibited a high degree of structural conservation between C. bracteatus and C. zingiberoides, despite their contrasting pollination syndromes, suggesting strong functional constraints on protein architecture.

The LNG-like protein was modeled for four species, Musa acuminata (not under positive selection), Costus zingiberoides (not under positive selection), Costus bracteatus (not under positive selection) and Zingiber officinale (under positive selection). The RMSD Cα values indicated moderate structural divergence among homologs (RMSD = 1.06; Figure 3). When comparing each Costus protein to the set of non-Costus homologs, both C. zingiberoides and C. bracteatus exhibited comparable RMSD values (1.23 Å for both species), suggesting a similar degree of structural conservation relative to other lineages. Structural alignment highlighted several regions of divergence between species, particularly in the IPR025486 domain of LNG-like protein and with unknown function identified by the Pfam webserver, spanning from AA 919 to AA 1100 for M. acuminata, AA 820 to 984 for Z. officinale, and from 890 to 1065 for C. bracteatus and C. zingiberoides (Figure 4).

In silico alanine scanning was performed to evaluate the impact of single-residue substitutions on protein stability, revealed that most alanine substitutions in LNG-like sequences are destabilizing for both proteins, mainly in residues of its domain (Figure 4). However, when matching the highlighted amino acid replacements with the alanine-scanning energy values, three specific sites in Z. officinale showed a stabilizing effect on the protein structure — Ser842 (+0.93 kcal·mol⁻¹, position 31), Asn926 (+0.93 kcal·mol⁻¹, position 146), and Thr929 (+0.93 kcal·mol⁻¹, position 149). In contrast, the residues at the corresponding positions in the Costus species generally exhibited destabilizing effects upon alanine substitution, suggesting that these stabilizing mutations may contribute to structural optimization unique to Z. officinale.

AGL6 proteins were modeled for Dendrocalamus latiflorus, Cymbidium faberi, and C. bracteatus, as well as C. zingiberoides. Although the former species are more distantly related, their proteins yielded higher-quality structural models, thereby enabling more robust comparative analyses. The overall RMSD-Cα among all modeled proteins indicated moderate divergence (0.858 Å; Figure 3). When compared individually to non-Costus homologs, C. zingiberoides (0.847 Å) and C. bracteatus (0.876 Å) displayed similar structural deviation, suggesting comparable levels of conservation within the genus. The MADS-box domain of proteins across species was shown to be a highly conserved structure with a RMSD- Cα of 0.512 (Å), while its K domain showed greater structural variation, with a RMSD- Cα of 1.126 (Å) (Figure 3).

Few conformational changes were detected between Costus species in proteins that were not inferred to have undergone positive selection. Structural and mutational alignment analyses in silico showed high conservation of protein domains between the two species (Supplementary Figures 5–13). RADIALIS stood out by presenting variation in total protein length, with the C. zingiberoides sequence being shorter than that of C. bracteatus. Despite this difference, both species exhibited a similar mutational profile under alanine scanning, with a higher proportion of substitutions causing strong destabilizing effects (Supplementary Figure 8). Another noteworthy case was the CHS-like protein, where specific point amino acid replacements distinguished the two species structurally, and alanine scanning predicted that several sites in both proteins could lead to considerable destabilization if replaced (Supplementary Figure 12).

In contrast to the conserved patterns described above, some genes exhibited lineage-specific structural features between the two Costus species, potentially reflecting divergent evolutionary trajectories.

Despite the structural conservation of the domains, at least 32 mutations occurred in the AGL6 MADS-Box domain sites of C. zingiberoides (start of the domain in 2–78 aa) when compared to the other three proteins, and 22 amino acid replacements are shown exclusively in the K domain of C. zingiberoides (domain starts at 97aa and ends at 167aa) (Figure 5).

The mutational analysis in silico revealed that most alanine substitutions in the AGL6 protein sequences of C. zingiberoides resulted in either moderate or significant destabilization, particularly within residues located in the MADS-box and K domains (Figure 5). Residues whose replacement did not significantly alter stability were classified as neutral, with 21 residues in the MADS-box and 37 in the K domain falling into this category (Supplementary Table 1). Conversely, ten residues were identified as contributing to enhanced protein stability, with nine categorized as slightly stabilizing and one as highly stabilizing. Notably, seven of the ten stabilizing residues were situated within the MADS-box domain, while only one was located in the K domain, highlighting the MADS-box’s central role in maintaining structural integrity. Among the mutations observed exclusively in the MADS-box domain of C. zingiberoides, three—Cys20 (-0.48 kcal·mol⁻¹), Ser33 (-0.57 kcal·mol⁻¹), and Gly44 (-1.88 kcal·mol⁻¹)—were predicted by alanine scanning to confer additional stabilization to the protein.

In Costus bracteatus, the AGL6 protein exhibited a reduced overall length compared to other species (Figure 3), with a shorter MADS-box domain composed of 29 amino acid residues and a K-box domain with 47 residues (Figure 5). The mutational effect analysis in silico for C. bracteatus revealed a predominance of slightly destabilizing substitutions distributed along the protein sequence (Figure 5).

For CUC-like, we modeled the NAC domain (start of the domain at 22aa and ends at 175aa) in six species, Triticum aestivum, Zingiber officinale, C. bracteatus, and C. zingiberoides, and in Oryza sativa and Setaria viridis. The RMSD Cα values were low across all species, indicating high conservation of the NAC domain. Furthermore, the structural alignment (Figure 6) showed a small amount of amino acid replacements occurring in this domain, reinforcing the conservative nature of this region. When compared individually to non-Costus homologs, C. bracteatus exhibited a lower RMSD value (0.378 Å) than C. zingiberoides (0.610 Å), suggesting that C. bracteatus retained a more conserved NAC domain conformation. Structural alignments revealed minimal variation in C. bracteatus, while C. zingiberoides displayed several lineage-specific amino acid substitutions within the NAC domain, reflecting a subtle yet potentially meaningful divergence in a functionally constrained region.

The alanine scanning analysis (Figure 6) revealed a consistent destabilizing pattern across both Costus species, with predicted amino acid replacements leading to substantial destabilization, comparable to the profiles observed in the two other species whose CUC-like genes are under positive selection (Supplementary Figure 4). Amid the observed amino acid replacements, three residues—Glu (position 5), Thr (position 147), and Ser (position 160)—were conserved between the two positively selected species, whereas Costus harbors different residues at these sites. The alanine scanning outcomes revealed that the substitution at position 5 generally reduces stability across species; position 147 is neutral in O. sativa, S. viridis, and C. bracteatus, but destabilizing in C. zingiberoides (-0.474 kcal·mol⁻¹); and position 160 shows a stabilizing effect in both Costus species, contrasting with neutral effects in the positively selected species. These findings suggest that while the NAC domain is generally conserved, specific residues may act as hotspots critical for the protein’s stability and functional dynamics in different monocot species.

AGAMOUS protein was modeled for Musa acuminata, C. bracteatus, and C. zingiberoides (not under positive selection), Zea mays, Lacandonia schismatica, Elaeis guineensis, Alpinia oblongifolia, and Lilium longiflorum, the latter of which were found under positive selection. Between the two Costus species, RMSD overall values were 1.561 Å (C. bracteatus) and 1.588 Å (C. zingiberoides), suggesting slightly higher structural flexibility relative to other genes analyzed. The global RMSD-Cα for the MADS-box domain in these all species was highly conserved, with a value of 0.378 (Å). Similarly, the K-domain also showed notable conservation with an RMSD-Cα of 0.767 (Å) (Figure 3). The structural alignment of the protein sites revealed that fewer amino acid substitutions occur in the MADS-box domain than in the K domain, which presents a large number of replacements (Figure 7). When comparing the MADS-box and K domains, the MADS-box consistently showed more destabilizing amino acid replacements across most species.

Nevertheless, the AGAMOUS proteins from C. bracteatus and C. zingiberoides presented a distinct structural profile. Both species lacked the MADS-box domain entirely and retained only the K-box domain (Figure 3), which was notably conserved across taxa (Figure 7). Despite the absence of the MADS-box, the alanine scanning results for both Costus species revealed a consistent pattern of slightly destabilizing effects distributed throughout the protein sequence (Figure 7).

One of the most remarkable discoveries was the identification of positive selection acting on specific branches of the phylogenies for four genes. In LONGIFOLIA-like—they play a role in controlling flower size (Disch et al., 2006), positive selection was detected in Z. officinale, whereas AGL6 emerged as the most compelling candidate for pollination-related adaptation, showing evidence of diversifying selection in C. zingiberoides, an ornithophilous species pollinated by hummingbirds (Maas, 1977).

AGL6 is a MADS-box family gene that plays a central role in floral organ identity, carpel fusion, and perianth development across angiosperms, with particularly important functions in monocots where petaloid organs often derive from modified androecial structures (Yockteng et al., 2013). Changes in AGL6 expression or protein structure have been linked to variation in floral symmetry, organ fusion, and reproductive morphology. These findings align with the hypothesis that shifts in pollination syndromes may be associated with substitutions in genes regulating floral traits, as proposed by Valderrama et al. (2022).

Interestingly, the branch leading to C. zingiberoides in the AGL6 phylogeny is markedly elongated relative to other taxa, despite comparable sequence lengths among homologous sequences. This expansion may suggest an episode of accelerated molecular evolution, potentially shaped by adaptive divergence associated with shifts in pollination strategy (Frachon and Schiestl, 2024).

The differential selection regime acting on Z. officinale and C. zingiberoides may reflect adaptive changes in flower size regulation and carpel fusion in response to the specific requirements of their pollinators, with bees and birds favoring different floral characteristics (Fenster et al., 2004; Fattorini and Glover, 2020).

The CUP-SHAPED COTYLEDON (CUC) genes belong to the NAC family of transcription factors and are primarily involved in establishing organ boundaries, maintaining meristem, and separating floral organs during development (Hibara et al., 2006). Alterations in the function of boundary genes can influence floral architecture, organ fusion, and overall flower shape.

CUC-like selection tests showed evidence of positive selection in O. sativa and S. viridis—species that are mostly self-pollinated but can also rely on wind pollination (Muhammad et al., 2022). Although wind-pollinated species usually face weaker selective pressure on floral traits (Friedman and Barrett, 2009; Fattorini and Glover, 2020), this pattern suggests that CUC-like may also contribute to broader developmental roles beyond floral morphology, possibly linked to boundary formation and other processes essential for grass fitness.

AGAMOUS is a canonical MADS-box gene that specifies stamen and carpel identity and contributes to floral meristem determinacy (Zhang et al., 2020). While its core functions are conserved, several monocot lineages exhibit structural or regulatory divergence in AGAMOUS-like genes (Almeida et al., 2015). We detected signatures of episodic positive selection in AG-like across multiple monocot lineages with contrasting reproductive strategies (Devos et al., 2005; Sakai et al., 2013). In Costus, the observed truncation and domain differences should be interpreted cautiously, as they may reflect either lineage-specific divergence or limitations of draft assemblies/annotation. Nonetheless, these features are consistent with previously reported cases of altered selective regimes and domain evolution in floral developmental regulators (Auge et al., 2019; Dellinger, 2020).

Taken together, these results indicate that divergence between C. bracteatus and C. zingiberoides is not characterized by widespread adaptive evolution across floral developmental genes, but instead by lineage-specific selection acting on a restricted subset of key regulators.

The structural comparison of LNG-like proteins revealed a high degree of conservation among the two Costus species, consistent with a stable evolutionary pattern previously described for LONGIFOLIA-like proteins in Arabidopsis (Lee and Kim, 2018). This structural stability suggests that LNG-like proteins in Costus have maintained a conserved role in regulating organ elongation and floral morphology, which may reflect the shared developmental constraints. In contrast, Z. officinale exhibited specific substitutions—Ser842, Asn926, and Thr929—that slightly stabilized its protein conformation, distinguishing it from the other species, whose equivalent residues were generally destabilizing upon alanine mutation in silico. While the biological significance of this subtle stabilization is uncertain, it may reflect lineage-specific variation in protein flexibility, potentially related to shifts in floral traits observed across Zingiberales (Yockteng et al., 2013; Almeida et al., 2018; Valderrama et al., 2022). These observations highlight structural differences among LNG-like proteins across taxa, but further experimental or comparative analyses would be required to determine their functional or evolutionary implications.

In C. bracteatus, the NAC domain remains highly conserved, with minimal divergence compared to other monocots not under positive selection, consistent with a functionally constrained developmental role. In contrast, C. zingiberoides harbors more divergent residues in this region, suggesting lineage-specific modulation of boundary gene function. Notably, alanine scanning revealed a predominantly destabilizing mutational profile in both species, resembling patterns seen in positively selected lineages and pointing to structural sensitivity at key residues.

Remarkably, three positions that are conserved in positively selected species—Glu3 (position 5), Thr153 (position 147), and Ser162 (position 160)—correspond to distinct residues in C. zingiberoides and C. bracteatus, with varying effects on stability. These changes may alter the interaction of CUC-like proteins with other NAC family members or downstream targets, thereby affecting traits such as floral architecture, organ fusion, and meristem development. Because pollination syndromes appear less directly implicated, it is plausible that adaptive modifications in CUC-like reflect broader developmental processes, including boundary formation and organogenesis. This interpretation is consistent with evidence that boundary genes like CUC-like act synergistically with other transcription factors to establish spatial patterns essential for plant development (Takada et al., 2001; Nakamura et al., 2024).

Analysis of the structure of AGL6 revealed notable differences in the MADS-box and K-domains across species. The MADS-box domain exhibited a high degree of conservation (RMSD-Cα 0.512 Å), whereas the K-domain showed greater variability (RMSD-Cα 1.126 Å), particularly in C. zingiberoides, which presented 32 unique amino acid replacements in the MADS-box and 22 in the K-domain when compared to C. bracteatus, D. latiflorus and C. faberi (Figure 5). In contrast, C. bracteatus shows slightly higher RMSD values (0.876 Å) but fewer unique residues, indicating a combination of overall structural conservation with lineage-specific divergence in C. zingiberoides. These differences could potentially arise from the presence of multiple gene copies or from the assembly process, as heterozygosity or incomplete gene assembly could contribute to these observed variations (Sun et al., 2022). Further investigation, such as using allele-specific strategies, would be necessary to confirm the presence of multiple copies or other sources of variation.

Notably, while the MADS-box domain is generally considered highly conserved due to its critical role in DNA binding and transcriptional regulation (Riechmann et al., 1996; Adhikari and Kasahara, 2024), the presence of exclusive amino acid replacements in C. zingiberoides might be associated with the pollination syndromes in this lineage (Valderrama et al., 2022). Furthermore, three amino acid replacements — Cys20 (-0.48 kcal mol⁻¹), Ser33 (-0.57 kcal mol⁻¹), and Gly44 (-1.88 kcal mol⁻¹) — were identified as stabilizing, suggesting compensatory effects that help maintain structural integrity while permitting functional flexibility (Dreni and Zhang, 2016). Although no site-specific selection test was performed, these residues occur at positions that are otherwise conserved across species, which may indicate diversifying selection within C. zingiberoides. Given that amino acid replacements in the MADS domain can alter oligomerization and DNA-binding specificity, such substitutions could influence the assembly or target recognition of AGL6 complexes (Lai et al., 2019; Käppel et al., 2021), potentially contributing to lineage-specific regulatory divergence.

The K-domain of AGL6, which mediates protein-protein interactions (Fan et al., 1997; Veron et al., 2007), exhibited considerable divergence, with 15 non-consensus sites across aligned sequences. Such variability supports the view that variability in this region can promote species-specific adaptations, particularly those influencing floral morphology and pollinator interactions (Almeida et al., 2015).

In C. zingiberoides, the high number of amino acid replacements in the MADS-box relative to the K-domain of AGL6 contrasts with the typical pattern in other species, where variability is concentrated in the K-domain (Yockteng et al., 2013). This inversion may reflect lineage-specific functional shifts linked to the recurrent transitions and reversions of pollination syndromes documented in Costus (Vargas et al., 2020; Valderrama et al., 2022). More broadly, these results provide a molecular basis for the hypothesis that structural mutations in conserved domains, although rare, can be functionally necessary and positively selected to allow conformational flexibility and adaptation to changing ecological pressures (Sen et al., 2011; Piot et al., 2017). This flexibility may be particularly important for traits like carpel fusion, a key determinant of reproductive success in monocots and regulated by AGL6 (Ruelens et al., 2013).

Interestingly, Costus bracteatus exhibited an AGL6 protein model with a reduced overall length and lacking the canonical MADS-box domain, a feature that is typically conserved across monocots. However, because gene models were derived from a draft genome using ab initio and homology-based annotation, this apparent truncation cannot be unequivocally distinguished from a potential annotation or assembly artifact. Despite this uncertainty, similar truncations in MADS-type transcription factors have been reported in other species, such as maize, where lineage-specific loss or reduction of MADS domains has been associated with neofunctionalization or altered regulatory functions (Liu et al., 2020). Despite the absence of the MADS-box, the K-domain remains relatively conserved, and the alanine scanning indicated predominantly slightly destabilizing substitutions, suggesting that the truncated protein retains a degree of structural stability (Piot et al., 2017). Although this feature could reflect a novel aspect of AGL6 evolution in C. bracteatus, particularly in light of its distinct floral morphology, any functional association with pollination syndrome should be interpreted cautiously, given the speculative nature of in silico predictions and the lack of experimental validation (Gupta et al., 2014).

The AGAMOUS proteins of C. bracteatus and C. zingiberoides lack the canonical MADS-box domain but retain a conserved K-domain, representing a striking structural simplification that may indicate a divergent evolutionary trajectory within Costus. Although this absence could reflect technical limitations inherent to in silico gene recovery, similar domain losses or truncations have been reported in other plant lineages, such as grasses, where MADS-type genes experienced partial domain loss associated with relaxed selective constraints and potential neofunctionalization (Schilling et al., 2015; Liu et al., 2020).

These examples suggest that the loss of the MADS-box in Costus AGAMOUS protein may not solely reflect annotation artifacts but could represent a genuine evolutionary modification, potentially linked to subfunctionalization or compensatory regulatory mechanisms among other MADS-box paralogs. Consistent with this, alanine scanning revealed only slightly destabilizing effects in both proteins, suggesting a residual structural stability that may permit limited or alternative functionality even in the absence of the canonical DNA-binding domain, as previously proposed for other truncated MADS-type proteins (Ferrario et al., 2004).

Taken together, the observed amino acid replacements and alanine scanning results suggest that AGAMOUS contributes to floral phenotypic diversity across monocots. The interplay between stabilizing and destabilizing amino acid replacements may represent an evolutionary strategy that preserves protein stability while allowing for adaptive change. Such modifications likely facilitate adaptation to specific pollination syndromes and underscore the pleiotropic role of AGAMOUS in both floral development and diversification (Zhang et al., 2024).

Overall, our results indicate that the divergent pollination syndromes of C. bracteatus and C. zingiberoides are shaped by selective fine-tuning of a small number of key developmental regulators rather than by generalized divergence across the floral gene network. While most floral genes remain under strong functional constraint, lineage-specific selection and structural modification of proteins such as AGL6 and CUC-like likely modulate floral architecture, organ fusion, and morphology in ways that affect pollinator interaction.

Although the C. zingiberoides genome shows moderate contiguity and the gene annotation relied on ab initio and homology-based approaches—well-established for non-model species and widely applied in in silico evolutionary studies—assemblies with moderate scores can still yield robust evolutionary inferences when analyses focus on well-characterized coding sequences (Simão et al., 2015; Veeckman et al., 2016; Kersey, 2019). This way, our study providing initial evolutionary inferences presented here. Nonetheless, future work incorporating long-read sequencing, improved genome assemblies, transcriptome-based gene models, molecular dynamics and experimental validation would allow a more comprehensive assessment of gene and protein evolution across Costus species.