P. asiatica JP233, isolated from the soil of vegetable greenhouse in Shouguang, Shandong Province, has been proved to be a potentially efficient Phosphate-solubilizing bacteria (PSB) (Wang F. et al., 2022; Wang L. L. et al., 2022; Wang Z. H. et al., 2022; Yu et al., 2022). Bacterial strains were streaked and purified on Luria-Bertani (LB) solid medium, single colonies were selected and transferred to LB liquid medium for shake culture at 180 rpm and 28°C for 24 h. The bacterial culture was then centrifuged, with the supernatants being discarded, and washed and re-suspended with sterile water as inoculant. Two kinds of soil with different phosphorus (P) (Supplementary Tables S1, S2), high-P soil and low-P soil, were used in the plant pot experiment, and Dutch hard powder tomato varieties were purchased from local suppliers.

The experiment involved 12 treatments, utilizing two different P levels soils [low-P (L) and high-P (H)]. Three levels of P application at 0, 50, and 100 mg kg⁻¹ were applied under varying P levels soils. The experimental group was inoculated with the JP233 strain, while the control group received an equal amount of water. Each treatment included 12 biological replicates. The seeds were disinfected with alcohol and sodium hypochlorite, then placed in sterile water-soaked petri dishes at 28°C to aid germination. Subsequently, the tomato seedlings were transplanted into planting cups. When the tomato reached true leaf stage, the experimental group was inoculated with JP233 bacterial solution by root irrigation, and each plant received a 6 mL bacterial suspension with 10⁸ CFU mL⁻¹. The control group received the same amount of sterile water. The pot experiment was run in a controlled environment at 25°C in full light in the greenhouse. Samples were collected at 7, 14, 21, 28 days post-inoculation of the bacterial solution to measure plant height, stem diameter, fresh weight/dry weight of the above/underground plant parts, total P content in the plant parts, total phosphorus (TP) content in the soil, and available phosphorus (AP) content in the soil. Plant total P was analyzed using the H₂SO₄-H₂O₂ digestion method and vanadium molybdate blue colorimetric method (Bao, 2000), soil TP was determined using the potassium persulfate digestion method (Su and Chen, 2010), and soil AP was assessed using the Olsen method (Bao, 2000).

The H0 (0 mg kg⁻¹ P added in high P soil) level, which is the most significant level of increase in soil AP content, was selected and the tomato rhizosphere soil was sampled for metagenomic sequencing analysis. The rhizosphere soil was collected by taking tomato root from its pot and gently shaking off the loosely attached bulk soil. The soil firmly adhering to the roots, with a thickness of approximately 1 mm, was deemed as rhizosphere soil and gathered through a washing and centrifugation process. DNA from rhizosphere soil samples of H0 tomato was isolated and purified with magnetic bead genomic DNA extraction kit (Beijing, Bioteke, China). Using the truSeq nano DNA LT library preparation kit (Illumina, United States), a DNA library was created. LC Bio Technology CO., Ltd. (Hangzhou, China) used the Illumina NovaseqTM 6000 platform to sequence DNA libraries. PE150 was the sequencing mode used. Using the fastqc program,¹ quality control of the original data was carried out. To eliminate host contamination, reads were aligned to the host genome using bowtie2.² To ensure that the subsequent assembly and analysis results were microbial sequences (Langmead and Salzberg, 2012). Valid data of each sample assembled by MEGAHIT software³ to acquire FASTA format file (Li et al., 2015). Coding Region (CDS) prediction was carried out using MetaGeneMark software to retain contigs sequences with a length greater than 500 bp. Sequences with CDS length less than 100 nt were filtered according to the predicted results (Zhu et al., 2010; Gurevich et al., 2013). To extract non-redundant unigenes, CD-HIT⁴ was used to cluster the CDS sequences of all samples. In the meantime, the longest sequence was chosen as the representative sequence to build the Unigenes collection, and identity 95% and coverage 90% were employed for clustering (Fu et al., 2012). The different annotation information taxonomic levels of species were obtained by using DIAMOND software⁵ to compare the Unigenes protein sequence and Non-Redundant Protein Sequence Database (NR) meta library. Compared with Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO), Evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG), functional databases to obtain annotation information of each functional database (Buchfink et al., 2015). Finally, the abundance informations of various species and functional taxonomic levels were obtained based on the abundance statistics of Unigenes. Statistical analysis was performed based on species, function and Unigenes among the comparison groups. The difference threshold was p < 0.05 and |log₂(fold change)| ≥ 1 (log₂FC), which was applied for subsequent analysis.

