A 10-year monitoring of soil properties dynamics and soil fertility evaluation in Chinese hickory plantation regions of southeastern China

Monitoring the temporal and spatial variation of soil properties is helpful to understand the evolution of soil properties and adjust the management method in time. Soil fertility evaluation is an urgent need to understand soil fertility level and prevent soil degradation. Here, we conducted an intensive field investigation in Chinese hickory (Carya cathayensis Sarg.) plantation to clarify the spatial and temporal variation of soil properties and its influencing factors, and to evaluate the change of soil fertility. The results showed that the soil pH and soil organic carbon (SOC) significantly increased from 2008 to 2018, while available nitrogen (AN) significantly decreased from 2008 to 2018. The semi-variance revealed that except available phosphorus (AP), the spatial dependencies of soil properties increased from 2008 to 2018. An increasing south-north gradient was found for soil AN, AP, available potassium (AK) and SOC and a decreasing south-north gradient was found for soil pH. The average soil fertility in the whole area was increased from 2008 to 2018. Our findings demonstrated that the changes of the management measures were the reason for the change of soil properties from 2008 to 2018. Therefore, rational fertilization strategies and sod cultivation are recommended to maintain the long-term development of the producing forest.


Results
Descriptive statistics. After Box-cox transformation, the soil properties all passed the K-S test (K-S P > 0.05) ( Table 1). The coefficient of variation (CV) values ranged from 10 to 130%. According to Fu et al. (2014) reported 47 , CV < 10%, between 10 and 90%, and > 90% indicate low, moderate and high variabilities, respectively. With the exception of AP concentrations for 2008 and 2018, which were highly variable, all other soil properties were moderately variable. The average pH was significantly higher in 2018 than that in 2008 (P < 0.05). From 2008 to 2018, the AN concentration declined by 38% ( Table 1). The concentration of AK in Table 1. Descriptive statistics of the soil attributes. AN available nitrogen, AP available phosphorus, AK available potassium, SOC soil organic carbon. Different letters in the same variable indicate significant differences among years at P < 0.05 level. CV: coefficient of variation; K-S p : significance level of Kolmogorov-Smirnov test for normality. The K-Sp values in brackets were calculated after transformation. Spatial cluster and spatial outlier analysis. The local indicators of spatial association (LISA) maps ( Fig. 2) indicated a significant positive spatial autocorrelations for all soil properties (P < 0.05). The local Moran's I result identified high-high spatial clusters of soil pH in the middle region of the study area, while low-low clusters of soil pH were distributed in the northwest region of the study area (Daoshi town) from 2008 to 2018 ( Fig. 2a-c). On the contrary, high-high clusters of AN, AP and SOC were mainly located in the northwest region of study area ( Fig. 2d-i, m-o). Meanwhile, high-high clusters of AK concentration shifted from northwest to northeast of the study area from 2013 to 2018 ( Fig. 2k-l).
To further describe the spatial structures of soil properties in 2008, 2013 and 2018, we calculated the semi-variances function of each study variable and selected the best-fitted models and their related parameters ( Table 2). The spatial dependencies (C 0 /C 0 + C) for soil pH, AN, AK and SOC were moderate and strong for AK in 2008 and 2013. The spatial autocorrelation for SOC was improved from 2008 to 2018 ( Table 2). The ranged of soil properties varied from 0.16 km (AP) to 40.7 km (pH) in 2008, and from 0.13 km (AN) to 23.73 km (AP) in 2018, respectively.
The spatio-temporal distribution maps of soil properties were revealed by the ordinary Kriging interpolation method based on the semi-variance models for 2008, 2013 and 2018. The concentrations of AN, AP, AK and SOC had similar spatial distribution patterns ( Fig. 3d-o), with high values mainly located in the northwest and northeast parts of the study area, while low values in the central and south regions. However, pH values showed an opposite spatial distribution pattern with a gradually increasing trend from north to south (Fig. 3a-c). Generally speaking, the spatial distributions of soil properties were similar to the above spatial clusters identified by local Moran's I (Fig. 2). Meanwhile, soil properties varied considerably from 2008 to 2018 (Fig. S3). The pH value generally increased, with the largest increase in the western part of the study area in 2008-2018 ( Fig. S3a-c). Contrary to soil pH, soil AN concentration decreased in study regions, and the dropped trend declining from 40 to 200 mg kg −1 (Fig. S3d-f).

