This study was conducted in native Mexican maize fields in the rural community of San Juan de las Manzanas, Ixtlahuaca de Rayón, Estado de México (−99.842025 and 19.556658). The experimental plots were established using a randomized complete block design (RCBD) with five replicates per treatment (Piepho et al., 2004). Each plot received one of the following treatments: no application (control) (NA), conventional fertilization backpack (MF) and drone-assisted (DF), and biofertilization with SynCom applied manually (MB) and via drone technology (DB) (Mogili and Deepak, 2018; Saleem et al., 2019). Comprehensive site characterization was performed for each field to account for environmental variability. This included soil type classification (IUSS Working Group WRB, 2015), collection of historical crop data, and recording climatic conditions throughout the growing season (Lobell et al., 2011). These environmental parameters were integrated into subsequent data analyses to assess treatment efficacy across diverse agroecological contexts.

Based on functional traits, SynCom enhanced plant growth, nutrient cycling, and pathogen suppression. The formulation process involved screening bacterial isolates for nitrogen fixation, phosphate solubilization, and production of phytohormones such as indole-3-acetic acid (IAA). The SynCom comprised seven beneficial bacterial strains isolated from teosinte seeds: Serratia nematodiphila EDR2, Klebsiella variicola EChLG19, Bacillus thuringiensis EML22, Pantoea agglomerans EMH25, Bacillus thuringiensis EBG39, Serratia marcescens EPLG52, and Bacillus tropicus EPP72 (Table 1) (De-la-Vega-Camarillo et al., 2023a). Each bacterial strain was cultured independently in R2A broth media at 28°C with orbital shaking at 150 rpm until reaching an optical density (OD600) of 1.0 (Reasoner and Geldreich, 1985). Cultures were then centrifuged at 13,000 rpm for 10 min, and the resulting pellets were resuspended in sterile saline solution (0.85% NaCl) to achieve a final concentration of 1 × 10⁹ CFU/ML (De-la-Vega-Camarillo et al., 2023b). The SynCom was prepared for field application by mixing equal volumes of each bacterial suspension. The viability of each strain in the final mixture was confirmed through colony-forming unit (CFU) counts on selective media (Timmusk et al., 2017).

The treatments were applied using drone-assisted and backpack methods. For drone application, a DJI Agras MG-1S drone equipped with four precision nozzles was utilized (Zhang and Kovacs, 2012). The drone was calibrated to deliver 50 mL of treatment solution per plant, with flight paths programmed to ensure uniform coverage of each plot (Tripicchio et al., 2015). Backpack applications were conducted using backpack sprayers, calibrated to deliver 50 mL per plant, following standard agronomic practices (Szilagyi-Zecchin et al., 2016). Applications were carried out at four key growth stages: pre-treatment at 0 days after sowing (DAS), mid-season at 25 DAS, early reproductive stage at 55 DAS, and harvest at 85 DAS, aligning with critical periods in maize development (Abendroth et al., 2011). Weather conditions, including temperature, humidity, and wind speed, were recorded during each application to account for potential variability in treatment efficacy (Sánchez et al., 2014). These environmental parameters were later incorporated into the data analysis to assess their impact on treatment performance.

Soil samples were collected from each experimental plot 5 days after the application to observe the changes in the microbial communities due to fertilization and biofertilization processes as follows: 0 DAS (pre-treatment), 30 DAS (mid-season), 60 DAS (early reproductive stage), and 90 DAS (post-harvest) (Ma and Biswas, 2015). Each sample comprised five subsamples taken from the corners and center of each plot at a depth of 0–15 cm using a sterilized soil drill (Peigné et al., 2013). The subsamples were homogenized to create a composite sample (~500 g) per plot, following standard soil sampling protocols (Larkin, 2015). Immediately after collection, soil samples were placed in sterile, airtight containers and transported on ice to the laboratory. To preserve microbial community structure and DNA integrity, samples were stored at −80°C until further processing (Rissanen et al., 2010).

