Selective logging enhances ecosystem multifunctionality via increase of functional diversity in a Pinus yunnanensis forest in Southwest China

Background: The impacts of selective logging on ecosystem multifunctionality (EMF) remain largely unexplored. In this study, we analyzed the response of nine variables related to four ecosystem functions (i.e. nutrient cycling, soil carbon stocks, decomposition, and wood production) to ve selective logging intensities in a Pinus yunnanensis-dominated forest. We included a control group with no harvest to evaluate the potential shifts in EMF of the P. yunnanensis forests. We also assessed the relationship between above- and belowground biodiversity and EMF under these different selective logging intensities. Additionally, we evaluated the effects of biotic and abiotic factors on EMF using a structural equation modeling (SEM) approach. Results: Individual ecosystem functions (EFs) all had a signicant positive correlation with selective logging intensity. Different EFs showed different patterns with the increase of selective logging intensity. We found that EMF tended to increase with logging intensity, and that EMF signicantly improved when the stand was harvested at least twice. Both functional diversity and soil moisture had a signicant positive correlation with EMF, but soil fungal operational taxonomic units (OTUs) had a signicant negative correlation with EMF. Based on SEM, we found that selective logging improved EMF mainly by increasing functional diversity. Conclusion: Our study demonstrates that selective logging is a good management technique from an EMF perspective, and thus provide us with potential guidelines to improve forest management in P. yunnanensis forests in this region. The functional diversity is maximized through reasonable selective logging measures, so as to enhance EMF.

indicate that selective logging could change the carbon production and allocation (Riutta et al. 2018), impact soil organic carbon and nutrient stocks (Lontsi et al. 2019), and in uence decomposition (Ojanen et al. 2017;Wang et al. 2018). The simultaneous occurrence of these different ecosystem functions (EFs) is de ned as ecosystem multifunctionality ('EMF') (Hector and Bagchi 2007). This concept was rst proposed in the study of seagrass (Duffy et al. 2003) and has been widely used in grasslands (Soliveres et al. 2016;Meyer et al. 2018), aquatic (Lefcheck et al. 2015;Perkins et al. 2015) and forest ecosystems (van der Plas et al. 2016;Huang et al. 2019;Luo et al. 2019).
We know little about the effect of selective logging on multiple EFs and EMF. Due to the functional uniqueness of the species, biodiversity was crucial for EMF. Studies on the relationship between biodiversity and EMF ('BEMF') have become a hotspot in ecology (Zavaleta et al. 2010;Gamfeldt and Roger 2017). Biodiversity consists of three components: taxonomic diversity (de ned as species abundance and richness), functional diversity (de ned as various resource utilization strategies and growth types), and phylogenetic diversity (de ned as the presence of different evolutionary lineages) (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, since it is a simple metric to study 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). Plant functional traits are the characteristics of morphological, phenological and physiological, which could re ect how plants respond to environmental changes, and in uence other trophic levels and different ecosystem processes (Pérez-Harguindeguy et al. 2013). Functional diversity is based on functional traits, which are closely related to resource utilization (Díaz and Cabido 2001). Therefore, functional diversity might be a better predictor of ecosystem function or multifunctionality than taxonomic diversity (Díaz and Cabido 2001;Gross et al. 2017). 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 or multifunctionality, comparing to taxonomic or functional diversity (Flynn et al. 2011;Srivastava et al. 2012;Bagousse-Pinguet et al. 2019). The effects of different biodiversity attributes on EMF might vary in magnitude (Bagousse-Pinguet et al. 2019). In addition to 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. In fact, biodiversity is not the only aspect of biological systems that affects EMF, the effects of various abiotic factors on EMF have also received much attention (Jing et al. 2015;Giling et al. 2019).
Pinus yunnanensis Franch. is the main pioneer species in forests of southwest China and is widely used in construction (Xu et al. 2016). The P. yunnanensis forests range 23° to 30° N and 96° to 108° E. They grow mainly in river valleys between 700 and 3000m above sea level (Li et al. 2020). The P. yunnanensis forests can provide a large amount of high-quality timber and industrial raw materials, which play an important role in maintaining species diversity, conserving water resources and retaining soil (Xu et al.

Page 4/29
2016; Li et al., 2020). However, due to excessive human disturbance activities and the invasion of pests, diseases and alien species, the EFs of P. yunnanensis forests are not fully developed (Xu et al. 2016;Liu et al. 2019). Therefore, how to improve the EFs of P. yunnanensis forests is an important task for forest management.
The effect of multiple plant diversity attributes, soil microbial diversity and abiotic factors on EMF under different selective logging intensities is unclear. To ll these knowledge gaps, we quanti ed four key ecosystem functions: nutrient cycling (including soil total nitrogen, soil hydrolysable nitrogen, soil total phosphorus and soil available phosphorus), soil carbon stocks (soil total carbon), decomposition (including β-glucosidase, urease and phosphatase activity) and wood production (aboveground woody biomass). In this study 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) in uenced EMF under ve different logging intensities applied to a P. yunnanensis forest. We predicted that: (1) in the absence of management (i.e. no logging) overall lower values of diversity and EMF. Thus, selective logging could increase EMF both directly and indirectly via taxonomic, functional and/or phylogenetic diversity; (2) effects of functional and phylogenetic diversity would be greater; (3) soil microbial diversity could improve EMF. The results of this study could provide new insight for forest management and conservation.

