Effects of Tree Diversity, Functional Composition, and Large Trees on the Aboveground Biomass of an Old-Growth Subtropical Forest in Southern China

Abstract

Forest aboveground biomass (AGB) plays an important role in regulating the global carbon cycle and is thus an essential component of ecosystem functioning. In the relationships between biodiversity and ecosystem functioning (BEF), studies have shown that many biotic factors (e.g., species, functional traits, and large trees) and abiotic factors have significant impacts on AGB. However, the relative strength of these affecting factors remains unclear. In this study, we analyzed woody plants (diameter at breast height [DBH] ≥ 1 cm) within a 1.6 ha plot in an old-growth subtropical natural forest in southern China. We used structural equation models to test the effects of tree diversity (species, phylogenetic, functional, and size inequality), functional composition, large trees, and environmental factors (topography, soil nutrients, and understory light) on AGB. Our results indicated that size inequality, the community-weighted mean of maximum DBH (CWM_MDBH), and large trees had significant, positive effects on AGB (p < 0.001), while lower soil phosphorus content was found to promote an increase in AGB. Furthermore, large trees, which were mostly composed of dominant tree species, were the main driver of AGB, and the effect of functional composition (e.g., CWM_MDBH) on AGB was substantially reduced by large trees. We argue that the selection effect plays a key role in regulating BEF relationships in subtropical natural forests and conclude that retaining large-diameter trees and dominant species, along with sustaining a complex stand structure, are key measures for improving productivity.

1. Introduction

Numerous experimental or monitoring studies have shown that changes in biodiversity can significantly affect ecosystem functioning and stability by increasing resistance and resilience to stresses such as extreme weather events and invasive species [1]. Unprecedented biodiversity losses over the past century may severely threaten key ecosystem functions that humans both benefit from and rely on for sustainable development. Therefore, the relationship between biodiversity and ecosystem functioning (BEF) has become increasingly central in ecological studies [2]. For the reason that forest ecosystems cover most of the earth’s surface and host the majority of terrestrial biodiversity [3], understanding the BEF relationships of forest ecosystems is critical for mitigating climate change and regulating the global carbon cycle. Forest aboveground biomass (AGB), which accounts for approximately 80% of total forest biomass, is the main source of terrestrial carbon stocks and is considered a fundamental indicator of forest ecosystem functioning [4]. However, the BEF relationships of natural forests remain disputed, particularly in subtropical forests where a high diversity of tree species (e.g., seasonal, evergreen, or deciduous; tropical vs. austral origin) often coexist in and co-dominate the same plant community [5]. The effects of biodiversity on AGB may be more complex in subtropical forests than in temperate forests due to the unpredictable effects of tree–tree interactions among a large number of species as well as uncertainty about the effects of abiotic conditions on BEF relationships.

The niche complementarity and selection effects hypotheses are widely used to explain BEF relationships. The niche complementarity hypothesis posits that niche partitioning among species may lead to greater resource use efficiency, thus promoting increased ecosystem stability and productivity (or biomass) [6]. Early studies quantified niche differentiation based on species diversity, assuming that the higher the species diversity, the higher the degree of ecological niche differentiation in a community [7]. Despite the convenience and simplicity of this approach, numerous studies have concluded that niche overlap or dissimilarity can be better captured using phylogenetic and functional traits, suggesting that distantly related or functionally different species compete less for resources and space [8]. For example, studies have shown that phylogenetic and functional diversity influence AGB more strongly than species diversity in forest communities with less abundant tree species [9,10]. While the effects of species diversity on ecosystem productivity have been widely studied in various vegetation types, a range of BEF relationships, including positive, negative, unimodal, or irregular curves, have been observed [11]. Recent studies postulate that tree size structure rather than inter-specific differences may underlie the positive effects of tree diversity on forest community productivity because size diversity promotes light use efficiency and thus may increase AGB by promoting the coexistence of large numbers of different sized trees [12,13,14,15]. Furthermore, species and size structural diversity are highly correlated, particularly in species-rich natural forests [16]. Studies based on the cumulative abundance profile (CAP) of DBH have found that the size structure determines the species richness of the stand [17]. Therefore, the relationship between species diversity and tree size variation should not be considered in isolation when exploring the effects of diversity on ecosystem functioning. An increase in species diversity may also promote tree size heterogeneity at both the local and community scales, implying that the simultaneous influences of different aspects of tree diversity should be considered in biomass forecasting models and that their interrelationships and relative contributions to AGB should be explicitly explored [18]. To date, however, our understanding of how tree diversity affects ecosystem functioning in species-rich forests via the niche complementarity effect is still restricted.