Soil physicochemical properties were analyzed following standardized protocols. Before analysis, samples were air-dried at room temperature (25 ± 2°C), ground, and sieved through a 2-mm mesh to remove coarse particles and plant debris (International Organization for Standardization, 2006). The following parameters were assessed: pH (1:2.5 soil/water ratio), electrical conductivity (EC), organic matter content (Walkley-Black method), total nitrogen (Kjeldahl method), available phosphorus (Olsen method), and exchangeable potassium (ammonium acetate extraction) (Bao, 2000; Jones, 2001).

Soil nutrient availability was assessed using a combination of satellite and drone-based multispectral imaging. Satellite data were obtained from Sentinel-2 MSI (10 m resolution) (Drusch et al., 2012). At the same time, high-resolution imagery (3 cm/pixel) was acquired using a DJI Phantom 4 Multispectral drone equipped with a six-band multispectral camera (blue, green, red, red edge, near-infrared, and RGB) (Hassan et al., 2019). Drone flights were conducted at 40 m altitude between 10:00 and 14:00 under clear sky conditions (cloud cover < 10%), coinciding with soil sampling dates (Assmann et al., 2019). Radiometric calibration was performed using a calibrated reflectance panel before each flight (Wang and Myint, 2015). Both satellite and drone imagery were atmospherically corrected using Sen2Cor (version 2.8) and Pix4D Mapper (version 4.6.4), respectively (Main-Knorn et al., 2017; Pix4D SA, 2020).

Spectral indices were calculated to estimate specific soil nutrients: Modified Soil-Adjusted Vegetation Index (MSAVI) for nitrogen content, Band Ratio Phosphorus Index (BRPI) for available phosphorus, and Potassium Abundance Index (KAI) for exchangeable potassium (Wang et al., 2014; Song et al., 2023). Additional indices included the Iron Oxide Ratio (IOR) and Normalized Difference Salinity Index (NDSI) for micronutrient availability assessment (Ge et al., 2011). Calibration models were developed using partial least squares regression (PLSR) to relate spectral indices to laboratory-measured nutrient concentrations (Viscarra Rossel and Behrens, 2010). Model validation was performed using a subset of soil samples (30%) that was not used in calibration.

For each experimental plot, spectral data were extracted from both satellite (10 m² pixels) and drone imagery (9 cm² pixels) using zonal statistics in QGIS 3.22 (QGIS Development Team, 2024). Temporal variations in nutrient availability were analyzed using repeated measures ANOVA with Bonferroni correction for multiple comparisons (p < 0.05). The relationship between spectral-derived and laboratory-measured nutrient values was assessed using Spearman correlation coefficients and root mean square error (RMSE) calculations (Gorelick et al., 2017).

Total genomic DNA was extracted from 0.5 g of soil using the DNeasy PowerSoil Kit (Qiagen, Hilden, Germany), following the manufacturer's instructions with an additional bead-beating step for enhanced cell lysis (Lim et al., 2010). The V3–V4 regions of the 16S rRNA gene were amplified using 341F/785R (CCTACGGGNGGCWGCAG/GACTACHVGGGTATCTAATCC) primers for bacteria (Klindworth et al., 2013), and the internal transcribed spacer (ITS) regions were amplified using ITS1F/ITS2 (CTTGGTCATTTAGAGGAAGTAA/GCTGCGTTCTTCATCGATGC) primers for fungi (Gardes and Bruns, 1993); both with unique 6-nucleotide barcodes for sample identification. PCR reactions were performed in 30 μL volumes containing 15 μL Phusion Master Mix (New England Biolabs, Ipswich, MA, USA), 0.2 μm of each primer, and 10 ng of template DNA. The thermal cycling conditions were as follows: initial denaturation at 98°C for 1 min, followed by 30 cycles of 95°C for 10 s, 50°C for 30 s (16S primers)/56°C for 30 s (ITS primers), and 72°C for 30 s, with a final extension at 72°C for 5 min (Wu et al., 2015). Amplified products were purified using the GeneJET Gel Extraction Kit (Thermo Fisher Scientific, Waltham, MA, USA) and verified by 2% agarose gel electrophoresis. Sequencing was performed on the Illumina MiSeq platform at Genome Science Core, Wayne State University (Detroit, USA), with a read length of 250 bp paired-end (Caporaso et al., 2011).

