Aboveground carbon storage is driven by functional trait composition and stand structural attributes rather than biodiversity in temperate mixed forests recovering from disturbances

Abstract

Key message

Functional trait composition and stand structural complexity rather than biodiversity substantially enhance aboveground carbon storage in temperate mixed forests, while accounting for the effects of disturbance intensity. This study provides a strong support to the mass ratio effect in addition to the niche differentiation and facilitation effects.

Context

The underlying mechanisms for the relationships between biodiversity and ecosystem function remain hotly debated for the last four decades.

Aims

We tested how do biodiversity, functional trait composition, stand structural attributes, and topographic variables explain aboveground C storage under different disturbance regimes.

Methods

We used linear mixed effects and structural equation models to simultaneously evaluate the effects of biodiversity, stand structure attributes, functional trait composition, and topographic variables on aboveground C storage while considering for the effects of disturbance intensity. We used biophysical data from 260 plots within 11 permanent temperate mixed forests in Northeastern China.

Results

Aboveground C storage was driven by stand basal area, individual tree size inequality, community-weighted mean of maximum height and wood density, and diversity (functional evenness and mean nearest taxon distance). The structural equation model showed that aboveground C storage was positively affected by individual tree size inequality and trait composition (i.e., CWM of maximum height), after accounting for the strongest negative direct and indirect effects of disturbance intensity.

Conclusion

Conserving functional identity of species and maintaining complex stand structure would be the alternative choices for higher aboveground C storage in temperate mixed forests.

1 Introduction

Forests harbor about two thirds of the terrestrial biodiversity and account for 44% of the world’s forest carbon pool (Pan et al. 2011). Therefore, understanding the relationships between biodiversity and ecosystem functions such as aboveground biomass or carbon (C) storage and productivity in natural forests is crucial for predicting the consequences of biodiversity loss on ecosystem functioning and services (Cardinale et al. 2012; Forrester and Bauhus 2016). Most of the previous studies have reported positive relationships between species diversity and ecosystem functions across forest types and biomes (Scherer 2014; Liang et al. 2016; Ali and Yan 2017b). Several studies have suggested that ecosystem functioning is related to either taxonomic, functional, or phylogenetic diversity (Cadotte et al. 2008; Cavanaugh et al. 2014). However, little is known about the relative importance of different metrics of biodiversity, stand structural attributes, and environmental factors on aboveground C storage while accounting for the effects of disturbance intensities in natural forests.

Table 1 Environmental and forest characteristics of the 11 studied sites included in this study

Table 2 The best multiple model obtained from a series of linear mixed-effect models for predicting aboveground C storage based on 260 plots

Table 3 The direct, indirect, and total standardized effects of disturbance intensity, stand structure attributes, functional trait composition, and diversity on the aboveground C storage based on structural equation model (SEM) in temperate mixed forests

Fig. 1

Distribution of the study sites in and around the Changbai mountain region in China. Site names are given in Table 1

Fig. 2

Boxplots of aboveground C storage and multiple metrics of biodiversity and stand structural attributes for different disturbance intensity levels. Different letters (a–c) indicate significant differences at p < 0.05 (Tukey’s test). Factors are aboveground C storage (AGC), species richness (S), Pielou’s evenness index (E), Faith’s phylogenetic diversity (PD_F), mean pairwise distance phylogenetic distance between all pairs of species on the phylogenetic tree (MPD), mean of the shortest distances of species on the phylogenetic tree (MNTD), functional richness (FRic), functional evenness (FEve), functional dispersion (FDis), coefficient of variation in DBH (CV_DBH), tree DBH diversity based on 6 cm width (H_d), stand basal area (BA), community-weighted mean of wood density (CWM_WD), and maximum height (CWM_MH)

Fig. 3

The best-fit structural equation model for linking disturbance levels, diversity (FEve), stand structure attribute (CV_DBH), functional trait composition (CWM_MH), and aboveground C storage in the temperate mixed forests. Solid arrows represent significant paths, and dashed arrows represent nonsignificant paths. For each path, the standardized regression coefficient is shown. R² indicates the total variation in a dependent variable that is explained by the combined independent variables. The detail model-fit statistics summary is provided in Table S6 (χ² = 1.31, P value = 0.25, RMSEA = 0.035, SRMR = 0.016, CFI = 1.00). *P < 0.05; **P < 0.01

