The study was performed in the Qinling National Forest Ecosystem Research Station (33°18′-33°28′N, 108°21′-108°39′E), which is located in Ningshaan County, Shaanxi Province, China. The area covers 22.25 km² and has a subtropical humid montane climate. The average annual temperature is 10.5°C, and the annual precipitation is 1000 mm, with an average frost-free period of 199 days and a growth period of 177 days (Pang et al., 2021b, 2022). The soil types in the sampling area are classified as Cambisols, Umbrisols, and Podzols according to the Food and Agriculture Organization (FAO) classification system. Across the site, the average soil layer thickness and humus matter thickness are 50 cm and 8 cm, respectively. The forests underwent rotational felling or firewood cutting during the 1960s and 1970s at this research station. After natural regeneration, secondary growth dominated the study area, and the dominant tree species are Quercus aliena var. acutiserrata, Quercus variabilis, Pinus armandii, Betula albosinensis, Picea asperata, and Populus davidiana. Additionally, shrubs including Lonicera tragophylla, Cerasus stipulacea, and Symplocos paniculate, as well as herbs including Lysimachia christinae, Rubus parvifolius, Saussurea mutabilis, and Rubia cordifolia occupy the understory space.

To guide the development of secondary forests towards a healthy stand structure, this study implemented a target tree operation system to transform the existing secondary forests based on the principles of near-natural forest management. Based on field surveys, we selected secondary stands with similar geographical and micro-topographic conditions. The whole-tree harvesting experiment with a complete randomized block design with 5 treatments and 4 replicates, making a total of 20 plots (20 m in length × 20 m in width), were established in the secondary forests in September 2013. At the same time, tree height and diameter at breast height (DBH ≥ 5 cm, measured at 1.3 m) were recorded for each plot. The treatments were 5 levels of thinning intensity: no thinning as a control with (CK), 15% removal of the stand volume (T1), 30% removal of the stand volume (T2), 45% removal of the stand volume (T3), 60% removal of the stand volume (T4). To avoid potential edge effects, a 5-m-wide buffer zone was established around each plot. All harvested materials, including leaves, branches, stems, bark, and twigs, were removed from the plots. We conducted the survey and sampling for this study in 2018, 5 years after the whole-tree harvesting treatment. Detailed setup information about the blocks and plots are presented in Supplementary Figure S1 .

Plant samples were collected in August 2018 during the peak growing season. For trees, the tree height and diameter at breast height (DBH ≥ 5 cm, 1.3 m) in each plot were measured again before sampling, and the trees were also classified and counted by species. Tree leaves, branches, stems, bark, and roots were then sampled by species through pruning, cutting, drilling, and digging, as detailed in our previous study (Pang et al., 2021a). We calculated the tree organ biomass based on the allometric equations developed for trees in or near the study area ( Supplementary Table S1 ). According to the biomass ratios of different species organs, the different organ samples of the trees were mixed evenly. To evaluate the species diversity of shrub and herb, five shrub subplots (2 × 2 m) and five herb subplots (1 × 1 m) were established in each plot, and a whole-plant sampling technique was used to collect shrub and herb biomass. After dividing the shrubs into leaves, branches, and roots, and the herbs into aboveground and belowground parts, then mixed into a uniform sample.

At the plot level, we used a root auger (90 mm diameter) to collect 9 fine root point samples at depths of 0–20 cm, 20–40 cm and 40–60 cm in September and November 2018, and April and June 2019. According to the method described by Brassard et al. (2013), the fine roots were divided into living and dead groups, including three diameter classes (<0.5 mm, 0.5–1mm, and 1–2 mm). For litter sampling, six litter traps (90 cm in diameter) were fixed 1 m above the ground to collect tree litter three times in one year, and all the ground litter was collected in five 1 × 1 m subplots. Litter was categorized into leaves, branches, and miscellaneous material. Eventually, all subsamples of plant and litter were transported to the laboratory and oven-dried at 70 °C to constant weight.

