The microbial agent contained two main components, Paecilomyces lilacinus DZ910 and Bacillus subtilis KC1723. P. lilacinus DZ910 was isolated from Zibo, Shandong Province, China. B. subtilis KC1723 was isolated from Jinan, Shandong Province, China. P. lilacinus DZ910 and B. subtilis KC1723 strains were deposited in the China General Microbiological Culture Collection Center under preservation numbers CGMCC 40079 and CGMCC 24397, respectively.

Watermelon was planted in a solar greenhouse located in Renfeng Town, Jinan City, China (37.16° N, 117.36° E) on January 13, 2022. The test site had a planting density of 36,000 plants per hm² and has been used for growing watermelons for more than a decade. The soil is considered a fluvo-aquic type. The characteristics of the soil are as follows: pH, 7.76 ± 0.12; soil organic matter, 23.7 ± 0.6 g/kg. Ridge planting, conventional drip irrigation. All the treatments were managed in the same way.

The experiment consisted of four treatments: CK (untreated control), DZ910 (P. lilacinus DZ910), KC1723 (B. subtilis KC1723), D_K (mixture of P. lilacinus DZ910 and B. subtilis KC1723).Fungal strains were grown in PDA medium (P. lilacinus DZ910) and LB medium (B. subtilis KC1723). Conidial suspensions were prepared in sterile distilled water. Fungal conidial suspensions of P. lilacinus DZ910 (1 × 10⁹/mL), B. subtilis KC1723 (1 × 10¹¹/mL), and a combined suspension of DZ910 and KC1723 (1 × 10⁸/mL + 1 × 10¹⁰/mL, respectively) were prepared and used for experiments. The fungal conidial suspension used in this study was 90 l/hm², and diluted 10 times when used. Before transplanting, approximately 50 ml diluted suspension solution was poured into each hole. Watermelons of uniform size were used for this experiment. The daily field management of all treatments was consistent. Each treatment was replicated three times.

Three months after application of different microbial agents, chlorophyll content and photosynthesis parameters were measured using a portable instrument (Handy PEA, Hansatech). Three partitioned blades were tested, the lower blade (LB) was 50 cm away from the ground, the middle blade (MB) was 100 cm away, and the upper blade (UB) was 150 cm away. The experiment was repeated three times independently with 20 replicates each.

Virulence assays were performed at the same time. Ten leaves per plant were examined from 50 cm above the ground. According to the proportion of diseased area to total area, the disease index (DI) was classified into the following grades: 1, ≤5%; 3, 6–10%; 5, 11–25%; 7, 26–50%; 9, ≥50%. DI = Σ(number of disease plants×corresponding disease grade value)/(total number of plants investigated×highest value) × 100. The experiment was repeated three times independently with 15 replicates each.

Ten plants from each treatment were randomly assigned to assess watermelon fruit on April 19, 2022. The fruit were weighed and then used for quality determination. Fruit firmness, soluble solid content and sugar content were detected using the appropriate instruments following the corresponding operation instructions.

Total nitrogen (TN) contents were measured using an automatic Kjeldahl distillation titration unit (Foss, Sweden). Total phosphorus (TP) was digested with H₂SO₄-HClO₄ (Miranda et al., 2001). Total potassium (TK) was digested with HNO₃-HClO₄-HF (Shen et al., 2013). Soil total organic matter (OM) was assayed according to the vitriol acid potassium dichromate oxidation method (Nelson and Sommers, 2018).

Soil total DNA was extracted from 500 mg of fresh soil using a PowerSoil® DNA Isolation Kit (MO Bio Laboratories, San Diego, CA, USA) according to the manufacturer’s protocol. The bacterial universal primer pair 338F (5′-GTACTCCTACGGGA GGCAGCA-3′) and 806R (5′-GTGGACTACHVGGGTW TCTAAT-3′) targeting the V3-V4 region of the 16S rRNA gene was used to characterize soil bacterial communities (Kim et al., 2011). The fungal universal primer pair ITS1-F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and ITS2-R (5′-TGCGTTCTTCATCGATGC-3′) targeting the internal transcribed spacer (ITS) region was used to identify soil fungal communities (Martin, 2011).

