The study site is R. soongorica shrubland that is located 5 km apart from Urat Desert Grassland Research Station of Chinese Academy of Sciences (106°58′ E and 41°25′ N, 1650 m elevation) in western Inner Mongolia, China (Figure 1). This area has been exclosured since 2010, covering an area of about 333.34 hm². The main plant species includes Stipa klemenzii, Allium polyrhizum, Salsola collina, and Peganum harmala, as well as afforested Caragana korshinskii (Zhang, et al., 2019). The main soil types are brown-calcium soil and gray-brown desert soil without litter layer (Shaokun et al., 2016; Zhang, et al., 2019).

This area is a typical continental monsoon climate, and is affected by the continental air mass of Mongolian high pressure (Zhang et al., 2019). The mean annual temperature of 5.3°C, with the highest temperature of 37°C, and the lowest temperature of −34.4°C. The annual accumulated temperature at 10°C is 2,000–3,000°C. The mean annual precipitation is 180 mm, with a large annual variation rate. The precipitation is mainly concentrated in summer (July-August), accounting for 70% of the annual precipitation. Annual evaporation is ten times more than the annual precipitation. The frost-free period is about 130 days a year. The mean annual wind speed is 5 m s⁻¹, and the number of days with wind speed greater than 10 m s⁻¹ per year is as high as 52 days.

In May 2019, five repeated plots were set up in the R. soongorica shrubland. The distance between the edge of each plot was 15–20 m, with the area of 50 × 50 m² each. Within each plot, there were 12 sampling points set up with four sampling points beneath the shrub canopy of respect R. soongorica and C. konshinskii, and with four sampling points in the open areas between shrubs. The distance between each sampling points beneath the shrub canopy of R. soongorica and C. konshinskii and in the open areas was 5–10 m. Totally, 60 sampling points were obtained, with (4 R. soongorica + 4 C. konshinskii + 4 open areas) × 5 replicate plots.