The soil samples were collected from 9 replicate points along an “S” shape in each plot at 0–20 cm, 20–40 cm and 40–60 cm depths by using a 40 mm diameter stainless-steel auger and mixing the samples into one composite sample per plot. The soil samples were then immediately sieved through a mesh < 2 mm, and stones, plant roots, fauna, and debris were removed. Each soil sample was divided into three parts, the first part was used to determine the soil moisture content; the second part was airdried for physicochemical analyses; and the other portion was stored at 4°C for microbial biomass analyses. Soil bulk density samples were obtained randomly from three points per plot by volumetric rings (100 cm³).

All the chemical indicators of plant, litter and soil samples were determined following a previously described method (Bao, 2000). The carbon (C) contents in plant, litter and soil (SOC) samples were analyzed using the K₂Cr₂O₇ oxidation method, while total nitrogen (TN) and total phosphorus (TP) contents were analyzed by colorimetry with an automatic discontinuous elemental analyzer (Clever chem200+, Germany) after wet digestion. Plant, litter and soil total potassium (TK) and soil available potassium (AK) were measured through a flame photometry. Soil available nitrogen (AN) and available phosphorus (AP) were determined by alkaline hydrolysis diffusion and colorimetry respectively. Soil pH was measured in a 1:2.5 soil:water suspension using a glass-electrode meter. Soil moisture content was determined by oven-drying fresh soil at 105°C for 48 hours to a constant weight. Similarly, the volumetric ring soil was dried at 105 °C to constant weight, and the ratio of soil mass to total volume (g·cm⁻³) was calculated to obtain the soil bulk density (De Vos et al., 2005). Microbial biomass C (MBC), N (MBN) and P (MBP) were quantified with a total organic C analyzer (TOC-LCSN, Shimadzu Co., Japan), flow analyzer (AutoAnalyzer 3, Seal Analytical Ltd., UK) and spectrophotometer (UV-1900, Shimadzu Co., Japan) after fumigation (Kim et al., 2018).

The biomass and productivity of tree and ecosystem, and ground litter biomass were significantly reduced with increasing thinning intensity, whereas shrub and herb biomass and productivity showed the opposite trend ( Supplementary Figures S2–S9 , p < 0.05). Although litter miscellaneous productivity was significantly higher in T2 and T3 than in other thinning treatments ( Supplementary Figure S8A , p < 0.05), total litter productivity did not respond significantly to the thinning treatments ( Supplementary Figure S8B , p > 0.05).

Overall, thinning significantly reduced tree C, N, P and K stocks and accumulation rates, ground litter C, N, P and K stocks, soil C, N and P stocks, ecosystem C, N and P stocks and accumulation rates and ecosystem K accumulation rate, and significantly increased the ecosystem K stock and the stocks and accumulation rates of shrub and herb C and N, shrub P and herb K ( Supplementary Table S3 and Supplementary Figures S10–S16 , p < 0.05). The shrub K and herb P stocks and accumulation rates were nonsignificant among the different thinning intensities ( Supplementary Table S3 and Supplementary Figures S11G, H , S15G, H , p > 0.05).

As shown in Supplementary Figure S17 , the tree litter C concentration was highest in winter, while N, P, and K concentrations peaked during the spring-summer period while the return amounts of tree litter C, N, P, and K were significantly higher in autumn than in other seasons (p < 0.05). However, total aboveground litter C, N, P, and K return amounts did not respond significantly to thinning treatments ( Supplementary Figure S18 , p > 0.05).

Although the <0.5mm fine root in T2 and T4 exhibited significantly higher productivity, turnover rate, and C, N, P, and K return amounts than in other treatments ( Supplementary Figures S19 , S20A–D , p < 0.05), no significant differences were observed in total (<2mm) fine root C, N, P, and K return amounts across thinning intensities ( Supplementary Figures S20E-H , p > 0.05).

