Multiple environmental stressors interactively affect soil phosphorus cycling microbiomes

Abstract

Soil phosphorus (P) cycling microorganisms (PCMs) facilitate organic P mineralization (represented by the phoD gene) and inorganic P dissolution (pqqC), driving terrestrial P cycling and enhancing P availability to plants. However, the interactive effects (i.e., additive, synergistic, and antagonistic interactions) of multiple environmental stressors (i.e., acidification, salinization, drought, metal pollution and human disturbance) on PCMs remain unclear. Combining a large-scale survey and controlled experiments, we found that the interactive effects of stressors with significance on PCM abundance and diversity were more ubiquitous than the effects of individual stressors. Specifically, acidification and drought had antagonistic and synergistic effects on the abundances of organic P mineralization microorganisms and inorganic P dissolution microorganisms, respectively, while antagonistic effects between drought and metal pollution were observed on the diversity of these microorganisms. The PCM changes were associated with soil P-related stoichiometric parameters under interactive effects, indicating that interactive stressors strongly influenced microbial-regulated P cycling and P availability in soil. These findings underscore the need to incorporate interactive stressor effects into strategies for sustainable soil P management and ecosystem resilience.

Interactive stressors influence microbial-regulated phosphorus cycling and availability in soil by changing the abundance and diversity of phosphorus cycling microorganisms and associated stoichiometric parameters, according to an analysis of soil samples from a large-scale survey across China

Introduction

Phosphorus (P) is a critical nutrient in ecosystems, essential for plant growth and ecosystem productivity¹. However, its availability is often limited due to strong fixation in soils. Phosphorus cycling microorganisms (PCMs) are essential for regulating P biogeochemical cycling and making phosphorus available to plants. Microorganisms containing phoD, pqqC, phnK, bpp, phoX, cphy and gcd genes increase soil P availability by mineralizing organic P (phytate, phospholipid, etc.) or dissolving inorganic phosphates, such as Fe-P, Al-P, etc^{2, 3–4}. In particular, the phoD gene, commonly found in phosphorus-cycling microorganisms, encodes alkaline phosphatase that hydrolyzes organic phosphorus (OP) compounds into inorganic phosphate, thereby making phosphorus accessible to plants^5,6. The pqqC gene encodes pyrroloquinoline quinone (PQQ) synthase, which governs the oxidation pathway of glucose to gluconic acid, facilitating mineral phosphate solubilization⁷.

The abundance and diversity of soil P cycling microorganisms vary across ecosystems^{8, 9–10}, and are influenced by edaphic, climatic factors and human disturbances^{11, 12, 13–14}. A higher abundance of P cycling microorganisms generally leads to greater overall P solubilization potential, and higher diversity provides functional redundancy and complementary mechanisms for P solubilization, ensuring resilience and adaptability to environmental changes¹⁵. Studies have shown that drought increases while soil acidification decreases the abundance of the gcd and phoD genes^4,16. Furthermore, acidification, salinization, and metal pollution significantly change the P-cycling microbial diversity, such as community composition^{9,17, 18–19}. However, these studies primarily address the individual effects of environmental factors, overlooking their interactive impacts on P-cycling microorganisms. As global climate change and human activities increasingly cause those environmental stressors (e.g., soil acidification, metal pollution, salinization, drought and human disturbances) to co-occur in soil ecosystems^20,21, their effects on P-cycling microorganisms are likely driven by complex interactions among those stressors²². Moreover, past studies failed to provide quantitative data on the extent or direction of the interactions of different stressors on P-cycling microorganisms^18,23,24.

Predicting the combined effects of environmental stressors based on the additive effects of single factors may overlook the complexity of interactive impacts, such as synergistic and antagonistic interactions, on phosphorus-cycling microorganisms (PCMs). Synergistic effects can lead to more pronounced changes in microbial communities than predicted by single-stressor studies, while antagonistic effects can mitigate the overall impact²⁵. Regarding interactive effects on PCM communities, three critical gaps remain: (1) interactive effects need to be quantified to identify dominant environmental stressors, as these may induce more significant impacts on PCM communities than individual stressors²⁶, 2) the direction of interactive effects (synergistic or antagonistic) on soil PCMs need to be clarified to determine whether stressor combinations counteract or amplify each other; and 3) the extent to which changes in interactive effects on PCMs are closely associated with soil P-related parameters remains unexplored in previous studies.

Here, we collected soil samples from 180 sites across diverse ecosystems encompassing various climatic and edaphic conditions. These sites make up environmental gradients for soil pH, salinity, aridity, metal pollution, and human disturbance intensity. Samples with various environmental stressors were categorized into no-stressed and stressed groups based on established thresholds²¹. Then, we investigated the responses of the abundance and diversity of soil organic P mineralization and inorganic P dissolution microorganisms to single stressors and their combinations. We also conducted controlled experiments to verify the results from the large-scale field study. The primary objectives were: 1) to evaluate the individual and interactive effects of environmental stressors on the abundance and diversity of soil P cycling microorganisms, and 2) to find the relationship between changes in PCM abundance and diversity with soil P-related stoichiometric parameters (e.g., total C/total P, total C/available P, etc.). We hypothesized that, in addition to individual effects, the interactive effects of environmental stressors play an important role in regulating soil P cycling microorganisms, which further contributed to changes in soil P-related stoichiometric parameters. Our study aims to improve the prediction and management of microbially-regulated soil P cycling under intensified global change and human activities.

Results

Overview of environmental stressors

Of the 180 soil samples collected, the threshold values at 50% of pH value, aridity, nemerow pollution index, human influence index and sodium adsorption ratio, which represented soil acidification, drought, metal pollution, human disturbances, and salinization were 6.81, 0.82, 1.85, 26, and 1.19, respectively (Fig. 1a and Table S1). The samples with the values (i.e., aridity, nemerow pollution index, human influence index and sodium adsorption ratio) higher than 50% thresholds were categorized into stress groups, except for the pH that had the opposite criterion. Variance partitioning analysis revealed that acidification accounted for the largest variance in the abundances of P cycling microorganisms (62.3%), organic P mineralizing microorganisms (57.9%), and inorganic P dissolving microorganisms (77.8%) (Fig. 1b). Drought also had a significant impact, explaining 21.2% of the variance for P cycling microorganisms, 23.9% for organic P mineralizing microorganisms, and 11.8% for inorganic P dissolving microorganisms. By comparison, salinization, human disturbance, and metal pollution contributed less to the abundance of P-cycling microorganisms than acidification and drought (Fig. 1b).

Fig. 1 The classification of environmental stresses and the effect of environmental stress on P cycling microorganisms. [Images not available. See PDF.]

a Frequency distribution of soil pH, aridity, sodium adsorption ratio (SAR), integrated Nemerow pollution index (NPI), and human influence index (HII) among the studied soil samples. b Contribution of environmental stressors to the abundances of P cycling microorganisms using a variance partitioning scheme. The pie charts indicate the percentage of explained variance for each stressor associated with the abundance of soil microorganisms. Total microbe: all microorganisms in the soil; Total PCM: total P cycling microorganisms; Mineralization: organic P mineralization by microorganisms; Dissolution: inorganic P dissolution by microorganisms.

Individual and interactive effects on the PCM abundance

