Thinning enhances ecosystem multifunctionality via increase of functional diversity in a Pinus yunnanensis natural secondary forest


 Background: The impacts of thinning on ecosystem multifunctionality (EMF) remain largely unexplored. In this study, we analyzed nine variables related to four ecosystem functions (nutrient cycling, soil carbon stocks, decomposition, and wood production) under five thinning intensities. We included a control group to evaluate the shift in EMF of a Pinus yunnanensis natural secondary forest. We also assessed the relationship between above- and belowground biodiversity and EMF under these different thinning intensities. Additionally, we evaluated the effects of biotic and abiotic factors on EMF with the structural equation model (SEM). Results: We found that EMF tended to increase with thinning intensity, and that thinning significantly improved EMF except the low intensity of thinning. Individual ecosystem functions (EFs) all had a significant positive correlation with thinning intensity. Different EFs showed different patterns with the increase of thinning intensity: the nutrient cycling and the soil carbon stock of thinning three times and five times were significantly greater than other thinning intensities and control group; decomposition correlated directly to the increase of thinning intensity; the wood production of the fourth thinning was greatest. Thinning intensity had a significant positive correlation with functional diversity and soil moisture. Both functional diversity and soil moisture had a significant positive correlation with EMF, but soil fungal operational taxonomic units (OTUs) had a significant negative correlation with EMF. Based on SEM, we found that thinning improved EMF mainly by increasing functional diversity. Conclusion: Our study both demonstrates that thinning is a good management technique from an EMF perspective, and provides an input to improve management of a P. yunnanensis natural secondary forest.


Background
Thinning is a broadly used and effective method in forest management that could significantly change ecosystem function by removing certain trees and reducing competition (Baena et al. 2013;Smolander et al. 2015;Trentini et al. 2017). Compared to traditional commercial logging, thinning is not only considered for timber production, but also for ecological restoration, reduction of fire risk, biodiversity conservation, the creation and maintenance of wildlife habitat, and so on (Keyser et al. 2012). Some studies indicate that thinning could promote productivity and soil nutrients (Dang et al. 2018), impact soil microbial activity and biomass (Baena et al. 2013), and increase decomposition of cellulose substrate in the organic horizon (Thibodeau et al. 2000). The simultaneous occurrence of these ecosystem functions (EFs) is defined as ecosystem multifunctionality ('EMF') (Hector and Bagchi 2007). We know little about the effect of thinning on multiple EFs and EMF. Because of the functional uniqueness of the species, biodiversity was crucial for EMF, and studies on the relationship between biodiversity and EMF ('BEMF') has become a hotspot in ecology (Zavaleta et al. 2010;Gamfeldt and Roger 2017).
Biodiversity consists of three components: taxonomic diversity, functional diversity and phylogenetic diversity (Bagousse-Pinguet et al. 2019;Luo et al. 2019). Most studies have focused on taxonomic diversity, e.g. species richness as a metric of biodiversity, because it was simple and could be extended to other diversity metrics (Isbell et al. 2011;Maestre et al. 2012;Gamfeldt and Roger 2017).
Scientists pointed out that species richness was not a particularly good metric for biodiversity because species richness does not provide direct information about traits like functional or phylogenetic diversity does (Gamfeldt and Roger 2017;Zirbel et al. 2019). Functional diversity might be a better predictor of ecosystem function or multifunctionality (Díaz and Cabido 2001;Gross et al. 2017). Functional diversity is based on functional traits, which are closely related to resource utilization (Díaz and Cabido 2001). Phylogenetic diversity could capture information on unmeasured traits as they relate to ecosystem function, so it might be a key component to predict ecosystem function (Flynn et al. 2011;Srivastava et al. 2012). The effects of different biodiversity attributes on EMF might vary in magnitude (Bagousse-Pinguet et al. 2019). Except aboveground plant biodiversity, a growing body of studies suggests that soil microbial diversity also plays a critical role in maintaining EMF (Jing et al. 2015;Delgado-Baquerizo et al. 2016;Luo et al. 2018). Many studies have shown that the loss of soil microbial diversity, including bacterial and fungal diversity, would threaten EMF (Wagg et al. 2014;Jing et al. 2015;Delgado-Baquerizo et al. 2016;Delgado-Baquerizo et al. 2017). The effect of multiple plant diversity attributes and soil microbial diversity on EMF under different thinning 4 intensities is unclear.
To fill these knowledge gaps, we quantified four key ecosystem functions: nutrient cycling (soil total nitrogen, soil hydrolysable nitrogen, soil total phosphorus and soil available phosphorus), soil carbon stocks (soil total carbon), decomposition (β-glucosidase, urease and phosphatase activity) and wood production (aboveground woody biomass). Using these EFs, under five different thinning intensities in a Pinus yunnanensis natural secondary forest, we analyzed how aboveground plant biodiversity (taxonomic, functional and phylogenetic diversity), soil microbial biodiversity (bacterial and fungal diversity), and abiotic factors (soil moisture and pH) influenced EMF. We predicted that: (1) thinning could either directly or indirectly improve EMF; (2) the effects of functional and phylogenetic diversity were greater than taxonomic diversity; (3) soil microbial diversity could improve EMF. The results of this study could provide new insight for forest management and conservation.

