To carry out this study, strains of the bacteria B. velezensis (5D5 and 6E9) and B. subtilis (1A11) were isolated and selected in vitro based on their potential to promote growth in plants, as well as their ability to survive under osmotic stress as described by de Oliveira-Paiva et al. (2024). The microorganisms were isolated from soil samples, starting with the screening of a total of 414 strains samples. In the subsequent stages of evaluation and screening, which included the capacity to grow in medium with a high concentration of sorbitol, assessing exopolysaccharide (EPS) production capacity, biofilm production, and a qualitative evaluation of siderophore and growth phytoregulators, the three strains identified as 5D5 (germplasm code BRM051757), 6E9 (germplasm code BRM051761), and 1A11 (germplasm code BRM051734) were selected for further analysis. The methods and protocols used during this stage of the study are outlined below.

The greenhouse experiments were carried out at Embrapa Maize and Sorghum, located in Sete Lagoas, MG, Brazil, in 2023. The study evaluated the dose-range effect of inoculants individually formulated with the strains 1A11, 5D5, and 6E9 on the growth of soybean seedlings. This assessment was carried out under conditions of adequate moisture availability, as well as during episodes of drought stress.

The inoculants were prepared with bacterial colony-forming unit concentration for each strain as follows: 1A11 = 8.67x10⁹ bacterials mL⁻¹; 5D5 = 1.62x10¹⁰ bacterial mL_-1 to CFU; and 6E9 = 9.41x10⁹ bacterial mL_-1 to CFU. Soybean seeds were treated with doses of 1, 2, 3, and 4 mL kg⁻¹ for all tested strains, with application occurring just before sowing. The soybean seeds were homogenized with the liquid inoculant by mixing them in a container until all seeds were thoroughly moistened and covered by the inoculant. As a standard practice for soybean cultivation to enhance symbiotic nitrogen fixation, it was applied to all treatments the liquid inoculant SimbioseNod™, formulated from strains SEMIA 5079 and SEMIA 5080 of Bradyrhizobium japonicum, at a concentration of 7.2 x 10⁹ bacterial mL_-1 to CFU of inoculant, at a dose of 3 mL kg⁻¹ of seed. The B. japonicum inoculant was applied separately, immediately after the tested Bacillus inoculants were administered. The control group received only the B. japonicum and this inoculant was mixed with the seed mass in a manner similar to the process used for the Bacillus inoculants.

For each experimental plot, 20 L pots filled with clayey red latosol, typical of the Brazilian Cerrado biome, were utilized. Soil fertilization and liming were consistent across all treatments, adhering to recommendations based on soil physicochemical analysis. Two soybean plants of the cultivar ‘BRS Valiosa RR’, were maintained in each pot. The experimental design was completely randomized, comprising two sets of 26 treatments with three replicates per treatment. The data were analyzed by the R statistical package.

Of the 78 experimental units, 39 were sustained under adequate irrigation, with soil field capacity continuously maintained at -18 kPa. The other 39 experimental units underwent drought stress conditions starting 27 days after germination, simulated by keeping the soil at a water tension of -138 kPa for 14 days. Daily monitoring of soil water availability for both irrigated and non-irrigated treatments was conducted using a tensiometer installed at a depth of 20 cm in each pot. At the conclusion of the drought stress period, the irrigated and non-irrigated plants were collected at the early flowering stage (R1), and the weights of fresh and dry matter for both the aerial parts and roots were evaluated.

To assess the impact of microbiological inoculants on soybean crops, field-trials were established in three different locations in Brazil: Birigui (geographical coordinate -21.318220, -50.355165), Itapira (-22.401927, -46.778187) and Piracicaba (-22.769027, -47.582643) and evaluated over two consecutive growing seasons (2022–2023 and 2023–2024). The soybean varieties was chosen based in its agronomic recommendation for each region and in its precocity of the cycle. In Birigui was cultivated the variety ‘BMX Nexus i2X’ of the maturity group 6.4, in Itapira was cultivated the variety ‘Monsoy 6601 i2X’ of the maturity group 6.6, and in Piracicaba was cultivated the variety ‘NS 5933 iPro’ of the group of maturity 6.1. To establish a stand of approximately 260,000 plants per hectare, about 13 seeds were sown per linear meter of furrow, maintaining a spacing of 50 cm between rows and the planting depth was set at 5 cm. During the planting of seeds was applied granulated fertilizer with a composition of 3-21-21 (NPK, YaraBasa™) at rates based on the soil’s chemical analysis for each location: 150 kg ha⁻¹ in Birigui; 130 kg ha⁻¹ in Itapira; and 180 kg ha⁻¹ in Piracicaba and the fertilizer was positioned next to the seed, 5 cm below it. During the phenological phases between V3 and V6, a foliar micronutrient fertilizer (Micromax Foliar Concentrate™, Aqua do Brasil) was applied at a rate of 500 g ha⁻¹.

