Spatial Association of Shrubs and Their Interrelation to Burrowing Site Preference of Subterranean Rodents on Dune Slope in the Otindag Sandy Land , China

Rangelands worldwide have more shrubs now, and subterranean rangeland rodents show close interaction to shrubs when choosing a burrowing site. The study was conducted in Otindag Sandy Land in Inner Mongolia, China with the objective of determining the effects of slope position on spatial pattern and interaction of shrubs; how rodents choose their habitat in different slope; and shrubs and rodents influence each other. To accomplish the objective set, we used three physiographic units: Plot 1 (upper slope), Plot 2 (middle slope), and Plot 3 (lower slope), and all individual woody plants and rodent holes in the three plots were mapped. The result of the study showed that: (1) two shrub species show a random distribution trend in all three plots except an aggregated trend only at the smaller scale on the upper slope; (2) the majority of subterranean rodents preferred to select their burrowing sites under the shrub crown, and these selected shrub individuals had generally larger crown length than those unselected individuals. At the same time, the majority of these burrowing sites were located on the lower right direction. (3) The distribution of rodents holes differ across the slopes in the study area. In the three samples, the relative locations of burrowing sites to shrubs are mostly distributed down slope of shrubs. From upper slope to lower slope, this trend gradually enhanced. Our conclusion is that the increase in shrubs represents a pioneer phase in the rehabilitation of degraded sandy land ecosystems, and colonization of subterranean rangeland rodents near the shrubs is a clear indicator of stabilization of sand dunes.

We suggest that the role of shrub species is important in the rehabilitation of degraded sandland ecosystems, rather than as a desertification indicator.Many shrub species are the preferred plant materials for large-scale stabilization of sand dunes.Some shrub species regenerate naturally on the bare sand dunes following cessation of grazing or other human disturbance so that the bare sand dunes became gradually stabilized by this vegetation [15,29,30].These shrubs often form shrubland or shrub-hummock on the slopes of semi-or completely stabilized sand dunes, which can provide shelter for the survival of other plants (i.e., nurse plant syndrome) [31,32] and eventually facilitate more complex vegetation restoration on sand dunes [33].
Shrubs can also provide a food source and potential microhabitat for some small animals, especially rodents [34].Rodents as imperative consumers of plant materials are definitely a significant component of desert ecosystems [35,36].For instance, subterranean rodents play a very important role in desert ecosystem sustainability processes through seed predation and dispersal, and nutrient cycling [37][38][39][40].Subterranean rodents can also increase resource and landscape heterogeneity through burrowing activity [41,42] and a great deal of research about habitat selection by subterranean rodents has been undertaken [43][44][45][46][47].However, burrowing site preferences in desert ecosystems have not received much attention [48].
In this study, a typical slope of a stabilized sand dune in the Otindag Sandy Land was selected for plant and rodent survey.The following questions are addressed based on different slope positions of this sand dune: (1) How are shrubs spatially distributed?(2) How is burrowing site preference of subterranean rodents related to shrub distribution pattern?
Our objective was to better understand the effects of slope position on spatial pattern and interaction of shrubs, and relation between shrubs and rodent selection on the slope of stabilized sand dunes, which will be helpful in managing rehabilitated ecosystems for stability and sustainability.The present paper summarizes our efforts over a one-year period in sites within the Otindag Sandy Land to quantify and explain the spatial distribution of rodent burrows and the implications of this information for rehabilitation of sandy land in Inner Mongolia.

Ethics Statement
This study was conducted in the Otindag Sandy Land, which belongs to the Xilinguole League (115 • 16 E, 42 • 50 N), in the desert steppe region in Inner Mongolia, China.The study plots were selected on private land.We obtained permission to do this study on private land by paying money to the landowner.The permission also requires that open flames cannot be used to prevent forest fires; tall trees cannot be cut; and the forest ecosystem must be protected from any damage.No rare or endangered wild animals or plants are involved in this experiment.Furthermore, the samples in this study only contain shrubs and burrowing sites by rodents, thus no samplings can directly impact vertebrate survival.There is no threat to the environment from this experiment, as no wild animals or plants were research objects.

