The organic fertiliser ORGAON® PK (OPK), previously described by Robas-Mora et al. (2023), was prepared from the leachate of horticultural residues supplemented with phosphorus (0.10 ± 0.03% w/w) and potassium (2.58 ± 1.52% w/w). The physicochemical properties of OPK are summarised in Supplementary Table S1. The product was applied at a 1: 512 dilution (OPK: H₂O).

To assess the effect of the organic matrix in the absence of its native microbiota, a sterilised variant (OPK_ST) was produced. Fifty millilitres of OPK were exposed to UV-C radiation (265 nm, Philips TUV 30 W G30T8) for 15 min in a laminar-flow cabinet; the fertiliser was distributed in Petri dishes to a depth < 5 mm to ensure uniform irradiation. Pre- and post-irradiation aliquots were analysed for humic extract, fulvic acids, DOC and inorganic N/P (Supplementary Table S2). No significant differences were detected (t-test, p > 0.05), confirming that the 15-min UV-C protocol removed viable microbiota while leaving the chemical profile essentially unchanged (Supplementary Table S3). Sterility was verified by plating 100 μL of the treated solution, in duplicate, on LB agar (Condalab®, Madrid, Spain) and incubating at 28 ± 1 °C for 72 h. Only samples showing no microbial growth in either replicate were deemed sterile. Sterilised solutions were stored at 4 °C and used within 24 h of preparation.

Prior to field application the residue-microbiome safety was screened. The crude residue (OPK) was assayed for potential human or environmental pathogens. Triplicate 10-mL aliquots were pelleted (10,000 × g, 10 min), the pellets were processed with the DNeasy PowerSoil Pro kit (QIAGEN) and the V3–V4 16S-rRNA region was sequenced as described above. After DADA2 denoising, amplicon sequence variants were classified against SILVA 138.99; additionally, all quality-filtered reads were mapped with Kraken2 (v2.1.3, max-hits = 1) against the NCBI RefSeq “bacteria, archaea & viruses” database to capture low-abundance taxa. Taxa were flagged as putative pathogens if they belonged to any of the following genera: Escherichia/Shigella, Salmonella, Klebsiella, Enterobacter, Citrobacter, Listeria, Staphylococcus, Pseudomonas aeruginosa group, Clostridium sensu stricto 1, Vibrio or Mycobacterium tuberculosis complex. None of these genera exceeded 0.01% of total reads in any replicate (Supplementary Table S1).

The OPK microbiome was dominated by soil saprotrophs—Rubrobacter (18%), Nocardioides (9%), Streptomyces (7%) and Sphingomonas (4%)—plus Planctomycetota lineages typical of lignocellulose-rich residues (e.g., Pirellula, 3%). Absolute read counts for the 25 most abundant genera are shown in Supplementary Figure S1; together they accounted for 72% of the library, indicating moderate evenness. The combined Kraken2 screen and genus-level inspection confirmed the absence of clinically relevant pathogens at detectable levels, supporting the biosafety of OPK for open-field use.

The inorganic reference fertiliser was a commercial NPK solution (Universal Fertilizer Complet®, COMPO Iberia S. L., Barcelona, Spain). Its guaranteed analysis is 7–5-6 (N-P₂O₅-K₂O) with micronutrients: 7% total nitrogen (2.1% nitrate-N, 1.3% ammoniacal-N, 3.6% urea-N), 6% water-soluble phosphorus (P₂O₅) and 5% water-soluble potassium (K₂O). The solution also supplies (all water-soluble): 0.01% boron (B), 0.002% copper (Cu), 0.02% iron (Fe), 0.01% manganese (Mn) and 0.002% zinc (Zn), each chelated with EDTA, together with 0.001% molybdenum (Mo).

The biofertiliser was prepared by enriching OPK with the selected PGPB strains (B. pretiosus—C1 or P. agronomica—C2) at the dilution recommended for fertigation of M. sativa (Robas-Mora et al., 2023). Twenty millilitres of a bacterial suspension adjusted to McFarland 4 were added per liter of OPK, yielding ≈ 10⁸ cfu mL⁻¹ at 10⁻¹ dilution. Fresh biofertiliser was prepared weekly and stored at 4 °C until use. The same inoculation protocol was applied to CF and to the water control (W).

Each fertiliser was applied once, as the initial watering at the start of the trial. Thereafter, plants were maintained under rain-fed conditions from September 2022 to May 2023; no supplementary irrigation was provided. Two inoculations were performed, one at the beginning of the experiment in September and the second one before spring in February.