Control factors for soil properties.
One way ANOVA analysis indicated that MAP and MAT had a significant influence on the change of pH (Table S2). What's more, altitude has a significant influence on AK and SOC. However, most of soil properties was significantly influenced by anthropogenic factors such as fertilization, weeding, and harvesting methods (Table S2). This showed that the difference of comprehensive management mode, such as management method and intensity, will lead to the change of soil properties.
Comprehensive appraisement of soil fertility. The improved Nemerow method was used to evaluate the soil integrated fertility index (IFI) of Chinese hickory plantation, and the results were shown in Fig. 4. The soil fertility of Chinese hickory plantation was at medium level, but the IFI value increased year by year, which was IFI = 1.096 in 2008, IFI = 1.097 in 2013 and IFI = 1.156 in 2018. In 2008, areas with IFI < 0.9 accounted for 21.5% of the whole study area. With the increase of years, areas with IFI < 0.9 gradually decreased, accounting for 11.5 and 8.6% of the whole study area in 2013 and 2018, respectively (Fig. 4). Compared with the previous two periods, IFI in most areas of the study area was between 1.1 and 1.3 in 2018. These results indicated that the soil fertility in the study area increased gradually with the increase of years, and the heterogeneity of soil fertility in different areas decreased (Fig. 4).

Discussion
In this study, we monitored the temporal and spatial variability of soil properties and soil fertility in hickory plantations under intensive management for 10 years. Geostatistical studies have determined the degree of variability of soil properties in different regions and the degree of spatial dependence of soil properties. Kriging interpolation and exploratory spatial analysis can provide data support for regional soil management policies. Soil pH is a fundamental property that has significantly influences on numerous soil physical, chemical, biological properties and processes. Therefore, it is considered to be a key soil variable in terrestrial ecosystem 48,49 . Long-term monitoring showed that the soil pH decreased significantly from 2008 to 2013 (Table 1). Previous study showed that long-term agricultural activities lead to the decreased of soil pH 50 . It is well known that excessive application of inorganic fertilizer is an important cause of soil acidification 51,52 . NH 4 + ions in the fertilizer replace the basic cations (Ca 2+ , Mg 2+ , K + , Na + ) bound on the soil surface, making it easier to leach from the soil and reducing its buffering effect on acidification. In addition, when an NH 4 + ion is absorbed by the crop, an H + ion is released into the soil solution, causing soil acidification. More serious, Al 3+ increases in concentration in more acidic soils 53 . For hickory, soil acidification increased the risk of canker diseased, leading to mass mortality of hickory plantations 54 . This may be the reason that the decrease of soil pH in 2013.The global problem was recognized by local governments, which instructed farmers drastically reduced the amount of inorganic fertilizer applied, and added soil conditioners (lime, bamboo char and potassium humate) to raise the pH of the hickory soil in 2010 55 . Our study also found that precipitation had a significant effect on soil pH (Table S2). The reason is that precipitation can accelerate the decomposition of litter and promote the return of more cations to the topsoil to supplement the soil exchangeable cationic pool, thereby increasing the soil pH 56,57 . This may be the reason that the increase of soil pH in 2018. Kriging interpolation showed that soil pH value in the northwest part of the study area increased most during this decade. (Fig. S3a-c). This may be due to the steep slope and thin soil layer thickness of hickory plantations in the northwest of the study area, and the rapid continuous and rapid weathering of parent material supplementing rich base cation exchange, thus improving the soil buffer capacity 52 .
Descriptive analysis showed that soil AN decreased significantly from 190 to 71 mg kg −1 in this decades. According to the classification levels of the State Soil Survey Service of China (SSSSC 1996) 58 , the concentration of AN more than 120 mg kg −1 can be considered as high level. At the same time, Kriging interpolation analysis showed that high levels of AN were observed in almost the whole study area in 2008 and 2013 (Fig. 3). Lu 59 concluded that soil AN in southern China is not a limiting nutrient, because the overuse of nitrogen fertilizer leads to the excess of soil AN. Geisseler and Scow's survey of 10 years long-term monitoring sites in China showed that soil pH was significantly reduced by 0.45-2.20 units during 8-25 years when inorganic N, P and K were applied 60 . This was consistent with the negative correlation between soil pH and AN, AP and AK. Zhao et al. 61 found that when the amount of nitrogen applied was greater than the amount required by plants, further nitrogen application would be unfavorable to plant growth and might even limit the uptake of nitrogen by plants. Lei et al. 62 showed that excessive accumulation of nitrogen in the soil would cause an increase in the incidence of canker diseased in Chinese hickory plants, which could lead to plant death (Fig. S4). In addition, in southeast of China, where the temperature is high and the rain is heavy, excessive nitrogen input significantly increased the leaching of reactive nitrogen and gas emissions, causing environmental pollution. Based on this, farmers reduced the amount of fertilizer they applied, and soil nutrients were taken away with the hickory harvest 57 . At the same time, the soil erosion of hickory plantations on the steep slope and exposed surface also caused the loss of nutrients. This significantly reduced the content of available nutrients, especially AN in 2018. www.nature.com/scientificreports/ Compared with pH, AN, AP and SOC, the concentration of AP showed a higher variability (90%-120%). According to LISA map, the distribution of AP presented a serious polarization phenomenon, while the highhigh clusters were mainly distributed in the northeastern and the low-low clusters were mainly distributed in the southeastern of study area (Fig. 2g-i). This was related to the relatively low phosphorus concentration derived from the parent materials. In addition, phosphorus was adsorbed and fixed by soil clay or calcium carbonate, so phosphate was less lost in soil runoff. However, in low soil pH, the active iron and aluminum in the soil increased, and the precipitation of iron and aluminum phosphate was easy to immobilized, resulting in the decrease of the availability of phosphorus 63 . Bruun et al. 64 suggested that additional P fertilizer was not needed because P fertilizer has a strong aftereffect and higher temperature will increase soil phosphorus availability when soil available phosphorus was higher than 10.0 mg kg −1 . Unfortunately, the AP concentration in the most of study area was less than 10 mg kg −1 from 2008 to 2018 indicating a severe phosphorus deficiency (Fig. 3g-i). Therefore, it is necessary to continue to apply phosphate fertilizer in regions with low P to improve or at least maintain hickory yields because of its poor mobility in soil 65 . Tong et al. 66 reported that when AK was higher than 100 mg kg −1 in soil, it could meet the requirements of carbohydrate transportation and fat synthesis during the growth period of Chinese hickory nuts. However, the AK < 100 mg kg −1 were widely found in Daoshi town, Qingliangfeng town and Tuankou town (Fig. 3l). The soil in these areas were still potassium-deficient because potassium fertilizer supply was still below consumption. More specifically, Chinese hickory growing areas in these towns shared a common characteristic: it used to grow on steep slopes greater than 40 degrees and had serious soil erosion (field survey data). This will cause potassium loss with water runoff or leaching to deep layers. Therefore, increasing the application of potassium fertilizer, improving the availability of phosphorus and adjusting the fertilization structure are the key points for the sustainable development of hickory plantations.
SOC is commonly considered as an important indicator for evaluating soil fertility [67][68][69][70] . It can not only maintain and release soil nutrients, but also improve the physical structure of soil [71][72][73] . Table 2 showed that compared with 2008 and 2013, low "Nugget-to-sill" ratio was found in 2018 because the intrinsic soil management factors probably enhanced the spatial correlation due to the weakening of farmer's intensive management after 2010. This was consistent with Chen et al. (2019) that long-term large-scale planting can reduce the spatial variability of soil properties and homogenize soil properties 37 . In our study, the concentration of SOC has significantly increased  (Fig. S2). Moreover, the increase of soil SOC was the largest in Daoshi town and Qingliangfeng Town. The average elevation of these towns is over 500, and some are as high as 1000 m. Bangroo et al. (2017) 74 reported that with the increase of altitude, the decrease of microbial activity resulted in the increase of SOC stability, which was one reason for the increased of SOC. The decrease of annual application of inorganic fertilizers under intensive management reduced the activity of microorganisms and the mineralization of SOC which resulted in the increased of SOC content. The traditional method of collecting Chinese hickory fruit required that the ground be kept clean, so herbicides have been the first choice for farmers to clear undergrowth in the past (Fig S1a,b). However, in recent years, due to the intensification of soil erosion, farmers have gradually changed their harvesting methods from knocking to laying nets (Fig. S1a-d). This method did not require the removal of undergrowth and thus reduced the application of herbicides. This enhanced the stability of soil aggregates 28,75 . Moreover, according to field surveys, white clover, ryegrass and milk vetch were grown in large quantities in these areas. The undergrowth's fine root, root hair, associated mycorrhizal and root secretion of C will directly into microaggregate, provide physical protection against SOC mineralization and effectively increase soil C content 76 . Qian et al. 77 study confirmed that soil organic carbon content in hickory plantations increased by 12.9% after 2 years of planting grass. There was a negative correlation between soil pH and SOC (Fig. 1), which was due to the reduced activity of fungi at lower soil pH, thus reducing the decomposition of organic matter 78 . Table 2 showed that forest age had a significant effect on soil pH and SOC. The reason may be that with the increased of forest age, soil surface litter increased and microorganisms promoted the decomposition of litters, thus the SOC content increase. However, the decomposition of litter would produce acid, which has a negative impact on soil pH. SOC is an important indicator of soil fertility, its depletion can trigger or increase soil degradation processes. Therefore, it is necessary to adopt sod cultivation, increase the application of organic fertilizer and carry out continuous monitoring in areas with low soil organic matter content. Previous study showed that soil fertility varied significantly under different cropping systems 37 . And there was a significant correlation between soil fertility and rice crop yield, suggesting that soil health was crucial for sustainable crop production 65 . In our study, the overall soil fertility of study area belonged to the middle level in 2008-2018. Though soil fertility increased slightly in 2018 due to changes in land management intensity since 2010 (Fig. 4). However, the uneven spatial distribution of soil properties exists (Fig. 3). This required variable fertilization rates to maintain soil nitrogen, phosphorus, potassium and organic carbon balance. Soil AP and AK can be increased by adding additional potassium fertilizer and phosphate fertilizer. Soil conditioner (biochar, potassium humate) are also effective in reducing soil acidification. Sod cultivation will be an effective means to increase SOC. At present, production forests (Chinese Torreya and bamboo) are being intensively managed, including heavy application of chemical fertilizers and herbicides such as glyphosate 19,21 . Such intensive management measures seriously harm the ecological balance and further threaten the health of plants. Therefore, rational fertilization strategies and sod cultivation are recommended to maintain the long-term development of the producing forest.

