The previously isolated PSB A. pittii gp-1 (accession number: MK641660) with indole acetic acid production ability was used in potted experiment. The strain gp-1 was inoculated to 200 ml of the National Botanical Research Institute’s phosphate (NBRIP) medium and incubated at 28°C with shaking of 180 rpm for 5 days. NBRIP medium contained 10 g/L of glucose, 5 g/L of Ca₃(PO₄)₂, 0.25 g/L of MgSO₄⋅7H₂O, 5 g/L of MgCl₂⋅7H₂O, 0.2 g/L of KCl, 0.1 g/L of (NH₄)₂SO₄, and 2 ml/L of trace element solution (EDTA, 10 g/L; MnSO₄⋅H₂O, 2.2 g/L; FeSO₄⋅7H₂O, 1.0 g/L; CuSO₄⋅5H₂O, 0.5 g/L; CoCl₂⋅6H₂O, 0.3 g/L; Na₂MoO₄⋅2H₂O, 0.2 g/L; and CaCl₂, 0.1 g/L) (Nautiyal, 1999). After incubation, bacteria were collected by centrifuging and washed three times with sterile water.

The experimental potted soil was collected from an uncultivated field in Wuhan, China (30°28′N, 114°21′E). The soil type is calcareous, with original pH, total carbon, total nitrogen, availability phosphorus, and total phosphorus of 6.9, 0.52, 0.68%, 0.22 mg/g, and 0.89 mg/g, respectively. These P-deficient soils were sieved through a 2-mm mesh to remove stones and plant residuals. TCP was applied as phosphorus source in plant growth promotion experiment as described in previous literatures (Yu et al., 2011; Liu et al., 2014). Four potted treatments were designed: 200 g of sieved soil + 100 ml of sterile water (CK treatment), 195 g of sieved soil + 5 g of TCP + 100 ml of sterile water (Tri treatment), 200 g of sieved soil + 10 ml of bacterial suspension (10⁷ cfu/ml) + 90 ml of sterile water (Sup treatment), and 195 g of sieved soil + 5 g of TCP + 10 ml of bacterial suspension (10⁷ cfu/ml) + 90 ml of sterile water (Bac treatment). Each treatment had five replications. Soybean seeds (Glycine max w82) were purchased from China National Seed Group, pre-cultivated in sterile nutritious soils, and allowed them grow to about 10-cm length of sprouts. Each sprout with same growth potential was transplanted to each plastic pot as described above, and the strain gp-1 was inoculated to soybean rhizosphere in Sup and Bac treatments. Each plot was covered with Nylon membrane. These pots were randomly placed in greenhouse and incubated at 25°C with the cycling treatment of 16-h light and 8-h dark for a total of 40 days.

Every 10 days, we used alcohol-wiped shovels and tweezers to collected about 5 g of bulk soils near soybean root from each pot. In the experiment of plate-colony counting for abundance of PSB, 1 g of freeze-dried soil was added to 10 ml of sterile water and shaken at 180 rpm for 30 min, and the mixture is allowed to stand for 10 min. Then 1 ml of soil suspension was diluted, 0.1 ml × 10^–6 of diluent was evenly spread on NBRIP solid medium containing 0.2 g/L of cycloheximide acting as fungicide and incubated at 28°C for 5 days. After incubation, the cfu in different plates were counted. We also estimated content of indole acetic acid by using Van Urk Salkowski reagent, and the standard approach has been described previously (Biswas et al., 2018).

After 40-day growth of soybean, we excluded pots with the best and worst soybean growth in each treatment, and then 12 pots were left. We scraped rhizosphere soils by using a brush. We measured soil physicochemical properties, including pH, total carbon, total nitrogen, and available P, based on standard methods (Wan et al., 2021a). Microbial biomass P was evaluated by chloroform fumigation extraction and was calculated as the difference between fumigated and non-fumigated subsamples and simultaneously revised for the incomplete recovery of a P spike (Roberts et al., 2013; Ragot et al., 2016).

Soil alkaline phosphatase activity and phytase activity were determined according to previous methods (Wan et al., 2020a). Phosphatase activity and phytase activity were expressed as μg pNPP produced per gram of freeze-dried soil in 1 h and μmol P produced per gram of freeze-dried soil in 1 day, respectively.

The pots in each group was kept to measure the plant height, plant fresh weight, plant dry weight, leaf number, leaf fresh weight, root length, and root fresh weight. Soybean shoots and roots were separated from plants and dried at 60°C. The clean and dried root and shoot were separately cut into small pieces and digested by concentrated H₂SO₄–H₂O₂. The digested solutions were applied for measuring the content of root P and shoot P (Fraser et al., 2017).