The abundance of P cycling microorganisms was influenced by both individual and interactive effects of environmental stressors (Fig. 2). Acidification and drought had the greatest number of individual effects on gene abundance, followed by salinization, with human disturbance and metal pollution showing the least impact (Fig. 2a). Drought exhibited positive individual effects on gene abundances. In contrast, acidification had negative effects (Fig. 2b).

Fig. 2 Individual and interactive effects of environmental stressors on the abundances of genes involved in microbial P cycling. [Images not available. See PDF.]

a Number of individual effects of environmental stressors on the abundances of the genes involved in microbial P cycling. b Boxplots of individual effects of acidification and drought on the abundances of genes involved in microbial P cycling. c Number of interactive effects of environmental stressors on the abundances of genes involved in microbial P cycling. d Interactive effects of acidification*drought on the abundances of genes involved in organic P mineralization (phoD, phnK, bpp, phoX and cphy) and inorganic P dissolution (gcd and pqqC). The error bars in the dots represent the 95% confidence interval (CI). If the 95% CI overlapped with zero, the interactive effect was considered additive. If the 95% CI did not overlap with zero, the synergistic or antagonistic interactive effects were suggested. The criteria for defining synergistic or antagonistic interactive effects are outlined in the Methods section. Total PCM: total P cycling microorganisms; Mineralization: organic P mineralization by microorganisms; Dissolution: inorganic P dissolution by microorganisms. The genes involved in microbial P cycling are listed in Table S1.

The combination of acidification and drought produced the highest number of interactive effects (i.e., 8) on gene abundance. In contrast, interactive effects from other stressor combinations were predominantly additive (Fig. 2c). Antagonistic interactions were observed between drought and acidification for the total abundance of all PCMs (effect size of 0.54), as well as for those involved in organic P mineralization (0.53) and inorganic P dissolution (0.55). Specifically, the interactive effects of acidification*drought were antagonistic for abundances of phoD (0.63), bpp (0.27), and gcd (0.58), but synergistic for phnK (0.62) and pqqC (0.47) genes (Fig. 2d).

Individual and interactive effects on PCM diversity

Alpha diversity and community composition of PCM harboring phoD and pqqC genes were influenced by both individual and interactive effects of environmental stressors (Fig. 3). Acidification and drought had individual effects on the alpha diversity, community dissimilarity and genera abundances of PCMs containing phoD and pqqC genes (Fig. 3a). Metal pollution also affected the community dissimilarity of PCMs containing phoD and pqqC genes.

Fig. 3 Individual and interactive effects of environmental stressors on the diversity and community composition of P cycling microorganisms. [Images not available. See PDF.]

a The number of individual effects of environmental stressors on alpha diversity (Chao 1 index), community dissimilarity and the abundances of P cycling microorganisms at the genus level. b The number of interactive effects of environmental stressors on alpha diversity (Chao 1 index), community dissimilarity and the abundances of P cycling microorganisms at the genus level. c Interactive effects of drought*metal pollution on alpha diversity (Chao 1 index) and community dissimilarity of P cycling microorganisms. d Interactive effects of drought*acidification and drought*metal pollution on the abundances of P cycling microorganisms containing pqqC and phoD genes at the genus level. The error bars in the dots represent the 95% confidence interval (CI). If the 95% CI overlapped with zero, the interactive effect was considered additive. If the 95% CI did not overlap with zero, the synergistic or antagonistic interactive effects were suggested. The criteria for defining synergistic or antagonistic interactive effects are outlined in the Methods section.

The combination of drought and metal pollution had antagonistic interactive effects on the community dissimilarity of PCMs containing phoD and pqqC genes (Fig. 3b, c). Synergistic effects between drought and acidification were found for phoD-harboring Mesorhizobium and Stenotrophomonas, and pqqC-harboring Mycolicibacterium. In contrast, antagonistic interactions of drought*acidification affected phoD-harboring Streptomyces, while the combination of drought and metal pollution antagonistically affected pqqC-harboring Methylobacterium (Fig. 4d).

Fig. 4 Relationship between P cycling microorganisms and soil P-related parameters caused by individual and interactive effects of stressors. [Images not available. See PDF.]

a Spearman correlations between soil P-related parameters and the abundance, alpha diversity and community dissimilarity of pqqC and phoD genes caused by individual effects and interactive effects of stressors. b The explanation of soil P-related parameters in the variances of the abundances of genes involved in P cycling caused by individual and interactive effects of stressors revealed by random forest analysis.

Relationships between PCMs and soil P-related stoichiometric parameters

The abundance, alpha diversity, and community dissimilarity of phoD and pqqC-harboring communities were not only correlated with soil P-related stoichiometric parameters under individual effects but also with those under interactive effects. The interactive effects had contrasting effects on soil P-related stoichiometric parameters compared with individual effects (Fig. 4). Specifically, the soil TC/TP ratio was positively correlated with the abundances of both phoD (p < 0.05, r = 0.66) and pqqC (p < 0.05, r = 0.73) microorganisms, and the soil TN/AP ratio was positively correlated with community dissimilarity for both phoD (p < 0.05, r = 0.88) and pqqC (p < 0.05, r = 0.83) microorganisms under interactive effects (Fig. 4a). However, the soil TC/TP and TN/AP showed no correlations (p > 0.05) with the corresponding gene abundances under individual effects.

Random forest models indicated that the individual effects of stressors on soil TC/TN, TC/TP, TN/TP, TC/AP, and TN/AP ratios were significant predictors of the abundance of microorganisms containing phoD and pqqC genes (p < 0.05, Fig. 4b). Under interactive effects, the soil TC/TP ratio was the strongest predictor of abundances of phoD and pqqC-containing microorganisms. For alpha diversity, soil TC/TN, TN/TP, TC/AP, and TN/AP under individual effects were significant predictors, whereas no parameters under interactive effects significantly predicted alpha diversity except the TC/TP and TN/TP for pqqC community (Fig. 4b). Soil TN/AP ratio under both individual and interactive effects was the strong predictor of community dissimilarity for phoD- and pqqC-harboring communities (Fig. 4b).

Results from controlled experiments

Controlled experiments showed trends consistent with results from the large-scale study, irrespective of soil type (Fig. 5). For both soil types, acidification stress reduced the abundances of the phoD and pqqC genes, while drought increased their abundances. The interactive effects of acidification and drought were antagonistic for phoD gene abundance but synergistic for pqqC gene abundance (Fig. 5a, b). Drought decreased the community dissimilarity of microorganisms containing the phoD and pqqC genes in the Red soil but did not affect them in the Brown soil. Metal pollution increased the community dissimilarity of microorganisms containing the phoD gene only. Consistent with the survey results, the interactive effects of metal pollution and drought were synergistic for microorganisms containing both the phoD and pqqC genes (Fig. 5c, d).

Fig. 5 Changes in the abundance and community dissimilarity of microorganisms containing phoD and pqqC genes under the interaction of stressors for two typical soils by verification experiments. [Images not available. See PDF.]

a Abundance of the phoD gene in the treatments of control, drought, acidification and drought + acidification for two soils, and their interactive effects. b Abundance of the pqqC gene in the treatments of control, drought, acidification and drought + acidification for two soils, and their interactive effects. c Community dissimilarity for microorganisms containing the phoD gene in the treatments of control, drought, metal and drought+metal, and their interactive effects. d Community dissimilarity for microorganisms containing the pqqC gene in the treatments of control, drought, metal and drought+metal, and their interactive effects. The error bars in the dots represent the 95% confidence interval (CI). If the 95% CI overlapped with zero, the interactive effect was considered additive. If the 95% CI did not overlap with zero, the synergistic or antagonistic interactive effects were suggested. The criteria for defining synergistic or antagonistic interactive effects are outlined in the Methods section.

