Soil Quality Assessment for Sustainable Management: A Minimum Dataset for Long-Term Fertilization in Subtropical Plantations in South China

1. Introduction

Subtropical forests comprise 11% of the world’s forest ecosystems [1]. In areas with increasing population density and rapid economic development, subtropical forest ecosystems are often transformed into alternative land-use practices that prioritize the production of wood and other forest products [2]. However, such forms of conversion are typically unsustainable, resulting in significant degradation of the land and forest ecosystems [3-5]. Protective forest plantations play a pivotal role in preventing soil erosion and ecosystem restoration. Evidence suggests that fast-growing trees can promote erosion control through mechanisms including soil improvement and water regulation [6, 7]. As a result, extensive plantations have been established in southern China since the 1980s. Afforestation efforts in this region frequently utilize fast-growing tree species such as Eucalyptus urophylla and Acacia auriculiformis [8].

Managing fast-growing monoculture plantations for long-term sustainability requires careful consideration of risks including potential nutrient limitation [9-11]. Subtropical forest soils, characterized by high weathering, often face phosphorus (P) limitation [12-14], which can hinder tree growth, reduce soil fertility, and potentially alter nutrient cycling processes [15-17]. Furthermore, subtropical areas in China have been regarded as nitrogen (N) deposition hotspots primarily attributed to the increased consumption of N fertilizers and fossil fuels in recent decades, which might experience even greater P limitation [18, 19]. The resulting imbalance between N and P supply in the subtropical plantations can have negative consequences for soil quality and the sustainability of plantation ecosystems [20-22]. Ren et al. [23] found that most of these monoculture plantations have low productivity, thereby offering insufficient forest goods and ecosystem services.

Given that both N and P are particularly important nutrients for ecosystem structure, processes and function, fertilization has been one of the intensive forest management practices in plantations [24-27]. Results from previous studies demonstrate a strong and consistent association between fertilization and changes in soil nutrients. Lu et al. [28] reported that N fertilization significantly increased soil carbon (C) levels in subtropical forests via a decade-long N addition experiment and meta-analysis. Meanwhile, P fertilization was found to improve N-induced negative effects including nitrate leaching and nitrous oxide emission. While combined N+P fertilization is supposed to have synergistic effects on soil microbial activity and nutrient dynamics [29]. According to Fan et al. [30], soil physicochemical indicators including available P and related enzyme activity were considerably changed under N and P fertilizers. Overall, these studies illustrate how different fertilizer applications affect soil nutrient pools and the effects can be site-specific. However, few studies have quantitatively and comprehensively examined the consequences of long-term fertilization practices on soil quality of subtropical plantations.

Soil quality, defined as “the capacity of a soil to function within ecosystem and land-use boundaries to sustain biological productivity, maintain environmental quality and promote plant and animal health” [31], is a multifaceted concept for assessing the impacts of management practices [32-34]. By focusing on a specific soil threat like nutrient imbalance and a key function like nutrient supply, soil quality assessment can provide more targeted insights into the impacts of management practices. Building on Bünemann et al. [35], who emphasized the importance of assessing changes in soil quality status to devise adaptive management strategies, this study aims to evaluate the effects of long-term fertilizer application on soil quality in subtropical plantation soils via soil quality assessment.

Although soil quality cannot be measured directly, certain soil physical, chemical and biological indicators offer valuable information about soil processes and functions [36-39], and a minimum dataset (MDS) is often preferred because it can significantly reduce laboratory analysis workload and associated costs [40-43]. An MDS represents a subset of carefully chosen soil indicators selected from the total dataset (TDS), offering sufficient information for comprehensive soil quality evaluations. Many multivariate techniques including principal component analysis (PCA), redundancy analysis and factor analysis can capture the most relevant information while minimizing redundancy and complexity [35]. Among them, PCA is a popular method to identify MDS and has been extensively applied in various soil quality assessment studies [44-47].

In this study, we took advantage of the long-term fertilization experimental plot in the Acacia auriculiformis and Eucalyptus urophylla plantations starting in 2010. We hypothesized that Acacia auriculiformis plantation would exhibit higher soil quality than that of Eucalyptus urophylla plantation, due to the N-fixing capabilities of leguminous species in the former. Previous studies have reported positive impacts on soil nutrient pools and microbial communities in plantations dominated by leguminous species such as Acacia auriculiformis [48-50]. Given the differential responses of soil indicators to fertilization, we refrained from making any a priori predictions regarding SQI variation among treatments. Our investigation focused on the key questions: 1) How do long-term fertilization regimes affect a range of potential soil quality indicators? 2) How does soil quality vary among different fertilization regimes and plantations as measured by an SQI? 3) Can an MDS be identified that effectively represents the overall soil quality with long-term fertilization ?

