Stand Structure and Abiotic Factors Modulate Karst Forest Biomass in Southwest China

: Understanding the driving factors of forest biomass are critical for further understanding the forest carbon cycle and carbon storage management in karst forests. This study aimed to investigate the distribution of forest aboveground biomass (AGB) and the e ﬀ ects of stand structural and abiotic factors on AGB in karst forests in Southwest China. We established a 25 ha plot and sampled all trees ( ≥ 1 cm diameter) in a subtropical mixed evergreen–deciduous broadleaf forest. We mapped the forest biomass distribution and applied a variation of partitioning analysis to examine the topographic, stand structural, and spatial factors. Furthermore, we used structural equation models (SEM) to test how these variables directly and / or indirectly a ﬀ ect AGB. The average AGB of the 25 ha plot was 73.92 Mg / ha, but that varied from 3.22 to 198.11 Mg / ha in the 20 m × 20 m quadrats. Topographic, stand structural, and spatial factors together explained 67.7% of the variation in AGB distribution. The structural variables (including tree density and the diameter at breast height (DBH) diversity) and topographic factors (including elevation, VDCN (vertical distance to channel network), convexity, and slope) were the most crucial driving factors of AGB in the karst forests. Structural equation models indicated that elevation, tree density, and DBH diversity directly a ﬀ ected AGB, and elevation also indirectly a ﬀ ected AGB through tree density and DBH diversity. Meanwhile, AGB was indirectly inﬂuenced by VDCN, convexity, and slope. The evaluation of stand structural and abiotic drivers of AGB provides better insights into the mechanisms that play a role in carbon storage in karst forests, which may assist in improving forest carbon management. a karst forest biomass and examined stand structural and abiotic factors driving AGB in a 25 ha dynamic forest plot in Southwest China, which provided important information on biomass distribution patterns for a better understanding of karst ecosystems. We found that the mean forest AGB was 73.92 ± 23.51 Mg / ha. The structural, topographic, and spatial factors predominantly contributed to the distribution of AGB. The stand structural complexity (tree density and DBH diversity) and topographic properties (elevation, VDCN, convexity, and slope) most strongly inﬂuenced AGB. The SEM suggested that AGB is directly driven by tree density, DBH diversity, and elevation. In addition, topographic factors mainly had indirect inﬂuence on AGB. Our work provides the links to understanding and establishing relationships between AGB and the driving factors on karst forests. These results will contribute to provide basic data for global and regional vegetation and carbon models. Further, regional forest carbon models need to consider community structural and topographic parameters.