Study sites and sampling
This study was conducted at the Yongren National Forest Farm (26°12´~26°2´N; 101°27´~101°37′E), located in the Jinsha River watershed of Yunnan Province, southwest China (Fig. 1). The site has a mean annual temperature of 17.8℃, and a mean annual precipitation of 840 mm, with an altitudinal range from 1, 800 m to 2, 500 m. The soil type is hilly red soil. The Yongren National Forest Farm has been producing timber since the 1950s. We obtained field data from 52 plots in April 2018. The plots were set as Li et al. (2020) described: thinning one time (T1, 10 plots); thinning two times(T2, 10 plots); thinning three times (T3, 10 plots); thinning four times (T4, 10 plots); thinning five times (T5, 9 plots); control group without thinning (CG, 3 plots). These plots (20m × 30m) were representative of the major vegetation types and management styles in the area. The thinning began in the 1970s, and occurred 10-15% of the interim total time. Selected species included P. yunnanensis and broadleaf species. Within the plots, we identified species and used a diameter tape and tree altimeter to measure diameter at breast height (DBH, cm) and height (H, m) for all individuals > 1cm DBH.
Following the standardized protocols of Cornelissen et al. (2003), we gathered 50 to 100 sun-exposed mature leaves from five to ten individuals of each species in the middle or upper section of the tree 5 per plot. We labeled the leaves by species and put them into the paper bags. For broad-leaf species, we used a laser area meter (LI-COR 3100C Area Meter, LI-COR, USA). For conifer species, such as P.
yunnanensis, we used a method described by Zhang et al. (2008) that saw the leaf of P. yunnanensis as a cylinder. The leaf area of P. yunnanensis was calculated by the equation shown below: LA is the leaf area of P. yunnanensis, d is the leaf diameter of P. yunnanensis and L is the leaf length of P. yunnanensis. The d and L were measured by a Vernier caliper (precision: 0.02 mm) and a ruler (precision: 1 mm), respectively.
The fresh weight of each leaf was measured immediately with an electronic balance (precision: 0.0001 g). We weighed them again after drying for 72 h at 60 °C. Afterward, we used a ball mill (NM200, Retsch, Haan, Germany) to grind all leaf samples into a fine powder to prepare for the measurement of leaf N concentration and leaf P concentration.
In each plot, we collected the wooden core for large trees and shrubs (DBH > 5cm) by increment borer. Then, the volume was measured using the diameter and the length of the wooden core. For other small shrubs (DBH ≤ 5cm), we collected 5-10 branches and removed the bark, then measured volume using the water replacement method. The wooden cores and the branches were weighed using an electronic balance (precision: 0.0001 g) after drying for 72 h at 60 °C.
We also collected 5 soil samples (0-20 cm) from each plot and mixed them evenly. Those samples were divided into three sections and placed in plastic bags. The first portion was stored in a portable refrigerator at 4 °C to maintain freshness and moistness for soil enzyme activity. The second portion was stored at -20 °C to allow for DNA extractions to assess soil bacterial and fungal diversity. When we got back to the lab, the third portion was sieved by a 2mm mesh and air-dried for one month for physiochemical analyses.