Raw sequencing data were processed using QIIME2 (version 2021.2) for demultiplexing and quality filtering (Bolyen et al., 2019). Reads with a Phred score below 20, chimeric sequences, and those shorter than 200 bp were discarded. Operational taxonomic units (OTUs) were delineated at 97% sequence similarity using DADA2 (Callahan et al., 2016). Taxonomic identities were attributed in comparison to the SILVA database (v13.8) for bacterial sequences (Quast et al., 2012) and the UNITE database (v132) for fungal sequences (Nilsson et al., 2019).

Microbial community diversity was assessed using both alpha and beta diversity metrics. Alpha diversity was quantified using Shannon, Simpson, Chao1, ACE (Abundance-based Coverage Estimator) indices, and observed species richness, calculated using the scikit-bio library (scikit-bio Development Team, 2020) in Python 3.8 (Van Rossum and Drake, 2009). Rarefaction curves were generated to evaluate sampling depth adequacy using the rarefaction_curve function from scikit-bio. Beta diversity was assessed using weighted and unweighted UniFrac distances (Lozupone et al., 2011), Bray-Curtis dissimilarity index (Bray and Curtis, 1957), Jaccard index (Jaccard, 1912), and Hellinger distance (Legendre and Gallagher, 2001), computed with the scipy.spatial.distance module (Virtanen et al., 2020). To visualize relationships between microbial communities, we performed Principal Coordinate Analysis (PCoA) and Non-metric Multidimensional Scaling (NMDS) using the scikit-learn library (Pedregosa et al., 2011). Differences in alpha diversity metrics between treatment groups were assessed using one-way ANOVA followed by Tukey's HSD post-hoc test, implemented with the scipy—stats module. For beta diversity, the statistical significance of community differences was tested using Permutational Multivariate Analysis of Variance (PERMANOVA) with 999 permutations, implemented through the skbio.stats.distance module. To elucidate potential microbial interactions, co-occurrence networks were constructed using three complementary approaches: the SparCC algorithm (Friedman and Alm, 2012) implemented in the FastSpar package (Watts et al., 2019), the CoNet method (Faust and Raes, 2016) using the RMT-based approach in the MENA package (Deng et al., 2012), and SPIEC-EASI (Kurtz et al., 2015) using the SpiecEasi Python package. Networks were visualized using NetworkX (Hagberg et al., 2008) and plotted with Matplotlib (Hunter, 2007). Network properties, including modularity, average path length, and clustering coefficient, were calculated using NetworkX. To identify microbial taxa significantly enriched or depleted between different treatments, we performed differential abundance analysis using DESeq2 (Love et al., 2014) through the DESeq2 Python package, considering taxa with an adjusted p-value < 0.05 and absolute log₂ fold change > 1 as significantly differentially abundant. All statistical analyses were performed in Python 3.8, and plots were generated using Matplotlib and Seaborn (Waskom, 2021).