Three rhizosphere soils from each group were used to extract total DNA using a DNA extraction kit (Mo Bio, Carlsbad, CA, United States) according to the manufacturer’s instruction. DNA concentrations were determined using a NanoDrop 2,000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). All extracted DNA samples were stored at –80°C.

The absolute abundances of phosphorus-cycling-related genes in soil bacteria were measured using qPCR with SYBR green mix. Primer sequences for amplifying P-cycling-related genes (i.e., phoD, bpp, gcd, and pstS) and quantitation PCR condition are summarized in Supplementary Method 1. Additionally, we used these primers to amplify bpp, phoD, gcd, and pstS from A. pittii gp-1.

The V3–V4 region of bacterial 16S rRNA gene was amplified using the primers 338F (5′-ACT CCT ACG GGA GGC AGC A-3′) and 806R (5′-GGA CTA CHV GGG TWT CTA AT-3′) (Mori et al., 2013). A PCR of 20 μl was performed in triplicate using a thermal cycler (ABI 9700, Thermo, United States) and conducted at the following conditions: an initial denaturation at 95°C for 3 min, 30 cycles of 95°C for 40 s, 58°C for 40 s, and 72°C for 50 s, and then a final extension at 72°C for 10 min. Sequencing was conducted on an Illumina MiSeq platform at Majorbio Bio-Pharm Technology Co., Ltd., Shanghai, China.

The raw reads were processed to gain purified sequences following the pathway of QIIME (Caporaso et al., 2010). We eliminated (1) sequences that did not exactly match barcodes and primers; (2) sequences with an average quality score < 20; (3) sequences with maximum homopolymers < 10 bp; and (4) sequences that contained ambiguous bases call. The purified sequences were clustered into operational taxonomic units (OTUs) at 97% identity against the SILVA v128 reference set.

Significant differences were calculated by the one-way analysis of variance with means compared using Tukey’s test in R. Venn diagram and non-metric multidimensional scaling (NMDS) plot were used to reflect bacterial community composition. Pairwise analysis of similarity (ANOSIM) was applied to quantitatively evaluate difference in bacterial community composition by using the ‘‘anosim’’ function in the ‘‘vegan’’ package of R. Permutational multivariate analysis of variance (PERMANOVA) was applied to evaluate pure effect of factors (e.g., physicochemical parameters and enzyme activity) on vegetation properties by using the ‘‘adonis’’ function in the ‘‘vegan’’ package of R. Linear discriminant analysis (LDA) effect size (LEfSe) statistical analysis was conducted on the online interface Galaxy¹ at a significant level of p < 0.05 and an LDA score > 4. Functional profiling of bacterial taxa was carried out by applying the “Tax4Fun2” package in R, and the functional redundancy index for each sample was calculated based on 16S rRNA gene similarity (Wemheuer et al., 2020). Canonical analysis of principal coordinates was applied to investigate influences of components including soil physicochemical parameters, gene abundance, cell exudates (include enzyme and indole acetic acid), and relative abundances of phylum bacteria on the vegetation properties. To identify core taxa, OTUs observed in more than 50% of all samples (> 6 samples, 875 OTUs) were applied to build a co-occurrence network. The co-occurrence network was visualized using Gephi v. 0.9.2² at a significant level of p < 0.01 and Spearman’s correlation coefficient higher than 0.67 (Wan et al., 2021b). Structural equation model was built to show relationships among vegetation properties, physicochemical properties, gene abundance, cell exudate, and bacterial community composition by using the packages of “sem” and “plspm” in R. The first principal component (PC1) value of soil physicochemical properties, P-cycling-related gene abundance, bacterial community composition, cell exudate, and vegetation properties accounting for 96.19, 85.19, 41.56, 98.99, and 96.37% of the total variances, respectively, were used as a proxy in structural equation model.

The PSB abundance represented by the number of cfu showed significant difference in four treatments (CK, Tri, Sup, and Bac) during 40-day growth of soybean (Figure 1A). The abundance of PSB in Bac treatment dramatically increased from 3.57 × 10⁷ cfu/g soil at day 10 to 6.96 × 10⁷ cfu/g soil at day 40 (p < 0.05). The population of PSB fluctuated in CK, Tri, and Sup treatments during 40 days but did not significantly ascend at day 40 than at day 10 (p > 0.05). The abundance of PSB in Bac treatment was significantly higher than that in other groups; this difference might be partially due to the input of A. pittii gp-1 and TCP. In addition, we randomly picked 10 colonies from the plate and found that 16S rRNA gene sequence of three bacterial colonies presented 100% similarity with that of A. pittii gp-1.