To avoid confounding factors such as rodent activity, wind or farm traffic, harvested biomass was dried ex situ. Material from each subplot was transferred to a secure shed and spread in single layers, preserving the spatial order of the field layout to retain sample traceability. Stacking was avoided to prevent fermentation and ensure uniform desiccation. Drying proceeded under ambient conditions (20–25 °C, natural ventilation, no direct sunlight) for 18 days. Weight stability was checked on days 12, 14, 16 and 18 by weighing three representative samples per treatment every 48 h; a change of < 1% between consecutive measurements was taken as the end-point. Total biomass per subplot was recorded on a platform balance (capacity 20–150 kg, accuracy ± 0.1 kg), noting both fresh and final dry weights to calculate moisture content and dry-matter yield.

Within 24 h of harvest, samples were analysed at Rock River Labs Spain (Lalín, Pontevedra) using near-infrared spectroscopy (NIRS DS2500, FOSS Analytical, Denmark). Plant material was dried, separated into leaf and stem fractions, and milled to 1 mm with a Cyclotec mill (FOSS). Spectra were collected from 400 to 2,498 nm and processed with proprietary calibrations validated by Rock River Laboratory. The following variables were determined: water-soluble carbohydrates (WSC), undigested neutral-detergent fibre after 240 h (uNDF₍₂₄₀₎), total tract fibre digestibility (TTFDND), soluble protein, crude protein, indigestible modified NDF fraction (aNDFₘₒ), and total amino acids. Each measurement was performed in triplicate to ensure robustness. Statistical analyses were conducted in R v.4.3.1. One-way ANOVA tested treatment effects; post-hoc differences were examined with Tukey’s HSD; effect sizes were expressed as Cohen’s d. Multivariate patterns were explored by principal-component analysis (packages stats, factoextra, ggplot2 and psych).

To evaluate antibiotic resistance, the supernatant was adjusted to an optical density equivalent to 0.5 McFarland with sterile saline (0.45% NaCl). Aliquots of 100 μL were spread aseptically onto Mueller–Hinton agar (Condalab®, Madrid, Spain) using sterile swabs. Etest® strips (bioMérieux®, France) were applied to determine community minimum inhibitory concentrations (MICs) for the following antibiotics: amoxicillin (AML), amoxicillin + clavulanic acid (AUG), cefotaxime (CTX), piperacillin (PP), cefepime (PM), piperacillin + tazobactam (TZP), imipenem (IMI), imipenem + EDTA (IMD), trimethoprim (TS), gentamicin (CN), nalidixic acid (NA) and ciprofloxacin (CIP). Plates were incubated at 30 °C for 24 h, or as recommended by the manufacturer.

MICs were read visually at the point where the inhibition ellipse intersected the strip’s scale (González-Reguero et al., 2023; Robas-Mora et al., 2023). Prior to statistical testing, MIC values were log₂-transformed to stabilise variances and approximate normality. A general linear model (GLM) ANOVA (p < 0.05) was used to detect significant differences among biological treatments (C0, C1, C2). Multivariate patterns were explored by principal component analysis (PCA) and correspondence analysis (CA), complemented with hierarchical clustering (Euclidean distance) and heat-map visualisation. A correlation matrix of MIC values was computed to identify potential cross-resistance. Post-hoc pairwise comparisons were performed with Tukey’s HSD test. All statistics were run in R (R Core Team, 2023).

Functional diversity of the rhizosphere communities was determined with Biolog EcoPlates® (Biolog Inc., Hayward, CA, USA), each containing 31 carbon sources in triplicate plus a negative control. The bacterial extract described above was adjusted to 0.5 McFarland with sterile saline (0.45% NaCl); 150 μL of this suspension were dispensed into each well, and plates were incubated at 25 °C for 7 days. Optical density (OD₅₉₅) was recorded every 24 h with a Multiskan FC microplate reader (Thermo Fisher Scientific). For each reading the corresponding blank value was subtracted, and triplicate means were calculated for each substrate. Average Well Colour Development (AWCD) was plotted against incubation time to obtain metabolic-activity curves for each community, and the time point of maximum AWCD was used to compute functional-diversity indices.

Metabolic diversity was expressed as the Shannon–Weaver index (Hₘ): Hm = −∑ qi log₂(qi); where qi is the proportion of corrected absorbance for well i: qi = Ai / ∑Ai; with Ai being the blank-corrected absorbance of well i and ∑Ai the summed absorbance for the plate at that time point. The index was calculated for each biological replicate, and the resulting values were subjected to statistical analysis to detect significant differences among treatments.

