The bulbs of Lilium brownii (lily) were collected from Lintao County, Gansu Province. First, the bulbs were rinsed with sterile water for 1 min, repeating the process two to three times, and allowed to air dry. Second, the bulbs underwent disinfection with 75% ethanol and 3% sodium hypochlorite, followed by three rinses with sterile water to eliminate any residual ethanol and sodium hypochlorite. Third, the disinfected bulbs were crushed into a fine powder using a sterilized mortar. One gram sample was mixed with 9 mL sterile water and vortexed to ensure homogeneity. Fourth, the mixture was serially diluted to achieve concentrations of 10⁻⁶, 10⁻⁷, 10⁻⁸, and 10⁻⁹, and 100 μL aliquots of each dilution was spread onto solid Luria-Bertani (LB) medium. Finally, these LB plates were incubated at 30°C for 48 h. Distinct single colonies were selected for further isolation and purification, culminating in the successful isolation of the B. velezensis BN.

The culture conditions for the B. velezensis BN were as follows: LB medium with the NaCl concentration adjusted to 0.17 M, the pH maintained at 7.0, and the temperature at 30°C. The cells were observed using a fluorescence microscope (Yongxing N-800F, China).

To evaluate the antifungal activity of B. velezensis BN against diverse fungal pathogens, a dual-culture experiment was conducted, as previously described (Liu et al., 2024). The study encompassed six distinct fungal species: Colletotrichum gloeosporioides (C. gloeosporioides), Fusarium oxysporum (F. oxysporum), Rhizoctonia solani (R. solani), Phytophthora infestans (P. infestans), Fusarium graminearum (F. graminearum) and Fusarium tricinctum (F. tricinctum). Following a 7 days incubation, the pathogens were inoculated at the center of the Potato Dextrose Agar (PDA) plates. The PDA medium was prepared with 200 g L⁻¹ potato, 20 g L⁻¹ glucose, and 20 g L⁻¹ agar, with the pH adjusted to 7.0. B. velezensis BN was inoculated ~1 cm from the edge of each plate, with control plates containing only the pathogens. Each experiment was replicated three times to ensure statistical robustness. All plates were incubated at 28°C until the fungal growth in the control plates extended to the edges.

