We followed the method described by Wang et al. (2021) to isolate the soil bacteria. Rhizosphere soil in which healthy tomato (Lycopersicon esculentum) plants grew was collected from areas infested with tomato root-knot nematode disease. Soil samples (5 g) were suspended in 50 ml of sterile distilled water and mixed on a table concentrator for 30 min. These soil samples were serially diluted (up to 10^–7-fold), plated on potato dextrose agar (PDA), and incubated at 30 ± 2°C for 2–3 days. Bacterial colonies growing on the plates were isolated according to their morphological characteristics for further study. Subsequently, the isolated strains were screened using in vitro tests, for which purpose, they were transferred to an LB solid medium plate. After activation at 30°C, purified single colonies were picked using a sterile inoculating loop and placed in a 250 ml conical flask containing 50 ml LB liquid medium. After 24 h of cultivation at 30°C, with stirring at 200 rpm, the inoculum was connected to the 50 ml LB culture medium at a ratio of 2% and cultured at 30°C and 200 rpm in a 250 ml conical flask for 48 h to prepare a fermentation broth. This broth was then centrifuged at 13,000 rpm for 10 min and filtered through a 0.22 μm sterile syringe filter. Lastly, 0.8 ml of the filtrates was placed in 1.5 ml centrifuge tubes before adding 100 s-stage juveniles (J2s) of M. incognita to each centrifuge tube. The corrected mortality of the J2s of M. incognita was calculated at 24 h and again at 48 h. Each treatment was repeated six times.

Strain JF-22 of B. licheniformis was preliminarily identified based on its morphological characteristics and then inoculated on PDA medium and observed for colony morphological characteristics after incubation at 32°C for 5 days. The strain was subsequently cultured for 24 h in LB liquid medium prior to Gram staining. Morphology was observed and photographed under a microscope. Further identification of JF-22 was achieved through the analysis of its 16S rDNA and gyrB gene sequences. The DNA of strain JF-22 was extracted with the Tiangen bacterial DNA extraction kit and stored at 4°C. 16S rDNA was amplified by the polymerase chain reaction (PCR) using the bacterial universal primers 27F (5′-AGAGTTTGATCCTGGCTCAG-3′) and 1492R (5′-GGTTACCTTGTTACGACTT-3′). The gyrB gene was amplified by PCR using specific primers G3F (5′-GAAGTCATCATGACCGTTCTGCAYGCNGGNGGNAARTTY GA-3′) and G3L(5′-AGCAGGGTACGGATGTGCGAGCCRTCN ACRTCNGCRTCNGTCAT-3′). The reaction conditions for 16S rDNA amplification were as follows: 95°C for 5 min; 33 cycles at 95°C for 1 min, at 58°C for 30 s, and at 72°C for 2 min, and a final extension at 72°C for 10 min. The reaction conditions of gyrB gene amplification were as follows: 96°C for 5 min, 32 cycles at 96°C for 30 s, at 59°C for 30 s, and at 72°C for 1 min, and a final extension at 72°C for 10 min. Then, the PCR products were purified and sequenced by a commercial sequencing company (Sangon Biotech, Shanghai, China). Nucleotide sequences were analyzed and compared with sequences in the GenBank database using BLAST. A phylogenetic tree of the JF-22 strain was constructed with 16S rDNA and gyrB gene sequences using the neighbor-joining method in Mega 7.0.21.

Meloidogyne incognita was isolated and cultured according to the method described by Jing et al. (2020). Tomato seeds of the tested variety were sown in a nursery tray, and when the seedlings reached the fourth true-leaf stage they were transplanted into a plastic pot 20 cm in diameter and infested with 300 J2s specimens of M. incognita per plant. Two days after transplanting, each tomato plant was watered with 100 ml of JF-22 bacterial suspension at a concentration of 10⁸ CFU/ml, and the same amount of sterile water was used as a control (CK). After inoculation, pots were randomly placed on the operating platform of a greenhouse. Each treatment was repeated 60 times. Tomato plants were uprooted separately to determine plant biomass indices. Disease severity was determined using the method of Lu et al. (2017), and the incidence index was determined. Rhizosphere soil around the tomato roots was collected and stored at −80°C for analysis of microbial community structure.

The volatile composition of strain JF-22 was determined using SPME–GC–MS as described by Wang et al. (2021). Strain JF-22 was inoculated on PDA slants placed in triangular flasks, which were then sealed with foil and incubated at 30°C. After 3 days, polydimethylsiloxane-divinylbenzene (65 μm) was used to collect volatile organic compounds (VOCs), and an RTX-5MS capillary column (60 × 0⋅25 μm ID × 0⋅25-μm thick film) was used to separate VOCs. Helium was used as the carrier gas at a rate of 1 ml min^–1 to analyze VOCs by GC–MS using the following procedure: desorb SPME-DVB fiber at 200°C for 3 min. For VOCs analysis, first, the oven temperature was raised to 40°C for 3 min, and then increased to 160°C at a rate of 8°C min^–1; it was maintained at that point for 2 min, and finally increased to 240°C at a rate of 15°C min^–1 and maintained at that point for 3 min. The mass spectrometer was set to 70 eV electron ionization mode, with a source temperature of 200°C. Continuous scanning from 45 to 500 m/z was used. VOCs were identified by comparison with data provided in the National Institute of Standards and Technology guidelines (number 14).

