Spatial distribution of soil organic carbon stock in Moso bamboo forests in subtropical China

Moso bamboo (Phyllostachys heterocycla (Carr.) Mitford cv. Pubescens) is an important timber substitute in China. Site specific stand management requires an accurate estimate of soil organic carbon (SOC) stock for maintaining stand productivity and understanding global carbon cycling. This study compared ordinary kriging (OK) and inverse distance weighting (IDW) approaches to study the spatial distribution of SOC stock within 0–60 cm using 111 soil samples in Moso bamboo forests in subtropical China. Similar spatial patterns but different spatial distribution ranges of SOC stock from OK and IDW highlighted the necessity to apply different approaches to obtain accurate and consistent results of SOC stock distribution. Different spatial patterns of SOC stock suggested the use of different fertilization treatments in Moso bamboo forests across the study area. SOC pool within 0–60 cm was 6.46 and 6.22 Tg for OK and IDW; results which were lower than that of conventional approach (CA, 7.41 Tg). CA is not recommended unless coordinates of the sampling locations are missing and the spatial patterns of SOC stock are not required. OK is recommended for the uneven distribution of sampling locations. Our results can improve methodology selection for investigating spatial distribution of SOC stock in Moso bamboo forests.


Results
Descriptive statistics. The summary of descriptive statistics for measured SOC stock is presented in Table 1.
Results show that SOC stock showed a decreasing trend with increasing soil depth; for 0-20 cm, SOC stock fell within 20.7-90.0 Mg ha −1 , and it was 10.5-72.2 Mg ha −1 for 20-40 cm, 6.6-80.1 Mg ha −1 for 40-60 cm and 52.8-229.7 Mg ha −1 for 0-60 cm. Mean SOC stock was 50.9, 42.6, 33.3 and 126.7 Mg ha −1 for soil layers 0-20 cm, 20-40, 40-60 cm and 0-60 cm, respectively. The lowest coefficient of variation (CV, 32%) was identified for SOC stock at the 0-60 cm layer while the highest CV (47%) was at the 40-60 cm layer. A CV value of 10% indicates low variability and values ranging from 10-90% indicate a moderate variability; CV values > 90% indicate high variability 31 . Therefore, SOC stock in our study area suggested a moderate variability. p values for the Shapiro-Wilk test ranged from 0.007 to 0.041, indicating a non-normal distribution of SOC stock for the different soil layers at the 0.05 level of significance. Therefore, before conducting spatial interpolation, a natural log-transformation was undertaken to meet the assumption of normal distribution.
The relationships between SOC stock and topographic variables. The relationships between topographic variables (elevation, slope and aspect) and SOC stock at the different soil layers were analysed using linear regression. Positive significant relationships between elevation and SOC stock (Table 2) were identified, indicating SOC stock increased with elevation. However, no significant relationship was found between slope, aspect and their interactions, and SOC stock. The interactions between elevation and slope and elevation and aspect led to significant effects on SOC stock, except for SOC stock at 20-40 cm and 40-60 cm layers for the interactions between elevation and slope. Spatial autocorrelation and trend surface analysis. Moran concurred with the assumption of OK, that the SOC stocks (regionalized variables) distributed in different locations were spatially correlated. Trend surface analysis indicated that there were significant spatial trends in SOC stocks (Table 3). Second order trend surface explained 27-39% of the variation of SOC stocks. The trend was removed during interpolation of the results.
Geostatistical analysis. Experimental semivariograms are presented in Fig. 1 and the parameters are shown in Table 4. Based on the selection criteria of highest determination coefficient and lowest residuals, spherical model can best describe the semivariograms for SOC stock at 0-20 cm and 40-60 cm layers, while Gaussian model best described the results for the 20-40 cm and 0-60 cm layers. The determination coefficient ranged from 0.55 to 0.75 with a residual range of 0.0141-0.109. The semivariogram of SOC stock at the 40-60 cm layer showed a larger Nugget effect than at other soil layers, followed by 20-40 cm, 0-60 cm and 0-20 cm layers. The ratio of Nugget/Sill was lowest for the 40-60 cm (27%) layer and highest for the 20-40 cm (42%) layer. Spatially, the variation ranges decreased from 30,900 m at 0-20 cm to 15,800 m at the 0-60 cm layer.