The field-trials of Birigui and Itapira consisted of 16 experimental plots for each location and twenty rows of soybean measuring 0.5 m between furrows and 100 m in length, while the field-trial in Piracicaba also consisted of 16 experimental plots but with 10 rows of soybean measuring 0.5 m between furrows and 25 m in length. For each soybean variety, there were 4 blocks per field-trial, totaling 48 experimental plots per season and 96 experimental plots considering two harvest seasons. The random block design with four blocks per field-trial was adopted in this study to enable the measurement of random environmental variation. All experimental fields followed management practices aligned with standard commercial soybean cultivation in Brazil, ensuring consistent agronomic standards across all varieties. Soybean plots were harvested and measured in tons of soy per hectare (t ha⁻¹) at the end of their growth cycle, once the seeds reached peak of maturation (13% of humidity based on the dry weight). Yield was assessed based on the average grain weight obtained from experimental replicates for each treatment across all locations (Birigui, Itapira and Piracicaba) over two consecutive growth cycles: 2022–2023 and 2023–2024.

To assess the impact of three microbiological inoculants on soybean crops, each bacterial strain—5D5, 6E9, and 1A11—formulated as liquid inoculant was applied individually at a rate of 3 mL of liquid inoculant per kilogram of soybean seeds at the same concentration of the greenhouse study. A control treatment without inoculant application was included for comparison.

The inoculants containing strains 5D5, 6E9, and 1A11 were applied to the seeds in a similar manner to the procedure described in the greenhouse experiment. After application, the seeds were allowed to dry for five minutes before sowing, making sure to avoid direct sunlight exposure. Additionally, after the Bacillus inoculants application, B. japonicum inoculant was applied to the seeds of all treatments, including the control, using the same commercial product at the same concentration and dosage as in the greenhouse experiment.

The climate data for the field-trials in Birigui, Itapira and Piracicaba were compiled from the NASA Power (NASA POWER Project, 2024) database as indicated by Monteiro et al. (2018). The dataset, which spans the entire soybean growing season for each location and harvest season, was sourced directly from the NASA POWER Project (2024).

The EVI2 VI (Jiang et al., 2008) was calculated for all field-trials during the phenological phase of vegetative growing. This VI can be used as an alternative to the traditional NDVI (Normalized Difference Vegetation Index) without changing the spectral bands used in its calculation (red and near infrared), which often reach saturation points in the most advanced stages of crop development, making it difficult to observe differences between treatments (Santos et al., 2024). For the fields of Birigui and Itapira were used satellite images from the PlanetScope PSB.SD instrument (Frazier and Hemingway, 2021), provided by RedeMAIS/MJSP source, with spatial resolution of 3 m/pixel and daily temporal resolution, regarding the following bands: blue (B): 465–515 nm; green (G): 547–583 nm; red (R): 650–680 nm; red edge (RE): 697–713 nm; and near-infrared (NIR): 845–885 nm. For the field of Piracicaba, due to the smaller size of the experimental area, were used images obtained from missions performed by a drone, model DJI Phantom 4 Pro equipped with RTK and two cameras: one RGB, for images in the visible spectrum; and another for multispectral array, covering the following bands: blue (B): 450 ± 16 nm; green (G): 560 ± 16 nm; red (R): 650 ± 16 nm; red edge (RE): 730 ± 16 nm; and near-infrared (NIR): 840 ± 26 nm. In this case, for each mission, orthomosaics were generated by the Agisoft Metashape software, considering the aforementioned bands with a ground sample distance (GSD) of approximately 10 cm. The satellite images or the orthomosaic assembled from drone images were then used to extract and calculate the VIs of the field-trials. For both cases, only the images or othomosaic that best represented the phenological stadium of vegetative growth (approximately 70 to 90 days) were effectively chosen to extract the EVI2 values. The average EVI2 value of the internal pixels of each experimental plot, excluding the pixels considered as borders, was then calculated and considered in this study.

