Functional and Comparative Phylogenomics of Pathogen-Inhibiting Native Mexican Bacillus velezensis Across Regional and Global Scales Reveal Open Pan-Genome Architecture and Adaptive Divergence

1. Introduction

Bacillus velezensis is one of the most prominent plant growth–promoting rhizobacteria (PGPR) recognized by its biocontrol activity against a broad spectrum of phytopathogens, including fungi, bacteria, and nematodes [1,2]. As a Gram-positive species it is characterized by its remarkable ability to form endospores and by harboring an extensive repertoire of biosynthetic gene clusters responsible for the production of antimicrobial secondary metabolites, which enable the inhibition of competing microorganisms [2,3,4]. In addition, B. velezensis produces molecules with phytohormonal activity, such as indole-3-acetic acid (IAA) [5,6], and cytokinins [7,8], among others, that facilitate close interaction with the host plant and exert significant effects on key physiological processes including cell division regulation, thereby promoting plant growth [7]. Furthermore, this species plays a relevant role in modulating plant immunity, notably through the activation of induced systemic resistance (ISR) which strengthens the plant’s defense responses against phytopathogens [9,10].

Beyond its biocontrol capabilities, B. velezensis also functions as a bio-organic fertilizer due to its ability to solubilize phosphorus [11,12] and potassium [13,14], enhance sulfur assimilation [15,16] and mediate iron chelation and uptake [17,18]. Numerous studies highlight Bacillus spp., and particularly B. velezensis, as effective rhizobacteria against fungal phytopathogens affecting plant roots, as well as bacterial pathogens impacting both root systems and foliar tissues [19,20]. Some strains have already been developed into commercial agrobiotechnilogical formulations, including RhizoPlus^®, RhizoVital^® (Andermatt Group AG, Suitzerland), Taegro^® (Novozymes, México, SA. de CV) Botrybel^® (Probelte, España), Fungifree AB^® (Agrobiotics, México), VELEZIGEN^® (Cergen SRL Bio Tecnología, Argentina), among others.

Whole sequencing of microorganisms, coupled with the rapidly expanding number of genome records available in public databases, has opened extensive opportunities for comprehensive and in-depth studies. These approaches enable the characterization of the full gene repertoire, or pangenome, that defines a given phylogenetic clade, as well as the identification of the core genome shared across all members of a microbial population. This core set comprises essential genes required for fundamental metabolic functions that underpin key phenotypic traits and are critical for microbial survival within their predominant ecological niche [21,22].

Isolates of B. velezensis have been reported across all continents, with the highest number of identified strains documented in Asia, particularly in China, where the largest number of genomes has been sequenced. Species within the genus Bacillus typically possess genomes ranging from approximately 3.73 to 5.79 Mb, with B. velezensis exhibiting an average genome size of around 3.99 Mb [23,24].

In north-central Mexico, particularly in the state of Zacatecas, rhizobacteria have been isolated for the inhibition and biocontrol of root pathogenic fungi in chili pepper plants, as well as root and foliar pathogens in common bean. Among these isolates, strains exhibiting the highest inhibitory potential against phytopathogens and plant growth-promoting traits have, for reasons yet unclear, consistently been identified as belonging to the species B. velezensis [25,26]. Four genomes from these strains have been sequenced, corresponding to strains 2A-2B, 3A-6A, 2A-10A and 3A-25B, with NCBI accession numbers NZ_MLCV00000000.1, MLCX00000000.1, MKWT00000000.1 and MLCW00000000.1 respectively.

The aim of the present study was to analyze the phenotypic variation among these four B. velezensis strains in terms of their ability to inhibit phytopathogens and their interaction with chili pepper plants following root inoculation, in relation to their genomic and genetic features directly or indirectly associated with plant interaction and plant growth promotion. Additionally, the study evaluated the flexibility of the pan-genome, core genome, and singleton gene content at regional, continental, and global scales.

2. Materials and Methods

2.1. Dual Culture Assay Between Bacterial Strains and Fungal or Oomycete Phytopathogens

Virulent strains of the pathogens Phytophthora capsici, Rhizoctonia solani, Fusarium solani, and Fusarium oxysporum on chili pepper plants of the cultivar Mirasol were cultured on PDA or KB agar plates; fungi were grown for 4 days and the oomycete for 8 days. Bacillus velezensis strains 2A-2B, 3A-6A, 2A-10A and 3A-25B were grown in liquid KB medium for 12 to 18 h and subsequently subcultured for an additional 40 to 8 h until reaching a OD₆₀₀ of 2.0, equivalent to 2x10⁷ UFC/mL. After confirming that both that fungi and the oomycete grew adequately on KB medium, dual confrontation assays between bacterial strains and pathogens were performed.

For the dual confrontation assay, 8 mL of bacterial suspension (2 × 10⁷ CFU/mL) obtained from 4-6 h subcultures, were inoculated at four equidistant points on KB agar plates, positioned 5 mm from the edge of the Petri dish. Subsequently, each plate was centrally inoculated with the corresponding fungal or oomycete pathogen by placing a 3.5 mm diameter mycelial plug, obtained with a cork borer from actively growing cultures previously established on PDA medium. As a control treatment, each pathogen was grown individually from the center of KB plates without bacterial confrontation. Plates were incubated at 28ºC in darkness for 5-12 days. Four biological replicates were included per treatment.

Pathogen growth and inhibition were measured every 3 days up to 12 days. The percentage of inhibition was calculated using the following formula:

$$\mathrm{Pathogen}\mathrm{inhibition}={\displaystyle \frac{\mathit{r}\mathit{c}-\mathit{r}\mathit{t}}{\mathit{r}\mathit{c}}}\mathit{x}100$$