Comparison of OK and IDW.
To check the interpolation performance of OK and IDW, the predicted values were plotted against the measured values ( Fig. 2) The linear model intersected 1:1 line for SOC stock. Before the intersection, the linear model (continuous line) overestimated SOC stock, after intersection SOC stock was underestimated. This result was due to the nature of the algorithms used for parameter estimation which aimed to achieve unbiased prediction of the mean values 32,33 . The values of AME, ME, RMSE and pseudo R 2 are shown in Table 5. The closer the AME, ME and RMSE values are to zero, the better the model performed. ME of SOC stock at different soil layers varied from − 0.94 to − 6.98 Mg ha −1 and determination coefficient of OK and IDW ranged from 0.35 to 0.46, indicating that predicted values derived from OK and IDW slightly underestimated SOC stock, but they were suitable for mapping SOC. Results for OK analysis had a higher determination coefficient and lower AME, ME and RMSE values than IDW for all soil layers, therefore having a better performance. This finding may be attributed to the sampling design and the nature of the algorithm of OK (see the discussion section for more details).

Spatial prediction of SOC stock.
To further compare the differences of OK and IDW for the spatial interpolation of SOC stock at different soil layers, spatial distribution maps were produced using both OK and IDW (Fig. 3); results of which showed a strong spatial variability of SOC stock across the whole study area. Generally, SOC stocks were highest in the northeast of the study area for all soil layers and lowest in the central areas, correlating to city centre locations. For 0-20 cm, 20-40 cm, 40-60 cm and 0-60 cm layers, SOC stock derived from OK fell by 30-77, 19-68, 12-62 and 61-200 Mg ha −1 , respectively. These ranges were lower than those derived from IDW. The distribution patterns for OK and IDW were generally similar, however the absolute values of SOC differed. For example, SOC stock for the 20-40 cm layer in the southern part of the study area varied from 45-60 Mg ha −1 for OK (Fig. 3c) and 35-45 Mg ha −1 for IDW (Fig. 3d).
Total SOC stock prediction. The SOC pool from the different soil layers are shown in Table 6