Discussion

The PCMs, which typically harbor phoD and pqqC genes involved in organic P mineralization and inorganic P dissolution, enhance P availability and drive P biogeochemical cycling⁴. In general, microorganisms with higher abundance directly increase P dissolution rates by organic acid production and enzyme activity. Microorganisms with higher diversity support multiple P-solubilization strategies, such as mineral dissolution and enzymatic hydrolysis, and ensure resilience against environmental stressors^27,28. While existing research primarily focuses on the effects of individual stressors on PCMs^29,30 or qualitatively describes the interactive effects^18,23,24, this approach overlooks the complex interactions between co-occurring environmental stressors in real-world conditions.

Based on both large-scale survey and controlled experiments, we found that most combinations of environmental stressors exhibited primarily additive effects on PCMs (Table S2), indicating that, in many cases, the impacts of multiple stressors can be estimated from the additive effects of single factors. However, the soil acidification + drought demonstrated significant interactive effects on PCM abundance, which indicates the capacity of microorganisms to release the fixed P in soil⁴, with greater interactive impacts on gene abundances than observed with additive effects of the individual stressors. This suggests that the interactive effects of environmental stressors are potent and must be accounted for to accurately assess microbial P solubilization capacity. Similarly, drought, acidification, and metal pollution showed interactive effects on PCM diversity, which indicates the complementary functional traits of P solubilization across microbial species^15,31, underscoring the importance of interactive effects of multiple stressors on P solubilization efficiency. We acknowledge that controlled experimental setups, such as adding soluble heavy metals, may artificially amplify metal toxicity to microorganisms compared to natural conditions. Nevertheless, the convergence of these experimental results with large-scale field survey data validates our core conclusions. These findings emphasize that studying the abundance and diversity of PCMs solely under single-stressor conditions may miss critical interactions that influence soil P availability and biogeochemical cycling.

Shifts in microbial P solubilization capacity are closely tied to soil P availability³², as also evidenced by positive correlations between gene abundance/diversity and P stoichiometric ratios (Figs. 4a and 4b). This suggests that shifts in microbial P solubilization capacity are closely tied to soil P availability and soil P cycling under combined stress conditions. Notably, correlations driven by interactive effects differed from those observed with individual effects, indicating that interactive effects, often overlooked, provide additional explanatory power for soil P-related parameter changes beyond what is accounted for by individual stressors alone. Our findings suggest that multiple stressors can amplify or mitigate microbial responses in complex ways, particularly regarding the community composition and function of microbial P cycling.

By quantifying the individual and interactive effects of environmental stressors, we identified acidification and drought as dominant factors significantly influencing the abundance of PCMs. The abundance of PCMs reflects the soil’s potential P-solubilizing capacity to release fixed P. As a stressor, acidification reduces PCM abundance (Fig. 2). Previous research indicates that low soil pH imposes physiological constraints on microorganisms³³, and the PCMs generally perform more effectively at higher pH levels³⁴. Thus, soil acidification may directly hinder the growth of PCMs. Additionally, moderate acidification increases P availability (Table S3), which could suppress the growth of phoD- and pqqC-harboring microorganisms, as these P-cycling microorganisms tend to remain dormant when P is sufficiently available³⁵. In contrast, drought increased the abundance of PCMs (Fig. 2). Typically, microorganisms release more hydrolytic phosphatases (for mineralization) and organic anions (for dissolution) when P levels are low³⁶. Drought reduces nutrient diffusivity and mass flow³⁷, slowing the conversion of recalcitrant P into more accessible forms and leading to P deficiency¹⁶, consistent with decreased soil available P under drought in our study (Table S3). This P scarcity likely stimulates P solubilization³⁸. While some studies suggest that reduced moisture availability can hinder microbial growth³⁹, the increased PCM abundance under the drought observed in this study may be due to 1) a moderate, rather than extreme, drought under which certain PCMs can thrive⁴⁰, and 2) the drought-induced P deficiency that could stimulate the P solubilizing activity of PCM.

The antagonistic effect of acidification and drought on phoD-harboring microorganisms (involved in organic P mineralization) and their synergistic effect on pqqC-harboring microorganisms (involved in inorganic P dissolution) suggest that the interactive effects are complex and non-additive, impacting microbial P solubilization capacity in distinct ways. The antagonistic effect on phoD microorganisms indicates that although acidification and drought each stress these microorganisms individually, their combined effect may alleviate the negative pressures on microbial potential of organic P mineralization. Conversely, the synergistic effect on pqqC microorganisms suggests that regions suffering from acidification and drought together would accelerate the suppression of microbial inorganic P dissolution, experiencing an extreme P deficiency. Those contrasting interactive effects on phoD- and pqqC-harboring microorganisms may stem from their different metabolic pathways⁴¹: phoD-harboring microorganisms produce extracellular enzymes to mineralize organic P, while pqqC-harboring microorganisms dissolve inorganic P through H⁺ and organic acid production^2,3. The antagonistic interaction for phoD microorganisms likely arises because acidification’s P-release partially counteracts drought’s P-deficiency stimulus, and the direct pH inhibition remains dominant. The synergistic suppression of pqqC microorganisms occurs as drought strongly stimulates microorganisms, but the combination of inherent soil acidification with their massive output of H⁺ and organic acids creates an extremely acidic environment. This combined acidity surpasses microbial tolerance thresholds, resulting in greater suppression than either stressor applied alone.

Here, we demonstrate that PCM abundance and diversity undergo dramatic changes only when soil pH shifts from alkaline (pH > 7.0) to acidic (pH < 7.0). Furthermore, at the lower pH thresholds tested (6.5, 6.0, and 5.5), which fall entirely within the acidic range, the changes in PCMs are less pronounced than those observed at the threshold of 6.81 derived from the 50% of maximum value criterion (Fig. S1). This strongly suggests that the transition from alkaline to acidic conditions represents a critical turning point for the PCM community shift during soil acidification.

Our study found that drought and metal pollution are dominant stressors influencing the diversity of PCMs. Drought was observed to reduce PCM community heterogeneity, leading to community homogenization due to the survival of a limited number of drought-tolerant species^16,39,42. This effect likely arises because only a subset of highly drought-tolerant PCMs can survive and thrive under restricted moisture conditions, thereby reducing overall community diversity⁴³. Conversely, metal pollution was found to increase community heterogeneity (Fig. 3c), suggesting that the microbial community structure shifted towards more diverse compositions under metal stress. This increase in diversity may be attributed to different microorganisms employing various strategies to tolerate or detoxify metals, such as sequestration, efflux pumps, or metabolic modifications, which support a wider range of species coexisting in metal-stressed environments⁴⁴. Furthermore, correlation analyses between total and bioavailable concentrations of heavy metals revealed strong, significantly positive relationships, observed for both raw total concentrations and normalized data (Fig. S2). This suggests that the biological stress exerted by bioavailable heavy metals on PCM diversity probably corresponds to that of total metal concentrations.