where rc is the radial growth in the control treatment and rt is the radial growth in the respective treatment. Differences in inhibition percentages for each pathogen in confrontation with each bacterial strain were analyzed. Statistical analyses were performed using PSPP software version 2.0.1. [27].

For these bacterial strains, three genomes had been previously sequenced, whereas the fourth was sequenced more recently using illumina (MiSeq) technology. Comparative genomics and subsequent bioinformatic analyses indicated that all four strains are consistently assigned to the species B. velezensis [25,26]. Three of the genomes had been previously deposited in the NCBI database, while the fourth genome was submitted more recently.

2.2. In Vitro Confrontation of B. velezensis Bacteria with Chili Pepper Plants cv. Mirasol

Fifteen-day-old chili pepper plants (cv. Mirasol) grown in vitro on MS medium in Petri dishes were root-inoculated with each B. velezensis strain independently, using a liquid suspension at 2x10⁷ CFU/mL. Inoculation was performed by inmmersing the roots in the bacterial suspension for 10 min with agitation every 3 min. After inoculation, the plants were transferred to fresh MS medium and maintained at 37ºC under a 12-h light/12-h dark photoperiod for 8 days of observation.

The experiment consisted of six treatments with six replicates each. T1 = Control (plants without bacterial inoculation), T2 = Inoculation with strain S50-2A-1 (positive controls, hostile strain), T3 - T6 = Inoculation with strains 2A-2B, 2A-10A, 3A-25B and 3A-6A, respectively. Plant damage was assessed using a five-level scale, based on Sunwoo et al. [28] with minor modifications: 0 = no visible disease symptom (healthy plant); 1 = initial brownish lesions on root and stems, with slight leaf wilting, 2 = necrosis affecting the main root system and stem base, with 30-50% of the plant affected, 3 = necrosis affecting the main root system and stem base, with 50-70% oh the plant affected, and 4 = necrosis affecting the main root system and stem base, with ≥ 70% of the plant affected. Bootstrap analysis was performed using the boot package in R, with a 95% confidence interval.

2.3. Review and Analysis of the Genomic Features of Native B. velezensis Strains 2A-2B, 3A-6A, 2A-10A and 3A-25B.

Genome assembly quality was assessed using QUAST algorithm, while genes and coding sequences were annotated with the Prokka algorithm on the Galaxy server. This analysis allowed us to confirm, quantify, and summarize the basic genomic features, including genome size, GC content, number of coding sequences (CDs), RNA genes, and putative pseudogenes. The genomes of the four strains (2A-2B, 3A-6A, 2A-10A and 3A-25B) under study have been registered in the NCBI database and assigned accession numbers NZ_MLCV00000000.1, MLCX00000000.1, MKWT00000000.1 and MLCW00000000.1 respectively. (The accession numbers for the genomes of strains 3A-6A, 2A-10A in the NCBI database will be released next week).

2.4. Comparative Genomic and Functional Analyses of Four Native B. velezensis Strains from North-Central Mexico, Including Reference Genome

For the analyses of the pan-genome, core genome, dispensable genome and singleton gene of the four native B. velezensis strains (2A-10A, 3A-6A, 2A-2B, and 3A-25B) the genome of strain JS25R (accession NZ_CP009679) was included as a reference. All the genomic datasets were analyzed using EDGAR platform [29]. Gene prediction and coding sequence (CDS) annotation were standarized within the platform to ensure comparability. The pan-genome, core genome, and singleton functions were computed using the corresponding EDGAR modules, and their proportions were determined based on the total number of gene families identified across the five genomes. Pan-genome development curves and core genome reduction trends were generated through iterative genome addition analysis. Gene cardinality and the distribution of shared and unique CDS among strains were assessed using the Geneset tools, including Venn diagram analyses.

Genomic similarity among strains was evaluated through Average Nucleotide Identity (ANI) analysis using FastANI algorithm implemented in EDGAR environment. Pairwise ANI values were calculated for all genome combinations to determine intra- and interspecific relationships. Additionally, the Porcentage of Conserved Proteins (POCP) was computed to assess evolutionary relatedness and funcional conservation at the proteome level.

Functional annotation and classification of genes were performed using both KEGG and COG framworks. KEGG-based analyses categorized genes into metabolic pathways, genetic information processing, environmental information processing, and cellular processes, using KEGG Sunburst tool [30]. In parallel, Clusters of Orthologous Groups (COG) classification was applied to assign genes into functional categories such as motility, secondary metabolite biosynthesis, and defense mechanisms [31]. The distribution and abundance of genes within each functional category were quantified for each genome and comparatively analyzed across strains.

All analyses, including pan-genome modeling, ANI, POCP, functional categorization, and gene distribution assessements, were performed under default settings unless otherwise specified.

2.5. Multi-Scale Comparative Pan-Genome and Phylogenomics Analyses of B. velezensis Across Regional, Continental, and Global Populations

