The study site is in Hongyuan County (32° 49.823′ N, 102° 35.237′ E), eastern Qinghai-Tibetan Plateau with the typical continental plateau climate. The mean annual temperature (MAT) is 1.1°C and the mean annual precipitation (MAP) ranges from 650 mm to 800 mm. The soil type is Felty soil following the Chinese Soil Classification System. The aboveground plant species mainly include Elymus nutans, Agrostis hugoniana, Kobresia setchwanensis, Kobresia humilis, Saussurea nigrescens etc.

The field experiment was initiated in 2019 and included the following seven treatments with three replicated plots: (1) P0, none fertilization; (2) P10, 10 g⋅m^–2⋅a^–1 (single superphosphate, containing P₂O₅ 16.0%); (3) P20, 20 g⋅m^–2⋅a^–1; (4) P30, 30 g⋅m^–2⋅a^–1; (5) P40, 40 g⋅m^–2⋅a^–1; (6) P50, 50 g⋅m^–2⋅a^–1; and (7) P60, 60 g⋅m^–2⋅a^–1. The P fertilizer was distributed evenly by hand on the soil surface in each plot once a year during natural rain event in early growing season on 6 May 2019 and 4 May 2020. For each replicate, a random quadrat (1 m × 1 m) was assigned as the sampling quadrat, and five scattered soil cores were collected from each quadrat with a 5 cm auger from the topsoil (0∼10 cm) after harvesting on August 15, 2020. All five soil cores from each quadrat were combined into one composite sample and mixed thoroughly to minimize sampling errors. All plant litters and rocks were removed by hand; 50 g of soil were stored at −80°C for DNA extraction; 350 g of soil were stored at 4°C for microbial biomass phosphorus (MBP) analysis within 1 week and the remaining sampled soil was air dried for physicochemical analysis.

Soil moisture content (MC) was determined gravimetrically at 105°C for 24 h and expressed as a ratio of soil water to dry soil weight. Soil pH was measured in the soil-water slurry (soil: water = 1: 2.5) using a pH meter (FE20-FiveEasyTM pH, Mettler Toledo, Germany). The total carbon (TC) and total nitrogen (TN) contents were measured using an elemental analyzer (Vario MAX, Elementar, Germany). Fe, Al, and Ca in acid digested extract were measured by ICP-OES (PerkinElmer Optima 8300, USA). Soil MBP was determined using the chloroform fumigation-extraction method (Wu et al., 2000).

The modified Hedley sequential fractionation method was used to measure soil P forms (Tiessen and Moir, 1993). Air dried soils were milled and sieved to 2 mm. Then 0.5 g soil was extracted sequentially with an anion exchange membrane in deionized water, 30 ml of 0.5 M NaHCO₃ at pH 8.5, 30 ml of 0.1 M NaOH and 30 ml of 1 M HCl in a 50 ml screw up centrifuge plastic tube. In each extraction step, samples were shaken continuously for 16 h at 25°C, centrifuged at 3,500 rpm for 15 min at 0°C and thereafter the supernatant was filtered using a 0.45 μm filter paper. The supernatant was stored at 4°C for phosphate determination. Soil residue in the centrifuge plastic tube was then re-suspended using concentrated H₂SO₄ and 30% H₂O₂ to extract more chemically stable P (Residual-P). The supernatant was analyzed for dissolved Pi and total P (Pt). Pi was determined directly in the supernatant through the molybdate-ascorbic acid method using a spectrophotometer at 712 nm (Murphy and Riley, 1962), Pt was measured through the same method after digestion with ammonium persulfate and H₂SO₄ at 121°C in an autoclave. The Po was then calculated as the difference between Pt and Pi. We omitted the addition of hot concentrated HCl from the original methodology. Soil Bio-P was the sum of Resin-Pi and NaHCO₃-Pi, which is the biologically most available Pi form for plants (Moro et al., 2015; Tian et al., 2022).