Total DNA extraction from the soil was performed using the BIO-TEK OMEGA Soil DNA Kit (Omega Bio-tek, Norcross, GA, United States). Concomitantly, the bacterial 16S rDNA V3-V4 region and the fungal rDNA-ITS gene were amplified. PCR amplification was performed as previously described (Sui et al., 2019). After mixing the PCR products of the same sample, they were separated using 2% agarose gels, and the recovered products were purified using AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, United States), and detected by 2% agarose gel electrophoresis. The recovered products were quantified using a Quantus fluorometer (Promega, United States). Furthermore, the NEXTFLEX Rapid DNA-Seq Kit was used for the DNA library construction. The library was constructed following these steps: (1) Adapter linking; (2) Magnetic bead screening to remove adapter self-ligating fragments; (3) PCR amplification to enrich library templates; (4) Magnetic beads recovery of PCR products to produce the final library. Sequencing was performed on the Illumina MiSeq PE300 platform. The raw sequences were quality controlled using the Trimmomatic software and spliced using the FLASH software as follows: (1) the 300 bp reads were truncated at any site receiving an average quality score of < 20 over a 50 bp sliding window, and the truncated reads shorter than 50 bp were discarded, reads containing ambiguous characters were also discarded; (2) according to the overlap relationship between PE reads, paired reads were merged into a sequence, and the minimum overlap length was 10 bp; (3) maximum mismatch ratio allowed in the overlap region of the spliced sequence was 0.2, and the non-conforming sequences were screened; and (4) samples were distinguished according to the barcodes and primers at the beginning and end of the sequence, and the sequence direction was adjusted. The number of mismatches allowed by the barcode was zero, and the maximum number of primer mismatches was 2. All reads were clustered with a 97% similarity cut-off using UPARSE (ver. 7.1)¹, and chimeric sequences were identified and removed using UCHIME (Zheng et al., 2015). The taxonomy of each 16S rDNA and ITS rDNA gene sequence was analyzed using the RDP Classifier against the Silva (SSU123) 16S rDNA database (Amato et al., 2013) and the UNITE 7.0/ITS database (Abarenkov et al., 2010), respectively, using a confidence threshold of 70%. Bacterial population functions were analyzed using the PICRUSt 2 database. Fungal ecosystem analyses were performed using the FUNGuild database (Nguyen et al., 2015). Raw reads were deposited in the NCBI Sequence Read Archive database (accession numbers: PRJNA803317 and PRJNA803323).

Data were analyzed using Microsoft Excel and SPSS 22.0 (SPSS Inc., Chicago, IL, United States). Duncan’s multi-range test and Student’s t-test were performed to compare the significance of the differences at p < 0.05.

For the selection of potential biocontrol agents against M. incognita, 15 bacterial strains were isolated, purified, and screened using in vitro assays to test the toxicity of the isolated fermentation supernatant on J2s of M. incognita. The results showed that the fermentation supernatant of four bacterial strains had toxic effects on M. incognita J2s. Among these, specifically, the fermentation supernatant of JF-22 had the strongest toxic effect on M. incognita J2s. When J2s of M. incognita were treated for 24 h with the fermentation supernatant of JF-22, the death rate of the J2s was 77%, and when they were treated for 48 h, mortality reached 89%. This indicated that JF-22 has a strong ability to kill M. incognita at the J2s developmental stage (Table 1).

The colony of strain JF-22 was round, white, and non-transparent, with a wrinkled surface, irregular edges, and a central bulge. Under the microscope, strain JF-22 was found to be a rod-shaped, gram-positive bacterium. Based on phylogenetic analysis of 16S rDNA and gyrB (Figure 1A), strain JF-22 was identified as B. licheniformis (Figure 1B).

The biological functions of B. licheniformis strain JF-22 were verified by a pot test. Compared with the control treatment CK, JF-22 significantly increased tomato biomass by 18.41%. The application of B. licheniformis strain JF-22 significantly reduced the severity of the root-nematode disease, indicating that B. licheniformis effectively controlled M. incognita (Table 2).

