In October 2019, approximately 200 kg of soil (0 to 20 cm deep), was collected from Jianning County, Sanming city, China (34°04’ N, 108° 10′ E) with no history of cucumber cultivation. The red soil type is classified as Ferrosols-Udic in the Chinese Soil Taxonomy (Huang et al., 2017). The soil samples were air-dried in a natural manner, sieved through a 5 mm mesh to eliminate plant debris and rocks, homogenized, and subsequently stored in plastic bags at room temperature prior to utilization. From 20th March 2020 to 20th May 2020, we carried out a pot experiment (Figure 1A) at the Nanjing Agricultural University Test Base (31°36’ N, 119°10′ E), Baima Town, Nanjing, Jiangsu Province. The climate type at test site is the transition zone from north subtropical to mid-subtropical. The climate is mild and humid, the annual average for temperature is 15.4°C, rainfall is 1087.4 mm, and the frost-free period is 237 d. Before the experiment, the soil properties were: the pH value is 4.65, soil organic matter content (SOM) is 18.2 g kg⁻¹, total nitrogen content (TN) is 0.96 g kg⁻¹, available potassium content (AK) is 122.7 mg kg⁻¹ and available phosphorus content (AP) is 74.3 mg kg⁻¹.

The seeds of cucumber cultivar (Cucumis sativus L.), “Feng-Lu,” a commercial variety in China, were purchased from Jiangsu Academy of Agricultural Sciences. On March 20, 2020, the cucumber seeds were surface sterilized with 75% ethanol for 30 s followed by 5% NaClO for 5 min before planting. The sterilized seeds were placed in Petri dishes with wet autoclaved filter paper in a growth chamber (25°C, 70% relative humidity in the dark). After 2 days of pregermination, the cucumber seedlings (two true-leaf stage) were transplanted to plastic pots (15 cm diameter, 20 cm height) with 3 kg of soil.

The three treatments were as follows: (1) CK, no fertilizer control; (2) CF, chemical fertilizer treatment; and (3) BIO, bioorganic fertilizer treatment. CF treatment was consisting of analytically pure urea (N 46%), calcium superphosphate (P₂O₅ 15%), and potassium sulphate (K₂O 50%). Bioorganic fertilizer (pH 6.6, N 4.77%, P₂O₅ 2.26%, K₂O 1.00%, organic matter content 40%) was purchased from Jiangsu Lianye Fertilizer Co., Ltd., (the colony count of Bacillus amyloliquefaciens SQR9 was 1 × 10⁸ CFU g⁻¹ dry weight). Based on the equal nutrient content and nutrient requirements of cucumber (N: P₂O_5: K₂O = 1:0.5:1), CF and BIO treatments were adjusted to the same amounts of N, P and K according to the nutrients of BIO treatment if necessary by urea, calcium superphosphate, and potassium sulphate, respectively. The amount of nitrogen, phosphorus and potassium in each fertilized pot was 2.86 g N, 1.36 g P₂O₅ and 2.86 g K₂O, respectively. Thus, each BIO chamber received 60 g (2% w/w) bioorganic fertilizer and 2.26 g potassium sulphate. Each CF chamber was supplied 6.22 g urea, 9.07 g calcium superphosphate, and 5.72 g potassium sulphate. Bioorganic fertilizer and calcium superphosphate were used as basal fertilizers on November 25, 2019. The urea and potassium sulphate was applied on March 28, 2020. There were 36 seedlings of each treatment, for a total of 108 pots (Figure 1A); the pots were then randomly placed in a growth chamber (28/26°C day/night cycle, 70% relative humidity, and 180 μmol light m⁻² s⁻¹) and irrigated as needed. On 20th May 2020, 6 pots were randomly pooled into one sample, and thus 6 replicates from each treatment. In total, 18 samples were collected. The comprehensive workflow of our research is depicted in Figure 1A.

On 20th May 2020, plant growth traits including height (cm), stem thickness (D, cm), relative growth rate (RGR), leaf area (LA), and root traits [e.g., average lateral root diameter (LRAD), total lateral root number (TLRN), total lateral root length (TLRL), mean inter-branch density (MID), maximum order of lateral roots (MaxO), and branching density (BI, ration of branches number to root length)] were measured following previous protocols (Kang et al., 2021b). The methods used for collecting the plant samples are detailed in the Supplementary Method S1.