At each sampling point, four soil cores from four directions (east, south, west, and north) were collected with a 100-cm³ metal cylinders (soil depth: 0–15 cm), and then mixed. All the composite soil samples were put in plastic sealed bags, and brought back to the laboratory. The soil sampling was carried on in spring (May), summer (August), and autumn (October) in 2019, when the annual precipitation and average temperature of was 190.8 mm and 6.46°C, respectively. The monthly precipitation and temperature in 2019 see the Figure 2 (Data sources: http://data.cma.cn/).

The soils were sampled with the exception of winter due to few soil microarthropods obtained since early November (the winter started with the temperature below 0°C).

One-third composite soil samples were used for soil microarthropod extraction by using Berlese-Tullgren funnels (Macdonald et al., 2009) in the laboratory at Ningxia University. The microarthropod extraction lasted for 48 h at room temperature using one 25-W bulb per funnel. Animals were collected and preserved individually in 75% ethanol (Macdonald et al., 2009). Extracted fauna were sorted and categorized into family level based on literatures including Pictorial Keys to Soil Animals of China (Yin et al., 2010), Acarology (Li, 1989), and Larvae Taxonomy (Ma, 1990). All arthropods were identified under the microscope (Nikon E200MV) and stereomicroscope (Olympus SZX 163).

The remaining 2/3 composite soil samples (i.e.,) were used for soil physical-chemical analysis. Some soil samples were used for the determination of soil moisture content (0–10 cm; SM, %) by oven drying (Institute of Soil Science, Chinese Academy of Sciences (ISSCAS, 1978). Other soil samples were passed through a 2 mm soil sieve to remove plant and others debris, and were then air-dried. Soil pH and electrical conductivity (EC, μS·cm⁻¹) were measured in a 1: 5 (v/v) soil-water extract, which was measured by PHS-3C acidity meter and portable electrical conductivity meter (DDSJ-308F), respectively.

Composed soil particles was measured by a Mastersizer3000 laser diffraction particle size analyzer. According to USAD classification, it can be divided into coarse sand (CS, 250–1000 μM), fine sand (FS, 50–250 μM), silt (S, 2–50 μM) and clay (C, < 2 μM) (Su et al., 2004). At each sampling points, soil temperature (ST, °C) was measured for 2 weeks by a portable thermometer with conductive wires (Sato KeiryokiMfg Co., Ltd., Japan); the averaged values by the 12 sampling points and 2 weeks represented soil temperature of each sampling seasons. The soil properties at each sampling point see the Table 1.

The number of microarthropods for each taxon groups collected from 12 sampling points within each plot was combined, in order to ensure the representativeness of microarthropod samples covering R. soongorica shrubland. We determined the abundance of each taxonomic group (individuals•200g). According to the proportion of each taxa occupying the total number of individuals, the taxonomical groups can be classified into dominant groups (accounting for more than 10% of the total individual numbers), common groups (1%–10%), and rare groups (< 1%) (Macdonald et al., 2009).

Seasonal diversity variation was examined based on Hill numbers (q = 0, 1, 2) (Chang et al., 2019). Calculated Hill number unifies the reported biological diversity parameter such that the Hill numbers 0, near 1, and 2 provide the Chao2 estimated species richness, exponent of Shannon index, and inverse of Simpson index, respectively (Chao et al., 2014; Osbrink et al., 2018). D q = ( ∑ i = 1 S P i q ) 1 / ( 1 − q ) (1) D 1 = lim q → 1 q D = e x p ( − ∑ i = 1 S P i log ⁡ P i ) (2) Where, i in the formula represents the ith species, and S is the total number of species. Pi is the percentage of the total number of individuals belonging to the group category i. When q = 0, ⁰D refers to group richness, which does not consider the relative abundance of species. When q = 1, Eq. 1 is meaningless, but when q approaches 1, the limit of this formula is the index of Shannon index (Eq. 2), which represents Shannon diversity. q = 2, indicating inverse of Simpson index, the role of other species except dominant species is underestimated, so it can be interpreted as the effective number of dominant species in the community (Chao et al., 2014; Osbrink et al., 2018).

All statistical analyses were carried out using SPSS 20.0 for Windows. Before applying parametric tests, we tested for normality and homogeneity of variances (Liu et al., 2019). One-way ANOVA and multiple comparison were used to compare the differences between different data groups. Figures were created by using Origin 2018 (Origin Lab, Northampton, United States). For all tests, statistically significant differences were assigned at p < 0.05.

We identified indicator taxon for each season using the labdsv package in R (R Development Core Team, 2009) using the duleg command (numitr = 1000, all season combinations considered). According to the results of the Indval () function, the taxa with the indicator value greater than 0.7 could be considered as the indicator taxon (Nakamura et al., 2007).

Redundancy analysis (RDA) by using Canoco Software 4.5 (Microcomputer Power, Ithaca, United States) was chosen to examine the environmental factors affecting the structure of the soil microarthropod community (Macdonald et al., 2009). The data firstly analyzed detrended correspondence analysis (DCA), which suggested RDA as an appropriate approach for further analysis (length of gradient <4 for soil microarthropod communities). RDA correlated each soil microarthropod group to the environmental factors by selecting linear combinations of environmental factors that gave the smallest residual sum of squares.

In spring, there were five dominant groups including Pygmephoridae, Bdellidae, Oppiidae, Onychiuriclae and Tomoceridae, according for 94.12% of total individuals. From spring toward summer, the number of dominant groups decreased to two groups: Pygmephoridae and Isotomidae, with the reduction of dominance to 60.56%. After summer, the composition of dominant groups shifted strongly alternative with the new groups including Stigmaeidae, Tydeidae and Oppiidae, according for 36.67%, 40.00% and 13.33% of total individuals, respectively.

We identified indicator taxon for each season (Table 2): there were two indicator taxon including Cryptopalmites (p = 0.0368) and Isotomidae (p = 0.0317) in summer; there was one indicator taxon including Stigmaeidae (p = 0.0350) in autumn. However, there was no indicator taxon observed in spring.

The total abundance and Chao2 estimated species richness were found to be significantly (p < 0.05) greater in summer than in other two seasons (Figure 3). However, there were no significant differences in exponent of Shannon index and inverse of Simpson index of soil microarthropod between the three seasons.

From the results of redundancy analysis (RDA) (Table 3), the eight variables (soil moisture, pH, electrical conductivity, soil temperature, soil clay, silt, fine sand and coarse sand) together explained 57.4% of the total variation in the data, with axes 1 and 2 explained 43.1% and 14.3% of the total variation, respectively. The species-environment relationship for axes 1 and 2 accounted for 84% of the total variance, indicating these axes accounted for the bulk of the variance in the soil microarthropod group data that could be attributed to environmental factors. Species-environment correlations for these axes were > 0.94, indicating that the soil microarthropod group data were strongly correlated with the environment factors. Monte Carlo significance tests revealed that both the first axis (p = 0.02) and all axes (p = 0.012) combined together explained a significant amount of the variation within the data.

From Figure 4, all the three seasons were sorted into two groups: (1) the summer season was on the right; (2) the spring and autumn seasons were on the left. The second groups could be sorted into two subgroups: (1) the spring was on the up; (2) the autumn was on the bottom.

Figure 4 also demonstrates that different groups within the microarthropod community responded differently to the environmental factors (Table 4). The microarthropod groups whose arrows pointed in approximately the same direction in the ordination as the arrows for the environmental factors were strongly correlated with those factors. For example, the Teneriffiidae were positively correlated with soil moisture content (p < 0.05). The Macrochelidae and Isotomidae were positively correlated with soil moisture content whereas negatively correlated with soil clay (p < 0.05). The Tydeidae was negatively correlated with soil fine sand and soil coarse sand whereas negatively correlated with soil silt (p < 0.05). The Tomoceridae was negatively correlated with soil moisture content (p < 0.05).

From the results of partial RDA (Table 5), soil temperature (F = 6.53, p = 0.002) and soil moisture content (F = 4.36, p = 0.002) were significant in the Monte Carlo permutation test, whereas those of the remaining variables were not significant under the same test (in all case p > 0.05). Soil temperature and soil moisture content appeared to be the most important environment factors driving the distribution of soil microhabitats across seasons.