In contrast to the niche complementarity hypothesis, which stresses resource use efficiency and positive interactions among trees, the selection (mass ratio) effect hypothesis highlights the critical role of a few dominant, productive species (or functional traits) in determining community biomass and postulates that an increase in species diversity will increase the probability of the occurrence of productive or high functioning species in the community, particularly in species-rich tropical and subtropical forests [19,20,21]. The community-weighted mean (CWM) of functional traits, i.e., mean trait values weighted by species’ abundances at the plot level, is often used to characterize the strength of selection effects in forest communities [22]. CWM–AGB relationships may be sensitive to the functional traits used, as different traits reflect different resource use strategies [23]. For example, high CWMs for specific leaf area (SLA) or leaf nitrogen or phosphorus concentrations in plant communities are commonly associated with high productivity, as these traits maximize resource acquisition rates and are often associated with acquisitive species [22]. By contrast, high CWMs for leaf dry matter content (LDMC) or wood density may be indicative of low-productivity communities dominated by conservative species with a high tolerance for resource limitations and stressful conditions [24,25]. While niche complementarity and selection effects typically operate simultaneously, the interactions among them remain unclear [26,27].

Large trees are also important direct contributors to AGB in natural forests [28]. Compared with medium-sized and small trees, large trees significantly impact tree size, structural heterogeneity, community functional composition, and the maintenance of species diversity at the local scale [21,29]. However, few studies have explored the indirect influence of large trees on AGB via their effects on tree diversity and functional composition. Moreover, the BEF relationships can be confounded by abiotic conditions since factors such as soil nutrients and light conditions directly influence the availability of resources for tree growth in natural forests, thus influencing diversity [30,31]. As such, the abovementioned processes may jointly determine forest productivity, and the complex relationships between AGB, species diversity, tree size inequality, functional composition, and large trees, along with environmental factors, remain unclear and require further study.

We used plot data from an old-growth subtropical forest to further clarify the BEF relationships in a species-rich forest ecosystem. Such forests are characterized by a large number of mature trees, complex stand structure, rich biodiversity, and carbon storage capacity [32]. We used structural equation models (SEMs), in which explanatory variables were added to the model in a stepwise approach, to explore the following three questions:

(1) How does tree diversity affect AGB? We hypothesize that size inequality is the most crucial factor contributing to AGB and that it weakens the contribution of other tree diversity indices to AGB. To this end, we constructed a SEM that included species, phylogenetic and functional diversity, and environmental factors as explanatory variables (Figure 1A). We then fitted a second model (Figure 1B), which included size inequality based on the first model, and compared the effect sizes of tree diversity and the AGB interpretability between the two models.

(2) What role do niche complementarity and selection effects play in BEF relationships? Here, we consider the effect of tree diversity on AGB as niche complementarity, while the effect of functional composition is a selection effect. We selected the tree diversity index with the highest explanatory power from the second model (Figure 1B) and then added functional composition to construct a new SEM (Figure 1C). We sought to clarify the relative strength of the two effects (niche complementarity and selection effect) by comparing the standardized path coefficients.

(3) How does the presence of large trees influence AGB either directly or indirectly through tree diversity and functional composition? We added the effects of large trees to the third model and identified the significance of tree diversity and functional composition (Figure 1D). This further helped us clarify the relative strength of the complementary niche effect and the selection effect on AGB.

2. Materials and Methods

2.1. Study Site

