Bowman-Birk protease inhibitors (BBIs) are a family of serine-type protease inhibitors that play a pivotal role in plant development and defenses (Bowman, 1946; Shan et al., 2008; Dramé et al., 2013; Malefo et al., 2020; Xie et al., 2021). Recent phylogenetic analyses suggest that angiosperms share a common ancestral BBI inhibitory loop with early vascular plants that appeared in the lycopod plant lineage 370 million years ago (James et al., 2017). In modern plants, BBI-type inhibitors are mainly distributed in the Poaceae and Leguminosae families. In legumes, they possess a double inhibitory-loop structure, which was lost in the monocot lineage due to the loss of two cysteine residues, resulting in a single-headed structure (James et al., 2017). In Fabaceae, BBIs are composed of about 70 amino acids that fold to expose hydrophobic patches but lack a hydrophobic core. This structure is highly stable under different conditions, such as high temperature and a wide pH range, due to seven disulfide bridges, and forms a double-headed configuration with two independent inhibitory domains (Odani and Ikenaka, 1973; Domoney, 1999; Mello et al., 2003; Barbosa et al., 2007; Losso, 2008). Classified as I12 family in the MEROPS database (Rawlings et al., 2014)¹, these inhibitors share a binding loop motif where P1 residues serve as primary contact residues defining the interaction with the target enzyme (Supplementary Figure S1) (Schechter and Berger, 1967; Grosse-Holz and van der Hoorn, 2016). Different classes of BBI proteins, which arose from an ancestral duplication of a single-headed inhibitor gene, have demonstrated the ability to inhibit trypsin (where both reactive sites inhibit trypsin) and chymotrypsin (where the first reactive site inhibits trypsin and the second inhibits chymotrypsin). Furthermore, some studies suggest that BBI also reduces the proteolytic activities of other proteases such as elastase, cathepsin G, and chymase (Clemente and Domoney, 2006; Losso, 2008; De Paola et al., 2012; Ferreira et al., 2019).

Legume seeds are a rich source of BBIs and their activity suppresses pest and pathogen digestion, playing a defensive role conferring resistance in response to wounding by insects (e.g. Coleoptera and Lepidoptera) in different plant species (Hilder et al., 1987; Rakwal et al., 2001; Qu et al., 2003) Moreover, BBIs exhibit in vitro antifungal activity against different pathogens such as Botrytis cinerea and Fusarium culmorum (Ye et al., 2001; Pekkarinen et al., 2007) and appear to be involved in diverse biological processes including drought stress (Yan et al., 2009; Malefo et al., 2020) and tolerance to salinity (Shan et al., 2008).

In the last decade, BBIs have received a growing interest, and epidemiological studies have shown that diets based on legumes have a positive impact on human health (Sánchez-Chino et al., 2015; de Freitas et al., 2020; Carbonaro and Nucara, 2022). Although the underlying mechanism is not yet fully understood, the capacity to control the cellular growth and proliferation of some cancers (e.g. colorectal, prostate, ovarian, oral cavity) as well as the ability to suppress chronic inflammations of different origin have been widely reported (Clemente et al., 2005; Chen et al., 2005; Clemente and del Carmen Arques, 2014; Juritsch and Moreau, 2018; Gitlin-Domagalska et al., 2020). Some studies suggest that their beneficial effects are due to the inhibition of the proteasome, related to anticancer activity, and to the inhibition of proteases produced by macrophages and mast cells during inflammation, such as chymase, neutrophil elastase, and cathepsin G, directly supporting their anti-inflammatory properties (Gitlin-Domagalska et al., 2020). BBI targets are primarily members of the S1 and S3 peptidase families, and interactions with these enzymes are likely responsible for many of their biological properties (Rawlings et al., 2014)².