Soil DNA was extracted from the 0.25 g soil using the PowerSoil DNA Isolation Kit (MOBIO Laboratories, Carlsbad, CA, USA) as directed by the manufacturer within 2 weeks after sampling. The quantity and quality of the extracted DNA were measured using NanoDrop ND-2000 (Thermo Fisher Scientific, Waltham, MA, USA). The primer set phoD-F733 (5′-TGGGAYGATCAYGARGT-3′)/phoD-R1083 (5′-CTGSGCSAKSACRTTCCA-3′) (Ragot et al., 2015) was employed to target the microbial phoD gene. The PCR mixture (25 μl) contained 2× Taq MasterMix (TaKaRa, Dalian, China) 12.5 μl, Forward primer 1 μl, Reverse primer 1 μl, ddH₂O 9.5 μl, DNA (10 ng⋅μl^–1) 1 μl. The PCR amplification program included an initial denaturation at 95°C for 3 min, followed by 35 cycles of 95°C for 40 s, 59°C for 60 s, and 72°C for 60 s, and a final extension at 72°C for 10 min. Triplicate reaction mixtures per sample were pooled and purified using DNA Gel Extraction Kit (Tsingke, Beijing, China). Finally, the amplicons of all samples were pooled at equimolar concentrations, and sequencing was performed using the Illumina MiSeq (Illumina, San Diego, CA, USA) sequencing platform.

The raw sequences of the phoD gene were merged using FLASH with 30 bp overlap (Magoc and Salzberg, 2011), and the joint sequences were processed using the Quantitative Insights Into Microbial Ecology (version 2.0, QIIME2) pipeline (Bolyen et al., 2019). Reads were shorter than 300 bp, reads containing more than two ambiguous bases, and reads with an average quality score lower than 30 were removed. The chimeras were removed and the FrameBot was adopted to correct frameshift mutation sequences in the Ribosomal Database Project (RDA) function gene pipeline (Wang et al., 2013). Using the deblurring algorithm, 10,739 amplicon sequence variants (ASVs) were obtained from the 1,238,300 successfully corrected sequences in 21 samples. The RDP FunGene database was used to assign taxonomic affiliations of phoD sequences (Cole et al., 2014). To eliminate the bias on diversity comparison as caused by uneven sequencing, we rarified sequences to 25,690 for each sample randomly.

One-way ANOVA was performed to examine the effect of P application on soil physicochemical properties, P forms, phoD-harboring microbial diversity, and the relative abundance of the dominant phyla. The comparison of significant means was carried out using Tukey’s Honest Significant Difference (HSD) test. To visualize the community structure of the phoD-harboring microbes among different treatments, the Non-metric multidimensional scaling (NMDS) analysis was constructed based on the Bray--Curtis distance. The significance of the phoD-harboring community structure was tested by Permutational multivariate analysis (PERMANOVA), analyses of similarities (ANOSIM) and Multi Response Permutation Procedure (MRPP) based on 999 permutations with the Vegan package in the R environment.¹

Differentially abundant taxa among the seven treatments were identified using the linear discriminant analysis (LDA) effect size (LEfSe) method using online Galaxy,² and the LEfSe algorithm uses a P-value < 0.05 and an LDA score > 2.0 for the non-parametric factorial Kruskal–Wallis test.

Network analyses were performed to better understand the species relations within the phoD-harboring microbial communities in R as described by Mendes et al. (2014). We calculated all possible Spearman’s rank correlation coefficients, and the high correlations with cutoff at r > 0.8 and statistically significant P-value < 0.01 were kept in the network construction. The networks were constructed using the igraph package. The nodes in the reconstructed network represent ASVs, and the edges represent high and significant correlations between nodes. The network graphs were visualized based on nodes, edges and taxonomic information using the ggraph package. Moreover, the within-module connectivity (Zi) and among-module connectivity (Pi) were calculated using ggClusterNet package. According to Zi = 2.5 and Pi = 0.62, all nodes were divided into four categories: peripherals (Zi < 2.5, Pi < 0.62), connectors (Zi < 2.5, Pi ≥ 0.62), module hubs (Zi ≥ 2.5, Pi < 0.62), and network hubs (Zi ≥ 2.5, Pi ≥ 0.62). From the ecological perspective, connectors, module hubs, and network hubs are the generalists while peripherals are the specialists (Deng et al., 2012; Feng et al., 2022).

