Spatial patterns and drivers of soil total nitrogen in anthropogenic shrub encroachment in desert steppe

Abstract

Nitrogen is the most important driving factor in primary production and decomposition in arid and semi-arid ecosystems. The effects of shrub encroachment on nitrogen cycling have been investigated at the site scale but seldomly conducted at the landscape scale. Here, we selected 43 shrubland sites distributing across 3000 km² area in temperate desert grassland in eastern Yanchi County of Ningxia Hui Autonomous. We investigated the spatial heterogeneity and driving factors of soil total nitrogen (STN) at the landscape scale by using geostatistical analysis and the geographical detector method. Our results showed that the average soil total nitrogen decreased in the order of 0–5 cm (0.21 g kg⁻¹) > 5–15 cm (0.19 g kg⁻¹) > 15–40 cm (0.18 g kg⁻¹). Geostatistical analysis showed that soil total nitrogen exhibited the strong spatial autocorrelation in the 0–5 and 5–15 cm soil layers and the moderate spatial autocorrelation in the 15–40 cm soil layer. Furthermore, the geographic detector method indicated that soil physicochemical properties exhibited the stronger effects than these of topographic and vegetation biomass in determining the spatial distribution of soil total nitrogen. Specifically, soil water content in the 0–20 cm soil layer explained 35% of variation in soil total nitrogen spatial pattern in the 0–5 cm soil layer, while soil organic carbon content in the 15–40 cm soil layer explained 64% and 45% of variation in soil total nitrogen spatial patterns in the 5–15 cm and 15–40 cm soil layers, respectively. It was concluded that soil water content and organic carbon content primarily drove the formation of soil total nitrogen spatial heterogeneity in shrubland at the landscape scale, indicating that anthropogenic shrub encroachment evidently affected soil water content and redistribution in dryland.

1 Introduction

Shrub encroachment occurs widely in 10%–20% of arid and semi-arid grassland ecosystems, which is usually accompanied by gradual and irreversible changes in plant community composition (Van Auken, 2009; Eldridge et al., 2011; D’Odorico et al., 2012; Nackley et al., 2017; Bestelmeyer et al., 2018; Travers et al., 2019). This phenomenon can be caused by anthropogenic factors (overgrazing of livestock or reduced fire frequency) and climate change (increased atmospheric carbon dioxide concentration or change in rainfall patterns) (Archer et al., 1995; D’Odorico et al., 2012; Devine et al., 2017; Venter et al., 2018). The microclimate formed in shrubland also promotes the further establishment and expansion of shrub (Thompson et al., 2017). In some regions, shrub encroachment can be consciously implemented to introduce native shrub species to overgrazed desert grassland for ecological and vegetation restoration (Zhao et al., 2018; Zhao et al., 2020; Du et al., 2021). Both natural and anthropogenic shrub encroachment have important and unpredictable effects on grassland ecosystem structure, function, and services (Nackley et al., 2017; Bestelmeyer et al., 2018; Travers et al., 2019; Li et al., 2021; Zhao et al., 2021). For example, the formation of excessive shrub in the temperate desert steppe can significantly deplete soil moisture and decrease herbaceous diversity (Yu et al., 2021; Zhao et al., 2021). Moreover, growth of woody plants is at the expense of herbaceous plants, which reduces the productivity of livestock grazing (cattle, sheep, and horses) (Archer et al., 2017; Pierce et al., 2019). Changes in plant communities also lead to distinctly different ecosystem processes at the structural, morphological, and biological levels and significantly alter nutrient fluxes, spatial distribution, and biogeochemical cycles in arid and semi-arid grassland (Throop and Archer, 2008; Ratajczak et al., 2012; Sala and Maestre, 2014; Bestelmeyer et al., 2018; Stanton et al., 2018; Zhou et al., 2019). Consequently, shrub encroachment has become an important social and ecological issue for land managers (Anadón et al., 2014; Belayneh and Tessema, 2017; Ayalew and Mulualem, 2018).

Nitrogen content alters both abiotic (soil nutrients and pH) and biotic (biomass) composition in soil and significantly influences ecosystem functioning (Turnbull et al., 2010). In arid and semi-arid ecosystems, nitrogen is the most important nutrient influencing soil primary production and decomposition (Gebauer and Ehleringer, 2000). Soil total nitrogen (STN) increase, decrease, or no change has been reported along the transition from grassland to shrubland (Boutton and Liao, 2010; Li et al., 2019; Turpin-Jelfs et al., 2019; Li et al., 2021). In addition, STN responded varyingly to vegetation patches in different geographical regions. In subtropical grassland in New Mexico, shrub encroachment drove STN to accumulate beneath shrub patches but not under grass patches, while STN tended to increase under both grass and shrub patches in shrub encroached Tallgrass Prairie and Mediterranean grassland (Hibbard et al., 2001; Lett et al., 2004; Báez and Collins, 2008; Maestre et al., 2009). The varying STN accumulation speeds under shrub and grass patches should also cause varying levels of STN distribution fragmentation at the regional scale or landscape scale. However, such information has been largely lacking as the majority of studies in STN have been conducted at the small scale or site scale. The enriched knowledge of STN spatial pattern at the landscape scale is essential for us to understand the regional flow of STN over large geographical area (Archer et al., 2017; Ward et al., 2018; Connell et al., 2021).