Plant height was measured from the base to the flag leaf collar at the silking stage using a measuring tape. A plant canopy analyzer assesses leaf area index (LAI) (Welles and Norman, 1991). Chlorophyll content was determined using a SPAD-502 meter (Uddling et al., 2007). Nitrogen, phosphorus, and potassium content in leaves were determined from samples collected at the tasseling stage. Dried and ground leaf samples were analyzed using the Kjeldahl method for nitrogen and inductively coupled plasma optical emission spectrometry (ICP-OES) for phosphorus and potassium (Jones et al., 1991). Water use efficiency (WUE) was calculated as the ratio of grain yield to total water use during the growing season (Passioura, 2006). Nitrogen use efficiency (NUE) was determined as the ratio of grain yield to total nitrogen applied (Moll et al., 1982). Maize was harvested at physiological maturity, ~120 days after sowing (DAS) (Abendroth et al., 2011). Yield was measured as tons per hectare (tons ha⁻¹) by weighing the total biomass of ears per plot after shelling and adjusting to 14% standard moisture content (Zia et al., 2013). The harvest index was calculated as the ratio of grain yield to total aboveground biomass (Hay, 1995; Badu-Apraku et al., 2012). Kernel weight was assessed by weighing 1,000 grains per plot (CIMMYT standard procedures). Cob length and diameter were measured using a digital caliper (Carcova and Otegui, 2001).

Statistical analysis of yield data was performed using one-way ANOVA, followed by Tukey's HSD post-hoc test for pairwise comparisons, with significance set at p < 0.05. Statistical analyses were performed using R software (version 4.1.0, R Foundation for Statistical Computing, Vienna, Austria) (R Core Team, 2021).

All statistical analyses were performed using Python (version 3.9.19) (Van Rossum and Drake, 2009). The following libraries were utilized: pandas (version 2.2.2) for data manipulation and preprocessing (McKinney, 2010), scikit-learn (version 1.5.1) for machine learning-based feature selection (Pedregosa et al., 2011), statsmodels (version 0.14.2) for inferential statistics and hypothesis testing (Seabold and Perktold, 2010), and SciPy (version 1.12.0) for additional statistical computations (Virtanen et al., 2020). Data visualization was performed using matplotlib (version 3.9.2) and seaborn (version 0.13.2) (Hunter, 2007; Waskom et al., 2020).

To assess the relationship between microbial diversity and crop performance, multiple statistical approaches were implemented. The normality of the data was tested using the Shapiro–Wilk test, while Levene's test assessed the homogeneity of variances. When assumptions of parametric tests were met, one-way ANOVA followed by Tukey's HSD post-hoc test was applied to compare means across different treatments. For non-parametric comparisons, the Kruskal–Wallis test, followed by Dunn's test with Bonferroni correction, was performed.

Multiple linear regression models were constructed using the ordinary least squares (OLS) method in the statsmodels library to evaluate the influence of microbial diversity and environmental factors (temperature, precipitation, and soil pH) on crop yield and health metrics. The regression model used was:

Where Y represents the dependent variable (e.g., crop yield or plant health index), X₁, X₂,..., X_n are independent variables (diversity indices, environmental covariates, and fertilization treatment), β₀, β₁, β₂,..., β_n are regression coefficients, and ε is the error term. Model assumptions were verified using quantile-quantile (Q-Q) plots, the Shapiro–Wilk test for normality of residuals, and the Breusch-Pagan test for homoscedasticity.

Significance levels were set at p < 0.05, and 95% confidence intervals were calculated for all parameter estimates. Effect sizes were reported where applicable using Cohen's f² for regression models and eta-squared (η²) for ANOVA. All statistical analyses were conducted in a fully scripted and reproducible manner, with code available upon request.

This analysis showed a remarkably high precision in the correlation between laboratory-measured soil physicochemical properties and spectral indices derived from georeferenced satellite imagery. A highly positive correlation (r ≈ 0.98, p < 0.001) was evident between laboratory analyses and spectral image analyses for corresponding sample points (e.g., S1 laboratory vs. S1 spectral, S2 laboratory vs. S2 spectral, and so on up to S50) (Figure 1).

Importantly, this high correlation was observed only between matching sample points and not between different field areas (see red diagonal series in Figure 1). For instance, the spectral data for sample S1 showed a near-perfect correlation with the laboratory data for S1 but not with laboratory data from S2, S3, or any other sample points. This pattern was consistent across all 50 sample points (S1–S50).