The correlations between phoD-harboring microbial community composition and environmental factors were detected by Mantel tests using the linkET package in R. Then, the structural equation modeling (SEM) framework was further applied to test the significance of hypothesized causal relationships between soil properties, P forms and phoD-harboring microbial community structure using AMOS 23.0 (Amos Development, Spring House, Armonk, NY, USA). To simplify the SEM model, the principal component analysis (PCA) was used to create multivariate functional indices for (NMDS1 and NMDS2) phoD-harboring microbial community structure. The first principal component axis (PCA1), explaining 71.24% of the total variance in phoD-harboring microbial community structure, were used in the SEM analysis. The best-fit model was derived using maximum likelihood based on model fit, i.e., chi-square test (χ²), goodness-of-fit index (GFI), and the root mean square error of approximation (RMSEA) (Hooper et al., 2008).

One-way ANOVA analysis showed that soil properties except for soil MC and Ca responded significantly to the P application (P < 0.05) (Table 1). The NaHCO₃-Po (Figure 1C) and pH (Figure 1I) declined with increasing P application rate, with the lowest values in the P60 treatment. However, the NaOH-Po (Figure 1E), Al (Figure 1N), Fe (Figure 1O), and MBP (Figure 1P) showed a rough increase pattern with increasing P application rate, with the highest value in the P60 treatment. Soil MC (Figure 1L) and Ca (Figure 1M) varied slightly with increasing P application rate. Soil TN (Figure 1J) and TC (Figure 1K) in the P0 treatment were significantly higher than that in P10∼P40 treatments. Additionally, soil Resin-Pi (Figure 1A), NaHCO₃-Pi (Figure 1B), NaOH-Pi (Figure 1D), HCl-Pi (Figure 1F), Residual-P (Figure 1G), and Bio-P (Figure 1H) all peaked at P30 treatment. This indicated that soil Resin-Pi, NaHCO₃-Pi, NaOH-Pi, HCl-Pi, Residual-P, and Bio-P showed non-linear relationships with increasing P application amount, instead, there might be a tipping point at the P application rate of 30 g⋅m^–2⋅a^–1.

Pearson analysis indicated that soil TC was the most relevant factor affecting soil Resin-Pi (r = −0.77, P < 0.01), NaOH-Pi (r = −0.78, P < 0.01), HCl-Pi (r = −0.74, P < 0.01), and Bio-P (r = −0.74, P < 0.01), while soil TN was the strongest correlate of NaHCO₃-Pi (r = −0.82, P < 0.01) (Figure 6). Soil NaOH-Po was negatively correlated with pH (r = −0.57, P < 0.01) but positively correlated with Al (r = 0.69, P < 0.01) and Fe (r = 0.73, P < 0.01). Soil Residual-P was also positively correlated with Al (r = 0.70, P < 0.01) and Fe (r = 0.69, P < 0.01). The links between phoD-harboring microbial community structure and environmental factors were detected by Mantel tests. The phoD-harboring microbial community structure was significantly correlated with NaOH-Po (r = 0.15, P < 0.05), pH (r = 0.18, P < 0.05), and Fe content (r = 0.15, P < 0.05).

Based on the most relevant factors in Mantel analysis and previous studies (Chen et al., 2019; Xu et al., 2021; Wang M. M. et al., 2022), the NaOH-Po, Fe, pH, TC, and Bio-P were included in the optimal SEM analysis. The model explained 78% of the variance in phoD-harboring microbial community structure (Figure 7A). The Fe (standardized direct effects: −0.61) and NaOH-Po content (standardized direct effects: −0.41) negatively and directly affected phoD-harboring microbial community structure, while Bio-P positively and directly affected phoD-harboring microbial community structure (standardized direct effects: 0.60). Standardized total effects derived from the SEM revealed that Bio-P had the largest positive effect (0.60) on soil phoD-harboring microbial community structure (Figure 7B). However, the standardized total negative effects on phoD-harboring microbial community structure were in the order of P application (−0.65) > TC (−0.55) > Fe (−0.37) > NaOH-Po (−0.31) > pH (−0.03).

