**1. Introduction**

Shrub encroachment in grasslands and the densification of woody plant cover in savannas have been widely documented across many arid and semiarid areas of the world [1–4], including in South America, Australia, and the warm deserts of the Southwestern United States. However, until now, there has been a clear division of opinion about its ecological indication. One view is that shrub encroachment is an indicator of land degradation that is often associated with ecosystem degradation [5], such as declines in forage productivity, biodiversity, and socioeconomic potential, as well as increased erosion [6]. However, an alternative viewpoint proposing that shrub emergence is a sign of the restoration of degraded ecosystems has emerged recently [7,8], considered to support the biodiversity and a variety of ecosystem services [9]. In addition, shrub encroachment has been described as an alternative stable state occurring several times during the last two millennia in African savannas [10]. Some studies also found that the role of shrub species was important in the rehabilitation of degraded sandy land ecosystems [11]. Furthermore, di fferent opinions exist as to the causes of shrub encroachment in di fferent types of ecosystems [12]: some studies found that underuse leads to shrub and subsequent tree encroachment and, finally, conversion to forest [12]; however, overgrazing usually causes shrub encroachment in Africa and the Mongolia Plateau [10]. This di fference in the perception of shrub encroachment could lead to completely di fferent judgments of the states and transitions of shrub-encroached ecosystems, which would further a ffect decisions about their conservation and managemen<sup>t</sup> [13]. An innovative study [14], combined with the landscape history and other environmental factors, quantified the process and indicative significance of shrub encroachment and provided us with new ideas. Therefore, more attention needs to be paid to ecosystem changes after shrub encroachment to clarify their ecological significance, especially the significance of the degree of shrub encroachment and changes in distribution patterns and structure. This information could then be used to quantify the ecological indications of shrub encroachment.

As a common tree species distributed widely throughout the forest-steppe ecotone on the Mongolian Plateau [15,16], *Ulmus pumila* forms a stable savanna-like woody-herbaceous complex ecosystem in association with grass in the Horqin Sandy Land, the Otindag Sandy Land, and the Hulunbuir Sandy Land in Northern China [17,18]. The sparse *U. pumila* trees on the temperate savanna-like ecosystem have ecological significance because they provide sand stabilization and also resting places for livestock [19]. In recent years, due to various changing conditions, such as climate change and land use change in the course of economic development [20–22], there has been severe damage to the *U. pumila*-dominated temperate savanna-like ecosystem. This damage has led to an increase in the number of shrubs and the densification of woody plant cover as well as outcomes such as decreasing temperate savanna-like ecosystem area, the loss of structural integrity, poor population regeneration, and changing spatial patterns [23,24]. Previous studies have focused mainly on the vegetation distribution, possible encroachment mechanisms, vegetation characteristics, species composition, and carbon budget in the *U. pumila*-dominated temperate savanna-like ecosystem after shrub encroachment [25,26]. It was found that shrub coverage in the grasslands of Inner Mongolia reached 12.8% and that the emergence of *Caragana microphylla* was an indicator of grassland degradation in Inner Mongolia [27]. However, the spatial patterns of woody plants and their interrelationships associated with shrub encroachment are often overlooked despite having important e ffects on ecosystem functions [28].

Spatial pattern analysis is an important method for studying the distribution and relationships of di fferent plants [29,30]. The spatial patterns of species and the spatial relationships among species significantly impact growth, reproduction, death, resource utilization, and gap formation among species [31,32]. The analysis of a species' spatial pattern helps us to understand both the ecological processes that form the pattern (such as intra- and interspecific competition, interference, and environmental heterogeneity) and the ecophysiological traits of the plant species, including the relationships between these plants and their environment [33,34]. Recently, a spatial pattern analysis method has been used to clarify the vegetation degradation processes underlying the patterns of individual species in semiarid and arid areas [35,36].

In this study, our objective is to explore how shrub encroachment a ffects the spatial distribution and interaction of woody plants in a temperate savanna-like ecosystem. We hypothesize that there are

di fferences in the relationships between *U. pumila* trees of di fferent DBHs and shrubs. The following questions will be addressed: (1) What are the distribution patterns of *U. pumila* trees and shrubs in a temperate savanna-like ecosystem? (2) What are the interspecific interactions between *U. pumila* trees and shrubs? and (3) Do tree–shrub spatial associations change as a function of tree size?

#### **2. Materials and Methods**

#### *2.1. Study Area*