The study site was located in Yachang National Nature Reserve (24°44′16″—24°53′58″ N, 106°11′31″—106°27′04″ E) in Guangxi Zhuang Autonomous Region, southern China (Figure 2). The reserve encompasses 22,062 ha, and historically, it belonged to the state-owned Yachang Forest Farm. The reserve starts at the southeast edge of the Yunnan-Guizhou Plateau, a transitional zone between plateau and hilly regions. Regionally important rivers, including the Beipan, Nanpan, and Hongshui, flow through the area, and the reserve is geographically unique and ecologically important. The site is characterized by a subtropical monsoon climate with distinct dry and wet seasons. The minimum, maximum, and annual average temperatures are 5.3 °C, 34 °C, 16 °C, and 23 °C, respectively, and the annual rainfall and annual evaporation are 1058 mm and 1484.7 mm, respectively. Precipitation falls on 210 days per year on average, and approximately 50%–60% of total precipitation falls in the summer. The average annual relative humidity is 82%. Soils are dominated by brown laterites in the Hongshui River valley below altitudes of 500 m, whereas mountain red soils occur at altitudes of 500–1000 m, and mountain yellow soils are common above 1000 m. The vegetation is characterized by the typical evergreen broad-leaved forests of the southern subtropical zone [33].

The sample plot was established in the Jiulong Forest Farm of Yachang National Natural Reserve. In July 2020, we used a total station instrument (Southern Mapping Company, Johannesburg, South China; precision: 2 mm + 2 ppm) and established a 1.6 ha plot (80 m × 200 m, datum point: 24°38′54″ N, 106°20′25″ E, altitude, 1700 m). We divided the plot into 40 small quadrats (20 m × 20 m), in which we georeferenced and labeled all woody plants with a diameter at breast height (DBH) ≥1 cm. We tagged each plant and recorded the species, DBH, and height. A total of 10,083 woody plants (including branches), representing 33 families, 55 genera, and 84 species (APG IV), were recorded in the plot. The dominant species were Liquidambar formosana Hance, Rhododendron cavaleriei H.Lév., Eurya impressinervis Kobuski, Eurya nitida Korth., and Ilex ficoidea Hemsl. et al. (Table S2).

2.2. Data Sources

2.2.1. Trait Data

We collected data on 10 key functional traits: leaf thickness (LT), chlorophyll relative content (SPAD), specific leaf area (SLA), leaf dry matter content (LDMC), leaf total nitrogen content (LNC), leaf organic carbon content (LCC), leaf total phosphorus content (LPC), branch wood density (WD), maximum DBH (MDBH), and maximum tree height (MH). These traits reflect plant life history strategies and substantially influence growth [34]. In August 2021, we collected and measured the leaf and branch traits of 76 species, following the protocols in a handbook of standardized approaches for measuring plant functional traits [35]. We randomly selected at least 10 mature, healthy individuals of each species, or 3–5 individuals for rare species (i.e., species represented by <10 individuals). We collected at least 20 mature, sunward leaves and 3 branches from each plant and sought to distribute the sampling throughout all the small quadrats to the greatest degree possible. We measured a total of 10,860 leaves and 1629 branches from 543 plants.

In the field, we wrapped the samples in wet paper towels and placed them in sealed plastic bags. Upon our return to the laboratory, we soaked the branches in water for at least 2 h, and measured the basic physical properties of the leaves within 12 h. Leaves were scanned using a scanner (Epson, Suwa, Nagano, Japan), and analyzed using ImageJ 18.0 software. LT was measured using micrometers (San Liang, JDE03, Dongguan, China). We measured the tip, middle, and base of each leaf and used the average value to represent LT. A SPAD-502 chlorophyll meter was used to determine SPAD. The saturated fresh weight of leaves was measured using a thousandth electron balance. We dried the leaves in an 80 °C oven for 48 h to a constant weight and recorded their dry weights. SLA and LDMC were calculated based on Formulas (1) and (2). The dried leaves were ground to a powder prior to measuring LNC, LCC, and LPC. We soaked branches to saturation and measured wood volume using the drainage method. Branch dry weight was determined after samples were oven-dried to a constant weight. WD was calculated based on Formula (3). The height of each tree was measured using a Haglöf Vertex IV (Haglöf Sweden AB, Långsele, Sweden), and the maximum height of each species was determined based on the Flora of Guangxi.

SLA = Leaf area/Leaf dry weight

(1)

LDMC = Leaf dry weight/Leaf fresh weight

(2)

WD = Branches dry weight/Wood volume

(3)

where Leaf area is in cm², Leaf dry weight, Leaf fresh weight and Branches dry weight is in g, Wood volume is in·cm³.

