1 Introduction
Phosphorus (P) is an essential nutrient that governs plant growth and primary productivity in terrestrial ecosystems
[1–
3]. Its availability is particularly critical in highly weathered tropical and subtropical soils, wherein a large fraction of P is immobilized through interactions with iron and aluminum oxides
[4,
5]. As a largely immobile nutrient derived mainly from mineral weathering, P availability in the rhizosphere is a key determinant of plant P acquisition and overall ecosystem function
[6,
7]. However, soil P exists in diverse inorganic and organic forms, and the mechanisms by which plants access these different pools remain poorly understood
[8,
9]. Clarifying how plant diversity influences the bioavailability of different P forms in the rhizosphere is important for developing sustainable strategies to mitigate P limitation.
The rhizosphere, defined as the soil zone influenced by plant roots, is a critical hotspot for nutrient acquisition, where rhizodeposition and root turnover provide substantial carbon inputs
[8,
10]. In species-rich plant communities, elevated productivity increases P demand
[11]. Concurrently, enhanced carbon availability stimulates microbial growth and drives compositional shifts toward P cycling-related bacterial phyla, such as Bacteroidetes, which facilitate organic P mineralization and inorganic P solubilization
[12,
13]. These interactions form a plant–microbe feedback loop that enhances P bioavailability in diverse ecosystems
[14,
15]. A mechanistic explanation for sustained P availability in species-rich systems requires closer examination of plant–microbe–soil interactions.
Plant species richness can enhance P acquisition through two key mechanisms: resource partitioning and direct facilitation
[16]. Resource partitioning occurs when co-occurring species employ divergent spatial, temporal, or chemical strategies to access different soil P pools, thus increasing overall P-uptake efficiency
[17]. Direct facilitation involves certain plant species that enhance P availability by exuding compounds that mobilize sparingly soluble P forms
[9]. However, interspecific variation in P-mobilization capacity limits the explanatory power of species richness alone
[18]. Functional traits and their diversity may offer stronger explanatory power, as they capture variation in resource-use strategies that directly influence P availability
[19,
20]. Specifically, the community-weighted mean (CWM) of leaf economic traits reflects the dominant resource-use strategy of a community, typically associated with high nutrient-uptake efficiency. In parallel, the functional diversity in these traits captures the breadth of niche differentiation among coexisting species
[21,
22]. Considering both aspects provides a more complete basis for explaining how plant communities regulate P availability.
In this study, we conducted a tree diversity experiment along a gradient of species richness (1, 2, 4, and 6 species) to investigate how tree species richness (SR), tree leaf functional trait diversity (FD), CWM values of resource-acquisitive traits, tree productivity, soil organic carbon (SOC), and microbial communities jointly influence rhizosphere P bioavailability. We tested the following hypotheses: (a) Rhizosphere P bioavailability will be enhanced by both tree species richness and functional diversity, because complementary resource-use strategies and facilitation among co-existing species increase the community-level efficiency of accessing soil P pools. (b) Rhizosphere P bioavailability will be positively associated with higher CWM values of resource-acquisitive traits (e.g., high specific leaf area (SLA), leaf nitrogen content (LNC), and leaf phosphorus content (LPC)), as these traits directly enhance the capacity of the plant community to acquire and mobilize P. (c) The positive effects of tree diversity and CWM values of resource-acquisitive traits on P bioavailability are mediated by rhizosphere microbial communities through shifts in community composition and activity, which promote P solubilization and mineralization.
2 Material and methods
2.1 Study site and climate
This study was conducted at the Experimental Center of Tropical Forestry, Chinese Academy of Forestry (22°10′ N, 106°50′ E) in Pingxiang City, Guangxi Zhuang Autonomous Region, China. This site experiences a subtropical monsoon climate with a mean annual temperature of 22.3 °C and mean annual precipitation of 1300 mm. Situated at 350 m above sea level, the experiment was established in 2015 on a cleared plantation that previously harbored a Pinus massoniana monoculture. From the regional native tree species pool, eight species were selected: P. massoniana (P), Erythrophleum fordii (E), Castanopsis hystrix (C), Dalbergia odorifera (D), Michelia macclurei (M), Mytilaria laosensis (L), Aquilaria sinensis (A), and Magnolietia glauca (G). Among these, E and D are nitrogen-fixing species. Detailed species combinations for each diversity level are provided in Table S1.
2.2 Experimental design
A randomized complete block design was used, comprising nine tree species compositions (representing richness levels of 1, 2, 4, and 6 species) replicated across four blocks. The design included monocultures (12 plots), two-species mixtures (eight plots), four-species mixtures (8 plots), and six-species mixtures (eight plots), yielding a total of 36 experimental plots (Table S1). The total experimental area was 40 ha. Within each of four terrain-based blocks, plots ranging from 0.27 to 1.2 ha were randomly arranged. A 20 m × 20 m subplot was randomly established within each plot. To achieve intimate species mixing, trees were planted in a checkerboard pattern with alternating rows and individuals at 2 m × 2 m spacing, ensuring that each tree was directly bordered by a heterospecific neighbor. The soil was classified as Ferralsol according to the World Reference Base for Soil Resources.
2.3 Soil sampling
In 2023, nine pits were excavated in the topsoil layer (0–10 cm) after removing the forest floor in each plot. Rhizosphere soil was collected by systematically sampling fine roots with adhering soil from all tree species present; soil was isolated by gentle manual shaking of the roots
[23]. Samples from the same plot were combined into a composite sample for subsequent analysis.
2.4 Rhizosphere soil phosphorus fractions
Plant P acquisition was assessed using four extractants that target distinct P pools: (a) CaCl
2 for soluble P accessed via root interception and diffusion; (b) citrate for mineral-sorbed P mobilized by organic acid complexation; (c) enzymes for labile organic P released through hydrolysis; and (d) HCl for occluded inorganic P dissolved by proton extrusion
[24–
28].
To obtain operationally defined proxies for plant-accessible P pools, we used the chemical fractionation scheme, which uses extractants that simulate rhizosphere processes
[29]. Four P fractions were quantified along a lability gradient: CaCl
2-P (0.01 mol·L
−1 CaCl
2-extractable P), representing soluble inorganic P accessible via root interception and diffusion; Citrate-P (10 mmol·L−1 citrate-extractable P), representing the active inorganic P pool solubilized by organic acids; Enzyme-P (0.2 enzyme-unit-extractable P), representing readily hydrolysable organic P mineralized by extracellular enzymes; and HCl-P (1 mol·L
−1 HCl-extractable P), representing more stable, mineral-adsorbed inorganic P accessed via proton extrusion.
For each fraction, 0.5 g of fresh soil was shaken with 10 mL of the respective extractant for 3 h at 180 r·min
−1. After centrifugation (4000 ×
g, 30 min, 25 °C), supernatants were collected. Citrate-P and HCl-P extracts were diluted 10-fold and 20-fold, respectively, whereas CaCl
2-P and Enzyme-P were analyzed undiluted
[29]. All P concentrations were determined colorimetrically using the malachite green method on a Multiskan Spectrum (Tecan Infinite 200 Pro, Switzerland). CaCl
2-P, Citrate-P and HCl-P were analyzed within 30 days of sampling, and Enzyme-P was analyzed within 7 days.
2.5 Rhizosphere SOC and microbial biomass
Rhizosphere SOC was measured using an elemental analyzer (LECO TruSpec CN, USA). Microbial biomass carbon (MBC) was determined using the chloroform fumigation-extraction method
[30].
2.6 Forest productivity
All woody plants with diameter at breast height (DBH) ≥1 cm within each 20 m × 20 m plot were tagged, mapped, and measured for DBH and height during censuses in 2020 and 2024. Aboveground biomass was estimated using species-specific allometric equations for Chinese tree species
[31](Table S2). Tree productivity (t·ha
−1·a
−1) was calculated as the sum of biomass increments of surviving trees and the biomass of newly recruited trees between the two censuses
[32].
2.7 Leaf functional traits
Following standardized protocols
[33–
35], 3–5 representative trees per species per plot were selected, and 5–10 mature sun-exposed leaves were collected from each tree. Leaves from conspecific individuals were pooled into a composite sample. Leaf area (LA) was measured using a Li-3000 leaf area meter; leaf fresh weight (LFW) and dry weight (LDW, after oven-drying at 65 °C for 48 h) were recorded. SLA was calculated as LA/LDW, and leaf dry matter content (LDMC) was calculated as LDW/LFW. LNC and LPC were determined using the semi-micro Kjeldahl method and the molybdenum-antimony colorimetric method, respectively
[36].
2.8 Tree leaf functional trait diversity and Functional identity
SR was the number of species per plot. FD was calculated as abundance-weighted functional dispersion using the ‘
dbFD’ function in the R package
‘FD’ package
[37]. Functional identity was expressed as the community-weighted mean (CWM) of each trait, weighted by the species basal area
[38].
where pi is the relative contribution of species i to the community, and traiti is the trait value of species i.
2.9 Microbial community profiling
DNA was purified using the PowerSoil DNA Kit (Qiagen, Germany). The content was assessed using a FLUOstar Optima fluorescence plate reader (BMG Labtech) using the PicoGreen dye. Extracted DNA was stored at –80 °C until used. The V3–V4 region of the 16S ribosomal RNA (rRNA) gene was amplified using the primers 338F (5ʹ-ACTCCTACGGGAGGCAGCAG-3ʹ) and 806R (5ʹ-GGACTACHVGG GTWTCTAAT-3ʹ) to evaluate the soil bacterial community composition
[39]. The polymerase chain reaction (PCR) mix (20 μL) comprised 10 ng of DNA template, 2 μL dNTPs, 0.8 μL for each primer, 4 μL 5× FastPfu Buffer, 0.2 μL BSA, and 0.4 μL FastPfu polymerase. The cycling conditions were 95 °C for 3 min; 27 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 45 s, and 72 °C for 10 min.
The internal transcribed spacer 1 (ITS1) region was amplified using the primers ITS1F (5ʹ-CTTGGTCATTTAGAGGAAGTAA-3ʹ) and ITS2R (5ʹ-GCTGCGTTCTTCATCGATGC-3ʹ) to assess the soil fungal community
[40]. The PCR mix (20 μL) included 10 ng of DNA template, 2 μL dNTPs, 0.8 μL of each primer, 2 μL 10× Buffer, 0.2 μL BSA, and 0.2 μL rTaq polymerase. The cycling conditions were 95 °C for 3 min; 35 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 45 s; and 72 °C for 10 min.
Both bacterial and fungal amplicon libraries were sequenced on an Illumina MiSeq platform at Majorbio Biopharm Technology Co., Ltd. (Shanghai) using the TruSeq Sample Prep Kit with 2 × 150 bp paired-end sequencing. Sequencing data were analyzed using Trimmomatic, Usearch, and FLASH (version 7.1). Operational taxonomic units with 97% similarity were clustered using UPARSE (version 7.1), and chimeric sequences were identified and removed.
2.10 Statistical analyses
All statistical analyses were performed in R 4.1.0. One-way ANOVA with Tukey’s Honestly Significant Difference (Tukey’s HSD) post-hoc tests was used to compare rhizosphere P fractions, productivity, SOC, and microbial biomass across richness levels. The community composition of fungi and bacteria was quantified using non-metric multidimensional scaling (NMDS) based on Bray–Curtis distances. Pairwise permutational multivariate analysis of variance (PERMANOVA) was performed with 999 permutations to test the effect of SR on microbial community composition using Bray–Curtis dissimilarity matrices. The community-level leaf economics spectrum, representing the acquisitive-conservative trait continuum, was quantified using the first axis of a principal component analysis (PCA) based on the CWM of SLA, LDMC, LNC and LPC (CWM.PC1)
[41]. PCA was performed to characterize the rhizosphere P bioavailability, incorporating four P fractions: CaCl
2-P, Citrate-P, Enzyme-P, and HCl-P. In this case, positive PCA scores corresponded to higher overall P bioavailability. Differences in PCA scores representing acquisition strategies and rhizosphere soil P bioavailability across species compositions were evaluated using one-way ANOVA followed by Tukey’s HSD post hoc multiple comparison tests. Pearson’s correlation coefficients between rhizosphere P bioavailability and the predictor variables (SR, FD, CWM.PC1, productivity, SOC, MBC, and Bacteroidetes abundance) were analyzed. Linear mixed-effects models (LMMs) were fitted to examine the relationships between rhizosphere P bioavailability and predictors, with species composition as a random effect. Residuals of LMMs were assessed for normality using the Shapiro–Wilk test and for homoscedasticity (Table S3). Prior to analysis, all variables were Z-score standardized to meet normality assumptions.
3 Results
3.1 Patterns of rhizosphere soil P bioavailability across tree species richness
Both rhizosphere soil available P (Fig. S1) and bioavailable P fractions (Fig. 1) varied significantly with tree species richness (p < 0.05). Available P was higher in two-species mixtures than in monocultures (Fig. S1). Among the sequentially extracted fractions, CaCl2-P was significantly higher in four-species mixtures than in monocultures and six-species mixtures (Fig. 1(a)). Citrate-P was lower in monocultures than in two-species mixtures (Fig. 1(b)). In contrast, enzyme-P and HCl-P did not differ significantly across the richness gradient (Fig. 1(c,d)).
3.2 Forest productivity and rhizosphere soil organic carbon
Forest productivity and rhizosphere SOC showed clear shifts along the diversity gradient (Fig. 2). Forest productivity was higher in two-species mixtures than in monocultures (Fig. 2(a)). SOC was greater in 6-species mixtures relative to that in four-species mixtures (Fig. 2(b)).
3.3 Rhizosphere soil microbial biomass and community composition
MBC did not differ significantly across the richness gradient (Fig. 3(a)). NMDS and subsequent correlation analysis indicated that tree species richness significantly influenced the composition of the rhizosphere bacterial community, whereas no such effect was observed for the fungal community (Fig. 3(b,c)). Pairwise PERMANOVA confirmed that the bacterial community structure differed significantly across the richness gradient (Table S4).
3.4 Associations between plant acquisition traits and rhizosphere P bioavailability
CWM.PC1 explained 64.7% of the trait variation (Fig. 4(a)). Compositions containing legumes (E, PE, PECD, and PECDMM) scored positively on PC1 and differed significantly from non-legume compositions (Fig. 4(a,b)). A separate PCA of the four P fractions showed that the first axis accounted for 52.4% of the variance in rhizosphere P bioavailability (Fig. 4(c)). Legume-included mixtures (PE, PECD, and PECDMM) again aligned with higher P bioavailability along this axis. Notably, legume monocultures showed lower P bioavailability than that with legume-based mixtures (Fig. 4(d)).
3.5 Drivers of rhizosphere P fractions
LMMs indicated that SR had no significant direct effect on any P fraction (Fig. 5(a), Table S5). In contrast, FD showed positive relationships with Citrate-P, Enzyme-P, and HCl-P (Fig. 5(b), Fig. S2). Similarly, CWM.PC1 was positively associated with CaCl2-P, Citrate-P, and HCl-P (Fig. 5(c), Fig. S2).
Forest productivity positively influenced CaCl2-P, Citrate-P, and HCl-P (Fig. 6(a), Fig. S2), whereas SOC was correlated with Citrate-P and HCl-P (Fig. 6(b), Fig. S2). MBC was positively related to Citrate-P and Enzyme-P (Fig. 6(c), Fig. S2), as was the relative abundance of Bacteroidetes (Fig. 6(d), Fig. S2).
4 Discussion
4.1 Role of tree diversity in rhizosphere P bioavailability
Our findings are consistent with previous studies reporting a positive association between tree diversity and rhizosphere P bioavailability
[12,
18,
42]. However, our results further show that the increase in P bioavailability was not directly driven by tree species richness, but instead responded significantly and positively to tree leaf functional trait diversity. This finding supports the view that functional diversity may better capture the niche complementarity in resource-acquisition strategies that enhance P mobilization in mixed communities
[43], consistent with previous observations in diverse stands
[18,
44]. Extending this perspective to managed subtropical plantations, functional diversity appears to be a key driver of rhizosphere P bioavailability, whereas species richness alone plays a more limited role.
The enhancement of P bioavailability by functional diversity operates through two interdependent pathways mediated by forest productivity and SOC. First, higher functional diversity promoted forest productivity, which may have stimulated root exudation of organic acids, phosphatases, and protons, processes known to mobilize sparingly soluble P
[11,
25]. Increased productivity was also linked to a higher abundance of Bacteroidetes, a phylum encoding P-solubilizing genes that play a key role in soil P cycling
[45,
46]. Second, functional diversity increased the SOC, which acted as a fundamental substrate for microbial metabolism and biomass synthesis
[12,
47]. The resulting boost in microbial biomass could further enhance P-mobilizing processes in the rhizosphere
[48]. These results point to a close linkage between C and P cycling under higher functional diversity, although the underlying physiological and microbial processes remain to be tested directly.
4.2 The influence of tree functional identity
Beyond tree diversity, soil nutrient cycling was also influenced by the functional identity of the dominant species
[9,
49], as reflected in the CWM of acquisitive traits
[50,
51]. Communities dominated by resource-acquisitive strategies showed significantly higher rhizosphere P bioavailability
[22,
52]. This pattern aligns with previous findings that acquisitive species enhance the exudation of carboxylates, protons, and extracellular enzymes, all of which contribute to mobilizing sparingly available P
[18]. Additionally, the presence of nitrogen-fixing tree species intercropped with non-fixing counterparts was associated with higher rhizosphere P bioavailability
[9,
12], supporting the hypothesis that facilitative interactions between P-mobilizing and non-mobilizing plants enhance community-level P availability
[9,
53].
The dominance of acquisitive traits, which was positively associated with forest productivity, was in turn correlated with rhizosphere P bioavailability. Physiologically, traits such as high leaf nitrogen content supported Rubisco synthesis and improved photosynthetic water-use efficiency, thereby sustaining rapid growth
[22,
54]. Previous studies have shown that increased productivity stimulates key microbial groups involved in P cycling, creating a positive plant–microbe feedback that reinforces P availability in the rhizosphere
[12,
46].
4.3 Limitations and perspectives
We acknowledge several limitations of this study. First, pooled rhizosphere sampling—although appropriate for community-scale assessment—limits inference on species-specific contributions; thus, the observed mixture effects may reflect both selection effects (e.g., nitrogen-fixers) and complementarity effects
[55,
56], which our data cannot disentangle. Second, our mechanistic interpretations (e.g., organic acid exudation and phosphatase activity) are largely inferred from correlational evidence. Although we observed consistent relationships among functional diversity, microbial attributes, and P bioavailability, causality cannot be established because we did not directly quantify root exudates, enzyme activities, or functional gene expression. Third, the measured phosphorus fractions measured (CaCl
2-P, Citrate-P, Enzyme-P, and HCl-P) are operationally defined by sequential extraction and not by the direct measurement of plant-available P. Although this scheme offers insights into P mobilization, these fractions indicate potential bioavailability rather than actual plant uptake. Future studies should therefore incorporate direct measurements (e.g., root exudates, enzyme activities, metagenomics, and
in-situ P uptake assays) to validate the proposed pathways.
5 Conclusions
Our study provides evidence that tree leaf functional trait diversity, rather than tree species richness, is closely associated with rhizosphere P bioavailability in subtropical plantations. This enhancement was found to be increased forest productivity and SOC, which were associated with shifts in microbial biomass and the abundance of Bacteroidetes. Resource-acquisitive tree species were also positively related to P bioavailability through increased productivity. Importantly, mixtures of nitrogen-fixing and non-fixing species proved particularly effective at enhancing rhizosphere P availability. These findings suggest that may benefit from focusing on tree leaf functional trait diversity rather than species richness alone, particularly by incorporating species with strong P-mobilizing capacity to support long-term productivity.
The Author(s) 2027. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)