Introduction
Forests account for 30% of the earth's total land area, and contribute to 75% of the terrestrial gross primary production and 80% of the total global plant biomass [1,2]. Previous studies have shown that forests play important roles in the carbon cycle and reducing the greenhouse effect by sequestering carbon dioxide [3][4][5]. Southwest China has the second largest forest region in China [6]. This region also contains the largest karst forest region in China, with an extremely fragile geological background, small environmental carrying capacity, and low tolerance to disturbance [7]. During the past century, a large part of the karst region in southwest China was severely degraded following the destruction of natural vegetation due to human disturbance. Thankfully, most of the degraded land has been undergoing ecological restoration, either through natural regeneration or afforestation, attributed to the implementation of ecological restoration projects (e.g., the "Grain for Green" project, the Karst Rocky Desertification Restoration Project) over the past two decades [8,9]. Therefore, a detailed study on karst forest biomass is important for improving the accuracy of estimating Chinese vegetation carbon storages and sequestration potential.
Southwest China contains one of the largest karst regions in the world [8]. Previous studies have assessed the biomass in this region on small and large scales. For example, the largest biomass increase in the world occurred in the South China karst region between 1999 and 2012 [10]. Other studies reported the biomass in typical karst ecosystems, or along vegetation restoration areas based on surveys of small plots [2,6,[11][12][13][14]. However, no biomass investigations have been conducted on a local scale in the karst region, especially in forests. Mapping forest biomass patterns with accurate spatial information for each recorded individual tree can specifically describe biomass variance along a topographic gradient and reveal the fine-scale relationships of potential factors affecting the distribution [15]. This is beneficial for improving the accuracy of the forest ecosystem carbon cycle model.
Forest biomass is a complex property affected by site conditions, forest structure, and ecological processes [1,16,17]. Abiotic factors, such as topography and edaphic conditions, greatly influence the forest biomass because they determine resource availability for plant growth and survival [18]. Topography, as a key factor, maintains the heterogeneity of forest landscapes. Topographic features (e.g., slope, curvature, and terrain relief) greatly influence local scale variation in soil chemistry, hydrology, and microclimate [19]. Not only abiotic factors, but also biotic factors can determine forest biomass. For instance, forest structural attributes (e.g., tree size variance and stem density) play a vital role in regulating the distribution of forest biomass [20]. Stand structural complexity can enhance aboveground biomass directly in tropical forests, and may be a result of the high efficient utilization of light, water, and soil nutrient resources under the action of niche differentiation and facilitation [21]. Additionally, higher stem density may increase tree biomass via greater canopy packing which leads to more light capture [22]. Similarly, tree size variation is linked to biomass in natural forests through better spatial distribution of different tree crowns, and thus, there is more efficient utilization of light [23,24]. However, the influence of abiotic and biotic factors on karst forest biomass has not been well examined.
The objectives of this study were: (1) To estimate the aboveground biomass (AGB) and map the spatial pattern of AGB in the Mulun 25 ha plot; (2) to quantify the contributions of topographic, stand structural, and spatial factors on biomass spatial distribution and select the main driving factors; and (3) to evaluate the direct and indirect effects of the main driving factors on AGB in the karst forest.

Study Site
Mulun National Natural Reserve (MNNR) and the adjacent Maolan National Natural Reserve in the Guangxi and Guizhou provinces, respectively, preserve the last undisturbed remnants of karst forests in Southwest China, and possibly in the world [7]. Our study was conducted at MNNR (107 • 54 01 -108 • 05 51 E, 25 • 07 01 -25 • 12 22 N), in the Huanjiang County, northwestern Guangxi Province, China. MNNR was set up in 1991 to protect the subtropical mixed evergreen-deciduous broadleaf forest ecosystem that developed on a limestone substrate. The topography in this reserve is characterized by steep hills separated by lowland depressions with numerous potholes and caves, and extensive underground streams [7]. The average annual temperature is 19.4 • C and the average annual precipitation is 1500 mm in this area [7].

Predictor Variables
The Mulun plot was divided into a grid composed of 625 20 m × 20 m quadrats. For each quadrat, seven topographic variables were recorded: elevation, convexity, slope, aspect, rock outcrop ratio (ROR), topographic wetness index (TWI), and vertical distance to channel network (VDCN). The first four topographic factors were computed following the method in Du et al. [7]. ROR is a ratio of the rock coverage area. TWI was calculated as the ratio of the area upslope from any given point on the landscape to the local slope at that point. VDCN was calculated as the vertical distance from the channel network. TWI and VDCN are commonly used to quantify topographical control on hydrological processes [25,26].
where Hd is Shannon's diversity of DBH, and pj is the relative basal area of the jth tree DBH classes in the given plot. We selected tree density and DBH diversity as the stand structure variables to examine the relationships with AGB.

Statistical Analysis
For the estimation of AGB of the Mulun plot, we used a mixed allometric equation based on DBH (cm) to calculate the AGB of the trees following Wang et al. [27]. The relative contribution of structure variables and topography in the spatial pattern of AGB were determined by variance partitioning based on multiple regressions [28]. The principal coordinate neighbor matrix (PCNM) eigenfunctions on a 20 m × 20 m scale were computed by performing a principal coordinate analysis

Predictor Variables
The Mulun plot was divided into a grid composed of 625 20 m × 20 m quadrats. For each quadrat, seven topographic variables were recorded: elevation, convexity, slope, aspect, rock outcrop ratio (ROR), topographic wetness index (TWI), and vertical distance to channel network (VDCN). The first four topographic factors were computed following the method in Du et al. [7]. ROR is a ratio of the rock coverage area. TWI was calculated as the ratio of the area upslope from any given point on the landscape to the local slope at that point. VDCN was calculated as the vertical distance from the channel network. TWI and VDCN are commonly used to quantify topographical control on hydrological processes [25,26].
where Hd is Shannon's diversity of DBH, and p i is the relative basal area of the ith tree DBH classes in the given plot. We selected tree density and DBH diversity as the stand structure variables to examine the relationships with AGB.

Statistical Analysis
For the estimation of AGB of the Mulun plot, we used a mixed allometric equation based on DBH (cm) to calculate the AGB of the trees following Wang et al. [27]. The relative contribution of structure variables and topography in the spatial pattern of AGB were determined by variance partitioning based on multiple regressions [28]. The principal coordinate neighbor matrix (PCNM) eigenfunctions on a 20 m × 20 m scale were computed by performing a principal coordinate analysis based on a truncated Euclidean distance matrix [29]. The negative eigenvalues were filtered, and the positive values were used as predictors. Forward selection was applied to select the significant (p < 0.05 after 999 simulations) structure variables, topographic, and PCNMs. The adjusted R 2 was used to explain pure effects and joint effects for each group of factors [15].
In order to assess the complex relationships among the variables, structural equation modeling (SEM) was used to explicitly test the direct and indirect effects of a combination of factors on the AGB including all significant variables from the PCNM analysis. Here, we designed a hypothetical model in which the topographic factors and structural variables could directly influence the AGB. The topographic factors could also directly influence tree density and DBH diversity. These correlations are supported by previous studies [21,30]. We selected the best SEMs based on the fit indices, including a nonsignificant Chi-square (χ 2 ) test statistic (p > 0.05), S SRMR < 0.08, and both GFI and CHI > 0.95 > through removing the nonsignificant variables with the lowest Akaike information criterion (AIC) [21,31].
All statistical analyses were performed in R 3.5.1. The predicted map of AGB was depicted using ordinary kriging in the geoR package [32]. The PCNM analysis was conducted using the PCNM package [33]. Variation partitioning analysis was computed using the 'varpart' function in the vegan package [34]. The SEM analysis was performed using the lavaan' package [35].

Results
According to the first census, the plot had 110,203 individuals belonging to 61 families, 147 genera, and 227 species. The average DBH of trees was 4.48 cm. The mean stem density was 4408 individuals per ha and the mean basal area was 17.00 m 2 /ha. The average forest AGB was 73.92 ± 23.51 Mg/ha (± standard deviation) in the 20 m × 20 m quadrats of the Mulun plot, although AGB ranged widely from 3.22 to 198.11 Mg/ha. AGB among different diameter ranges was submitted to a normal distribution (Shapiro-Wilk test; p = 0.0526), and the biomass of trees in the 5-10 cm DBH range was higher than in the other diameter classes (Figure 2). Trees from 1 to 5 cm DBH accounted for only 14.89% of the biomass, although they accounted for 72.22% of the stems.
Forests 2020, 11, x; doi: FOR PEER REVIEW www.mdpi.com/journal/forests based on a truncated Euclidean distance matrix [29]. The negative eigenvalues were filtered, and the positive values were used as predictors. Forward selection was applied to select the significant (p﹤0.05 after 999 simulations) structure variables, topographic, and PCNMs. The adjusted R 2 was used to explain pure effects and joint effects for each group of factors [15]. In order to assess the complex relationships among the variables, structural equation modeling (SEM) was used to explicitly test the direct and indirect effects of a combination of factors on the AGB including all significant variables from the PCNM analysis. Here, we designed a hypothetical model in which the topographic factors and structural variables could directly influence the AGB. The topographic factors could also directly influence tree density and DBH diversity. These correlations are supported by previous studies [21,30]. We selected the best SEMs based on the fit indices, including a nonsignificant Chi-square (χ 2 ) test statistic (p > 0.05), S SRMR < 0.08, and both GFI and CHI > 0.95，through removing the nonsignificant variables with the lowest Akaike information criterion (AIC) [21,31].
All statistical analyses were performed in R 3.5.1. The predicted map of AGB was depicted using ordinary kriging in the geoR package [32]. The PCNM analysis was conducted using the PCNM package [33]. Variation partitioning analysis was computed using the 'varpart' function in the vegan package [34]. The SEM analysis was performed using the lavaan' package [35].

Results
According to the first census, the plot had 110,203 individuals belonging to 61 families, 147 genera, and 227 species. The average DBH of trees was 4.48 cm. The mean stem density was 4408 individuals per ha and the mean basal area was 17.00 m 2 /ha. The average forest AGB was 73.92 ± 23.51 Mg/ha (± standard deviation) in the 20 m × 20 m quadrats of the Mulun plot, although AGB ranged widely from 3.22 to 198.11 Mg/ha. AGB among different diameter ranges was submitted to a normal distribution (Shapiro-Wilk test; p = 0.0526), and the biomass of trees in the 5-10 cm DBH range was higher than in the other diameter classes (Figure 2). Trees from 1 to 5 cm DBH accounted for only 14.89% of the biomass, although they accounted for 72.22% of the stems. The forest AGB showed obvious spatial heterogeneity in the Mulun plot (Figure 3) where the topography, and especially the elevation, played a vital role in regulating the spatial pattern of AGB. The mountain peak almost had a higher biomass than the low-lying areas. Variance partitioning of AGB showed that 67.7% of the variation was explained by stand structural, topographic, and spatial factors together. Specifically, stand structural complexity explained 56.8% (b+ab+bc+abc) of the variation, topographic variables explained 8.0% (a+ab+ac+abc) of the variation, and spatial variables (70 PCNMs) explained 49.1% (c+ac+bc+abc) of the variation (Figure 4). Among The forest AGB showed obvious spatial heterogeneity in the Mulun plot (Figure 3) where the topography, and especially the elevation, played a vital role in regulating the spatial pattern of AGB. The mountain peak almost had a higher biomass than the low-lying areas. Variance partitioning of AGB showed that 67.7% of the variation was explained by stand structural, topographic, and spatial factors together. Specifically, stand structural complexity explained 56.8% (b+ab+bc+abc) of the variation, topographic variables explained 8.0% (a+ab+ac+abc) of the variation, and spatial variables (70 PCNMs) explained 49.1% (c+ac+bc+abc) of the variation (Figure 4). Among seven topographic variables, elevation and VDCN were significant (p < 0.01), but only explained 5.4% and 1.1% of the AGB spatial variance, respectively. Tree density and DBH diversity were significant (p < 0.001) structural variables that explained 18.2% and 38.3 of the AGB spatial variance, respectively (Table 1).   . Variance partitions of topographic, stand structural, and spatial variables account for the AGB variation (%). The symbols a, b, c represent the independent effects of topography, stand structure, and selected spatial variables, respectively; ab is the interactive effect of topography and stand structure; ac, the interactive effect of topography and selected spatial variables; bc, the interactive effect of stand structure and selected spatial variables; and abc, the interactive effect of topography, stand structure and selected spatial variables.    . Variance partitions of topographic, stand structural, and spatial variables account for the AGB variation (%). The symbols a, b, c represent the independent effects of topography, stand structure, and selected spatial variables, respectively; ab is the interactive effect of topography and stand structure; ac, the interactive effect of topography and selected spatial variables; bc, the interactive effect of stand structure and selected spatial variables; and abc, the interactive effect of topography, stand structure and selected spatial variables. . Variance partitions of topographic, stand structural, and spatial variables account for the AGB variation (%). The symbols a, b, c represent the independent effects of topography, stand structure, and selected spatial variables, respectively; ab is the interactive effect of topography and stand structure; ac, the interactive effect of topography and selected spatial variables; bc, the interactive effect of stand structure and selected spatial variables; and abc, the interactive effect of topography, stand structure and selected spatial variables. According to the results of PCNM, the topographical properties (elevation, VDCN, convexity, and slope) and structural properties (tree density and DBH diversity) served as two variables for carrying out SEM ( Figure 5). This SEM explained 74.6% of the AGB (p < 0.001) in the Mulun plot as follows: elevation (β = 0.15), DBH diversity (β = 0.68), and tree density (β = 0.71) positively influenced AGB. However, DBH diversity was directly affected by VDCN (β = 0. 19), and also negatively affected by elevation (β = −0.21), convexity (β = −0.16), and slope (β = −0.26). Tree density was directly affected by elevation (β = 0.36) and slope (β = 0.27), and negatively affected by elevation (β = −0.08). According to the results of PCNM, the topographical properties (elevation, VDCN, convexity, and slope) and structural properties (tree density and DBH diversity) served as two variables for carrying out SEM ( Figure 5). This SEM explained 74.6% of the AGB (p﹤0.001) in the Mulun plot as follows: elevation (β = 0.15), DBH diversity (β = 0.68), and tree density (β = 0.71) positively influenced AGB. However, DBH diversity was directly affected by VDCN (β = 0.19), and also negatively affected by elevation (β = −0.21), convexity (β = −0.16), and slope (β = −0.26). Tree density was directly affected by elevation (β = 0.36) and slope (β = 0.27), and negatively affected by elevation (β = −0.08).

Discussion
In this study, the average AGB (73.92 Mg/ha) of trees in a mixed evergreen-deciduous broadleaf forest was lower than the AGB (136.4 Mg/ha) estimated in the same type of karst forest in the neighboring Guizhou province [13]. The difference may be attributed to differences in the plot areas studied. An investigation using a large plot may allow for greater consideration of spatial variability of community structure. More robust estimates can be achieved with large sampling [15,36]. We also found that the karst forest had lower AGB than typical forests in non-karst regions in the same climate zone [15,36]. However, the highest value of AGB (198.11 Mg/ha) was 22.96 times as large as the lowest one. Therefore, this karst forest still has a certain level of potential for productivity due to the large proportion of small trees. Thus, the specialty of the nonzonal vegetation should be considered in estimating the regional carbon pools. The AGB has great variability with a variation coefficient of 30.80%. Other studies also showed a high variation in AGB

Discussion
In this study, the average AGB (73.92 Mg/ha) of trees in a mixed evergreen-deciduous broadleaf forest was lower than the AGB (136.4 Mg/ha) estimated in the same type of karst forest in the neighboring Guizhou province [13]. The difference may be attributed to differences in the plot areas studied. An investigation using a large plot may allow for greater consideration of spatial variability of community structure. More robust estimates can be achieved with large sampling [15,36]. We also found that the karst forest had lower AGB than typical forests in non-karst regions in the same climate zone [15,36]. However, the highest value of AGB (198.11 Mg/ha) was 22.96 times as large as the lowest one. Therefore, this karst forest still has a certain level of potential for productivity due to the large proportion of small trees. Thus, the specialty of the nonzonal vegetation should be considered in estimating the regional carbon pools. The AGB has great variability with a variation coefficient of 30.80%. Other studies also showed a high variation in AGB at local scales [15,36,37]. The difference in topography, forest structure, light condition, and their interaction are likely to contribute to this variation [37].
The variance partition results demonstrated that DBH diversity and tree density are the crucial stand structural driving factors for AGB, which is consistent with Xu et al. who reported that biomass of a subtropical mountain moist forest was significantly affected by stem density [15]. Ali et al. also reported that stand structural complexity (tree DBH diversity) may affect AGB [21]. We further found that tree density and DBH diversity had significantly positive effects on AGB in a karst forest based on SEM ( Figure 5). Many studies indicated that both stem density and tree size variation have stronger direct effects on AGB [17,21,24,38,39]. In the Mulun plot, AGB increased with increasing tree density or DBH diversity, and then decreased ( Figure 6). Hui et al. and McEwan et al. also found this unimodal relationship between density and biomass [17,39]. Within a certain range, AGB was also found to increase with DBH variation or diversity [21,24]. Niche differentiation and facilitation allow species complementarity due to variations in individual tree sizes [21], and may improve light capture and light-use efficiency through greater canopy densities, as well as enhancing AGB. Therefore, maintaining the appropriate levels of DBH diversity is an efficient way to improve the carbon sequestration potential of the forest ecosystem [21].
Forests 2020, 11, x; doi: FOR PEER REVIEW www.mdpi.com/journal/forests at local scales [15,36,37]. The difference in topography, forest structure, light condition, and their interaction are likely to contribute to this variation [37]. The variance partition results demonstrated that DBH diversity and tree density are the crucial stand structural driving factors for AGB, which is consistent with Xu et al. who reported that biomass of a subtropical mountain moist forest was significantly affected by stem density [15]. Ali et al. also reported that stand structural complexity (tree DBH diversity) may affect AGB [21]. We further found that tree density and DBH diversity had significantly positive effects on AGB in a karst forest based on SEM ( Figure 5). Many studies indicated that both stem density and tree size variation have stronger direct effects on AGB [17,21,24,38,39]. In the Mulun plot, AGB increased with increasing tree density or DBH diversity, and then decreased ( Figure 6). Hui et al. and McEwan et al. also found this unimodal relationship between density and biomass [17,39]. Within a certain range, AGB was also found to increase with DBH variation or diversity [21,24]. Niche differentiation and facilitation allow species complementarity due to variations in individual tree sizes [21], and may improve light capture and light-use efficiency through greater canopy densities, as well as enhancing AGB. Therefore, maintaining the appropriate levels of DBH diversity is an efficient way to improve the carbon sequestration potential of the forest ecosystem [21]. In the Mulun plot, of the seven topographic variables, elevation and VDCN explained the majority of variation, followed by convexity and slope. Generally, it is difficult to establish causal In the Mulun plot, of the seven topographic variables, elevation and VDCN explained the majority of variation, followed by convexity and slope. Generally, it is difficult to establish causal relationships between topography and biomass because topography is a composite variable, which can directly and indirectly affect the forest biomass. On the basis of the result of variance partitioning, we attempted to test the direct effects of topography on AGB and the indirect effects via structural attributes. The SEM results showed that elevation had a direct positive effect on AGB. In the Mulun plot, the relative elevation difference reached 209 m from depression to peak. The karst forest AGB increased with elevation ( Figure 6), which is consistent with Hui et al. and de Castilho et al. [17,40]. Additionally, McEwan et al. also found that elevation was significantly related to AGB in a Fushan 25 ha forest plot [39]. However, in some other large plots, AGB seemed to be insensitive to elevation, for example, in the Nouragues and Piste de Saint-Ellie 70 ha plot [41], Barro Colorado Island 50 ha plot [42], and Lienhuachih forest 25 ha plot [39]. This may be due to the differences in species diversity or climatic region. Further work is needed to investigate the causes.
The effect of topography on biomass is largely indirect, with little independent influence on biomass. In our study, topography also had an indirect effect on AGB via tree density and DBH diversity, consistent with previous studies [15,43]. Li et al. noted that environmental conditions strongly influence plant community compositions, and therefore, ecosystem functioning such as biomass [20]. Elevation and convexity had a positive influence on tree density in the Mulun plot. Consistent with our results, Xu et al. suggested that convexity was an important variable in regulating tree density [15]. Yuan et al. also found that altitude significantly correlated with the individual tree number in karst mountain forests in Southwest Guangxi, China [44].

Conclusions
We investigated a karst forest biomass and examined stand structural and abiotic factors driving AGB in a 25 ha dynamic forest plot in Southwest China, which provided important information on biomass distribution patterns for a better understanding of karst ecosystems. We found that the mean forest AGB was 73.92 ± 23.51 Mg/ha. The structural, topographic, and spatial factors predominantly contributed to the distribution of AGB. The stand structural complexity (tree density and DBH diversity) and topographic properties (elevation, VDCN, convexity, and slope) most strongly influenced AGB. The SEM suggested that AGB is directly driven by tree density, DBH diversity, and elevation. In addition, topographic factors mainly had indirect influence on AGB. Our work provides the links to understanding and establishing relationships between AGB and the driving factors on karst forests. These results will contribute to provide basic data for global and regional vegetation and carbon models. Further, regional forest carbon models need to consider community structural and topographic parameters.