2.2.2. Environmental Variables

To calculate topographic factors (i.e., elevation, slope, aspect, and convexity), we created a digital elevation model (DEM) using ArcGIS 10.2. The DEM was based on the three-dimensional coordinates of the geo-referenced trees and had a spatial resolution of 0.8 m [36]. We then used this model to calculate the topographic attributes of each small quadrat. Elevation was determined based on the average elevation values of the four vertices in each quadrat. Convexity was calculated based on the elevation of each quadrat minus the mean elevation of the eight neighboring quadrats. The slope was assessed based on the mean angular deviation from the horizontal plane of each of the four triangular planes formed by connecting three corners of the quadrat [37]. Finally, the aspect was determined based on the mean angle between the true north and the orientation of the same four planes [38].

Five soil samples, ranging in depth from 0–20 cm, were collected along the diagonal and at the center of each quadrat. Soil samples were used to measure the following physical and chemical soil properties: PH, total nitrogen, organic carbon, total phosphorus, total potassium, alkaline hydrolysis nitrogen, available phosphorus, available potassium, and soil moisture content. To reduce the number of soil variables in our analysis and to avoid strong correlations among variables, we standardized each variable (mean = 0, standard deviation = 1) and used principal component analysis (PCA) to reduce the dimensionality of the data. The first PCA axis (soil PCA1, 47.85%) was driven primarily by pH and total potassium. The second axis (soil PCA2, 28.69%) was driven primarily by total phosphorus and available soil nutrients (alkaline hydrolysis nitrogen, available phosphorus, and available potassium). These two axes were used to represent soil physicochemical properties (

Table 1. The loading values of soil factors on the first two principal components.

Soil Factors	Soil PCA1	Soil PCA2
Soil PH (PH)	0.39	0.24
Soil total nitrogen (TN)	−0.44	0.17
Soil organic carbon (OC)	−0.45	0.11
Soil total phosphorus (TP)	0.08	0.57
Soil total potassium (TK)	0.33	0.38
Soil alkaline hydrolysis (AHN)	−0.41	0.25
Soil available phosphorus (AP)	−0.14	0.27
Soil available potassium (AK)	−0.07	0.54
Soil moisture content (SMC)	−0.38	−0.09
Explained variance proportion	47.85%	28.69%
Cumulative proportion	47.85%	76.54%

).

The understory light conditions in each quadrat were measured using an EOS 80D single lens reflex camera (Canon, Tokyo, Japan) with a Sigma 4.5 mm F2.8EXDC fish-eye lens converter (Sigma-Aldrich, St. Louis, MO, USA). We positioned the camera 2 m above the ground using a tripod and rotated the tripod to set the magnetic needle to the north. We took three hemispherical canopy images at each quarter point and the midpoint on the diagonal line of each quadrat, for a total of 600 images. Photographs were taken at sunrise or sunset on sunny days to avoid direct sunlight [39]. We used Hemiview2.1 software to analyze the images and calculate canopy opening (Co) and leaf area index (LAI) [40]. The understory light of each subplot was the average LAI or Co from 15 images taken within a quadrat.

2.2.3. AGB Calculation

Following Chave et al. [41], we used the general allometric equation model to calculate the AGB of all woody plants (DBH ≥ 1 cm) in the sample plot.

$$\mathrm{A}\mathrm{G}\mathrm{B}=0.0673\times {\left(\rho {\mathrm{D}}^2\mathrm{H}\right)}^{0.976}$$

(4)

where D is in cm, H is in m, and $\rho$ is in g·cm⁻³.

2.2.4. Multivariate Diversity Metrics

• Species diversity;

Species richness (SR), the Shannon index, and the Simpson index were calculated for measuring species diversity [42,43].

$$Species richness=S$$

(5)

$$Shannon=-\sum P_{i}ln{P}_i$$

(6)

$$Simpson=1-\sum P_i^2$$

(7)

where S is the number of species and P_i is the proportion of species abundance belong to total individuals. These indices were calculated using the R package “vegan”.

• Phylogenetic diversity;