Shrub invasion can affect a broad spectrum of ecological processes at the landscape scale (Pickett and Cadenasso, 1995; Baveye and Laba, 2015). In terms of vegetation, Zuo et al. (2021) found that shrub invasion simplified community composition and amplified the spatial heterogeneity of vegetation. Regarding soil, Eldridge et al. (2015) found that shrub encroachment changed the spatial pattern of soil infiltration, decreased interstitial infiltration connectivity, and reduced infiltration under shrub patches compared to grass patches. Meanwhile, woody plant encroachment can change soil organic carbon, total nitrogen, total phosphorus, and δ¹⁵N within the 1.2 m soil profile at the landscape scale (Zhou et al., 2017a; Zhou et al., 2018a; Zhou et al., 2018b; Zhou et al., 2018c). In the temperate desert steppe, the introduced shrubs accelerated water depletion and drought in soil, resulting in the enhanced spatial heterogeneity and fragmentation of soil moisture in the 0–200 cm soil layer (Zhao et al., 2018; Zhao et al., 2020). Shrub encroachment was also observed to aggravate the spatial heterogeneity and fragmentation of soil organic carbon (SOC) in the soil surface layer (Zhao et al., 2019). However, few studies have been conducted to investigate STN in grassland-shrubland transition, and our understanding of STN spatial pattern in shrub encroached grassland is still largely limited.

Semi-variogram and kriging methods have been frequently used to study spatial autocorrelation and heterogeneity of soil properties at the landscape scale (Western et al., 1998; Maroufpoor et al., 2020). By using these methods, we can predict that spatial pattern of a soil indictor is affected by factors of structure, randomness, or both. However, a major limitation of these two methods is that the interaction between spatial properties cannot be predicted (Western et al., 1998; Yang et al., 2017). To complement this, a new tool, GeoDetector, can be used to not only measure, mine, and exploit spatial heterogeneity but also estimate the degree of interpretation of independent spatial variables and interactions between dependent spatial variables (Wang et al., 2016; Wang and Xu, 2017). Here, all these methods were recruited to analyze and predict the spatial distribution pattern of STN and its corresponding driving factors in grassland-shrubland transition at the landscape scale in desert steppe. Specifically, our objectives were to 1) examine how STN content varied at different soil depths in shrubland; 2) estimate the spatial heterogeneity of STN at different soil depths in shrubland; and 3) determine how topography, vegetation, and soil physicochemical properties determine spatial pattern formation of STN at different soil depths in shrubland.

2 Materials and methods

2.1 Site description

Fieldwork for this study was conducted in eastern Yanchi County of Ningxia Hui Autonomous Region in China (Figure 1; 106°30′-107°39′E, 37°5′−38°10′ N). Situated on a typical desert steppe region on Loess Plateau comprising ≈3,000 km², eastern Yanchi county ranges in elevations of 1,283 to 1,700 m a.s.l. The mean monthly maximum and minimum temperatures in this region were −6.9°C in January and 23.1°C in July from 1986 to 2017, respectively, with a mean annual temperature of 8.1°C. Mean annual precipitation for the same period was 301.2 mm, with 80% of it being received during from July to September. The mean annual sunshine duration, mean annual frost-free period, and mean annual evaporation between 1986 and 2017 were 2,863 h, 160 days, and 2,139 mm, respectively. Soil types in the study area are identified as the sierozem and aeolian sandy soil, which is described as having a low level of soil fertility with loose structure according to Chinese Soil Taxonomy (Zhao et al., 2018). Vegetation is dominated by several native desert grassland species, such as Agropyron mongolicum, Lespedeza potaninii, Artemisia scoparia, Stipa breviflora, Stipa bungeana, Leymus secalinus, and Pennisetum centrasiaticum. This region has been characterized by low vegetation coverage, land desertification, and soil erosion due to overgrazing over the past decades. As a result, since the 1970s and 1980s, a native xerophytic shrub (Caragana korshinskii) has been introduced to this area for vegetation reconstruction (Zhao et al., 2018; Zhao et al., 2019; Du et al., 2021). Other species including Artemisia ordosica and Medicago sativa have been planted in the area but they occupy very little area. To date, shrub cover in this area has anthropogenically increased to more than 48.5% of the total land area over the past 40 years. These shrubs have formed a large number of grass-shrub mosaics (Zhao et al., 2021).

FIGURE 1

2.2 Field sampling and data collection

From July 18^th to August 17^th in 2017, soil was sampled over the whole region. We selected 43 sites in desert grassland where C. korshinskii were planted (Figure 1). The planting volume, planting interval, and age of C. korshinskii were 1,665–3,330 clusters/hm², 2–40 m, and 3–29 years, respectively (Zhao et al., 2018). Three representative plots (10 × 10 m for each plot) were randomly and homogeneously assigned in each of 43 sites. Six soil columns (3 columns in each of shrub patches and shrub interspace) with a diameter of 5 cm and a depth of 40 cm (0–5, 5–15, and 15–40 cm soil layers) without the litter layer were collected in a homogeneous manner in each plot. These soil samples (n = 2,322) were then taken to the lab to be homogenized, air-dried, and sieved through with a 2-mm mech to remove gravel and litter for proper STN measurement. STN content was then determined using the Kjeldahl digestion procedure (Bremner and Tabatabai, 1972). The Digital Elevation Model (DEM) data comes from the Geospatial Data Cloud (http://www.gscloud.cn/) with a resolution of 30 m. The elevation, slope, and aspect data of 43 sites were calculated and derived from the DEM data (Supplementary Table S1). Soil organic carbon content, soil water content, and aboveground vegetation biomass data were obtained from Zhao et al. (2018), Zhao et al. (2019), Zhao et al. (2021). It is worth noting that the soil organic carbon, soil moisture, and aboveground vegetation biomass previously examined by Zhao et al. (2018), Zhao et al. (2019), Zhao et al. (2021) were reused as potential STN influencing factors for a different purpose in the present study (Supplementary Table S1).

2.3 Research methods

2.3.1 Descriptive statistics

The datasets of STN in the 0–5, 5–15, and 15–40 cm soil layers were analyzed to determine the classical statistical parameters, i.e., maximum, minimum, mean, median, standard deviation, and coefficient of variation (CV = mean/standard deviation) (Yang et al., 2017). These parameters are generally used to indicate the midpoint and spread of data. The Kolmogorov-Smirnov method, skewness, and kurtosis were used to evaluate the normality of STN datasets. Log₁₀ transformation was performed on non-normally distributed data for normality correction.

2.3.2 Geostatistical statistics

Based on geostatistical analysis and theory of regionalized variables, we used a semi-variogram with a best-fit model for STN to quantify spatial structure to provide parameters for kriging optimal spatial interpolation (Matheron, 1963; Western et al., 2004; Webster et al., 2007; Maroufpoor et al., 2020). Semi-variograms describe the ratio change in STN between sampling sites using the following equation:where for each site 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 h) is the number of pairs of sample points separated by h. The measured data were used to calculate the experimental semi-variogram and were then fitted by authorized theoretical models, i.e., linear, Gaussian, spherical, and exponential models (Goovaerts, 1999). The best-fitted model was considered to be the one with the smallest residual sum of squares (RSS) and the largest coefficient of determination (R²) (Maroufpoor et al., 2020). Four major parameters could be derived from the fitted model, i.e., the nugget (C₀), the structural C), the sill (C + C₀), and the spatial autocorrelation range (A₀), which could identify the spatial structure of variables at the given scale. The total variance (Sill, C + C₀) is expressed as the sum of structural and nugget variances. To determine the magnitude of spatial dependence, the percentage of the total variance (sill) explained by structural variance C) was calculated using a structural variance ratio [C/(C + C₀)] (Eldridge et al., 2015). A variable is considered to have weak, moderate, or strong spatial dependence if the structural variance ratio is equal to or lower than 25%, between 25% and 75%, or greater than 75%, respectively (Cambardella et al., 1994). The spatial autocorrelation range (A₀) represents the maximum distance within which variables show spatial dependence (Cambardella et al., 1994).