The 16S rRNA gene and ITS regions in the DNA samples were amplified by polymerase chain reaction (PCR) before sequencing. The PCR was performed in a total volume of 50 μl, containing 2 μl of template DNA, 5 μl of 10 × buffer, 2 μl of 10 μM forward primer, 2 μl of 10 μM reverse primer, 4 μl of 2.5 μM dNTPs, 0.3 μl of DNA polymerase (2.5 U μL⁻¹), and 34.7 μl of ddH2O. Reaction conditions were 5 min at 95°C, followed by 30 cycles (for the bacterial 16S rRNA gene) or 35 cycles (for the fungal ITS region) of 30 s at 95°C, 30 s at 56°C and 40 s for elongation at 72°C, followed by a final extension at 72°C for 10 min and a hold at 10°C.

All PCR amplicons were examined using 2% (w/v) agarose gels. Purified DNA was quantified using a QuantiFluor™-ST (Promega, USA) and subsequently sequenced (PE300) using the Illumina MiSeq platform (Illumina Inc., USA) by Shanghai Majorbio Bio-pharm Technology Co., Ltd. Operational taxonomic units (OTUs) were assigned at a threshold of 97% similarity level and clustered using UPARSE (version 7.1) (Edgar, 2013). Chimeric sequences were identified and removed using UCHIME (Knight, 2011). Representative sequences of each OTU were taxonomically classified using the Silva and Unite database (Quast et al., 2012). Sequences matching the archaea, mitochondria and chloroplasts were removed and the number of sequences were rarefied to 32,650 reads per sample. The sequencing datasets supporting the current study are available in the National Center for Biotechnology Information (NCBI) under BioProject accession number PRJNA908896.

Data were analyzed by one-way analysis of variance (ANOVA) using SPSS 20.0 (SPSS Inc., USA). Alpha-diversity analyses, including ACE index, Chao index and Shannon index, were calculated using Mothur software (Quast et al., 2012). Taxonomic dissimilarity analysis between samples was based on the PCoA method with Bray-Curtis distances (beta diversity) (Lozupone et al., 2011). The R language vegan package (version 3.3.1) was used to draw community heat maps. The relationship between genus diversity and soil physicochemical properties was analyzed using redundancy analysis (RDA) by Canoco for Windows (ver. 4.5). The significance level for differences was set at p ≤ 0.05.

We measured chlorophyll content, photosynthetic efficiency and disease index of leaves to assess the growth of plants under different treatment (Figures 1A,B; Supplementary Figure S1a). Microbial agents had no obvious effect on chlorophyll content (Figures 1A). The photosynthetic efficiency of upper leaves (PE-up) and lower leaves (PE-low) were significantly higher under P. lilacinus (DZ910), B. subtilis (KC1723) and mixed microbial agent (D_K) treatments than in the control (CK) (p < 0.05) (Supplementary Figure S1a). The disease index of leaves was significantly lower under DZ910 and D_K treatments than under CK and KC1723 treatments (p < 0.05) (Figure 1B).

The single fruit weight (SWF), soluble solids and sugar degree were determined to assess the growth of watermelon fruit (Figures 1C–E; Supplementary Figure S1b). The D_K treatment produced the highest SWF, which was significantly higher than that in other treatments (Figure 1C; Supplementary Figure S1b). Sugar degree was significantly higher in the D_K treatment compared with the other treatments (p < 0.05) (Figure 1E). However, there were no significant differences in the soluble solids of fruit among different treatments (Figure 1D).

The effects of using P. lilacinus (DZ910), B. subtilis (KC1723) and the mixed microbial agent (D_K) on soil physicochemical characteristics are shown in Figure 1. Contents of TN and TK did not differ significantly between treatments (Figures 1F,G). TP in the D_K treatment was significantly higher than that in the other treatments (Figure 1H). Likewise, soil OM in the DZ910 treatment was significantly higher than that in other treatments (Figure 1I).

We evaluated and compared the richness and diversity of bacterial and fungal communities in the different soil samples using Chao1, Shannon and ACE indices (Figure 2). As shown in Figures 2A,D, the Chao1 estimator for the bacterial community showed no significant difference between samples, whereas a lower Chao1 estimator for the fungal community was observed in the soil treated with KC1723. The Shannon index of the bacterial community was higher in the sample under DZ910, KC1723 and D_K treatments than under CK, with the highest value in the soil treated with DZ910 (Figure 2B). However, the Shannon index for the fungal community showed a lower value in the soil under KC1723 and D_K treatment (Figure 2E). The ACE index value for the bacterial community was significantly higher under DZ910 treatment than under other treatments, while the ACE index value for the fungal community was significantly lower under KC1723 treatment than under other treatments (Figures 2C,F).