Study Site
The research area is located in the Otindag Sandy Land, which is administered by the Xilinguole League (115 • 16 E, 42 • 50 N), in the desert steppe region in Inner Mongolia, China.The elevation of this research area is about 1320 m.This region has a continental climate: the mean annual precipitation is 250-350 mm, most of which falls from June to August; there are hot summers; and long and cold winters.With the extreme minimum temperature of −38 • C, the annual mean temperature in this area is 1.7 • C. The annual sunshine duration amounts to more than 1000 h.With 4 m −s of annual mean wind speed, the wind level exceeding Beaufort scale value of 8 is 90 days per year, and the dominant wind direction is from the northwest.The main soil type is classified as an aeolian sandy soil with a mean depth of 200 cm.The field investigation shows that the calcic horizon occurs at 30-100 cm.The horizon in this area is so hard that it is difficult for plant roots to penetrate.
A typical stabilized sand dune was selected for investigation.Its leeward slope with a length 300 m and average gradient 15 • runs in a southeasterly direction.This dune had been recorded as "shifting" early in the 1980s and grazing was excluded at the end of that decade to allow revegetation.Two shrub species have become established on the slope: Spiraea aquilegifolia Pall and Caragana microphylla Lam.
S. aquilegifolia as a long-lived xero-mesophytic shrub widely distributed over forest, steppe and desert along the whole Inner Mongolia Plateau from east to west.The flowering of this shrub mainly happens around June and seeding appears in August or September.C. microphylla, as another dominant species in this area, is also a long-lived xero-psammophytic shrub, scattering in steppe and desert communities along the whole Inner Mongolia Plateau from east to west.It flowers in May and then lasts almost 20 days; the seeds ripen in July; and summer rainfall gives rise to seed germination [49][50][51].Two small subterranean rodent species occur in the research area, i.e., Microtus gregalis Pall and Citellusdauricus (Spermophilus dauricus) Brandt [52], which are similar in body size.

Data Collection
Three physiographic units were setup down the slope of the sand dune, namely Plot 1 (upper slope), Plot 2 (middle slope) and Plot 3 (lower slope).Each plot (50 m × 40 m) was divided into 80 contiguous 5 m × 5 m sub-plots to provide the basis of the vegetation survey.We recorded location coordinates of each shrub, stem position and orientation, crown diameter, height, and health status of all shrubs in all three plots.At the same time, the coordinates of burrowing sites of subterranean rodents were also recorded.The locations of all shrubs and digging sites by rodents in three plots were recorded and mapped using the total station transit (modelGTS-3B, Topcon, Paramus, NJ, USA) within an accuracy range of approximately 1 cm.
To clarify the intraspecific interaction of S. aquilegifolia shrubs, we classified them into two categories, namely, larger shrubs (≥40 cm) and smaller shrubs (<40 cm) (Table 1), according to (height (H) + long crown (LC) + short crown (SC))/3.In addition, rodent burrowing sites were recorded as under a shrub crown or in bare soil surfaces.Therefore, the shrubs were divided into those with burrowing sites (SwBS) or those without burrowing sites (Sw/oBS) for characterizing the relation between shrubs and rodents.