Likewise, thinning did not significantly impact total plant C, N, P, and K return amounts at the ecosystem level ( Supplementary Figure S21 , p > 0.05).

Significant differences were observed among ecosystem N, P, and K cycling parameters, which showed varied responses to thinning treatments ( Figure 1 ). The ecosystem N total return amounts, absorption amounts, absorption coefficient, and utilization coefficient were the highest ( Figures 1A–D , p < 0.05), while the N surplus amounts and surplus rate were the lowest ( Figures 1F, J , p < 0.05). The ecosystem P productivity, retention coefficient, and plant nutrient turnover time were at their highest levels ( Figures 1G, I, K , p < 0.05). The ecosystem N and K cycle coefficient, recycling coefficient and plant nutrient accumulation coefficient were comparable and significantly higher than those parameters of P ( Figures 1E, H, L , p < 0.05).

The ecosystem N, P, and K total return amounts and productivity had no significant responses to the thinning treatments ( Figures 1A, G , p > 0.05). The N, P, and K absorption amounts, retention coefficients, and plant turnover nutrient time significantly decreased as thinning intensities increased ( Figures 1B, I, K , p < 0.05), while the N, P and K cycle coefficient, recycling coefficient, surplus amounts, and plant nutrient accumulation coefficient significantly increased ( Figures 1E, H, J, L , p < 0.05). Thinning significantly increased the ecosystem N absorption coefficient and utilization coefficient, and showed opposite effects on the P and K absorption coefficients and utilization coefficients ( Figures 1C, D , p < 0.05). The N surplus rate in T2 was significantly lower than other treatments, while there was a significant increasing trend in ecosystem P and K surplus rates following thinning intensities ( Figure 1F , p < 0.05).

Compared with CK, thinning improved the ecosystem N, P, and K cycling function indexes by 5%~232%, 32%~195%, and 104%~233% respectively ( Figure 2 , p < 0.05).

The ecosystem N, P, and K cycling function indexes were significantly negatively correlated with SOC, TN, TP, AN, AP, tree productivity, tree ecosystem level K, and herb ecosystem level N, while there were significant positive correlations between those nutrient cycling function indexes and the shrub and herb productivity, and the herb ecosystem level C:N ( Supplementary Tables S2–S4 and Figure 3A , p < 0.05). The ecosystem N and P cycling function indexes were significantly negatively related to soil C:P, N:P ratios and ground litter biomass ( Supplementary Tables S2–S4 and Figure 3A , p < 0.05). The tree ecosystem level P had negative correlations with the ecosystem P and K cycling function indexes ( Supplementary Tables S2–S4 and Figure 3A , p < 0.05). The ecosystem K cycling function index was positively correlated with the tree ecosystem level C, tree ecosystem level C:N, tree ecosystem level C:P and the herb ecosystem level C:P, but had the opposite effect on the herb ecosystem level P ( Supplementary Tables S2–S4 and Figure 3A , p < 0.05). The optimal predictors for ecosystem N, P, and K cycling function indexes, as identified by the Boruta method, aligned with with Pearson’s significantly correlated results ( Supplementary Tables S2–S4 and Figures 3A–D , p < 0.05).

Based on the classification of selected predictors, the VPA analyses indicated that these soil nutrient, litter biomass, plant productivity and plant nutrient collectively explained the majority of variation associated with the ecosystem N, P, and K cycling functions ( Figure 4 ). The soil nutrient, litter biomass, plant productivity and plant nutrient accounted for 67.1%, 37%, 77.2%, 27.5% of the variation of N cycling functions ( Figure 4A ), 71.3%, 33.5%, 79.4%, 28.4% of the variation of P cycling functions ( Figure 4B ), and 58.3%, 26.6%, 58.8%, 87.5% of the variation of K cycling functions ( Figure 4C ), respectively.

