Knowledge on the fine-scale spatial structure is very important for understanding plant community dynamics, as interactions among sessile plants mainly occur in the neighborhood. However, the spatial assemblage mechanism and interrelationships with productivity remain poorly understood. Here, species segregation index (Ms_com) was used to quantify the interspecific segregation (or intraspecific aggregation) of the neighborhood in plant communities from three contrasting plateaus. We employed the coordinates for 125,726 individual plants (within 1 m² plots) sampled in natural grassland communities across Loess Plateau (LP), Mongolian Plateau (MP), and Tibetan Plateau (TP) to elucidate the assembly mechanisms of Ms_com and interrelationships with above-ground net primary productivity (ANPP). We found that these grassland communities tended to assemble in intraspecific aggregations; the phenomena were more apparent in desert grassland communities. Apart from the controls of self-organization, climatic heterogeneity also plays an important role in determining Ms_com. In grassland communities on LP and MP, Ms_com was predominantly affected by reproductive allocation (RA); whereas on TP, Ms_com was more responsive and sensitive to climatic heterogeneity (especially precipitation and low temperatures). We also found that Ms_com better explained the overall variations in ANPP than species diversity, and mediated the response of ANPP to environmental variations; Ms_com and environmental effects jointly explained 22.7, 19.1, and 49.8% of the variations in ANPP on LP, MP, and TP, respectively. Our findings provide new insights into the underlying mechanisms of spatial assemblage of grassland communities, and emphasize the critical role of spatial structure in predicting the impacts of environmental change.

 Data collection

Transects setup and field survey

We tested our hypothesis by the experimental design of multiple transects comparison. We set up three transects with reference to the International Geosphere-Biosphere Programme (IGBP) single-transect investigation (Koch et al. 1995). Each transect spanned 10 sites and more than 1000 km along the precipitation gradient, including three transects for meadow grassland, four for typical grassland, and three for desert grassland (Figure S1). At each site, eight 1 × 1 m² plots were established at intervals of 50 m for a total of 240 plots. Field sampling was conducted from mid-July to late August 2018.

A 1 m² sample box that was equally divided into 100 grid cells size of 10 × 10 cm² was used for the survey. We considered the lower left corner of each cell as the coordinate origin, and identified the vascular plants and their respective coordinates. The number of stems or tillers associated with an individual was also recorded. For those plants that were larger than a single grid cell and for which the root and stem positions could not be determined, such as Androsace tapete Maxim., the center of the grid was defined as its coordinates.

In addition, for plants with both vegetative and reproductive organs, their height was measured separately, and the coverage of each species within the plots was recorded. We also harvested all plants at ground level in all plots. The samples were oven-dried to constant weight at 65 °C in laboratory to calculate the above-ground biomass (AGB) (g/m²). AGB was considered as above-ground net primary productivity (ANPP). For July/August, this AGB approximates the ANPP in grasslands of the Northern Hemisphere (Ma et al. 2010; Guo et al. 2012). The important values (IV) used to calculate Shannon and Simpson indices were calculated using the following formula: IV = (relative abundance + relative frequency + relative coverage) / 3 (Fang et al. 2009).

Soil sampling and measurements

We randomly collected one 10-cm-deep soil sample from each 1 m² plot, and collected a total of 240 samples. The soil samples were air dried and then sieved (< 2 mm) to homogenize the soil samples. We measured soil organic carbon (SOC), carbon to nitrogen ratio (C:N), and soil bulk density (BD) and pH, which have been widely used in previous studies and are easily interpretable in a biological context (Chen et al. 2018). Soil organic carbon (SOC) was analyzed using the H₂SO₄–K₂Cr₂O₇ oxidation method. Soil pH was measured by a pH electrode (Leici). Soil bulk density was determined using the volumetric ring method. The C and N content in soil samples was measured by an elemental analyzer (Vario MAX CN, Elementar, Germany).

Climatic data

The MAT and growing season temperature (T_season) (from April to August) were used to represent heat conditions. MAP and growing season precipitation (P_season) were used to represent water conditions. Because the communities of plateau ecosystem are also the result of adaptation to the characteristically stressful climates of the biomes (Ding et al. 2020), we also considered minimum temperature in January (T_min) and ultraviolet radiation (UR).

MAT, T_season, T_min, MAP, and P_season were downloaded from the WorldClim database (www.worldclim.org) at a spatial resolution of 30 arc seconds. The annual average UR for the years 2000–2014 was obtained from the Science Data Bank (http://www.dx.https://doi.org/10.11922/sciencedb.332).

Calculation of spatial structure parameters

We determined the location coordinates of 125,726 plant individuals in all plots and determined the number and identity of neighbors of each plant individuals by constructing a Voronoi diagram based on the neighborhood analysis tool of ArcGIS (Liu et al. 2017a) (Figure 1b). Species segregation index (Ms_com) (Gadow 1993; Hui & Gadow 2003) was selected to quantify the spatial structure. The value range of Ms_com is 0–1, and higher Ms_com indicated a higher interspecific segregation (or smaller intraspecific aggregation). Here, Ms_com was calculated as follows:

(1) Building structural units and removing edge effects

Each plant individual (central plant i) and its neighboring plants (j) constitute a structural unit. The number of polygon sides was the number of neighboring plants (n). The central plants (i) of adjacent polygons are neighbors (j) to each other. The neighbors of plants on the edge of the quadrat may be outside the plots, therefore we performed a de-edge effect on the constructed Voronoi diagram. Inside the blue wire frame were the spatial structural units after removing the edge effect (Liu et al. 2017a) (Figure 1b).

(2) Calculating Ms_com

Here, Msi represents the degree of species segregation of a structural unit. The minimum value was 0, indicating that the plants in the spatial structural unit were of the same species, the maximum was 1, indicating that the plants in the structural unit were of all different species. The Msi was calculated using Eq. 1. In Figure 1a,b, only the neighbor is 4 (n = 4) as an example. Further, the community-scale Mscom was calculated using Eq. 2: where, si is the number of species contained in each structural unit; n is the number of neighboring plants (n = 0, 1, 2, 3...); N is the number of structural units in each quadrat; i is the central plant, and j is the adjacent plant. When i and j are the same species, vij = 0. When i and j are different species, vij = 1.