Importantly, the world’s forest area is declining due to anthropogenic disturbances (e.g. logging), thus raising important concerns for biodiversity loss and ecosystem functioning (Asner et al. 2009; Thom and Seidl 2015). For example, disturbance intensities have a stronger effect on the relationships between species diversity and aboveground biomass across different biomes (Grace et al. 2016; Yeboah and Chen 2016; Sanaei et al. 2018). The intermediate disturbance hypothesis predicts that moderate intensities of disturbance maintain the highest diversity, resulting in a hump-shaped relationship between species diversity and ecosystem function (Connell 1978). Although multiple drivers affect aboveground C storage (Paquette and Messier 2011; Lasky et al. 2014; Ali et al. 2016), few studies have tested whether the importance of these drivers might change across disturbance intensities (Kröber et al. 2015; Jucker et al. 2016).

Two nonmutually exclusive hypotheses have been proposed to explain the relationships between biodiversity and ecosystem functions: (1) the niche complementarity hypothesis and (2) the mass ratio hypothesis. The first hypothesis suggests that high species or functional trait diversity can enhance ecosystem functions, and hence biomass productivity increases through resource-use efficiency by component species or interacting individuals due to the niche partitioning and/or facilitation (Tilman et al. 1997; Zhang et al. 2012; Forrester and Bauhus 2016). Under this expectation, recent studies have reported that functional trait and/or phylogenetic diversity indices are better predictors for aboveground biomass or productivity than species richness, because they better capture the degree of functional redundancy and niche overlap (Paquette and Messier 2011; Lasky et al. 2014; Ruiz-Benito et al. 2014). The second hypothesis suggests that ecosystem functions, at a given time, is mainly determined by the trait(s) of most dominant species (Grime 1998). Under this expectation, higher aboveground C storage could be closely related to the functional traits of the dominant species, which can be quantified using the functional trait composition or identity such as community-weighted mean (CWM) of a certain functional trait (Cavanaugh et al. 2014; Tobner et al. 2016; Ali et al. 2017). Understanding how the relative contributions of these two mechanisms change according to the context may inform about new strategies for the forest management and biodiversity conservation (Reich et al. 2001; Cavanaugh et al. 2014; Lin et al. 2016; Ratcliffe et al. 2016). For example, the higher relative importance of the mass ratio effect on tree growth at latitudinal extremes of the European continent suggests the key role of certain traits in harsh-climate forests, whereas the niche complementarity in Mediterranean water-limited forests strongly promotes ecosystem functioning (Ratcliffe et al. 2016).

Stand structural attributes are often characterized using stand-level metrics of tree size inequality or diversity based on tree heights or diameters at breast height (DBH) (Pretzsch 2014; Ali et al. 2016). In natural forests, stand structure is strongly shaped by inter- and intra-specific competition for resources (Clark 2010; Zhang and Chen 2015). Consequently, it can also be considered a measure of the realized degree of species complementarity through niche differentiation and facilitation (Yachi and Loreau 2007; Clark 2010; Zhang and Chen 2015). Recent studies have suggested that individual tree size inequality or stand structural diversity (i.e., tree DBH and height diversity) rather than species diversity substantially enhances aboveground C storage, biomass, and productivity in natural forests (Ali et al. 2016; Danescu et al. 2016). However, other studies found that increasing the heterogeneity of stand structure reduces aboveground biomass or productivity (Binkley et al. 2010; Ryan et al. 2010; Soares et al. 2017), but the direction of these relationships may vary across natural forests, plantations, and manipulated experimental forests.

In this study, we tested how aboveground C storage was driven by biodiversity (i.e., taxonomic, functional, and phylogenetic), functional trait composition, and stand structural attributes across several disturbance intensities (i.e., timber harvesting or logging) in temperate mixed forests while accounting for the effects of topographic factors as covariates. Specifically, we addressed the following objectives by analyzing biophysical data from 260 subplots within 11 permanent temperate forests in Northeastern China: (1) to identify differences in diversity, structural attributes, and carbon storage depending on disturbance; (2) to select the best combination of multiple metrics of biodiversity, stand structural attributes, and topographic factors explaining aboveground C storage; and (3) to test the direct, indirect, and total effects of disturbance intensities on aboveground C storage via biodiversity, stand structure, and functional composition.