PLS-PM identified direct and indirect effects of forest thinning, plant productivity, litter biomass, soil nutrient and plant nutrient on the ecosystem N, P, and K cycling functions ( Figures 5 – 7 ). For N, P, and K cycling functions, thinning positively affected the plant productivity (0.84 of the direct effects) and negatively affected litter biomass (-0.76 of the direct effects) ( Figures 5 - 7A ). For N, and P cycling functions, plant productivity further negatively regulated plant nutrient and soil nutrient (-0.51 and -0.53 of the direct effects), and litter biomass positively regulated soil nutrient (0.33 of the direct effects) ( Figures 5 , 6A ). Plant and soil nutrients jointly exerted a negative direct effect on the ecosystem N (-0.28 and -0.74) and P (-0.35 and -0.68) cycling function ( Figures 5 , 6A ). For K cycling function, plant productivity directly positively regulated plant nutrient (0.51) and negatively regulated soil nutrient (-0.64) ( Figure 7A ). And plant nutrient and soil nutrient have a direct positive effect (0.56) and a direct negative effect (-0.52) on the ecosystem K cycling function ( Figure 7A ). Overall, the total positive effects of forest thinning on ecosystem N, P, and K cycling functions were 0.639, 0.626, and 0.585 respectively ( Figures 5 - 7B ).

Nutrient cycling is a key ecosystem process closely related to forest ecosystem structure and function and is also the basis for supporting organic matter production in forest ecosystems (Vrignon-Brenas et al., 2019; Perron et al., 2021). Although ecosystem nutrient cycling is crucial to the healthy and stable development of communities, there have been few evaluations of nutrient cycling characteristics at the ecosystem level. In our study, the ecosystem N, P, and K cycling parameters were significantly different from each other, which was consistent with previous research results on the key cycling processes of different nutrient elements (absorption, accumulation, distribution, utilization efficiency, etc.) (Rodríguez-Soalleiro et al., 2018; Schlesinger, 2021). A possible explanation is that the cycling process of different nutrient elements will be corporately affected by the physical and chemical properties of the elements, biological effects, environmental conditions (Vitousek, 2018; Johnson and Turner, 2019; Sayer et al., 2020), exhibiting different nutrient cycling characteristics. Specifically, the ecosystem’s high N return, absorption, and utilization rates, along with low N surplus, indicate severe nitrogen limitation for plant growth in this area, which conformed with our previous results regarding nutrient stoichiometry and nutrient reabsorption in multi-plant organs (Pang et al., 2021b). Additionally, compared with N and K, the ecosystem’s higher P productivity, retention coefficient, and plant nutrient turnover time and lower P cycle coefficient, recycling coefficient, and plant nutrient accumulation coefficient implied that although P productivity is high in this study area, the cycling intensity is the lowest. On the one hand, ecosystem P is mainly derived from the weathering of soil phosphate minerals, which is a relatively slow process (Zhou et al., 2018; Koester et al., 2021). On the other hand, ecosystem P cycle also depends on the amount of P returned from plant organisms (Whitehead, 2000; Mengel and Kirkby, 2012), but in this study, plant P showed higher retention and higher turnover time, which also led to lower ecosystem P cycling intensity. As for the N cycle, ecosystem N sources are diversified. In addition to soil N nutrient mineralization and plant return, N can also be input through atmospheric nitrogen deposition and biological nitrogen fixation, so it has a higher cycling intensity relative to P (Wang et al., 2007; Nevison et al., 2022). Unlike N and P, which often form complex molecular structures, K mainly exists as soluble ions. These ions are more easily leached, resulting in a faster nutrient cycle (Qiao et al., 2023). Moreover, in order to maintain the K balance, plants in the forest ecosystem have a “pumping” mechanism, which further promotes the K circulation rate in the ecosystem (Sardans and Peñuelas, 2015; Schlesinger, 2021). Overall, ecosystem N, P, and K cycles exhibited heterogeneous patterns, supporting our hypothesis 1.