The genomic DNA of B. velezensis BN was extracted using a bacterial DNA extraction kit (General, Shanghai, China). The 16S rRNA gene was amplified using primers 27F (5′-AGAGTTTGATCMTGGCTCAG-3′) and 1492R (5′-GGYTACCTTGTTACGACTT-3′). Sequencing of the PCR products was performed by Sangon Biotech Co., Ltd., (China). The phylogenetic analysis of the 16S rRNA was conducted using the Type (Strain) Genome Server (TYGS) (https://tygs.dsmz.de/) (Le Han et al., 2022; Meier-Kolthoff and Göker, 2019).

Whole-genome sequencing of B. velezensis BN was performed at Majorbio Laboratory (Shanghai, China) and involved both Scaffolding and completing genome assemblies. For scaffolding, an Illumina paired-end (PE) library with an average insert size of ~400 bp was constructed and sequenced using the Illumina HiSeq platform. This process included de novo genome assembly to obtain preliminary genomic sequences and functional annotation. To refine the genome assemblies, a PacBio RS II library with a read length of ~10 kb was generated and sequenced on the PacBio RS II platform, followed by comprehensive bioinformatics analysis. Scaffolding was performed using SOAPdenovo2 software for the assembly of multiple short reads, while completion of genome assemblies were achieved with Unicycler v0.4.8 and Pilon v1.22 for sequence polishing.

Phylogenetic analysis of the gyrA gene and whole-genome sequences were conducted using the TYGS server (Le Han et al., 2022; Meier-Kolthoff and Göker, 2019). Average nucleotide identity (ANI) analysis was performed using the JSpeciesWS53 web server (https://jspecies.ribohost.com/jspeciesws/) (Richter et al., 2016). dDDH analysis was conducted using the Genome-to-Genome Distance Calculator 3.0 (GGDC) web server (https://ggdc.dsmz.de/ggdc.php#) (Meier-Kolthoff et al., 2013). Collinearity analysis of the four strains was carried out with Mauve software. Core and pan-genome analyses were performed using BPGA v1.3 software.

Functional annotation of the predicted coding genes was accomplished through a dual-pronged approach. Initially, we aligned the predicted genomic data with the genomes of the B. velezensis species, selected as reference genomes based on a phylogenetic analysis of 16S rRNA. This analysis facilitated the identification of homologous genes and their subsequent functional annotation. Subsequently, a database comparison method was employed, leveraging six major databases (NR, Swiss-Prot, Pfam, EggNOG, GO, and KEGG) to enhance the precision of functional descriptions for gene annotation.

Carbohydrate active enzymes were identified by the CAZy database (http://www.cazy.org/). Genes implicated in plant growth-promoting traits were analyzed by PGPT-Pred (https://plabase.cs.uni-tuebingen.de/pb/form.php?var=PGPT-Pred) (Patz et al., 2021). Furthermore, AntiSMASH 7.0 (https://antismash.secondarymetabolites.org/) (Blin et al., 2023) was employed to identify and characterize the biosynthetic gene clusters of secondary metabolites with molecular weights <2.5 kDa.

A bacterial strain isolated from the rhizosphere of Lilium brownii (lily) was designated as BN. Lily is widely recognized for its various health benefits, such as having calming properties, skin soothing effects, moisturizing dryness, and relieving cough. Additionally, research has shown that Lilium polysaccharides exhibit a variety of pharmacological activities, including antitumor, immunomodulatory, hypoglycemic, and radioprotective properties (Guo et al., 2024). The optimal culture conditions of BN were determined to be a temperature of 30°C, pH 7.0, and 0.17 M NaCl. BN is a Gram-positive bacterium, rod-shaped bacterium with length of 1–3 μm (Supplementary Figure S1A). It is characterized the formation of a substantial biofilm, which may potentially facilitate the colonization in the plant rhizosphere (Paula et al., 2020) (Supplementary Figure S1B).

Our study shown that BN exhibits a broad-spectrum antagonistic activity against several fungal plant pathogens, including C. gloeosporioides, F. oxysporum, R.solani, P. infestans, F. graminearum and F. tricinctum (Figure 1). C. gloeosporioides is a globally widespread plant pathogen that causes anthracnose in various plants, including wolfberry, legumes, and strawberries, leading to lesions on fruit and leaf surfaces (Dean et al., 2012). F. oxysporum is a destructive soil-borne pathogen that causes fusarium wilt and ceitocybe bescens in over 100 crops, including those in the Musaceae, Rosaceae, and Cucurbitaceae families (Berrocal-Lobo and Molina, 2008; Edel-Hermann and Lecomte, 2019). R. solani is a challenging pathogen to control which causes significant damage to a wide range of vegetables, fruits, and crops like potatoes, rice, corn, and cotton, leading to substantial yield and economic losses (Singh A. et al., 2023). P. infestans is responsible for causing the serious solanaceous plant disease known as late blight, particularly affecting potatoes, and was a major culprit in the 1840s European (Nowicki et al., 2012). F. graminearum is the primary causative agent of Fusarium head blight (FHB), one of the most destructive diseases affecting cereal crops worldwide (Dash et al., 2012; Lu et al., 2022). F. tricinctum is a ubiquitous and significant plant pathogen and mycotoxin producer that affects a variety of crops, notably causing Fusarium head blight (FHB) and ceitocybe bescens. In lilies, F. tricinctum can result in root rot, bulbs rot, and the yellowing and shedding of basal leaves (Li et al., 2013; Wang Y. et al., 2022). In summary, the BN strain shows promising potential as a highly effective broad-spectrum biocontrol agent.

We performed a phylogenetic analysis on the 16S rRNA gene sequence of the BN strain (Figure 2). The result revealed that BN exhibited the highest similarity with B. velezensis FZB 42, with a 16S rRNA sequence identity of 99.74%. B. velezensis FZB 42 serves as the model strain for promoting plant growth and biocontrol rhizosphere (Fan et al., 2018). Additionally, BN shared a sequence identity of 98.84% with B. velezensis B19. In summary, based on the sequence similarity of the 16S rRNA genes, coupled with the clustering patterns observed on the phylogenetic tree, BN was preliminarily identified as B. velezensis. The complete genome of BN was sequenced to further explore the antimicrobial capabilities and identify additional resistance genes, laying a foundation for developing engineered biocontrol strains.

The circular genome map of BN was constructed using the GC view program, as shown in Figure 3A (Grant and Stothard, 2008; Le Han et al., 2022). The whole genome sequence of BN consists of a circular chromosome spanning 3,929,791 bp, with a G+C content of 46.5%, and no plasmid sequences were detected. A total of 3,747 protein-coding genes (CDS) were predicted, accounting for 88.44% of the genome. These CDS genes have an average G+C content of 47.34% and an average length of 928 bp, inclusive of 281 pseudogenes. Additionally, the genome contains 81 small RNAs (sRNAs), 27 rRNAs (nine copies each of 5S, 16S, and 23S rRNAs), and 86 tRNAs.

The ANI and dDDH analyses are widely used to evaluate the similarity across bacterial strains based on whole-genome sequences. According to studies by Zhang et al. (2016), the compared strains exhibiting ANI values of ≥96% and DDH values of ≥70% are typically classified within the same species. In this study, the ANI values between BN and B. velezensis FZB 42, B. velezensis B19, and B. velezensis JS25R were 98.3%, 98.1%, and 98.2%, respectively (Figure 3B). Correspondingly, the DDH values were 85.4%, 83.4%, and 84.5% (Figure 3B).

The phylogeny based on the gyrA gene sequence from whole-genome sequences (Supplementary Figure S2), revealing that BN exhibited 98.76% and 98.33% sequence identity with the gyrA sequences of B. velezensis FZB 42 and B. velezensis B19, respectively. The gyrA gene, a highly conserved housekeeping gene in the Bacillus genus, exists as a single copy across all bacterial species, providing a higher resolution for Bacillus identification than the 16S rRNA gene, which may exhibit heterogeneity due to the multicopy (Liu et al., 2022). Furthermore, a phylogenomic tree was constructed based on the complete genome of BN and the published whole-genome sequences of other strains (Figure 3C). Collectively, these data fully confirmed the classification of the BN strain as a member of B. velezensis (Figures 2, 3, Supplementary Figure S2).

Comparative genomic analysis using whole genome sequencing with the Mauve program demonstrated that B. velezensis BN shares a high degree of homology with the majority of genes from B. velezensis FZB 42, B. velezensis B19, and B. velezensis JS25R. However, compared to B. velezensis FZB 42 and B. velezensis JS25R, certain genes were found to be absent in B. velezensis BN. Additionally, in comparison to B. velezensis B19, both deletion and inversion of specific genes were observed in B. velezensis BN (Supplementary Figure S3A). Core pan-genomic analysis revealed that the four strains collectively possess 3,266 core genes, constituting 87.2% of the total gene count of B. velezensis BN, further substantiating the high similarity between B. velezensis BN and the other three strains. Additionally, B. velezensis BN contains 111 unique genes (Supplementary Figure S3B).

COG (Clusters of Orthologous Groups) gene annotation reveals that the largest proportion of genes fall into three categories: Amino acid transport and metabolism, Transcription, and Carbohydrate transport and metabolism (Supplementary Figure S4A). Notably, genes associated with secondary metabolite biosynthesis, transport, and catabolism constitute 2.1% of the total coding genes.

Carbohydrate active enzymes (CAZy) are a crucial class of enzymes, classified into glycoside hydrolases (GHs), glycosyl transferases (GTs), polysaccharide lyases (PLs), carbohydrate esterases (CEs), Carbohydrate-Binding Modules (CBMs) and Auxiliary Activities (AAs). These enzymes are responsible for degrading, modifying, and generating glycoside bonds. The ability of CAZy to degrade carbohydrates not only provides a competitive advantage against other bacteria but also activates the immunity of host cell (Bu et al., 2023; Wardman et al., 2022). For instance, endophyte-derived α-mannosidase can activate the host immunity by degrading the rice cell wall and releasing oligosaccharides as DAMPs (damage-associated molecular patterns), thereby enhancing the rice's disease resistance (Bu et al., 2023). The function of the B. velezensis BN gene was annotated based on the CAZy database (Supplementary Figure S4B). The analysis identified 128 of CAZy, including 41 GHs, 42 GTs, two CBMs, 31 CEs, three PLs, and nine AAs.

PGPT-Pred is an analytical tool designed to predict plant growth-promoting traits within single bacterial genomes (Patz et al., 2021). As of the analysis conducted in August 2024, PGPT-Pred indicates that the gene distribution within B. velezensis BN is as follows: 26% are associated with plant colonization systems, 22% with competitive exclusion, 22% with stress control mechanisms, 12% with bio-fertilization, 10% with phytohormones, 7% with bioremediation, and 2% with plant immune response stimulation (Figure 4A). Notably, genes involved in plant colonization, competitive exclusion, and stress control, as well as plant immune response stimulation, indirectly influence plant growth, while those associated with bio-fertilization, phytohormone, and bio-remediation have a direct effect on promoting plant growth (Figure 4B, Supplementary Figure S5).

Biocontrol strains can produce a variety of secondary metabolites with diverse properties and structures, exhibiting a broad range of activities. These secondary metabolites include antibiotics, pigments, plant growth promoters, effectors of ecological competition, enzyme inhibitors, and pheromones (Chaabouni et al., 2012). Among them, compounds with a molecular weight below 2.5 kDa, produced during the stationary phase, represent the most important class of secondary metabolites for the control of plant diseases (Jayakumar et al., 2011; Ruparelia et al., 2022). They exhibit a wide range of antibacterial properties, triggering systemic acquired resistance, and promoting biofilm formation which facilitates the colonization for the producing strains. The antiSMASH platform is the leading resource for identifying and analyzing the biosynthetic gene clusters (BGCs) responsible for compounds with a molecular weight below 2.5 kDa in both archaea and bacteria (Blin et al., 2023).

The antiSMASH was employed to predict the BGCs in B. velezensis BN, aiming to assessing its capacity to suppress pathogenic activity. The analysis revealed a total of 12 BGCs dedicated to the production of compounds with a molecular weight below 2.5 kDa (Table 1), spanning a combined length of 754.8 kb, which constitutes a significant 19.2% of the genome's total length. Among the identified BGCs are two encoding nonribosomal peptide synthetases (NRPS), two trans-Acyl transferase polyketide synthases (transAT-PKS), one type III polyketide synthase (T3PKS), two transAT-PKS-NRPS hybrid systems, one lanthipeptides, one putative polyketide synthase (PKS), two terpene clusters, and one cluster of an undetermined type (Table 1).

Further analysis revealed that seven BGCs in B. velezensis BN exhibit a high degree of similarity to those found in the model strains B. velezensis FZB 42 and B. subtilis 168. The sequences of macrolactin H, bacillaene, fengycin, difficidin, bactillihactin, and bacilysin are 100% homologous with the reference strains, while surfactin show 82% similarity based on alignment with the public database (Table 1). Additionally, we compared the sequence differences of seven compounds among four strains: B. velezensis BN, B. velezensis FZB 42, B. velezensis JS25R, and B. velezensis B19. The results demonstrate that the four strains possess similar core biosynthetic genes, but B. velezensis B19 exhibits some differences compared to the other three (Figure 5), which aligns with the findings of the collinearity analysis (Supplementary Figure S3A). These differences may be attributed to the distinct environmental in which these strains evolved, leading to the development of unique metabolic pathways as adaptive mechanisms.

In this study, we conducted an antiSMASH analysis to provide a detailed summary of several key compound biosynthetic pathways (Figure 5). The biosynthetic gene cluster for surfactins is composed of four open reading frames (ORFs): licA, licB, licC, and srfATE, with a total length of 26.1 kb. In other species, the gene cluster ORFs related to surfactin biosynthesis are named srfAA, srfAB, srfAC, and srfAD (Ongena and Jacques, 2008). The synthesis of macrolactin H is catalyzed by a PKS encoded by seven genes, pksE, mlnB, mlnC, mlnD, mlnE, mlnF, and mlnG, totaling 48.2 kb. The biosynthesis of bacillaene involves a PKS composed of pksC, pksD, pksE, acpK, pksG, pksH, pksI, and pksL, along with a NRPS composed of pksJ, pksN, and pksR, resulting in a total gene cluster length of 70.1 kb. Fengycin is catalyzed by a NRPS composed of ppsA-E, with a total gene cluster length of 37.7 kb. Difficidin is catalyzed by a PKS composed of dfnD-J, totaling 60 kb. Bactillihactin comprises five ORFs, entA, entC, entE, entB, and dhbF, with a total length of 11.7 kb. In addition, bacillibactin is a siderophore that can absorb iron elements at extremely low iron concentrations. Given that iron atoms in nature primarily exist as ferric iron, which has low bioavailability at physiological pH, the role of bacillibactin is particularly crucial (Jia et al., 2021; Ma et al., 2017). The bacilysin is primarily catalyzed by a series of enzymes encoded by bacA-G, each performing distinct functions. The biosynthetic pathways of these natural active substances not only reveal their roles in microbial metabolism but also provide a potential for constructing artificial biocontrol engineering strains.