In recent years, a growing number of legumes have been analyzed with encouraging results showing that BBIs extracted from cowpea (Vigna unguiculata L.) seeds are good candidates for the treatment of various pathological states such as neurodegeneration, hypertension and different types of cancers such as breast, colorectal, and prostate (Joanitti et al., 2010; Souza et al., 2014; Panzeri et al., 2020; Gitlin-Domagalska et al., 2020). However, analyses of BBI expression levels as well as a general characterization of BBI genes in the genus Vigna have not yet been performed. This lack of information is surprising, given both the extent to which many communities around the world are reliant on both domesticated and wild lines of Vigna (Maxted et al., 2004; Nair et al., 2023), and the potential health benefits already observed from BBIs isolated from cowpea (Joanitti et al., 2010; Mehdad et al., 2016; Panzeri et al., 2020).

The genus Vigna (Fabaceae) includes over 100 species (Plants of the World Online, 2025), many of which are cultivated in tropical and subtropical regions (Tomooka et al., 2002; Maxted et al., 2004), with cowpea (V. unguiculata) having the largest cultivation area (Kebede and Bekeko, 2020). These species are highly adaptable to extreme environments, such as drought, poor soils, and high temperatures, making them valuable for subsistence agriculture (Nair et al., 2023). They provide food, forage, and other benefits, while wild relatives serve as genetic resources for improving stress tolerance, disease resistance, and agronomic traits, supporting the development of resilient varieties under climate change (Sharma et al., 2013; Dwivedi et al., 2016; Kumari et al., 2021). Furthermore, beans are widely consumed by African and Asian communities, who integrate these legumes into an expanding variety of traditional and modern dishes known for their affordability and sustainability (Wang et al., 2022; Singh et al., 2022; Kanishka et al., 2023; Harouna et al., 2024).

To further our understanding of the distribution of BBI genes within the genus Vigna, we used a multidisciplinary genetic, phylogenetic and computational approach and investigated their expression levels in a number of carefully selected domesticated and wild lines that are representative of their geographic and environmental distribution. Given the high adaptability of this genus to a wide variety of edaphic conditions (Nair et al., 2023) and the role of BBI proteins in plant defenses to both biotic and abiotic stresses (Clemente and Domoney, 2006), we hypothesized that there would also be a high degree of genetic variability in terms of BBI gene expression across the genus.

BBIs have been isolated and described in various legumes, however, the expression levels as well as the characterization of BBI genes in Vigna genus is lacking. In this study, we retrieve all the BBIs genes present in 12 Vigna species genomes. Then we conducted a preliminary analysis of the expression of BBI genes in five Vigna species in both leaves and seeds, as these tissues are highly relevant to both plant fitness and human nutrition. The two BBI genes mainly expressed in seeds were successively sequenced in 42 Vigna taxa, including both crops and wild relatives. The levels of interaction energy with the main human targets (trypsin and chymotrypsin) are calculated for each mature protein and the results are discussed in an evolutionary context. This study allowed us to gain deeper insight into the expression, natural variability, evolution and property of BBI genes in the genus Vigna, particularly in seeds, providing valuable information for the development of future crops.

BBIs are an important family of serine protease inhibitors, and their high rate of duplication events combined with the forces of natural selection drive the formation of divergent expression profiles and the development of new functional roles (Xie et al., 2021). Our results revealed a common expression pattern across all the seeds analyzed, evidencing that BBI1 and BBI2 genes have the highest expression level when compared with other BBI genes. This difference was confirmed in all species studied (Figure 1). The BBI1 and BBI2 genes correspond to the two genes located on chromosome 4 in the Vigna unguiculata (cowpea) reference genome available in the NCBI database⁶. In this species, their predicted inhibitory targets are trypsin–trypsin and trypsin–chymotrypsin, respectively. However, due to the uncertainty regarding target specificity across different species, in this work they are referred to as BBI1 and BBI2.