The strength and specificity of these correlations underscore the accuracy of spectral image analysis in capturing soil physicochemical properties. Critical parameters such as pH, organic matter (OM), total nitrogen (TN), and available phosphorus (AP) all showed correlations of r > 0.95 (p < 0.001) between laboratory and spectral measurements for each respective sample point.

The application of different fertilization strategies significantly influenced soil properties, nutrient uptake, plant physiological traits, maize yield, and disease incidence. The heatmap analysis (Figure 2a) showed strong correlations between fertilization treatments and key soil and plant parameters. Biofertilization treatments, particularly drone biofertilization, were associated with higher soil microbial biomass (0.72 ± 0.05), improved soil aggregate stability (0.78 ± 0.03), and increased soil organic matter content (0.86 ± 0.04). In contrast, conventional fertilization resulted in lower microbial biomass (0.46 ± 0.04) and aggregate stability (0.48 ± 0.02), suggesting a reduced impact on soil structure and long-term fertility. Soil pH was significantly higher in biofertilization treatments (6.8 ± 0.1) compared to conventional fertilization (6.3 ± 0.2) and untreated controls (5.9 ± 0.3) [ANOVA, F(4, 95) = 21.7, p < 0.001].

Biofertilization also enhanced plant nutrient uptake and physiological performance. Nitrogen content in leaves was highest under drone biofertilization (0.86 ± 0.02%), followed by backpack biofertilization (0.80 ± 0.03%) and conventional drone fertilization (0.74 ± 0.02%). In contrast, conventional fertilization and untreated controls had significantly lower nitrogen content (0.63 ± 0.04% and 0.41 ± 0.03%, respectively) (p < 0.05). Similar phosphorus and potassium uptake trends were observed, with the highest concentrations detected in biofertilization treatments. Chlorophyll content was also significantly higher in drone biofertilization (48.2 ± 1.8 SPAD) compared to conventional fertilization (42.7 ± 2.1 SPAD) and untreated plants (35.4 ± 2.5 SPAD) (p < 0.001), indicating improved photosynthetic capacity.

Plant morphological responses aligned with these physiological changes. Plants in drone biofertilization treatments exhibited greater height (2.1 ± 0.1 m) and larger leaf area index (3.9 ± 0.2) compared to conventionally fertilized (1.8 ± 0.1 m, 3.2 ± 0.1) and untreated plants (1.4 ± 0.2 m, 2.5 ± 0.2) (p < 0.05). These growth improvements translated into enhanced cob development, with drone biofertilization resulting in the most considerable cob length (19.4 ± 0.7 cm) and diameter (5.6 ± 0.2 cm), whereas untreated plants produced significantly smaller cobs (11.3 ± 0.9 cm length, 4.1 ± 0.3 cm diameter) (Figure 2c).

Yield outcomes reflected these physiological advantages. Drone biofertilization resulted in the highest mean yield (7.4 ± 0.03 tons/ha), significantly outperforming backpack biofertilization (7.2 ± 0.14 tons/ha), conventional drone fertilization (6.7 ± 0.19 tons/ha), conventional fertilization (6.3 ± 0.33 tons/ha), and untreated controls (2.5 ± 0.21 tons/ha) [ANOVA, F(4, 95) = 42.6, p < 0.001]. Post-hoc Tukey's HSD tests confirmed that drone biofertilization yielded significantly higher than all other treatments (p < 0.05) (Figure 2b).

Disease incidence varied significantly across treatments, with biofertilization demonstrating a strong protective effect. Kruskal–Wallis tests, followed by Dunn's post-hoc comparisons, showed that drone biofertilization significantly reduced the incidence of Fusarium wilt (20% [IQR: 15–25%]) compared to conventional fertilization (55% [IQR: 50–65%]) and untreated controls (80% [IQR: 70–85%]) (p < 0.05). Similar reductions were observed for root rot, stalk rot, and head smut, indicating that biofertilization, particularly when delivered via precision technology, enhances plant resilience against pathogenic infections (Figure 2b).