Assessing aboveground diversity
Species richness was calculated for taxonomic diversity. As a classical measure of functional diversity, Rao's quadratic entropy (RaoQ) could reflect both functional richness and functional divergence and has been applied to the study of the relationship between biodiversity and ecosystem function (Botta-6 Dukát 2005; Wasof et al. 2018). RaoQ was calculated using specific leaf area (SLA = leaf area / dry weight, mm 2 ·mg -1 ), leaf dry mass content (LDMC = leaf dry weight / leaf saturated fresh weight, %), wood density (WD = wood dry weight / wood volume, g·cm -3 ), leaf N concentration (LNC, mg·g -1 ) and leaf P concentration (LPC, mg·g -1 ). Faith's PD was calculated as a proxy of phylogenetic diversity (Faith 1992).

Assessing belowground diversity
The number of soil bacterial and fungal operational taxonomic units (OTUs) was used to assess belowground diversity. According to manufacturer's instructions, we used the PowersoilRDNA Isolation and FITS7/ITS4 primer sets, respectively (Ihrmark et al. 2012). We defined an OTU as having similarity above 97% (Öpik et al. 2009).

Measurement of physiochemical data
Soil total nitrogen, leaf N concentration and soil total carbon were determined by combustion in a CHN elemental analyser (2400Ⅱ CHN elemental analyser, PerkinElmer, Boston, MA, USA). Soil hydrolysable nitrogen was determined through alkaline hydrolysis diffusion. We used a molybdenum antimony blue colorimetry method to analyze soil total phosphorus and leaf phosphorus concentration. We used the Olsen method to analyze soil available phosphorus. For the three activities of enzymes, β-glucosidase activity was analyzed as described by Tabatabai (1982), urease activity was analyzed as described by Nannipieri et al. (1980), and phosphatase activity was analyzed as described by Tabatabai and Bremner (1969). For soil bacterial and fungal biodiversity, we used a suite of molecular techniques as described by Jing et al. (2015). We used a pH meter (Thermo Orion-868) to analyze soil pH in a 1:5 ratio of fresh soil to deionized water slurry. The soil moisture was calculated by the ratio of soil dry weight to soil fresh weight, and was measured by an electronic 7 balance (precision: 0.0001 g).

Individual ecosystem functions
We used an approach described by Lucas-Borja and Delgado-Baquerizo (2019) in which EMF was assessed by nutrient cycling, soil carbon stocks, decomposition and wood production (Table 1). Soil total nitrogen, soil hydrolysable nitrogen, soil total phosphorus and soil available phosphorus were used as the proxy of nutrient cycling. The content of soil total carbon was used as the proxy of soil carbon stocks. The three enzyme activities (β-glucosidase, urease, phosphatase) were used as the proxy of decomposition. Finally, the aboveground woody plant biomass was used as the proxy of wood production. We used the general allometric equation to estimate the aboveground biomass (AGB, kg/hm 2 ) of individual trees. For some species without allometric equation, we chose the equation of a similar species instead (Paul et al., 2014). We chose the allometric equation as Ali et al. (2015) described to calculate the biomass of the less frequently occurring species. The total AGB per plot was calculated as the sum of the AGB of all trees. The allometric equations for suitable species are presented in Additional file: Table S1. All variables were standardized using the Z-score transformation (Bagousse-Pinguet et al. 2019). Then, these four components of EMF were calculated as the average of the standardized values across the corresponding variables.

Assessing ecosystem multifunctionality
We used the averaging approach (Maestre et al. 2012) to obtain a multifunctionality index. Firstly, we calculated the standardized Z-score value for each of the nine variables per plot. Then the EMF index was obtained by averaging the Z-score value of each of the nine ecosystem variables across plots. We also used this approach to calculate the index of the sets of functions related to nutrient cycling, soil carbon stocks, decomposition and wood production. The averaging approach has been widely used in the study of BEMF (Valencia-Gómez et al. 2015;Berdugo et al. 2017).