To maintain the forest ecosystem sustainable development, in addition to focusing on forest productivity sustainability, the importance of nutrient sustainability cannot be ignored (Haberl et al., 2010; Rodríguez-Soalleiro et al., 2018). It is well-known that the whole-tree harvesting method permanently removes a large amount of wood and plant residues from the nutrient cycling system of plant growth, and this high nutrient output rate has a negative impact on the ecosystem nutrient budget (Smolander et al., 2010; Thiffault et al., 2011; Pang et al., 2021b). Therefore, understanding the effects of whole-tree harvesting on the forest ecosystem nutrient cycling characteristics and functions at the ecosystem level is of great significance for explaining the forest ecosystem nutrient balance mechanism after thinning and formulating scientific forest management strategies (Jordan, 1985; Pei et al., 2017). Contrary to our initial hypothesis 2, whole-tree harvesting can significantly enhance the ecosystem N, P, and K cycling functions. The results of this study were inconsistent with previous researchers who suggested that whole-tree harvesting aggravated nutrient loss and had a negative impact on nutrient cycling in forest ecosystems (Kaarakka et al., 2014). From the details, there was no difference in the nutrient concentration and returned biomass of above-ground litter and underground fine roots under different thinning intensities ( Supplementary Figures S8 , S17 , S19 ; Supplementary Table S3 ), and the nutrient input from atmospheric precipitation was also consistent. So, the total nutrient return amounts and nutrient productivity did not respond to thinning treatments, and had no effects on the ecosystem nutrient cycling functions ( Figure 1 ). However, N, P, and K cycle coefficients, recycling coefficients, surplus rates, and plant nutrient accumulation coefficients significantly increased with the increasing thinning intensity, while N, P, and K absorption amounts, retention coefficients, and plant nutrient turnover time significantly decreased ( Figure 1 ). Collectively, these results proved that the N, P and K cycling intensity in the ecosystem increased following thinning, promoting the improvement of N, P and K cycling functions at the ecosystem level. Simultaneously, there are coupling effects and mutual influences in the multi-nutrient element cycling process (Zhang et al., 2019). In this study area, due to N limitation, the ecosystem N absorption coefficient and utilization coefficient increased significantly with the increased thinning intensity, and the N surplus amounts showed a significantly low value at T2, indicating that the ecosystem N utilization and consumption increased significantly, the N cycling speed was accelerated, and which may also have led to the improvement of the ecosystem P and K cycling functions. Furthermore, previous studies on the effects of whole-tree harvesting on nutrient cycling in forest ecosystems mostly focused on the nutrient budget of plant or soil pools (the single ecosystem component) (Yan et al., 2017; Pang et al., 2021b), and lacked understanding of the nutrient cycling function at the ecosystem level. In this study, the ecosystem function index was employed to systematically extract comprehensive information regarding ecosystem nutrient cycling characteristics, providing in a more representative assessment of the impact of whole-tree harvesting on nutrient cycling functions at the ecosystem level.