Discussion
Soil is the largest C pool and accurate estimation of SOC is therefore necessary to assess C sequestration or emission potential caused by global environmental changes 1 34 . This result indicates that the Moso bamboo forest soil is a larger carbon pool compared to the average soil results of other forest types. In south China, driven by fast and high economic benefits due to annual harvests, the area of Moso bamboo is increasing by the rate of 3% per year 35 . The increase of Moso bamboo forest area indicates an increase of SOC stock in Moso bamboo forests under the intensive management. In addition, the carbon stock of (including above-and below-ground) Moso bamboo is dynamically balanced between annual harvest and growth of new bamboos, and one sixth of aboveground biomass is harvested for timber every year 36 , thus leading to a higher ecosystem production than the fast-growing Chinese fir plantations 28 . Our results suggest that Moso bamboo forest soils play an important role in alleviating future climate change. In contrast, deforestation of Moso bamboo forests could lead to significant losses of carbon to the atmosphere, thus causing a negative feedback to climate change. SOC stock results were also higher than those from Moso bamboo forests in Hubei Province (65-99 Mg ha −1 ) 27 , Zhejiang Province (94 Mg ha −1 ) 29 and Jiangxi Province (111 Mg ha −1 ) 28 . The higher SOC results of our   investigation are mainly attributed to fertilization treatment in the majority of Moso bamboo forests in our study area to improve stand production and bamboo shoots; the other three comparison bamboo forests did not receive fertilization treatments [27][28][29] . SOC stock in this study was also higher than the average SOC stock of plantations (104 Mg ha −1 ), and it was comparable to natural forests across China (129 Mg ha −1 ) 37 . This indicates that Moso bamboo forests act as an important role in global C cycling with storing more C in soil. SOC stock results decreased as soil depth increased, a finding which is in general accordance with the majority of previous investigations: SOC in Moso bamboo forests showed a diminishing trend 38,39 . However, some investigations also reported a positive trend of SOC stock with increasing soil depth 40 ; a finding which we attribute to an unusual trend due to mixing effects of tillage 40 .
Many approaches have been developed in the framework of geostatistical analysis techniques for spatial interpolation, such as OK, co-kriging and universe kriging. As the OK technique normally performs better and is easy to apply, it has been widely used in spatial interpolation of SOC stock 3,29 and soil physical properties 30 . To describe the spatial variability of soil properties, semivariograms are fitted and used 3,13,29 .
Nugget values present undetectable experimental errors, field variation within the minimum sampling space and inherent variability 41 . Nugget values in this study were small and varied from 0.0427 to 0.0947, indicating a positive nugget effect, a weak field or random variability 41 . The nugget-to-sill ratio represents a spatial dependence. The ratios < 25%, 25-75% and > 75% suggest a strong, moderate and weak spatial dependence 42 . A strong spatial dependence of soil properties is attributed to soil intrinsic properties, such as soil parent materials, soil texture, topography and vegetation 29,41 . On the other hand, a weak spatial dependence of soil properties indicates that the spatial variability is mainly regulated by extrinsic variations, such as soil fertilization and cultivation practices 41,43 ; moderate spatial dependence is controlled by both intrinsic and extrinsic factors. In this study, the nugget-to-sill ratios ranged from 27% to 42%, demonstrating a moderate spatial dependence of SOC stock for the different soil layers which was controlled by both intrinsic and extrinsic factors. This was evidenced by local, intensive stand management practices, such as tillage and fertilization, widely used in the study area to improve the stand productivity and bamboo shoot output 23,44 , and the complex local topography. A previous study also showed that intensive management led to a decrease of soil organic matter content in Moso bamboo forests 19 . This study also detected a significant positive relationship between SOC stock and topographic variables, notably elevation ( Table 2). This result also supported findings from previous studies which showed that SOC stock increased with elevation due to a reduction in human disturbance and slower decomposition rates of soil organic matter 19,45 . In areas with high rates of disturbance, for example city centres and surround areas, SOC stock was notably lower than areas without disturbance (Fig. 3). Soil temperature, a key determinant of soil organic matter decomposition, decreases with increased elevation which therefore decreases the loss of soil carbon to the atmosphere 45 .
Ranges (A 0 ) indicate different influence zones of environmental factors at different scales, beyond the range, the measured data are not spatially dependent 41,46 , which means that the sampling points cannot be applied for spatial interpolation. This suggests that A 0 can therefore be an effective criterion for evaluating sampling design and SOC stock mapping. The lowest A 0 (15,800 m) was found at the 0-60 cm layer and the largest A 0 (30,900 m) was found for the 0-20 cm layer. These ranges were larger than the sampling intervals (Fig. 1), which suggested that the sampling strategy in this study was appropriate for studying the spatial pattern of SOC stock. However, A 0 differed from soil layers, indicating that the influence zones of SOC stock were not uniform and were depth dependent. This result was mainly associated with intensive human disturbance during timber harvest, digging bamboo shoots and applying fertilizers. These processes lasted for more than six months per year and mainly occurred within the 0-40 cm soil layers where bamboo roots and shoots are widely distributed 27 . ME values in our analysis were negative, indicating both OK and IDW underestimated SOC stock. OK analysis had higher R 2 values and lower AME, ME and RMSE results suggesting this analysis had a more suitable fit than IDW, a finding which is supported by previous studies 14,32,47 . The R 2 values, although being relatively low (0.35 to 0.46) were consistent with previous reports on the spatial interpolation of soil properties 29,48 . This problem is associated with the dataset that soil samples were not collected using a probabilistic sampling design 48 . In our study, due to uneven distribution of bamboo forest, the majority of sampling plots were located in the southern and eastern areas of Yong'an City (Fig. 4). It is also recommended that probabilistic-based design and depth function should be applied to further study of the vertical distribution of SOC stocks since the fertilization treatment could considerably affect the nutrient supply and availability, especially for the top layers, where most of roots are distributed. Another possible explanation for low R 2 values might be the strong local variation of SOC stock due to the variability in the environmental conditions in the study area. However, the low correlation between the predicted and measured data indicated that a better methodology such as Artificial Neural Network  Table 5. Cross-validation indices for ordinary kriging (OK) and inverse distance weighting (IDW) methods. AME = absolute mean error; ME = mean error; RMSE = root mean square error; Pseudo R 2 = pseudo determination coefficient.
Scientific RepoRts | 7:42640 | DOI: 10.1038/srep42640 or Random Forest may improve the accuracy of spatial interpolation of SOC stock, including both intrinsic and extrinsic factors. Generally, OK and IDW produced similar results for the spatial pattern of SOC stock, together with comparable total SOC stock derived from OK and IDW for different soil layers, demonstrating the suitability of OK and IDW for spatial interpolation of SOC stock. However, a larger range was produced by IDW compared to OK, indicating the necessity of using different approaches to study the spatial SOC stock, particularly areas with a complex topography. This issue is related to the different algorithms and computational efficiencies for the spatial interpolation of OK and IDW, thus selection of the appropriate approach is important to improve the interpolation accuracy and efficiency. Theoretically, OK can provide the best linear unbiased estimations and information on the spatial patterns of estimation errors 32 . However, it is important to note that the assumption of stationarity may be not appropriate in practice 49 . The IDW method involves a simple and quick calculation and does not require assumptions about the data. However, IDW does not have the statistical advantages compared to OK (Table 4). The IDW formula has the effect of giving data points close to the interpolation point relatively large weights, while points further away from the interpolation point exert little effect 50   used, the more influence the point close to the estimation value is given. Therefore, as a result of the irregular distribution of the sampling locations, combined with better statistical performances, OK is recommended as a more suitable approach for similar studies in the future. Although OK and IDW generated a similar result of SOC pool of Moso bamboo forests for the study area, SOC pool derived from CA was higher than that of OK and IDW, demonstrating the importance of selecting the appropriate approach to estimate SOC stock. Although CA is a simple approach once the mean SOC stock per unit and area of the study site are known, CA is unable to provide information about the continuous mapping of SOC stock and therefore cannot test the accuracy of spatial distribution of SOC stocks. CA can only provide limited information for optimizing stand management to improve stand productivity, thus making CA of limited use in studying the spatial pattern of SOC stock. Therefore, CA is only recommended to be used when the coordinates of the sampling locations are missing and the spatial patterns of SOC stock are not required; different geostatistical approaches are recommended to be used to obtain accurate and consistent spatial patterns of SOC stock and regional SOC pools.
Compared to other studies on SOC stock in Moso bamboo forests, such as Zhang, et al. 45 and Fu,et al. 29 , this study is the first to attempt to (1) compare different spatial interpolation approaches; and (2) compare geostatistical approaches and CA for regional SOC stock estimates. These results can improve the methodology selection of studying spatial distribution of SOC stocks. In addition, scientific management of Moso bamboo forests requires site specific maps of SOC stock to improve stand productivity. Regarding stand management, this study further proposed that rather than using a consistent treatment of fertilizers across the whole study area, different distribution patterns of SOC stocks indicated different fertilizer treatments should be conducted in different sites since SOC is an important indicator for soil fertility. For example, organic fertilizer treatment could be applied in the centre areas of the southern study regions (Fig. 3). Together with organic fertilizers, the addition of other nutrients, such as nitrogen, phosphorus and potassium, should be added since, nitrogen and phosphorus are the most important limited nutrients in Moso bamboo forests in south China 51 .
Regarding to uncertainties, analysis of random soil cores for the presence of stones and rocks suggested low contents (< 5%), thus we did not correct for gravel content. This could be a source of uncertainty for regional SOC stock estimates, especially in nutrient poor soils. However, data from poor soils is sporadic as the stands were fertilized and managed every year. The estimates of SOC stock from OK and IDW differed by 13% and 16%, respectively, compared to CA. This difference is due to the lack of the relative area weighted mean of CA. Although model efficiency of different models fitted for semivariograms ranged from 55% to 75% (Table 5), R 2 of the model values and predicted values varied from 0.35 to 0.46. This could be an important uncertainty of regional estimates of SOC stock. However, despite the uncertainty in model efficiency, both geostatistical interpolation (OK) and deterministic (IDW) approaches compared and produced similar estimates of regional SOC stock (3.8% difference, Table 6). This result highlighted that the estimates of total SOC stock were accurate.