The high level of concentration of these two genes expressed in seeds, combined with their ability to inhibit different serine proteases, enhances defenses during germination and seedling growth, protecting the plant from potential pathogen attacks during its most vulnerable developmental stage. However, some authors have also interpreted that the abundance of serine protease inhibitors in seeds participate in the regulation of endogenous protease activity and serve as sulfur storage proteins because of their particularly high cysteine content (Pusztai, 1972; Tsybina et al., 2001). Several protease families are required in seeds and storage organs for protein turnover and their activity is regulated by endogenous inhibitors (Vorster et al., 2023). Therefore, the abundant concentration of BBIs in seeds might serve a triple role. Beyond providing protection against non-plant organisms, they could also regulate endogenous proteases during the mobilization of reserve proteins. Our results also revealed that the expression levels of BBIs in leaves are highly variable. This result is consistent with their ability to exhibit differential expression during various life cycles of the plant or in response to different environmental conditions (Eckelkamp et al., 1993; Malefo et al., 2020). Due to the involvement of BBIs in responses to abiotic and biotic stress, the sensitivity of BBI gene expression in leaves to environmental variations could have significant implications for agriculture (Dramé et al., 2013; Chen et al., 2024). However, additional experiments are required to evaluate the gene expression levels of BBIs across various developmental stages and environmental conditions, as well as to comprehend how to optimize the production of these inhibitors. It is also important to acknowledge the limitations of relying on a single reference gene. Although actin is commonly used, its stability may fluctuate under different conditions (Ruan and Lai, 2007). Therefore, incorporating multiple candidate reference genes is preferable for identifying the most stable one under the experimental conditions, especially in detailed studies examining gene expression across various stages and environments.

To better understand the variability and functionality of BBI1 and BBI2 genes, we have analysed two datasets of DNA sequences including their homologs identified in 42 Vigna species. A phylogenetic analysis demonstrated that the amino acid residue P1 of the first domain, which is expected to confer inhibitory specificity for trypsin, is under positive selection in both genes. The interaction with protease is primarily determined by the P1 residue, and alterations to this residue significantly influence the inhibitor’s specificity (Clemente and Domoney, 2006). In general, positive selection promotes genetic variants that increase the organism’s fitness, favouring changes in amino acid residues to optimize adaptation to environmental pressures. At the P1 residue of the first domain, in addition to the amino acids arginine and lysine which commonly target trypsin, we identified alanine (in BBI1) which emerged in some African species, as well as histidine and glutamic acid (in BBI2) found in V. racemosa and V. stipulacea, respectively. According to the literature, when alanine is at P1, the BBI is known to inhibit elastase, a protease associated with various medical conditions, including inflammatory processes and their corresponding immune responses (Kawabata et al., 2002; Polverino et al., 2017; Ferreira et al., 2019). Evidence of elastase inhibitory activity, potentially linked to the BBI family, has also been recently reported in cowpea and it was suggested that natural inhibitors targeting elastase could hold potential for medical trials and therapeutic applications (Ferreira et al., 2019). The interaction between alanine and elastase has been widely predicted in the literature and confirmed in studies on other species (Clemente and Domoney, 2006; Rocco et al., 2011; De Paola et al., 2012). Meanwhile, the molecular targets of the other two amino acids, glutamic acid and histidine, remain largely unknown and have been much less investigated. While some studies described glutamic acid as inactive in this context (Clemente and Domoney, 2006), we hypothesize that histidine exhibits affinity for trypsin due to its positive charge. Therefore, based on all these findings, we propose that the P1 residue of the first domain has undergone continuous changes throughout evolution and that the residues defining the interaction have attempted to optimize their affinity with different proteases. As confirmation of this, BBI1 and BBI2 sequenced from different species exhibited varying energy requirements, but interaction with chymotrypsin in the second domain of BBI2 proved to be generally favored (Figures 4–6). The fact that this latter interaction generally requires less binding energy is attributed to the presence of amino acid residues at the P1 position which confer a higher specificity toward chymotrypsin. However, two important aspects must be considered for a correct interpretation of our results. Firstly, amino acids variability observed in the P1 residue seem insufficient to fully explain varying energy requirements. Additional residues, particularly those within the binding loops, play a significant role in determining inhibitory specificity (Supplementary Figures S5, S6, S7 and S8). Unfortunately, our current findings do not permit quantification of their individual contributions. Therefore, more in-depth analyses are needed to precisely characterize these interactions and the specific roles of these residues. Secondly, the interaction values obtained in this study should be interpreted as a preference of a specific domain for a particular target. We cannot rule out the possibility that, in nature, the domain may also bind to targets other than the one defined by P1.