The visual assessment of cob quality and field performance further supported these quantitative findings. Cobs from drone biofertilization treatments exhibited more uniform kernel filling, greater grain weight, and lower levels of kernel abortion compared to conventionally fertilized and untreated plants (Figure 2c). Additionally, maize plants under biofertilization showed increased vigor, with more robust stems and healthier foliage at the flowering stage compared to the control (Figure 2d).

The metabarcoding sequencing of soil samples yielded 2.5 million raw reads across 40 libraries, representing five fertilization treatments: no application, backpack fertilization, drone fertilization, backpack biofertilization, and drone biofertilization, with four sampling points and two replicates each. After quality filtering and chimera removal, 1.8 million high-quality sequences were retained, clustering into 57 distinct genera (27 bacterial and 30 fungal OTUs). Taxonomic classification revealed four major phyla across 14 classes, 25 orders, and 34 families, with Ascomycota (≈65%) and Basidiomycota (≈25%) dominating the fungal community. Z-score analysis revealed distinct microbial distribution patterns across treatments, significantly enriching beneficial microorganisms in biofertilization treatments. Control soils (NA) showed a higher abundance of Cladosporium (Z = 1.55), Setophoma (Z = 1.33), and Enterobacter (Z = 1.32), while biofertilized soils, particularly under drone application (DB), exhibited significant enrichment of plant growth-promoting microorganisms including Rhizobium (Z = 1.51), Mortierella (Z = 1.53), and Penicillium (Z = 1.53) (Figure 3A).

Bacterial communities exhibited complementary patterns, with plant growth-promoting rhizobacteria showing significant enrichment in biofertilization treatments. Rhizobium abundance increased markedly in DB (Z = 1.51) and MB (Z = 1.00) treatments, mainly from time points 2 to 4. Similarly, Pseudomonas and Bradyrhizobium showed progressive enrichment in biofertilization treatments (Z-scores increasing from 0.54 to 0.98 and 0.32 to 1.01, respectively). The drone application methods (DB and DF) showed more stable community compositions across replicates than backpack applications (average standard deviation: drone = 0.24, backpack = 0.41). Temporal analysis revealed distinct succession patterns: initial time points (1–2) showed moderate shifts from control conditions (average ΔZ = 0.45), while later time points (3–4) exhibited more pronounced community restructuring (average ΔZ = 0.89). This temporal progression was particularly evident in biofertilization treatments, where beneficial microorganisms showed consistent enrichment patterns (temporal correlation coefficient r = 0.78, p < 0.001). Post-hoc Dunn's tests with Benjamini-Hochberg correction confirmed that biofertilization treatments, particularly drone-assisted application, resulted in significantly higher microbial richness than other treatments (p < 0.05). The most pronounced differences were observed between DB4 and NA1 treatments (average ΔZ = 2.14 for fungi, ΔZ = 1.89 for bacteria). The progressive increase in Z-scores from NA to DB treatments, particularly evident in time points 3 and 4, suggests a cumulative positive effect of biofertilization on microbial community structure (Figure 3B).

Bacterial and fungal communities analysis revealed significant differences in taxonomic composition and abundance patterns across treatments and time points (Figure 4). Kruskal–Wallis tests showed substantial variations among treatments for both fungi [χ²(4) = 38.2, p < 0.001] and bacteria [χ²(4) = 42.7, p < 0.001]. For fungi, drone biofertilization (DB) showed the highest enrichment in beneficial genera, with Mortierella (Z = 1.53), Penicillium (Z = 1.53), and Fusarium (Z = 1.56) increasing significantly from time points 1 to 4. Backpack biofertilization (MB) showed similar but less pronounced trends (average Z-scores rising from 0.68 to 1.21 across time points). In contrast, chemical fertilization treatments (MF and DF) showed moderate enrichment (Z-scores ranging from 0.29 to 0.84), with drone application (DF) showing more consistent patterns than backpack application (MF) (coefficient of variation: DF = 18.2%, MF = 27.4%).