Rhizosphere soil was collected and subjected to the collection of root samples on May 20, 2020. In summary, the loosely attached soil on the cucumber roots was removed and discarded, while the soil that remained tightly adhered was collected as rhizosphere soil (Wen et al., 2020). Concurrently, bulk soil samples were collected. A portion of approximately 180.0 g of soil was air-dried for the determination of soil physicochemical properties, whereas another portion of approximately 10.0 g was stored at 4°C for the analysis of soil extracellular enzyme activities. The remained part of about 10.0 g was stored at −80°C and used for the microbial analyses. Detailed procedures for determining those properties are described in Supplementary Method S2.

About 0.5 g rhizosphere soil was used to extract the Genomic DNA by using the PowerLyzer PowerSoil DNA Isolation Kit (Qiagen, Germany) in accordance with the manufacturer’s instructions. Bacterial V4 region for amplicon sequencing, and the primers 515F: GTGYCAGCMGCCGCGGTAA and 806R: (GGACTACNVGGGTWTCTAAT) (Walters et al., 2015) to yield an amplicon of 292 bp. Detailed procedures for PCR amplification are described in Supplementary Method S3.

In this study, the “rhizosphere metabolome” refers to the collective pool of small molecule metabolites in the rhizosphere, including those released by cucumber roots and produced by the microbial community. Soil samples were divided into two portions (0.2 g each) and transferred into 2 mL EP tubes. Subsequently, 24 μL of Adonitol solution (1 mg mL⁻¹ in dH₂O) was added as an internal standard for rhizosphere metabolome extraction. The detailed extraction procedures are provided in Supplementary Method S4. The GC-TOF-MS (gas chromatography-time of flight mass spectrometry) analysis and the analysis of the raw peak performed as reported by Wen et al. (2020).

In a soil application trial, seven low molecular weight metabolites (ST; glycerol, Sorbitol, Phytol, 1,2,4-Benzenetriol, succinate semialdehyde, alpha-ketoglutaric acid and D-Glyceric acid) were tested for their effects on the assembly of soil microbial communities. Detailed procedures for incubating are described in Supplementary Method S5. Finally, 12 samples for the two treatments (2 treatments × 6 soil samples) were obtained and stored at −80°C.

The primers 341F/806R (F: CCTAYGGGRBGCASCAG; R: GGACTACHVGGGTWTCTAAT) were used to amplify the V3–V4 region of the 16S rRNA gene of the bacterial communities from the 12 collected samples with an amplicon size of 465 bp. The target sequences were identified from the raw sequences by aligning them with the 515F/806R primers, ensuring that they corresponded to the same region of the 16S rRNA gene as observed in the rhizosphere samples (Wen et al., 2020). Afterward, the PCR amplification and sequencing of the 16S rRNA were conducted in a manner consistent with the procedures outlined previously.

The alpha diversity of microbial communities was assessed through statistical analysis and visualization, utilizing metrics such as Richness, Shannon, and Pielou evenness. This analysis was conducted on the Amplicon Sequence Variant (ASV) table with the minimum read number across all samples, and the results were visually represented using boxplots. Next, the relative abundance of each ASV was standardized using the normalize_table.py script in Qiime (version 1.9.1) (Caporaso et al., 2010), followed by the preparation of Bray–Curtis similarity matrices using the beta_diversity.py script. Principal coordinate analysis (PCoA) plots were produced by utilizing Bray–Curtis similarity matrices generated with the R package ggplot2 (Wen et al., 2020). For beta diversity, PERMANOVA (Adonis, transformed data by Bray-Curtis, permutation = 999) (Gower, 1966) was used to determine the similarity of bacterial communities among treatments. Network analysis utilizing the Spearman correction and its associated properties were computed using the “igraph” package (Wen et al., 2022). All analyses and graphical visualizations were performed in R Environment using libraries dplyr, reshape2 and ggplot2. A nonparametric t-test was used to determine significant differences in alpha diversity, metabolites and soil properties among treatments using the “EasyStat” package in R with a false discovery rate (FDR).

Here, we used two approaches to examine bacterial community assembly. First, we applied the neutral model described by Merifield and Sloan (2006). Second, we estimated the influence of deterministic processes on community assembly by calculating the β-nearest taxon index (βNTI) between pairs of samples as described in Stegen et al. (2012). Detailed procedures for analyzing the bacterial community assembly are described in Supplementary Method S6.