Statistical analyses
Firstly, we used "LSD" test in one-way analysis of variance (ANOVA) to test the difference of multiple All statistical analyses were conducted in R version 3.3.2. RaoQ was calculated in the "FD" package (Laliberté et al. 2014). For phylogenetic diversity, we first built the list of species found in our field survey, and then we used Phylomatic software (Webb et al. 2008) to create a phylogeny. The Phylomatic software used the Zanne et al. (2014) tree, including a dated molecular tree of > 32, 000 angiosperm species. Then we used the Phylocom software to calculate Faith's PD (Webb et al. 2008).
The SEM was conducted in the "lavaan" package.

Effect of thinning intensity on EFs
Thinning had different effects on different EFs (Table 2). Nutrient cycling and soil carbon stocks responded similarly to increasing thinning intensity, appearing "N" shaped. The values of the T3 of nutrient cycling and soil carbon stocks were significantly greater than other thinning and control groups (P < 0.05), except for the T5. Decomposition responded directly to increasing thinning intensity, with the values of the T4 and T5 being significantly greater than the control group (P < 0.05), and the T5 showing the highest value. Wood production displayed an "M" pattern with increasing thinning intensity, and the value of the T4 was significantly greater than other thinning and control groups (P < 0.05), except for the T2. Different lower case letter indicating significant differences at P < 0.05.

EMF under different thinning intensities
We found that, excepting the low levels of thinning (one times and two times), other thinning intensities could significantly improve EMF (p < 0.05) (Fig. 3). Except for T2, there was a significant difference in EMF between T1 and other thinning intensities (p < 0.05) (Fig. 3).

Relationships between thinning intensity and above-, belowground biodiversity and abiotic factors
There was a significant positive correlation between thinning intensity and functional diversity (p < 0.01, r = 0.369) ( Table 3). Thinning intensity had a very significant positive correlation with soil moisture (p < 0.01, r = 0.524) (Table 3). However, thinning intensity had no significant correlation with taxonomic diversity, phylogenetic diversity, soil fungal OTUs, soil bacterial OTUs or soil pH (p > 0.05) ( Table 3).

Effect of abiotic factors on EMF and multiple EFs
Soil moisture was significantly positively correlated with EMF (P < 0.01), and explained 13.2% of the variation in EMF (Fig. 5a). There was no significant correlation between soil pH and EMF (P > 0.05) (Fig. 5b). Soil moisture was significantly positively correlated with nutrient cycling (P < 0.05), soil carbon stocks (P < 0.05) and wood production (P < 0.05) (Additional file: Figure S2a, c, g). Soil pH had no significant effect on multiple EFs (P > 0.05) (Additional file: Figure S2b, d, f, h).

Effect of thinning intensity, biotic and abiotic factors as predictors of EMF
Thinning might directly and indirectly through functional diversity impact EMF (P < 0.05) (Fig. 6a). The variations of EMF, functional diversity, soil moisture and soil fungal OTUs explained by this model were 42.6%, 13.7%, 27.5% and 5.7%, respectively (Fig. 6a). Soil fungal OTUs could also directly impact EMF (P < 0.05), but there was no significant relationship between thinning and soil fungal OTUs (P > 0.05) (Fig. 6a). The strongest relationship in the SEM analysis was between thinning and soil moisture, but we also found a weaker positive relationship between soil moisture and EMF (P > 0.05) (Fig. 6a). The total effects of functional diversity, soil moisture, soil fungal OTUs and Thinning intensity were 0.292, 0.364, -0.192 and 0.558, respectively.