To evaluate pan-genome structure, genomic diversity, and phylogenetic relationships of B. velezensis across increasing geographic scales, a hierarchical comparative genomics approach was implemented using EDGAR platform. Genome datasets were structured into four levels: (i) regional (north-central Mexico, including four native strains plus one reference genome), (ii) North America (19 genomes from México, the United States, and Canada, (iii) intercontinental comparison (North America vs. South America; 9 genomes each), and (iv) global scale (30 genomes from Americas, Eeurope, Asia and Africa). Genome selection was based on availability of complete genomes in the NCBI database and geographic representativeness.

All genomes were processed within EDGAR under standarized annotation pipelines to ensure consistency in coding sequences (CDS) prediction and ortholog identification. Pan-genome, core genome, dispensable genome, and singleton fractions were computed using EDGAR modules. Pan-genome development and core genome decay were assessed through sequential genome addition analyses, allowing evaluation of genome openness on Heaps’law [32]. Gene distribution patterns and cadinality were analyzed to determine shared and strain-specific genomic content across datasets.

Comparative analyses between geographic regions were conducted by normalizing genome samples sizes where necessary (e. g., equal number of genomes for North and South America) to reduce sampling bias. Pan-genome structure and dynamics were evaluated independently for each region and comparatively interpreted across scales.

Phylogenetic relationships were inferred usin EDGAR phylogenomics tools based on conserved gene sets. Phylogenetic trees were constructed to assess clustering patterns, lineage divergence using ANI and POCP, as implemented in EDGAR. These metrics provided quantitative support for intra-species relationships and protein-level conservation across strains.

Finally, comparative analyses of pan-genome structure across scales (regional, continental, and global) were performed by integrating absolute gene counts and relative proportions of core, dispensable, and singleton fractions. Observed patterns were interpreted in the context of genome plasticity, geographic structure, and evolutionary processes sucha as gain and horizontal gene transfer. Default parameters were used in all EDGAR analyses unless otherwise specified.

2.6. Analysis of Bacterial Gene Groups Associated with Plant-Microbe Interaction and Plant Growth Promotion

To analyse genes involved in plant-bacteria interactions and plant growth promotion, the PLaBAse platform wass used. Within PLaBAse the PGPT-Pred tool was employed to predict plant growth-promotion traits, enabling the identification of gene groups organized into eight functional categories: colonizing planr system, competitive exclusion, stress control/biocontrol, biofertilization, phytohormone and plant signal production, bioremediation, plant immune response stimulation, and putative functions [33,34].

For the analysis of genes associated with plant-bacteria interactions, the PIFAR-Pred tool was utilized. This tool classifies genes into 13 functional categories, including exopolisaccharide synthesis, hormone biosynthesis, antibiotic biosynthesis, detoxification, motility, volatile compound production, adhesion, plant cell wall-degrading enzymes (PWDE), proteases, among others [34].

Both PGPT_Pred and PIFAR-Pred operate through the identification of marker genes and are supported by a database comprising approximstely 6,900 plant growth-promoting ttraits (PGPTs) associated with 6,965,955 protein sequences [33,34].

3. Results

3.1. Phenotypic Variation Among Bacterial Strains in Phytopathogen Inhibition and Root Interaction in Chili pepper cv. Mirasol

In a screening of 150 rhizobacterial isolates obtained from wild plants in forest soils and cultivated plants in the state of Zacatecas, Mexico, 21 strains exhibited ≥35% inhibition of each of the pathogens Rhizoctonia solani, Phytophthora capsici, Fusarium solani, and Fusarium oxysporum in in vitro dual confrontation assays. Among the most effective rhizobacteria, isolates 2A-2B, 3A-6A, 2A-10A and 3A-25B consistently stood out, showing inhibition levels ranging from 38% to 83% against all four phytopathogens. (Figure 1 and Figure 2). In the comparative phytopathogen inhibition analysis, strain 3A-25B inhibited R. solani in a 37.8%, while a 40.8%, 44.2% and 45.5% inhibition was achieved with 3A-6A, 2A-10A and 2A-2B strains respectively, however, it exhibited a major inhibiting effect over P. capsici than 3A-6A strain, with an 81.2 and 69.7% respectively, and of similar manner over F. solani, with a 67.9 and 61.7% respectively. On the other hand, strain 2A-10A briefly surpasses the other 3 strains in pathogens P. capsici and F. solani inhibition. No statistically significant differences were observed in any case (Figure 1 and Figure 2).

In root inoculation assays evaluating bacterial patogenicity in chili pepper plants under in vitro conditions, observations at 8 days post-inoculation showed that bootstrap analysis (10,000 iterations) identified the treatment inoculated with strain S50-2A-1 (T2), corresponding to the positive control (pathogenic bacterium causing necrosis in the primry and lateral roots), as having a disease index of 87.5%. This value was significantly higher than those of the remaining treatments, which ranged from 6.25% to 18:75%. Confidence imtervals showed substantial overlap among these latter treatments, indicating no significant differences between them. Bootstrap analysis was performed using the boot package in R, with 95% confidence intervals (Figure 3).

3.2. Genomic Features of Four Bacillus velezensis Strains Native to North-Central Mexico

Genome sequencing achieved approximately 99% coverage of the total genome length. Genome sizes ranged from 3.93 to 4.01 Mb, values very close to the reported average for this species (~3.9 Mb). GC content, gene number, and coding sequences were higly similar across the four strains. However, variation was observed in the number of rRNAs, tRNAs, miscellaneous RNAs (miscRNAs), and tmRNAs. Strain 3A-25B exhibited the greater divergence, particularly in miscRNAs and tmRNAs, where reduced counts were detected, whereas the other three strains showed high similarity in these features (Table 1).

3.3. Pan-Genome of Four B. velezensis Strains from North-Central Mezico with a Reference Genome

The pan-genome of the four B. velezensis strains native to north-central México, including strain JS25R (accession NZ_CP009679) as the reference genome, comprises 4,258 genes. Of these, 77.9% correspond to the core genome (3,345 genes), 14.4% to the dispensable genome (582 genes), and 7.7% to singleton genes (331 genes) (Figure 4-A). In the pan-genome development, the sequential addition of genomes showed that the inclusion of strain 3A-25B resulted in a marked increase in pan-genome size, while the core genome exhibited a gradual decline (Figure 3-B). The core genome encompasses higly conserved genes present in most of the analysed genomes, presumably associated with essential cellular functions, including central metabolic processes, maintenance and processing of genetic information, and other components of the fundamental bacterial cellular machinery.

In terms of gene cardinality, the five genomes displayed a relatively balanced number of coding sequences (CDS), forming a conserved core genome composed of 3,337 CDS shared among all strains. Each genome contributed a variable number of unique genes, with strain 3A-25B being the most prominent, harboring 149 singleton genes, exceeding the combined total of the three strains (138 singleton genes) (Figure 4-C).

Genomic similarity among these bacterial strains was assessed using Average Nucleotide Identity (ANI) calculated with the FastANI allgorithm, which is considered accurate for both complete and draft genomes. Consistent with established tresholds, ANI values >95% indicate intraspecific relationships, wheras values <83% suggest interspecific divergence [35]. The analysis revealed that strain 3A-25B exhibits the highest similarity to strain JS25R, with an ANI of 99.15%, while ANI values between 3A-25B and the other three native strains ranged from 97.8% to 97.9%. In contrast, ANI values among strains 2A-2B, 2A-10A and 3A-6A ranged between 99.6% and 99.8% (Figure 5-B).

To further assess evolutionary relatedness and potential phenotipic differences in terms of protein conservation, ultimately responsible for cellular funcitonal processes, the Percentage of Conserved Proteins (POCP) was calculated [36]. POCP analysis performed on EDGAR server showed that strain 3A-25B is evolutionarily closer to both the reference strain JS25R and native strain 2A-2B, with values of 95.6% and 95.5%, respectively, whereas POCP values among the other three native strains ranged from 96.9% to 98.8% (Figure 5-C).

Functional characterization based on metabolic pathways, signaling, and other biological processes using KEGG, together with lineage and protein function–based analysis using COG, revealed that the genomes of the four B. velezensis strains, in comparison with the reference genome, display highly similar functional profiles. According to KEGG classification, the most abundant gene categories correspond to genetic information processing (534–551 genes) and signaling and cellular processes (559–574 genes), whereas the least represented categories include terpenoid and polyketide metabolism, as well as organismal systems. Notably, a substantial number of genes remain unannotated, ranging from 1,353 to 1,472, predominantly in strain 2A-2B and, to a lesser extent, in strain 3A-6A (Figure 6). Comparative analysis of gene counts across functional categories among the four native strains, as well as among all five strains including the reference, indicates an overall high degree of similarity (Figure 6).

Based on Clusters of Orthologous Groups (COG) classification, these bacteria possess 46-48 genes associated with motility, an important trait for chemotaxis during rhizosphere colonization. They also harbor 87-91 genes involved in secondary metabolite biosynthesis and 105-113 genes related to defense mechanisms, both relevant for antibiosis and inhibition of competing microorganisms. Additionally, the native B. velezensis strains from Mexico contain beween 28 and 29 transposase genes, in contrast to the reference strain JS25R, which contains only 15. As observed with KEGG, COG analysis also revealed a high number of unannotated genes, ranging from 609 to 684, again predominantly in strain 2A-2B and secondarily in strain 3A-6A (Figure 6, Table 2).

3.4. Pan-Genome and Phylogenetic Placement of Four B. velezensis Strains Native to North-Central Mexico at the North American Scale

Pan-genome analysis of the native strains fron north-central Mexico, expanded to include 19 B. velezensis strains fron North America, comprising 5 strains from Canada, 6 from United States, and 9 from Mexico (the four under study plus 5 additional strains) with strain FZB24 as he reference genome, revealed a total pan-genome of 5,759 genes. Of these, 2,584 correspond to the core genome (44.9%), 2,183 to the dispensable genome (37.9%), and 992 to singleton genes (17.2%) (Figure 7-A). The high proportion of accesory and singleton genes (55.1%) indicates substantial genomic plasticity in this B. velezensis species, at least within the northern region of the American continent.

Pan-genome development analysis showed that the inclusion of the ninth genome, originating from Montreal, Canada, resulted in a marked reduction of the core genome from 3,197 to 2,699 genes. Thereafter, the core genome remained relatively stable with the addition of the remaining genomes, ultimately reaching 2,583 genes. This suggests that approximately 2,580-2,700 core genes are sufficient for these bacteria to maintain standard growth under natural conditions. In contrast, the pan-genome exhibited a gradual upward trend with the sequential addition of genomes. According to Heaps’law, the continuous increase in pan-genome size indicates an exponent α < 1, characteristic of an open pan-genome [37] (Figura 6-B).

Phylogenetic analysis based on the genomes of the same 19 bacterial strains revealed that the native Mexican strains do not form a single clade. While 2A-2B, 3A-6A, and 2A-10A cluster closely together, indicating high genomic similarity, strain 3A-25B is positioned on a separate phylogenetic branch (Figure 8).

3.5. Comparative Distribution, Pan-Genome Dynamics, and Phylogeny of B. velezensis Strains Across North and South America

In this comparative analysis, it should be noted that the number of complete genomes available in the NCBI database for B. velezensis from South America is limited. Therefore, nine accessions, mostly from Brazil, were included, and, for balance, nine genomes from North America were included, comprising three from Mexico (including strain 2A-2B), three from the United States and three from Canada.

The North American group exhibits a predominant core genome representing 58.6% of the pan-genome (2,762 CDS), along with smaller fractions of dispensable genes (27.2%) and singletons (14.2%). The pan-genome accumulation curve shows a continuous increase without reaching a plateau, while the core genome gradually declines and stabilizes at approximately 2,700 CDS (Figure 9-A). This pattern indicates an open pan-genome alongside a high degree of genomic conservation, consistent with a relatively homogeneous population.

In contrast, South American strains display a higher proportion of variable (dispensable) genes. The core genome accounts for only 32.3% of the pan-genome (1,615 CDS), whereas dispensable genes represent 50.3% and singletons 17.5%, with the dispensable fraction being the largest component. As additional genomes are incorporated, the pan-genome shows sustained growth, accompanied by a pronounced reduction in the core genome, which stabilizes at approximately 1,600 CDS (Figure 9-B), also consistent with an open pan-genome.

In the phylogenetic analysis of strains from the South American region, including the four native strains from north-central Mexico, strains 2A-10A, 3A-6A and 2A-2B from this region clustered together, whereas strain 3A-25B is positioned in a separate subclade alongside Brazilian strains LABIM22, FTCO1 and 3TSA-3 (Figure 10).

3.6. Comparative Analysis of B. velezensis Pan-Genome Structure Across Increasing Levels of Geographic Sampling and Population Size

Considering 30 genomes of B. velezensis strains at a global scale, originating from the American and Afro-Eurasian continents, including the four strains from north-central Mexico, the pan-genome comprises 6,535 genes. Of these, the core genome consists of 2,462 genes (37.7%), the dispensable genome includes 2,614 genes (40.0%), and the singleton fraction accounts for 1,459 genes (22.3%). These results indicate marked intraspecific genetic diversity (Figure 11-A).

In pan-genome development, the initial gene set comprises 3,187 genes. The core genome shows a sharp decline upon incorporation of the third genome, followed by a more gradual decrease from the fourth genome onward, as each of the remaining 27 genomes is sequentially added, ultimately stabilizing at 2,462 genes. In contrast, the pan-genome trajectory exhibits several modest yet noticeable increments, reaching a final total of 6,535 genes. Overall, the progressive addition of genomes results in a continuous increase in the total number of genes without approaching saturation, while the core genome gradually declines and stabilizes at approximately 2,500–2,700 genes (Figure 11-B). This pattern is indicative of an open pan-genome, in which each newly incorporated genome contributes additional genes.

A comparative assessment of B. velezensis pan-genome structure was conducted across hierarchical levels of geographic sampling and population size, ranging from lower to higher complexity. The lowest level (4CNM+1ref) represents north-central Mexico (5 genomes), intermediate levels correspond to North America (14 genomes) and South America (11 genomes) analyzed separately, and the highest level represents the global dataset (30 genomes). Both absolute gene counts and relative proportions of genomic fractions were evaluated (Figure 12).

In terms of absolute gene numbers (Figure 12, left panel), a progressive increase is observed as both the number and geographic diversity of genomes expand. The 4CNM+1ref dataset exhibits the smallest pan-genome, whereas North and South America show intermediate expansion, and the global dataset displays the largest pan-genome. This pattern is characteristic of an open pan-genome, in which the incorporation of additional genomes continuously contributes new genes. The core genome correspondingly decreases as sampling expands, a predictable mathematical effect whereby increasing genome numbers reduce the conserved fraction. In contrast, the dispensable genome and singleton fractions increase with broader geographic scale, indicating greater genetic diversity and potential regional specialization.

In relative terms (Figure 12, right panel), the pan-genome is normalized to 100% as a reference. The core genome proportion decreases from regional to continental and global scales, whereas the dispensable and singleton fractions increase, particularly at the global level and in the ProPan database (prokaryotic pan-genome database) for B. velezensis species. This trend reflects a lower proportion of conserved genes and a higher contribution of accessory and unique genes. According to ProPan post-processed datasets [34], the reported distribution for B. velezensis consists of 8.13% core genome, 55.36% dispensable genome, and 36.51% singletons. The results obtained in the present study deviate somewhat from these values; however, the South American dataset shows the closest correspondence to the ProPan distribution (Figure 12).

Overall, the regional, continental, and global pan-genomes of B. velezensis analyzed in this study consistently exhibit an open structure, with genetic diversity strongly influenced by geographic structuring, high evolutionary plasticity, and expansion of accessory genome associated with ecological niche adaptation.

With respect to the global evolutionary history of the four B. velezensis strains native to north-central Mexico, phylogenetic relationships inferred from a dataset of 30 representative strains across continents, including 13 from North America (among them the four Mexican strains), 8 from South America, 4 from Asia, 3 from Europe, and 2 from Africa, indicate that strains 2A-10A, 3A-6A and 2A-2B consistently cluster together, closely related to a strain from Brazil. In contrast, strain 3A-25B is more distantly positioned and shows closer evolutionary affinity to strains from Europe and Asia, specifically from Italy and China (Figure 13).

This pattern may be associated with genomic differences, particularly the presence of a transposase gene unique to strain 3A-25B. A comparative search for transposase genes revealed that all four strains share up to six transposase genes, including types such as ISAau1, ISSsu9, ISCac2, and ISBsu1, among other variants. However, strain 3A-25B harbors a transposase gene associated with the Tn916 transposon, which is absent in the other three strains. Tn916 belongs to a group of mobile genetic elements capable of integrating into the genome, mediating conjugative transfer to recipient bacterial cells, and undergoing intracellular transposition to different genomic loci. Additionally, these elements can readily acquire accessory genes and promote genome rearrangements [34]. This distinctive feature of strain 3A-25B may underlie genome restructuring over time, contributing to its phylogenetic divergence from the other native B. velezensis strains.

3.7. Bacterial Gene Groups Associated with Microbe-Plant Interaction and Plant Growth Promotion

For bacteria intended for biocontrol of phytopathogens, it is essential to assess their genomic composition in terms of gene involved in plant interaction and plant growth promotion. A useful resource for this purpose is the PLaBAse platform, which provides tools such as PGPT-Pred for the prediction of plant growth-promoting traits (PGPTs) and PIFAR-Pred for the identification of genes involved in plant-bacteria associations. Both tools rely on marker gene identification and are supported by a database comprising approximately 6,900 PGPTs associated with 6,965,955 protein sequences [33,34].

Key biological functions involved in interaction with plant tissues in the rhizosphere include chemotaxis, adhesion, synthesis of exopolysaccharides and lipopolysaccharides, plant cell wall-degerading enzymes, hormone biosynthesis, and volatile compound production, among others. Analysis using PIFAR-Pred tool across the four genomes under study revealed a high degree of similarity among these B. velezensis strains. All strains possess gene repertoires associated with exopolysaccharide synthesis, hormone biosynthesis, and antibiotic production, incluing compounds such as syringomycin, syringopeptin, toxoflavin, and tagetitoxin, as well as genes related to motility, adhesión, and detoxification.

The main differences among the four strains were observed in toxin-related categories. Notably, strain 2A-10A exhibited a slight increase (~1%) in genes associated with biosynthesis of three antibiotic types. Strain 3A-25B also showed moderate enrichment in genes related to metabolism, pigment production, motility and antibiotic biosynthesis (Figure 14).

To assess similarities and differences in plant interaction-related bacterial factors between these strains and those form other geographic regions, strain JS25R was included as reference genome; this strain originates from Yancheng province, Jiamgsu, China. The results indicate that the native Mexican strains exhibit broadly similar profile of plant interaction factors compared to the reference strain. However, differences were observed in specific functional categories: the reference strain showed reduced representation of genes associated with metabolism and plant cell-wall degrading enzymes (PCWDE), whereas it displayed comparatively higher proportions of genes related to exopolysaccharide production and toxin synthesis.

Additionally, native strain 3A-25B showed greater similarity to the reference strain JS25R than to other three native strains, particularly in the proportion of genes associated with exopolysaccharide production, metabolism, and motility (Figure 15).

With respect to plant growth-promoting traits (PGPTs), analysis of the four B. velezensis strains from north-central Mexico revealed a high degree of similarity among them. An exception was strain 3A-25B, which showed a higher number of genes associated with bioremediation and biofertilization, but a lower representation of genes related to competitive exclusion (Figure 16). Notably, all four strains harbor the rbcL and cbbR genes, which are associated with the regulation of CO₂ fixation through alternative pathways or those related to RUBISCO.

Within the biofertilization category, the strains possess approximately 40 genes involved in iron acquisition, 51 in nitrogen acquisition, more than 100 in phosphate solubilization, fewer than 100 in potassium solubilization, and around 14 in sulfur assimilation. Gene-level comparisons indicate the highest similarity between strains 2A-2B and 3A-6A. In the phytohormone signaling category, approximately 40 genes are associated with indole-3-acetic acid (IAA), 7 with cytokinins, and smaller numbers with abscisic acid (ABA) and gibberellins, along with 13 genes related to acetoin/butanediol production. Approximately 825 genes are involved in plant tissue colonization, with slightly higher representation in strains 2A-2B and 2A-10A, although gene-level similarity is again highest between strains 2A-2B and 3A-6A.

In competitive exclusion category, around 630 genes were identified, with a slightly lower number in strain 3A-25B; gene-level similarity was again highest between strains 2A-2B y 3A-6A. For plant immune response induction, 66-67 genes were identified, and for stress control/biocontrol, between 585 and 599 genes, with the highest similarity consistently observed between strains 2A-2B and 3A-6A, including three genes with insecticidal function.

Across the seven functional PGPT categories, the total number of genes identified was 4,095, 4,070, 4,043, and 4,111 for strains 2A-2B, 2A-10A, 3A-6A and 3A-25B, respectively. These values are slightly higher than the total gene counts obtained from genome annotation, likely because annotation relies on predictive approaches that do not always capture all functional elements and are primarily based on coding sequences (CDS).

To assess the degree of similarity or divergence among native strains from north-central Mexico in the proportion of plant growth-promotion traits, both among themselves and relative to strains from other geographic regions, strain JS25R fron China was included as a reference. The proportional distribution of genes associated with plant growth-promoting functions is highly consistent among the ntive Mexican strains. Again, strain 3A-25B represents an exception, exhibiting a reduced proportion of genes associated with competitive exclusion (21% vs. 22%) and a slightly higher proportion associated with bioremediation (7% vs. 6%).

In comparison with the reference strain, the proportions of genes related to bioremediation, plant colonization, biocontrol, and stress response are identical, at least for strains 2A-2B, 2A-10A and 3A-6A. However, strain 3A-25B shows a 1% reduction in the competitive exclusion category, while simultaneously matching the reference strain with a 1% higher proportion in the bioremediation category (Figure 17).

3.8. Strain-Specific Genes (Singletons) in B. velezensis Across Regional and Global Scales

The presence of singleton genes was analyzed in bacterial genomes across regional (North America and South America) and global scales. In the North American comparison, which included 20 genomes from Mexico, the United States, and Canada, the four strains from north-central Mexico collectively harbored 70 singleton genes distributed as follows: 40 in strain 3A-25B, 12 in strain 2A-2B, 10 in strain 3A-6A, and 8 in strain 2A-10A. In strains 2A-2B, 3A-6A and 2A-10A, all singleton genes were annotated as hypothetical proteins. In contrast, among the 40 singleton genes in strain 3A-25B, 8 correspond to proteins with known functions and 32 to hypothetical proteins; among the characterized genes, three are associated with phage structure or recombination. As a reference genome, strain FZB42, native to Wisconsin, USA, was included; this strain contains only three singleton genes, all encoding hypothetical proteins (Table 2).

For South America, where genome availability in public databases is more limited, a total of 13 genomes were included, In this comparison, the four strains from north-central Mexico collectively accounted for 99 singleton genes, with the majority (72) found in strain 3A-25B, followed by 10 in strain 2A-2B, 10 in strain 3A-6A, and 7 in strain 2A-10A. In strains 2A-2B, 3A-6A and 2A-10A, all singleton genes were hypothetical proteins, whereas in strain 3A-25B, 50 were hypothetical and the remainder corresponded to proteins with known functions, including roles in phage-related processes, manganese transport, and DNA packaging. The reference genome for this dataset was strain TrigoCor1448, native to Brazil, which contains 13 singleton genes, of which 5 encode hypotheticl proteins and 8 have known biological functions. Among these are members of the TetR/AcrR family, which act as transcriptional regulators involved in gene expression associated with antibiotic resistance and responses to environmental stresses such as osmotic and oxidative stress [38].

At the global scale, based onthe analysis of 30 genomes distributed across North America, South America, eurasia, and Africa, the four strains from north-central Mexico collectively harbor 51 singleton genes: 5 in strain 2A-2B, 7 in strain 2A-10A, 10 in strain 3A-6A, and 29 in strain 3A-25B. Most are annotated as hypothetical proteins, except for five genes with known functions identified in strain 3A-25B, including genes involved in phage structural assembly, viral genome packaging into capsids, autolysin activity, and DNA packaging (Table 2). The reference genome included at this scale was strain JS25R originating from Jiangsu Province, China, which contains 13 singleton genes, 4 encoding hypothetical proteins and 9 with known biological functions. Among these is a CbrC family protein, associated with bacterial tolerance to colicin E2, an antibiotic-like toxin [39], among other functions.

4. Discussion

The four B. velezensis strains (3A-25B, 3A-6A, 2A-10A and 2A-2B) native to north-central Mexico consistently inhibited the four major chili root pathogens, with inhibition percentages ranging from 38% to 80%. Rhizoctonia solani was the least affected pathogen, whereas Phytophthora capsici was the most susceptible, with no statistically significant differences among bacterial strains (Figure 1 and Figure 2). These values are consistent with reports for B. velezensis strains from the Amazon region, where inhibition levels of 73–100% have been observed against phytopathogenic fungi such as Botrytis cinerea, Colletotrichum acutatum, and Penicillium expansum [40], as well as strains from Morocco showing inhibition of 67–78.4% against Phytophthora infestans [41]. This suggests that core antifungal mechanisms in B. velezensis may be conserved across geographically distant lineages.

In plant interaction assays conducted in vitro, strains 2A-2B, 2A-10A, 3A-6A, and 3A-25B exhibited low pathogenic impact, with the highest mean disease index reaching only 18.7%, in contrast to the pathogenic strain S50-2A-1, which reached 87% (Figure 3). These findings align with numerous studies demonstrating that B. velezensis strains are generally plant-friendly and capable of promoting plant growth [42,43].

At the genomic level, the four strains showed high similarity in genome size, GC content, gene number, and coding sequences (CDS), but differed in non-coding RNA elements, including rRNAs, tRNAs, miscRNAs, and tmRNAs. One strain exhibited nearly half the number of tRNAs compared to the others, while another contained only five miscRNAs and lacked tmRNAs enterily (Table 1). These differences in non-coding RNAs do not appear to directly influence antifungal activity, as all four strains remained statistically equivalent in this trait.

Pan-genome analysis of the four native strains plus reference strain JS25R revealed that strain 3A-25B harbors 149 singleton genes, exceeding both the 110 singletons in strain JS25R and the 22-25 singletons found in strains 3A-6A, 2A-10A and 2A-2B. Approximately 80% of shared genes constitute the core genome, with 14% classified as dispensable and 8% as singletons, indicating a moderately variable pan-genome with a large conserved functional backbone. Comparative genomic analyses based on average nucleotide identity (ANI) and percentage of conserved proteins (POCP) showed that strain 3A-25B is more divergent from the other native strains and clusters more closely with the reference strain JS25R from China. This pattern suggests non-uniform divergence among native strains, potentially reflecting distinct evolutionary trajectories, differential gene acquisition, or similar ecological selection pressures acting on strains 3A-25B and JS25R.

Functional annotation based on KEGG and COG categories revealed highly similar gene distributions among the four Mexican strains and the reference strain JS25R. The most represented categories corresponded to genetic information processing, signaling, and cellular processes, whereas terpene/polyketide metabolism and organismal systems were the least represented. Notably, approximately 17% of genes remained unannotated (~1,400 genes) (Figure 6), indicating a substantial proportion of “functional dark matter” with potential biotechnological relevance [44]. COG analysis further showed enrichment in amino acid transport, transcription, cell wall/membrane biogenesis, signal transduction, and general functional categories, alongside a similarly high proportion of uncharacterized genes (610–680) (Figure 6).

At the North American scale (19 genomes), the pan-genome exhibitedd continuous increase, while the core genome declinend notably after incorporation of the ninth genome (from Montreal, Canada), consistent with an open pan-genome structure. Phylogenetic analysis revealed that strains 2A-10A, 3A-6A and -2A-2B cluster together, whereas strain 3A-25B occupies a more distant subclade alongside other Mexican strains (Figure 7). This configuration underscores genetic differentiation among bacteria sharing a common geographic origin and suggests partially independent evolutionary trajectories, potentially associated with distinct rhizosphere niches. Notably, strain 3A-25B harbors a transposase gene associated with the Tn916 transposon, absent in the other three strains, belonging to a group of mobile genetic elements capable of genome integration, conjugative transfer to recipient cells, and intracellular transposition across genomic loci, as well as acquisition of accessory genes and promotion of genome rearrangements [34]. This distinctive genomic feature may have facilitated structural genome remodeling over time, thereby contributing to the phylogenetic divergence of strain 3A-25B from the other tree B. velezensis strains.

In contrast, the South American pan-genome (13 genomes) displayed a markedly different structure, with a reduced core genome (29.9%) and an expanded dispensable genome (53%), while singleton proportions remainded similar (~17%). This inverse relationship between core and accessory fractions suggests that geographic separation has contributed to increased evolutionary divergence and genomic plasticity. Phylogenetically, strain 3A-25B clustered closely with a Brazilian strain, whereas strains 2A-2B, 3A-6A, and 2A-10A remained grouped together (Figure 9). These results suggest greater intraspecific genomic diversity in South America, potentially associated with higher environmental heterogeneity.

At the global scale (30 genomes), the pan-genome structure resembled that observed for North America, with 38% core, 40% dispensable, and 22% singleton genes. The continuous expansion of the pan-genome and gradual decline of the core genome further support an open pan-genome model. Notably, South American strains exhibited a distinct genomic distribution, with a higher proportion of dispensable genes, aligning more closely with pan-genome structures reported in the ProPan database (Figure 12). The relatively high proportions of dispensable genes and singletons suggest considerable genomic plasticity, likely linked to environmental interactions, adaptation to local ecological niches, microbial competition in the rhizosphere, and specific selective pressures.

Globally, phylogentic relationships revealed that strains 3A-6A, 2A-10A y 2A-2B cluster with strains from diverse geographic origins, whereas strain 3A-25B form a distinct subclade with strains from Italy, China, and northern Mexico. This pattern may reflect the high dispersal capacity of B.velezensis, likely facilitated by endospore formation, enabling long-term survival under desiccation, radiation, and extreme temperatures [45], and long-distance transport via air currents [46], water systems [47], migratory organisms, or human-mediated activities [48]. Additionally, some Bacillus lineages exhibit relatively conserved pan-genomes with a stable core genome, supporting the occurrence of shared evolutionary origins among geographically distant strains [49]. In contrast, strain Sam8H1 forms a distinct phylogenetic branch; this strain originates from Botswana, Africa, where semi-arid to arid conditions (30–35 °C, high evaporation, nutrient-poor sandy soils) may have driven its marked evolutionary divergence and unique phylogenetic placement. Collectively, these findings suggest that the B. velezensis species possesses a conserved genetic backbone sufficient to maintain functional identity, complemented by a substantial proportion of variable genes that enhance its adaptive capacity and ecological success at a global scale. The significant contribution of accessory and unique (singleton) genes underscores the role of gene gain, including horizontal gene transfer, as a key evolutionary driver of diversity in this rhizobacterial species.

The analysis of plant growth-promoting traits (PGPTs) revealed a largely homogeneous functional profile among the four native strains, with genes associated with phytohormone biosynthesis (notably indole-3-acetic acid), motility, adhesion, detoxification, and plant cell wall degradation. Overall similarity with the Asian reference strain was high, although differences were observed in metabolic and plant cell wall-degrading enzyme categories. Notably, strain 3A-25B again showed greater similarity to the reference strain than to the other native strains.

Strain 3A-25B was distinguished by a higher abundance of genes related to biofertilization (iron acquisition, phosphate solubilization, sulfur assimilation) and bioremediation, suggesting enhanced ecological versatility. However, this was accompanied by a reduced number of genes associated with competitive exclusion, indicating a potential functional trade-off favoring host interaction over microbial competition.

A notable feature across all four strains was the presence of genes related to alternative CO₂ fixation pathways linked to Rubisco-like mechanisms, suggesting metabolic plasticity that may enhance survival under carbon-limited conditions [50]. Additionally, the high abundance of genes for nutrient acquisition (iron, nitrogen, phosphorus, potassium) and phytohormone signaling (IAA, cytokinins, ABA, gibberellins), along with volatile compound production (acetoin and butanediol), supports their classification as multifunctional plant growth-promoting rhizobacteria (PGPR).

Overall, the larger number of genes associated with plant colonization (~825) and microbial competitiveness (>600) indicates that plant-microbe interactions are governed by complex and redundant functional networks. Within this framework, strain 3A-25B stands out as functionally specialized, combining enhanced biofertilization and colonization capacity with reduced competitive exclusion, reinforcing its distinct adaptive trajectory among the four analyzed strains.

5. Conclusions

The four Bacillus velezensis native to north-central Mexico exhibit a functional and genomic profile that supports their strong potential as biocontrol agents and plant growth promoters, maintaining consistentt antifungal activity against chili phytopathogens without causing significant adverse effects on the host plant. The absence of statistically significant differences in antagonistic acctivity, despite underlying genomic variation, suggest that key antifungal mechanisms are highly conserved within the species, even among lineages with distinct evolutionary histories and broad geographic distribution.

At the genomic level, these strains share a robust and stable core genome, while variation in dispensable and singleton genes, particularly pronounced in strain 3A-25B, reflects ongoing microevolutionary diversification and differential gene acquisition. Although these differences do not translate into major changes in antifungal capacity, they are associated with ecological specialization, especially in traits related to plant colonization, biofertilization, and metabolic adaptation to environmental conditions.

Pan-genome analyses at regional and global scales indicate that B. velezensis possesses an open pan-genome, consistent with its wide geographic distribution and adaptability to diverse environments. Differences observed between North and South America, together with phylogenetic proximity among geographically distant strains, support the notion of a highly dispersible species, likely facilitated by endospore formation and multiple natural and anthropogenic vectors. Within this framework, strain 3A-25B emerges as a more genetically divergent lineage within the Mexican group, showing evolutionary affinities with strains from other regions, suggesting shared selective pressures or convergent evolutionary trajectories.

From a functional perspective, the abundance of genes associated with biofertilization, nutrient acquisition, phytohormone signaling, and volatile metabolite production underscores the multifunctional nature of these strains as plant growth-promoting rhizobacteria (PGPR). The consistently high proportion of unannotated genes across all scales highlights the presence of substantial “functional dark matter” with unexplored biotechnological potential. Particularly noteworthy is the identification of genes linked to alternative Rubisco-like carbon fixation pathways, suggesting a previously underappreciated metabolic plasticity in B. velezensis that may confer adaptive advantages under carbon-limited conditions.

Overall, the integratioin of phenotypic, genomic, and phylogenetic evidence positions these strains as valuable biological resources for sustainable agricultural applications, while also emphasizing the evolutionary complexity of B. velezensis and need for further functional studies to elucidate the ecologial and biotechnological roles of its variable and uncharacterized genes.