Fertilization is the main way of exogenous P input to soil, and the content of soil P forms is directly affected by P application (Mahmood et al., 2020). Soil Pi (Resin-Pi, NaHCO₃-Pi, NaOH-Pi, HCl-Pi, and Bio-P) increased firstly and then decreased with increasing P application rate (0∼9.6 g P₂O₅⋅m^–2), which verified our hypothesis. The study in Loess Plateau pointed out that NaHCO₃-Pi, NaOH-Pi, and HCl-Pi increased with P application rate (0∼9.2 g P₂O₅⋅m^–2), which was inconsistent with our result (Mahmood et al., 2020). The decreasing soil pH caused by P application could dissolve Ca-P in calcareous soils of the Loess Plateau, thus resulting in increasing soil Pi. However, most added Pi was absorbed and precipitated with Fe, Al, and Ca in acid soils, and the lower soil pH in high-P treatments would promote the absorption and precipitation (Wang et al., 2014; Maranguit et al., 2017). The increasing trend of NaOH-Po with P application rate may owe to the decrease in soil pH (Figures 1E, 6), because NaOH-Po is often associated with humic and fulvic acids (Wang et al., 2006). Moreover, NaOH-Po also contains the Po absorbed on Al and Fe hydrous oxides (Moro et al., 2015), and the increasing soil NaOH-Po is coupled well with Al and Fe content in current study (Figure 6). Soil Residual-P increased firstly and after that it had an approximate constant with the increase of P application rate, suggesting that soil P coated by Fe/Al hydrous oxides tends to be saturated when P application rate > 4.8 g P₂O₅⋅m^–2.

The hydrolysis of organic P is predominantly mediated by phosphatase, especially acid and alkaline phosphatases (Fraser et al., 2015b). The phoD alkaline phosphatase is solely of microbial origin and widespread in soil (Ragot et al., 2015). Agricultural practices such as fertilization could affect the soil properties and thus influence the phoD-harboring microbial community (Fraser et al., 2015b; Liu et al., 2020). Previous studies have indicated that P fertilization could affect specific taxa related to the phoD-harboring microbial community (Liu et al., 2020; Wan et al., 2020; Chen et al., 2021; Wang M. M. et al., 2022). For example, the study in Loess Plateau found that the relative abundance of Proteobacteria group in high-P treatments (9.2 g P₂O₅⋅m^–2) was higher than that in low-P treatments (0∼6.9 g P₂O₅⋅m^–2) (Liu et al., 2020). Similarly, we found that the relative abundance of Proteobacteria in high-P (9.6 g P₂O₅⋅m^–2) treatments was higher than that in low-P (0∼4.8 g P₂O₅⋅m^–2) treatments (Figure 2I). Proteobacteria was reported to have copiotrophic lifestyle and thrive better under higher nutrient conditions (Veraart et al., 2015; Li et al., 2020), and the higher abundance of Proteobacteria was accordant with the higher soil TC in our study. Though the relative abundance of Actinobacteria varied slightly among the seven treatments (Figure 2B), we found four of six biomarkers in P30 derived from Actinobacteria (Figure 3), and soil TC in P30 was the lowest (Figure 1K). This was consistent with the oligotrophic lifestyle of Actinobacteria (Ho et al., 2017). The Deinococcus_Thermus group could resist the ultraviolet light, and the lower abundance of Deinococcus_Thermus (Figure 1D) in high-P treatments might be owed to the higher plant coverage (Griffiths and Gupta, 2007).