The parameters of the best-fitted model based on semi-variance analysis were used in the ordinary kriging method, and the interpolation map was used to test the predicted distribution of STN (Yang et al., 2017). In addition to testing R² and RSS of the best-fitted model, another way of testing the accuracy of kriging interpolation maps is to calculate the diagnostic statistics of the cross-validation approach. Generally, three statistics: the mean error (ME), the mean square deviation (MSE), and the root-mean-squared error (RMSE) were considered as follows (Webster et al., 2007):where n is the number of the locations; P_i and M_i are the predicted and measured values at location i, respectively, and δ_i is the kriging standard error at location i. The model was optimal if these statistics met the following criteria: ME and MSE were close to 0 and RMSE was close to 1 (Webster et al., 2007).

2.3.3 Geographical detector method

The geographical detector method (GDM) is a new statistical method to detect spatial stratified heterogeneity and reveal the corresponding driving factors (Wang et al., 2010). GDM has no linear hypothesis, while its idea assumes that if an independent variable X) has an essential impact on a dependent variable Y), the spatial distribution of that independent variable X) and that dependent variable Y) should be similar (Wang et al., 2010; Wang et al., 2016). Different from traditional methods, GDM can detect the influencing factors (numerical variables or categorical variables) accounting for STN distribution and determine how any of two factors interactively determine STN distribution (Wang and Xu, 2017). GDM consists of four parts: factor detector, ecological detector, risk detector, and interaction detector (Wang and Xu, 2017). As the core part of the GDM, the q-value is used to check the influencing factors and compare the dispersion variance between the observed values of the whole study area along with the strata variables (Wang et al., 2016; Wang and Xu, 2017). The power of influencing factors on STN distribution can be calculated as follows (Wang and Xu, 2017):Where h = 1, …, L is the number of stratums affecting STN (Y) and its influencing factors(X); N_h and N are the number of units in strata h and all samples in the study area; and ² are the variance of strata h and Y value of the whole region; SSW is the within the sum of squares and SST is the total sum of squares. The q-value indicates that the independent variable X explains 100 × q% of the dependent variable in the range of [0,1] (0 if there is no stratification heterogeneity and 1 if the overall stratification is entirely) (Wang and Xu, 2017). If the value of the q tends toward 1, the better the explanatory power of the variable in explaining the distribution of STN.

Interaction detection was used to test the effect of interaction between different influencing factors (X_S). In other words, this method was used to evaluate if the interaction of two factors (independent variables X1 and X2) will enhance or weaken the explanatory power of the dependent variable Y), or if the factor is simply independent of Y (Wang and Xu, 2017). Specifically, we utilized the analysis method described by Wang and Xu (2017), which begins by calculating the q-values of two factors X1 and X2 to generate q (X1) and q (X2) representing that how much variation that X1 and X2 can explain Y. Moreover, the interaction of the two factors yields the q-value: q (X1∩X2). Finally, the values of q (X1), q (X2) and q (X1∩X2) are compared. Using this method, the results of the interaction between the two factors can be divided into five possible relationships (Wang and Xu, 2017).

2.4 Statistical analysis

