Differences in spatial patterns and driving factors of biomass carbon density between natural coniferous and broad-leaved forests in mountainous terrain, eastern Loess Plateau of China

Mountain forests in China are an integral part of the country’s natural vegetation. Understanding the spatial variability and control mechanisms for biomass carbon density of mountain forests is necessary to make full use of the carbon sequestration potential for climate change mitigation. Based on the 9th national forest inventory data in Shanxi Province, which is mountainous terrain, eastern Loess Plateau of China, we characterized the spatial pattern of biomass carbon density for natural coniferous and broad-leaved forests using Local Getis-ord G* and proposed an integrative framework to evaluate the direct and indirect effects of stand, geographical and climatic factors on biomass carbon density for the types of forests using structural equation modeling.


Results
There was no signi cant difference between the mean biomass carbon densities of the natural coniferous and broad-leaved forests. The number of spots with a spatial autocorrelation accounted for 51.6% of all plots of the natural forests. Compared with the broad-leaved forests, the hot spots at the 1% signi cance level for the coniferous forests were distributed in areas with higher latitude, higher elevation, lower temperature, and lower precipitation. Geographical factors affected biomass carbon density positively and indirectly, via the stand and climatic factors, with larger effects for the natural coniferous than broad-leaved forests. Latitude and elevation are the most crucial driving factors for coniferous forests, but stand age and forest coverage are for broad-leaved forests. Climatic factors had weaker effects than other factors, with negative effects of temperature for coniferous and no effects for broadleaved forests.

Conclusions
The effects of stand, geographical and climatic factors on biomass carbon density are different between natural coniferous and broad-leaved forests, respectively. Employing the integrative framework can improve the prediction of the impact of stand, geographical and climatic factors on natural forests in mountainous areas.