Faith’s phylogenetic diversity (PD), mean pairwise distance (MPD), and mean nearest taxon distance (MNTD) were used as proxies for phylogenetic diversity. PD is calculated based on the sum of phylogenetic branch lengths weighted by species abundance [44]. We generated a phylogenetic tree (APG IV) and calculated the PD of each quadrat using the R package “V.PhyloMaker” [45]. MPD represents the average phylogenetic distance between all individuals in a community, and MNTD is calculated based on the average phylogenetic distance between individuals and their closest non-conspecific relative [46]. The phylogenetic indices were calculated using the R package “vegan” [47].

• Functional composition;

Functional composition was quantified based on the community-weighted means of the ten functional traits [48].

$$CWM=\sum _{i=1}^n{A}_{ij}{T}_{ij}$$

(8)

where A_ij is the abundance of species i in subplot j, and T_ij is the mean trait value of species i in subplot j. The number of species on the subplot was measured by n.

• Functional diversity;

Functional diversity was quantified using Functional richness (FR_ci), Functional evenness (FE_ve), Functional divergence (FD_iv) and Functional dispersion (FD_is) [8,49].

$${FR}_{ci}=\frac{{SF}_{ci}}{R_c}$$

(9)

$${FE}_{ve}=\frac{\sum _{i=1}^{S-1}min\left({PEW}_i,\frac{1}{S-1}\right)-\frac{1}{S-1}}{1-\frac{1}{S-1}}$$

(10)

$${FD}_{iv}=\frac{2}{\pi }arctan\left\{5\times \sum _{i-1}^N\left[{\left(\ln C_i-\overline{\ln x}\right)}^2\times A_i\right]\right\}$$

(11)

$${FD}_{is}=\frac{\sum a_j{z}_j}{a_j}$$

(12)

where FR_ci is the functional richness of plant functional trait c in community i, SF_ci is the niche space occupied by species in the community i, and R_c is the absolute value range of plant functional trait c. FE_ve describes the evenness of abundance distribution in a functional trait space [8], S is the species richness, and PEW_i is the weighted evenness of species i. FD_iv is a functional separation index that contains multiple functional traits, C_i is the value of the ith functional trait, Inx is the natural logarithm of traits, A_i is the abundance proportion of the ith functional trait, and N is the number of species in the community. FD_is is the average distance from the weighted abundance of species to the centroid of all species in the community, a_j is the relative abundance of species j, z_j is the distance from species j to the weighted centroid. Functional composition and functional diversity indices were calculated using the R package “FD” [50].

• Size structure diversity;

Tree size inequality was quantified using the coefficient of DBH variation (CV) and the Gini coefficient (Gini) of individual trees.

$$CV=100\%\frac{\sqrt{\frac{1}{N}{({DBH}_k-\mu )}^2}}{\mu }$$

(13)

$$Gini=\frac{{\sum }_{i=1}^N{\sum }_{j=1}^N\left|x_i-x_j\right|}{2N^2\overline{x}}$$

(14)

where DBH_k is the DBH of the kth individual in a subplot, μ is the average DBH of all individuals in a subplot. N is the total number of individuals in a subplot, $\left|x_i-x_j\right|$ is the absolute value of the DBH difference between any two trees, and $\overline{x}$ is the average DBH in each subplot. The CV index was calculated using the R package “raster”, and the Gini index was calculated using the R package “ineq” [51,52].

2.2.5. Defining Large Trees

Criteria for defining large-diameter trees may vary by forest type and successional stage [22]. Lutz et al. [28] proposed three methods for identifying large trees. The first is to define the largest 1% of trees in a community as large trees, based on DBH. The second is based on a fixed threshold (DBH ≥ 60 cm), and the third defines large trees based on 50% cumulative biomass. We applied the first method to our data using thresholds of 1%, 2%, and 3%, yielding large tree datasets comprising 98, 196, and 294 individuals, respectively. The DBHs of these groups were concentrated in the range of 40–50 cm, 30–40 cm, and 25–35 cm, respectively. The most abundant species included in all three groups were Liquidambar formosana, Schima argentea, Clethra kaipoensis and Rhododendron cavaleriei (Figure S1). We defined the largest 2% of trees, based on DBH, as large trees, as few individuals had DBHs ≥ 60 cm. However, we also tested the 1% and 3% thresholds in our AGB models (Figure S5).

2.3. Data Analysis