Microbial diversity, including alpha diversity and community dissimilarity, is indicative of the ecosystem’s capacity to release P from complex P compounds, as higher diversity allows a broader array of microorganisms to contribute to different P releasing pathways⁴⁵. The interactions between drought and metal pollution had antagonistic effects on the community dissimilarity of phoD- and pqqC-harboring microorganisms, suggesting that changes in microbial community composition occur in counterbalancing ways. This antagonistic interaction may help preserve functional diversity even as species richness declines, thus maintaining the soil’s capacity for P solubilization^29,46,47. This outcome is likely due to the counteraction between metal-induced microhabitat heterogeneity and drought-driven homogenization: drought primarily favors drought-tolerant microorganisms, reducing diversity due to limited resources and habitat homogeneity, whereas metal pollution favors metal-tolerant species employing a range of tolerance mechanisms. Their opposing effects may counterbalance dominance by any single tolerance strategy, thereby preserving microbial functional diversity. Additionally, combined stressors, such as drought and metal pollution, or drought and acidification, affect the abundance of specific genera, including Mesorhizobium, Stenotrophomonas, and Mycolicibacterium. These taxa appear to be the primary contributors to shifts in community composition and are especially responsive to combined stress. For instance, Mesorhizobium is known to benefit leguminous plants, and Stenotrophomonas has been associated with plant growth promotion and pathogen suppression^48,49. Therefore, the interactive effects on the abundances of these genera could extend beyond microbial dynamics to potentially influence plant growth and health.

Conclusions

Understanding the individual and interactive effects of environmental stressors on the abundance and diversity of soil PCMs provides valuable and comprehensive insights into the regulation of soil P availability and cycling in a rapidly changing world. Acidification, drought, and metal pollution affected the abundance and diversity of PCMs compared to salinization and human disturbances by causing interactive effects (Fig. 6), underlining the importance of managing soil acidity, aridity and preventing pollution to restore P solubilizing capacity if they occurred simultaneously. Future studies should investigate how shifts in the soil PCM community caused by the interactive effects of stressors impact aboveground productivity, especially given that P is often a limiting nutrient for plant growth. Overall, this study highlights the underappreciated role of interactive environmental stressors in regulating soil microbial functions essential for ecosystem nutrient cycling.

Fig. 6 Conceptual diagram illustrating the individual and interactive effects of five environmental stressors on the abundance and community of microorganisms involved in organic P mineralization (phoD gene) and inorganic P dissolution (pqqC gene). [Images not available. See PDF.]

OP: organic P; IP: inorganic P; AP: available P.

Methods and Materials

Site description and sampling

We collected soil samples from 180 sites across China (Fig. S3), with those sites representing a broad range of environmental conditions of precipitation, temperature, soil type and human influence. At each site, we collected five soil cores (0–20 cm depth), which were then combined to form a composite sample. After sampling, the composite samples were transported on ice to the laboratory. Each sample was homogenized and divided into two subsamples: one was air-dried for soil chemical analysis, and the other was frozen at −80 °C for DNA extraction.

Analysis of soil chemical properties

Soil pH was measured with a pH meter in a soil-to-water suspension of 1:2.5. Soil total carbon (C) and nitrogen (N) were determined by an elemental analyzer (Vario EL II, Germany). Total phosphorus (TP) was heat-digested using concentrated H₂SO₄-HClO₄, and available phosphorus (AP) was extracted from the soil using 0.5 mol/L NaHCO₃ (adjusted to pH 8.5) at a soil-to-solution ratio of 1:20, the extracted solution was then filtered through Whatman Grade No. 42 quantitative filter paper. Soil samples (200 mg) were extracted with 1 mol/L HCl (20 ml, 16 h), and the supernatant was separated from the residue by centrifugation (20 min, 4530 rpm) to assess the concentration of inorganic phosphorus (IP). The residue was calcined for 3 h at 450 °C, and the ash was extracted for OP using 1 mol/L HCl⁵⁰. The P concentration in all extracts was measured by the molybdenum-blue colorimetry method. To measure the total heavy metal (Cd, As, Cr, Cu, Mn, Ni, Pb, and Zn) concentrations, soil samples were digested by the HF-HNO₃-H₂O₂ (1:2:1, v/v/v) at 210 °C for 60 min using a microwave. Soil available heavy metals were extracted using a solution containing 0.005 mol/L DTPA, 0.01 mol/L CaCl₂, and 0.1 mol/L triethanolamine (pH 7.30). The extracts were then filtered through Whatman Grade 42 filter paper⁵¹. Ammonium acetate (1 mol/L, pH=7.0) was used as an ion-exchange extractant to extract soil-exchangeable Na, Ca, and Mg. Finally, inductively coupled plasma optical emission spectrometry (ICP-OES) (iCCP 6300, Thermo Scientific, USA) was used to determine the concentrations of soil heavy metals and exchangeable bases.

DNA extraction and gene quantification

Soil microbial DNA was extracted using the FastDNA SPIN kit (MP Biomedicals, Solon, OH, USA) according to the manufacturer’s protocol. Absolute abundances of genes involved in organic P mineralization (e.g., bpp, phnK, cphy, phoX, phoD) and inorganic P dissolution (gcd and pqqC) were quantified using a high-throughput qPCR-based chip (QMEC) on the SmartChip Real-Time PCR System (WaferGen Biosystems, Fremont, USA). Briefly, the qPCR amplification protocol began with a denaturation step at 95 °C for 10 min, followed by 40 cycles, each consisting of 30 s of denaturation at 95 °C, 30 s of annealing at 58 °C, and 30 s of extension at 72 °C. Each DNA sample was amplified in triplicate. Results were excluded from analysis if amplification efficiencies were below 1.8 or above 2.2, if amplification occurred in the negative control, or if the threshold cycle (CT) value exceeded 31⁵². The specific functions of genes involved in P cycling are detailed in Fig. S4 and Table S4.

Sequencing and diversity estimates

Soil DNA samples were then sent for the specific phoD and pqqC gene (key functional genes involved in organic P mineralization and inorganic P dissolution, respectively) sequencing on Illumina NovaSeq 6000 (Illumina Inc., San Diego, CA). The phoD gene primers for organic P mineralization were phoD-F733 (TGGGAYGATCAYGARGT) and phoD-R1083 (CTGSGCSAKSACRTTCCA)¹¹. The pqqC gene primers for inorganic P dissolution were: pqqC2-F (AACCGCTTCTACTACCAG) and pqqC2-R (GCGAACAGCTCGGTCAG)³⁴. A Qubit® 2.0 Fluorometer was used to monitor the concentration and quality of DNA samples (DNA concentration > 10 ng/µL).

Raw data quality was ensured using Fastp v. 0.20.0 with sliding window parameters (-W 4, -M 20), and primers were eliminated using Cutadapt⁵³. Paired-end reads were merged, and a feature table was generated using the Uparse, with chimeric reads trimmed⁵⁴. Taxonomy assignment to sequences was performed using FunGene v1.0 database at the confidence threshold of 0.8. Alpha diversity, i.e., Chao 1 index, of microorganisms containing phoD and pqqC genes was obtained before rarefying all sequences at a minimum number of reads per sample, 93674 and 89695, respectively. Community dissimilarity of pqqC and phoD microbial communities were estimated based on Bray–Curtis distances between samples, and the distances were calculated using OTU tables resampled to a minimum number of sequences from each sample.

Quantification of environmental stressors