Based on the genes related to P cycle in soil microorganisms reported in the literature (Dai et al., 2020; Liang et al., 2020; Huang K. et al., 2024; Huang Y. L. et al., 2024; Wang et al., 2024), a total of 33 genes involved in P cycle were collected by comparing the gene name and functional classification with KEGG annotation information of metagenomic clustering Unigene set, excluding genes involved in intracellular phosphatase production in microbial metabolic activities. According to the function of genes, they can be divided into three categories: (1) Genes involved in P-starvation response regulation; (2) Genes involved in P-uptake and transport system; (3) Genes involved in inorganic P-solubilization; (4) Genes involved in organic P-mineralization. The names, functions and classifications of genes related to P cycle in soil microorganisms are shown in (Supplementary Table S3). Based on the corresponding relationship between species and function in the samples, the correlation between species abundance and functional abundance was analyzed. DIAMOND software was used to compare Unigenes protein sequence with NR database (blastp, evalue≤1e−5). Each Unigene with the best comparison result index was selected as the species classification (Buchfink et al., 2015), and Unigene, species (all hierarchical species) and function (all hierarchical annotation results) were annotated accordingly.

Student’s t test (T-test) and homogeneity of Variance Test were performed using IBM SPSS 25.0 (IBM Corporation, Armonk, NY, United States) to examine the significant differences between different P treatments. Metagenomic data using OmicStudio tools for bioinformatics analysis.⁶ The functional abundance tables of all Unigenes and eight functional databases were statistically analyzed by Mann Whitney U test, and the threshold value of |log₂FC| ≥ 1, (p < 0.05) was used to determine the significant differences. The four types of genes involved in P-cycling were mapped with Unigenes and visualized in the volcano map. The results of KEGG Pathway Definition functional nonparametric difference analysis were visualized in STAMP plot. The gene abundance of four major types of microorganisms involved in P-cycling, as identified in metagenomic data sets, was normalized using Z-score transformation in R (Wu et al., 2022). The differences were then evaluated using a two-sample T-test, and GraphPad Prism 10.0 (GraphPad Software, La Jolla, CA, United States) was used to show the results. Spearman was used for correlation calculation, and the correlation threshold was |rho| > 0.5 and p < 0.05. The correlation network map was used to compare the complexity of gene association between the two groups.

To investigate the impacts of P status on PSB, strain JP233 was inoculated to tomato plants cultivated in two kinds of soil with low (L) and high (H) P background. Different P fertilization levels were also included for both soils by amending 0, 50, 100 mg kg⁻¹ of KH₂PO₄ (designated as L0, L50, L100, H0, H50, H100). In low P soil, JP233 could enhance soil AP within 28 days post inoculation (dpi) across all three P fertilization levels (Figure 1A). Generally, the growth of tomato was much better under L50 and L100 compared to L0. Although JP233 could also promote plant growth, the increase of aboveground dry weight was not statistically significant after 14 dpi (Figure 1B; Supplementary Table S4). Strain JP233 could increase the P content of aboveground plant, and the increment was significant within 7 dpi (Figure 1C). Similarly, it could also increase the P content of tomato root, albeit the increment was not significant (Figure 1D). Notably, the root P content under L0 was much higher than the roots under L50 and L100 after 14 dpi, potentially due to the stressed growth conditions of the tomato plants under L0 (Figures 1B,D).

In high P soil (Supplementary Table S4), strain JP233 could also enhance soil AP during 28 dpi across all three P fertilization levels (Figure 2A). It is worth noting that extra P fertilization was not beneficial for tomato plant growth, as showed by the reduced aboveground dry weight under H50 and H100 compared with H0 (Figure 2B; Supplementary Table S5). Similar to in low P soil, JP233 could also promote plant growth and enhance the P content of aboveground plant parts and root, albeit the increment was not significant after 7 dpi (Figures 2B–D; Supplementary Table S5). After a comprehensive comparison, the H0 treatment was selected for further metagenomics analysis, as JP233 demonstrated effective P solubilization and plant growth promotion, and the influence of exogenous P could also be minimized under H0 conditions.

As indicated by Chao1, Shannon and Simpson indices, JP233 inoculation did not significantly change the α-diversity of microbial community (Supplementary Figure S1). However, the Anosim analysis revealed that the JP233 inoculation significantly altered the β-diversity (p < 0.05). The Principal Coordinate Analysis (PCoA) showed that JP233-treated soil samples clustered separately from the control (CK), with PCo1 explained 30.11% variance (Figure 3A). Species annotation based on Unigene sequence assembly discovered that the majority of Unigenes in tomato rhizosphere soil (>83%) originated from bacteria. Following JP233 inoculation, the relative abundance of bacteria was enhanced, while the abundance of viruses and eukaryotes were both significantly depressed (p < 0.05) (Figures 3B,C). The top 20 abundant species that were altered by JP233 are shown in Figure 3D. As expected, the relative abundance of P. asiatica was found enhanced significantly (p < 0.05). Pseudomonas putida was also found to be elevated, but there is reason to speculate that the reads may have originated from JP233, given the very close phylogenetic relationship between P. putida and P. asiatica (Wang F. et al., 2022; Wang L. L. et al., 2022; Wang Z. H. et al., 2022). It is noteworthy that the relative abundance of Bdellovibrio bacteriovorus significantly increased after JP233 inoculation (Figure 3D). Given that B. bacteriovorus is a predatory bacterium that invades and preys on other Gram-negative bacteria (Strauch et al., 2007), the increase has a detrimental effect on the survival of JP233.