2. Materials and Methods

2.1. Site Description

Soil quality indicators were assessed in two forest plantations located at the Heshan National Field Research Station of Forest Ecosystems (112°50′E, 22°34′N), Guangdong Province, South China (Figure 1A). The station is situated within a humid subtropical monsoon zone. The annual precipitation is about 1300 mm, with most rainfall occurring between May and September. The mean annual temperature is 21.7 °C. The terrain is hilly, with slopes ranging from 8° to 20°. Both plantations were established in 1984, covering areas of 4.6 ha and 1.9 ha, respectively [51]. Historically, unsustainable exploitation practices and inadequate land-use management have resulted in the degradation of the original forest ecosystem (Figure 1B, C). To address this issue, restoration efforts began in 1980s with the introduction of various tree species. The plantations investigated in this study were dominated by two different vegetation types: leguminous plants and non-leguminous plants. The leguminous plantation was primarily composed of Acacia auriculiformis, while the non-leguminous plantation was dominated by Eucalyptus urophylla (Figure 1D, E). These are typical fast-growing species for restoration in southern China and both were 40 years old. Details of the two plantations are presented in Tables S1 and S2.

2.2. Experimental Design and Sampling

In August 2010, an N and P fertilization study was established in the Acacia auriculiformis and Eucalyptus urophylla plantation as a randomized complete block design with three blocks, and three fertilization treatments varying in application dose. The field experiment included a control plot without fertilizer application (Control), an N fertilizer of 50 kg N ha^-1 yr^-1 (N fertilizer), a P fertilizer of 50 kg P ha^-1 yr^-1 (P fertilizer) and an N+P fertilizer of 50 kg N ha^-1 yr^-1 + 50 kg P ha^-1 yr^-1 (NP 50 fertilizer). Within each plantation, twelve 10 m × 10 m plots were established, separated by 10 m wide buffer zones to minimize interference between different fertilization treatments. Ammonium nitrate (NH₄NO₃) and/or dihydrogen phosphate (NaH₂PO₄) dissolved in 10 L water was used for the application. The fertilizer solution was then sprayed evenly on each plot every two months starting in August 2010. The control plots received 10 L of water during each fertilizer application.

The soil was sampled in September (wet season) and December (dry season) of 2023, respectively in the two plantations. Five soil samples were randomly collected from the 0-20 cm soil layer in each plot and mixed into a single composite soil sample. In this way, a total of 48 (4 treatments × 3 replicates × 2 plantations × 2 seasons) soil samples were collected. The soil samples were brought to the lab, immediately sieved in a 2 mm mesh, and then divided into two parts. One part was air-dried at room temperature and the other part was stored at -4 °C for further analyses. We collected litter within each plot during the wet (July-September 2023) and dry (October-December 2023) seasons, using a 100 × 100 cm sampling frame. Three bags of litter samples were collected from each subplot and then dried at 105 °C to determine the litter mass.

2.3. Laboratory Analysis

A total of eighteen physical, chemical, and biological indicators were measured in each collected soil sample. The soil physical indicators included soil water content (SWC) and soil structure. SWC was determined by gravimetric method with oven drying [52]. Particle size distributions of sand (0.05–2.0 mm), silt (0.002–0.05 mm) and clay (< 0.002 mm) were determined by the laser diffractometry method [53]. The soil chemical indicators included pH, soil organic C (SOC), total N (TN), inorganic N (IN), total P (TP) and cation exchange capacity (CEC). Soil pH was measured using a combination glass electrode meter (Mettler Toledo, Zurich, Switzerland) in a 1:2.5 soil-to-water slurry. SOC and TN were measured using an elemental analyzer (vario Max cube, Elementar, Germany). Soil IN, including ammonium (NH₄⁺) and nitrate (NO₃^-), were extracted by potassium chloride at a 1:2.5 soil–water ratio and determined using a continuous segmented flow analyzer (AA3, SEAL Analytical, Germany). Soil TP was measured using an X-ray fluorescence spectrometer (CNX-808, NCS, China). The exchangeable cations (K⁺, Ca²⁺, Na⁺, Mg²⁺, Al³⁺, Fe³⁺, and Mn²⁺) were determined following protocols as described previously [54], and CEC was calculated as the sum of the exchangeable cations. The soil biological indicators included microbial biomass C (MBC), microbial biomass N (MBN) and potential enzyme activities related to C, N and P cycle including β-glucosidase (BG), xylanase (XYL), β-D-cellobiohydrolase (CBH), β-N-acetylglucosaminidase (NAG), phosphodiesterase (PDE) and phosphomonoesterase (PME). MBC and MBN were determined by the chloroform fumigation method [55, 56]. Potential enzyme activities were measured using the microporous plate fluorescence method based on Bell et al. [57].