The indole acetic acid content was significantly higher in Bac treatment than in other treatments in each period (p < 0.05; Figure 1B). Additionally, the indole acetic acid content noticeably increased in Bac and Sup treatments during 40 days (p < 0.05), while in CK and Tri treatments, it did not (p > 0.05). Linear regression indicated that abundance of PSB was significantly correlated with content of indole acetic acid (Supplementary Figure 1). This suggests that PSB could produce and release indole acetic acid, which in turn might promote soybean growth.

After 40 days’ growth, the soybean presented erect leaves that became dark green (Supplementary Figure 2). Differences in vegetation properties were found in four treatments, including the plant height, plant fresh weight, plant dry weight, leaf number, leaf fresh weight, root length, root fresh weight, shoot P, and root P (Table 1). Plant dry weight, root length, shoot P, and root P were significantly higher in Bac group than in other three groups (p < 0.05). More importantly, the plant length, plant fresh weight, plant dry weight, leaf number, leaf fresh weight, root length, root fresh weight, root P, and shoot P were dramatically higher in Sup treatment than in CK treatment (p < 0.05). This suggests that the inoculation of PSB A. pittii gp-1 promotes soybean growth.

The soil pH (6.1–7.0) was significantly lower in Bac and Sup treatments than that in CK and Tri treatments (Table 1). Total carbon (0.43–2.04%), total nitrogen (0.07–0.24%), available P (0.19–0.99 mg/g), microbial biomass P (0.08–0.24 mg/g), alkaline phosphatase activity (3.79–10.25 μg/g/h), and phytase activity (0.66–1.63 μmol/g/day) were remarkably higher in Bac treatment than in other treatments (p < 0.05). These results indicate that the inoculation of PSB A. pittii gp-1 increases P availability and microbial activity.

Basically, the abundances of bpp, phoD, gcd, and psts genes were higher in Bac and Sup treatments than in CK and Tri treatments (Figure 1B). Expectedly, the abundance of gcd was significantly higher in Bac treatment (3.36 × 10⁷ copies/g soil) than that in CK treatment (1.86 × 10⁶ copies/g soil), Tri treatment (4.68 × 10⁶ copies/g soil), and Sup treatment (9.62 × 10⁶ copies/g soil). Linear regressions reflected significantly positive correlations between phoD gene abundance and alkaline phosphatase activity (R² = 0.585, p < 0.01), between bpp gene abundance and phytase activity (R² = 0.892, p < 0.001), and between gcd gene abundance and indole acetic acid content (R² = 0.854, p < 0.001) (Supplementary Figure 3). Additionally, gcd and pstS could be amplified from strain A. pittii gp-1 using primers described above, while bpp and phoD did not. These results might imply that the addition of A. pittii gp-1 could increase the abundances of organic P-cycling-related bacterial abundance.

A total of 2,829 OTUs were found across 12 soil samples. The CK, Tri, Sup, and Bac treatments possessed 1,670, 1,556, 1,413, and 906 OTUs, respectively; and they shared 181 OTUs (Figure 2A). A total of 39 phyla were observed, and 11 phyla with relative abundance > 0.01% were found across these 12 samples (Figure 2B). Proteobacteria, Chloroflexi, Actinobacteria, and Firmicutes were the first level dominant bacteria, with corresponding relative abundance from 7.94 to 48.22%, from 0.92 to 46.58%, from 8.86 to 44.47%, and from 1.13 to 25.25%, respectively. Acidobacteria, Bacteroidetes, Cyanobacteria, Deinococcus–Thermus, Gemmatimonadetes, Nitrospirae, and Saccharibacteria were the secondary dominant bacteria. The NMDS result showed that distinct difference in bacterial community composition among four treatments (Figure 2C). ANOSIM confirmed further the significant difference (R = 0.6451, p < 0.001). According to LEfSe result, bacteria including Bacillus and Acinetobacter were dramatically abundant in Bac treatment, while bacteria including Acinetobacter, Nitrospira, and Rhodobacter were significantly abundant in Sup treatment (Supplementary Figure 4). According to PERMANOVA results, the application of TCP explained 23.40% of the total variation in community composition (F = 5.27, p < 0.01), and the application of A. pittii gp-1 explained 29.11% of the total variation in community composition (F = 5.75, p < 0.001).

The bacterial community diversity represented by the Shannon–Wiener index (3.85–6.16) and community richness represented by Chao1 index (687–1303) were significantly lower in Bac treatment than in other treatments (p < 0.05; Supplementary Figure 5). This suggests that the addition of A. pittii gp-1 and TCP decreased rhizosphere bacterial diversity.