The inoculation of JP233 induced modifications in the rhizosphere soil gene profiling, resulting in a total of 134,184 unigenes with increased abundance and 164,550 unigenes with decreased abundance (Supplementary Figure S2). According to KEGG annotation, these unigenes with significant abundance changes can be categorized into 33 pathways (Figure 4A). Among them, biosynthesis of cofactors, fatty acid metabolism, and ribosome, among others, were enhanced, whereas pathways such as monobactam biosynthesis, lysine biosynthesis, ascorbate and aldarate metabolism exhibited decreases. Enrichment analysis revealed that the unigenes exhibiting significant abundance alterations were predominantly enriched in pathways including the biosynthesis of secondary metabolites, amino acids, carbon metabolism, and others (Figure 4B). Given the advantageous roles of secondary metabolites, amino acids, and organic carbon in fostering plant health and growth, these alterations in gene profiling could potentially be attributed to the enhancement of plant growth.

The processes of the soil P cycle are primarily governed by four distinct gene groups, which are associated with the regulation of P-starvation response, P-uptake and transport, inorganic P-solubilization, and organic P-mineralization (Liu et al., 2023). Based on KEGG annotation, our focus was on 33 genes involved in P cycle (Supplementary Table S3), excluding genes related to intracellular phosphatase production during microbial metabolism (Bergkemper et al., 2016). Notably, the inoculation of JP233 resulted in a significant 3.4% increase in the total relative abundance of P cycle genes (p < 0.05). Furthermore, approximately 79% of these P cycle genes exhibited a heightened abundance. According to PCoA analysis, there were significant difference in the profiling of soil P cycle genes attributable to the inoculation of JP233 (Figure 5A).

The most prominent P cycle genes identified in the soil samples were phoR, pstB, ppk, pstS and gcd (Figure 5B). Specifically, phoR, which encodes the phosphate regulon sensor histidine kinase, is associated with the P-starvation response. The pstB and pstS genes encode components of the high-affinity phosphate-specific transporter PstSABC, while the ppk gene, which encodes polyphosphate kinase, is also involved in P-uptake and transport. Additionally, the gcd gene, encoding quinoprotein glucose dehydrogenase, is linked to inorganic P-solubilization through gluconic acid production. Despite JP233 possessing the gcd gene, it did not significantly increase its abundance in the soil (p > 0.05). However, JP233 did enhance the abundance of the ppx gene, which encodes exopolyphosphatase (PPX), playing a crucial role in the degradation of inorganic polyphosphate into phosphate (Song et al., 2020). Additionally, JP233 increased the abundances of two genes related to organic P-mineralization: ugpQ and phnO, which encode glycerophosphodiester phosphodiesterase (GDPD) and a C-P lyase subunit, respectively (p < 0.05). Conversely, JP233 also decreased the abundances of two P cycle genes, ugpC and phoN, which encode components of the Ugp transport system and C-P lyase, respectively (p < 0.05).

In addition to the individual P cycle genes, the inoculation of JP233 also modified their co-occurrence network. Specifically, JP233 boosted the connectivity among soil P cycle genes. Comparing with the JP233-treated soil samples, the P gene network in control (CK) soil exhibited lower complexity and sparser linkage density (Figure 5C). The inoculation of JP233 significantly improved key topological parameters, including the average degree and the interconnections among P cycle genes. Notably, genes implicated in P-uptake and transport (phnCDE, pstB), inorganic P-solubilization and organic P-mineralization (gcd, phnJLHMNP), as well as the regulation of P-starvation response (phoB) demonstrated a substantial increase in node degree and interaction patterns (Figure 5D).

To ascertain the origins of the alterations in the abundance and co-occurrence network of soil P cycle genes following the inoculation of JP233, their taxonomic compositions were traced. Figure 6 displays the top 20 most abundant genera of the P cycle genes that underwent significant alterations due to JP233. Notably, the abundances of the ppx gene from Sphingopyxis, Sphingosinicella, and Bdellovibrio increased significantly (p < 0.05). Similarly, the abundances of the ugpQ gene from Sphingopyxis and Marinobacter also showed a significant increase (p < 0.05). However, for the genes phnO, ugpC, and phoN, no individual genus within the top 20 showed significant alterations. While the abundances of other P cycle genes remained largely unchanged, their taxonomic compositions did shift, as evidenced by PCoA analysis (Figure 6F).