Study area
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 classi ed as hilly red soil. The main forest type in this area is P. yunnanensis forest. P. yunnanensis forest is mainly divided into two layers, the upper layer is P. yunnanensis, which occupies the absolute dominant position, the main species of the lower layer is yunnanensis forests is the selective logging for valuable timber. The Yongren National Forest Farm has been producing timber since the 1950s (Li et al. 2020).

Sampling and experimental design
We obtained eld data from 52 plots in April 2018. These plots (20m × 30m) were representative of the major vegetation types and management styles in the area. The plots were set as Li et al. (2020) described: selectively logged once (SL1, 10 plots, Fig. 2a), logged twice (SL2, 10 plots, Fig. 2b), logged three times (SL3, 10 plots, Fig. 2c), logged four times (SL4, 10 plots, Fig. 2d), logged ve times (SL5, 9 plots, Fig. 2e), control group without selective logging (CG, 3 plots, Fig. 2f). The forest stand characteristics (stand density, basal area and height) of these treatments are presented in Additional le: Table S1. The selective logging began in the 1970s, the interval of selective logging is about 7-8 years.
Selected species included P. yunnanensis and broadleaf species. On average, 10% of the number of individuals was selectively harvested at one time.
Follow the standard eld protocol, within the plots, we identi ed 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 (Condit 1998;Yuan et al. 2020). A total of 13 species were investigated, belonging to 7 families and 11 genera (Table S2).
To obtain plant functional traits to calculate functional diversity, following the standardized protocols of Cornelissen et al. (2003), we gathered 50 to 100 sun-exposed mature leaves from ve to ten individuals of each species per plot. We labeled the leaves by species and put them into 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 ne powder to prepare for the measurement of leaf N concentration and leaf P concentration.
In each plot, we collected the wooden core for woody plants (DBH > 5cm) by increment borer. Then, the volume was measured using the diameter and the length of the wooden core. For other woody plants (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.
To measure soil enzyme activity, belowground diversity and soil physiochemical properties, 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 rst 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. All measurements of physiochemical data are presented in Supplementary Methods 1.

Assessing aboveground and belowground diversity
Species richness was calculated for taxonomic diversity. As a classical measure of functional diversity, Rao's quadratic entropy (RaoQ) is a multivariate measure of functional divergence which could increase ecosystem function due to e cient use of resources (Mason et al. 2005) and has been applied to the study of the relationship between biodiversity and ecosystem function (Botta-Dukát 2005; Wasof et al. 2018). RaoQ was calculated using speci c 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 phylogenetic diversity (PD) was calculated as a proxy of phylogenetic diversity (Faith 1992). The number of soil bacterial and fungal operational taxonomic units (OTUs) was used to assess belowground diversity. Details on the several molecular techniques employed for assessing belowground diversity are presented in Supplementary Methods 2.

Individual EFs
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 (Lucas-Borja and Delgado-Baquerizo 2019). The content of soil total carbon was used as the proxy of soil carbon stocks (Gamfeldt et al. 2013). The three enzyme activities (βglucosidase, urease, phosphatase) were used as the proxy of decomposition (Lucas-Borja and Delgado-Baquerizo 2019). Finally, the aboveground woody plant biomass was used as the proxy of wood production (Mensah et al. 2020). We used the allometric equation form Li et al. (2020) to estimate the aboveground biomass of individual trees of P. yunnanensis forest. A list of the allometric regressions used in our study is included in Table S3. 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. 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 rst built the list of species found in our eld 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 in R.

Effect of selective logging intensity on EFs
Selective logging had different effects on different EFs (Table 2). Nutrient cycling and soil carbon stocks responded similarly to increasing selective logging intensity. It showed an increasing trend from SL1 to SL3, and a decreasing rst and then increasing trend from SL3 to SL5. The values of the SL3 of nutrient cycling and soil carbon stocks were signi cantly greater than other selective logging intensities and CG (P < 0.05), except the SL5. Decomposition responded directly to increasing selective logging intensity, with the values of the SL4 and SL5 being signi cantly greater than the CG (P < 0.05), and the SL5 showing the highest value. For wood production, the value of the SL4 was signi cantly greater than other selective logging intensities and CG (P < 0.05), except the SL2. The value of the SL1 was signi cantly lower than SL2, SL3 and SL4 (P < 0.05), except the CG and SL5. We found a signi cant positive correlation between selective logging intensity and nutrient cycling (p < 0.01) (Table S4), soil carbon stocks (p < 0.01) (Table S4), decomposition (P < 0.01) (Table S4) and wood production (p < 0.01) (Table S4). With the increase of selective logging, nutrient cycling (P < 0.01) ( Fig. 3a), soil carbon stocks (P < 0.01) (Fig. 3b), decomposition (P < 0.001) (Fig. 3c) and wood production (P < 0.05) (Fig. 3d) all showed an increasing trend. Selective logging intensity explained 17.1% of variation of nutrient cycling, 14.7% of variation of soil carbon stocks, 33.2% of variation of decomposition and 11.4% of variation of wood production (Fig. 3).

EMF under different selective logging intensities
We found that that except for SL1, the rest of treatments signi cantly improved EMF when compared to the control treatment (CG), showing a signi cant difference between SL1 and other logging intensities (p < 0.05) (Fig. 4). Interestingly, no signi cant differences among SL2, SL3, SL4 and SL5 were found with regards to EMF (p > 0.05) (Fig. 4).
Relationships between selective logging intensity and above-, belowground biodiversity and abiotic factors There was a signi cant positive correlation between selective logging intensity and functional diversity (p < 0.05, τ = 0.217) (

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

Effect of selective logging, biotic and abiotic factors as predictors of EMF
Selective logging might directly and indirectly through functional diversity impact EMF (P < 0.05) (Fig.   7a). 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. 7a). Soil fungal OTUs could also directly impact EMF (P < 0.05), but there was no signi cant relationship between selective logging and soil fungal OTUs (P > 0.05) (Fig. 7a). A strong relationship in the SEM analysis was found between logging and soil moisture, but we also found a weaker positive relationship between soil moisture and EMF (P > 0.05) (Fig. 7a). The standardized total effects of functional diversity, soil moisture, soil fungal OTUs and selective logging intensity derived from SEM were 0.292, 0.364, -0.192 and 0.558, respectively (Fig. 7b).

Discussion
Our study explored the relationship between above-and belowground biodiversity and abiotic factors and EMF under different logging intensities. The results basically met our rst prediction: selective logging could directly signi cantly improve EMF except logging one time also indirectly signi cantly improve EMF through functional diversity. Logging intensity in uenced positively nutrient cycling, soil carbon stocks, decomposition and wood production. Different EFs showed different response patterns to logging intensity. Among biodiversity elements, logging intensity only in uenced positively functional diversity.
Similarly, logging intensity only in uenced positively soil moisture between the two abiotic factors. Contrary to our second prediction, the effect of functional diversity, not phylogenetic diversity, on EMF was more signi cant compared to taxonomic diversity. Neither phylogenetic diversity nor taxonomic diversity had a signi cant effect on EMF. In addition, the relationship between belowground biodiversity and EMF was also contrary to our third prediction: only soil fungal OTUs, not soil bacterial OTUs negatively in uenced EMF.
Selective logging is a widespread forest management technique that could alter forest microclimate. Studies show that air temperature, soil temperature, relative air humidity and soil moisture increased as canopy density and tree density decreased after selective logging (Masyagina et al. 2006;Oldén et al. 2019). Variations of temperature and moisture in forest microclimate directly affect some key processes in forest ecosystems (Vanwalleghem and Meentemeyer 2009). Selective logging 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 logging intensity promoted multiple EFs. Consistent with previous studies (Hwang and Son 2006;Kim et al. 2016), we observed that soil carbon stocks increased with logging intensity. Changes in soil carbon stock were caused by root mortality and the amount of dead woody vegetation entering the soil at different logging intensities (Kim et al. 2016). After logging, changes in vegetation structure and species composition, together with the addition of necromass, increased the decomposition rate of organic carbon and thus increased the soil carbon stocks (Barros and Fearnside 2016). Wang et al. (2018) indicated that the amount of upper canopy removal would affect decomposition rates. Selective logging reduced canopy density and allowed more sunlight and precipitation to reach the ground (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 logging is one of the main explanations for the increase in decomposition. In our study, logging intensity led to an increase in nutrient cycling, a relationship caused by the nutrient pool formed by logging residues entering the soil; the loss of individuals, reducing competition for nutrients; the increase of leaching loss (Carlyle 1995). Nutrient cycling implies a positive feedback between logging intensity and decomposition (Kim et al. 2016). In our study, wood production increased with logging intensities because logging 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 logging on these four EFs leads to the direct positive effect of selective logging on EMF.
There was no signi cant correlation between logging intensity and taxonomic diversity. This result was consistent with the results of many previous studies (Ares et al. 2010;Nunes et al. 2014). Selective logging might change the structure of the community but not the number of species. In the eld, we did not nd any invasive species. Phylogenetic diversity, supposed to be the most expected predictor of EMF, was not signi cantly related to EMF, perhaps because the phylogenetic signal of functional traits related to each ecosystem function might not always be present (Davies et al. 2016). Zirbel et al. (2019) suggested that if the functional traits that determined an ecosystem function were not phylogenetically conserved, phylogenetic diversity might not be a good predictor of that ecosystem function. We concluded that some species, functional traits, and phylogenetic signals were lost during logging. The disappearance of a species means the loss of evolutionary history (Ribeiro et al. 2016). In addition, phylogenetic diversity is a measure of biodiversity that combines species richness and evolutionary history (Cadotte 2013). The fact that taxonomic diversity was not related to EMF in our study might also lead to the failure of phylogenetic diversity prediction. Compared to taxonomic and phylogenetic diversity, functional diversity partly in uenced EMF under selective logging intensity. Moreno-Gutiérrez (2011) reported that the changes of soil conditions and light after logging intensity potentially in uence plant functional traits. Our result that logging improves functional diversity aligns with previous studies (Ares et al. 2010;Nunes et al. 2014). It has been proven that functional diversity was the main driver of EMF in a P. yunnanensis natural secondary forest (Huang et al. 2019). Increased functional diversity would aid forest ecosystems in responding to global climate change (Ares et al. 2010). In the forest ecosystem with high functional diversity, greater ecological differentiation enables coexisting species to survive through niche partitioning and complementarity of resource use among species, thus improve community resilience to maintain more EFs (Valencia et al. 2015;Mensah et al. 2020). Contrary to other ndings (Wagg et al. 2014;Jing et al. 2015;Delgado-Baquerizo et al. 2016;Delgado-Baquerizo et al. 2017), this study found that soil biodiversity (mainly soil fungal OTUs) had a signi cant negative correlation with EMF and showed similar negative trends with multiple EFs. This nding was interesting, though the speci c causality is unknown. One potential answer is that the average of multiple functions included in the EMF index might be a poor metric to accurately re ect the multiple functions performed by soil fungi and the complex interactions between soil and plant fungi (Jing et al. 2015).
Our ndings suggest that selectively logged twice or more is necessary if EMF is to be enhanced for forest managers. During selective logging, species with large differences in traits should be retained as much as possible. In the current study, although the expected results have been achieved, four EFs we used can't represent all functions performed by ecosystems. The number of measured functions may in uence the results of the study (Meyer et al. 2018). Gamfeldt and Roger (2017) pointed out that identifying mechanisms was the key to the success of BEMF eld. Therefore, future studies need to identify underlying mechanisms behind the effects of biotic and abiotic factors on EMF.

Conclusions
Our results suggest that selective logging of P. yunnanensis forests was a signi cant driver of EMF when it is applied at least twice or more. However, this does not mean that more selective logging is better. For

Availability of data and materials
The datasets analysed during the current study are available from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests.  The location of the study at Yongren National Forest Farm (black squares, n=52). CG: control group; SL1: selective logging one time; SL2: selective logging two times; SL3: selective logging three times; SL4: selective logging four times; SL5: selective logging ve times.

Figure 2
Sampling was implemented under different selective logging intensities in the study location. SL1: selective logging one time (a); SL2: selective logging two times (b); SL3: selective logging three times (c); SL4: selective logging four times (d); SL5: selective logging ve times (e); CG: control group (f). Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

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