The study sampled from a maize field reported an increasing trend of phoD-harboring microbial diversity with P application rate (Chen et al., 2019). Another study in pasture soils also found PhoD gene diversity increased with P application rate (Tan et al., 2013). In the present study, we found that phoD-harboring microbial diversity also increased with increasing P application rate when P application rate < 4.8 g P₂O₅⋅m^–2 (Figure 4), which disagreed with our second hypothesis. The increasing trend of phoD gene diversity indicated that a small amount of P input could meet the P demand of phoD-harboring microbes and further promote the growth of phoD-harboring microbes (Chen et al., 2019). However, the decreasing trend of phoD-harboring microbial diversity for P application rate > 4.8 g P₂O₅⋅m^–2 (Figure 4) suggested that excessive P inputs might suppress the growth and reproduction of phoD-harboring microbes (Wei et al., 2019). Pimm and Pimm (1982) pointed out that sufficient resource could help to improve the propagation of other competitive species, thus phoD-harboring microbes might lose the advantage to outcompete other competitive microbes after high-P application. The opposite relationship between phoD-harboring microbial diversity and P application amount for below and above P application amount = 4.8 g P₂O₅⋅m^–2 suggested a P fertilization threshold for optimistic growth of phoD-harboring microbes in the alpine grassland on the Qinghai-Tibetan Plateau. Similarly, Liu et al. (2020) conducted a P fertilization experiment (0∼9.2 g P₂O₅⋅m^–2) in Loess Plateau and they also found a P fertilization threshold (6.9 g P₂O₅⋅m^–2) for optimistic growth of soil microbes. Those P fertilization thresholds in the Qinghai-Tibetan Plateau and Loess Plateau further proved the “intermediate disturbance hypothesis” that local species diversity to be maximal at an intermediate level of disturbance (Bongers et al., 2009). Considering the higher percentage of the dominant phyla (e.g., Actinobacteria, Cyanobacteria, and Proteobacteria) in P60 treatment (Figures 2B,C,I), the lower phoD-harboring microbial diversity in P60 treatment likely depended on the growth of dominant taxa, but not rare ones.

The soil phoD-harboring microbes interact with each other and form a connected network to adapt to their habitat. In the current study, the interactions between ASVs (nodes and edges), degrees, connectance, and transitivity were higher in the high-P network (P40, P50, and P60, Figure 5B) than that in low-P network (P0, P10, P20, and P30, Figure 5A), suggesting a larger and more complex network in high-P treatments. The complexity of a network is one of the important traits of ecosystem stability and functioning (Wei et al., 2019). The higher network complexity in the high-P group benefits the stability of the phoD-harboring microbial community, suggesting that phoD-harboring microbes in the high-P group could mediate internal interactions to resist to stresses. Those microbes driving the same biochemical process can compete for similar resources and resulting in conflicting interactions (Le Gac et al., 2012). The more conflicting interactions in the phoD-harboring microbial network under low-P fertilized conditions (percent of the negative edges: 48.23%) than under high-P fertilized conditions (percent of the negative edges: 34.63%) indicated a stronger competition within the phoD-harboring microbes in low-P fertilized condition, and P fertilization could alleviate this competition. Furthermore, the generalists in the phoD-harboring microbial communities in the low-P (Figure 5C) and high-P (Figure 5D) networks were different. In the low-P network, the generalist Rhizobiales is a nitrogen-fixing taxon, but it also can solubilize phosphate (Halder et al., 1990), suggesting that Rhizobiales could increase soil available P in the P deficient condition. Apart from Rhizobiales, Burkholderiales, and Streptomycetales were the generalists in the high-P network. Burkholderiales is well-known as pathogens (Banerjee et al., 2018) while Streptomycetales could excrete antibiotics (Harir et al., 2018). Different functions of the generalists suggested a versatile phoD-harboring microbial community in high-P fertilized conditions.