One-way analysis of variance (ANOVA) and Duncan’s comparisons test (p < 0.05) were used to determine whether the mean of STN significantly differed among different soil layers. We tested the normal distribution of STN in different soil layers using the Kolmogorov-Smirnov test at p = 0.05 significance level. Non-normally distributed STN was log₁₀ transformed to normal distribution in SPSS 21.0 software. Spatial heterogeneity of STN content at different soil depths was then determined by geostatistical methods. The spatial distribution map of STN content in the 0–40 cm soil layer at the landscape scale was determined using the ordinary kriging tool from ArcGIS 10.2 software. To determine the driving factors of the spatial distribution of STN in the 0–40 cm soil layer, the factor detector and interaction detector were used through a geographic detector method (GDM) within Excel-GeoDetector (http://www.geodetector.org) (Wang and Xu, 2017). Since this method requires the independent variables to be numerical variables and the dependent variables to be categorical variables, potential influencing factors (i.e., soil organic carbon content, soil water content, elevation, slope, and vegetation aboveground biomass) were divided into categorical variables by a classification algorithm. In ArcGIS, those potential influencing factors are then subjected to ordinary kriging interpolation and divided by the natural break classification method. Thereafter, the corresponding coordinate information can be obtained by overlapping analysis using the geographic coordinates of 43 sample sites. Finally, STN content in the 0–40 cm soil layer was imported into the GDM (Wang and Xu, 2017) to show the influences of a single factor and interaction of any of two factors on STN distribution.

3 Results

3.1 Descriptive statistics

Descriptive statistics, including overall variability and distribution of STN in the 0–5, 5–15, and 15–40 cm soil layers are presented in Figure 2. The average STN decreased in the order of 0–5 cm (0.21 g kg⁻¹) > 5–15 cm (0.19 g kg⁻¹) > 15–40 cm (0.18 g kg⁻¹) soil layers, though no statistically significant difference was found for STN among these soil layers (p > 0.05). Figure 2 showed that STN had the moderate level of variability in all soil layers, with variation values varying from 24.46% to 46.98%.

FIGURE 2

3.2 Geostatistical analysis

The STN content in the 0–5, 5–15, and 15–40 cm soil layers conformed Gaussian, exponential, and spherical models as best-fitted models, respectively (see Table 1 and Figure 3). Each best-fitted model has a large coefficient of determination (R²) of 0.557–0.774 and a small residual sum of squares (RSS) of 0.002–0.011 (Table 1). The results of the cross-validation analysis showed that the mean error and mean standard error were close to 0, while the root-mean squared-error was close to 1 (Table 2). These results indicated that these models fitted well (Table 2). In addition, the structural variance ratios [C/(C₀+C)] in the 0–5 and 5–15 cm soil layers were greater than 99.8%, demonstrating the strong spatial autocorrelation. Moreover, the structure variance ratio in the 15–40 cm soil layer was 50.6%, which belonged to the moderate spatial autocorrelation (Table 1). Furthermore, the spatial autocorrelation range (A₀) showed an increasing trend as the soil layer became deeper (Table 1).

FIGURE 3

The spatial distribution pattern of STN content was higher in the north than in the south after shrub encroachment in the desert steppe of eastern Ningxia (Figure 4). In the 0–5 cm soil layer, STN distribution was highly fragmented, reflecting in complex spatial distribution exhibiting dense and irregular patches of different sizes. In the 5–15 cm soil layer, STN distribution was also patchy but showed an evident trend that was high in the north and low in the south. In contrast, STN in the 15–40 cm soil layer was mainly distributed in a strip-shaped pattern, as the contour lines here were sparser and smoother than those in other soil layers. The distribution of STN in this soil layer tended to be higher in the north than in the south, confirming the continuity of STN distribution across different soil layers.

FIGURE 4

3.3 Geographical detector method

Nine potential driving factors of STN were analyzed by GDM (Figure 5). The q-value of all driving factors in the 0–5, 5–15, and 15–40 cm soil layers ranged from 0.06 to 0.35, 0.08 to 0.64, and 0.08 to 0.45, respectively. SWC_0-20 (q = 0.35), SOC_15-40 (q = 0.64), and SOC_15-40 (q = 0.45) had the highest q-value than other factors, meaning that they were the major contributors to spatial heterogeneity of STN in the 0–5, 5–15, and 15–40 cm soil layers. The top potential drivers that influenced the spatial distribution of STN in the 0–40 cm soil layer were soil water content (q = 0.35–0.37) and soil organic carbon (q = 0.27–0.64). The q-value of AGB, elevation, aspect, and slope were below 0.23, so these factors had less impact on STN distribution in the shrub-introduced temperate desert steppe.

FIGURE 5

The non-linear and double factor enhancements caused by two-way interaction of factors on STN distribution were shown in Table 3. For example, the q-value of elevation in the 0–5 cm soil layer as a single factor was 0.06, but q-value could be larger than 0.46 when elevation interacted with one of other factors. Among all potential factors, the interaction of any two factors displayed the higher q-values than those calculated from a single factor, indicating that spatial distribution of STN in the 0–40 cm soil layer was primarily determined by multiple factors other than any single factor. Not surprisingly, the interaction between aspect and soil organic carbon or soil moisture explained major (q = 0.52–0.99) effect on STN distribution, while aspect itself only explained minor effect (q = 0.08–0.11) on STN distribution. In most cases, aspect showed the higher interpretation rates when it interacted with other factors, indicating that aspect was an important topographic factor that should not be ignored.

4 Discussion

Aboveground vegetation can affect spatial heterogeneity of soil properties (Schlesinger et al., 1995; D’Odorico et al., 2010; Pellegrini et al., 2021; Zhou et al., 2021). In grassland-shrubland transition, the aboveground vegetation is swiftly switching from grass to shrub. The current study investigated how invasion of woody vegetation affected spatial heterogeneity of STN in different soil depths in grassland, aiming to enrich our understanding of STN dynamics during such an ecological transition and provide land restoration guidance for land managers.

The STN contents of 43 sites in shrubland in temperate desert steppe decreased along with the increased soil depth after anthropogenic shrub encroachment. This finding is consistent with our expectation that large amount of shrub litter is decomposed in the top soil layer under shrub patches (Van Auken, 2009; Archer et al., 2017; Zhou et al., 2019). Moreover, STN in the 0–5 and 5–15 cm soil layers exhibited the strong spatial autocorrelation (Table 1), which was controlled by structural factors (i.e., meteorology, topography, soil texture) (Zhou et al., 2017b). The surface soil (0–15 cm) is vulnerable to human and animal activities, but this is not the case in the study region after 2003 when this region has been fully enclosed to prohibit human and animal entering (Li et al., 2021). In this area, rainfall is the only water source for soil water supply. In this regard, rainfall became a potential major driver to influence STN (Zhao et al., 2020). The STN in the 15–40 cm layer had the moderate spatial autocorrelation (Table 1), which would be affected by structural and random factors (i.e., microhabitat, soil parent material, texture, or mineralogy, and land management) (Cambardella and Karlen, 1999; D’Odorico et al., 2010; Wang et al., 2009; Liu et al., 2013). Consistent with our findings, other studies have also reported the moderate to strong spatial autocorrelation of STN in the 0–40 cm soil layer (Wang et al., 2009, 2017; Liu et al., 2013). Land use/cover change may be a potential major factor affecting spatial variation of STN, reflecting as transition from desert grassland to shrubland in the present study (Li et al., 2017; Zhao et al., 2017). Our results found that the spatial autocorrelation range of STN increased from 4.12 km in the 0–5 cm soi layer to 45.16 km in the 15–40 cm soi layer, which was evidently smaller than the spatial autocorrelation range of 374–461 km in 382 sampling sites across the Loess Plateau (620,000 km²) (Liu et al., 2013). This is likely due to differences in rainfall, elevation, and slope among different regions, which significantly affect the spatial variation of STN (Liu et al., 2013). A study case in subalpine forest reported a spatial autocorrelation range of 1.08 km in an area of 5.64 km², with topographic factors and normalized difference vegetation index (NDVI) being the major factors explaining more than 40% of STN variation (Wang et al., 2017). These varying results suggest that the spatial heterogeneity of STN is highly dependent on geographic location, land use type, and study scale.

Soil moisture is an important limiting factor, especially in semi-arid desert steppe ecosystems, as it affects vegetation development, carbon and nitrogen cycling, and microbial processes (Li et al., 2021; Yu et al., 2021; Zhao et al., 2021). Shrub encroachment has negative effects on soil moisture, leading to soil desiccation, herbaceous species reduction, and possible further changes of soil nutrient cycling or availability (Eldridge et al., 2011; Archer et al., 2017; Zhou et al., 2019; Zhao et al., 2020). Soil water and soil organic carbon were found to be the main influencing factors controlling the spatial variability of STN after shrub planting, as suggested by the factor detector in GDM (Figure 5). The strong coupling relationship or interaction between soil organic carbon and soil total nitrogen means that a stable stoichiometric equilibrium exists within the carbon-nitrogen cycle (Wang et al., 2017; Zhou et al., 2018a; Chai et al., 2019). Previous studies have shown that topography, vegetation, and NDVI are the main factors determining the spatial heterogeneity of soil properties (Xin et al., 2016). In our study, topographic factors (elevation, slope, and aspect) and aboveground vegetation biomass were not dominant factors for spatial differentiation of STN in the topsoil (all q were less than 0.20) (Figure 5). Effects of two-way interaction of any two of nine factors were greater than these calculated from any single factor (Table 3). These findings indicated that the spatial pattern of STN after shrub introduction in temperate desert steppe was caused by the combination of multiple factors (Wang et al., 2016; Wang and Xu, 2017). Among them, the interaction between aspect and soil moisture or soil organic carbon greatly enhanced the explanatory power of STN spatial distribution (q > 0.76) (Table 3). This result implied that aspect was an essential factor that should not be ignored when it comes to investigate the effects of spatial drivers on spatial distribution of soil properties at the landscape scale.

5 Conclusion

References

• 1

  AnadónJ. D.SalaO. E.TurnerB. L.BennettE. M. (2014). Effect of woody-plant encroachment on livestock production in North and South America. Proc. Natl. Acad. Sci. U.S.A.111, 12948–12953. 10.1073/pnas.1320585111

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  ArcherS.SchimelD. S.HollandE. A. (1995). Mechanisms of shrubland expansion: Land use, climate or CO₂?Clim. Change29, 91–99. 10.1007/BF01091640

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  ArcherS. R.AndersenE. M.PredickK. I.SchwinningS.SteidlR. J.WoodsS. R. (2017). “Woody plant encroachment: Causes and consequences,” in Rangeland systems. Springer series on environmental management. Editor BriskeD. D. (Cham: Springer International Publishing), 25–84. 10.1007/978-3-319-46709-2_2

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  AyalewS.MulualemG. (2018). A review on bush encroachment effect on cattle rearing in rangelands. J. Rangel. Sci.8, 403–415.

  • Google Scholar
  • View reference in article

• 5

  BáezS.CollinsS. L. (2008). Shrub invasion decreases diversity and alters community stability in northern chihuahuan desert plant communities. PLoS ONE3, e2332. 10.1371/journal.pone.0002332

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BaveyeP. C.LabaM. (2015). Moving away from the geostatistical lamppost: Why, where, and how does the spatial heterogeneity of soils matter?Ecol. Modell.298, 24–38. 10.1016/j.ecolmodel.2014.03.018

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  BelaynehA.TessemaZ. K. (2017). Mechanisms of bush encroachment and its inter-connection with rangeland degradation in semi-arid african ecosystems: a review. J. Arid. Land9, 299–312. 10.1007/s40333-016-0023-x

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BestelmeyerB. T.PetersD. P. C.ArcherS. R.BrowningD. M.OkinG. S.SchooleyR. L.et al (2018). The grassland–shrubland regime shift in the southwestern United States: misconceptions and their implications for management. BioScience68, 678–690. 10.1093/biosci/biy065

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BouttonT. W.LiaoJ. D. (2010). Changes in soil nitrogen storage and δ¹⁵N with woody plant encroachment in a subtropical savanna parkland landscape. J. Geophys. Res.115, G03019. 10.1029/2009JG001184

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  BremnerJ. M.TabatabaiM. A. (1972). Use of an ammonia electrode for determination of ammonium in Kjeldahl analysis of soils. Commun. Soil Sci. Plant Anal.3, 159–165. 10.1080/00103627209366361

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  CambardellaC. A.KarlenD. L. (1999). Spatial analysis of soil fertility parameters. Precis. Agric.1, 5–14. 10.1023/A:1009925919134

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  CambardellaC. A.MoormanT. B.NovakJ. M.ParkinT. B.KarlenD. L.TurcoR. F.et al (1994). Field-scale variability of soil properties in central Iowa soils. Soil Sci. Soc. Am. J.58, 1501–1511. 10.2136/sssaj1994.03615995005800050033x

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  ChaiQ.MaZ.AnQ.WuG.ChangX.ZhengJ.et al (2019). Does Caragana korshinskii plantation increase soil carbon continuously in a water‐limited landscape on the Loess Plateau, China?Land Degrad. Dev.30, 1691–1698. 10.1002/ldr.3373

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  ConnellR. K.O’ConnorR. C.NippertJ. B.BlairJ. M. (2021). Spatial variation in soil microbial processes as a result of woody encroachment depends on shrub size in tallgrass prairie. Plant Soil460, 359–373. 10.1007/s11104-020-04813-9

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  D’OdoricoP.FuentesJ. D.PockmanW. T.CollinsS. L.HeY.MedeirosJ. S.et al (2010). Positive feedback between microclimate and shrub encroachment in the northern Chihuahuan desert. Ecosphere1, 1–11. 10.1890/ES10-00073.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  D’OdoricoP.OkinG. S.BestelmeyerB. T. (2012). A synthetic review of feedbacks and drivers of shrub encroachment in arid grasslands. Ecohydrol5, 520–530. 10.1002/eco.259

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  DevineA. P.McDonaldR. A.QuaifeT.MacleanI. M. D. (2017). Determinants of woody encroachment and cover in African savannas. Oecologia183, 939–951. 10.1007/s00442-017-3807-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  DuL.ZengY.MaL.QiaoC.WuH.SuZ.et al (2021). Effects of anthropogenic revegetation on the water and carbon cycles of a desert steppe ecosystem. Agric Meteorol300, 108339. 10.1016/j.agrformet.2021.108339

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  EldridgeD. J.BowkerM. A.MaestreF. T.RogerE.ReynoldsJ. F.WhitfordW. G. (2011). Impacts of shrub encroachment on ecosystem structure and functioning: towards a global synthesis: Synthesizing shrub encroachment effects. Ecol. Lett.14, 709–722. 10.1111/j.1461-0248.2011.01630.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  EldridgeD. J.WangL.Ruiz-ColmeneroM. (2015). Shrub encroachment alters the spatial patterns of infiltration. Ecohydrol8, 83–93. 10.1002/eco.1490

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  GebauerR. L. E.EhleringerJ. R. (2000). Water and nitrogen uptake patterns following moisture pulses in a cold desert community. Ecology81, 1415–1424. 10.2307/177218

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  GoovaertsP. (1999). Geostatistics in soil science: state-of-the-art and perspectives. Geoderma89, 1–45. 10.1016/S0016-7061(98)00078-0

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  HibbardK. A.ArcherS.SchimelD. S.ValentineD. W. (2001). Biogeochemical changes accompanying woody plant encroachment in a subtropical savanna. Ecology82, 1999–2011. 10.2307/2680064

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  LettM. S.KnappA. K.BriggsJ. M.BlairJ. M. (2004). Influence of shrub encroachment on aboveground net primary productivity and carbon and nitrogen pools in a mesic grassland. Can. J. Bot.24, 1363–1370. 10.1139/b04-088

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  LiZ.LiuC.DongY.ChangX.NieX.LiuL.et al (2017). Response of soil organic carbon and nitrogen stocks to soil erosion and land use types in the Loess hilly–gully region of China. Soil Tillage Res.166, 1–9. 10.1016/j.still.2016.10.004

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  LiH.ShenH.ZhouL.ZhuY.ChenL.HuH.et al (2019). Shrub encroachment increases soil carbon and nitrogen stocks in temperate grasslands in China. Land Degrad. Dev.30, 756–767. 10.1002/ldr.3259

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  LiZ.WangH.SunZ.ZhaoY.LiX. (2021). Responses of soil nitrogen to the transition from desert grassland to shrubland in eastern Ningxia, China. Chin. J. Appl. Ecol.32, 1230–1240. 10.13287/j.1001-9332.202104.008

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  LiuZ. P.ShaoM. A.WangY. Q. (2013). Spatial patterns of soil total nitrogen and soil total phosphorus across the entire Loess Plateau region of China. Geoderma197–198, 67–78. 10.1016/j.geoderma.2012.12.011

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  MaestreF. T.BowkerM. A.PucheM. D.Belén HinojosaM.MartínezI.García-PalaciosP.et al (2009). Shrub encroachment can reverse desertification in semi-arid Mediterranean grasslands. Ecol. Lett.12, 930–941. 10.1111/j.1461-0248.2009.01352.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  MaroufpoorS.Bozorg-HaddadO.ChuX. (2020). “Geostatistics,” in Handbook of probabilistic models (Elsevier), 229–242. 10.1016/B978-0-12-816514-0.00009-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  MatheronG. (1963). Principles of geostatistics. Econ. Geol.58, 1246–1266. 10.2113/gsecongeo.58.8.1246

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  NackleyL. L.WestA. G.SkownoA. L.BondW. J. (2017). The nebulous ecology of native invasions. Trends Ecol. Evol.32, 814–824. 10.1016/j.tree.2017.08.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  PellegriniE.BoscuttiF.AlbertiG.CasoloV.ContinM.De NobiliM. (2021). Stand age, degree of encroachment and soil characteristics modulate changes of C and N cycles in dry grassland soils invaded by the N₂-fixing shrub Amorpha fruticosa. Sci. Total Environ.792, 148295. 10.1016/j.scitotenv.2021.148295

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  PickettS. T. A.CadenassoM. L. (1995). Landscape ecology: spatial heterogeneity in ecological systems. Science269, 331–334. 10.1126/science.269.5222.331

  • CrossRef
  • Google Scholar

• 35

  PierceN. A.ArcherS. R.BestelmeyerB. T.JamesD. K. (2019). Grass-shrub competition in arid lands: an overlooked driver in grassland–shrubland state transition?Ecosystems22, 619–628. 10.1007/s10021-018-0290-9

  • CrossRef
  • Google Scholar

• 36

  RatajczakZ.NippertJ. B.CollinsS. L. (2012). Woody encroachment decreases diversity across North American grasslands and savannas. Ecology93, 697–703. 10.1890/11-1199.1

  • CrossRef
  • Google Scholar

• 37

  SalaO. E.MaestreF. T. (2014). Grass-woodland transitions: determinants and consequences for ecosystem functioning and provisioning of services. J. Ecol.102, 1357–1362. 10.1111/1365-2745.12326

  • CrossRef
  • Google Scholar

• 38

  SchlesingerW. H.RaikesJ. A.HartleyA. E.CrossA. F. (1995). On the spatial pattern of soil nutrients in desert ecosystems: ecological archives e077-002. Ecology77, 364–374. 10.2307/2265615

  • CrossRef
  • Google Scholar

• 39

  StantonR. A.BooneW. W.Soto-ShoenderJ.FletcherR. J.BlaumN.McCleeryR. A. (2018). Shrub encroachment and vertebrate diversity: A global meta-analysis. Glob. Ecol. Biogeogr.27, 368–379. 10.1111/geb.12675

  • CrossRef
  • Google Scholar

• 40

  ThompsonJ. A.ZinnertJ. C.YoungD. R. (2017). Immediate effects of microclimate modification enhance native shrub encroachment. Ecosphere8, e01687. 10.1002/ecs2.1687

  • CrossRef
  • Google Scholar

• 41

  ThroopH. L.ArcherS. R. (2008). Shrub (Prosopis velutina) encroachment in a semidesert grassland: spatial-temporal changes in soil organic carbon and nitrogen pools. Glob. Change Biol.14, 2420–2431. 10.1111/j.1365-2486.2008.01650.x

  • CrossRef
  • Google Scholar

• 42

  TraversS. K.EldridgeD. J.ValJ.OliverI. (2019). Rabbits and livestock grazing alter the structure and composition of mid-storey plants in a wooded dryland. Agric. Ecosyst. Environ.277, 53–60. 10.1016/j.agee.2019.03.003

  • CrossRef
  • Google Scholar

• 43

  TurnbullL.WainwrightJ.BrazierR. E.BolR. (2010). Biotic and abiotic changes in ecosystem structure over a shrub-encroachment gradient in the southwestern USA. Ecosystems13, 1239–1255. 10.1007/s10021-010-9384-8

  • CrossRef
  • Google Scholar

• 44

  Turpin-JelfsT.MichaelidesK.BiedermanJ. A.AnesioA. M. (2019). Soil nitrogen response to shrub encroachment in a degrading semi-arid grassland. Biogeosciences16, 369–381. 10.5194/bg-16-369-2019

  • CrossRef
  • Google Scholar

• 45

  Van AukenO. W. (2009). Causes and consequences of woody plant encroachment into Western North American grasslands. J. Environ. Manage.90, 2931–2942. 10.1016/j.jenvman.2009.04.023

  • CrossRef
  • Google Scholar

• 46

  VenterZ. S.CramerM. D.HawkinsH.-J. (2018). Drivers of woody plant encroachment over Africa. Nat. Commun.9, 2272. 10.1038/s41467-018-04616-8

  • CrossRef
  • Google Scholar

• 47

  WangJ.XuC. (2017). Geodetector: Principle and prospective. Acta Geogr. Sin.72, 1116–1134. 10.11821/dlxb201701010

  • CrossRef
  • Google Scholar

• 48

  WangY.ZhangX.HuangC. (2009). Spatial variability of soil total nitrogen and soil total phosphorus under different land uses in a small watershed on the Loess Plateau, China. Geoderma150, 141–149. 10.1016/j.geoderma.2009.01.021

  • CrossRef
  • Google Scholar

• 49

  WangJ.LiX.ChristakosG.LiaoY.ZhangT.GuX.et al (2010). Geographical detectors‐based health risk assessment and its application in the neural tube defects study of the Heshun region, China. Int. J. Geogr. Inf. Sci.24, 107–127. 10.1080/13658810802443457

  • CrossRef
  • Google Scholar

• 50

  WangJ. F.ZhangT. L.FuB. J. (2016). A measure of spatial stratified heterogeneity. Ecol. Indic.67, 250–256. 10.1016/j.ecolind.2016.02.052

  • CrossRef
  • Google Scholar

• 51

  WangT.KangF.ChengX.HanH.BaiY.MaJ. (2017). Spatial variability of organic carbon and total nitrogen in the soils of a subalpine forested catchment at Mt. Taiyue, China. CATENA155, 41–52. 10.1016/j.catena.2017.03.004

  • CrossRef
  • Google Scholar

• 52

  WardD.TrinoggaJ.WiegandK.du ToitJ.OkubamichaelD.ReinschS.et al (2018). Large shrubs increase soil nutrients in a semi-arid savanna. Geoderma310, 153–162. 10.1016/j.geoderma.2017.09.023

  • CrossRef
  • Google Scholar

• 53

  WebsterR.OliverM. A.OliverM. (2007). Geostatistics for environmental scientists. 2. ed.Hoboken, NJ: Wiley.

  • Google Scholar

• 54

  WesternA. W.BlöschlG.GraysonR. B. (1998). Geostatistical characterisation of soil moisture patterns in the Tarrawarra catchment. J. Hydrol.205, 20–37. 10.1016/S0022-1694(97)00142-X

  • CrossRef
  • Google Scholar

• 55

  WesternA. W.ZhouS. L.GraysonR. B.McMahonT. A.BlöschlG.WilsonD. J. (2004). Spatial correlation of soil moisture in small catchments and its relationship to dominant spatial hydrological processes. J. Hydrol.286, 113–134. 10.1016/j.jhydrol.2003.09.014

  • CrossRef
  • Google Scholar

• 56

  XinZ.QinY.YuX. (2016). Spatial variability in soil organic carbon and its influencing factors in a hilly watershed of the Loess Plateau, China. CATENA137, 660–669. 10.1016/j.catena.2015.01.028

  • CrossRef
  • Google Scholar

• 57

  YangY.DouY.LiuD.AnS. (2017). Spatial pattern and heterogeneity of soil moisture along a transect in a small catchment on the Loess Plateau. J. Hydrol.550, 466–477. 10.1016/j.jhydrol.2017.05.026

  • CrossRef
  • Google Scholar

• 58

  YuL.WangH.GuoT.YangQ.SunZ. (2021). Bistable-state of vegetation shift in the desert grassland-shrubland anthropogenic Mosaic area. Acta Ecol. Sin.41, 9773–9783. 10.5846/stxb202012023080

  • CrossRef
  • Google Scholar

• 59

  ZhaoB.LiZ.LiP.XuG.GaoH.ChengY.et al (2017). Spatial distribution of soil organic carbon and its influencing factors under the condition of ecological construction in a hilly-gully watershed of the Loess Plateau, China. Geoderma296, 10–17. 10.1016/j.geoderma.2017.02.010

  • CrossRef
  • Google Scholar

• 60

  ZhaoY.ZhouY.WangH. (2018). Spatial heterogeneities of soil water content under anthropogenic introduced shrub (Caragana Korshinskii) encroachment in desert grassland of the eastern Ningxia, China. Chin. J. Appl. Ecol.29, 3577–3586. 10.13287/j.1001-9332.201811.001

  • CrossRef
  • Google Scholar

• 61

  ZhaoY.DuY.MaY.ZhaoY.ZhouY.WangH. (2019). Soil organic carbon dynamics and the prediction of their spatial changes in response to anthropogenically introduced shrub encroachment in desert steppe of the Eastern Ningxia, China. Chin. J. Appl. Ecol.30, 1927–1935. 10.13287/j.1001-9332.201906.001

  • CrossRef
  • Google Scholar

• 62

  ZhaoY.YuL.ZhouY.MaQ.LuQ. (2020). Soil moisture dynamics and deficit of desert grassland with anthropogenic introduced shrub encroachment in the eastern Ningxia, China. Acta Ecol. Sin.40, 1305–1315. 10.5846/stxb201812152735

  • CrossRef
  • Google Scholar

• 63

  ZhaoY.ZhaoY.WangH.MaY.LiZ. (2021). Response of spatial heterogeneity and threshold value for soil water and aboveground biomass of desert grassland-shrubland anthropogenic transition in desert steppe of Ningxia, China. Sci. Silvae Sin.57, 1–12. 10.11707/j.1001-7488.20211201

  • CrossRef
  • Google Scholar

• 64

  ZhouY.BouttonT. W.WuX. B. (2017a). Soil carbon response to woody plant encroachment: Importance of spatial heterogeneity and deep soil storage. J. Ecol.105, 1738–1749. 10.1111/1365-2745.12770

  • CrossRef
  • Google Scholar

• 65

  ZhouY.BouttonT. W.WuX. B.YangC. (2017b). Spatial heterogeneity of subsurface soil texture drives landscape-scale patterns of woody patches in a subtropical savanna. Landsc. Ecol.32, 915–929. 10.1007/s10980-017-0496-9

  • CrossRef
  • Google Scholar

• 66

  ZhouY.BouttonT. W.WuX. B. (2018a). Soil phosphorus does not keep pace with soil carbon and nitrogen accumulation following woody encroachment. Glob. Change Biol.24, 1992–2007. 10.1111/gcb.14048

  • CrossRef
  • Google Scholar

• 67

  ZhouY.BouttonT. W.WuX. B. (2018b). Woody plant encroachment amplifies spatial heterogeneity of soil phosphorus to considerable depth. Ecology99, 136–147. 10.1002/ecy.2051

  • CrossRef
  • Google Scholar

• 68

  ZhouY.MushinskiR. M.HyodoA.Ben WuX.BouttonT. W. (2018c). Vegetation change alters soil profile δ¹⁵N values at the landscape scale. Soil Biol. biochem.119, 110–120. 10.1016/j.soilbio.2018.01.012

  • CrossRef
  • Google Scholar

• 69

  ZhouL.ShenH.ChenL.LiH.ZhangP.ZhaoX.et al (2019). Ecological consequences of shrub encroachment in the grasslands of northern China. Landsc. Ecol.34, 119–130. 10.1007/s10980-018-0749-2

  • CrossRef
  • Google Scholar

• 70

  ZhouY.TaylorR. J.BouttonT. W. (2021). Divergent patterns and spatial heterogeneity of soil nutrients in a complex and dynamic savanna landscape. J. Geophys Res. Biogeosci126. 10.1029/2021JG006575

  • CrossRef
  • Google Scholar

• 71

  ZuoX.MaoW.QuH.ChenM.ZhaoS.LiuL.et al (2021). Scale effects on spatial heterogeneity of herbaceous vegetation in desert steppe depend on plant community type. Ecol. Indic.127, 107769. 10.1016/j.ecolind.2021.107769

  • CrossRef
  • Google Scholar