Nutrient cycling in forest ecosystems involves complex processes driven by various influencing factors (Johnson and Lindberg, 2013; Vitousek, 2018). Whole-tree harvesting drastically changes the biotic and abiotic factors in the forest ecosystem, thereby most likely significantly affecting the ecosystem nutrient cycling function (McDaniel et al., 2013; Sicard et al., 2016). In the present study, PLS-PM revealed that whole-tree harvesting mainly promoted ecosystem N and P cycling functions through two pathways: (1) whole-tree harvesting regulated litter biomass and affected soil nutrient characteristics, then regulating the ecosystem N and P cycling functions; (2) whole-tree harvesting adjusted ecosystem plant productivity, influenced both plant and soil nutrient characteristics, and then adjusting ecosystem N and P cycling functions. First of all, early studies pointed out that litter nutrient return is an important source of soil nutrients (Chao et al., 2019; Keller and Phillips, 2019), and litter biomass has a positive impact on soil nutrients ( Figures 5 , 6 ). Our results showed that whole-tree harvesting significantly reduced litter biomass return and had a negative impact on soil nutrients ( Supplementary Figure S3 ; Supplementary Tables S2 , S3 ), in agreement with previous results (Clarke et al., 2021; Smith et al., 2022; Vos et al., 2023). Secondly, the ecosystem plant productivity represents the ability of plant systems to compete and consume nutrients (Liu, 2009). We found that although the productivity of the tree layer was significantly reduced by whole-tree harvesting, the productivity of the understory vegetation increased rapidly after thinning ( Supplementary Figure S6 ). At the same time, the VPA analysis results showed that the ecosystem plant productivity was the main controlling factor of the ecosystem N and P cycling functions while the PLS-PM path coefficient showed that the whole-tree harvesting had a significant positive effect on the ecosystem plant productivity. All these results indicated that the productivity of understory vegetation after thinning played a more significant role in regulating the ecosystem N and P cycling functions and confirmed findings from previous studies that thinning promotes understory vegetation growth, with its productivity making an essential contribution to ecosystem nutrient reserves and nutrient cycling (Ma et al., 2007). Therefore, when whole-tree harvesting caused a large amount of nutrient loss, the rapidly increased understory vegetation productivity intensified nutrient consumption and competition, having negative effects on plant and soil nutrient characteristics (Elser et al., 1996; Bauhus et al., 2001). Ultimately, the degradation of plant and soil nutrient characteristics had direct negative impacts on the ecosystem N and P cycling functions, supporting the previous conclusion that intensity disturbance negatively affects the nutrient budget and nutrient cycling function (Kaarakka et al., 2014).

Compared with the ecosystem N and P cycling functions, PLS-PM results showed that whole-tree harvesting mainly affected the plant and soil nutrient characteristics by adjusting the ecosystem plant productivity, thereby enhancing the ecosystem K cycling function. Moreover, within the regulation pathway, it was noteworthy that the positive effect of ecosystem plant productivity on plant nutrient characteristics, as well as the positive effect of plant nutrient characteristics on K cycling function, differed from the effects of plant productivity on N and P cycling functions, where it exerts an indirect positive effect through its negative impact on plant nutrient characteristics. It is well known that removing trees through thinning can reduce competition pressure for nutrient resources in plant communities (Hertel et al., 2009). Relative to the important regulatory role of herb nutrients (plant nutrient characteristics) on the ecosystem N and P cycling functions, tree nutrients were included in the plant nutrient characteristics in the K cycling function PLS-PM model ( Figures 5 - 7C ). Therefore, in the K-cycling process, the direct positive effects of ecosystem plant productivity on plant nutrient characteristics (herbs and trees) may be due to the fact that the nutrient competition pressure released by thinning exceeded the nutrient consumption and competition pressure caused by increased understory vegetation productivity, thus inducing the positive effects of ecosystem plant productivity on plant nutrient characteristics. Additionally, the differential regulation patterns of N, P, and K cycling functions may depend on synergistic effects of multiple factors not monitored in this study, such as climate drivers, elemental heterogeneity or biocatalytic effects (Vitousek, 2018; Johnson and Turner, 2019; Sayer et al., 2020). And unexplained parts of the VPA and PLS-PM models require further investigation ( Figures 6 , 7 ). Furthermore, the direct positive effect of plant nutrient characteristics on K cycling function was consistent with the positive correlation between plant nutrient stoichiometric ratios and K cycling function index ( Figures 3A , 7A ), suggesting the importance of plant nutrient stoichiometric ratios in the ecosystem K cycling function regulation. Similar findings were also reported elsewhere in studies of forest or non-forest ecosystem nutrient cycling (Tripler et al., 2006; Chen B. et al., 2023). Overall, whole-tree harvesting showed inconsistent regulatory mechanisms for N, P and K cycling functions, supporting our hypothesis 3.