2 Materials and methods

2.1 Study area, sites, and forest plots

The datasets used in this study originate from the “Project of coniferous-deciduous mixed temperate forest diversity survey in northeast China,” which was launched to understand the effects of environmental changes and disturbance intensities on biodiversity and ecosystem functions and services (Song et al. 2014). This study was based on the 11 permanent sites in temperate mixed forests (site size varied between 0.6 and 1 ha) on the Changbai Mountain (40° 54′ to 44° 03′ N, 124° 47′ to 130° 09′ E), located in Liaoning and Jilin Provinces in Northeastern China (Table 1; Fig. 1). The region is characterized by a temperate continental climate with long, cold winters and warm summers. The mean annual temperature is 2.8 °C, and the temperature of the warmest and coldest months are 19.6 and − 13 °C, respectively (Zhang et al. 2005; Hao et al. 2007). The annual precipitation is 700 mm (mostly fall between June and September), the annual evaporation is 1250.9 mm, and the relative humidity is 72%. The altitude of the study area varies between 650 and 1020 m a.s.l (Zhang et al. 2005). The zonal vegetation is characterized as a broad-leaved Korean pine (Pinus koraiensis) mixed forest (Hao et al. 2007), and the soils are classified as the dark-brown soil according to the FAO soil classification system (Yang and Li 1985).

Besides the core area of Changbai Mountain, the natural broad-leaved Korean pine forest is distributed in some residual fragments with unequal areas in Changbai Mountains and Small Hinggan Mountains. In order to monitor this forest type on different sites, we first approximated its distribution area using satellite imagery, and then established permanent plots in suitable positions from July 2012 to July 2013 (Song et al. 2014). The basic sampling information, geographic location and forest structure of the 11 study sites are described in Table 1. Within each plot, all individual trees with a DBH ≥ 1 cm in contiguous 20 × 20 subplots were tagged, identified, and measured, and their geographic coordinates were recorded following a standard field protocol (Hao et al. 2007). A total of 22,766 stems were recorded belonging to 81 species, 46 genera, and 26 families.

2.2 Assessments of disturbance intensities

The study area has been subjected to both natural disturbances and variable intensities of human disturbances. All studied forest sites have been protected from anthropogenic disturbance since the implementation of Natural Forest Protection Project in 1998 (Dai et al. 2004). Forests in the study area contain stands with different successional stages (Chen et al. 2014).

The disturbance intensity of each plot was assessed by counting the number of tree stumps that had been removed in the field using a chainsaw (Kahl and Bauhus 2014). In addition, the official records of the Local Forestry Bureau, Jilin and Liaoning Provinces, were reviewed to extract selective logging data, by following the method of Ali et al. (2017). After that, plots were categorized into three disturbance intensity levels based on the percent of removed stems: low (< 10%), medium (10–30%), and high (> 30%; Table 1). Plots with a low level of disturbance were located in intact and relatively undisturbed forests within the core zone of the Changbai Mountain Nature Reserve (Fig. 1). Established in 1960, this reserve is part of the World Biosphere Reserve Network under the Man and the Biosphere Project in 1980 (Yang and Li 1985).

2.3 Estimation of aboveground carbon storage

The aboveground biomass of each individual tree was estimated from DBH using species-specific available allometric equations (covering ~ 80% of the studied tree species) and a generic allometric equation for the remaining species (Table S1 in the Electronic supplementary material). We summed individual aboveground biomass of all trees to obtain total aboveground biomass within each subplot (converted to Mg ha⁻¹). Finally, we multiplied aboveground biomass by a factor of 0.5 to derive the aboveground C storage for each subplot (Guo et al. 2010).

To test the potential effects of other variables such as stand density and climate on biomass estimations, we applied an alternative biomass expansion factor method through converting timber volume derived from DBH to forest biomass (Fang et al. 1998, 2001). The results revealed a strong correlation with species-specific allometric regression equations (Pearson’s r = 0.92). Therefore, we reported only the results based on allometric equations.

2.4 Quantification of multiple metrics of biodiversity