In order to enhance the prediction model’s stability and address potential multicollinearity concerns, all predictor variables were incorporated into the initial analysis. The assessment of the Variance Inflation Factor (VIF) (Belsley et al., 1980) was performed to detect and eliminate variables with high collinearity, leading to a more streamlined and reliable selection of variables. The Shapiro-Wilk test was used for normality (Shapiro and Wilk, 1965) and the Breusch-Pagan test for homoscedasticity (Breusch and Pagan, 1979).

To investigate the determinants of soybean productivity, we define a set of explanatory variables that comprise both categorical and continuous covariates. The response variable is the grain yield, denoted by y_i , measured in tons per hectare (t ha⁻¹). Each observation i = 1,…,96 corresponds to a unique combination of experimental conditions. The variables are structured as follows:

Categorical variables (x_{i 1} to x_{i 4}) were transformed into binary indicators through dummy encoding to ensure compatibility with linear modeling frameworks. Continuous variables (x_{i 5} and x_{i 6}) were retained in their original form to preserve their numeric information. This structure enables the application of linear regression models and facilitates the interpretation of individual effects on soybean yield.

A linear regression model was used to evaluate the effect of explanatory variables on soybean grain-yield, measured in tons per hectare (t ha⁻¹). The model includes both categorical and continuous covariates previously described and is given by:

where y_i represents the response variable in observation i. The term β ₀ corresponds to the model intercept, which indicates the value of the response variable when all the explanatory variables are zero. The coefficients β_j are the regression parameters associated with the explanatory variables x_ij , for j = 1,2,…,p, which quantify the influence of each of these variables on the response variable. The variables x_ij are the explanatory variables that can directly influence the response variable. Finally, the term ϵ_i represents the random error associated with observation i, which captures the variations unexplained by the variables included in the model. It is assumed that the error ϵ_i is a random variable with a normal distribution, mean zero, and constant variance, i.e., ϵ_i ∼ N(0,σ ²). Thus, the equation expresses how the response variable is affected by the explanatory variables, incorporating an error component that reflects the unexplained variability.

The model was fitted using the lm() function from the R software environment (R Core Team, 2024), which estimates the regression coefficients using the Ordinary Least Squares (OLS) method. This method aims to minimize the sum of squared residuals:

where y ^ i is the fitted value for observation i. In matrix form, the estimator is given by:

where X is the design matrix, y is the response vector, and β ^ is the vector of estimated coefficients. The formulation and estimation are based on standard theory for linear models as described in Chambers (1992).

The selection of Bacillus strains from the Active Germplasm Bank of Multifunctional and Phytopathogenic Microorganisms at Embrapa Maize and Sorghum focused on laboratory tests to assess their ability to grow in vitro under high osmotic potential conditions ( Table 1 ). This included evaluating colony growth in a medium with high sorbitol concentration and low water activity. Additionally, the selected strains were examined for mechanisms that may help mitigate plant drought stress, such as the ability to exude biofilm and exopolysaccharides. Some characteristics associated with growth promotion of roots and canopy, including the synthesis of indole-3-acetic acid (IAA) and siderophores, were also evaluated. Based on the results obtained from strains 5D5 (B. velezensis), 6E9 (B. velezensis), and 1A11 (B. subtilis), the next step was to evaluate the behavior of these three strains in a biological system involving both plant and soil in a controlled growth environment. For this purpose, soybean was selected due to its agronomic significance and because it is a crop with sensitivity to moderate and extreme variations in temperature and soil humidity—conditions that are commonly observed in many agricultural systems in Brazil.

The identification of microorganisms that can thrive in low humidity and simultaneously produce biologically active substances holds significant promise for various applications, including enhancing plant tolerance in marginal environments. The bacterial strains were initially isolated from arid soils and the strains capable of growing in culture medium containing sorbitol at concentrations of 520 g L⁻¹ and/or 780 g L⁻¹ proceeded to the next stages of selection and evaluation. Further, a series of screenings was carried out to identify multifunctional mechanisms for plant growth-promoting (PGP) related to drought stress tolerance, particularly those related to exopolysaccharide (EPS) and biolfim production and capacity of growing in a low water activity. Water activity (aw) refers to the amount of water available for biological processes, and a value of 0.807, where 5D5 and 6E9 strains could grow, indicates a reduced moisture availability in the soil. For example, most microbes cease cell division below 0.9, and few prokaryotic strains have been shown to continue cell division at aw smaller than 0.755 (Amina and Lotfi, 2024). However, the strains 5D5, 1A11, and 6E9 screened in water-restricted culture media were efficient in producing EPS in vitro.

EPS are high-molecular-weight carbohydrates attached to the external surface of bacteria, related to the formation of biofilms and the attachment of bacterial cells to surfaces, including plant roots and soil particles. Since EPS are hydrated compounds, including 97% water in a polymeric matrix, they can enhance growth and ensure the survival of plants under drought stress (Vu et al., 2009; Nocker et al., 2012). In addition, they can protect the plant against desiccation due to the formation of hydrophilic biofilms on the root surface (Rossi et al., 2012). Although EPS is the main constituent of the biofilm, the 1A11 strain did not show biofilm production in vitro. Since the presence of EPS and biofilm are closely related, the difference between the number of bacteria producing these two compounds can be explained by the need for other substances to be present for biofilm formation, such as extracellular proteins or molecules related to hydrophobicity, which makes it more or less soluble (Vlamakis et al., 2013; Zhang et al., 2015). The production of EPS by bacteria is important because it increases soil permeability and aggregation, maintaining high water potential near the roots (Alami et al., 2000), which facilitates the diffusion of nutrients and the mass flow of soluble substances, increasing plant absorption and therefore, promoting the mitigation of water deficit. Inoculation experiments of maize seeds with EPS-producing bacterial strains, in combination with their respective purified EPS, resulted in increased soil moisture content, plant biomass, root and shoot length, and leaf area (Zhang et al., 2015). The microenvironments provided by EPS retain water and dehydrate more slowly than the surrounding environment, thus protecting the rhizosphere from desiccation.

The Bacillus isolates in this study also produced Indole-3 acetic acid (IAA), ranging from 14.19 to 29.12 µg mL⁻¹. IAA is a member of the group of phytohormones normally considered the most important native auxins. It has been frequently reported as an important mechanism for plant growth promotion, improving root size and distribution, and ensuring an increase in nutrient uptake from the soil. The production of IAA by bacteria differs between species and can vary between strains, but it is also influenced by culture conditions, availability of vitamins, salts, oxygen, pH value, temperature, growth rate, and available nitrogen source (Duca et al., 2014), while the tryptophan, an amino acid used as an IAA precursor, can be provided by the plant roots.

To assess the ability of B. velezensis strains 5D5 and 6E9, and B. subtilis strain 1A11 to enhance growth in soybean plants in comparison to a control inoculated only with B. japonicum, an experiment using pots were carried out in a greenhouse with controlled environmental conditions. The experiments consisted of dose-range of the three inoculants, evaluating the effects of 1, 2, 3, and 4 mL of inoculant per kg of seeds. The experiment was divided into two groups: one with adequate moisture for optimal plant development and the other subjected to water restriction, simulating a drought episode.

Table 2 displays the response of soybean plants in terms of Shoot Fresh Weight (SFW), Shoot Dry Weight (SDW), Root Dry Weight (RDW), Number of Nodules per Root (NNR), and Mass of Functional Nodules per Root (MFNR) under irrigated and non-irrigated conditions. In both conditions, it was observed that the average treatment effects for certain inoculant doses were higher than the average of the control groups (T13 and T26 for irrigated and non-irrigated, respectively).

Under irrigated conditions (T1 to T13), treatment T4 (1A11, 4 mL·kg⁻¹) exhibited the highest mean shoot fresh weight (SFW), not differing significantly from treatments T5 (6E9, 1 mL·kg⁻¹), T6 (6E9, 2 mL·kg⁻¹), T7 (6E9, 3 mL·kg⁻¹), T8 (6E9, 4 mL·kg⁻¹), T9 (5D5, 1 mL·kg⁻¹), and T10 (5D5, 2 mL·kg⁻¹). All these treatments outperformed the non-inoculated control (T13). A similar pattern was observed for shoot dry weight (SDW) and root dry weight (RDW), where the same treatments, except T5 for SDW, showed significantly higher values compared to the control. In contrast, no significant improvements were observed in the number of nodules (NNR) and the mass of functional nodules (MFNR) in inoculated treatments compared to the control (T13).

Under non-irrigated conditions (T14 to T26), for the variable shoot fresh weight (SFW), the treatments T17 (1A11, 4 mL·kg⁻¹), T21 (6E9, 4 mL·kg⁻¹), and T22 (5D5, 1 mL·kg⁻¹) exhibited the highest shoot fresh weight (SFW) means, significantly surpassing the non-inoculated control (T26). Similarly, treatments T18 (6E9, 1 mL·kg⁻¹), T23 (5D5, 2 mL·kg⁻¹), T24 (5D5, 3 mL·kg⁻¹), and T25 (5D5, 4 mL·kg⁻¹) showed intermediate performance, but still significantly higher than the control. For shoot dry weight (SDW), treatments T17 (1A11, 4 mL·kg⁻¹), T18 (6E9, 1 mL·kg⁻¹), T21 (6E9, 4 mL·kg⁻¹), and all doses of 5D5 (T22 through T25) recorded superior means compared to the control. In contrast, no significant differences were observed for root dry weight (RDW), suggesting that root development was uniformly constrained under water-limited conditions. Regarding the number of nodules (NNR), none of the treatments outperformed the control. However, for the mass of functional nodules (MFNR), treatment T18 had the highest mean, statistically comparable to T17, T22, and T23, all of which were significantly superior to the control.

The results revealed a strain- and dose-dependent effect of inoculation on plant growth and nodulation parameters, with notable distinctions between irrigated and non-irrigated conditions. Under irrigation, strains 1A11 and 6E9 significantly enhanced shoot fresh weight (SFW) and the mass of functional nodules per root (MFNR), particularly at moderate to high doses. Strain 1A11 reached the highest SFW (72.07 g) at 4 mL·kg⁻¹, while 6E9 showed consistent performance across doses. In contrast, strain 5D5 exhibited reduced performance at higher doses, suggesting a possible inhibitory effect, which may be attributed to microbial competition or metabolic burden on the host plant (Etesami, 2025). It is also important to highlight that each microbial strain possesses distinct mechanisms for promoting drought tolerance, such as the production of phytohormones, enzymes, and biofilm formation. Additionally, the applied dose influences the number of viable cells and spores available to colonize the seed, thereby affecting both the survival and colonization efficiency of the inoculant. Notably, very high or very low inoculum doses may not necessarily benefit the plant, as the concentration of bioactive compounds, such as phytohormones, may become either excessive or insufficient, potentially leading to suboptimal plant responses (Etesami, 2025).

Under non-irrigated conditions, although the overall performance of inoculated plants was reduced compared to irrigated treatments, certain inoculant-strain combinations still demonstrated potential to confer drought resilience. For instance, 6E9 at 1 mL·kg⁻¹ improved SFW (35.10 g) and MFNR (0.65 g) relative to the non-inoculated control (27.38 g and 0.27 g, respectively), and 5D5 at the same dose yielded the highest SFW (36.42 g) under drought. These results suggest that low-dose inoculation with well-adapted strains can mitigate drought stress effects, likely through improved symbiotic efficiency and modulation of root architecture, as previously documented in studies on plant-microbe interactions under water-limited conditions (Timmusk et al., 2014; Ngumbi and Kloepper, 2016; Vurukonda et al., 2016). In contrast, 1A11 enhanced SFW (38.03 g) and SDW (8.71 g), relative to the non-inoculated control under drought condition (27.38 g and 6.48 g, respectively).

Among the evaluated variables, SFW, SDW, and MFNR were the most positively responsive to inoculation, highlighting their utility as indicators of inoculant performance under contrasting moisture regimes. These findings reinforce the importance of selecting compatible strain-dose combinations and support the strategic use of plant growth-promoting bacteria to enhance legume productivity in both irrigated and drought-prone environments.

Abrupt and extreme climate variations have become one of the main challenges for agricultural practices, particularly in Brazil, a tropical country with a strong agricultural tradition. Regular planting and cultivation seasons are changing due to shifts in rainfall patterns, as well as variations in average maximum and minimum temperatures (Dos Santos et al., 2023). These changes increase the level of risk and make agricultural activities more unpredictable.

To mitigate these effects and identify their occurrence, farmers are adopting technological alternatives as part of their strategy to reduce losses and ensure economic returns. In this context, microbiological inoculants that improve plant resilience during mild to moderate water deficits—such as the bacterial strains evaluated in this study—are fundamental for maintaining crop grain yield. Some products are already available on the market, including those aimed at increasing phosphorus solubilization (de Oliveira-Paiva et al., 2024) and microbial-based solutions to reduce the effects of drought (Kavamura et al., 2013).

Additionally, incorporating technologies like remote sensing with time-series satellite or drone images enables producers to track in advance the environmental impacts on crop development. This allows them to take proactive palliative or corrective measures.

During the field-trials with soybeans, significant variations were observed in rainfall patterns and temperature peaks ( Figure 1 ). In the 2022–2023 season, rainfall was more regularly distributed and of a greater volume compared to the 2023–2024 season. The maximum temperatures recorded during 2022–2023 generally stayed below 35°C, whereas during 2023–2024, temperatures frequently exceeded 40°C, negatively affecting the development of soybean crops.

As a result, the 2022–2023 harvest yielded an average of nearly one ton more than the 2023–2024 harvest average across the three locations involved in the study ( Figure 2c ). This clearly indicates the detrimental effects of climate extremes on soybean production. The analysis also revealed that Itapira had slightly higher average grain yield than Birigui and Piracicaba, likely due to better rainfall averages, especially durin7g the 2023–2024 season ( Figure 2b ). Beyond climate variations across locations ( Figure 1 ), the differences in environmental conditions—particularly soil physicochemical properties—and the distinct genetic backgrounds of the three soybean varieties may also influence the results. This is worth mentioning that the three soybean varieties utilized in this study are categorized into maturity groups ranging from 6.1 (‘NS 5933 iPro’) and 6.4 (‘BMX Nexus i2X’) to 6.8 (‘Monsoy 6601 i2X’). Therefore, these varieties are suitable for the southeastern region of Brazil, where the field trials were carried out, and it also provided a wider planting window for the locations. This flexibility allowed to plant later in the season, exposing the crops to a period of lower rainfall and reduced soil humidity, and thus evaluate the effect of water shortage on the Bacillus/soybean crop association under conditions that are appropriate to water deficit.

Figure 2a shows that microbial strain 1A11 achieved the highest soybean yield across the two growing seasons and three tested locations, producing an average of 620 kg ha⁻¹ more than the control (an 18% increase). Strains 5D5 and 6E9 also led to yield improvements over the control, with 5D5 yielding 390 kg ha⁻¹ (11.3% increase) and 6E9 yielding 440 kg ha⁻¹ (12.8% increase). These results suggest that the application of these treatments resulted in yield gains relative to the control.

Among the locations, Itapira achieved the highest average yield (4.66 t ha⁻¹), followed by Birigui (3.69 t ha⁻¹) and Piracicaba (3.08 t ha⁻¹), as shown in Figure 2b . Furthermore, significant differences in yield were observed between the two growing seasons, with 2022–2023 (4.28 t ha⁻¹) outperforming 2023–2024 (3.34 t ha⁻¹), as shown in Figure 2c . According to Figure 2 , increased cumulative precipitation ( Figure 2d ) and higher EVI2 ( Figure 2e ) values are associated with indications of greater soybean grain yield, suggesting that these variables may have a positive influence on crop yield.

Table 3 summarizes the mean and standard deviation of grain yield (t ha⁻¹) and the EVI2 index for different locations, treatments, and years. The data cover two seasons (2022–2023 and 2023–2024) and three locations: Birigui, Itapira, and Piracicaba. For each location and year, values for four treatments—1A11, 5D5, 6E9, and control—are provided. In 2022–2023 harvest, grain yield means ranged from 3.58 t ha⁻¹ in the control treatment at Birigui to 5.14 t ha⁻¹ in treatment 1A11 at the same location. EVI2 values were relatively stable across all treatments and locations in 2022–2023 season, with mean values close to 0.95. In 2023–2024 harvest, a noticeable decrease in grain yield was observed at most locations and treatments, particularly at Piracicaba, where the control treatment showed the lowest grain yield (1.78 t ha⁻¹). EVI2 values remained consistent across treatments, with some variations, especially in Itapira and Birigui. These results reflect both grain yield levels and vegetation dynamics (as indicated by EVI2) across the two seasons.

Additionally, Table 3 highlights that treatments with strains 1A11, 5D5, and 6E9 performed better in terms of EVI2 compared to the control. However, in Birigui during the 2022–2023 season, all treatments showed the same average EVI2 values, which were close to the maximum of the EVI2 scale (0.95), suggesting potential saturation under these specific conditions.

The high EVI2 values in Birigui during 2022–2023 may be attributed to the greater canopy development of soybean plants in all treatment groups, compared to other locations or years. Environmental variations in open field conditions likely contributed to these discrepancies.

Based on the multicollinearity analysis, the explanatory variable Year (x_{i 4}) was excluded due to its high VIF in relation to the other covariates. The final model, without this variable, was fitted to represent the relationship between the dependent variable y_i and the selected explanatory variables. The model equation is given by:

This revised model, excluding the variable with high multicollinearity, results in a more robust and precise selection of explanatory variables.

The results of the statistical model ( Table 4 ) indicated that treatments 1A11, 5D5, and 6E9 showed significant differences compared to the control (represented by the intercept), with statistical significance determined by p-values less than 0.10, however, these treatments did not differ from each other ( Table 5 ). These findings highlight the efficacy of strains 1A11, 5D5, and 6E9 in enhancing soybean development under non-irrigated, open-field conditions. The observed grain yield gains were statistically significant compared to the not inoculated control, underscoring the potential of microbial inoculants to improve crop performance even in challenging agricultural environments.

Regarding location, Itapira showed significantly higher average grain yield compared to the other sites, as illustrated in Figure 2b . The bar plot indicates that the mean yield in Itapira was greater than in other locations, particularly when compared to Birigui. This difference may reflect environmental and management variability between locations and the availability of rainwater throughout the crop cycle, with the rainfall regime being more favorable in some of the locations where the tests were carried out than in others ( Figure 1 ).

As for cumulative precipitation, the model results indicate that this variable was statistically significant. However, excessive precipitation appears to be associated with a slight reduction in soybean grain yield. This trend is evident in Figure 2d , which shows a small decrease in grain yield at the highest levels of accumulated rainfall during the early stages of the experiment.

Finally, the EVI2 index showed a significant positive effect on grain yield: according to the fitted model, for each one-unit increase in EVI2, soybean yield increased by an average of 4.45 t/ha. It is worth noting that, as the crop approaches its final growth stages, a natural decline in vegetation indices is expected due to plant maturation and senescence.

The half-normal plot (Moral et al., 2017) with simulated envelope ( Figure 3a ) showed that all standardized residuals fell within the envelope. This suggests that the linear model provides an adequate fit to the data, with residuals behaving as expected under the assumptions of normality and homoscedasticity.

The Shapiro-Wilk test yielded W = 0.9898 with a p-value of 0.6756, and the Breusch-Pagan test produced BP = 12.508 with df = 10 and a p-value of 0.2525. Both tests were non-significant, indicating that the residuals are normally distributed and exhibit constant variance, thus supporting the adequacy of the linear model.

Additionally, the model achieved a coefficient of determination (R ²) of 0.71, indicating that approximately 71% of the variability in soybean grain yield was explained by the included explanatory variables.

Figure 3b illustrates the comparison between observed and predicted grain yield values obtained from the multiple linear regression model. The close alignment between observed and predicted points demonstrates that the model effectively captures the variability in soybean grain yield. In the figure, blue triangles represent the observed values, while red circles indicate the predicted values. The overlapping of these points reflects the natural variability in the data and supports the model’s accuracy in predicting grain yield.

Research utilizing remote sensing tools to study microbiological inoculants is surprisingly limited in the literature, revealing a significant gap but also presenting opportunities to improve the assessment of these inoculants’ effectiveness on agronomic crop performance. For example, de Souza et al. (2022) conducted a study on common beans, using UAVs to evaluate and monitor the biofertilization process. They applied VIs to map spectral changes in common beans after the application of phosphate-solubilizing bacteria in a remote and non-destructive manner. Similarly, Santos et al. (2024) investigated phosphate-solubilizing bacterial inoculants in sugarcane crops developed a predictive statistical model that accurately forecasted sugarcane grain yield in fields treated with these microbiological inoculants, producing predictions months ahead of harvest that closely aligned with actual field data.