Five environmental stressors, i.e., acidification, metal pollution, salinization, drought and human disturbances, which are common environmental factors caused by global change and human activities²¹, were used in this study. These stressors have been chosen by their prevalence across the world and their reported ecological significance on soil ecological services and microbial communities^20,21. The aridity index indicating drought degree was estimated using the Global Aridity Index and Potential Evapotranspiration Climate Database v3 (https://csidotinfo.wordpress.com/2019/01/24/global-aridity-index-and-potential-evapotranspiration-climate-database-v3/). The human influence index indicating human disturbance intensity was obtained from the Global Human Influence Index Dataset of the Last of the Wild Project, Version2, 2005 (LWP-2) (https://sedac.ciesin.columbia.edu/data/set/wildareas-v2-human-influence-index-geographic). Soil SAR (sodium adsorption ratio), a key index indicating soil salinization degree, is the ratio of Na⁺ relative to Ca²⁺ and Mg²⁺ in the water extract (solution phase) from a saturated soil paste⁵⁵:

1

$$SAR=\frac{{Na}^{+}}{\sqrt{\frac{{Ca}^{2+}+{Mg}^{2+}}{2}}}$$ where the base cation concentrations were in mmol(c) L^–1 as measured in solution. We used NPI (the integrated Nemerow pollution index) to evaluate the degree of soil heavy metals. The PIi (single pollution index) was calculated first⁵⁶:

2

$${PI}_i=\frac{C_i}{S_i}$$ where C_i is the actual concentration of heavy metal element i in soils (mg kg⁻¹). S_i is the soil background concentration of element i (mg kg⁻¹), and the background values were referred to Chen et al.⁵⁷. The integrated Nemerow pollution index (NPI) was calculated as:

3

$$NPI=\sqrt{\frac{MaxP{I_i}^2+AveP{I_i}^2}{2}}$$ $Max{PI}_i$ and $Ave{PI}_i$ are the maximum and average values of the single pollution indexes⁵⁸. We selected total metal concentrations as the primary stress metric because (1) background data for soil bioavailable fractions are lacking, preventing application of the Nemerow Pollution Index related to bioavailable data; and (2) total concentrations provide a more robust indicator of long-term cumulative effects. Soil pH was directly used to indicate the degree of acidification.

Individual and interactive effects calculation

For further analysis of interactive effects, we categorized the 180 samples into non-stress and stress groups using threshold methods⁵⁹. Samples were considered impacted by an environmental stressor (e.g., metal pollution) if the stressor’s level exceeded 50% of the maximum observed value. Using this 50% threshold, we classified the 180 samples with five environmental stressors into non-stress ( ≤ 50% threshold) and stress groups ( > 50% threshold). We further divided the samples into categories based on the presence of none, one, or two environmental stressors. Specifically, samples without any environmental stressors (i.e., all stressor values ≤ 50% threshold) were designated as the control (CK). Samples with only one environmental stressor above the 50% threshold were designated as Group A or Group B, depending on which stressor exceeded the threshold. Samples with both stressors exceeding the 50% threshold were classified as group A + B (Table S5). To ensure comparability, we first standardized all environmental stressor values to a range between 0 and 1 within each dataset. To maintain methodological consistency, we applied the same threshold criterion ( > 50% of the maximum observed value) to five environmental stressors. However, for soil pH, the threshold derived from this criterion (6.81) aligns more closely with alkaline soil conditions. Therefore, we additionally calculated effects using alternative pH thresholds (7.0, 6.5, 6.0, and 5.5). Our analysis revealed that the results obtained with a pH threshold of 7.0 were identical to those using the 50% maximum value criterion (Fig. S1). Consequently, in this study, soil acidification stress is defined as the process of converting alkaline soils to acidic soils (reaching a pH ≤ 7.0).

The individual and interaction effects of environmental stresses on the abundance and diversity of soil phosphorus-cycling microorganisms were calculated using Hedge’s d. The individual effects refer to the net effect of a specific environmental stressor. The individual effect of stressor A ( $d_A$ ), stressor B ( $d_B$ ) and their interaction ( $d_I$ ) were calculated as²:

4

$$d_A=\frac{\left(\overline{X_A}+\overline{X_{AB}}\right)-\left(\overline{X_B}+\overline{X_C}\right)}{2s}\times J\left(m\right)$$

5

$$d_B=\frac{\left(\overline{X_B}+\overline{X_{AB}}\right)-\left(\overline{X_A}+\overline{X_C}\right)}{2s}\times J\left(m\right)$$

6

$$d_I=\frac{\left(\overline{X_{AB}}-\overline{X_B}\right)-\left(\overline{X_A}-\overline{X_C}\right)}{s}\times J\left(m\right)$$ where $\overline{X_C}$ , $\overline{X_A}$ , $\overline{X_B}$ and $\overline{X_{AB}}$ are means of a variable in the control and treatment groups of A, B, and their combination (A + B), respectively. The standard deviation (s), degree of freedom (m) and correction term J were estimated by equations (7), (8) and (9), respectively, for the individual and interactive effects.

7

$$s=\sqrt{\frac{\left(n_C-1\right){\left(s_C\right)}^2+\left(n_A-1\right){\left(s_A\right)}^2+\left(n_B-1\right){\left(s_B\right)}^2+\left(n_{AB}-1\right){\left(s_{AB}\right)}^2}{n_C+n_A+n_B+n_{AB}-4}}$$

8

$$m=n_C+n_A+n_B+n_{AB}-4$$

9

$$J=1-\frac{3}{4\times m-1}$$ According to the methods of Folt et al.⁶⁰, we classified the interactions between two factors into three types: additive, antagonistic, and synergistic. The variances ( $v$ ) of $d_i$ were calculated according to Gurevitch et al.⁶¹. The two-tailed critical value of the standard normal distribution (Cα/2, i.e., 1.96) and standard error of $d_i$ (st) were used to calculate the 95% CI as Hedges’ dt ± Cα/2 × st⁶¹. If the 95% CI did not overlap with zero, the interactive effect was considered to be additive. For factor pairs whose $d_A$ and $d_B$ exhibit the same sign if $d_I$ and $d_A$ have the same sign, the interaction type is classified as synergistic. Conversely, if $d_I$ differs in sign from $d_A$ , the interaction is classified as antagonistic. When $d_A$ and $d_B$ possess opposing signs, the interaction type is determined as follows: if $d_A+d_B$ , aligns with the sign of $d_I$ , the interaction is labeled synergistic. Alternatively, if the $d_A+d_B$ differs in sign from $d_I$ , the interaction is classified as antagonistic²⁵. The additive effect means that the combined effect of two environmental stressors is simply the sum of their individual effects. The synergistic effect means that the combined effect of two stressors is significantly greater than the sum of their individual effects. The antagonistic effect means that the combined effect of two stressors is less than the sum of their individual effects.

Verification experiments

To verify the significant interactive effects of drought and acidification on phoD and pqqC gene abundances from large-scale survey, as well as the effects of drought and metal pollution stresses on the diversity of microbial communities with phoD and pqqC genes, we selected two types of soils: red soil (19.89° N, 110.59° E) and brown soil (35.91° N, 116.19° E). These soils are widely distributed across China and represent two major global soil types. Red soil exhibits a pH of 5.7, a total carbon content of 2.0%, and a texture of 45 clay, 35 silt, and 20% sand. Brown soil has a pH of 7.4, a total carbon content of 3.3%, and a texture of 30 clay, 48 silt, and 22% sand. For acidification stress, a mixture of HCl and H₂SO₄ (1:1, v:v) with a pH of 3.5 was gradually added to the soil using a peristaltic pump until the soil pH was reduced by one unit. The acidified soils were then incubated at 20 °C for 7 days at 60% water holding capacity (WHC). To introduce drought stress, soils at 60% WHC were air-dried at room temperature until the moisture content reached 5% WHC, followed by incubation at 20 °C for 7 days. For metal stress, the concentrations of As, Cd, Cr, Cu, Ni, Pb, and Zn were adjusted to 10, 0.3, 72, 22, 36, 21, and 81 mg kg⁻¹ (which were the median values in the stress group of metal pollution) respectively, by adding Na₂HAsO₄·7H₂O, CdCl₂·5/2H₂O, K₂Cr₂O₇, Cu(NO₃)₂·3H₂O, N₂O₆·Ni·6H₂O, HNO₃·1/2Pb, and Zn(NO₃)₂·6H₂O to the soils, which were then aged for two weeks to achieve realistic metal bioavailability and ecological relevance in experiments. Post-aging, available metal concentrations were measured and confirmed to belong in the range of metal contaminated soil group. Untreated soils served as controls, while soils with combined acidification+drought and drought+metal stress represented the interactive effects of dual stress factors.

The verification experiment included six treatments, each with four replicates: 1) control, 2) acidification stress, 3) drought stress, 4) metal stress, 5) acidification + drought stress and 6) drought+metal stress. The prepared soils were incubated in the dark at 20 °C, with deionized water periodically added by weight to maintain soil moisture. After 30 days of incubation, samples from control, acidification, drought, and acidification+drought treatments were collected for measurement of phoD and pqqC gene abundances by qPCR. Samples from control, metal, drought, and drought + metal treatments were collected to assess the diversity of microorganisms with phoD and pqqC genes using 16S sequencing. The qPCR and 16S sequencing procedures followed previously described methods.