We first used a generalized least squares (GLS) model to assess spatial autocorrelation among neighboring quadrats [53]. The spatial model with the lowest Akaike information criterion (AIC) score was selected from among five candidate spatial autocorrelation models. We then compared this model to a non-spatial GLS model. We found that the model without spatial autocorrelation consistently had a lower AIC score (Table S3), indicating an absence of spatial autocorrelation among the quadrats. We then used linear regressions to assess the relationships between AGB and each explanatory variable. All data were log-transformed and normalized (mean = 0, standard deviation = 1) prior to analysis. We fitted linear models for AGB and each explanatory variable to obtain standardized coefficients. Forest plots of the normalized coefficients for each variable were generated using the R package “ggplot2” [54].

Due to limitations in the relationship between sample size and path coefficients in SEMs analysis, it is necessary to screen explanatory variables. Random Forest is a machine-learning algorithm that is relatively insensitive to multicollinearity and overfitting, allowing for the inclusion of numerous explanatory variables [18]. We used the random forest to assess the importance of each explanatory variable and selected the 15 most important abiotic and biotic variables that most strongly influenced AGB (Figure S3) for correlation analysis. The importance of the selected variables was quantified based on the percent increase in mean square error. Random Forest models were run in the R package “randomForest” [55]. To avoid overfitting of the SEMs due to high collinearity among variables, we performed a correlation analysis of the 15 selected variables. If the correlation coefficient of a pair of candidate variables was >0.70, the first variable was retained and the second variable was excluded (Figure S4). The remaining variables were included in the SEMs. In addition, to evaluate the conciseness and suitability of the selected variables, we used the dredge function in the R package “MuMIn” [56] to generate a full set of models based on maximum likelihood estimation and then used AIC to select the best model (Table S4). Models were considered equally supported if the difference in AIC (delta AIC) scores was <2.

Finally, we used linear SEMs that gradually added variables to assess the effects of the biotic and abiotic variables on AGB based on a priori conceptual models (Figure 1A–D). Standardized path coefficients can be used as a measure of the sensitivity of response variables to explanatory variables [57]. To improve model fit, we conducted a stepwise removal of non-significant paths [58]. The best-fit SEM was identified based on several parameters, including the chi-square test (${\chi }^2$), Akaike information criterion (AIC), mean square and square root of asymptotic residuals (RMSEA), comparative fit index (CFI), and standardized root mean square residuals (SRMR) [59]. The SEMs were fitted using the “lavaan” package [60]. To aid in model interpretation, we calculated the direct, indirect, and total standardized effect sizes of each explanatory variable. All calculations were conducted in R4.1.1 [61].

3. Results

3.1. Bivariate Relationships between AGB and Biotic Variables

There was no significant correlation observed between AGB and species diversity or phylogenetic diversity. Notably, the community-weighted mean of leaf phosphorus concentration (CWM_P) and functional dispersion (FDis) exhibited significant negative impacts on AGB. Conversely, community-weighted means of maximum breast-height diameter (CWM_MDBH), Gini index, coefficient of variation (CV), and large trees were associated with significant positive impacts on AGB. The relative effect sizes were ranked in the following order: Large trees > Gini > CV > CWM_MDBH (Figure 3).

3.2. The Effects of Biotic and Abiotic Factors on AGB

The random forest analysis indicated that, among tree diversity indices, FDis impacted AGB more than species diversity (Figure S3A). However, size structural diversity was far more important than other diversity indices if tree size inequality was considered in the model (Figure S3B). Among functional compositions for different traits, CWM_MDBH had the second highest importance value following tree size inequality (Figure S3C). Large trees had the strongest effect among the variables included in the full model (Figure S3D,F). Strong correlations were detected between the diversity indices as well as between the functional composition indices (Figure S4). By contrast, neither large trees nor tree size inequality were strongly correlated with the remaining explanatory variables (correlation coefficient < 0.7).