Phosphorus fertilization can affect the soil phoD-harboring community by changing soil physiochemical properties (Ragot et al., 2016; Hu et al., 2018; Chen et al., 2019; Wei et al., 2019; Liu et al., 2020; Wan et al., 2020). In the present study, P application affected soil phoD-harboring community structure through Bio-P, TC, Fe, NaOH-Po, and pH. The study in a maize field reported that the shift in PhoD gene community structure was related to P availability (Chen et al., 2019). Similarly, soil Bio-P held the largest positive effect on soil phoD-harboring microbial community structure in the present study (Figure 7B). Bio-P is widely known to be a key factor affecting soil phoD-harboring microbes, as P is the essential component of soil phoD-harboring microbial nutrition (Chen et al., 2019; Liu et al., 2020; Wan et al., 2020; Wei et al., 2021; Wang M. M. et al., 2022). The higher Bio-P could stimulate the growth of soil phoD-harboring microorganisms, which could enhance the mineralization of organic P, releasing Pi and supplementing Bio-P in return (Chen et al., 2019). However, the positive feedback could not last forever. Another study reported a negative relationship between P availability and phoD microbial diversity in wheat rhizosphere soil (Xie et al., 2020). This discrepancy may be explained by the different age of fertilization. The P application in the wheat field lasted 7 years, which was longer than that in our study (2 years), and the long-term P application could suppress the reproduction of the phoD microbes (Wang M. M. et al., 2022; Zeng et al., 2022). Previous studies reported the significant effect of soil carbon on phoD-harboring microbial community structure (Ragot et al., 2016; Luo et al., 2017; Hu et al., 2018; Liu et al., 2020; Wang M. M. et al., 2022). Carbon as an energy supply for microbes was demonstrated to promote the propagation of microbial community, and the phoD-harboring microbial respiration could lead to carbon loss (Luo et al., 2017). The lowest soil TC, the highest phoD-harboring microbial diversity and the highest Bio-P in P30 suggested that phoD-harboring microbes likely depleted carbon resource to release Pi. Based on the SEM, the total effect size of Bio-P (0.60) on PhoD gene community structure was higher than that of TC (−0.55), reflecting the more important role of Bio-P. Moreover, the different effects of Bio-P and TC on phoD-harboring microbial community structure suggested that P mineralization is decoupled with C mineralization in the alpine grassland of the Qinghai-Tibetan Plateau.

In addition, the effect of soil Bio-P on soil phoD-harboring community could not be distinguished from that of Fe, given that Fe is commonly correlated with soil available P in acid soils (Wang et al., 2006). Soil Fe associated P acts as a buffer in P supply when the soil P is sufficient, and it could be released into soil solution and transformed into Bio-P when the soil Bio-P is insufficient (Maranguit et al., 2017). Another study in Karst soils also reported that metal ions affected soil phoD microbes through available P (Hu et al., 2018). The negative relationship between NaOH-Po and soil phoD gene (Figure 7) was consistent with the result in subtropical Chinese fir plantations (Wang et al., 2021). NaOH-Po is a kind of organic P that can be slowly mineralized into inorganic P when the soil available P is inadequate (Cao et al., 2022). Conversely, the inorganic P could be transformed into organic P by plant and soil microbial assimilation (Giles et al., 2014). The microbial cell is a kind of NaOH-Po source (Gu et al., 2017), and the high MBP under high-P treatment (Figure 1) was conductive to NaOH-Po accumulation. Moreover, soil phoD-harboring microbes compete with other microbes for resources, the low phoD-harboring microbial diversity and high MBP suggested a strong competition between phoD-harboring microbes and other microbes under high-P conditions. Similarly, P fertilization could promote the growth of the plant (Mahmood et al., 2020), which may enhance the P demand of plants and then intensify the competition between plants and phoD microbes for P. The competitions between soil phoD-harboring microbes, other soil microbes, and plants might be the reason why NaOH-Po negatively correlated with soil phoD microbial community structure (Wang et al., 2021).

Fertilization has been reported to decrease soil pH, and soil pH has also been shown to drive the soil microbial community (Liu et al., 2018; Chen et al., 2019; Xu et al., 2021). Mantel test indicated that soil pH was significantly correlated with phoD-harboring microbial community structure (Figure 6), and this was also confirmed by the negative effect of soil pH on phoD microbial community structure in the SEM model (Figure 7), suggesting that phoD microbial communities are less similar in soils with greater difference in pH. Soil pH has also been found as the important factor affecting phoD microbial communities in soil under grass (Ragot et al., 2015, 2016), Luvisol under maize (Chen et al., 2019), and Karst soil under soybean (Hu et al., 2018). The decreased soil pH enhanced the NaOH-Po content, as NaOH-Po is associated with humic and fulvic acids (Cao et al., 2022), further led to the change in phoD-harboring microbial communities as discussed above.