Finally, observations of affinity values along the phylogenetic trees did not reveal any evolutionary pattern within the Vigna genus. Although our sampling was primarily concentrated in Africa and Asia, the ability of BBI1 and BBI2 to interact with the proteases analyzed in this study does not appear to be correlated with particular phylogenetic groups. Given the relatively recent origin of the Vigna genus (likely about ten million years ago), this timeframe may have been too short to create lineages with particularly favored interaction activity. However, the modern protein variants observed across different species have exhibited a wide range of interaction values, with differences that are not negligible. The inclusion of many taxa in the sampling allowed to increase the capability of exploration and new protein variants were found. The wide availability of wild Vigna species is an important resource for African and Asian countries and the findings obtained may prove useful for breeding programs or genome editing efforts aimed at selecting more efficient cultivars. However, to fully exploit the potential of BBIs for the benefit of both agriculture and human health, a nuanced approach is needed, that takes into account the differential specificity and the application of advanced breeding and genetic engineering techniques to tailor BBI profiles. Indeed, the complex and seemingly paradoxical nature of BBIs, which act as antinutrients by inhibiting digestive proteases in humans and pests, but also possess health-promoting properties, poses a difficult challenge. Nevertheless, there is now a range of strategies that can be deployed to address this issue. For example, precision genetic editing technology, such as CRISPR/Cas9, and RNA interference (RNAi) technology has been successfully applied to soybeans, resulting in a substantial reduction in trypsin and chymotrypsin inhibitory activity and improving digestibility without adversely affecting other nutritional characteristics (Kim et al., 2024, 2025). Similar strategies could also be applied in Vigna. In order to enhance specific human health benefits, it is necessary to identify or genetically manipulate BBI variants that interact more strongly with disease-related human proteases, but which bind less strongly to human digestive proteases. The “double-headed” structure of BBIs, prevalent in legumes like Vigna, with two distinct inhibitory domains that simultaneously and independently inhibit two different protease molecules provides a powerful platform for engineering. Using site-directed mutagenesis, one inhibitory site could be modified to reduce its binding to human digest enzymes and the other site (or a different BBI variant) could be modified to enhance its binding to specific human proteases involved in disease pathways. Concurrently, conventional breeding programmes, augmented by molecular markers and advanced phenotyping, could identify existing germplasm for BBI variants that inherently exhibit a more favourable balance between lower antinutritional activity for humans and enhanced activity against specific insect proteases and/or proteases that are beneficial to humans.

In this context, our findings offer a useful starting point for addressing this challenge. We highlight that docking simulations provide a valuable preliminary assessment of protein-target interaction potential, though biological validation is still required. Docking simulations applied to new genes encoding specific proteins represent an inexpensive and efficient method for the preliminary evaluation of protein variants hidden in wild and domesticated accessions.

Our findings on the high expression of BBI1 and BBI2 genes in Vigna seeds align with previous hypotheses suggesting a triple role in plant defense, regulation of endogenous protease activity, and storage protein function. We acknowledge that our expression data, derived from a representative subset of species and potentially influenced by non-specificity of some primer, provides only a preliminary view; confirming these patterns across a broader phylogenetic spectrum is an essential next step. Despite these considerations, our comprehensive evolutionary analysis across 42 species identifies the P1 residue of the first domain under positive selection. However, we interpret this selection signature with caution, recognizing that non-adaptive forces, such as genetic drift during population bottlenecks from distinct domestication events on different continents, can also elevate dN/dS ratios, thus creating signatures that mimic positive selection. This selective pressure, whatever its origin, appears to have driven a dynamic functional diversification. Our exploratory in silico docking, which provided an initial assessment of this diversity, was based on human protease models. Future studies must therefore validate these interactions using ecologically relevant targets, such as modeled insect or fungal enzymes, to truly understand their co-evolutionary context. Ultimately, this rich reservoir of natural variation presents a valuable genetic toolkit for future agricultural applications. We therefore present our conclusions as a robust working hypothesis, where the compelling genetic and initial computational insights now require large-scale transcriptomic and direct biological validation to fully harness their potential.