Statistical analyses

We used the ‘relaimpo’ package in R to quantify the variance explained by each individual environmental stressor through a variance partitioning approach. In this context, predictors accounting for a significant portion of variability in the response variables are deemed more influential⁶². Spearman correlations were applied to assess the relationships between soil P-related parameters and the abundance and diversity of phosphorus-cycling microorganisms, using both the individual and interactive effects datasets. This analysis tests whether often-overlooked interactive effects on PCMs also contribute to changes in soil P cycling. To determine the contributions of individual and interactive effects of soil P properties on the abundance and diversity of P-cycling microorganisms, we performed a random forest analysis using the ‘rfPermute’ package in R⁶³.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Acknowledgements

This study was supported by the National Key Research and Development Program of China (2023YFD1501700), the National Science Foundation of China (41721001) and Fundamental Research Funds for the Central Universities in China.

Author contributions

Z.D., X.T. and J.X. developed the original idea, designed and supervised the study. X.T., Y.C. and R.Z. performed the experimental studies. X.T., Y.C., R.Z. and D.H. carried out the analysis. Z.D., X.T., S.C. and J.X. wrote and revised the manuscript.

Peer review

Peer review information

Communications Earth & Environment thanks Jiang Tian and Kazuki Fujita, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: D’Arcy Meyer-Dombard, Somaparna Ghosh [A peer review file is available].

Data availability

All raw sequences have been deposited in the Sequence Read Archive (SRA) with the accession number PRJNA1234860 (phoD) and PRJNA1235101 (pqqC). The source data used in this study are available in the Zenodo database with the https://doi.org/10.5281/zenodo.17046095.

Competing interests