*U. pumila*-dominated temperate savanna-like ecosystem, which is distributed widely throughout the forest-steppe ecotone on the Mongolian Plateau, China, is quite di fferent from North American prairies and African savannas in terms of climate, soils, and dominant plant functional types; however, there are some taxonomical similarities among these ecosystems at the genus and family levels [37,38]. The study area experiences a typical temperate continental semiarid climate, with warm summers and cold winters. The mean annual precipitation is approximately 350 mm, 70% of which occurs from June to August. The mean annual pan evaporation is 1900 mm. The mean annual air temperature is 1.6 ◦C, with a maximum monthly temperature of 17.8 ◦C (July) and a minimum monthly temperature of −17.6 ◦C (January) [37,38]. The main soil types are aeolian sandy soils with a mean thickness of 200 cm and a calcic horizon occurring at 30–100 cm depth. In this forest-steppe ecotone, *U. pumila* is dominant; the other species that occur include other woody species, including *C. microphylla*, *Spiraea aquilegifolia*, *Ribes diacanthum* and *Salix linearistipularis*, and herbaceous plants, including *Leymus chinense* and *Cleistogenes squarrosa* [39]. The main land use in the area is grazing, and the livestock are primarily cattle. The local herders have implemented alternating seasonal grazing (winter and spring) and mowing (autumn) regimes since 2000 [40]; before that time, annual grazing was the local practice. Shrub encroachment, such as the spread of *C. microphylla*, has been clearly observed in grasslands and savannas in this area, resulting in a landscape characterized by a mosaic of shrub and grass patches.

#### *2.2. Research Site and Data Collection*

Most work about the spatial distribution and interaction of woody plants had been done in single hectare or smaller plots, but the relative rarity of many species in forests necessitated large-scale census plots. Thus, plots usually with more than 5 hectares, named large plots, are considered as representative of local vegetation, which could cover the local typical vegetation and topography in a region. By a large plot, the interference of scale and environmental heterogeneity could be avoided for the study on vegetation composition, pattern, and biodiversity [41]. This method has been widely used in global forest dynamic monitoring.

In this study, this large plot method was applied and one large plot with an area of 25 ha (500 m × 500 m) was established following the plot standards of the Center for Tropical Forest Science (CTFS) network [41]. It was located in the Otindag Sandy Land in Inner Mongolia, China (115◦16 E, 42◦50 N), in a typical area of *Ulmus pumila*-dominated temperate savanna-like ecosystem (Figure 1).

Then, this large plot was further divided into 625 subplots (20 m × 20 m). All free-standing trees with stem diameters at breast height (DBHs) of more than one centimeter and all shrubs were tagged, mapped, and identified to species during the summer of 2013–2014. The coordinates of woody plants in each subplot were recorded using an Electronic Total Station, with the southwestern corner of the subplot as the origin. Additionally, in the center of each subplot, a small plot with the size of 1 × 1 m was setup. All herbaceous species were identified and their cover ratios and numbers were recorded.

All *U. pumila* trees were classified into three categories according to their DBH [42], namely, old trees (DBH ≥ 20 cm), medium trees (5 cm ≤ DBH < 20 cm), and juvenile trees (DBH < 5 cm) (Table 1).

**Figure 1.** Map of the Otindag Sandy Land and the location of the study region (**a**) location; (**b**) sample map; (**<sup>c</sup>**,**d**) pictures of the landscape.


**Table 1.** Basic parameters of *U. pumila* trees and shrubs.

#### *2.3. Data Analysis*

The pair-correlation function g (r) is a statistical method used to estimate the number of points within concentric rings at a distance r rather than within a certain radius and is especially sensitive to small-scale effects [43,44]. There are two g (r) functions, i.e., univariate and bivariate g (r). In this study, the univariate g (r) function was used to analyze the spatial distribution patterns within woody plants, and the bivariate g (r) function to quantify the both intra- and interspecific spatial association among tree and shrub plants [43,44].

Firstly, the univariate g (r) function was used to analyze three tree categories (old trees, medium trees, and juvenile trees) and all shrub species. The null model of complete spatial randomness (CSR) as a null hypothesis was used for all the univariate analyses. Secondly, for the bivariate analyses, two cases were considered. One case was that the relationship between small and large trees was considered. Since large trees may impact the distribution pattern of small trees within their area of influence (competition), a bivariate g function analysis was conducted for these two size classes using both the toroidal shift and the antecedent condition null model options [43]. This tests whether the patterns of distribution of small and large trees were generated by independent processes. The antecedent condition model tests whether one pattern (small trees) is influenced by a second

pattern (large trees), assessing whether there are more (or fewer) small trees in the neighborhood of large trees than expected under a random distribution of small trees [43]. The second case concerns the interaction between trees and shrubs. Because the spatial distributions of plants in plots seem to be affected significantly by drought stress and habitat heterogeneity (e.g., soil patch and microtopography), we examined the spatial association between the two species with the independent null model [44].

All analyses were conducted with Programita software (2016). To assess the significance of the data by comparing them with the null model, 95% confidence intervals were obtained by running 199 Monte Carlo simulations [45].