The final SEMs showed that species, phylogenetic, and functional diversity had no significant effect on AGB (Figure 4A). Rather, 7.9% of the AGB variation was explained by the mean nearest taxon distance (MNTD), functional divergence (FDiv), and functional dispersion (FDis). The inclusion of tree size inequality increased the explanatory power of the model to 47.8% (Figure 4B). Among them, Gini had the strongest positive direct effect on AGB (β = 0.679, p < 0.001), but Simpson’s index had a significant negative total effect on AGB (Figure 5B). We observed a weak relationship between species diversity and tree size inequality. Soil nutrients had direct and significant effects on tree diversity but no significant effects on AGB (Figure 4A,B). When both the functional composition and tree size inequality were considered, the explanatory power of the model reached 57.6% (Figure 4C). CWM_MDBH (β = 0.484, p < 0.001) and CV (β = 0.679, p < 0.001) had strong, positive, direct effects on AGB, whereas CWM_LDMC (β = −0.454, p < 0.001) had significant, negative effects. Soil PCA2 had a significantly positive direct effect on AGB (β = 0.383, p < 0.01) but no total effect on AGB (Figure 4C and Figure 5C).

Following the addition of large trees, the model explained >60% of the variation in AGB (Figure 4D and Figure S5A,B). Large trees not only had significant and positive effects on AGB but also had a strong impact on functional composition-AGB or tree size inequality-AGB relationships. The effects of functional composition and tree size inequality on AGB weaken following the addition of large trees. As the threshold for large trees increased, the direct and total effects of CWM_MDBH on AGB weakened, whereas the direct and total effects of tree size inequality became increasingly strong (Figure 5D and Figure S5C,D).

4. Discussion

Assessing the main drivers of AGB in natural forests is critical for forest management, biodiversity conservation, and forest carbon sequestration [62]. Our analysis demonstrates that different factors, including tree diversity, functional composition, large trees, and environmental variables, jointly influence AGB in an old-growth subtropical natural forest. Our results confirm that CWM_MDBH, tree size inequality, and large trees significantly impact AGB, thus providing strong evidence of the importance of both selection effects and niche complementarity processes in BEF relationships. However, a large proportion of the AGB variance is explained by tree size inequality and large trees jointly. Notably, large trees completely weakened the effects of CWM_DBH on AGB, indicating that selection effects play a stronger role in regulating AGB formation in natural subtropical forests.

4.1. The Effects of Tree Size Inequality on AGB Were Stronger Than Those of Other Tree Diversity

Our model indicates that tree size inequality positively affects AGB, which is consistent with the results of numerous other studies conducted in subtropical forests [63]. Size inequality directly reflects the horizontal variations of individual trees while also indirectly indicating differences in habitat requirements or resource demands, particularly for light, among species [64]. Therefore, complex size structure improves forest canopy growth and stratification and enhances niche complementarity through efficient light capture or spatial stratification [65]. Although numerous studies have shown that AGB in subtropical forests increases with increasing biodiversity [65,66,67], our results suggest that species, phylogenetic, and functional diversity have negligible, or even negative, effects on AGB. One possible explanation for this may be that the competitive effects of dominant species inhibit resource utilization by other species, resulting in decreased species diversity as productivity increases [68,69]. Therefore, asymmetric inter-specific competition for available resources may result in a negligible relationship between species diversity and AGB. Furthermore, our results were consistent with other studies based on global forest databases that indicate that the effects of functional diversity on AGB are negligible [21]. Our focal community was in the late stages of succession, which is frequently characterized by niche saturation or functional redundancy, which in turn may have led to the observed weak association between AGB and functional diversity [26,70]. In theory, size structure diversity is a key mechanism underlying the positive relationship between species diversity and AGB [65]; however, our results show that the effects of tree size inequality on species diversity are not significant. Ren et al. [62] came to a similar conclusion in a study conducted in a natural subtropical forest. In summary, our results emphasize the importance of size structure in forest monitoring and management, which is consistent with the findings of Angiolini et al., who supported the role of size structure in classifying European forest habitat types [71].

4.2. The Relative Importance of Selection Effect and Niche Complementary Effect