The authors declare no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s43247-025-02772-6.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Cleveland, CC et al. Relationships among net primary productivity, nutrients and climate in tropical rain forest: a pan-tropical analysis. Ecol. Lett.; 2011; 14, pp. 939-947. [DOI: https://dx.doi.org/10.1111/j.1461-0248.2011.01658.x]

2. Rodrı́guez, H; Fraga, R. Phosphate solubilizing bacteria and their role in plant growth promotion. Biotechnol. Adv.; 1999; 17, pp. 319-339. [DOI: https://dx.doi.org/10.1016/S0734-9750(99)00014-2]

3. Gyaneshwar, P; Naresh Kumar, G; Parekh, LJ; Poole, PS. Role of soil microorganisms in improving P nutrition of plants. Plant Soil; 2002; 245, pp. 83-93.1:CAS:528:DC%2BD38XnvVCitrs%3D [DOI: https://dx.doi.org/10.1023/A:1020663916259]

4. Dai, ZM et al. Long-term nutrient inputs shift soil microbial functional profiles of phosphorus cycling in diverse agroecosystems. ISME J.; 2020; 14, pp. 757-770.1:CAS:528:DC%2BC1MXitl2lt73J [DOI: https://dx.doi.org/10.1038/s41396-019-0567-9]

5. Fraser, TD; Lynch, DH; Bent, E; Entz, MH; Dunfield, KE. Soil bacterial phoD gene abundance and expression in response to applied phosphorus and long-term management. Soil Biol. Biochem.; 2015; 88, pp. 137-147.1:CAS:528:DC%2BC2MXot1Wmu7c%3D [DOI: https://dx.doi.org/10.1016/j.soilbio.2015.04.014]

6. Rodríguez, H; Fraga, R; Gonzalez, T; Bashan, Y. Genetics of phosphate solubilization and its potential applications for improving plant growth-promoting bacteria. Plant Soil; 2006; 287, pp. 15-21. [DOI: https://dx.doi.org/10.1007/s11104-006-9056-9]

7. Meyer, JB; Frapolli, M; Keel, C; Maurhofer, M. Pyrroloquinoline quinone biosynthesis gene pqqC, a novel molecular marker for studying the phylogeny and diversity of phosphate-solubilizing pseudomonads. Appl. Environ. Microbiol.; 2011; 77, pp. 7345-7354.1:CAS:528:DC%2BC3MXhsVKnu7zN [DOI: https://dx.doi.org/10.1128/AEM.05434-11]

8. Dai, ZM et al. Elevated temperature shifts soil N cycling from microbial immobilization to enhanced mineralization, nitrification and denitrification across global terrestrial ecosystems. Glob. Change Biol.; 2020; 26, pp. 5267-5276. [DOI: https://dx.doi.org/10.1111/gcb.15211]

9. Liu, JL et al. Precipitation and soil pH drive the soil microbial spatial patterns in the Robinia pseudoacacia forests at the regional scale. Catena; 2022; 212, 106120.1:CAS:528:DC%2BB38XmtVaqsbk%3D [DOI: https://dx.doi.org/10.1016/j.catena.2022.106120]

10. Pastore, G; Kernchen, S; Spohn, M. Microbial solubilization of silicon and phosphorus from bedrock in relation to abundance of phosphorus-solubilizing bacteria in temperate forest soils. Soil Biol. Biochem.; 2020; 151, 108050.1:CAS:528:DC%2BB3cXitFyhtrzL [DOI: https://dx.doi.org/10.1016/j.soilbio.2020.108050]

11. Ragot, SA; Kertesz, MA; Bünemann, EK. phoD alkaline phosphatase gene diversity in soil. Appl. Environ. Microbiol.; 2015; 81, pp. 7281-7289.1:CAS:528:DC%2BC28XhslamtLw%3D [DOI: https://dx.doi.org/10.1128/AEM.01823-15]

12. Zhu, SG et al. Soil phosphorus availability and utilization are mediated by plant facilitation via rhizosphere interactions in an intercropping system. Eur. J. Agron.; 2023; 142, 126679.1:CAS:528:DC%2BB38XivVahsrvF [DOI: https://dx.doi.org/10.1016/j.eja.2022.126679]

13. Xu, HD et al. Long-term straw return promotes soil phosphorus cycling by enhancing soil microbial functional genes responsible for phosphorus mobilization in the rice rhizosphere. Agr. Ecosyst. Environ.; 2025; 381, 109422.1:CAS:528:DC%2BB2cXislelt7rN [DOI: https://dx.doi.org/10.1016/j.agee.2024.109422]

14. Wang, RB et al. Effects of short-term drought, nitrogen application and their interactions on the composition and functional genes of soil microbial communities in alfalfa grassland on the Loess Plateau. Front. Sustain. Food Syst.; 2023; 7, 1332683. [DOI: https://dx.doi.org/10.3389/fsufs.2023.1332683]

15. Allison, SD; Martiny, JB. Resistance, resilience, and redundancy in microbial communities. Proc. Natl. Acad. Sci. USA; 2008; 105, pp. 11512-11519.1:CAS:528:DC%2BD1cXhtVSjtrzO [DOI: https://dx.doi.org/10.1073/pnas.0801925105] supplement_1

16. Gao, YJ et al. Soil microbial functional profiles of P-cycling reveal drought-induced constraints on P-transformation in a hyper-arid desert ecosystem. Sci. Total Environ.; 2024; 925, 171767.1:CAS:528:DC%2BB2cXmtlygsrY%3D [DOI: https://dx.doi.org/10.1016/j.scitotenv.2024.171767]

17. Wan, WJ et al. Spatial differences in soil microbial diversity caused by pH-driven organic phosphorus mineralization. Land Degrad. Dev.; 2021; 32, pp. 766-776. [DOI: https://dx.doi.org/10.1002/ldr.3734]

18. Zhang, QF et al. Are the combined effects of warming and drought on foliar C: N: P: K stoichiometry in a subtropical forest greater than their individual effects? For. Ecol. Manage.; 2019; 448, pp. 256-266. [DOI: https://dx.doi.org/10.1016/j.foreco.2019.06.021]

19. Hu, MJ et al. Moderate salinity improves the availability of soil P by regulating P-cycling microbial communities in coastal wetlands. Glob. Change Biol.; 2023; 29, pp. 276-288. [DOI: https://dx.doi.org/10.1111/gcb.16465]

20. Rillig, MC et al. The role of multiple global change factors in driving soil functions and microbial biodiversity. Science; 2019; 366, pp. 886-890.1:CAS:528:DC%2BC1MXitFOnu7rE [DOI: https://dx.doi.org/10.1126/science.aay2832]

21. Rillig, MC et al. Increasing the number of stressors reduces soil ecosystem services worldwide. Nat. Clim. Change; 2023; 13, pp. 478-483. [DOI: https://dx.doi.org/10.1038/s41558-023-01627-2]

22. Eisenhauer, N; Cesarz, S; Koller, R; Worm, K; Reich, PB. Global change belowground: impacts of elevated CO₂, nitrogen, and summer drought on soil food webs and biodiversity. Glob. Change Biol.; 2012; 18, pp. 435-447. [DOI: https://dx.doi.org/10.1111/j.1365-2486.2011.02555.x]

23. Wang, XW et al. Microbial phosphorus-cycling genes in soil under global change. Glob. Change Biol.; 2024; 30, 17281. [DOI: https://dx.doi.org/10.1111/gcb.17281]

24. Yue, K et al. Individual and combined effects of multiple global change drivers on terrestrial phosphorus pools: a meta-analysis. Sci. Total Environ.; 2018; 630, pp. 181-188.1:CAS:528:DC%2BC1cXjt1Sjur8%3D [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.02.213]

25. Zhou, LY et al. Interactive effects of global change factors on soil respiration and its components: a meta-analysis. Glob. Change Biol.; 2016; 22, pp. 3157-3169. [DOI: https://dx.doi.org/10.1111/gcb.13253]

26. Crain, CM; Kroeker, K; Halpern, BS. Interactive and cumulative effects of multiple human stressors in marine systems. Ecol. Lett.; 2008; 11, pp. 1304-1315. [DOI: https://dx.doi.org/10.1111/j.1461-0248.2008.01253.x]

27. Shrivastava, M., Srivastava, P. C., & D’souza, S. F. Phosphate-solubilizing microbes: diversity and phosphates solubilization mechanism. Role of rhizospheric microbes in soil: vol. 2: nutrient management and crop improvement, 137-165 (Springer, 2018).

28. Louca, S et al. Function and functional redundancy in microbial systems. Nat. Ecol. Evol.; 2018; 2, pp. 936-943. [DOI: https://dx.doi.org/10.1038/s41559-018-0519-1]

29. Berlinches de Gea, A; Hautier, Y; Geisen, S. Interactive effects of global change drivers as determinants of the link between soil biodiversity and ecosystem functioning. Glob. Change Biol.; 2023; 29, pp. 296-307.1:CAS:528:DC%2BB38XislWhsLvM [DOI: https://dx.doi.org/10.1111/gcb.16471]

30. Zhou, X; Gu, X; Smaill, SJ. Rethinking experiments that explore multiple global change factors. Trends Ecol. Evol.; 2023; 38, pp. 399-401. [DOI: https://dx.doi.org/10.1016/j.tree.2023.01.009]

31. Naik, PR; Raman, G; Narayanan, KB; Sakthivel, N. Assessment of genetic and functional diversity of phosphate solubilizing fluorescent pseudomonads isolated from rhizospheric soil. BMC Microbiol.; 2008; 8, pp. 1-14. [DOI: https://dx.doi.org/10.1186/1471-2180-8-230]

32. Jones, D. L., & Oburger, E. Solubilization of phosphorus by soil microorganisms. Phosphorus in action: biological processes in soil phosphorus cycling, 169-198 (Springer, 2011).

33. Zhou, Z; Zheng, M; Xia, J; Wang, C. Nitrogen addition promotes soil microbial beta diversity and the stochastic assembly. Sci. Total Environ.; 2022; 806, 150569.1:CAS:528:DC%2BB3MXitFGhsLjI [DOI: https://dx.doi.org/10.1016/j.scitotenv.2021.150569]

34. Zheng, BX et al. Long-term nitrogen fertilization decreased the abundance of inorganic phosphate solubilizing bacteria in an alkaline soil. Sci. Rep.; 2017; 7, 1:CAS:528:DC%2BC2sXisVenur8%3D [DOI: https://dx.doi.org/10.1038/srep42284] 42284.

35. Yao, QM et al. Community proteogenomics reveals the systemic impact of phosphorus availability on microbial functions in tropical soil. Nat. Ecol. Evol.; 2018; 2, pp. 499-509. [DOI: https://dx.doi.org/10.1038/s41559-017-0463-5]

36. Luo, G et al. Soil carbon, nitrogen, and phosphorus cycling microbial populations and their resistance to global change depend on soil C: N: P stoichiometry. Msystems; 2020; 5, pp. 10-1128. [DOI: https://dx.doi.org/10.1128/msystems.00162-20]

37. He, M; Dijkstra, FA. Drought effect on plant nitrogen and phosphorus: a meta-analysis. New Phytol.; 2014; 204, pp. 924-931.1:CAS:528:DC%2BC2cXhvVOgurvL [DOI: https://dx.doi.org/10.1111/nph.12952]

38. Raymond, NS et al. Phosphate-solubilising microorganisms for improved crop productivity: a critical assessment. New Phytol; 2021; 229, pp. 1268-1277.1:CAS:528:DC%2BB3MXjsVemt7c%3D [DOI: https://dx.doi.org/10.1111/nph.16924]

39. Schimel, JP. Life in dry soils: effects of drought on soil microbial communities and processes. Annu. Rev. Ecol. Evol. S.; 2018; 49, pp. 409-432. [DOI: https://dx.doi.org/10.1146/annurev-ecolsys-110617-062614]

40. Kour, D. et al. Drought-tolerant phosphorus-solubilizing microbes: biodiversity and biotechnological applications for alleviation of drought stress in plants. Plant growth promoting rhizobacteria for sustainable stress management: Volume 1: Rhizobacteria in abiotic stress management, 255-308 (Springer, 2019).

41. Shi, W; Xing, Y; Zhu, Y; Gao, N; Ying, Y. Diverse responses of pqqC-and phoD-harbouring bacterial communities to variation in soil properties of Moso bamboo forests. Microb. Biotechnol.; 2022; 15, pp. 2097-2111.1:CAS:528:DC%2BB38XhtFKrs77O [DOI: https://dx.doi.org/10.1111/1751-7915.14029]

42. Yang, L; Du, L; Li, W; Wang, R; Guo, S. Divergent responses of phoD-and pqqC-harbouring bacterial communities across soil aggregates to long fertilization practices. Soil Tillage Res.; 2023; 228, 105634. [DOI: https://dx.doi.org/10.1016/j.still.2023.105634]

43. Preece, C; Verbruggen, E; Liu, L; Weedon, JT; Peñuelas, J. Effects of past and current drought on the composition and diversity of soil microbial communities. Soil Biol. Biochem.; 2019; 131, pp. 28-39.1:CAS:528:DC%2BC1MXhtlegtw%3D%3D [DOI: https://dx.doi.org/10.1016/j.soilbio.2018.12.022]

44. Xie, Y et al. Effect of heavy metals pollution on soil microbial diversity and bermudagrass genetic variation. Front. Plant Sci.; 2016; 7, 755. [DOI: https://dx.doi.org/10.3389/fpls.2016.00755]

45. Fierer, N. Embracing the unknown: disentangling the complexities of the soil microbiome. Nat. Rev. Microbiol.; 2017; 15, pp. 579-590.1:CAS:528:DC%2BC2sXhtlGktbfN [DOI: https://dx.doi.org/10.1038/nrmicro.2017.87]

46. Li, YQ et al. Effects of multiple global change factors on soil microbial richness, diversity and functional gene abundances: a meta-analysis. Sci. Total Environ.; 2022; 815, 152737.1:CAS:528:DC%2BB38XpvV2rsA%3D%3D [DOI: https://dx.doi.org/10.1016/j.scitotenv.2021.152737]

47. Gineyts, R; Niboyet, A. Nitrification, denitrification, and related functional genes under elevated CO₂: a meta-analysis in terrestrial ecosystems. Glob. Change Biol.; 2023; 29, pp. 1839-1853.1:CAS:528:DC%2BB3sXhvVGkurs%3D [DOI: https://dx.doi.org/10.1111/gcb.16568]

48. Laranjo, M; Alexandre, A; Oliveira, S. Legume growth-promoting rhizobia: an overview on the Mesorhizobium genus. Microbiol. Res.; 2014; 169, pp. 2-17. [DOI: https://dx.doi.org/10.1016/j.micres.2013.09.012]

49. Ryan, RP et al. The versatility and adaptation of bacteria from the genus Stenotrophomonas. Nat. Rev. Microbiol.; 2009; 7, pp. 514-525.1:CAS:528:DC%2BD1MXntlGns7s%3D [DOI: https://dx.doi.org/10.1038/nrmicro2163]

50. Waterlot, C. Alternative approach to the standard, measurements and testing programme used to establish phosphorus fractionation in soils. Anal. Chim. Acta; 2018; 1003, pp. 26-33.1:CAS:528:DC%2BC2sXhvV2mu7jF [DOI: https://dx.doi.org/10.1016/j.aca.2017.11.059]

51. de Santiago-Martín, A et al. Improving the relationship between soil characteristics and metal bioavailability by using reactive fractions of soil parameters in calcareous soils. Environ. Toxicol. Chem.; 2015; 34, pp. 37-44. [DOI: https://dx.doi.org/10.1002/etc.2772]

52. Zheng, BX; Zhu, YG; Sardans, J; Penuelas, J; Su, JQ. QMEC: a tool for high-throughput quantitative assessment of microbial functional potential in C, N, P, and S biogeochemical cycling. Sci. China Life Sci.; 2018; 61, pp. 1451-1462.1:CAS:528:DC%2BC1cXhsFOgsL7P [DOI: https://dx.doi.org/10.1007/s11427-018-9364-7]

53. Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. J.; 2011; 17, pp. 10-12. [DOI: https://dx.doi.org/10.14806/ej.17.1.200]

54. Bolyen, E et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol.; 2019; 37, pp. 852-857.1:CAS:528:DC%2BC1MXhsVeksr%2FO [DOI: https://dx.doi.org/10.1038/s41587-019-0209-9]

55. Rietz, DN; Haynes, RJ. Effects of irrigation-induced salinity and sodicity on soil microbial activity. Soil Biol. Biochem.; 2003; 35, pp. 845-854.1:CAS:528:DC%2BD3sXktFCqtrs%3D [DOI: https://dx.doi.org/10.1016/S0038-0717(03)00125-1]

56. Wu, W et al. Assessment of heavy metal pollution and human health risks in urban soils around an electronics manufacturing facility. Sci. Total Environ.; 2018; 630, pp. 53-61.1:CAS:528:DC%2BC1cXjtFClt7c%3D [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.02.183]

57. Chen, J; Wei, F; Zheng, C; Wu, Y; Adriano, DC. Background concentrations of elements in soils of China. Water Air Soil Pollut.; 1991; 57, pp. 699-712. [DOI: https://dx.doi.org/10.1007/BF00282934]

58. Ma, WC et al. Contamination source apportionment and health risk assessment of heavy metals in soil around municipal solid waste incinerator: a case study in North China. Sci. Total Environ.; 2018; 631, pp. 348-357. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.03.011]

59. Byrnes, JE et al. Investigating the relationship between biodiversity and ecosystem multifunctionality: challenges and solutions. Methods Ecol. Evol.; 2014; 5, pp. 111-124. [DOI: https://dx.doi.org/10.1111/2041-210X.12143]

60. Folt, CL; Chen, CY; Moore, MV; Burnaford, J. Synergism and antagonism among multiple stressors. Limnol. Oceanogr.; 1999; 44, pp. 864-877. [DOI: https://dx.doi.org/10.4319/lo.1999.44.3_part_2.0864]

61. Gurevitch, J; Morrison, JA; Hedges, LV; Morin, PJ. The interaction between competition and predation: a meta-analysis of field experiments. AM. NAT.; 2000; 155, pp. 435-453. [DOI: https://dx.doi.org/10.1086/303337]

62. Grömping, U. Relative importance for linear regression in R: the package relaimpo. J. Stat. Softw.; 2007; 17, pp. 1-27.

63. Archer, E. Package ‘rfPermute’. CRANhttps://doi.org/10.32614/CRAN.package.rfPermute (CRAN, 2016).