Materials and methods
Study area. The study area is located in Lin'an city (29°-31° N, 118°-120° E), Zhejiang province, southeastern China (Fig. 5). It is the largest production area of Chinese hickory, accounting for approximately 51% of the nationally planted areas 79 . The Chinese hickory planting densities range from 300 to 375 stems ha −1 , with an average diameter at breast height (DBH) of 12 cm and an average tree height of 8 m 80 . The area is characterized by subtropical monsoon climate with four distinct seasons, with the annual average temperature of 16 °C and annual precipitation of 1628 mm. The annual average daylight hours are 1774 h with 235 d frost-free 81 . It has undulating topography with an elevation range of 150-1000 m 82 . The forest ages range of hickory plantations from 30-100 years. According to Geology online map of Zhejiang Province, China (https:// www. osgeo. cn/ map/ m02cd) and Dong et al. 83 , the soil is derived from 7 major types of parent material, which include sandstone, sand shale, slate, phyllite, rhyolite, granite and quartz porphyry. During the period of 2008-2010, 600-800 kg ha −1 of a compound fertilizer (N:P 2 O 5 :K 2 O, 15:15:15) was applied every year 41 . And herbicide application and mowing were the main methods for controlling understory vegetation. In order to reduce agricultural non-point source pollution, reduce soil erosion and phosphorus pollution, improve soil fertility, the local government issued a document, the application of compound fertilizer was reduced to 150-300 kg ha −1 per year, herbicide was banned, and grass was planted under the forest after 2010.