While studies have shown that niche complementarity and mass ratio effects both drive AGB, the relative importance of the two mechanisms in natural forests remains disputed [72]. We found that both stand-level DBH variations and the functional component of maximum DBH had significant, positive effects on AGB, suggesting that dominant species with potentially larger DBHs and greater size variation drive AGB at our study site. Thus, our study demonstrates the presence of both niche complementarity and mass ratio effects. In models that only included size structural diversity and functional composition (Figure 4C, Figure 5C and Figure S3C), CWM_MDBH contributed less to AGB than size structural diversity (CV). Based on these findings, it can be preliminarily inferred that niche complementarity is more important than mass effects in determining AGB in natural subtropical forests. Many other studies have used the CWM of maximum tree height (CWM_MH), which reflects competition for light, to explore the relationship between functional composition and AGB and have demonstrated that CWM_MH is the main driver of AGB [27,73]. We opted to use DBH instead, as measurements of tree height are prone to substantial error and DBH is a better indicator of competitive advantages among plants [74]. In addition, we observed a negative relationship between AGB and the conservative trait CWM_LDMC. This may be because slow-growing, conservative species (e.g., species located in understory strata) are often suppressed by functionally dominant species [22].

4.3. The Role of Large Trees on AGB Reflects the Selection Effect

It is generally acknowledged that large trees contribute substantially to AGB, but little is known about their unique role in driving AGB in forest communities relative to other biotic and abiotic factors [22,75]. Our models demonstrate that among AGB, large trees, size inequality, and functional composition, large trees have the greatest impact on AGB, and as the threshold for large trees increases, so does the contribution of large trees to AGB (Figure 4D and Figure S5A,B). However, large trees also reduce the direct positive effects of size-structural diversity and CWM_MDBH on AGB, which is consistent with the results of previous studies on stand structural complexity and aboveground carbon stocks [76]. Therefore, it is possible that the effects of biodiversity on AGB may have been largely misattributed to stand structural diversity in previous studies [15,62]. However, the significant effect of tree size inequality was not completely weakened by large trees. Interestingly, as the threshold for defining large trees increased, the effect of diameter-at-breast-height dominant species on AGB gradually weakened and eventually became insignificant, while size structural diversity became increasingly important. This may be due to the fact that higher thresholds of large trees capture more dominant species, which account for most of the AGB (Figure S2). Thus, the effects of large trees on AGB belong to the selection effect to some extent. Furthermore, the relative importance of niche complementarity and mass ratio effects may differ among strata. Ali et al. [23] found that AGB was driven by conservative species in the understory but by dominant species in the canopy. Stratification may be an important direction for future research to further explore the relationship between biodiversity and ecosystem functioning.

4.4. The Direct and Indirect Effects of Environment Conditions on AGB

Compared with the bivariate BEF relationships (Figure 3), the significance of functional diversity on AGB changed when both biotic and abiotic factors were included in our models (Figure 4). This suggests that environmental factors have non-negligible effects on BEF relationships. Among the environmental factors included in our study (topography, soil nutrient, and understory light), only soil PCA2 had a positive direct effect on AGB; this variable also had significant positive effects on species, phylogenetic, and functional diversity. These positive effects may be attributable to soil nutrient levels. Other studies have shown that forest stands with nutrient-poor soils exhibit stronger species diversity effects and lower tree size variation [77,78,79,80,81]. Total phosphorus, which was generally low at our site, had the strongest contribution to soil PCA2, suggesting that phosphorus use by plants promoted niche differentiation (Table 1 and Table S1). In addition, soil PCA2 had a positive effect on CWM_LDMC, supporting the hypothesis that nutrient-poor soils are advantageous to species with conservative life history strategies [22,78]. By contrast, lower pH or potassium levels may limit the availability of soil nutrients [79]. In our study, soil PCA1 mostly reflected variations in pH and total potassium. Moreover, the soils at our study site were acidic, and the total potassium content was lower than average for evergreen broad-leaved forests in southern China [82]; therefore, soil PCA1 had a negative effect on tree size inequality.

5. Conclusions

Our study demonstrated that niche complementarity and mass ratio effects jointly maintain BEF relationships in subtropical forests, but that the relative importance of the selection effect was greater. Furthermore, tree size inequality exerted stronger positive effects on AGB, independent of species, phylogenetic, or functional diversity. Species with large DBH inhibited the contribution of conservative species to AGB. More importantly, large trees were the main drivers of AGB as well as a key driver of size diversity and functional composition. Concurrently, lower soil phosphorus content promoted an increase in the AGB of conservative species. Therefore, the retention of large-diameter trees and the maintenance of size structural complexity are conducive to improving productivity in subtropical natural forests.