In order to investigate the impact of dispersal on community assembly processes, we assessed ASV turnover using the abundance-weighted Raup-Crick metric (RCbray), following the methodology outlined by Stegen et al. (2013). RCbray assesses whether ASV turnover among sites differs from what would be expected based on ecological drift alone. Initially, the Bray–Curtis dissimilarity was calculated for each sample pair. Subsequently, a null-model community with equivalent size and richness was randomly generated for each sample. Every null model is created by randomly selecting ASVs (weighted by frequency across all samples) and assigning their relative abundance based on their relative abundance in the metacommunity. The Bray–Curtis dissimilarity between all pairs of simulated communities is subsequently calculated. This procedure is iterated 999 times to produce a null model distribution. The Relative Contribution of Bray–Curtis dissimilarity (RCbray) was determined through a calculation involving the summation of the number of simulated communities exhibiting a Bray–Curtis dissimilarity greater than the observed dissimilarity (Nsim>obs), half of the number of simulated communities with a Bray–Curtis dissimilarity equal to the observed dissimilarity (Nsim = obs), and subsequent division by the total number of simulations (999). Sample pairs with a |βNTI| < 2 and an |RCbray| > 0.95 suggest that bacterial community turnover is primarily driven by stochastic processes. Conversely, the bacterial community turnover is primarily driven by deterministic processes (Stegen et al., 2012). To determine the relative proportions of community assembly governed by deterministic and stochastic processes within each sample, we combined the results of the βNTI and RCbray analyses.

Additionally, the raw sequence data reported in this paper have been deposited in the Genome Sequence Archive of the BIG Data Center, Chinese Academy of Sciences, under accession code PRJCA034325.

A total of 969,634 high-quality sequences were obtained across 36 soil samples, including both bulk soil and rhizosphere samples. The average read count for each sample was 35,824 (standard deviation (SD) 2,388). All sequences were clustered into 6,284 ASVs. The majority of ASV belonged to the phyla Proteobacteria (57.99% with sd 0.065), Actinobacteria (30% with sd 0.051), Bacteroidetes (5.13% with sd 0.02), Firmicutes (2.01% with sd 0.014), Unassigned (1.71% with sd 0.005), Chloroflexi (1.68% with sd 0.008), Acidobacteria (0.51% with sd 0.003) Verrucomicrobia (0.32% with sd 0.002), Planctomycetes (0.28% with sd 0.002) and Spirochaetes (0.13% with sd 0.001).

The following section will examine two approaches to the study of bacterial community assembly. The initial step was to assess the extent to which stochastic processes are involved in the assembly of soil and rhizosphere bacterial communities. The neutral model, which accounts for 80.9% of the variation in ASV detection frequency, was estimated to have an estimated migration factor (m) of 0.26 in bulk soil (Figure 1B). In rhizosphere soil, the neutral model accounted for 50.5% of the variation in ASV detection frequency, with an estimated migration factor (m) of 0.54 (Figure 1C). The first step in the study was to assess the phylogenetic signal in habitat preference based on root traits, plant growth traits, rhizosphere metabolites, and soil chemical properties. A notable positive correlation was found between phylogenetic distance and habitat preference in closely related taxa (Supplementary Tables S2, S3), thus validating the assumptions of the βNTI metric. What’s more, deterministic processes (βNTI >2) was observed in the rhizosphere, while stochastic process (|βNTI| < 2) dominate phylogenetic turnover in bulk soil (Figures 1D,E).

Examination of the rhizosphere metabolome using GC-TOF-MS resulted in the detection of 323 chromatographic peaks, with 172 compounds identified in all samples. These included 14 amino acids and amides, 7 amines, 2 Pyridines, 26 alcohols, 61 acids, 4 phenols, 1 Alloxazine, 1 Purine, 4 aldehyde, 1 Alkaloids, 27 sugars, 34 sugar acids, 11 ketones, 1 Adenosines, 8 esters and 3 others (Supplementary Figure S1). Moreover, RDA/CCA showed that soil chemical properties (ANOVA, F = 1.53, p = 0.04) and important rhizosphere metabolites (ANOVA, F = 1.04, p = 0.001) were significantly affected bacterial composition at the ASVs level, while enzyme activities, root traits and indicators of aboveground growth were not significantly (Figures 2A–E). Mantel tests showed root exudates were a significant correlation (R = 0.397, p = 0.001) with rhizosphere microbiome. In contrast, no correlation was found between chemical properties, enzyme activities and root traits with rhizosphere microbiome (Supplementary Figure S2). Based on the random forest model, intriguingly, the rhizosphere microbiome was significantly correlated with glycerol (RET40), 7,8-Dimethylalloxazine (RET 319) and Capric Acid (RET80) (Figure 2F).