We measured three components of biodiversity within each plot: taxonomic, phylogenetic, and functional trait diversity. Taxonomic diversity indices were computed based on tree species richness and evenness (Pielou’s), using the vegan package (Oksanen et al. 2013) in R 3.3.3 (Team RDC 2017).

Functional trait diversity indices were calculated based on two functional traits, i.e., tree maximum height and wood density, which are closely related to the life-history strategies of the species. Tree maximum height provides a proxy for the potential tree height, which is considered an important indicator of light capture strategy (Kunstler et al. 2016). Wood density reflects the strategy of a species to allocate resources into fast growth and early reproduction or slow growth and resistance to environmental hazards (Kunstler et al. 2016). Typically, fast-growing species have a higher light demand and lower wood density than slow-growing conservative species (Spasojevic and Suding 2012). For a given species, the maximum height was estimated by the largest DBH value within the dataset. Specifically, we first selected the top 10 individuals with largest DBH from the dataset and measured their tree height using a laser rangefinder (Laser Technology Inc.). Based on the recommendations of Pérez-Harguindeguy et al. (2013) and Paquette and Messier (2011), we measured wood density values for 80% of the dominant species across the plots. For the rest of species (n = 9), wood densities were collected from a global wood density database (Chave et al. 2009). For the missing species (n = 7) in the global database, we used the mean family wood density value of the specific species.

Three multidimensional functional trait diversity indices are commonly used in biodiversity and ecosystem functioning studies to explore different facets of functional diversity in multidimensional trait space (Fotis et al. 2017; van der Sande et al. 2017). Functional richness (FRic) is the amount of multivariate trait space filled by the community. Functional evenness (FEve) indicates how species’ basal area spread over the multivariate-trait space and the index is higher when basal area distribution is homogeneous across trait space. Functional dispersion (FDis) is the average distance of the species to the basal-area-weighted centroid of all species in community trait space (Laliberté and Legendre 2010).