Discussion
Our study explored the relationship between above-and belowground biodiversity, abiotic factors and EMF under different thinning intensities. The results basically met our expectations: thinning could significantly improve EMF except the low levels of thinning. Thinning intensity was significantly positively correlated with nutrient cycling, soil carbon stocks, decomposition and wood production.
Different EFs showed different response patterns to thinning intensity. Among biodiversity elements, only functional diversity exhibited a significantly positive correlation with thinning intensity. Similarly, among abiotic factors, only soil moisture displayed a significantly positive correlation with thinning intensity. Based on SEM, we found that thinning could not only directly but also indirectly impact EMF through functional diversity. Contrary to our predictions, the effect of functional diversity, not phylogenetic diversity, on EMF was stronger compared to taxonomic diversity. Neither phylogenetic diversity nor taxonomic diversity had a significant effect on EMF. In addition, the relationship between belowground biodiversity and EMF was also contrary to our predictions: soil microbial diversity was negatively correlated with EMF, and even soil fungal OTUs and bacterial OTUs differed in significance.
Thinning is a widespread forest management technique that could alter forest microclimate. Studies show that soil temperature and moisture increased as canopy density and tree density decreased after thinning (Masyagina et al. 2006). Thinning could also alter soil conditions and nutrient contents (Smolander et al. 2013;Kim et al. 2015). Soil microclimate, soil conditions and nutrient contents could affect the enzyme activities related to decomposition (Adamczyk et al. 2015). In our study, an increase in thinning intensity promoted multiple EFs. Consistent with previous studies (Hwang and Son 2006;Kim et al. 2016), we observed that soil carbon stocks increased with thinning intensity.
Changes in soil carbon stock were caused by root mortality and the amount of felled woody vegetation entering the soil at different thinning intensities (Kim et al. 2016). Wang et al. (2018) indicated that the amount of upper canopy removed would affect decomposition rates. Thinning reduced canopy density and allowed more sunlight and precipitation to reach the ground (Thibodeau et al. 2000;Simonin et al. 2007). In addition, plant residues entering the soil could also increase decomposition (Adamczyk et al. 2015). The increase of soil moisture caused by thinning is one of the main explanations for the increase in decomposition. In our study, thinning intensity led to an increase in nutrient cycling, a relationship caused by: the nutrient pool formed by thinning residues entering the soil; the loss of individuals, reducing competition for nutrients; the increase of leaching loss (Carlyle 2011). Nutrient cycling implies a positive feedback between thinning intensity and decomposition (Kim et al. 2016). Wood production increased with thinning intensities because thinning promotes the radial growth of remaining trees due to the weakening of competition and the removal of water restriction (Sohn et al. 2016;Carbon et al. 2018). The positive effect of thinning on these four EFs leads to the direct positive effect of thinning on EMF.
There was no significant correlation between thinning intensity and taxonomic diversity. This result was consistent with the results of many previous studies (Ares et al. 2010;Nunes et al. 2014).
Thinning might change the structure of the community but not the number of species. In the field, we did not find any invasive species. Phylogenetic diversity, the strongest expected predictor of EMF, was not significantly related to EMF, perhaps because the phylogenetic signal of functional traits related to each EF might not always be present (Davies et al. 2016). Zirbel et al. (2019) suggested that if the functional traits that determined an EF were not phylogenetically conserved, phylogenetic diversity might not be a good predictor of that EF. We suspected that some species, functional traits, and phylogenetic signals were lost during thinning. In addition, phylogenetic diversity is a measure of biodiversity that combines species richness and evolutionary history (Cadotte 2013

Consent for publication
Not applicable.

Figure 1
The location of the study at Yongren National Forest Farm (black squares, n=52) 22 Figure 2 Relationship between thinning and multiple EFs. Multiple EFs including nutrient cycling (a), soil carbon stocks (b), decomposition (c), wood production (d).
23 Figure 3 EMF index with different thinning intensities. Different lower case letter represents significant difference at p < 0.05 level. CG: control group; T1: thinning one time; T2: thinning two times; T3: thinning three times; T4: thinning four times; T5: thinning five times.
25 Figure 5 Relationship between abiotic factors and EMF. Abiotic factors include soil moisture (a), soil pH (b).

Figure 6
Structural equation models (SEM) of the thinning intensity, biotic and abiotic factors as predictors of EMF. Solid red arrows represent significant positive effects (P < 0.05), solid black arrows represent significant negative effects (P < 0.05), and dotted grey arrows represent non-significant paths (P 0.05). Path standardized coefficients (β) are on the arrows.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.
supplementary information.doc