In order to ascertain the extent to which these properties account for community assembly patterns, we conducted analyses to assess the correlations between βNTI and various chemical properties, enzyme activities, root traits, and rhizosphere metabolites. There was no significant correlation between each index of chemical properties, enzyme activities and root traits with βNTI (Supplementary Figure S3). Subsequently, a correlation analysis was conducted between 11 groups of rhizosphere metabolites and βNTI, which identified a significant correlation between Phytol and βNTI (R = 0.25, p = 0.009) (Supplementary Figure S4). Furthermore, 26 compounds belonging to Phytol were analyzed for correlation with βNTI, revealing that only three compounds—glycerol (RET40), sorbitol (RET46), and Phytol (RET44)—showed significant correlations (Supplementary Figure S5). In order to ensure comprehensive coverage of potentially influential compounds in the microbiome assembly process, compounds exhibiting a coefficient of variation (CV) greater than 0.5 were selected for correlation analysis with βNTI. The results indicated that seven compounds (1,2,4-Benzenetriol, succinate semialdehyde, alpha-ketoglutaric acid, D-Glyceric acid, Glucoheptonic acid, 2-amino-2-methylpropane-1,3-diol, and 2-Butyne-1,4-diol) displayed significant correlations with βNTI (Figure 3A; Supplementary Figure S6).

In the investigation of core taxa during microbial community assembly, genera with a relative abundance greater than 0.01 were selected for correlation analysis with βNTI. The results indicated a significant correlation between the genus Rhizobium and βNTI (Supplementary Figure S7), prompting the removal of false positive correlations between taxa with high non-zero values and βNTI. Subsequent recalculation of βNTI and RCbray after excluding the genus Rhizobium revealed that stochastic processes (|βNTI| < 2) predominantly govern phylogenetic turnover (Figure 3H). Finally, in order to enhance the validity of the findings, βNTI and RCbray were recalculated iteratively by systematically excluding each Genus (except Rhizobium) individually, with a focus on instances where the mean value of βNTI >2 (Figures 3I,J).

Seven small molecule metabolites, including glycerol, sorbitol, phytol, 1,2,4-benzenetriol, succinate semialdehyde, alpha-ketoglutaric acid, and D-glyceric acid, were ultimately chosen for the soil culture experiment based on a selection criterion of R > 0.2. We successfully acquired 1,032,205 high-quality sequences from a combined 12 soil samples, encompassing both control (CK) and treatment (ST) groups. The average read count for each sample was determined to be 86,017 (SD 20476). All sequences were clustered into 15,126 ASVs. The majority of ASV belonged to the phyla Actinobacteria (38.12% with sd 0.05) Proteobacteria (32.33% with sd 0.06), Acidobacteria (12.92% with sd 0.021), Firmicutes (5.49% with sd 0.024) Unassigned (3.62% with sd 0.038) Gemmatimonadetes (1.56% with sd 0.005), Chloroflexi (0.75% with sd 0.002), Nitrospirae (0.47% with sd 0.001), Planctomycetes (0.36% with sd 0.003) and Bacteroidetes (0.3% with sd 0.001). However, there were no significant differences in α-diversity, as indicated by Chao1, Shannon, and Pielou evenness values, between the CK and ST (Figure 4A). A higher prevalence of Planctomycetes was observed in the ST treatment, whereas Actinobacteria were more abundant in the CK treatment (Figure 4B). The Principle coordination analysis (PCoA) illustrated the bacterial community structures dissimilarities (Adonis, p = 0.001, R = 0.463 PERMANOVA) between CK and ST treatments (Figure 4C). Moreover, network analysis identified more connections in the CK than in the ST (Figure 4D; Supplementary Table S5).

By comparing the βNTI and RCbray metrics, respectively, we further investigated whether and how community assembly processes differ between CK and ST. We found that variable selection (βNTI >2) dominates phylogenetic turnover in ST, whereas stochastic processes (|βNTI| < 2) dominate phylogenetic turnover in CK (Figure 4E). The Mantel tests revealed a significant correlation (R = 0.682, p = 0.001) between the ST and the composition of the rhizosphere microbiome, whereas no significant correlation (R = 0.375, p = 0.194) was found between the CK and the rhizosphere microbiome (Figure 4F). Furthermore, the presence of Rhizobium, a keystone taxa in the assembly process of rhizosphere bacterial communities, was more abundant in samples from the ST compared to those from the CK (Figure 4G).