The types and relative contents of VOCs produced by the JF-22 fraction were analyzed using SPME–GC–MS. Figure 2 and Table 3 show SPME–GC–MS results for VOCs produced by JF-22 strain. Volatiles produced by JF-22 were identified using NIST17 and NIST17s (National Institute of Standards and Technology) standard mass spectrometry libraries. Among VOCs produced by JF-22, acetoin and 2,3-Butanediol, [R-(R*,R*)]- were the main components, with the relative peak areas of 64.2 and 7.89%, respectively, followed by cyclotrisiloxane, hexamethyl-; cyclotetrasiloxane, octamethyl-; 3-Hexanol, 2-methyl-; cyclopentasiloxane, decamethyl-; and 3-Pentanol and silanediol, dimethyl-.

After sequence optimization, there were 252,803 sequences in the two treatment groups, with an average of 42,134 in each soil sample of bacteria. And there were a total of 234,239 ITS sequences, with an average of 39,040 in each soil sample. The sequence generated by sequencing was clustered and segmented at a similarity level of 97%, and corresponding bacterial and fungal operational taxonomic units (OTUs) were generated, respectively. The bacterial community structure diversity indexes Shannon and Chao in the CK and JF-22 groups didn’t change significantly. However, in the richness index of fungal community, the fungal diversity indexes Shannon and Chao of the JF-22 group decreased significantly compared to CK (Table 4).

Using high-throughput sequencing, we analyzed the effects of JF-22 inoculants on the tomato–rhizosphere microbial communities. PCoA analysis based on unweighted UniFrac distance analysis showed that B. licheniformis strain JF-22 changed the composition of the microbial community in the tomato rhizosphere, particularly, that of the fungal communities (Figures 3A,B). Further, the bacterial community in the tomato rhizosphere soil across treatments was mainly composed of 25 families whose members were the same across treatments, although the relative abundance of each member differed. Thus, for example, abundance of Norank_c_Subgroup_6 and Nocardioidaceae differed greatly among treatments. In the CK treatment, the relative abundance of norank_c_Subgroup_6 was 3.6%, but in the JF-22-treated pots, the relative abundance of norank_c_Subgroup_6 was 6.6%, and the relative abundance of Nocardioidaceae was 6.3, and 3.9% in the CK and the JF-22 treatment groups, respectively (Figure 3C).

The fungal community in the tomato rhizosphere was mainly composed of 14 families across treatments. The overall composition of the rhizosphere fungal communities was similar in all treatments, but the abundance of each member was different. Specifically, Chaetomiaceae showed the largest difference between treatments. In CK pots, the relative abundance of Chaetomiaceae was 12.6%, while it was 39.7% in JF-22-treated pots. Similarly, Mortierellaceae and unclassified _o_Sordariales differed between treatments. In CK, the relative abundance of Mortierellaceae was 25.2%, while it was 13.6% in JF-22-treated pots. Furthermore, the relative abundance of unclassified _o_Sordariales was 4.3 and 18.2% in the CK and JF-22-treated pots, respectively (Figure 3D).

The number of unique and common OTUs in the multiple samples is shown in a Venn diagram. In the case of bacteria, 2,918 OTUs were identified in CK and 2,850 OTUs were found in JF-22, among which, 2,436 were common to both CK and JF-22-treated pots, while 482 were unique to CK and 414 to JF-22-treated pots (Figure 3E). As for fungi, 456 OTUs were identified in CK and 289 in JF-22-treated pots. Among them, 249 were common to both CK and JF-22-treated pots (Figure 3F). These findings indicated that the application of JF-22 can change the composition of the tomato rhizosphere–soil microbial community.

FUNGuild was used to analyze the functional guild annotation of the fungi present across treatments. In the CK group, the number of plant pathogens was 2,516, and the number of Dung Saprotrophs and Soil Saprotrophs was 1,221; meanwhile, in the JF-22 group, the number of plant pathogens was 739, and the number of Dung Saprotrophs and Soil Saprotrophs was 737, indicating that B. licheniformis strain JF-22 may reduce the number of plant pathogens, pathogenic bacteria, dung saprotrophs, and soil saprotrophs (Figure 4A). Compared to the results of FUNGuild, the results of the CoG function classification did not differ much between the experimental treatment groups (Figure 4B). All samples were rich in energy production and conversion (relative abundance 7.41–7.46%); amino acid transport and metabolism (8.85–8.95%); nucleotide transport and metabolism (2.36–2.40%); carbohydrate transport and metabolism (6.31–6.37%); coenzyme transport and metabolism (4.12–4.13%); lipid transport and metabolism (4.63–4.81%); translation, ribosomal structure, and biogenesis (5.13–5.14%); transcription (6.82–6.96%); replication recombination and repair (5.40–5.46%); cell wall/membrane/envelope biogenesis (5.98–6.14%); post-translational modification, protein turnover, chaperones (3.96–4.00%); inorganic ion transport and metabolism (5.69–5.69%); secondary metabolites biosynthesis, transport, and catabolism (2.69–2.79%); general function prediction only (8.71–8.74%); and signal transduction mechanisms (6.35–6.50%). These COG functional classification results indicate that the microbial communities of tomato rhizosphere–soil bacteria maintained similar biological functions in all treatment groups.