Spatial Association of Shrubs
Within two common spatial patterns analysis techniques, Ripley's K-function and the pair correlation function g(r), we selected g-function in this study, as, compared to the K-function, g-function is more sensitive to small-scale effects.The g-function, a derivative of Ripley's K-function, replaces estimating the number of points within a radius, and analyzes the mean number of neighbors within concentric rings at a distance r, therefore isolating specific distance classes [53,54].Another advantage of the pair correlation function g(r) illustrated by Jian Zhang et al. is that "with increasing circles, g-function is more easily and intuitively analyzes the spatial patterns derived from ecological processes than K-function, by basing on the frequency of points co-occurring at a given distance" [54].K-function is an accumulative measure that cannot distinguish the effects accurately from large scales to those small scales.
To quantify the spatial patterns of individuals within communities in research area, univariate pair correlation functions g(r) was used to detective the single shrub species, while bivariate pair correlation functions g 12 (r) was used to analyze the interactions between two shrub species [53].The following formula is used for univariate analysis [53,[55][56][57]: where A represents the plot area, n means the total number of plants, and w ij means the weighting factor correcting of edge effects.A kernel function, k h , which is also called bandwidth parameter, is used for applying maximum weight to point pairs within a distance t.At a given distance, r, g(r) > 1 indicates the points within distance r are relatively more frequent than expected under complete spatial randomness (CSR), which reflected these points have a typical clustering trend [53].Otherwise, g(r) < 1 means that the points within distanced r are relatively less frequent than expected under CSR and show a regularity pattern.
The following formula is used for bivariate analysis [55][56][57]: where x i (i = 1, . . ., n 1 ,) and y j (j = 1, . . ., n 2 ) represent the points of Groups 1 and 2, respectively, with the weights w ij and kernel function k h as mentioned above.Under the given distance r, g 12 (r) = 1 means there is no interaction between Species 1 and 2. g 12 (r) <1 means Species 2 has a negative associated relationship with Species 1 at given distance.On the contrary, g 12 (r) > 1 indicates Species 2 is positively related to Species 1 at given distance r.
In this study, after comparing and analyzing the observed summary statistics, we chose the null model of complete spatial randomness (CSR) as an appropriate null hypothesis for the univariate analyses of two shrub categories (S. aquilegifolia, and C. microphylla), along age classes (larger and smaller).
In terms of bivariate analyses, the interaction between larger and smaller shrubs was tested.Since smaller shrubs would be influenced by distribution pattern of the larger shrubs in the same area, such as interspecies competition for the same resource in ecological system, two size classes of shrubs were analyzed using bivariate g-function patterns analysis through the toroidal shift and the antecedent condition null model.Firstly, the antecedent condition model can show whether the distribution of smaller shrubs occurs more or less frequently depends on their larger neighbor shrubs or it is an independent process.Secondly, this function also detects the interaction between different shrubs categories.In this research area, it seems that drought stress and habitat heterogeneity are two key elements that affect the spatial distributions of different shrub species.Thus, independent null model [56] was conducted to examine the spatial association between two shrub species.
To test the significance of the point pattern departures from the null model, the Monte-Carlo approach was used to test the significance and generate a summary statistic.Each summary statistic was generated using an approximate (two-sided) 95% simulation envelope after calculating for each distance r the 5th lowest and highest values from 199 simulations of a point process.All analysis processes underlying the null model were finished using the software Programita [58].

Burrowing Site Preference of Subterranean Rodents to Shrubs
To analyze the burrowing site preference of subterranean rodents to shrubs, we firstly compared the growth parameters of SwBS to Sw/oBS; then, for SwBS, we mapped the relative location of burrowing sites to shrubs with a relative origin (0, 0) for each shrub to analyze location preference; and, finally, for Sw/oBS, shrubs with nearest neighbor distance to burrowing sites (NN-Sw/oBS) were selected and their growth parameters were compared to the rest of Sw/oBS (R-Sw/oBS) and to SwBS.We conducted the ANOVA analyses with SPSS v17.0 software.All means are reported with ±SE and the rejection level for H 0 was set at p ≤ 0.05.

Univariate Analysis of Two Shrub Species
As shown in Figure 1 Thus, independent null model [56] was conducted to examine the spatial association between two shrub species.
To test the significance of the point pattern departures from the null model, the Monte-Carlo approach was used to test the significance and generate a summary statistic.Each summary statistic was generated using an approximate (two-sided) 95% simulation envelope after calculating for each distance r the 5th lowest and highest values from 199 simulations of a point process.All analysis processes underlying the null model were finished using the software Programita [58].

Burrowing Site Preference of Subterranean Rodents to Shrubs
To analyze the burrowing site preference of subterranean rodents to shrubs, we firstly compared the growth parameters of SwBS to Sw/oBS; then, for SwBS, we mapped the relative location of burrowing sites to shrubs with a relative origin (0, 0) for each shrub to analyze location preference; and, finally, for Sw/oBS, shrubs with nearest neighbor distance to burrowing sites (NN-Sw/oBS) were selected and their growth parameters were compared to the rest of Sw/oBS (R-Sw/oBS) and to SwBS.We conducted the ANOVA analyses with SPSS v17.0 software.All means are reported with ±SE and the rejection level for H0 was set at p ≤ 0.05.

Univariate Analysis of Two Shrub Species
As shown in Figure 1

Univariate Analysis of S. aquilegifolia
As shown in Figure 2, larger S. aquilegifolia had significantly aggregated trend at 0-1 m scale in Plot 1, and uniform distribution at 0-1 m scale in Plot 2 and at 0-1.5 m scale in Plot 3, while randomly distributed at other scales in all three plots.Smaller S. aquilegifolia were significantly aggregated at 0-2 m scale in Plot 2 and at 0-1.75 m scale in Plot 3; randomly distributed at other scales in Plot 2 and Plot 3; and randomly distributed at all scales in Plot 1.

Univariate Analysis of S. aquilegifolia
As shown in Figure 2, larger S. aquilegifolia had significantly aggregated trend at 0-1 m scale in Plot 1, and uniform distribution at 0-1 m scale in Plot 2 and at 0-1.5 m scale in Plot 3, while randomly distributed at other scales in all three plots.Smaller S. aquilegifolia were significantly aggregated at 0-2 m scale in Plot 2 and at 0-1.75 m scale in Plot 3; randomly distributed at other scales in Plot 2 and Plot 3; and randomly distributed at all scales in Plot 1.

Bivariate Analysis of Intra-and Inter-Specific Shrub Distribution
As shown in Figure 3

Bivariate Analysis of Intra-and Inter-Specific Shrub Distribution
As shown in Figure 3

The Analysis of Burrowing Site Preference of Subterranean Rodents
On the upper slope, there are 14 burrowing sites, among which 42.8% were found under shrub crown and 57.2% found in bare surface; on the middle slope, there are 21 burrowing sites, among which 80.9% were found under shrub crown, and only 19.1% found in bare surface; and, on the lower slope, there are 39 burrowing sites, among which 89.7% were found under shrub crown and only 10.3% found in bare surface (Table 1 and Figure 5).

The Analysis of Burrowing Site Preference of Subterranean Rodents
On the upper slope, there are 14 burrowing sites, among which 42.8% were found under shrub crown and 57.2% found in bare surface; on the middle slope, there are 21 burrowing sites, among which 80.9% were found under shrub crown, and only 19.1% found in bare surface; and, on the lower slope, there are 39 burrowing sites, among which 89.7% were found under shrub crown and only 10.3% found in bare surface (Table 1 and Figure 5).
In all three plots, the relative locations of burrowing sites to shrubs are mostly distributed down slope of shrubs, with 67.7%, 88.8% and 96.5% of burrowing sites down slope and 50%, 85.7% and 87.2% located at lower right direction in Plot 1, Plot 2 and Plot 3, respectively (Figure 6).

The Analysis of Burrowing Site Preference of Subterranean Rodents
On the upper slope, there are 14 burrowing sites, among which 42.8% were found under shrub crown and 57.2% found in bare surface; on the middle slope, there are 21 burrowing sites, among which 80.9% were found under shrub crown, and only 19.1% found in bare surface; and, on the lower slope, there are 39 burrowing sites, among which 89.7% were found under shrub crown and only 10.3% found in bare surface (Table 1 and Figure 5).In all three plots, the relative locations of burrowing sites to shrubs are mostly distributed down slope of shrubs, with 67.7%, 88.8% and 96.5% of burrowing sites down slope and 50%, 85.7% and 87.2% located at lower right direction in Plot 1, Plot 2 and Plot 3, respectively (Figure 6).
In the analysis of three types of shrubs, the height, and long and short crown showed no differences in Plot 1 (Figure 7); the height showed no difference in Plot 2 and Plot 3; and, compared to the other shrubs, the long and short crown (SwBS and SwoBS) showed significant difference in Plot 2 and Plot 3 (p ≤ 0.05).For comparison of growth parameters of SwBS to Sw/oBS, there were significant differences for the long and short crown in Plot 2 and Plot 3 (p 0.05).
Sustainability 2017, 9, 1729 8 of 14 In the analysis of three types of shrubs, the height, and long and short crown showed no differences in Plot 1 (Figure 7); the height showed no difference in Plot 2 and Plot 3; and, compared to the other shrubs, the long and short crown (SwBS and SwoBS) showed significant difference in Plot 2 and Plot 3 (p ≤ 0.05).For comparison of growth parameters of SwBS to Sw/oBS, there were significant differences for the long and short crown in Plot 2 and Plot 3 (p ≤ 0.05).In the analysis of three types of shrubs, the height, and long and short crown showed no differences in Plot 1 (Figure 7); the height showed no difference in Plot 2 and Plot 3; and, compared to the other shrubs, the long and short crown (SwBS and SwoBS) showed significant difference in Plot 2 and Plot 3 (p ≤ 0.05).For comparison of growth parameters of SwBS to Sw/oBS, there were significant differences for the long and short crown in Plot 2 and Plot 3 (p ≤ 0.05).

Alternative Stable State Formed by Spatial Self-Organization of Shrub Species
Spatial pattern of vegetation is determined by a variety of environmental factors at different scales [59,60], the zonal climate determines the vegetation distribution pattern at regional and global scale [61][62][63], while the azonal factor is one of the leading factors of vegetation distribution pattern at the landscape or finer scales [64].The slope position, as an important azonal factor, on the one hand, has direct effects on vegetation through the geomorphologic processes; on the other hand, controls space redistribution of resources factor through change of forms (ups, downs, etc.), thus indirectly affects the distribution of vegetation [65,66].The influence of slope position on vegetation pattern is receiving growing attention [67][68][69][70][71][72][73].In the mountains and hilly regions, the terrain controls the redistribution of solar radiation and precipitation, which can well indicate the microclimate condition of local habitat, and reflect the spatial differences on the thickness and soil nutrient [74-76], hence, the research on the relation between plants and topography has focused on the mountains [77][78][79][80] and hilly regions [65,[69][70][71]81,82].However, fewer researchers focus on the deserts and sand dunes [83,84].In our study, two shrub species show a random distribution trend in all three plots, except an aggregated trend only at the smaller scale on the upper slope (Plot 1), which results from runoff and sediment redistribution driven by microtopography in water-limited semi-arid ecosystem, and corresponding self-organization of vegetation pattern [69,[85][86][87][88]. From the upper slope to the lower slope, the individual number of S. aquilegifolia plants increased and that of C. microphylla decreased (Table 1), while the competition trend either intraspecific or interspecific increased (Figures 3 and 4), which would facilitate formation of spatially periodic vegetation patterns that are well known in arid and semi-arid regions around the world [89,90].C. microphylla, as an encroaching shrub species into xero-mosophytic shrub communities [16,91], can be an indicator of long-term climate change but its amount and its relationship with S. aquilegifolia show that an alternative stable state for S. aquilegifolia-dominated community has not changed due to encroachment of C. microphylla [92,93]: a stabilized sand dune cannot be destabilized due to vegetation self-patterning.

The Nurse-Protégé Interactions between Shrubs and Rodents
Nurse plant syndrome is a commonly accepted concept for positive interaction in plant communities [94].Extensively, some plant recruitment is also facilitated by elements of the microrelief such as stones, hollows or crevices, which are called as "nurse objects".It is well known that microclimatic amelioration and escape from seed and seedling predators are the most common mechanisms underlying plant-plant or object-plant facilitation [95][96][97].However, facilitation by nurse plants to rodents is poorly understood, especially on rodent selection in sand ecology system.Here, we used "nurse-protégé interaction" [31,98] to describe plant-animal relationships between shrub species and rodents.
In our study, the majority of subterranean rodents preferred to select their burrowing sites under the shrub crown, and these selected shrub individuals had generally larger crown length than those unselected individuals (Figures 6 and 7).As Figure 6 shows, the majority of these burrowing sites were located on the lower right direction (i.e., in a northeastern direction).As a protégé species, this preference for burrowing site selection was apparently to avoid the negative impact of rainfall and runoff from upslope on burrowing site and even underground holes [99,100], and to avoid higher surface temperature for rodents in the hot summers [101,102].This preference will also provide more hidden sites for rodents to reduce predation risk [103,104].Obviously, a typical nurse-protégé interaction between shrubs and rodents was formed in the alternative stable community.In addition, recovery of rodent populations after natural or human disturbance could be an indicator for degraded ecosystem restoration [105,106]; hence, nest building of rodents around the shrub species is also a strong evidence for alternative stable state on this stabilized sand dune slope.
In conclusion, this S. aquilegifolia-dominated community on the stabilized sand dune slope maintains an alternative state, which has not changed due to the encroachment of C. microphylla, and unidirectional facilitation from shrub species to rodents (nurse-protégé association) is formed.Further research should focus on the impacts of competition via seed feeding by rodents and facilitation via seed dispersal by rodents on ecosystem self-organization and resilience in semi-arid sandland, which will be more helpful for us to understand the processes of sand dunes, and to formulate scientific restoration strategies for degraded ecosystem.
, S. aquilegifolia had a significantly aggregated trend at 0.75-3.25 m scale in Plot 1, and uniform distribution at 0-1 m scale in Plot 2 and at 0-1.5 m scale in Plot 3. C. microphylla was significantly aggregated at 0-2 m scale in Plot 1, at 0-1 m scale in Plot 2 and at 0-1.5 m scale in Plot 3. The two shrub categories were randomly distributed at other scales in all three plots.Sustainability 2017, 9, 1729 5 of 14 , S. aquilegifolia had a significantly aggregated trend at 0.75-3.25 m scale in Plot 1, and uniform distribution at 0-1 m scale in Plot 2 and at 0-1.5 m scale in Plot 3. C. microphylla was significantly aggregated at 0-2 m scale in Plot 1, at 0-1 m scale in Plot 2 and at 0-1.5 m scale in Plot 3. The two shrub categories were randomly distributed at other scales in all three plots.

Figure 1 .
Figure 1.Univariate spatial patterns of two shrub species on different slope position.

Figure 1 .
Figure 1.Univariate spatial patterns of two shrub species on different slope position.

Figure 2 .
Figure 2. Univariate spatial patterns of S. aquilegifolia on different slope position.
, larger S. aquilegifolia and C. microphylla had significantly negative correlations at 1.75-4.5 m, 0-12.75 m and 0.25-3 m scales in Plot 1, Plot 2 and Plot 3, respectively; smaller S. aquilegifolia and C. microphylla had no correlation at all scales in Plot 1, and had a significantly negative correlations at 1.75-9.25 m scale and at 1.25-8.75m scale in Plot 2 and Plot 3, respectively.As shown in Figure 4, larger shrubs and smaller shrubs show a significantly positive correlation at 1-2 m scale in Plot 1, and significantly negative correlation at 3.25-7.75m scale and at 0-2 m scale in Plot 2 and Plot 3, respectively.

Figure 2 .
Figure 2. Univariate spatial patterns of S. aquilegifolia on different slope position.
, larger S. aquilegifolia and C. microphylla had significantly negative correlations at 1.75-4.5 m, 0-12.75 m and 0.25-3 m scales in Plot 1, Plot 2 and Plot 3, respectively; smaller S. aquilegifolia and C. microphylla had no correlation at all scales in Plot 1, and had a significantly negative correlations at 1.75-9.25 m scale and at 1.25-8.75m scale in Plot 2 and Plot 3, respectively.As shown in Figure 4, larger shrubs and smaller shrubs show a significantly positive correlation at 1-2 m scale in Plot 1, and significantly negative correlation at 3.25-7.75m scale and at 0-2 m scale in Plot 2 and Plot 3, respectively.

Figure 4 .
Figure 4. Bivariate spatial association for intra-specifics.S. aquilegifolia shrubs with both the toroidal shift and antecedent condition null models.

Figure 4 .
Figure 4. Bivariate spatial association for intra-specifics.S. aquilegifolia shrubs with both the toroidal shift and antecedent condition null models.

Figure 4 .
Figure 4. Bivariate spatial association for intra-specifics.S. aquilegifolia shrubs with both the toroidal shift and antecedent condition null models.

Figure 5 .
Figure 5.The distribution of burrowing site of subterranean rodents differs across the slopes in the study area.

Figure 6 .Figure 7 .
Figure 6.Map of the relative location of burrowing sites to shrubs with a relative origin (0, 0) for each shrub

Figure 5 .
Figure 5.The distribution of burrowing site of subterranean rodents differs across the slopes in the study area.

Figure 5 .
Figure 5.The distribution of burrowing site of subterranean rodents differs across the slopes in the study area.

Figure 6 .Figure 7 .
Figure 6.Map of the relative location of burrowing sites to shrubs with a relative origin (0, 0) for each shrub

Figure 6 .
Figure 6.Map of the relative location of burrowing sites to shrubs with a relative origin (0, 0) for each shrub.

Figure 6 .Figure 7 .Figure 7 .
Figure 6.Map of the relative location of burrowing sites to shrubs with a relative origin (0, 0) for each shrub