Based on functional profiling results, 3,113 functions at Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway level 3, including carbon-, nitrogen-, phosphorus-, and sulfate-cycling-related enzymes or proteins, displayed a higher functional redundancy in CK + Tri (without A. pittii gp-1 addition), whereas 3,772 functions had higher redundancies in Sup + Bac (with A. pittii gp-1 addition) (Figure 3A). It was worth noting that 206 functions representing P-cycling-related enzymes or proteins were higher in Sup + Bac than in CK + Tri, such as phosphoglycerate dehydrogenase (EC: 1.1.1.95) and phosphoglycerate kinase (EC: 2.7.2.3). Additionally, 35 functions [(e.g., L-iduronidase (EC: 3.2.1.76), dCTP deaminase (EC: 3.5.4.30), and phloroglucinol synthase (EC: 2.3.1.253)] were unique in CK + Tri, while 198 functions [e.g., phosphotransferase (EC: 2.7.1.-), neamine phosphoribosyltransferase (EC: 2.4.2.49), 5-phosphoribostamycin phosphatase (EC: 3.1.3.88), and uracil phosphatase (EC: 3.1.3.104)] were exclusive in Sup + Bac. At KEGG level 2, some functions (e.g., metabolism of cofactors and vitamins, energy metabolism, and translation) were significantly higher in Sup + Bac than in CK + Tri (p < 0.05), but some functions were not (Figure 3B).

A co-occurrence network was constructed to reveal the relationships among bacterial taxa (Figure 4A). We found 50,510 positive edges (represent significantly positive correlation) and two negative edges (denote dramatically negative correlation), suggesting that rhizosphere bacteria presented a less conflicting interaction. We also clarified the top 20 core nodes; i.e., those with the highest betweenness centrality were affiliated with Acidobacteria (e.g., OTU522), Actinobacteria (e.g., OTU947), Chloroflexi (OTU67), Firmicutes (e.g., OTU601), Gemmatimonadetes (OTU1967), and Proteobacteria (e.g., OTU1813) (Figure 4B).

According to PERMANOVA results, the application of TCP could explain 13.69% of the total variation (F = 20.23, p < 0.01) in vegetation properties, while the application of A. pittii gp-1 could explain 72.41% of the total variation (F = 107.05, p < 0.001). According to results of canonical analysis of principal coordinates, soil physicochemical properties (Figure 5A), gene abundance (Figure 5B), cell exudates (Figure 5C), and relative abundances of bacterial phyla (Figure 5D) explained more than 80% of the total variation in vegetation properties. Physicochemical parameter, gene abundance, enzyme activity and IAA, and relative abundance of bacterial phylum showed significantly pure effects on vegetation properties based on PERMANOVA (Figure 5).

Additionally, we also found that bacterial functions based on functional profiling were responsible for vegetation properties (Figure 3C). The function of metabolism of cofactors and vitamins (R² = 79.75%, F = 39.38; p < 0.01) showed greater effect on vegetation properties than other functions according to PERMANOVA results. The core taxa identified from co-occurrence network also have significant effects on vegetation properties based on PERMANOVA (Figure 4B). The OTU1813 regarded as Acinetobacter genus presented higher influence (R² = 52.08%, F = 10.87; p < 0.01) than other core taxa.

Ultimately, we used structural equation model to reveal interconnections among soil physicochemical properties, P-cycling-related gene abundance, bacterial community composition, enzyme activity, and vegetation properties (Figure 6). The model presented a good fit to our data, as indicated by the non-significant χ²-test (N = 12, χ² = 0.707, d.f. = 1, p = 0.400). On the one hand, bacterial community could affect soil physicochemical properties and P-cycling-related gene abundance, which in turn affect vegetation properties; on the other hand, soil physicochemical properties and P-cycling-related gene abundance could influence enzyme activity, which in turn influences vegetation properties. These results indicated that soil, plant, and bacteria presented close relationships.

Applying PSB can increase soil available P content (Maougal et al., 2014) and promote vegetation growth (Yu et al., 2011; Biswas et al., 2018; Hansen et al., 2020). However, it should be considered whether the PSB could maintain their activity, function, and abundance after inoculation. In this study, the inoculation of PSB A. pittii gp-1 significantly promoted the growth of soybean represented by better vegetation properties, which is in accordance with prior findings describing that PSB can enhance the growth of legume plant (Collavino et al., 2010; Bononi et al., 2020; Cumpa-Velásquez et al., 2021) and other kinds of plants (Yu et al., 2011; Liu et al., 2014). In these studies, the increase in the content of available P or small molecular organic acid is closely correlated with the growth of plants. The PSB Acinetobacter genus is reported to have the ability to release small molecular organic acid (e.g., indole acetic acid, gluconic acid, oxalic acid, and citric acid) (Marra et al., 2012; Marwa et al., 2019; Rasul et al., 2019). Interestingly, we found the A. pittii gp-1 could produce indole acetic acid detected by using the Van Urk Salkowski reagent. Therefore, the inoculation of the A. pittii gp-1 might increase the content of soil organic acid, which in turn increased the content of available P. Additionally, we detected Acinetobacter genus in Bac treatment by using simple 16S rRNA gene sequencing for single colony. Illumina MiSeq sequencing result also reflected that Acinetobacter dominated in Sup and Bac groups. These results suggest that A. pittii gp-1 could survive after inoculation and could promote the growth of soybean.

In addition, we used four pairs of primers as described above to amplify bpp, phoD, gcd, and pstS genes from A. pittii gp-1. Unexpectedly, only gcd and pstS genes could be amplified. Previous studies have reported that Acinetobacter genus harbors gcd and pstS gene (Marra et al., 2012; Farrugia et al., 2015; Wan et al., 2020b), and almost no study has reported that Acinetobacter genus possesses bpp and phoD genes. However, the abundances of P-cycling-related genes including bpp, phoD, gcd, and pstS were higher in Sup and Bac treatments. These results and findings suggest that the inoculation of PSB A. pittii gp-1 might significantly increase both inorganic and organic P-cycling-related gene abundance of soil indigenous bacteria. This phenomenon might be due to the solubilization of inorganic P by added PSB A. pittii gp-1 via releasing small molecular organic acid. Consequently, part of soluble P was assimilated by native bpp-harboring bacteria and phoD-harboring bacteria and in turn enriched the abundances of bpp and phoD genes and released more phosphatase and phytase. In addition, a part of inoculated A. pittii gp-1 might die; thus, the cell residues could be treated as nutrient for indigenous microorganisms. Previous literatures have reported that gcd-harboring bacteria can produce and release small organic acid to solubilize insoluble inorganic P, thus promoting the growth of plant (Wagh et al., 2014; Rasul et al., 2019). The bpp-harboring bacteria and phoD-harboring bacteria are reported to be responsible for the turnover of soil organic P by releasing extracellular enzyme, which in return promotes the growth of vegetation (Maougal et al., 2014; Ragot et al., 2016; Hu et al., 2018; Zhang et al., 2021). Therefore, the application of PSB A. pittii gp-1 could enhance utilization potentials of both inorganic and organic P.

Considering community diversity is closely correlated with soil ecosystem functions (Wan et al., 2021c), it is important to decipher effects of the application of PSB on plant rhizosphere bacterial community. We found significant decrease in rhizosphere bacterial diversity and distinct change in bacterial community composition, which is similar to findings in published literatures (Estrada-Bonilla et al., 2017; Wei et al., 2017; Widdig et al., 2019). In addition, earlier studies have reported that vegetation also affects the composition of bacterial community (Xue et al., 2017; Yang et al., 2018; Campos-Herrera et al., 2019). To the best of our knowledge, this is the first report that the addition of PSB A. pittii could promote the community function of rhizosphere bacteria especially phosphorus-cycling-related functions. This phenomenon might be due to elevated nutrient caused by inoculation of PSB A. pittii, which in turn affected rhizosphere bacterial community composition and function. An earlier study has reported that dead bacteria can be treated as available nutrient to affect growth of other microorganisms (Hanajima et al., 2019). Additionally, microbial biomass P contributes to P solubility in riparian vegetated buffer strip soils (Roberts et al., 2013).

Based on these results and findings, we raised one question of whether there were close relationships among plant, soil, and rhizosphere. The structural equation model reflected stronger interconnections among vegetation properties, soil physicochemical properties, P-cycling-related gene abundance, cell exudates, and bacterial community composition. This result is similar to our prior finding (Wan et al., 2021a). The co-occurrence network also showed that core taxa belonging to Acidobacteria, Chloroflexi, Gemmatimonadetes, and Proteobacteria presented significant effects on vegetation properties. Previous literature has reported that some specific phylum bacteria, such as Acidobacteria, Actinobacteria, and Proteobacteria, are responsible for vegetation growth under different P conditions (Bergkemper et al., 2016). Vegetation properties and microbes could also affect each other (Neal et al., 2017; Yang et al., 2018; Muñoz et al., 2021), suggesting that soil, plant, and bacteria have close relationships. In the future, we will explore molecular mechanisms to reveal close interconnections among soil, plant, and bacteria.