Conclusions
In this study, OK and IDW were applied to study the spatial interpolation of SOC stock at 0-20 cm, 20-40 cm and 40-60 cm using the measured data from 111 plots in Moso bamboo forests in Yong'an City, subtropical China. OK, IDW and CA were applied to estimate the regional SOC pool. These results can facilitate the accurate estimation of spatial distribution of SOC stock and regional SOC pool.
Spherical and exponential models were selected to describe the spatial pattern of SOC stock. A moderate spatial dependence of SOC stock was observed, indicating that SOC stock was controlled by both intrinsic factors (e.g. soil parent materials and soil texture) and extrinsic factors (e.g. application of fertilizers and tillage treatment).
OK and IDW produced similar spatial patterns of SOC stock, together with comparable SOC pool, indicating the suitability of both approaches in studying the spatial interpolation of SOC stock. However, OK produced a smaller distribution range of SOC stock compared to IDW, highlighting that it is essential to apply different approaches to obtain accurate and consistent results of SOC stock distribution. SOC pool derived from CA was higher than that from OK and IDW, thus CA is not recommended unless coordinates of the sampling locations are missing and the spatial patterns of SOC stock are not required.

Materials and Methods
Study area. The study area was located in the Yong'an City, Fujian Province, China (117°56′ -117°47′ E, 25°33′ -26°12′ N, Fig. 4). The area is characterized by a subtropical southeast monsoon climate with an average annual temperature of 19.3 °C (ranging from − 11 °C to 40 °C) and precipitation of 1600 mm 44,52 . Elevation in the study area spans 580 m to 1605 m above sea level 44,52 . The accumulated temperature of ≥ 10 °C is 4,520-5,800 °C, lasting for 225-250 days and relative humidity is about 80% 44 . The forest cover is 82% with 5.85 × 10 4 ha of Moso bamboo forests 52 . Moso bamboo forests are mainly distributed below 800 m, most of which are pure stands and are seldom mixed with Keteleeria cyclolepis, Cunninghamia lanceolata, Myrica rubra, Choerospondias axillaris, Liriodendron chinense, Schima Superba, etc. To improve the stand output and increase income, fertilizers have been widely applied to most of the Moso bamboo forests.
Soil sampling. Soil samples were collected from the sub-compartment of the forest resource management of Fujian province, China, an area which was established by the local Forest Bureau for soil mapping (Fig. 4). In the targeted sub-compartment, a cluster of three circular plots with a size of 33.3 m 2 were established, and 138 clusters were created in total. However, due to soil sample damage during transportation, soil samples from 111 plots were used for spatial interpolation of SOC stock. In each plot centre, soil samples from three layers (0-20 cm, 20-40 cm and 40-60 cm) were collected. Soil samples were air-dried at room temperature in the laboratory and prepared for sieving through 2-mm and 0.15-mm meshes for SOC content analysis. Identifiable plant residues and root materials were removed during sieving. As the majority of bamboo roots were distributed within the top 40 cm 27 , soil samples to a depth of 60 cm was deemed suitable to meet the research aims of this investigation. To determine bulk density, a cutting ring approach was used in the soil cores 53 . During fertilizer treatment in the plots, identical stones and rocks were removed from the Moso bamboo forests. This resulted in few stones and rocks being found in the cores, therefore correction for gravel content was not undertaken. Information about sample elevation, coordinates, soil depth, soil type, nitrogen content, phosphorus content and bamboo diameters were recorded and determined according to State Forestry Adiministration 53 .
SOC content was analysed using the K 2 Cr 2 O 7 -H 2 SO 4 wet oxidation method 53 . Specifically, 0.1-0.5 gram of air-dried soil was passed through a 0.15-mm sieve and digested with 5 mL 0.8 mol L −1 K 2 Cr 2 O 7 and 5 ml concentrated H 2 SO 4 (1.84 g mL −1 ) for 5 min at 170-180 °C. The digested solutions were then titrated using standardized 0.2 mol L −1 FeSO 4 solution mixed with 15 ml concentrated H 2 SO 4 per liter to prevent oxidization 53 . SOC stock was calculated as 39,54 : where, SOC is the soil organic C concentration (g kg −1 ); BD is bulk density (g cm −3 ); and D is the depth of the soil layer (cm).