Canonical Correlation Analysis (CCA) of treatment groupings based on environmental variables (Figure 5) revealed that the first two canonical axes explained 68 and 17% of the total variance, respectively. Permutation tests (999 permutations) confirmed the significance of the canonical correlations (p < 0.001). Microbial diversity and availability of nutrients, particularly nitrogen (N) and phosphorus (P), were positively associated with biofertilization. The CCA plot shows a clear separation of treatment groups, with drone biofertilization (DB) clustering distinctly from other treatments and associated with higher nutrient availability.

The weighted UniFrac Principal Coordinate Analysis (PCoA) revealed distinct clustering patterns for fungal and bacterial communities across different fertilization treatments. For fungal communities (Figure 6A), the first three principal coordinates explained 99.9% of the total variation (PC1: 95.9%, PC2: 3.3%, and PC3: 0.8%). The no-treatment control was separated from the fertilization treatments, with biofertilization treatments (both drone and conventional applications) forming a distinct cluster. Conventional fertilization and drone fertilization treatments showed intermediate positioning, suggesting a gradual shift in community composition.

The bacterial community structure (Figure 6B) showed a different pattern, with the three principal coordinates explaining 99.8% of the total variation (PC1: 97.6%, PC2: 1.6%, and PC3: 0.6%). The treatments exhibited clear spatial separation, with drone biofertilization and conventional drone fertilization clustering distinctly from conventional fertilization and no-treatment controls. Notably, the biofertilization treatment showed an intermediate position between conventional and drone-based applications, suggesting a gradient of community composition changes influenced by the fertilizer and application methods.

Network analysis revealed complex interactions among microbial guilds in response to drone-assisted biofertilization based on correlated changes in abundance throughout the crop cycle (Figure 7). Bacteria with potential plant growth-promoting abilities showed strong positive correlations among themselves (r > 0.8), particularly between Bacillus, Pseudomonas, and Acinetobacter suggesting a synergistic relationship in their response to biofertilization (Figure 7A). Potential nitrogen-fixing bacteria, represented by Enterobacter and Klebsiella, displayed positive correlations (r ≈ 0.6 to 0.8) with several PGPB (Figure 7A), suggesting a cooperative response to the treatment.

Interestingly, phytopathogenic fungi showed varied responses (Figure 7B). While some, like Pythium, exhibited negative correlations with PGPB (r ≈ −0.7), others, like Botrytis, showed weaker negative correlations (r ≈ −0.4), suggesting differential suppression. Organic matter-degrading fungi, including Mucor and Rhizopus, demonstrated positive correlations (r > 0.7) with both PGPB and nitrogen-fixing bacteria, indicating facilitation of nutrient cycling by these fungi in response to biofertilization (Figure 7). Plant growth-promoting fungi, particularly Penicillium and Trichoderma, showed strong positive correlations (r > 0.8) with PGPB (Figure 7D), suggesting a complementary response to the treatment. The network involving organic matter-degrading bacteria (Figure 7E) suggests a strong positive correlation (r > 0.8) with PGPB and nitrogen-fixing bacteria, indicating an integrated response in nutrient mobilization and plant growth promotion. Environmental fungi displayed varied correlations with other microbial groups (Figure 7F). Some genera, like Naganishia, showed moderate positive correlations (r ≈ 0.5 to 0.7) with PGPB, while others exhibited weak to moderate negative correlations (r ≈ −0.3 to −0.6), suggesting diverse ecological roles in the biofertilized soil.

These results highlight the complex, interconnected response of the soil microbiome to drone-assisted biofertilization and reveal potential synergies and antagonisms that emerged throughout the treatment.