2.4. The Minimum Dataset

A literature study was conducted using Core Collection databases from the Web of Science®, - Clarivate Analytics platform (https://webofknowledge.com/) to identify suitable indicators to be included in the MDS of fertilization experiments. We set the topic retrieval as (forest OR plantation) and (fertiliz*) and ("soil quality" OR "soil health"), and only research articles written in English were included. After downloading all lists of articles, a screening process for article eligibility was performed (Figure S1). MDS were collected and summarized from 32 publications to identify the most suggested indicators (Table S4). Three main steps were followed to conduct soil quality assessment in this study: (i) apply PCA to identify MDS, (ii) apply standard scoring functions to normalize the MDS indicators, and (iii) integrate the indicator scores into a soil quality index for quantification and comparison [32, 58]. PCA helps in reducing the dimensionality of the soil indicators measured by transforming the original data into a smaller set of uncorrelated variables called principal components (PCs).

The normalized data of soil physical, chemical and biological properties underwent PCA separately [59]. PCs with eigenvalues ≥ 1 were selected according to Kaiser's criterion to ensure that each retained PC explains at least as much variance as a single original variable. For each PC, only the indicators with an absolute loading value within 10% of the highest loading were retained. Pearson's correlation coefficient was used to assess the relationships between variables when more than one indicator was retained in each PC. If highly loaded indicators have a high correlation (p < 0.05), it indicates redundancy, and so only a representative indicator was chosen for MDS. Otherwise, each indicator was retained in the MDS. After determining the MDS for the SQI, each soil indicator was transformed into a unitless score ranging from 0.00 to 1.00 using the standard scoring function [32, 42]. The indicators were considered good (i.e., more is better) when increasing values indicate positive soil quality or bad (i.e., less is better) when decreasing values reflect better soil quality. “More is better” (Equation 1) and “less is better” (Equation 2) functions were used as follows:

$$S={\displaystyle \frac{X}{X_{max}}}$$

(1)

$$S={\displaystyle \frac{X_{min}}{X}}$$

(2)

where S is the score of the soil indicator, X is the soil indicator value, and X_max, X_min are the maximum and minimum values of each soil indicator observed among the two plantations.

The transformed indicator scores were integrated into an SQI using the following formula (Equation 3):

$$SQI={\sum }_{i=1}^n{W}_i·S_i$$

(3)

where W_i and S_i are, respectively, the weight and the indicator score. For validation of the relationship between TDS and MDS in SQI, W_i is derived from the communality of the indicator in each PCA for TDS and the variation in each PC for MDS, respectively. Communality refers to the proportion of variance in a soil indicator that is accounted for by all the retained principal components. For indicators in MDS, the total weight (i.e., 1) is evenly distributed among three categories of physical, chemical and biological indicators. Specifically, for uncorrelated indicators within a PC, weights were assigned equally, corresponding to the PC's proportion of the total. Conversely, correlated indicators within a PC had their weights determined by dividing the PC's explained variance proportion among them [33, 45, 60, 61].

2.5. Data Analysis

Data analysis was carried out using R (version 4.3.2). After confirming the normality of data and homogeneity of variance, one-way analysis of variance (ANOVA) was conducted to test the effect of fertilizer applications on individual soil indicators, and Tukey's Honestly Significant Difference (HSD) test was used to compare the mean difference of soil indicators among fertilizer treatments (p < 0.05). A Pearson correlation analysis was conducted to assess the correlations between soil indicators to avoid redundancy. PCA was performed using the “FactoMineR” package [62].

3. Results

3.1. Effect of Fertilization on Soil Indicators

The soil was acidic, with a pH ranging from 3.74 to 4.10 and classified as loam texture according to United States Department of Agriculture (USDA, 2017) classification system. The physical indicators were relatively stable compared to the chemical and biological indicators; the latter ones showed large variation across different fertilizer treatments (Table S3). Notably, fertilization practices significantly impacted the soil nutrient pool. P treatment resulted in an average decrease of IN by 25% and 33% in the Acacia auriculiformis and Eucalyptus urophylla plantations, respectively. Similarly, NP led to IN reductions of 27% and 21% on average in Acacia auriculiformis and Eucalyptus urophylla plantations, respectively. Furthermore, it was noteworthy that both P and NP treatment significantly increased TP while decreasing enzyme activities associated with the P cycle (i.e., PME and PDE). In the Acacia auriculiformis plantation, TP increased by an average of 101% under both P and NP treatment. The Eucalyptus urophylla plantation displayed similar trends, with TP increases of 64% and 92% on average under P and NP treatment, respectively.

3.2. Principal Component Analysis

The results of PCA for physical indicators showed that the first 2 PCs had eigenvalues > 1.0, explaining 49.8% and 33.1% of the total variance, respectively (Table 1). Together, they explained a cumulative variance of 82.9%. The variable vector for SWC and silt pointed in a positive direction on PC1, indicating a positive correlation between these variables and PC1 (Figure 2). The variable vector for sand pointed in a negative direction on PC1, likewise, PC1 was negatively associated with sand (the associated loading was -0.87). PC2 on the other hand, was largely influenced by clay (loading: -0.95).

The results of PCA for soil chemical indicators showed that the first 2 PCs had eigenvalues > 1.0, explaining 45.7% and 26.2% of the total variance, respectively (Table 1). In total, the first 2 PCs captured 71.91% of the total variation. PC1 was highly dependent on SOC, TN, and IN (loading: 0.90, 0.86, and 0.75). PC2 had a strong positive correlation with TP (loading: 0.88). With the variable vector for TP pointed in the direction of the cluster of P and NP points and away from the cluster of C and N points, it represented that both the P and NP treatment affected soil TP. Furthermore, indicators including IN, TN, and SOC were positively correlated because they clustered together closely.

The results of PCA for biological indicators showed that the first 2 PCs had eigenvalues > 1.0 and accounted for 66.9% and 15.0% of the total variance, respectively (Table 1). Together, they explained a cumulative variance of 81.9%. All the variable vectors pointed in a positive direction on PC1, indicating a positive correlation between these variables and PC1. Meanwhile, these variables themselves were positively correlated because the variable vectors were close together and they are all highly contributing variables in PC1. PC2 was significantly impacted by PDE (loading: 0.74) compared with other variables.

3.3. Selection of the MDS

While the measurement methods for these soil indicators may vary, they reflect the content or level of certain properties in the soil based on our literature study (Figure 2). The most used soil physical indicators are soil water storage and bulk density (> 20%), followed by porosity and soil texture. Soil organic carbon/matter (19%), and various forms of macronutrients including N, P, K, Ca, and Mg are the most frequently used chemical indicators. Additionally, soil pH is another commonly measured chemical indicator. The most utilized biological indicators comprise soil enzymes (35%), along with MBC and MBN.

Our analysis of soil physical properties revealed that factors contributing most to the overall composition (PC1) were silt and sand, however, silt and sand were significantly correlated (r = -0.79, p < 0.05) (Figure 3A). Additionally, the adsorption of organic matter onto silt particles is critical for stabilizing soil organic matter [63], therefore only silt is retained in MDS. While SWC did not have a high weight in PC1, it was the most frequently proposed physical indicator as revealed in the literature study. Given its important role in nutrient cycling and storage, SWC was also retained in MDS. For PC2, clay had the highest loading and was retained in MDS.

In the PCA for chemical indicators, the highly weighted indicators for PC1 were SOC and TN. They were significantly correlated with each other. SOC is the main source of nutrients and the most common chemical indicator used in MDS as revealed before, therefore it was not excluded from MDS (Figure 3B). Since TN (r = 0.94, p < 0.05) had a stronger correlation with SOC than IN (r = 0.55, p < 0.05) and IN was a more sensitive indicator in the ANOVA (Table S3), only IN was retained in MDS. TP had the highest loading in PC2 and was therefore retained in MDS.

In the PCA for biological indicators, BG, CBH, XYL, NAG and MBC had the highest loadings in PC1 and demonstrated significant correlations with each other. To reduce redundancy, only MBC was retained in MDS (Figure 3C). PDE as the highly weighted indicator in PC2 was retained in MDS. Finally, based on PCA results, literature study and the management goal, a total of eight soil indicators were chosen as the MDS, including three physical indicators (silt, clay, SWC), three chemical indicators (SOC, IN, TP) and two biological indicators (MBC, PDE).

The soil indicator values were normalized between 0 and 1 using a standard scoring function. The “more is better” function was used for all the indicators except for silt and sand. High sand content weakens soil's ability to retain nutrients while high silt content can increase susceptibility to topsoil erosion. Therefore “less is better” function was applied to silt and sand. Weights were assigned to the indicators and standardized to unity (Table 1). For the indicators in TDS, weights were assigned by the communality of each indicator. For the indicators in MDS, weights were assigned by the variation in each PC (Figure 4).

The SQI in TDS was calculated using the following formula (Equation 4):

$$\begin{gathered} SQI-TDS=0.038·S_{WC}+0.054·S_{clay}+0.06·S_{silt}+0.065·S_{sand}+ \\ 0.052·S_{pH}+0.064·S_{SOC}+0.064·S_{TN}+0.59·S_{TP}+0.054·S_{IN}+0.048· \\ {S}_{CEC}+0.053·S_{MBC}+0.042·S_{MBN}+0.056·S_{BG}+0.056·S_{CBH}+0.055· \\ {S}_{XYL}+0.061·S_{NAG}+0.06·S_{PDE}+0.057· \\ {S}_{PME} \end{gathered}$$

(4)

The SQI in MDS was calculated using the following formula (Equation 5):

$$\begin{gathered} SQI-MDS={\displaystyle \frac{1}{3}}\times \left({\displaystyle \frac{0.249·S_{silt}+0.249·S_{WC}+0.330·S_{clay}}{0.828}}+{\displaystyle \frac{0.228·S_{SOC}+0.228·S_{IN}+0.262·S_{TP}}{0.719}}+ \\ {\displaystyle \frac{0.669·S_{MBC}+0.150·S_{PDE}}{0.819}}\right)=0.1·S_{silt}+0.1·S_{WC}+0.133·S_{clay}+ \\ 0.106·S_{SOC}+0.106·S_{IN}+0.121·S_{TP}+0.272·S_{MBC}+0.061· \\ {S}_{PDE} \end{gathered}$$

(5)

where S_i is the score of the soil indicators and the coefficients are the weights.

3.4. Soil Quality Index

The SQI values ranged from 0.63 to 0.72 in the Acacia auriculiformis plantation and from 0.60 to 0.66 in the Eucalyptus urophylla plantation (Figure 5A). On average, the SQI was 6% higher in the Acacia auriculiformis plantation compared to the Eucalyptus urophylla plantation, except under the NP treatment (Figure 5B). However, statistical analysis indicated that this difference was not significant. Across different fertilization regimes, the SQI of all treatments were at similar levels, suggesting comparable soil quality under different treatments (Figure 5C). Nonetheless, the observed changes in specific soil indicators such as TP, PDE, IN, and SOC implied that fertilizer treatments may have influenced certain aspects of soil function and nutrient dynamics.

3.5. SQI Evaluation

Validation of the MDS is important [40], therefore a linear model was used to assess the relationship between MDS and TDS (Figure 6), revealing a strong correlation (R² = 0.81). This indicates the MDS effectively captured key soil quality attributes. Additionally, a linear regression analysis demonstrated a positive relationship between litter mass and SQI (R² = 0.37), suggesting that organic matter inputs positively influence soil quality (Figure 7).

4. Discussion

This study investigated a total of eighteen soil indicators which can generally represent soil nutrient pools, focusing on observing the response of plantation soil to different fertilization practices. Generally, fertilization had a significant impact on soil indicators. Soil biological, chemical and physical soil indicators interact in intricate ways, measuring the status of these indicators provides valuable insights into these complex interactions, a critical step towards achieving sustainable forest management practices.

4.1. Effects of Fertilization on Soil Indicators

The clear response of the soil TP and P-acquiring enzyme activities following long-term P treatment strongly suggests potential P limitation in these subtropical plantations. Meanwhile, P treatment resulted in a decrease of IN in both plantations, with a more significant reduction observed during the wet season. This decrease in IN could be due to several factors, including increased N demand by microbes and plants stimulated by P fertilization, as well as potential N losses through leaching processes in these plantations [64, 65]. It was worth noting that P treatment decreased SOC by an average of 28% in the Eucalyptus urophylla plantation. This aligns with the findings of Xia et al. [66] who found accelerated SOC mineralization after P input, likely due to enhanced microbial activity breaking down more complex SOC fractions.

Compared to P, N is less likely to be a limiting element in subtropical plantations due to N deposition, which might explain the weaker influence of N fertilizer on soil indicators. However, a noteworthy observation was a 35% decline in SOC in the Acacia auriculiformis plantation during the wet season while in the dry season, the SOC level stabilized. This decrease might be due to faster decomposition of surface litter during the wet season leading to a temporary decline in SOC, while deeper SOC layers remain relatively stable [30].

NP treatment can have complex effects on the soil nutrient pool in subtropical plantations. Due to the biological N fixation of leguminous plants, the background TN in the Acacia auriculiformis plantation was 1.2 g/kg higher than that of Eucalyptus urophylla plantation. NP treatment increased TN by an average of 9%, whereas N and P treatment decreased TN by an average of 18% and 27% respectively in the Eucalyptus urophylla plantation. NP treatment may have synergistic interactions and facilitate optimal N utilization by the non-leguminous plants, resulting in an overall increase in TN levels [67-69]. Application of N alone (N treatment) without proportional P availability, however, could lead to nutrient imbalances and N leaching contributing to the observed decrease in TN levels [70].

4.2. The Minimum Dataset of Two Plantations

Readily measurable and management-sensitive soil indicators are crucial for effective soil quality assessment. Existing frameworks like Soil Management Assessment Framework (SMAF) and Comprehensive Assessment of Soil Health (CASH), while valuable for farmlands, may not fully capture the complexities of forest ecosystems [58, 71]. Forest soils often have unique characteristics and challenges compared to farmlands, necessitating site-specific assessment methods. Focusing solely on soil physical and chemical properties can overlook the dynamic aspects of soil quality while relying solely on biological indicators might not capture long-term changes. Some studies assign weights based on specific soil functions, but the interconnected nature of soil indicators makes it challenging to definitively link each indicator to a single function [72]. To capture a comprehensive picture of soil quality and the interplay between these aspects, we assigned equal weight to all three types of indicators (physical, chemical and biological) in the soil quality assessment. The literature review (Table S5) showed that although some indicators like water storage, SOC and potential soil enzyme activity are frequently used in MDS, most selection processes consider site-specific properties such as soil texture, climatic conditions and management goals [73-75]. In this study, we chose three physical indicators (silt, clay, SWC), three chemical indicators (SOC, IN, TP) and two biological indicators (MBC, PDE) as the MDS, allowing us to evaluate the impact of fertilization practices on the soil quality in the plantations. These indicators were chosen because they represent key aspects of the overall nutrient pool and their responses to management practices allow us to track changes in soil quality over time.

SWC and soil physical structure play the important roles in nutrient delivery for biological activity and plant growth. Clay and silt content are key factors in the stabilization of soil organic matter through mechanisms including physical protection [63, 76-78]. Notably, these indicators exhibited negligible seasonal fluctuations, suggesting relative stability of soil physical properties under fertilization.

SOC, a primary source of energy and substrate for nutrient mineralization, exhibited surprisingly different responses to fertilization depending on the season and plantation. Fertilization significantly impacted SOC, with an average decrease of 32% observed in the Acacia auriculiformis plantation during the wet season. In the Eucalyptus urophylla plantation, both N and P treatments led to SOC reductions (16% and 33%, respectively), but these declines occurred during the dry season. These contrasting seasonal variations across plantations highlight the complexity of SOC dynamics. MBC, an important indicator of soil biological activity [79-81], showed no significant changes across fertilization treatments. TP was consistently higher under the P and NP treatments, indicating an increase in soil P pool, which may potentially alleviate P limitation for plants and microorganisms. Conversely, PDE activity was consistently higher under C and N treatments compared to P and NP treatments. This aligns with previous studies demonstrating that P fertilizer consistently suppresses P-acquiring enzyme activities, possibly by reducing the need for plants and microbes to produce these enzymes when sufficient inorganic P is readily available [82-84]. Although enzymes are known to be affected by factors other than nutrients, our study did not reveal strong effects of fertilization on soil pH. This suggests a potential response by plants and microbes to acquire relatively scarce P resources through the upregulated production of P-acquiring enzymes.

SQI, calculated based on the MDS, did not show any significant difference (p > 0.05), indicating that all the fertilization practices seem to maintain a relatively stable soil nutrient provision. However, analysis of individual soil indicators demonstrated that the nutrient pool has been changed. The observed stability of the SQI across fertilization treatments might be attributed to opposing trends in TP and PDE. While P and NP fertilization significantly increased TP (suggesting an alleviation of P limitation), control plots exhibited the highest levels of PDE activity (potentially indicating a microbial response to acquire scarce P). These contrasting changes in TP and PDE might have offset each other in the SQI calculation, masking the individual effects of these indicators on soil quality. As hypothesized, Acacia auriculiformis plantations generally exhibited higher SQI values compared to Eucalyptus urophylla plantations, although the ANOVA test indicated that the difference was not statistically significant (p > 0.05). This disparity can be attributed to the unique N-fixing capabilities of Acacia species. By establishing a symbiotic relationship with N-fixing bacteria, Acacia auriculiformis can enhance nutrient recycling [85-87]. The SQI in the two plantations became comparable under NP treatment, indicating that NP fertilization may have influenced biological N fixation processes in the Eucalyptus urophylla plantation. The external addition of N and P could diminish the natural advantage of Acacia auriculiformis in enhancing soil quality through N fixation, leading to similar SQI levels between the two plantations [69].

4.3. P Limitation and Soil Quality in Subtropical Plantations

Human activities are causing a continuous increase in N availability within ecosystems, mainly through various anthropogenic inputs. However, this rise is not matched by a corresponding increase in P inputs. This growing imbalance suggests a potential exacerbation of the existing nutrient limitation, potentially impacting ecosystem functioning and services [88-92]. Previous studies have confirmed the pervasive P limitation in subtropical areas [84, 93, 94]. Therefore, nutrient imbalance or limitation, particularly P deficiency, is an important aspect to consider in restoration efforts. Fertilization, as an intensive forest management strategy, is often applied to plantation forests established for productive purposes or ecosystem restoration. Plantation ecosystems rely heavily on their understory vegetation for structural integrity and ecological function.

This understory layer, particularly the dominant herb layer species such as Dicranopteris dichotoma (Thunb.) Bernh., Ottochloa nodosa var. micrantha (Balansa) Keng f. and Blechnopsis orientalis, exhibited sensitivity to fertilization practices (Figure S2). The observed responses in biomass across these representative species highlight the need for careful management of fertilization regimes. This balanced approach is crucial to maximize the potential benefits of fertilization while minimizing any potential negative consequences, such as disruptions to essential ecosystem functions, especially under conditions of P limitation. While our study found no significant changes in SQI with long-term fertilization, analysis of individual indicators, such as TP and PDE which significantly changed in fertilized plots, suggests potential alterations in nutrient cycling processes that warrant further investigation. Therefore, appropriately balanced fertilization regimes and long-term monitoring of soil quality, including analysis of individual indicators alongside SQI, are crucial for sustainable forest restoration.

5. Conclusions

To conclude, SQI can facilitate comparisons between different treatments and sites, offering valuable insight into soil quality and management practices. In this study, the SQI did not reveal significant differences among treatments. However, the indicators in the MDS provided a clear and comprehensive vision of the overall soil nutrient pool. These findings highlight the importance of analyzing individual soil indicators alongside SQI to gain a more comprehensive understanding of soil quality. SQI, while a valuable tool, may not always capture the complexities of nutrient cycling processes due to potential masking effects from opposing trends in MDS indicators. Future research will connect these findings to practical applications for sustainable forest management. This study utilized relatively small experimental plots, which may not fully capture the spatial variability of soil properties. Including results from a larger number of plots in future studies with diverse fertilization regimes would strengthen the generalizability of the findings to other similar plantations. However, despite these limitations, this study provides a valuable framework for assessing soil quality changes under long-term fertilization practices in subtropical plantations.