Field sampling and laboratory analysis.
According to the f distribution characteristics of China hickory plantations of Lin 'an district, a grid of 1 km × 1 km were set up. And the grids containing the hickory plantations was taken as sample plots. A total of 209 sample sites were established in 2008. For each site, a 20 m × 20 m plot was set in the center of each sample plot. Five subsoil samples with a depth of 0-30 cm were collected at 2-2.2 m from the tree trunk (the main distribution area of hickory root system) according to the "Z" shape, which were further mixed into one soil sample with a weight of 1 kg, and brought back to laboratory. For each sampling year, 209 soil samples were collected in situ in July of 2008, 2013 and 2018, respectively (Fig. 5). A portable global positioning system (GPS) was used to record the coordinates and altitude of each sampling location. Information on parent materials and forest age were recorded in 2008. The survey related to management measures (including fertilization, weeding and harvesting methods) of the Chinese hickory plantation regions was also carried out every 5 years. Mean average precipitation (MAP) and mean average temperature (MAT) information of sample plots come from the weather forecast network.
After all soil samples were air-dried, the visible root debris was removed and sieved through 2-mm nylon mesh. A portion of each soil sample was ground with an agate mortar to pass the 0.149 mm nylon mesh, and sealed in an enclosed polyethylene bag. The air-dried samples were prepared for physicochemical analysis. All  59 . Soil pH was determined using a pH meter at soil/water (w/v) ratio of 1:2.5. The soil AN was measured by a diffusion method. Soil AK was extracted with 1 mol L −1 NH 4 Ac, and then measured by an atomic absorption spectrometer. The AP was extracted using 0.5 mol L −1 NaHCO 3 (pH = 8.5) and the P concentration in the extract was determined using the molybdenum-blue method. SOC was determined by the K  84 . None of the data were normally distributed. Therefore, the Box-Cox transformation of soil properties were performed to meet the assumption of normality using Matlab r2019a software. Kernel density estimation was used to estimate the distribution of the soil properties in all sample plots (package stats in R statistical software 4.0.0). One-way ANOVA was used to compare the differences in soil properties in 2008, 2013 and 2018. Pearson correlation analysis was used to identify the correlations between soil pH, AN, AP, AK and SOC 47 . An alpha level of 0.05 for significance testing was used in all statistical analyses, unless mentioned otherwise.
Spatial autocorrelation analyses. Spatial autocorrelation analysis is a statistical method to measure the cluster degree of spatial variables 85 . Moran's I is a commonly used index of spatial autocorrelation, which reflects the similarity between adjacent samples 86 . The global Moran's I was used to describe the soil properties autocorrelation feature over the entire regions (see Supplementary material, Text S1 for detailed information).
Geostatistical analysis. The semi-variance (or variogram) is widely used in geostatistics to quantitatively describe the spatial variability of environmental variables, and this relationship was expressed through an effective variogram model, which can further provide input parameters for spatial interpolation of kriging 87 (see Supplementary material, Text S2 for detailed information). The ordinary Kriging method can be used to derive the optimal linear unbiased estimate of spatial variables 88 . The models that fit the semivariogram best according to the regression coefficient were determined. For the kriging interpolation, the transformed soil properties data were used. The ordinary kriging method was used to www.nature.com/scientificreports/ draw a spatial distribution map of soil properties and soil fertility index with ArcGIS desktop 10.7 (https:// www. esri. com/ en-us/ arcgis/ produ cts/ arcgis-deskt op/ resou rces).
Soil fertility evaluation. The calculation of Integrated Fertility Index (IFI) comprises three steps: (1) the selection of indicators, (2) the calculation of the individual fertility index (IFI i ), and (3) the calculation of IFI. The soil pH, AN, AP, AK and SOC were used in the calculations. To calculate the IFI the following equation was used: where IFI i is the individual fertility index; x is the measured value of each attribute 89 ; x a x b and x c are the upper and lower limits of each classification standard based on forest soils in Zhejiang. In the process of carrying out the cultivated land soil survey and evaluation project, soil scientists developed the soil property classification standard of Zhejiang Province. Soil properties such as pH, SOC, AP and AK are the main controlling factors of soil fertility, while AN is not considered to be the main controlling factor 59 . And according to the requirements of nut species on soil acidity, soil nutrient. The classification criteria and separate weights were determined as shown in Table S1 90 . The final step was to calculate IFI using the improved Nemerow Quality Index equation: where IFI is the soil integrated fertility index; IFI iave is the average values for the individual fertility indices; IFI imin is the minimum value for the individual fertility indices; n is the number of soil properties 91 . The degree of IFI was classified as follows: IFI < 0.9; low, 0.9 ≤ IFI < 1.8; moderate, 1.8 ≤ IFI < 2.7; high, and IFI ≥ 2.7; very high. (1) x > x c (2) IFI = 1 2