Total DNA was extracted from rhizosphere soil using the QIAsymphony PowerFecal Pro kit (QIAGEN®, Germany) in accordance with the manufacturer’s protocol. The V3–V4 region of the 16S rRNA gene was amplified with universal primers 341F (5′-CCTACGGGNGGCWGCAG-3′) and 805R (5′-GACTACHVGGGTATCTAATCC-3′). Libraries were prepared following Illumina’s standard workflow and quantified with the Quant-iT™ PicoGreen™ dsDNA Assay (Thermo Fisher Scientific). Equimolar pools were sequenced on an Illumina MiSeq® platform (2 × 250 bp), including 20% PhiX as an internal control.

Raw reads were processed in QIIME2 (version <insert>): adaptors and low-quality ends were trimmed, noise was removed with DADA2, and amplicon-sequence variants (ASVs) were inferred. Taxonomic assignment was performed with VSEARCH against the SILVA 138 99 database. The output comprised a BIOM ASV-frequency table and a rooted phylogenetic tree.

To characterise within-sample diversity, the following α-diversity metrics were calculated: Shannon index, Faith’s Phylogenetic Diversity, Pielou’s evenness and the number of observed taxa.

For partitioning of β-diversity, Jaccard dissimilarity was decomposed into turnover and nestedness components with the betapart package (v1.6). Pairwise distances from each treatment to the water control (WC0) were summarised in box-plots and compared by Kruskal–Wallis (α = 0.05). To gauge inoculant persistence, per-sample ASV counts assigned to Bacillus or Pseudomonas were summed and plotted against total reads. No further scaling was applied, as sequencing depth was normalised in previous steps. Differential abundance (ANCOM-BC). Bias-corrected log-fold changes were computed with ANCOM-BC v2.2 (ancombc) using the pseudo-count ADD = 1, structural zero detection ON and FDR adjustment (method = “BH”). Analyses were run separately at family and genus rank; species-level results were ignored owing to low 16S resolution. Only taxa with FDR-adjusted p < 0.05 were exported. Log₂FC and standard errors were visualised with ggplot2 (Figure 1). Absence of a panel denotes no significant taxa at that taxonomic level. All R analyses were scripted in R v4.3.2 under Ubuntu 22.04.

A total of 730,983 16S reads (mean length 301 bp; ≥ 14,171 valid sequences per sample; BioProject PRJNA 1162555) were obtained, and the rarefaction curves plateaued, confirming adequate sequencing depth.

Shannon indices ranged from 8.9 to 9.3, with no significant differences among biological treatments (C0 vs. C1 vs. C2); Mann–Whitney tests yielded p ≥ 0.32 and q ≥ 0.70 (Supplementary Table S3). Thus, inoculation with B. pretiosus or P. agronomica did not alter within-sample diversity.

Three pairwise comparisons showed no change (CF vs. OPK, CF vs. Water, OPK vs. Water; p ≥ 0.39), whereas the sterilised residue (OPK_ST) displayed a moderate decline in diversity: CF vs. OPK_ST, p = 0.02, q = 0.07; OPK vs. OPK_ST, p = 0.008, q = 0.03; OPK_ST vs. Water, p = 0.002, q = 0.01. These results indicate that UV sterilisation of the residue reduces bacterial richness/evenness relative to both its crude counterpart and the water control. No significant interactions between PGPB and matrix were detected (p > 0.05), so alpha diversity remains largely stable across the biological treatments.

PCoA, Bray–Curtis, accounted for 32.9 and 8.1% of the variance on the first two axes, respectively (Figure 5A). A test for homogeneity of multivariate dispersions detected no significant differences among treatments (p > 0.05), indicating that the observed separation reflects compositional shifts rather than changes in within-group variance. Centroids for the uninoculated controls (C0) clustered near the origin, whereas addition of B. pretiosus (C1) displaced centroids along axis 1, consistent with a moderate enrichment of rare taxa; the P. agronomica treatment (C2) produced a milder shift, remaining largely within the C0 cloud.

Ward’s hierarchical clustering based on Euclidean distance (Figure 5B) corroborated the PGPB-driven segregation. Uninoculated plots (WC0, CFC0, OPKC0, OPKSTC0) formed a single clade, while plots receiving OPK® plus PGPB (OPK-C1/C2 and OPKST-C1/C2) grouped on distinct branches, underscoring the combined influence of the organic matrix and the inoculants. A reference line corresponding to the mean distance from the water control (WC0) shows that all inoculated treatments exceed this threshold, whereas non-inoculated chemical and organic controls fall below it.

Although the global PERMANOVA was non-significant (p > 0.05), the centroid positions and dendrogram topology indicate that B. pretiosus—and to a lesser extent P. agronomica—induces a detectable re-ordering of community composition, particularly when applied to the OPK® matrix.

PERMANOVA comparisons between each treatment and the water control (WC0) found no compositional effect attributable to the chemical matrix alone (p > 0.05 for CF0, OPK0 and OPKSTC0). By contrast, inoculation with B. pretiosus (C1) shifted the centroid of every community—whether applied with mineral fertiliser, crude residue or sterilised residue—and generated Bray–Curtis distances significantly greater than those to WC0 (Figure 6). P. agronomica (C2) produced an intermediate response, whereas centroid distances in uninoculated plots remained close to the control. PERMDISP (F = ∞, p = 0.532) confirmed homogeneity of within-group variance, indicating that the PERMANOVA signal reflects genuine compositional change rather than increased dispersion among replicates.

Absolute read counts (Figure 7) confirm successful establishment of the inoculated strains. In every C1 treatment the proportion of the genus Bacillus (dark-green block) increased markedly, whereas in C2 treatments a parallel rise was observed for Pseudomonas (yellow block). This pattern was consistent across all chemical matrices—water, mineral fertiliser, crude residue and sterilised residue—demonstrating that both PGPB strains persisted in the rhizosphere irrespective of background nutrition. Abundances of the other dominant genera (e.g., Arthrobacter, Rubrobacter, Streptomyces) remained comparatively stable, indicating that the inoculants did not displace the resident microbiota but rather added an extra functional layer to the community.

Figure 8 displays three parallel box-plots in which each point is the pairwise distance of a sample to every water-treated community (WC), calculated with the Jaccard index and then decomposed into its turnover (taxon replacement) and nestedness (species-loss/gain) components. Distances are grouped by fertiliser: Water, mineral fertiliser (Chemical Fertilizer) and the valorised organic residue (ORGAON PK). Total Jaccard (red boxes, left panel) show median dissimilarities cluster around 0.55 for all three fertilisers and do not differ statistically (Kruskal–Wallis χ² = 3.50, p = 0.17), indicating a comparable overall share of genera between each treatment and the water reference. Nestedness (green boxes, centre panel) contributions are an order of magnitude smaller (medians < 0.15). Nonetheless, a clear upward gradient is evident: Water < Chemical < ORGAON PK, and the difference is highly significant (χ² = 24.0, p = 2.1 × 10⁻⁶). This pattern denotes a modest but measurable introduction of genera unique to the organic residue.

Finally, turnover (blue boxes, right panel) dominates the dissimilarity structure. Median values rise from ≈ 0.30 (Water) to 0.55 (ORGAON PK), with the mineral fertiliser intermediate; the Kruskal–Wallis test confirms heterogeneity among groups (χ² = 9.0, p = 0.011). Hence, most compositional change stems from replacement of resident genera rather than simple gains or losses. Collectively, these results show that switching from water to either fertiliser reshuffles the bacterial community, but the organic residue drives a markedly higher taxon replacement while adding only a small nestedness component. The mineral fertiliser elicits a weaker, intermediate response, reinforcing the view that ORGAON PK induces the most pronounced restructuring of the soil microbiome.

Finally, to pinpoint the taxa underlying the fertiliser-driven β-diversity shifts (Figure 8), we applied ANCOM-BC—a compositionality-aware test that returns bias-corrected log-fold changes (LFC) and controls the false-discovery rate (FDR). Only features with FDR < 0.05 are displayed in the Supplementary Figure S3. At the family level (Supplementary Figure S3a), the sterilised residue (OPK_ST) was significantly enriched in five families—an uncultured lineage (t__uncultured), Prolixibacteraceae, Labraceae, Frankiaceae and Rikenellaceae—all showing LFC ≈ +0.8 to +1.0 relative to the water control. The crude residue (OPK) exhibited a single significant depletion of Xanthobacteraceae (LFC ≈ −0.4; Supplementary Figure S3b). At the genus level (Supplementary Figure S3c), OPK_ST displayed a broader enrichment, with 10 genera rising by ≥0.6 log₂ units: Steroidobacteraceae (unclassified genus), Chitinophaga, Jatrophihabitans, JG1-0001001-H03, Corallococcus, Nitrosomonadaceae (unclassified genus), Variovorax, Sandaracinus, WCHB1-32 and Labrys. No significant changes were detected for the mineral fertiliser (CF) or for either inoculated matrix at the family or genus rank, indicating that OPK_ST alone drove the strongest compositional signal. These results demonstrate that the sterilised organic residue selectively recruits polymer-degraders (e.g., Chitinophaga), nitrifiers (Variovorax, Nitrosomonadaceae) and other specialised taxa without broad gains or losses in richness, complementing the turnover patterns shown in Figure 8.