Background
Human activities, such as the combustion of fossil fuels and deforestation, increase the concentration of atmospheric carbon dioxide that is one of the main causes of global climate change (IPCC, 2014). The role of forests as terrestrial sinks of atmospheric carbon dioxide has received increasing attention since the late 1990s (FAO, 2015). A series of forest management activities, such as LULUCF (land use, land-use change, and forestry) and REDD+ (reducing emissions from deforestation and forest degradation and other activities), have been implemented to mitigate the effects of climate change. Assessment of forest carbon storage and studying its driving factors of regional forests, including tropical (Chave et al., 2003;Slik et al., 2010), subtropical (L. ; M. , temperate and boreal forests (Wei et al., 2013), will improve our knowledge of terrestrial forest carbon to assist forest management practices.
Mountain forests cover about 900 million hectares of the world's land surface, constituting 24.3 percent of the world's forest cover (Price et al., 2000). Over 400 million hectares of the world's mountain forests are coniferous forests. The remainder is broad-leaved (Price et al., 2000). In China, forests cover about 195 million hectares (20.36%) of the national land surface, and the coniferous and broad-leaved forests accounted for 39.86% and 45.99% of the total forests area, respectively . Most of these forests are distributed in mountainous areas (Investigation and Planning Institute of the Ministry of forestry of China, 1981). Obviously, the mountain forests are a very important part of the forests whether in the world or China. The mountain forests are considered to be terrestrial sinks of carbon dioxide (Chang et al., 2015;Zald et al., 2016). However, because of their steep slopes, frequent extreme weather events, the ecosystem of mountainous areas is very fragile. To more effectively manage mountain forests and avoid the destruction of mountain forests, spatial distribution of forest carbon density and its in uencing factors is needed to be explored including terrain, site and environmental parameters, etc (FAO, 2015). However, there is little information on the spatial patterns of forest carbon densities and on how the forest carbon densities are driven by various factors for natural coniferous and broad-leaved forests in temperate mountainous areas.
The spatial patterns of forest carbon densities could not only play an important role in the evaluation of forest carbon sequestration potential and forest management practices (Fu et al., 2015;Lin et al., 2017) but also re ect the relationship between the geographic environment and forest carbon storage to some extent. Currently, many studies on the spatial pattern of carbon stocks (density) in large regions were reported in China. For instance, Ren et al. (Ren et al., 2013) examined the spatial pattern of carbon density in forest ecosystems in Guangdong Province and found that the spatial distribution of forest carbon density was uneven, and it was a re ection of the difference in forest management and economic and social development. Dai et al. (Dai et al., 2018) found that the forest carbon density decreased from southwest to northeast in Zhejiang Province, roughly in line with the topographic feature across this province by using Anselin Local Moran's I and geostatistical interpolation. Lin et al. (Lin et al., 2017) analyzed the spatial variability of forest carbon density in Jiangle County, Fujian Province by using Anselin Local Moran's I, Local Getis-Ord G* and geostatistical interpolation. However, little attempt has been made to compare the spatial patterns of carbon densities between the natural coniferous and broad-leaved forests.
The forest carbon density and its spatial pattern may be the result of complex interactions among stand, geographical and climatic factors. In mountain areas, the effect of elevation on forest carbon density has attracted great attention. China. This pattern would be attributed to the changes in temperature, precipitation, plant community type, and stem density with the elevation (Fisk et al., 1998;Sanaei et al., 2018). Similarly, latitude is also a critical geographical factor in determining forest carbon density. Wen et al. (Wen and He, 2016) found that the forest carbon density decreased with increasing latitude along the north-south transect of eastern China, while the biomass carbon density of Moso bamboo forests linearly increased with latitude (L. Xu et al., 2018; M. . In addition, Fan et al. (FAN et al., n.d.) also reported that carbon storage of Phyllostachys edulis forest was signi cantly affected by slope aspect and slope position. Meanwhile, the variation in forest carbon density had been related to changes in climatic condition, mainly including temperature and precipitation, which could affect the plant productivity, soil moisture availability, species distribution, and stand structure. For example, the annual mean temperature had negative indirect effects on carbon storage in dry, moist, and wet old-growth tropical forests across the globe based on multigroup structural equation modeling (Durán et al., 2015), and annual mean precipitation had a positive relationship with aboveground carbon density in temperate mature forests (Stegen et al., 2011). Moreover, stand factors were strong predictors of forest carbon density. For instance, Xu et al. (Xu et al., 2010) reported that forest age had a signi cant effect on biomass carbon density for most forest types in China. The average above-ground biomass carbon density of mature forests increased with stand age for the forests younger than 450 years . Up to now, our knowledge on the concomitant effects of various factors above mentioned on the carbon density of natural forests in temperate mountain areas had been far from complete.
Structural equation modeling (SEM) is a powerful tool for unraveling the structure linking variables that are correlated in a multivariate way, allowing the disentanglement of direct and indirect effects of a predictor variable (Shipley, 2016). It has been increasingly being used in researches in natural systems and ecology and contributes to theory maturation (Ali et  Shanxi Province is a mountainous terrain in north China, characterized by six mountain ranges, namely the Hengshan and Wutai Mountains in the north, the Zhongtiao Mountains in the south, the Lüliang Mountains in the west, the Taihang Mountains in the east, and the Taiyue Mountains at the center ( Fig. 1a). The forests which mainly consist of coniferous and deciduous broad-leaved forests are widely distributed in mountainous areas in the province. Currently, few studies have examined the carbon density (stock) of the forests in some mountainous areas of Shanxi Province. For example, Liu and Nan (Liu and Nan, 2018) examined the carbon stocks of three natural coniferous forests (Larix principisrupprechtii forest, Picea meyerii forest, and Pinus tabulaeformis forest) along an altitudinal gradient from 1200 to 2700m in the Guandi Mountain. Wang et al. (Wang et al., 2018) analyzed the spatial patterns of carbon densities of natural and planted forests in the Lüliang Mountains. Besides, Wang et al. (Wang, 2014) studied the carbon densities (stocks) of seven forest types in Shanxi, including Populus davidiana forest, Pinus tabuliformis forest, Platycladus orientalis forest, Quercus variabilis forest, Larix principisrupprechtii forest, Betula platyphylla forest, and Robinia pseudoacacia forest, and evaluated the carbon density basing on the allometric growth model of single-tree biomass built by themselves and analyzed the distribution patterns of the carbon density in the various layers (including arbor layer, shrub layer, grass layer, litter layer and soil layer) for each forest ecosystem. However, there has been no complete report on the carbon density for all the natural forests across the entire province. Almost all of natural forests in Shanxi are distributed in the mountainous areas (Fig. 1b), so it provides a good chance to study the spatial patterns and driving factors of the forest carbon density in the mountainous areas. What is more, there has been very little information on the comparative study between natural coniferous and broad-leaved forests on the spatial patterns of carbon densities and on how the carbon densities are driven by various factors.
In view of this, we rstly estimated the biomass carbon densities of natural coniferous and deciduous broad-leaved forests in Shanxi Province, China, using the method of biomass expansion factor based on the national forest inventory data in 2015. We then characterized the spatial variation of biomass carbon density of natural forests using Local Getis-Ord G*. Finally, we established an integrated model which revealed all possible connections between all variables including stand factors, geographical factors, climatic factors, and biomass carbon density. We further used multi-group structural equation modeling to evaluate the relative importance of driving factors in explaining variations of biomass carbon density, to estimate the direct and indirect effects of each driving factor, and to reveal the similarities and differences between the models for natural coniferous forests and natural deciduous broad-leaved forests.

Study region
Shanxi Province, also known as "Shanxi Plateau", is located in the east of the Loess Plateau in north China, and in the middle reach of the Yellow River, at 34°34′N-40°43′N and 110°14′E-114°33′E, with a total area of 156,271 km 2 (Fig 1a). The area of mountains and hills accounts for 80.1% of the total area of the province. The highest peak is Wutai Mountain in northeastern Shanxi with an altitude of 3,058 m (Ma, 2001), which is the highest peak in north China. Climatically, it belongs to the temperate continental monsoon climate and is in uenced by the summer and winter monsoons annually, with a well-de ned rainy season and a very dry winter (Fan and Wang, 2011). In terms of vegetation division, the south of Hengshan Mountain belongs to the warm temperate deciduous broad-leaved forest zone (Fig 1a) (Ma, 2001). The total area of forests is about 32,100 km 2 , accounting for 20.05% of the province's land area. Natural forests account for nearly 60% of all the forests in the province.

Data collection and preprocessing
Forest inventory data The data from a total of 878 eld sample plots were used in this study. Of all the plots from the 9th national forest inventory data in 2015 in Shanxi Province, 345 and 533 are natural coniferous forests and natural deciduous broad-leaved forests plots, respectively (Tables 1-2, Fig 1b). The permanent plots (each with an area of 667 m 2 ) were established systematically based on a 4km×4km grid (Xiao, 2005). The data collected in each plot included forest origin, dominant tree species, forest stand factors (stand age, age group, diameter at breath height of 1.3m (DBH), tree height and forest coverage), geographic location (latitude, longitude), and topographic factors (elevation, slope aspect, and slope position). For trees with DBH≥5cm, their DBHs were recorded. According to the dominant species, we classi ed all the natural coniferous forests into six categories, namely, Pinus tabuliformis forests, Pinus bungeana forests, Picea wilsonii forests, Larix principis-rupprechtii forests, Platycladus orientalis forests, and coniferous mixed forests. Correspondingly, we grouped the natural deciduous broad-leaved forests mainly into the following classes: Quercus liaotungensis forests, Quercus variabilis forests, Quercus acutissima forests, Betula platyphylla forests, Populus davidiana forests, deciduous broad-leaved mixed forests, and other forest types ( Fig S1). The shrub and herb characteristics were measured from three randomly selected 2m*2m quadrats in each plot.
The biomass carbon density (BCD) (Mg/ha) of individual living trees for each forest plot was estimated using the biomass expansion factor (BEF) method (Sun et al., 2020;Wang et al., 2018). The BCD of the shrub, herb, or litter layer for the forest-present plot was obtained by multiplying the biomass (Mg/m) by carbon content, which was 0.4627, 0.3270, and 0.4700 for the shrub, herb, and litter layers, respectively (Wang et al., 2018).
The biomass carbon density in this study represented the total carbon density in living trees, shrubs, herbs, and litter. We computed the mean BCD of each forest type (each forest age group) by dividing the sum of biomass carbon densities of the forest plots by the number of forest plots for the forest type (the forest age group).

Climate data
There are 108 surface meteorological stations in Shanxi Province. We collected the dataset of annual mean temperature (TEMP) and annual mean precipitation (PRCP) for the period 1981-2015 at all these stations from the meteorological database of the scienti c data platform of Shanxi Province. We then used the regression-kriging interpolation method (Lamsal et al., 2011;Plouffe et al., 2015;Sun et al., 2020) to derive the climatic data for each plot from that for each meteorological station, which combined regression of the climatic variable on topography variables (slope degree, slope aspect, and slope position) with kriging of the regression residuals (Eldeiry and Garcia, n.d.). We nally further corrected the derived TEMP for each sample plot according to the differences between the plot elevation and the interpolated elevation using the temperature lapse rate of 4.89°C km −1 as the correction factor (Sun et al., 2020; Wang, 2014).

Data analyses Statistical analyses
On the basis of passing the tests of normal distribution and homogeneity of variance, we used One-way ANOVA with Duncan's multiple tests to examine the difference in mean biomass carbon density among forest types or forest age groups.

Spatial hot and cold spots analyses
We used the local Getis-Ord G* to characterize the spatial hot and cold spots for forest BCDs, which works by looking at each feature within the context of neighboring features. When similar values are clustered and signi cantly larger (lower) than the mean, it is identi ed hot (cold) spots area. In the process of the local Getis-Ord G* test, the G* statistic is a Z score. For statistically signi cant positive Z scores, the larger the Z score is, the more intense the clustering of high values (hot spot). For statistically signi cant negative Z scores, the smaller the Z score is, the more intense the clustering of low values (cold spot). The "not signi cant" positive or negative Z scores represent no spatial autocorrelation (Manepalli et al., 2011). Getis-Ord Gi* can be described as: where x j is the value of attribute for observations j, and is the mean of the corresponding values of the attribute. The w i,j is the spatial weight between observations i and j, which can be de ned as the xed distance, which is given the same weight within the distance, while those outside the distance band are given the weight of 0. The symbol of n is the total number of observations.
We used the module of Hot Spot Analysis (Getis-Ord Gi*) of Spatial Statistics Toolbox within software ARCGIS 10.2 to compute the Local Getis-Ord Gi* index.

Multi-group SEM
We used SEM to assess the relative contribution of various in uencing factors to BCD. Multi-group SEM was used to compare the direct and indirect effects of these factors on BCD between natural broadleaved forests and natural coniferous forests. An exploratory SEM was supposed (Fig S2). We assumed that stand factors (stand age, stand coverage, and forest type) had direct effects on BCD. We also hypothesized that climatic (temperature and precipitation) factors had direct effects on BCD and indirect effects through direct effects on the stand factors. Moreover, geographic factors (latitude and elevation) affected BCD indirectly via the effects on stand and climatic variables. After testing the collinearity of all factors, we began with the construction of a prior model that included all hypothesized causal relationships between these variables based on the known effects of driving factors on BCD and relationships among these factors. As shown in Fig S3, seven driving factors were included in the prior model: stand age (AGE), stand coverage (COV), forest type (TYPE), latitude (LAT), elevation (ELE), mean annual precipitation (PRCP), and mean annual temperature (TEMP). Then we compared the actual and estimated covariance matrices of the SEM to determine whether the model was acceptable by using the indices of chi-square (χ2) test, comparative t index (CFI) and root square mean error of approximation (RMSEA). If the model is not acceptable, the model needed to be revised by removing or adding a connection between variables according to the modi cation indices and corresponding parameter changes expected. The process is performed continuously until nding a good-tting model ( Fig S3).
Upon nding a good-tting model, multi-group SEM was used to test whether the paths of the model differ statistically between natural coniferous and broad-leaved forests. In detail, we rstly started with a full free multi-group model with no between-group equality constraints to test that the same qualitative structure could apply to the two groups (P > 0.05) (Shipley, 2016), and then forced all the path coe cients of multi-group SEM to be equal between natural broad-leaved forests and coniferous forests, and judged whether the most constrained multi-group SEM were acceptable. If this fully constrained model was rejected (P<0.05), indicating that there was at least one equality constraint between the two groups was unreasonable (Shipley, 2016). We secondly tted a series of nested models to detect which paths were equal or which paths differ across groups by constraining only one path to be equal between groups at a time ( Table 3). The difference in the maximum likelihood χ2 statistics between models before and after adding an equal constrained path was used to test if this change was acceptable. Finally, we tted the multi-group SEM in which those acceptable path coe cients were forcing to be equal between the two forests type. The SEM and multi-group SEM were implemented using the R version 3.4.4 and lavaan package version 0.6-3 (Rosseel, 2012).

Biomass carbon density
The BCDs of all natural mountain forests in Shanxi varied from 4.79 to 209.16 Mg/ha, with a median value of 31.53 Mg/ha. There was no signi cant difference between the mean BCDs of coniferous and broad-leaved forests (P > 0.05), but the coe cient of variation of BCDs was obviously larger for coniferous forests (72.93%) than broad-leaved forests (57.35%) (Fig. 2). For all the natural forests, the average BCD of live trees was 30.06 Mg/ha, accounting for 82.4% of the total biomass carbon density ( Table 2). The mean BCD of live trees also showed no signi cant difference between the broad-leaved and coniferous forests (P > 0.05). However, the mean BCD of undergrowth including shrubs and herbs was larger for the coniferous than broad-leaved forests.
Pinus tabulaeformis forests and Quercus liaotungensis forests were the dominant forest types for the coniferous and broad-leaved forests, respectively, accounting for 51.01% and 38.27% of the corresponding total area ( Table 2). The mean BCDs of the arbor layer, undergrowth layer, and the total of them for Pinus tabulaeformis was signi cantly lower than those for Quercus wutaishansea Mary forests (P < 0.01).  (Over-matured). Across all the ve age groups but the middle-aged group, BCD of coniferous forests was lower than that of broadleaved forests (P < 0.05) (Fig. 3). The BCDs increased with forest age for both the coniferous and broadleaved forests, which were maximum near 85 year and 95 year, respectively (Fig. 3).

Spatial hot and cold spots analyses for biomass carbon densities
Anselin Local Moran's I analysis showed that 283 hot-spots and 170 cold-spots were detected by Getisord G* for BCDs of natural mountain forests in our study area (Fig. 4). The total number of spots with spatial autocorrelation accounted for 51.6% of all plots of the natural forests. The hot spots of biomass carbon density were distributed in the north of 36. Hot spots at the 1% signi cance level (HS99s), which have the highest average biomass carbon density of 47.86 Mg/ha, were distributed in the areas with an elevation of 1500-2000m, the mean annual temperature of 8 ~ 10.8 o C, and the mean annual precipitation of 450-490mm (Fig S4). Compared with the deciduous broad-leaved forests, HS99s for the coniferous forests were distributed in the areas with higher average latitude, higher average elevation, lower mean temperature, and lower mean precipitation ( Fig   S4). Overall, HS99s were mainly located in Luya Mountain, Guandi Mountain, and Taiyue Mountain from north to south (Fig. 4). In Luya Mountain, the number of HS99s for coniferous forests (18) dominated by Picea wilsonii forests and Larix principis-rupprechtii forests was more than that for deciduous broadleaved forests (11) dominated by Betula platyphylla forests and Quercus liaotungensis forests. In Guandi Mountain, the number of HS99s for coniferous forests (42) dominated by Pinus tabuliformis forests and coniferous mixed forests was more than that for deciduous broad-leaved forests (27) dominated by Quercus liaotungensis forests. In Taiyue Mountain, the number of HS99s for deciduous broad-leaved forests (39) mainly distributed in the west and dominated by Quercus liaotungensis forests and broadleaved mixed forests was more than that for coniferous forests (27) mainly distributed in the east and dominated by Pinus tabuliformis forests. Generally, the number of HS99s for coniferous forests was more than that for deciduous broad-leaved forests in the northern mountainous areas, but it was contrary in the southern mountainous areas of Shanxi province.
In uences of various factors on biomass carbon density of natural mountain forests The multi-group SEM analysis revealed variation in path coe cients between coniferous and deciduous broad-leaved forests (Fig. 5). Except for the ve paths including the effects of AGE and COV on BCD, the impact of PRCP on BCD, and the in uences of LAT on PRCP and AGE, the remaining paths were not equal between coniferous and deciduous broad-leaved forests. The seven factors included in our model could better explain the variation in BCD for coniferous forests (62.0%) than deciduous broad-leaved forests (51.6%) (Fig. 5). AGE and COV showed a direct and positive association with BCD consistently for coniferous and deciduous broad-leaved forests (Fig. 6). The effect of TYPE on BCD was stronger for the coniferous forests than the deciduous broad-leaved forests. LAT (or ELE) had only indirect effects on BCD, with the total positive effect of LAT (or ELE) stronger for the coniferous forests than the deciduous broad-leaved forests. TEMP had both direct and indirect effects on BCD for coniferous forests but no signi cant effects for deciduous broad-leaved forests. PRCP had weak negative effects on BCD for the coniferous forests and the deciduous broad-leaved forests. In short, the results of standardized total effects (Fig. 6) indicated that, for the coniferous forests, geographical factors (LAT and ELE) and biotic factors (AGE, COV, and TYPE) were important driving factors of BCD, and the former had stronger effects than the latter. While for the deciduous broad-leaved forests, biotic factors (AGE and COV) were the most crucial in uencing factors. No matter for the coniferous forests or the deciduous broad-leaved forests, climatic factors, speci cally temperature and precipitation, were fewer correlates of BCD than other factors.

Spatial patterns of biomass carbon density
The results of spatial hot and cold spots analyses showed that the areas in which the plots with high biomass carbon density distributed continuously in space were identi ed as hot spots areas (Fig 4). These hot spots were widely distributed in the mountainous areas with higher elevation (Fig S4), in most of which natural reserves were established, for instance, Luyashan National Natural Reserve in Luya Mountain, Pangquangou National Natural Reserve in Guandi Mountain, Lingkongshan National Natural Reserve in Taiyue Mountain. In these reserves, the forests were less disturbed by human activities and would be relatively older. Meanwhile, the coniferous forests in hot spots areas consisted of Picea wilsonii forests, Larix principis-rupprechtii forests, coniferous mixed forests, and Pinus tabulaeformis forests, and deciduous broad-leaved forests consisted of Quercus liaotungensis forests, Betula platyphylla forests, and deciduous broad-leaved mixed forests. Compared with the total average biomass carbon density of natural forests, these forests had higher mean biomass carbon densities ( Table 2). This promoted the formation of hot spots in these areas. Meanwhile, continuous distribution of the natural coniferous (deciduous broad-leaved) forests also helped to form hot spots areas in these areas (Fig 1b). In addition, the composition of tree species in hot spots areas was also relatively simple, and every tree species was widely distributed and concentrated in groups (Fig S1). This also reduced the inhomogeneity of spatial distribution of biomass carbon density in space due to the differences among forest types in these areas (Figs 5 and 6), which was unfavorable for the information of the hot spots.
In contrast to the hot spots, the cold spots for biomass carbon density of natural mountain forests were mainly detected in or around the Coal Mining areas. Shanxi province being abundant in coal resources is one of the best important providers of coal in China ( Shanxi Geological Exploration Bureau, 2018.). For instance, cold spots at the 1% signi cance level (CS99s) in the southernmost tip of the Lüliang Mountains are located in Xiangning colliery, which is part of Hedong coal eld (Shi, 2017; Shanxi Fenwei Energy Development Consulting Co., Ltd, 2007). Another CS99s region detected in the south of Taihang Mountains lied to the east of Lingchuan colliery or in the Jincheng colliery, which belongs to Qinshui coal eld (Shi, 2017; Shanxi Fenwei Energy Development Consulting Co., Ltd, 2007). Higher impacts of long-term mining activities would lead to low forest biomass carbon densities and the formation of cold spots in these areas.
The Zhongtiao Mountain, known as "Shanxi natural botanical garden" (Zhang and Wu, 2015), is a key region for the woody ora of North China, in which many rare and endangered plants were distributed. In our study, the entire Mountain was identi ed as an area of random distribution for the biomass carbon densities of natural mountain forests (Fig 4). The most likely explanation for this is that rich and varied plant species are randomly distributed in space (Fig 1b). On the other hand, it's well known that the Mountain is one of the birthplaces of Chinese civilization, and the emperors of Yao, Shun, Yu, and Tang had been here many times (Zhang and Wu, 2015). With the continuous expansion of agricultural production and the increasing population density around the Mountain, human activities had a more and more profound impact on natural forest vegetation in the course of history. Meanwhile, in ancient times the area was a major center for copper mining (Chenw. et al., 1996) and is today one of the operational bases of Zhongtiao Mountain Non-ferrous Metals Group Co., Ltd, one of China's largest metal processing companies. Therefore, long-term mining and processing of mineral resources had also been producing a negative effect on natural forests in the mountains. Besides, the continuous development and utilization of Tourism Resources in the Mountain could also be a reason for the random distribution of biomass carbon densities in Zhongtiao Mountain. In a word, frequent human disturbance should also be a decisive factor which cannot be neglected.

The relationships between stand factors and biomass carbon density
The results of SEM analysis showed that AGE and COV were important predictors in modulating the BCD for mountainous coniferous and deciduous broad-leaved forests. Especially for the deciduous broadleaved forests, the effect of AGE or COV was stronger than that of the other factors (Fig 6). The ndings were in agreement with Xu et al. (L. Xu et al., 2018) who found that canopy density and forest age were the dominant determinants of vegetation carbon density in subtropical forest ecosystems in Zhejiang Province, China. The strong positive effect of AGE on BCD, possibly because the stand age controls the duration of forest carbon accumulation (Pregitzer and Euskirchen, 2004) and biomass carbon density increases with stand development (Cheng et al., 2014;Fonseca et al., 2011). While owing to nutrient limitation, stomatal constraint and decline in photosynthesis during the stand development, stand net primary productivity (NPP) declines along with the increase of tree age (Gower, 2003;Mcdowell and Harireche, 2002;Tang et al., 2014). Consequently, with the increasing stand age, the biomass carbon accumulation reaches a steady state when the carbon sequestration rate of forests approaches zero (Goulden et al., 2011). We found that equilibrium points were about 95yr and 85yr for mountainous broad-leaved and coniferous forests in our study area, respectively (Fig 3b). Meanwhile, TYPE was another stand factor in modulating the BCD for the mountainous forests. The effect of TYPE on BCD was greater for the coniferous forests than the deciduous broad-leaved forests (Fig 6). The result could be due to the differences of averaged BCDs of different forest types between the coniferous and broad-leaved forests. The further result showed that the difference of mean BCDs of different forest types for the coniferous forests was weaker signi cantly larger than that for the broad-leaved forests ( Table 2) (t=0.973, P=0.058).

The effects of geographical factors on biomass carbon density
Compared with other geographical factors (longitude, slope aspect, and slope position), the introduction of ELE and LAT into our model could be mainly determined by the special geographical environment conditions of mountainous forests in our study. Deciduous broad-leaved forests and coniferous forests ranged from 520m to 2309m and from 483m to 2560m in elevation, respectively. Elevation gradient is associated with changes in temperature and precipitation, and forest type (Sanaei et al., 2018). By regulating moisture and soil water availability (Fisk et al., 1998), elevation can affect forest canopy, stem density, and stand basal area, and therefore affect aboveground biomass (L. Xu (Fig 6). Similar results were also reported by Liu et al. (Liu and Nan, 2018), who found a positive linear relationship between vegetation carbon stock and altitude across three forests (Larix principis-rupprechtii (LP) forest, Picea meyerii (PM) forest, and Pinus tabulaeformis (PT) forest) on Loess Plateau in an altitudinal range of 1200-2700m. In addition, a similar relationship between biomass carbon density and elevation was reported by Sumeet Gairola et al. (Gairola et al., 2011) in moist temperate forests of the Garhwal Himalaya. As for LAT, it also had a signi cant effect on BCD in our study due to that the distribution of natural forests also had a larger latitude span. Deciduous broad-leaved forests and coniferous forests ranged from 34.79° to 39.89° and from 34.97° to 39.75° in elevation, respectively. The larger latitude span would lead to great changes in the environmental conditions along latitude gradients from the north to the south, including the light, heat, and moisture, and would affect plant growth and therefore in uence the biomass carbon accumulation in the natural forest from the north to the south. It is noteworthy that in the areas of natural forests distribution in our study, the higher the latitude, the higher the elevation (r=0.708, P<0.001). So the positive effects of ELE on BCD enhanced the positive effects of LAT (Fig 6).

The effects of climatic factors on biomass carbon density
The total effects of climatic factors (TEMP and PRCP) were lower than those of geographical factors (ELE and LAT) and biotic factors (AGE, COV, and TYPE) in temperate mountainous forests in Shanxi. According to the SEM, we found a negative effect of TEMP on BCD for the natural coniferous forests, indicating that the lower the temperature, the higher the BCD of the forest. It could be the result of natural selection of the coniferous tree species for the site conditions, and the spatial distribution of these tree species in our study. For example, the Picea and Larix principis-rupprechtii forests with the higher mean biomass carbon densities were generally distributed in the regions with lower temperatures, but Platycladus orientalis forests with the lowest mean biomass carbon density were generally distributed in warmer areas (Table 1). A similar result was also found by Fehse et al. (Fehse et al., 2002), who demonstrated that the forests under favorable site conditions at high altitudes with low temperature were not inferior in biomass accumulation and productivity compared with the forests at low altitude with high temperature. Meanwhile, the smaller leaf surface area of coniferous trees which makes the water and heat not easy to lose is bene cial to cold resistance. And the lower temperature for mountainous coniferous forests in our study not only could not be lower to limit the photosynthetic carbon sequestration rate of these coniferous trees but also might decrease the amount of carbon released by respiration because the respiration generally required a higher temperature than photosynthesis.
Consequently, lower temperatures for coniferous forests probably increased the biomass carbon accumulation. Besides, for the temperate deciduous broad-leaved forests, which were the zonal vegetation in our study area, TEMP had no signi cant effect on BCD. Similarly, Liu et al.  also found a nonlinear relationship exists between the above-ground biomass carbon density and mean annual temperature in the temperate mountain system for mature forests on the global scale.
Integrative framework of biomass carbon density for the mountainous forests factors on forest carbon storage. In our study, an integrative framework of biomass carbon density for the mountainous forests was proposed, which combined the three kinds of driving factors, including biotic (AGE, COV, and TYPE), climatic (TEMP and PRCP), and geographical factors (ELE and LAT), and then the direct and indirect effects of these factors on biomass carbon density were derived from it. For example, we found that geographical factors had affected biomass carbon density indirectly in our study (Fig S2).
In detail, ELE had affected BCD through the direct effects of ELE on TYPE, AGE, and TEMP ( Fig S2).
Meanwhile, we also found a positive relationship between ELE and AGE (Fig 5). This was probably because the forests at higher elevations had usually experienced little disturbance, and could grow steadily and presented larger average stand age until the mature stage. Moreover, ELE had effects on the distribution of the coniferous forest types for the coniferous forests. Noticeably, for coniferous forests, the Picea and Larix principis-rupprechtii forests were generally distributed at a higher elevation than Pinus tabuliformis and Platycladus orientalis forests across the entire region (Table 1), and the BCDs for the former two was generally greater than for the latter two (Table 2). In short, the indirect effects were helpful for us to understand why elevation was an important driving factor of biomass carbon density for mountain forests and how geographical factors affect biomass carbon density in our study.
Meanwhile, multi-group SEM as an integrative approach also was used to compare differences between groups in our study. It could be revealed whether the effects of these factors were consistent between the coniferous forests and broad-leaved forests. For instance, the effect of AGE (COV) on BCD for the coniferous forests was equal to that for deciduous broad-leaved forests, but the effects of climatic and geographical factors varied greatly between them. For the deciduous broad-leaved forests, the positive connection between ELE and AGE was the main reason for the positive effects of ELE on BCD. But for the coniferous forests, ELE had effects on TEMP and TYPE besides the positive effect on AGE (Fig 5).
Therefore, the total effect of ELE on BCD was obviously different between the deciduous broad-leaved and coniferous forests (Fig 6). Furthermore, if we did not use the multi-group SEM, the model we proposed in our study, the qualitative structure of which was correct, would be incorrectly rejected because the coniferous forests and deciduous broad-leaved forests differ in the numerical strength of these causal relationships.

Conclusions
In our study, we found that there was no signi cant difference between the mean biomass carbon density for the natural coniferous forests and that for the natural deciduous broad-leaved forests in Shanxi, while the coe cient of variation of biomass carbon densities was obviously larger in the coniferous forest than in deciduous broad-leaved forests. The spatial patterns of biomass carbon densities for the natural forests were associated with the geographical conditions, forest type distribution, human interference density, and the implementation of natural forest resources protection measures in Shanxi Province.
Compared with the deciduous broad-leaved forests, hot spots at the 1% signi cance level (HS99s) for the coniferous forests, which have the highest average biomass carbon density of 47.86 Mg/ha, were distributed in the areas with higher average latitude, higher average elevation, lower mean temperature, and lower mean precipitation.
Multi-groups SEM that combined stand, geographical, and climatic factors and was used to compare the effects of these factors on biomass carbon densities between the natural coniferous forests and the natural deciduous broad-leaved forests in mountainous areas. These factors, including latitude (LAT), elevation (ELE), mean annual temperature (TEMP), mean annual precipitation (PRCP), stand age (AGE), coverage (COV), and forest type (TYPE) could explain 62.0% of the variation in biomass carbon density for coniferous forests and 51.6% for deciduous broad-leaved forests. The geographical factors (LAT and ELE) only had an indirect effect on BCD, and the effect for the coniferous forests was obviously larger than that for the deciduous broad-leaved forests. The indirect effects of geographical factors also revealed how they affected the biomass carbon density via other factors. Meanwhile, the effects of stand age and coverage showed no differences between the coniferous forests and the deciduous broad-leaved forests, but the effect of forest type varied between them. The effect of forest type was greater for the coniferous forests than that for the deciduous broad-leaved forests, probably because biomass carbon densities varied more greatly among different forest types for the former than the latter. Climatic (TEMP and PRCP) factors had weaker effects on biomass carbon density than other factors in mountainous areas on a provincial scale. This study enhances the understanding of variation in the biomass carbon density in natural coniferous forests and deciduous broad-leaved forests in the mountainous region and should be bene cial for more effective implementation of afforestation and forest protection programs to enhance biomass carbon sequestration and mitigate global climate change.

Declarations Acknowledgments
We acknowledge the Carbon Sequestration Monitoring Center of Forestry Investigation and Planning Institute in Shanxi for providing us national forest inventory data as well as useful advices for their manipulation.

Authors' contributions
Conceived and designed the study: SL and FX. Analyzed the data: SL, WQ and FX. Wrote the paper: SL, WQ and FX. All authors read and approved the nal manuscript.

Availability of data and materials
The dataset used and/or analyzed during the current study are available from the corresponding author on reasonable request.  The rst row shows the maximum likelihood χ2 estimates (MLχ2) when all parameters are free. The remaining rows show the result by constraining one parameter at a time. The difference between the fully free model and the rest is given as ΔMLχ2, and the P-value indicates the probability that the constraint of that parameter changes the model. Bonferroni-corrected P-value threshold, 0.05/15 = 0.0033. BCD, biomass carbon density; AGE, stand age; COV, stand coverage; TYPE, forest type; LAT, latitude; ELE, elevation; PRCP, mean annual precipitation  The kernel density estimation of natural broad-leaved and coniferous forests biomass carbon density in Shanxi, China.

Figure 3
Average biomass carbon density of different age groups (a) and different age intervals (b) of natural broad-leaved and coniferous forests in Shanxi, China.

Figure 4
Map showing the spatial cold and hot spots of biomass carbon density in Shanxi, China. The number in parentheses are the number of sample plots.