Since phylogenetic diversity (PD) reflects the evolutionary history of species, it is increasingly recognized as a predictor of the relationship between biodiversity and ecosystem function (Cadotte et al. 2008; Kembel et al. 2010). The relationship lies on the assumption that evolutionary diversification has generated trait diversification, which in turn could result in greater niche complementarity (Cadotte et al. 2008). Phylogenetic diversity indices were computed using the informatics tool Phylomatic (http://www.phylodiversity.net). Phylomatic utilizes the Angiosperm Phylogeny Group III (APG III 2010) phylogeny as a backbone. This supertree was then assigned branch length estimated from multi-gene molecular and fossil data implemented in Phylomatic, which provides the largest and most up-to-date time calibrated species-level phylogeny of seed plants (Zanne et al. 2014). From these phylogenies, we then calculated several measures of phylogenetic diversity including Faith’s phylogenetic diversity (PD_F, the total length of all branch lengths separating taxa in a community; the mean pairwise phylogenetic distance (MPD) and mean nearest taxon distance (MNTD) within a community. PD_F take into the consideration of the phylogenetic relations between the species in the community but not the relative abundance of species. Mean pairwise phylogenetic distance (MPD) and mean nearest taxon distance (MNTD) represent a different measurement of evolution information among species, which are commonly used in many recent studies (Erikson et al. 2014; Tucker et al. 2017). PD indices were calculated with the picante package in R 3.3.3 (Kembel et al. 2010).

Functional trait composition was quantified as the CWM of functional trait values. The CWM of a single trait was calculated as the mean trait value in each subplot (CWM_MH and CWM_WD), weighted by the species’ relative basal area. This metric represents the expected functional trait value of a specific community (i.e., subplot) (Díaz et al. 2007). Trait values were standardized before the calculation of functional trait diversity and CWM indices. All indices were calculated using the FD package in R 3.3.3 (Laliberté and Legendre 2010).

2.5 Quantification of stand structural attributes

For the quantification of stand structural attributes, we used stand basal area (BA), Shannon’s diversity index-based tree DBH diversity, and the coefficient of variation of DBH within each subplot. We calculated stand basal area through the summation of basal area of all individual trees within each subplot. It is worth noting that aboveground biomass or C storage of each individual trees is calculated based on tree DBH, and hence individual tree biomass scales closely with the basal area of the individual tree rather than stand basal area (Poorter et al. 2015). As recommended by Ali et al. (2016), tree DBH diversity (H_d) (Eq. 1) was calculated by considering several class widths for DBH (i.e., 2, 4, 6, and 8 cm). The species’ relative basal area was used to weight the proportion within each DBH class because the basal area is more related to biomass than relative abundance (Zhang and Chen 2015). Tree DBH diversity based on different discrete classes may influence the variation in aboveground C storage, and therefore, we evaluated the bivariate relationship between aboveground C storage and each of tree DBH diversity for each class and selected the class that conducted to the model with the lowest AIC (Table S2 in the Electronic supplementary material).

$$ {H}_d=-{\sum}_{i=1}^d{p}_i\times \ln \left({p}_i\right) $$

(1)

where p_i is the proportion of the basal area of ith DBH classes while d is the number of DBH classes within each subplot.

The coefficient of variation (CV) was used to quantify DBH variation within each subplot as a proxy of individual tree size inequality (Zhang and Chen 2015), because the overall variation in DBH reflect the extent of the realized niche differentiation via positive plant-plant interactions (Yachi and Loreau 2007; Ali et al. 2016).

2.6 Assessment of topographic variables

We initially assessed three topographic variables for each subplot including mean elevation, slope, and convexity. Mean elevation (Elev) was obtained as the mean value of the elevations of the four corners of a subplot, using an electronic distance measuring device. Following Harms et al. (2001), convexity was computed as the difference between the elevation of the focal subplot and the mean elevation of the eight surrounding subplots. For the edge subplot, convexity was taken as the elevation of the center points minus the mean elevation of the four corners. Each subplot was divided into four triangular planes, each formed by joining three corners of the subplot. The slope was the mean angular deviation from horizontal of each of the four triangular planes formed by connecting three of its corners (Harms et al. 2001).

2.7 Statistical analyses

We first applied Tukey’s test with one-way analysis of variance (ANOVA) in order to check for significant differences among three levels of disturbance intensities. Prior to the following analyses, aboveground C storage was natural-log transformed and all explanatory variables were standardized to obtain a mean of 0 and a standard deviation of 1, which aimed to improve the interpretability of regression coefficients (Schielzeth 2010). All statistical analyses were conducted in R 3.3.3 (Team RDC 2017).

Then, we tested for the spatial autocorrelation in aboveground C storage among subplots in each site by fitting generalized least-square (GLS) models without and with spherical autocorrelation structure for the relationship between each predictor and aboveground C storage and compared the GLS models (spatial and nonspatial) using the Akaike Information Criterion (AIC). Our results showed that models without spatial autocorrelation structures always had the lowest AIC values (Table S3 in the Electronic supplementary material), as similarly reported by previous studies in natural forests (Ali and Yan 2017b).

We did not find any strong confounding effect of spatial autocorrelation. Then, we tested the bivariate relationships between aboveground C storage and the predictors using ordinary least square regressions analysis in pooled data. We aim to preselected noncollinear variables for the mixed-effects models. After, we fitted linear mixed-effect models in order to simultaneously assess the effects of multiple predictors on aboveground C storage (second objective). Here, we included only those predictors which had significant bivariate relationships with aboveground C storage. In addition, we removed highly correlated predictors of aboveground C storage (i.e., r > 0.6), such as species richness and MPD (see the correlations among candidate predictors in Table S4 in the Electronic supplementary material). The full multiple linear mixed model included four biodiversity indices (such as PD_F, MNTD, FEve, and FDis), two functional trait composition indices (CWM_MH and CWM_WD), three stand structure attributes (BA, H_d, and CV_DBH), two topographic variables (elevation and slope), and disturbance intensity as fixed factors, while site was included as a random factor. The disturbance intensity was an ordinal categorical variable and was coded as 1 (low), 2 (medium), and 3 (high), being treated as a regular numeric variable as recommended (Zhang and Chen 2015). The best model was selected by considering the lowest AICc and number of predictors (Table S5 in the Electronic supplementary material) (Burnham and Anderson 2003; Bartoń 2016). The conditional (c) and marginal (m) R² for the multiple linear mixed models were calculated (Nakagawa and Schielzeth 2013). Conditional R² (R²_c) indicates the variance explained by both fixed and random factors, whereas marginal R² (R²_m) indicates the variance explained by fixed factors only (Prado-Junior et al. 2016). Multiple linear mixed models were conducted using the “nlme” package, and model selection was performed using the dredge function of “MuMIn” package (Bartoń 2016).

Finally, we applied the structural equation model (SEM) to test our third objective of direct and indirect effects of disturbance affecting through the modification of stand structure, diversity, and functional composition. Here, we used partially confirmatory SEM model due to the structure of conceptual model being fixed as we knew disturbance, abiotic, and biotic factors can influence aboveground C storage (Ali and Mattsson 2017a; van der Sande et al. 2017), but we did not know the relative importance of different combinations of variables within each category (such as trait composition) to aboveground C storage. According to the most parsimonious multiple linear mixed model (Table 2), we selected two variables related to diversity indices (MNTD and FEve), stand structure attributes (BA and CV_DBH), and trait composition (CWM_MH and CWM_WD) as the complete model. The best-fit SEM was evaluated based on chi-square (χ²) test (i.e., P value > 0.05) with evaluation of additional important indices or values such as highest Bentler’s Comparative Fit Index and Goodness of Fit Index (i.e., CFI and GFI > 0.9), coefficient of determination (R²), and standardized root mean square residual (SRMR ≤ 0.05) and lowest AIC value (Table S6 in the Electronic supplementary material). In the main result, we only reported the results derived from the selected best-fitted SEM. The SEM analysis was implemented using the “lavaan” package (Rosseel 2012).

Data availability

The datasets generated and/or analyzed during the current study are available from the corresponding author on reasonable request.

3 Results

3.1 Aboveground C storage and biodiversity indices across disturbance intensities

The mean aboveground C storage in highly disturbed plots was 55.6 Mg C ha⁻¹, which was significantly lower than plots with medium (88.6 Mg C ha⁻¹) and low disturbance levels (128.2 Mg C ha⁻¹; Fig. 2). In contrast, mean species richness increased from 30 to 42 with increasing disturbance intensity (Table 1). Tukey’s test indicated that plots with a low disturbance level had a low species richness and Faith’s phylogenetic diversity than plots with medium and low disturbance levels (Table 1; Fig. 2), but they had higher functional trait diversity indices (FEve and FDis), tree size inequality (CV_DBH), and CWM of tree maximum height (Fig. 2). Other two phylogenetic diversity indices (MPD and MNTD) were lowest in highly disturbed plots compared with plots with two other disturbance levels (Fig. 2).

3.2 Bivariate relationships of aboveground C storage with each of multiple predictors

Bivariate relationships, using pooled data, indicated that aboveground C storage significantly increased with phylogenetic diversity indices (i.e., MPD and MNTD) and functional trait diversity indices (i.e., FEve and FDis), and it decreased with species richness and Faith’s phylogenetic diversity (Fig. S1). The CWM of tree maximum height had the strongest positive relationship with aboveground C storage (R² = 0.35, P < 0.001) but a nonsignificant relationship with CWM of wood density. Aboveground C storage was positively related to stand structure (CV_DBH and BA). There was no significant association of aboveground C storage with species evenness (E) and functional richness (FRic). Among the three topographic variables, mean elevation of subplot had a positive relationship with aboveground C storage, but slope had a negative relationship, whereas convexity was not significantly related (Fig. S1).

3.3 Global multiple linear mixed effects and structural equation models: direct and indirect effects of disturbance on aboveground C storage

The best multiple linear mixed model showed that aboveground C storage was predicted by stand basal area, individual tree size inequality (CV_DBH), CWM of maximum height (CWM_MH), CWM of wood density (CWM_WD), diversity indices (MNTD and FEve), and disturbance level and accounted for 72% of the variation (Table 2). This result indicates that communities having high stand basal area and tree size inequality, dominated by tall and conservative species and having a high functional evenness, a high mean nearest taxon distance, and little disturbed are more likely to have high aboveground C storage.

The SEM results showed that aboveground C storage was positively affected by stand structure attributes (i.e., CV_DBH) and trait composition (i.e., CWM_MH), after accounting for the strongest negative direct effects of disturbance intensity (Fig. 3; Table 3).The disturbance intensity had a strongest total negative effect on aboveground C storage through the summation of direct (ß = −0.47, P < 0.001) and indirect (ß = −0.17, P < 0.001) negative effects via trait composition (i.e., CWM_MH) (Table 3). The stand structure attribute had a significant direct effect on C storage (ß = 0.10, P = 0.026) but a nonsignificant indirect positive effect via CWM_MH (ß = 0.03, P = 0.164) and FEve (ß = 0.02, P = 0.114). Disturbance intensity had also a negative direct effect on functional evenness but a nonsignificant indirect effect on aboveground C storage via functional evenness, probably due to a nonsignificant direct effect of functional evenness on aboveground C storage (Fig. 3; Table 3).

4 Discussion

In this study, we assessed how multiple metrics of biodiversity (taxonomic, functional, and phylogenetic), stand structural attributes, and trait composition drive aboveground C storage in temperate mixed forests while accounting for the direct, indirect, and total effects of disturbance on C storage. Based on the multiple linear mixed models, we found a strong support for the mass ratio effects based on CWM of maximum height and wood density, as well as support for the niche differentiation and facilitation based on stand structural complexity, functional trait, and phylogenetic diversity indices (Table 2). The topographic variables such as elevation and slope did not retain in the best multiple linear mixed model, suggesting that abiotic versus biotic factors may play a negligible role in driving aboveground C storage. However, our analysis of bivariate relationships revealed that phylogenetic (i.e., MPD and MNTD), as well as functional trait diversity indices (i.e., FEve and FDis), also enhanced aboveground C storage in the studied forests. In addition, we also found that aboveground C storage had negative relationships with Faith’s phylogenetic diversity and species richness but nonsignificant relationships with species evenness and functional richness (Fig. 3). These inconsistent relationships between multiple metrics of biodiversity and aboveground C storage might be attributable to the explanatory power of different metrics of biodiversity to gauge the actual range of strong interactions in the studied forests, rather than an intrinsic ecological mechanism (Yuan et al. 2016; Poorter et al. 2017).

Interestingly, multiple linear mixed model suggested that multiple metrics of biodiversity were weaker predictors of aboveground C storage as compared with functional trait composition (CWM indices) and stand structural attributes (Table 2), suggesting that the mass ratio and stand structural complexity effects outperform biodiversity for driving aboveground C storage. These conflicting results from bivariate relationships and multiple linear mixed model might happen due to the significant differences for predictors and aboveground C storage across levels of disturbance intensity (Fig. 2). Here, we found that functional trait composition (CWM_MH) and stand structure complexity (CV_DBH) decreased with disturbance intensity, and hence aboveground C storage increased with individual tree size inequality probably due to the light availability in the unequal stand dominated by large tree species (Fig. S1).

In this study, SEM results indicated that disturbance had a negative direct effect on CWM of maximum height and aboveground biomass, functional evenness, and individual tree size inequality (Fig. 3). Lower tree height and size inequalities were observed in high and medium disturbance plots, which suggest that past logging activities have decreased the proportion of large individuals. Therefore, low disturbance plots (close to natural old growth forest) are generally dominated by late successional species associated with greater maximum tree height (Chave et al. 2009). Thus, the negligible effects of disturbance on tree DBH diversity implied that the removal of a few large stems did not significantly influence the stand structural complexity (Fig. 2). In accordance with the expectation of intermediate disturbance hypothesis (Connell 1978), the forests with high and medium disturbance levels maintained more species richness than forests with a low disturbance level (Fig. 2). However, a disturbance-induced increase in species richness did not lead to higher functional trait diversity or stand structure complexity, and consequently, the increased tree size inequality as forests recovering from disturbance may result from the intra-specific rather than inter-specific differentiations (Paquette and Messier 2011; Ali et al. 2016; Tobner et al. 2016; Ali and Mattsson 2017a). As such, we found a negative indirect effect of disturbance on aboveground C storage via CWM of maximum height but negligible indirect effects via individual tree size inequality and functional evenness.

The strong positive effect of functional composition on aboveground C storage is consistent with recent several findings that aboveground biomass or productivity is driven by the trait(s) of dominant species in early successional stages of forests or young forests (Kröber et al. 2015; Lin et al. 2016; Prado-Junior et al. 2016; Tobner et al. 2016). This result suggests that large trees contributed disproportionally to stand-level carbon storage compared with small trees (Stephenson et al. 2014; Ali and Yan 2017b). For example, Wu et al. (2015) found that forest biomass is strongly positively correlated with forest height across a wide latitudinal gradient, ranging from tropical to boreal forests. Our previous studies in this area revealed that the relative abundance of larger trees (DBH > 60 cm) was about 1.2% but accounted for about 20% of aboveground biomass (Yuan et al. 2016), further demonstrating the importance of large trees in standing aboveground biomass (Ali and Yan 2017).

Stand structural attributes have been theorized to capture the degree of tree height variation of coexisting individuals among component species, and hence represent stronger drivers of aboveground C storage or productivity (Zhang and Chen 2015; Ali et al. 2016; Danescu et al. 2016). It is plausible that greater tree DBH diversity and individual tree size inequality might promote a vertical filling of available space within a site due to vertically stratified canopy, which in turn allows greater packing densities and enhances aboveground light capture and utilization (Yachi and Loreau 2007; Zhang and Chen 2015; Ali et al. 2016). As an added benefit, structurally complex stands also favor the attenuating fluctuations in temperature and increasing soil moisture, which might explain faster decomposition and nutrient cycling rates in mixed forests by coupling with likely increases in leaf litter production (Schwarz et al. 2014; Crockatt and Bebber 2015). However, other studies found that plots with diverse stand structure, in some cases, might have lower light use efficiency than those plots with uniform stand structure, resulting from the smaller average tree sizes (Binkley et al. 2010; Soares et al. 2017).

In this study, we also found that stand structural attributes outperformed taxonomic, functional trait, and phylogenetic diversity indices for driving aboveground C storage. In accord with the expectation of intermediate disturbance hypothesis, the forest that experienced high and medium disturbance levels maintained higher species richness than a forest with a low disturbance intensity (Fig. 2). However, forests with a low disturbance level (old growth) had lower species richness but higher stand structural diversity and functional trait diversity than highly and medium-disturbed forests (Fig. 2), and hence indicated the stronger influences of the niche complementarity on ecosystem functions due to reduced functional redundancy or competition (Loreau and Hector 2001). Consequently, in contrast with conventional view that inter-specific trait differences and the subsequent vertical stratification as the main drivers of higher productivities in mixtures, our study indicated that increased stand structural complexity may result from intra-specific tree size differentiation rather than inter-specific as forests recovering after disturbance activities (Paquette and Messier 2011; Ali et al. 2016; Tobner et al. 2016). However, recent studies found that stand structural attributes and species diversity was direct and independent drivers of aboveground biomass or productivity in both natural (temperate) forests and managed forests (Danescu et al. 2016; Ali and Mattsson 2017a). In partial contrast, empirical studies found that stand structural complexity act as a mechanism for the positive biodiversity and aboveground biomass relationships in natural forests (Zhang and Chen 2015; Ali et al. 2016). Yet, the interrelationships between multiple metrics of biodiversity and stand structural attributes for driving ecosystem function are still debated, possibly due to the confounding effects of other well-documented factors such as stand age, disturbance, and environmental (abiotic) factors in both natural forests and agroforests as well as in experimental plantations (Zhang and Chen 2015; Ali et al. 2016; Bourdier et al. 2016; Ali and Mattsson 2017a).

5 Concluding remarks

This study provides new evidence that multiple biotic factors, including stand structure complexity, functional trait diversity and composition (CWM of a trait values), and phylogenetic diversity, determine aboveground C storage in temperate mixed forests. However, past logging activities have significantly decreased the proportion of large trees, stand structure complexity, and functional diversity, and hence strongly influenced the aboveground C storage. As such, stand structural complexity and functional trait composition were the independent drivers of aboveground C storage, but stand structural complexity has outperformed the contribution of biodiversity to aboveground C storage. Collectively, from a practical point of view, this study suggests that conserving functional identity of species and complex stand structure would be the alternative choices for maintaining higher aboveground C storage in temperate mixed forests recovering from disturbances.