Extraction of topographic variables from a Digital Elevation Model (DEM).
A DEM with a resolution of 90 m was obtained from Geospatial Data Cloud (http://www.gscloud.cn/). Mean values of aspect, elevation and slope were extracted for each sample plot in ArcGIS 10.2 (http://www.esri.com/). Further details of the calculation of aspect, elevation and slope are described by Pierce, et al. 55 .
Statistical and geostatistical analyses. Traditional statistical analysis, such as mean, standard deviation and coefficient of variation, were calculated to describe the original data. The relationships between topographic variables (elevation, slope and aspect) were analysed using linear regression. Before starting geostatistical analysis, raw data was initially tested for normality using the Shapiro-Wilk test in R 56 . Instances where the data did not meet the assumption of normal distribution, the raw data was log-transformed and then transformed back using weighting mean in GS + 10.0 (www.gammadesign.com). In this study, OK and IDW were applied to estimate the spatial distribution of SOC stock.
Spatial autocorrelation and trend surface analyses. Moran 58 . In this study, second degree of polynomial surface was used because the increase of degree did not result in a significant increase of determination coefficient and F ratio.

Ordinary Kriging (OK).
Kriging is based on the theory of regionalized variables which are spatially distributed and autocorrelated 59 . The spatial autocorrelation can be indicated by Moran's I (see above). Spherical, exponential and Gaussian models are commonly used to calculate experimental semivariograms using the observed data 49 . The semivariograms are expressed as a function of distance between sampled points and calculate the integrity of spatial continuity in one or multiple directions using the following expression 13 : where, i, z(x i ) and z(x i+h ) are values of z at locations x i and x i+h , respectively; h is the lag and N(h) is the number of pairs of sample points separated by h. In this study, spherical, exponential, linear and Gaussian models were used to describe the semivariograms of SOC stock at 0-20 cm, 20-40 cm and 40-60 cm layers. The models with highest coefficient of determination and the smallest of residuals were chosen. These models where then applied to analyse spatial structure and to provide the input parameters for interpolation. There are three major parameters derived from the fitted models to identify the spatial structure of SOC stock for a given scale. The parameters are nugget (C 0 ), the sill (C + C 0 ) and the range (A 0 ).
The sill (C + C 0 ) represents total variation, and the ratio of nugget and sill is considered as a criterion to classify spatial dependence 3 ; A 0 represents the separation distance, beyond which the measured data are not spatially dependent 46 . More details about the semivariograms and kriging can be found in Goovaerts 49 . The most likely value R(x), which was expected to be encountered in a particular grid cell when using m neighbouring observations, was defined as: j m j i 1 Inverse distance weighting (IDW). Similarly, IDW is another important approach for spatial interpolation which assumes that each point influences the resulting surface to a finite distance 60 . IDW calculated an unsampled point as a weighting average of a known data point within local surroundings. The formula can be Scientific RepoRts | 7:42640 | DOI: 10.1038/srep42640 expressed by Eq. (4) 50 . In this study, data points of 16 without a fixed radius and the weight power of one were used; the weight power of one was found to perform better than the weight powers of two, three and four if the skewness is below one 32 where, r is the weight and d ij is the distance, which is the distance between the estimation point and the measured point.
Data validation. The prediction accuracy of SOC stock was evaluated using the leave-one-out cross-validation techniques 3,33,47 . In the determination of errors, one point was omitted and this value was estimated by the remaining values. Afterwards, the estimated value was compared with the real value in the situation of omitted point 3 . This process was repeated for all the observations. Four commonly used indices, i.e., absolute mean error (AME), mean error (ME), root mean square error (RMSE) and model efficiency (R 2 ), were used to compare the interpolation accuracy for OK and IDW. These indices were calculated as follows: