Dissecting the role of soybean rhizosphere-enriched bacterial taxa in modulating nitrogen-cycling functions

Abstract Crop roots selectively recruit certain microbial taxa that are essential for supporting their growth. Within the recruited microbes, some taxa are consistently enriched in the rhizosphere across various locations and crop genotypes, while others are unique to specific planting sites or genotypes. Whether these differentially enriched taxa are different in community composition and how they interact with nutrient cycling need further investigation. Here, we sampled bulk soil and the rhizosphere soil of five soybean varieties grown in Shijiazhuang and Xuzhou, categorized the rhizosphere-enriched microbes into shared, site-specific, and variety-specific taxa, and analyzed their correlation with the diazotrophic communities and microbial genes involved in nitrogen (N) cycling. The shared taxa were dominated by Actinobacteria and Thaumarchaeota, the site-specific taxa were dominated by Actinobacteria in Shijiazhuang and by Nitrospirae in Xuzhou, while the variety-specific taxa were more evenly distributed in several phyla and contained many rare operational taxonomic units (OTUs). The rhizosphere-enriched taxa correlated with most diazotroph orders negatively but with eight orders including Rhizobiales positively. Each group within the shared, site-specific, and variety-specific taxa negatively correlated with bacterial amoA and narG in Shijiazhuang and positively correlated with archaeal amoA in Xuzhou. These results revealed that the shared, site-specific, and variety-specific taxa are distinct in community compositions but similar in associations with rhizosphere N-cycling functions. They exhibited potential in regulating the soybean roots’ selection for high-efficiency diazotrophs and the ammonia-oxidizing and denitrification processes. This study provides new insights into soybean rhizosphere-enriched microbes and their association with N cycling. Key points • Soybean rhizosphere affected diazotroph community and enriched nifH, amoA, and nosZ. • Shared and site- and variety-specific taxa were dominated by different phyla. • Rhizosphere-enriched taxa were similarly associated with N-cycle functions. Supplementary Information The online version contains supplementary material available at 10.1007/s00253-024-13184-5.


Introduction
The rhizosphere is the soil surrounding plant roots, which is a microbial hotspot and plays a pivotal role in soil nutrient cycling (Reinhold-Hurek et al. 2015).The rhizosphere microbiota contains a subset of the microbes in bulk soil that are recruited by root exudates and has been well characterized in several crops, including maize, rice, barley, and soybean (Bulgarelli et al. 2015;Edwards et al. 2015;Mendes et al. 2014;Walters et al. 2018).Many beneficial microbes involved in soil nitrogen (N) cycling have been identified in the microbiota enriched in the rhizosphere, including diazotrophs (e.g., Rhizobium, Pseudomonas), nitrifiers, and denitrifiers (Trivedi et al. 2020).In addition to the microbes that possess beneficial functions, some rhizosphere-enriched taxa also affect nutrient cycling functions through complex interactions within the community (Chepsergon and Moleleki 2023).Understanding the link between rhizosphere-enriched taxa and nutrient cycling function could provide insights into increasing soil health while meeting crop nutrient demands.
Crop roots recruit varied microbial communities under the effects of soil properties and crop genotype, which can be attributed to the differences in native microbial communities in the bulk soil pool and the composition of root exudations, respectively (Kelly et al. 2022;Zhalnina et al. 2018;Zhang et al. 2018).Although the rhizosphere microbial communities may differ in their compositions and functions under these effects, some shared taxa are found across different sites and crop genotypes, which are often called the "core microbiota" (Jiao et al. 2022b;Yeoh et al. 2017).These shared taxa are considered to possess essential functions such as N fixation for crop growth or health maintenance and are less influenced by environmental variability (Castellano-Hinojosa and Strauss 2021;Chang et al. 2022;Xu et al. 2020;Ling et al. 2022).As the rhizosphere effect is under the influence of local soil conditions (Berg and Smalla 2009;Xu et al. 2009), some rhizosphereenriched taxa are "site-specific," which are shared across different crop genotypes only in specific locations (Jin et al. 2017).These "site-specific" taxa can be more sensitive to soil and environmental properties and contribute more to crop adaptation to the local environment when compared to the "core microbiota" (Thiergart et al. 2020).Despite the similarly enriched microbial taxa among sites or genotypes, the unique microbial taxa enriched in different varieties are likely linked to the heritable physiological properties of the varieties, such as fungal resistance and N use (Mendes et al. 2018;Zhang et al. 2019).Therefore, within the rhizosphere-enriched microbiota, microbial recruitment may be controlled through different mechanisms and contribute to community function differentially.
However, the differences in the taxonomy and function of the rhizosphere-enriched taxa shared across sites and genotypes and those specific to sites or varieties have not been well studied.
Soybean (Glycine max) is a major legume crop that can symbiose with rhizobia and fix N from the air (Schultze and Kondorosi 1998).Soybean relies less on synthetic N fertilizer when compared to nonlegume crops and inputs a considerable amount of natural N into the agricultural system, which affects the nutrient cycling function, especially the N-cycling function, of the rhizosphere microbiota (Sánchez and Minamisawa 2019;Tsiknia et al. 2020).In return, the abundance of the taxa involved in the N cycle in the rhizosphere community can also affect the efficiency of N transformation processes in the rhizosphere and thereby influence crop N use efficiency (Ren et al. 2023).Moreover, it was reported that the rhizosphere microbiota showed a close association with rhizobia, for example, the genus Bacillus in the soybean rhizosphere was associated with the symbiotic efficiency of rhizobia (Han et al. 2020).Therefore, the N-cycling bacteria of the soybean rhizosphere may be closely associated with rhizosphere-enriched taxa, but information about this association is still limited.
Accordingly, five soybean varieties with different variety types, growth habits, and stem terminations were planted in two fields with different soil types and environmental conditions, and the bacterial and diazotrophic communities and genes involved in the N cycle in rhizosphere and bulk soils were studied.The objectives of this study were (1) to determine the rhizosphere effect on the bacterial community and the microbes involved in the N cycle in different soybean varieties and field conditions; (2) to characterize the rhizosphere-enriched taxa that are shared among sites and varieties (shared taxa), shared among varieties but not sites (site-specific taxa), and unique to each variety (varietyspecific taxa) in the soybean rhizosphere; and (3) to explore the potential relationship between the N-cycling function and the shared and site-and variety-specific taxa enriched in the soybean rhizosphere.The results of this study will provide insights into the regulatory mechanisms of the nutrient cycling function in the soybean rhizosphere.
The study region has a warm temperate continental monsoon climate.The annual mean temperature and rainfall were 11.6 and 14.0 °C and 500 and 860 mm in Shijiazhuang and Xuzhou, respectively.The soil contains 72.5% sand (0.02-2 mm), 15.5% silt (0.002-0.02 mm), and 12.0% clay (< 0.002 mm) in Shijiazhuang and 49.7% sand, 30.6% silt, and 19.7% clay in Xuzhou.For both sites, wheat (Triticum aestivum L.) (sown in October and harvested in June) was rotated with soybean (sown in June and harvested in October) in a double-cropping system for more than 5 years before the experiment.
The field experiment was laid out following a randomized complete block design with five soybean varieties as treatments and four replicated for each treatment at both experimental sites.Each plot was 1.6 m wide and 5.5 m long.The five soybean varieties selected had different types, growth habits, and stem terminations according to the SoyFGB v2.0 platform (Supplemental Table S1) (Zheng et al. 2022).Soybean seeds were sown manually by 0.4 * 0.11 m 2 spacing on 16 and 12 June in Shijiazhuang and Xuzhou, respectively.The fertilizer was applied after sowing at the rate of 22.5 kg N ha −1 (applied as urea, 46% N), 34.5 kg P 2 O 5 ha −1 (applied as calcium superphosphate, 12% P 2 O 5 ), and 34.5 kg K 2 O ha −1 (applied as potassium sulfate, 50% K 2 O) for each site and each plot.Pesticides were applied the same as local cropping management.No irrigation was applied and weeds were controlled manually.The soybeans were not inoculated with rhizobia prior or during planting.
Five bulk soil samples were collected between soybean plant rows at approximately 20 cm away from soybean plants and 0-20 cm depth following the five-point sampling method at each site.One plant from each plot (four plants for each variety at each site) except the side row was randomly selected and excavated with roots and the surrounding soil at each site at the full bloom stage.The plants were shaken vigorously to remove the loosely adhered soil on the roots, and the roots together with rhizosphere soil were collected into an ethylene-oxide-sterilized ziplock bag (BKMAM, Changde, China).All samples were immediately transferred to the laboratory with ice packs.The rhizosphere soil was manually separated from the roots using a small brush according to Zhang et al. (2017).After removing fine residues and stones, each soil sample was divided into three parts.One part was stored at − 20 °C for DNA extraction; one was stored at 4 °C for the determination of pH and the contents of water, ammonia, and nitrate; and the remaining soil was air-dried, ground, and sieved to determine the contents of soil organic carbon, total nitrogen, and available phosphorous and potassium.

Measurement of soil properties
Soil pH was measured using a pH meter (PB-10, Sartorius, Göttingen, Germany) at a soil to distilled water ratio of 1:5 (weight/volume).The soil water content was measured by oven drying at 105 °C until a constant weight was achieved.Soil organic carbon and total N contents were measured using an elemental analyzer (Vario MACRO Cube, Elementar, Hanau, Germany).Soil nitrate and ammonia were extracted using 1 M KCl, and their contents were measured using a continuous flow autoanalyzer (Seal Auto Analyser 3, Seal Analytical, Norderstedt, Germany).The available phosphorous was extracted using 1 M NaHCO 3 , and its content was measured following the molybdenum-blue colorimetry method.The available potassium was extracted using 1 M NH 4 Ac, and its content was measured using a flame photometer (FP640,INASA,Shanghai,China).

Soil DNA extraction and amplicon sequencing
Total soil DNA was extracted from 0.5 g of each bulk and rhizosphere soil sample using the Fast DNA SPIN kit (MP Biomedicals, Irvine, CA, USA) following the manufacturer's instructions.The integrity and concentration of the extracted DNA samples were evaluated using 1% agarose gels.The V4 hypervariable region of the 16S rRNA gene and the nifH gene fragment were amplified using general universal primer sets 515F/806R (5′-GTG YCA GCMGCC GCG GTAA-3′/5′-GGA CTA CNVGGG TWT CTAAT-3′) and polF/polR (5′-TGC GAY CCSAARGCBGACTC-3′/5′-ATSGCC ATC ATY TCR CCGGA-3′), respectively, with Phusion high-fidelity DNA polymerase (New England Biolabs, Ipswich, MA, USA).The polymerase chain reaction (PCR) conditions consisted of an initial denaturation at 98 °C for 1 min, 30 cycles of denaturation at 98 °C for 10 s, annealing at 50 °C for 30 s, elongation at 72 °C for 30 s, and a final extension at 72 °C for 5 min.The PCR products were assessed on a 2% agarose gel, pooled in equidensity ratios, and purified with a GeneJET Gel Extraction kit (Thermo Scientific, Waltham.MA, USA).The sequencing libraries were constructed using a TruSeq DNA PCR-Free Library Preparation kit (Illumina, San Diego, CA, USA).The quality of the libraries was assessed on the Qubit 2.0 Fluorometer (Thermo Scientific, Waltham.MA, USA) and Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).The libraries were sequenced on an Illumina NovaSeq platform, and 250-bp paired-end reads were generated.The quality of the paired-end reads was assessed with FastQC v.0.11.5 (Andrews 2014).The raw sequences were deposited in the NCBI repository under BioProject accession number PRJNA1045824.

Quantitative PCR (qPCR) of N-cycling genes
The abundances of seven genes involved in N cycling were quantified by qPCR in this study.The name and function of the genes and the primer sets used are shown in Supplemental Table S2.The fragments of each gene were amplified from a random total soil DNA sample, ligated to the pGEM-T vector (Promega, Madison, WI, USA), and then transformed into Escherichia coli strain DH5α (Tiangen, Beijing, China).The cloning vector of each gene was screened by culturing on LB medium with IPTG (isopropyl ß-D-1-thiogalactopyranoside) and X-gal according to the method described by Green and Sambrook (2019) and confirmed by PCR amplification with the corresponding primer sets and sequencing analysis.The concentration of the extracted cloning vector of each gene was quantified using a NanoDrop spectrophotometer (NanoDrop, PeqLab Biotechnologie, Erlangen, Germany) and then calculated using the following equation: where N is the number of base pairs contained in the cloning vector.The quantified cloning vector of each gene was diluted to a tenfold dilution series from 10 to 10 7 copies μL −1 , which were later used as standards.The total soil DNA samples were quantified and then diluted to a concentration of approximately 100 ng µL −1 before being used as qPCR templates.
qPCR was performed on an Applied Biosystems 7500 Real-Time System (Life Technologies, Carlsbad, CA, USA) in a 20-μL mixture containing 10 μL of TB Green Premix Ex Taq II (Tli RNaseH Plus) (TaKaRa, Otsu, Japan), 0.4 μM of each forward and reverse primer for each gene, 0.4 μL of ROX reference dye II (50 ×), 2 μL of DNA template, and 6 μL of ddH 2 O.The qPCR conditions consisted of an initial denaturation at 95 °C for 30 s, 40 cycles of 95 °C for 5 s and 60 °C for 34 s, and a dissociation stage.All of the total soil DNA samples and standards were run in triplicate for technical replicates.Standard curves were generated by regressing the standard's copy number with its cycle threshold (Ct) value.The copy number of each gene in each soil sample was determined by the mean Ct value of the three (1) Copies of the cloning vector per L = ng cloning vector per L × 6.02 × 10 14 ∕N × 660 technical replicates using the corresponding standard curves.The absolute abundance of a gene refers to the copy number per gram soil dry weight.

Statistical analysis
The statistical analyses were performed using R software (version 4.0.2,https:// www.r-proje ct.org/).Two-way ANOVA followed by Tukey's HSD test was performed to determine the effects of soybean rhizosphere and variety on soil properties, the alpha diversity indices, and the relative abundances of the top 10 most abundant phyla and top 20 most abundant orders of bacterial communities.For alpha diversity indices of bacterial communities, richness was calculated as the number of OTUs observed in each sample, the Shannon index was calculated using the diversity function in the vegan package (https:// cran.r-project.org/ web/ packa ges/ vegan/ index.html), and the phylogenetic diversity (PD) index of the whole tree was calculated using the pd function in the picante package (https:// cran.r-project.org/ web/ packa ges/ pican te/ index.html).For the beta diversity of bacterial and diazotrophic communities, Bray-Curtis dissimilarity among samples was calculated using the vegdist function in the vegan package (https:// cran.r-project.org/ web/ packa ges/ vegan/ index.html), and the effects of the experimental site, compartment, and soybean variety were visualized by principal coordinate analysis (PCoA) and constrained PCoA (CPCoA).A permutational multivariate ANOVA (PERMANOVA) with 999 random permutations was performed using the Adonis function in vegan to determine the significance of the effects.The differences between the relative abundance of each OTU in the bulk and rhizosphere soil bacterial communities were assessed using the DESeq2 package (https:// bioco nduct or.org/ packa ges/ relea se/ bioc/ html/ DESeq2.html) at FDR (false discovery rate)adjusted p < 0.05 and estimated log2-fold change > 1.The Venn diagrams of the enriched OTUs were generated online at https:// bioin forma tics.psb.ugent.be/ webto ols/ Venn/ to identify the shared, site-specific, and variety-specific taxa.The phylogenetic trees were generated using the online tool iTOL (https:// itol.embl.de/) (Letunic and Bork 2021).The co-occurrence networks between the diazotrophic orders and the OTUs in the bacterial community were constructed using the ggClusterNet package (Wen et al. 2022) with r and p thresholds of 0.6 and 0.05, respectively.Student's t test was performed to test the difference between the absolute abundance of N-cycling genes in bulk and rhizosphere soils at each site at p < 0.05.The random forest model was applied to assess the importance and significance of the effects of each OTU on the community N-cycling function using the rfPermute package (https:// cran.r-project.org/ web/ packa ges/ rfPer mute/ index.html).Spearman correlation and Mantel tests were performed to analyze the correlation of the shared, site-specific, and variety-specific taxa with soil properties and the abundance of N-cycle-related genes using the correlate and Mantel_test functions in the linkET package (https:// github.com/ Hy4m/ linkET), respectively.

Soil properties
Soil properties were different between Shijiazhuang and Xuzhou, bulk and rhizosphere soil, and among soybean varieties (Table 1).The soil water content, organic carbon, and total N were significantly (p < 0.05) lower, while ammonium and available potassium were significantly higher in both bulk and rhizosphere soil in Shijiazhuang than in Xuzhou.The soil pH and organic carbon, total N, nitrate, and available potassium contents were significantly higher, and the water content was significantly lower in the soybean rhizosphere than in the bulk soil at both sites.The soil ammonium content was significantly lower in the rhizosphere than in the bulk soil in Xuzhou, and the available phosphorus content in the rhizosphere was significantly higher in Shijiazhuang but significantly lower in Xuzhou than in the bulk soil.Although all soil properties except for soil pH were significantly different among the rhizosphere soils of the five soybean varieties in at least one site, none of them exhibited a consistent pattern of variation between the two sites.

Bacterial communities
The soybean rhizosphere effect strongly influenced the diversity and composition of soil bacterial communities in Shijiazhuang and Xuzhou (Fig. 1).The PD whole tree index of the rhizosphere soil for the five cultivars was significantly (p < 0.05) lower than that in the bulk soil at both sites, while the richness and Shannon indices in rhizosphere soil except for QX were similar to those in bulk soils (Fig. 1a).The richness in the rhizosphere of QX was significantly lower than that in bulk soils in Shijiazhuang.The PCoA and CPCoA based on Bray-Curtis dissimilarity showed significant effects of the experimental site, sample compartment, and soybean variety on the bacterial communities (Fig. 1b).When both bulk and rhizosphere soil samples were included, the sample compartment explained 44% of the variance (p < 0.001), followed by the experimental site (R 2 = 0.14, p < 0.001), and their interactive effect was very small (R 2 = 0.04, p < 0.01).CPCoA was further performed on the rhizosphere samples, which showed small but significant effects of soybean variety (R 2 = 0.11, p < 0.05) and its interaction with the experimental site (R 2 = 0.13, p < 0.01) on the rhizosphere bacterial communities.
The bacterial community composition was significantly affected by the experimental site, compartment, and soybean variety (Fig. 1c, Supplemental Table S3).At the phylum level, Proteobacteria (24.4%) and Firmicutes (28.7%) were the most dominant in bulk soil, and Proteobacteria (29.7%) and Actinobacteria (19.5%) were the most dominant in rhizosphere soil at both sites.Among the top 10 most abundant phyla, Proteobacteria, Actinobacteria, Acidobacteria, Thaumarchaeota, Chloroflexi, Gemmatimonadetes, and Planctomycetes were significantly more abundant in rhizosphere soil than in bulk soil, and Firmicutes, Bacteroidetes, and Verrucomicrobia were significantly more abundant in bulk soil than in rhizosphere soil.Bacteroidetes, Chloroflexi, and Planctomycetes were significantly affected by site, and Chloroflexi in the rhizosphere soil was significantly affected by variety.Among the top 20 most abundant orders, nine orders, including Rhizobiales, Bacillales, Acidimicrobiales, Rhodospirillales, Deltaproteobacteria GR-WP33-30, Propionibacteriales, Solirubrobacterales, Gaiellales, Myxococcales, and Micrococcales, were significantly enriched in the soybean rhizosphere for all the five soybean varieties at both sites (Supplemental Table S4).

Shared and site-and variety-specific taxa enriched in the soybean rhizosphere
The rhizosphere-enriched taxa were identified by comparing the bacterial community present in the rhizosphere of each soybean variety with that in the bulk soil at each site.Based on their intersections, these OTUs were then classified into three categories: shared, site-specific, and variety-specific taxa (Fig. 2).The number of OTUs enriched in the soybean rhizosphere ranged from 443 to 613, which was lower than the number of depleted OTUs (Fig. 2a).The 49 enriched OTUs that were shared across all sites and varieties were classified as shared taxa; the 126 and 70 enriched OTUs that were shared among soybean varieties but not sites were classified as site-specific taxa in Shijiazhuang and Xuzhou, respectively; and the rhizosphere-enriched OTUs that were unique to each soybean variety were classified as varietyspecific taxa (Fig. 2b).
The shared, site-specific, and variety-specific taxa showed differences in the composition of both OTU taxonomy and relative abundance, although they were similarly dominated by Actinobacteria and Proteobacteria (Fig. 3).Among the shared taxa, most of the OTUs belonged to Actinobacteria (63.3%) and Proteobacteria (16.3%), and Thaumarchaeota and Actinobacteria were the phyla with the highest relative abundances (Fig. 3a).In Shijiazhuang, Nitrospirae was the third most dominant phylum by OTU numbers in both site-specific and variety-specific taxa, but it was relatively low in abundance (Fig. 3b, c).Thaumarchaeta was the third most abundant phylum after Actinobacteria and Proteobacteria among site-specific taxa, while Proteobacteria and Firmicutes were the most abundant among variety-specific taxa.In Xuzhou, Nitrospirae was the most dominant phylum by both OTU number and relative abundance among site-specific taxa (Fig. 3d).Thaumarchaeota and Acidobacteria were the most abundant phyla with higher relative abundances than Actinobacteria and Proteobacteria among variety-specific taxa (Fig. 3e).Notably, the variety-specific taxa were widely distributed in 21 phyla in both Shijiazhuang and Xuzhou, and 391 out of 480 OTUs in Shijiazhuang and 358 out of 460 OTUs in Xuzhou were rare taxa with relative abundances less than 0.1‰ (Supplemental Table S5).The diazotroph community composition was different between rhizosphere and bulk soil and between Shijiazhuang and Xuzhou (Fig. 4a).At the order level, Rhizobiales was the most abundant at both sites.Myococcales was more abundant in Xuzhou, while Pseudomonadales and Micrococcales were more abundant in Shijiazhuang.At Shijiazhuang, Pseudomonadales was enriched in the rhizospheres of four out of five varieties (HJ, MC, CX, and QX), and Micrococcales was enriched only in the rhizosphere of HJ.At Xuzhou, the top five most abundant orders, including Rhizobiales, Myococcales, Pseudomonadales, Micrococcales, and Propionibacteriales, were more abundant in the soybean rhizosphere than in the bulk soil.Unlike the total bacterial community, the diazotrophic community was mainly affected by the experimental site (PERMANOVA: R 2 = 0.47, p < 0.001), and the effect of compartment (PERMANOVA: R 2 = 0.08, p < 0.001) was smaller than that of site (Fig. 4b).The diazotrophic community was separated by experimental site at PCo 1, while the communities in bulk soil and rhizosphere soil were separated at PCo 2 only in Xuzhou.
Co-occurrence interaction networks were constructed between the diazotroph orders and the OTUs in bacterial communities (Fig. 4c).In the both-site network, the shared and site-and variety-specific taxa contributed 294 out of the 888 total edges.The average degrees of the shared and site-and variety-specific taxa were 2.84, 2.90, and 3.40, respectively.More negative correlations (652 edges) than positive ones (236 edges) were found in the both-site network, but more positive than negative correlations were found for Propionibacterials, Burkholderiales, Chlorobiales, and Rhizobiales.In the Shijiazhuang network, a total of 502 edges (263 negative and 139 positive) were found, and shared, site-specific, and variety-specific taxa were linked with diazotrophic orders through 41, 70, and 19 edges, respectively.They were negatively correlated with Chlorobiales, Nostocales, Methylococcales, and Bacillales and positively correlated with Propionibacteriales and Pseudomonadales.In the Xuzhou network, a total of 969 edges (645 negative and 324 positive) were found, and shared, site-specific, and variety-specific taxa were linked with diazotrophic orders through 70, 91, and 79 edges, respectively.They were positively correlated with Propionibacteriales, Enterobacteriales, and Rhizobiales and negatively correlated with the remaining diazotrophic orders.The N-cycling function and its relationship with shared and site-and variety-specific taxa enriched in the soybean rhizosphere Seven genes involved in N cycling were quantified to determine their absolute abundance in bulk soil and soybean rhizosphere soil, and six of them were enriched in the soybean rhizosphere (Fig. 5a).The abundances of nifH (involved in N fixation), archaeal amoA (involved in ammonia oxidation), and nosZ (involved in N 2 O reduction) were significantly (p < 0.05) higher in the soybean rhizosphere than in the bulk soil at both sites.The abundances of bacterial amoA in Xuzhou and nitrite reductase genes nirS and nirK in Shijiazhuang were significantly higher in the soybean rhizosphere than in the bulk soil.The abundance of narG (involved in nitrate reduction) was significantly lower in the soybean rhizosphere than in the bulk soil in Shijiazhuang.A differentiation between bulk and soybean rhizosphere soil was also observed in the predicted N-cycling functions (Supplemental Fig. S1).The predicted N-cycling functions were generally enriched in the soybean rhizosphere at both sites, except for nitrate denitrification and nitrite denitrification and respiration in Shijiazhuang and nitrate reduction at both sites.
To determine the contribution of the shared, site-specific, and variety-specific taxa to the community N-cycling functions, the function of each OTU was predicted, and their correlation with the abundance of N-cycle-related genes was examined (Fig. 5b, c).Only 19 OTUs were annotated with N-cycle-related functions, and they belong to the phyla Proteobacteria, Actinobacteria, Thaumarchaeota, Nitrospirae, and Firmicutes (Fig. 5b).Among the 19 OTUs, 11 OTUs had read counts larger than 50 (relative abundance > 1‰), and they possessed the functions of ureolysis, nitrate reduction, ammonia oxidation, nitrite oxidation, nitrification, and N fixation.According to the random forest analyses and the cross-validation curves, the 13 and 12 most important OTUs at Shijiazhuang and Xuzhou were correlated with the abundance of each N-cycle-related gene, respectively (Fig. 5c).In Shijiazhuang, four shared, three site-specific, and six variety-specific OTUs were included.Most of the OTUs were negatively correlated with bacterial amoA and narG, and less than half of the OTUs were positively correlated with nifH, archaeal amoA, nirS, nirK, and nosZ.In Xuzhou, three shared, two site-specific, and seven variety-specific OTUs were included.More significant (p < 0.05) positive correlations were found than in Shijiazhuang, where more than half of the OTUs were positively correlated with nifH and bacterial amoA, and negative correlations were found only between nirS and two OTUs.
To further explore the interplay among soil properties, shared and site-and variety-specific taxa, and community Fig. 4 The diazotrophic community and its correlation with shared and site-and variety-specific taxa enriched in the soybean rhizosphere.a The relative abundance of the top 20 orders in the diazotroph community in bulk soil and the rhizosphere of five cultivars in Shijiazhuang and Xuzhou; b PCoA of the diazotroph community in Bray-Curtis distance; c the co-occurrence network between diazotroph orders and site-and variety-specific taxa enriched in the soybean rhizosphere.Orders with a degree higher than 10 were labeled Fig. 5 The abundance of N-cycle-related genes and the relationships between the rhizosphere-enriched taxa and N-cycle-related genes.a The absolute abundance of N-cycle-related genes in the bulk soil and soybean rhizosphere.Error bars represent standard error; b the read count of the rhizosphere-enriched OTUs with predicted N-cycle functions; c random forest models of the rhizosphere-enriched OTUs predicting N-cycle function and Spearman's correlation between the OTUs and the abundance of N-cycle-related genes.The OTUs with the top 13 and 12 increases in mean squared error (MSE) were visualized at Shijiazhuang and Xuzhou, respectively, as determined by the tenfold cross-validation error (Supplemental Fig. S2).*, **, and ** indicate significant differences at p < 0.05, 0.01, and 0.001, respectively Page 12 of 17 N-cycle functions, their relationships were examined (Fig. 6).When both sites were included, shared taxa were correlated with archaeal and bacterial amoA and narG, site-specific taxa were correlated with bacterial amoA and narG, and varietyspecific taxa were correlated with bacterial amoA.Shared and site-and variety-specific taxa were all correlated with pH and soil organic carbon, while only site-specific and variety-specific taxa were correlated with soil water content and total N.When the two sites were analyzed separately, the shared and site-and variety-specific taxa did not show much difference between each other except for the correlation with available P. They correlated with bacterial amoA and narG at Shijiazhuang and with archaeal amoA at Xuzhou, with pH and soil organic carbon at both sites, with available K at Shijiazhuang, and with soil total N and NH 4 at Xuzhou.Available P was correlated with variety-specific taxa at Shijiazhuang and shared and sitespecific taxa at Xuzhou.

Soybean rhizosphere effect on the composition and N-cycling function of the bacterial community
The rhizosphere effect was found to be a major determinant in modulating the bacterial community in this study.The soybean rhizosphere had a substantially reduced PD whole tree index, and 1097-1471 OTUs were significantly depleted while only 443-613 were enriched (Figs.1b and 2a).These results confirmed the strong selective role of the crop rhizosphere, which leads to a less diverse but more specialized community (Berendsen et al. 2012;Ling et al. 2022).Among the orders enriched in all soybean rhizosphere at both sites, Rhizobiales and Bacillales possess potential functions in N cycling, growth promotion, and disease suppression AOB bacterial amoA, SWC soil water content, SOC soil organic carbon, STN soil total nitrogen, AP available phosphorus, AK available potassium (Han et al. 2020;Xing et al. 2022).Furthermore, the small but significant effects of experimental sites and soybean varieties (Fig. 1b) together with the site-specific and variety-specific enriched orders (Supplemental Table S3) suggested a complex interaction between soybean root exudate profiles and the local microbiome (Zhalnina et al. 2018).
The differential composition of diazotrophs and abundances of N-cycle-related genes revealed that the soybean rhizosphere effect significantly modulates the functional potential of N cycling (Figs.4a and b and 5a and Supplemental Fig. S1).The significant enrichment of N-cycling genes such as nifH, archaeal amoA, and nosZ indicated higher potential activity of nitrogen fixation, nitrification, and denitrification, suggesting that soybean rhizosphere created a microenvironment that promotes higher potential activity for N-cycling functions.These results are consistent with the metagenomics results reported by Mendes et al. (2014), which revealed a high abundance of genes associated with N fixation, denitrification, and nitrate/nitrite ammonification in the soybean rhizosphere.However, they partially deviate from the generalized rhizosphere effect on N-cycling functions demonstrated by a synthesis analysis of published sequences worldwide, which includes the enrichment of N fixation and denitrification microbes and the strong depletion of nitrification microbes (Ling et al. 2022).This inconsistency may be attributable to the influence of local soil conditions such as soil N availability, water content, and texture (Séneca et al. 2020;Lian et al. 2023).The influence of these local soil conditions could impact the rhizosphere's ability to shape N-cycling microbial consortia (Ma et al. 2023), which is also revealed by the profound effects of the experimental site (explaining 47% of the total variation) on the diazotroph community (Fig. 4b) and the variable rhizosphere effects on the abundances of bacterial amoA, narG, nirS, and nirK (Fig. 5a).Furthermore, the rhizospheres of soybean varieties exhibited variations in the enrichment of the diazotroph orders Rhizobiales and Pseudomonadales (Fig. 4a) and the N-cycling-related functions nitrate/nitrite ammonification, N/nitrate respiration, and denitrification (Supplemental Fig. S1).These results suggested that the soybean genotype affected the recruitment of microbes with specific N-cycle-related functions in the soybean rhizosphere.This might be attributed to the differences in the growth and nutrient requirements of these soybean varieties, which may result in the differences in their root exudates.However, the regulatory mechanisms by which soybean genotypes control the recruitment of N-cyclerelated microbes in the rhizosphere remain to be elucidated, potentially through detailed profiling of the root exudates of different soybean varieties.
The taxonomy of the shared and siteand variety-specific taxa enriched in the soybean rhizosphere Among the soybean rhizosphere-enriched OTUs, the shared taxa across all samples and those that were specific to sites and soybean varieties were identified, and they exhibited distinct compositions and correlated differently to soil properties (Figs. 2, 3, and 6).
In the shared taxa, most OTUs belonged to Actinobacteria, and the most abundance belonged to Thaumarchaeota (Fig. 3a).Actinobacteria are adept at decomposing complex organic compounds (Goodfellow and Williams 1983), which has also been previously reported to be enriched in the soybean rhizosphere (Mendes et al. 2014).Thaumarchaeota includes numerous ammonia-oxidizing archaea (Leininger et al. 2006) and contributes to maintaining microbial stability under continuous cropping of soybeans (Liu et al. 2022).The enrichment of shared taxa was correlated with the higher pH and soil organic carbon in the rhizosphere (Fig. 6a), which was in line with previous studies highlighting the vital role of soil pH and the organic carbon provided by soybean roots in attracting specific bacterial taxa (Bulgarelli et al. 2013).
Site-specific taxa exhibited significant variations between Shijiazhuang and Xuzhou, and they were dominated by Actinobacteria in Shijiazhuang but by Nitrospirae in Xuzhou (Fig. 3b, d).This differentiation could be linked to the different alterations of soil properties in the soybean rhizosphere between the two sites (Table 1).In Shijiazhuang, soybean rhizosphere exhibited an increase in the concentration of available phosphorus, which contrasted with a decrease in Xuzhou.Meanwhile, the rhizosphere showed a reduction in ammonium content in Xuzhou, which was not observed in Shijiazhuang.Notably, both available phosphorus and ammonium contents were found to correlate significantly with the respective site-specific taxa at each site (Fig. 6b, c).Actinobacteria is known as a phylum that includes many taxa with phosphorus-solubilizing capabilities (Mitra et al. 2022), which could account for the elevated available phosphorus levels in soybean rhizosphere in Shijiazhuang.Meanwhile, Nitrospirae comprises nitrifying bacteria, and some strains within this phylum have even been reported to be capable of completing the entire nitrification process independently, converting ammonium to nitrate (Daims et al. 2015), which might explain the reduction of ammonium content in the rhizosphere in Xuzhou.Therefore, the differences in these soil properties may not be causal factors but rather the consequences of the distinct site-specific taxa of each site.
Variety-specific taxa exhibited a more even distribution across several phyla (Fig. 3c, e), revealing the specific preferences of different soybean varieties when selecting microbes.Interestingly, a large proportion of these taxa had low relative abundances, with 391 out of 480 OTUs in Shijiazhuang and 358 out of 460 OTUs in Xuzhou presenting a relative abundance of less than 0.1‰ (Supplemental Table S5).As reported by Chen et al. (2020), rare microbial taxa were the major drivers of soil multifunctionality, and these low-abundance taxa may play specialized but potentially crucial roles in N cycling.Besides the soil properties that were also correlated with the shared and site-specific taxa, the variety-specific taxa were uniquely correlated with available phosphorus in Shijiazhuang.As phosphorus is known as an immobile and limiting nutrient that can influence N fixation and other N-cycling processes (Tang et al. 2020), this correlation indicated that the variety-specific taxa might contribute to the N-cycling functions by enhancing the availability of necessary nutrients such as phosphorus.
The relationship between N-cycle-related functions and the shared, site-specific, and variety-specific taxa enriched in the soybean rhizosphere The shared and site-and variety-specific taxa were closely related to the potential N-cycling function of the bacterial community in the soybean rhizosphere, as indicated by the network, random forest, and Mantel test in this study (Figs. 4c,5c,and 6).This is consistent with previous studies reporting that rhizosphere-enriched taxa are intricately linked to several processes that are pivotal for crop growth, including the N cycle (Schmidt et al. 2019;Wei et al. 2015).However, taking the shared and site-and variety-specific taxa enriched in the soybean rhizosphere together, only 19 out of 1090 OTUs were annotated with N-cycle-related functions (Fig. 5b).This suggested that the rest of the OTUs in shared, site-, and variety-specific taxa might contribute to the N cycle in the soybean rhizosphere through interacting with these functional microbes even though they did not participate in N transformations themselves.
In the network with diazotrophs, positive correlations of the shared and site-and variety-specific taxa were found with only eight out of 39 diazotroph orders, including Rhizobiales (Fig. 4c).More negative correlations than positive correlations indicated complex competition between diazotrophs and other bacterial taxa (Faust and Raes 2012), which may help soybean plants recruit a microbiome that can enhance N fixation while limiting energy costs to support their nutrient needs.These results suggest that shared and site-and variety-specific taxa may contribute to this process by competing with diazotrophs that are less beneficial to soybean growth (Chepsergon and Moleleki 2023;Wheatley et al. 2020), potentially improving the efficiency of N fixation and reducing the need for chemical fertilizers.Furthermore, the variety-specific taxa exhibited a higher average degree than shared and site-specific taxa, indicating a strong, fine-tuned interaction with specific soybean genotypes, potentially facilitating optimized N fixation and utilization that is unique to each variety.Nevertheless, direct evidence for the improved N fixation of soybean varieties due to rhizosphere-enriched taxa is still needed.
The shared, site-specific, and variety-specific taxa did not show substantial differences in their impacts on the abundances of N-cycle-related genes (Fig. 5c).In contrast, their association with N-cycle-related genes was distinct between Shijiazhuang and Xuzhou, where they were negatively correlated with bacterial amoA and narG in Shijiazhuang but positively correlated with archaeal amoA in Xuzhou (Figs. 5c and 6).Despite the shared taxa being consistently enriched in the soybean rhizosphere across both sites, the ratio of ammonium (the substrate for ammonia-oxidation) to nitrate (the substrate for denitrification) was much lower in Shijiazhuang than in Xuzhou.This site-specific correlation indicated that they interacted with ammonia-oxidizing and denitrifying microbes differentially according to the relative levels of ammonium and nitrate.This result suggested their potential role in modulating N-cycle processes, which aligned with previous studies highlighting the role of the core microbiota in sustaining the stability of rhizosphere functionality (Jiao et al. 2022a, b).Site-specific taxa showed similar correlation patterns with N-cycle-related genes as shared taxa, suggesting that these taxa may act synergistically or complement the functions of the shared taxa, contributing to the stability and efficiency of N cycling within the rhizosphere (Saleem et al. 2019).A greater number of OTUs in variety-specific taxa were found in the most influential OTUs that affected the community N-cycling function in the random forest model at both sites, indicating that variety-specific taxa may play a relatively more important role in enhancing the N cycle.This underscored the influence of crop genotype on N-cycle processes (Favela et al. 2021), which may recruit unique taxa with redundant functions with those of the shared taxa, ensuring the robustness of important functions such as those involved in the N cycle.These results revealed the complex interplay between the shared, site-specific, and variety-specific taxa and N cycling.While the shared taxa provide a foundational contribution to maintaining N-cycle functionality, site-specific and varietyspecific taxa provide supplementary resilience and adaptability to environmental variability.
Overall, this study of the bacterial community and N-cycling function in the rhizosphere of five typical soybean varieties grown in Shijiazhuang and Xuzhou demonstrated distinct compositions but similar associations with the potential N-cycling function of the shared, site-specific, and variety-specific taxa recruited by soybean roots.The results advance our understanding of the resilience and adaptability of the rhizosphere microbiota and its functions and provide potential insights for developing strategies to harness the rhizosphere microbiome for sustainable agriculture.The underlying mechanisms of the complex interactions between rhizosphere-enriched taxa and the N-cycle-related microbes need further investigation.

Fig. 1
Fig. 1 Bacterial communities in bulk soil and the rhizosphere of the five soybean varieties at Shijiazhuang and Xuzhou.a The richness, Shannon, and PD whole tree indices.Different letters indicate significant differences at p < 0.05; b unconstrained PCoA of the whole

Fig. 2
Fig. 2 OTUs significantly (adjusted p < 0.1) enriched or depleted in the rhizosphere of each soybean variety compared to the bulk soil at each site (a) and Venn diagrams of the enriched OTUs (b)

Fig. 3
Fig. 3 Phylogenetic relationships and relative abundances of rhizosphere-enriched taxa shared across varieties and sites (a), site-specific taxa at b Shijiazhuang and d Xuzhou, and variety-specific taxa at c Shijiazhuang and e Xuzhou

Fig. 6
Fig. 6 Relationships of the shared, site-specific, and variety-specific taxa with soil properties and N-cycle-related genes at a both sites, b Shijiazhuang, and c Xuzhou.*, **, and ** indicate significant differences at p < 0.05, 0.01, and 0.001, respectively.AOA archaeal amoA,

Table 1
Selected soil properties in the bulk soil and the rhizosphere of the five soybean varieties at Shijiazhuang and Xuzhou