The results of principal coordinate analysis (PCoA) of different soil samples are also shown in Figure 2. Variation in bacterial and fungal community composition was affected by 11 principal coordinate components, respectively (Figures 2G,H). In terms of bacteria, the cumulative contribution rates of the two maximally reflected differences were 42.09 and 14.51% and the cumulative explanatory variation was 56.6% (Figure 2G). The distance between CK and D_K treatments was small, and they grouped together. DZ910 and KC1723 treatments were both far away from CK and D_K treatments. In fungi, the two main coordinates extracted explained 63.25% of the variation, of which 39.19% was explained by PC1 and 24.06% by PC2 (Figure 2H). Similarly, the results showed that CK, DZ910 and D_K treatments gathered together. However, the KC1723 treatment was far away from CK, DZ910 and D_K treatments.

The bacterial and fungal populations in the different treatments were analyzed by high-throughput sequencing. In this study, a total of 26 bacterial genera and 20 fungal genera (relative abundance >1%) were detected in the different soil samples (Supplementary Tables S1, S2). To further analyze the relative abundance of the microbial community, we found that 3 identified bacterial genera, including Dongia, Pseudomonas and Flavobacterium, and five identified fungal genera, Aspergillus, Myceliophthora, Trichoderma, Humicola and Neocosmospora, showed significant differences in the soil under different treatments (Supplementary Tables S1, S2; Figure 3). In the bacterial community, OTUs belonging to Dongia were decreased in DZ910 and KC1723 treatments, but showed no significant change in D_K compared with those in CK (Figure 3). The abundance of Pseudomonas increased significantly in KC1723 and D_K, and the abundance of Flavobacterium increased significantly in all microbial agent treatments (DZ910, KC1723 and D_K) compared with those in CK (Figure 3). Within the fungal community, the abundance of Aspergillus decreased significantly in KC1723 and the abundance of Myceliophthora and Neocosmospora decreased significantly in KC1723 and D_K compared to those in CK (Figure 3). However, the abundance of Trichoderma and Humicola increased significantly in the KC1723 treatment compared to CK (Figure 3).

We constructed a correlation heatmap between microbial community at the genus level and physicochemical properties using Spearman’s correlation analysis (Figure 4). Among the top 50 most abundant bacterial genera, 16 genera were affected by at least one physicochemical property (Figure 4A). Nitrospira, Chryseolinea and norank_f_67–14 were significantly negatively correlated with TN. Norank_f_norank_o_norank_c_KD4-96 was significantly negatively correlated with TP. Nine genera were affected by TK, one of which was negatively affected by TK, while the other eight genera, Dongia, unclassified_k_norank_d_Bacteria, norank_f_JG30-KF-CM45, Nonomuraea, norank_f_AKYG1722, Streptomyces, Steroidobacter and norank_f_norank_o_Rokubacteriales, were positively correlated with TK. Four genera, norank_f_TRA3-20, norank_f_norank_o_Actinomarinales, MND1 and norank_f_67–14, were significantly negatively correlated with OM.

In the top 50 fungal genera, 17 genera were affected by at least one physicochemical property (Figure 4B). Kernia, Conocybe, Chordomyces and Cladosporium were significantly positively correlated with TN, and Coprinus was significantly negatively correlated with TN. Cheilymenia, Gibellulopsis, Madurella and unclassified_c_Eurotiomycetes were significantly positively correlated with TP. Three genera, unclassified_c_Eurotiomycetes, Stachybotrys and Amesia, were significantly positively affected by TK, and two genera, Mortierella and Trichoderma were significantly negatively affected by TK. Six genera were affected by OM: Kernia, Cladosporium, unclassified_o_Sordariales, Apiotrichum, Cyberlindnera and unclassified_o_Hypocreales.

We conducted network analysis among physicochemical properties, major microorganisms and agronomic variables to assess their relationships. There were 22 nodes and 44 edges in the network (Figure 5). TK was the physicochemical characteristic most related to microorganisms. It was positively related to Dongia, Amesia and Retroconis, and negatively related to Trichoderma and Humicola. TN was strongly associated with Nitrospira, which was similar to the results of the correlation heatmap (Figure 5). Interestingly, two bacterial genera, Pseudomonas and Flavobacterium, were positively related to SWF, PE-up and PE-low. At the same time, a positive correlation was observed between them. The fungal genus Myceliophthora was positively related to Neocosmospora and Aspergillus, and they were all negatively related to SWF and PE-up.