Carabid community structure in northern China grassland ecosystems: Effects of local habitat on species richness, species composition and functional diversity

Study areas and sampling design

The study was carried out in the Ningxia region, northern China. We selected three sampling areas representing the three main ecosystems in the region: desert steppe, typical steppe and meadow steppe (Figs. 1 and 2). In each area, we identified different habitats to reflect within-ecosystem variability. We selected the study sites to be representative of the variability of environmental conditions within and between Chinese grassland ecosystems. We adopted a stratified sampling design, with a different number of trapping sites among grassland types to reflect their within-ecosystem variability. To make results comparable, we used the same number of traps (15) for each habitat within each ecosystem.

The desert steppe area (Fig. 2A) is located in eastern Ningxia, Yanchi county (37°59′13″N–107°05′42″E). This area is characterized by a cold, semi-arid continental monsoonal climate zone (Liu et al., 2015), with an average annual temperature of 8.3 °C (−8 °C in January, 22 °C in July), and average annual precipitation of around 200 mm (Kang et al., 2007). The vegetation is characterized by typical drought-tolerant plant species, such as Agropyron mongolicum, Artemisia desertorum, Artemisia blepharolepi and Stipa spp. The typical steppe area (Fig. 2B) is located in southern Ningxia, Guyuan County, near the Natural Reserve of the Yunwu Mountain. This area is characterized by a continental monsoon climate. Average annual temperature is 5.7 °C (−22 °C in January, 28 °C in August) and annual rainfall is 350 mm (Kang et al., 2007). The top of the mountain (36°12′16″N–106°24′37″E) is characterized by grass vegetation crossed by patches of cut grasses that serve as fire belts. The natural vegetation on the top of the mountain includes Stipa bungeana, S. grandis, Artemisia frigida, Thymus mongolicus and Heteropappus altaicus. The bottom of the mountain (36°15′6″N–106°23′5″E) is occupied by crop fields and natural vegetation, including S. bungeana, Artemisia frigida, T. mongolicus and Potentilla acaulis. In this area, we selected three sectors; the first and second sectors were located at the top of the mountain, in natural patches of grass vegetation and in fire belts, respectively; the third sector was selected at the bottom of the mountain. The meadow steppe area (Fig. 2C) is located in western Ningxia, Haiyuan county, near the Nanhua mountain. This area is characterized by a semi-humid climate, with an annual temperature of 7 °C (−7 °C in January, 20 °C in July), and average annual precipitation of 450 mm (Kang et al., 2007). Within this area, we selected two sampling sectors to reflect different vegetation types. The first sector (36°26′50″N–105°38′24″E) was located at the south-west side of the mountain peak (2,600 m) and is dominated by several species of the genus of Festuca, principally the alpine fescue Festuca brachyphylla. The second sector (36°25′13″N–105°36′41″E) was located at the bottom of the mountain peak (1,800 m) and its vegetation is dominated by S. bungeana, Artemisia frigida and Achnatherum splendens.

In total, we selected 15 sampling sites in the desert steppe, 45 sites (15 sites per sector) in the typical steppe and 30 sites in the meadow steppe (15 sites per sector).

Sites were separated by at least 150 m from each other. At each sampling site, five pitfall traps were placed at a distance of least five m from each other. Pitfall traps were made of plastic cups (diameter: 7.15 cm, depth: nine cm) dug into the ground and filled with 60 ml of an attractant solution (vinegar, sugar, 70% alcohol and water in the following proportion: 2:1:1:20). We used pitfall traps with a diameter slightly over than seven cm because this size meets Luff’s (1975) suggested optimal diameter (seven cm) for carabids. Pitfall traps were put down once a month in mid-month, from May to September 2017, and collected 3 days after. A period of 3 days was chosen because many traps were found completely full of beetles (especially Tenebrionidae) in 2 or 3 days. In total, we used 2,250 pitfall traps (90 sampling sites × 5 pitfall traps × 5 sampling dates). Prior to analyses, we pooled the data from the five pitfall traps of each sampling, because soil and vegetation characteristics were observed at the sampling site level. Trap content was sorted in the laboratory and carabids identified to species level and assigned to trophic categories (herbivores vs. predators). All material is preserved in the insect collections of the School of Agriculture of Ningxia University.

Vegetation and soil characteristics

At each sampling site, we set up one quadrat frame of 0.25 m² to record plant dry biomass (PB, g/m²), cover (PC, % of soil covered by plants), density (PD, number of plants per m²), height (PH, average, cm) and species diversity, expressed as richness (PSD). Near the quadrat frames we collected samples of above-ground litter to measure litter dry mass (SL, g/m²) and samples of soil (10 cm depth) to measure SM (%) and SBD (g/cm³). SM was estimated using the thermogravimetric method also known as the oven dry method (Majumdar, 2001): SM = [(W₂−W₃)/(W₃−W₁)] × 100 where W₁ is the weight of the empty aluminum box (g); W₂ is the weight of the box + soil sample (g) and W₃ is the weight of the box and oven dry soil (g). SBD was estimated using the ring knife method (Bi, Zou & Zhu, 2014), with the formula SBD = (W_r × 100)/(V_r × (100 + SM)) where W_r is the weight of the soil in the ring knife; V_r is the ring knife volume. We also measured ST (10 cm depth) using a portable multiparametric probe TRS-II (Zhejiang Tuopu Instrument Co. Ltd., Hangzhou City, China; accuracy of ± 0.5 °C). Monthly mean values of humidity (Hum), precipitation (Prec) and temperature (Temp) were recorded from the meteorological stations in each county. Range, mean and standard deviation for the aforementioned environmental variables are given in Table S1.

Species characteristics

We expressed functional diversity with reference to two aspects: dispersal power and body size. For dispersal (FD-movement), we used the following morphometric traits, under the assumption that longer and more robust legs facilitate beetle movements on the ground: (1) length (from apex to coxa) and (2) maximum width of metafemurs, (3) length of metatibiae, (4) length of metatarsi and (5) presence of wing. For body size (FD-size), we used: (1) width of the head, (2) pronotum maximum width and (3) pronotum maximum height, (4) elytral length and (5) elytral width. All these traits were used to compute a total functional diversity (FD-total). To calculate FD-total, we also considered the following traits: (1) length of antennae, (2) feeding habits (predators vs herbivores) and (3) five characteristics in the mandibles (density of ventral groove; roughness of dorsal crenulations; sharpness of incisor ridge; ratio width/length of left and right mandible). Measurements were done using a digital caliper (precision to 0.01 mm, Stainless Hardened). For species with more than 50 collected individuals, we measured 50 specimens (14 species); for other species, we measured a variable number of individuals depending on their abundance (11 species).

Data analysis

We assessed species richness using the individual-based rarefaction method implemented in the “iNEXT” library of R (Hsieh, Ma & Chao, 2016). The input matrix was based on species abundance (total number of individuals from each sample site). We rarefied data to the smallest number of collected individuals. Since the number of traps in the different sites was different, in all analyses dealing with species abundance we used species’ activity density, calculated as the number of individuals from each species divided by the number of traps used in each site.

For each trait, we calculated functional diversity indices using Rao’s quadratic entropy, which expresses the sum of the dissimilarities in the trait space among all possible pairs of species weighted by species-relative abundances (Rao, 1982; Botta-Dukát & Wilson, 2005). High functional divergence should indicate a high degree of niche differentiation (Mason et al., 2005). To express functional diversity based on multiple traits, we calculated species dissimilarities with the commonly used Gower distance (Pavoine et al., 2009; Laliberté & Legendre, 2010). Calculations were done using the “StatMatch” package and the “melodic” function in R (De Bello et al., 2016).

Differences in species richness and functional diversity between the three grassland types were tested using a Nested analysis of variance (Nested ANOVA, with type of grasslands as fixed effect and sub-types of grassland as random effect), followed by post-hoc Tukey tests using the “multcomp” package in R (Hothorn, Bretz & Westfall, 2008). The effects of vegetation, soil and climate characteristics on beetles rarefied richness and functional diversity were investigated using a random-effect eigenvector spatial filtering (RE-ESF) approach (Murakami & Griffith, 2015) with the “spmoran” package in R (Murakami, 2018) to take into account spatial dependence.

A Moran test showed a spatial autocorrelation at regional level (Moran’s I = 0.024, P-value < 0.001) but not at the grassland scale (Desert steppe: Moran’s I = −0.008, P-value = 0.151; Typical steppe: Moran’s I = −0.004, P-value = 0.397; Meadow steppe: Moran’s I = −0.007, P-value = 0.645). However, we decided to use the RE-ESF in all models to make the results comparable.

Preliminary to RE-ESF analysis, variance inflation factors (VIF) were calculated using the “usdm” package (Naimi et al., 2014) to detect possible collinearity between explanatory variables and determine the stability of models. A high VIF (>10) indicates that the predictor is strongly dependent on others and does not carry independent information. No collinearity was found between the variables, with all being VIF <10 (Table S2).

Using the matrix of geographical coordinates, Moran’s eigenvectors and their corresponding eigenvalues were calculated using the “meigen” function in the “spmoran” package. The resulting eigenvectors are used as synthetic explanatory variables in regression analysis (Griffith & Peres-Neto, 2006).

The effects of vegetation, soil and climate characteristics on beetle community structure were investigated using Canonical correspondence analysis (CCA) with abundance data. This technique was particularly appropriate to our data because it addresses with the double-zero problem which characterizes community compositional data (Legendre & Gallagher, 2001) and does not try to display all variation in the data, but only the part that can be explained by the constraints (Oksanen et al., 2015). Permutation tests (999 permutations) were run to assess model significance. The sum of the canonical eigenvalues was used as a measure of the variability in the response variables explained by predictors. The importance of predictors was assessed beforehand by using the VIF (Table S3; Oksanen et al., 2015). Analyses were conducted in R using the “vegan” package (Dixon, 2009; Oksanen et al., 2015). We used the “step” function to determine the best model and the most important predictors in each CCA.

The variables used in the RE-ESF and CCAs were not redundant and represent different, but not mutually exclusive, hypotheses and each hypothesis has been evaluated individually with a selection procedure. We think that a further adjustment of P-values would result in a higher risk of pruning variables that are important. Thus, we did not adjust the P-values of variables selected as significant, but focused on the magnitude of the P-values and the consistency of results (see Moran, 2003).

Finally, we investigated patterns of ß-diversity, that is, species variations among habitats. We used the approach of Baselga, Jiménez-Valverde & Niccolini (2007) and Baselga (2010, 2012) for partitioning the overall ß-diversity (ßsor, Sørensen coefficient) among habitats into true species-replacement or pure turnover (ßsim, Simpson coefficient) and nestedness (ßnest = ßsor − ßsim) components. In this respect, nestedness quantified the part of compositional change caused by ordered species loss, whereas pure turnover was related to the exchange in species composition. Relationships among assemblages were investigated by cluster analysis using the UPGMA (unweighted pair-group method, arithmetic average) amalgamation rule. Calculations were done with PAST v.3 (Hammer, Harper & Ryan, 2001).

Differences in richness and functional diversity

We collected a total of 6,873 individuals belonging to 25 carabid species (Table 1). Overall, 18 species were predators and seven were herbivores (six herbivores and six predators in the desert steppe, four herbivores and 14 predators in the typical steppe, and four herbivores and 15 predators in the meadow steppe). Range, mean and standard deviation of rarefied richness, FD-total, FD-movement and FD-size are given in Table S4. The desert steppe was the grassland type with the lowest values of rarefied species richness and the three measures of functional diversity, whereas no significant differences were found between the meadow and the typical steppe (Fig. 3; Table S5).

Influence of environmental variables on richness and functional diversity

Climate predictors had strong effects on carabid species richness (Table 2). Hum was positively related to species richness in the desert and typical steppes, as well as at the regional scale. Prec had a positive effect on richness at regional scale, but a negative effect in the desert steppe. Temp had a positive effect on richness at regional scale, in the meadow and in the typical steppes. At regional scale SBD, ST and SM had negative effects. SBD and ST also had negative effects in the typical steppe, and SBD had a positive effect in desert steppe, whereas none of the soil predictors were significant for the meadow steppe. None of the vegetation predictors were important in explaining species richness in the desert or in the meadow steppe. PC and PSD had positive effects in the typical steppe and at regional scale, respectively.

At regional scale (Fig. 4; Table 3), the total variance explained by CCA was 24%, with the first two axes accounting for 57% of the explained variance. Constraints were significant (F = 9.53, P < 0.001) and all predictors were retained in the selection procedure (Table S6). However, the first axis showed that SM and ST were the most important variables, acting in opposite directions. In the desert steppe community (Fig. 5; Table 3), CCA explained 40% of total variance, with the first two axes accounting for 56% of the explained variance. Constraints were not significant (F = 1.53, P = 0.07), with only humidity and temperature being included in the final model after the selection procedure (Table S6). In the typical steppe community (Fig. 6; Table 3), CCA explained 37% of total variance, with the first two axes accounting for 68% of the explained variance. Constraints were significant (F = 9.74, P = 0.01), and the final model included PB (as the most important variable), PC, PD, SL, Hum, Prec and Temp (Table S6). In the meadow steppe community (Fig. 7; Table 3), CCA explained 25% of total variance, with the first two axes accounting for 65% of the explained variance. Constraints were significant (F = 3.38, P = 0.01), with PB, PH, SL, Hum and Temp being retained in the final model. PH and Temp were the most important variables (Table S6).

The response of functional diversity to habitat characteristics (vegetation, soil and climate) differed among grassland types (Table 4). At regional scale, five predictors showed important effects on FD-total: PB (negative, marginally non-significant), ST (negative, marginally non-significant), Hum, Prec and Temp (all positive). In the desert steppe only two predictors had important effects: PD (negative, marginally non-significant) and Hum (positive). In the typical steppe, only the climate characteristics Hum and Temp were important and both had positive effects. In the meadow steppe, PC (negatively) and ST (positively) were important predictors of meadow FD-total (temperature had a marginally non-significant positive effect).

The response of FD-movement (Table S7) at regional level was similar to that of FD-total, with the exception of Hum, which was not significant. By contrast, Hum was the only important predictor of FD-movement in the desert steppe. As at the regional scale, ST, Hum and Temp were important predictors of FD-movement in the typical steppe. In the meadow steppe, PB and PC had negative effects, whereas ST and Prec had positive effects.

At regional scale, four predictors showed important effects on FD-size (Table S8): ST (negatively) and Prec, Hum and Temp (positively). In the desert steppe, two variables influenced FD-size: Hum (positively) and PD (negatively and marginally non-significant). In the typical steppe, only Hum and Temp showed significantly positive effects on FD-size. In the meadow steppe, important predictors where PC and PH (negatively), and Prec (marginally non-significant) and Temp (both positively).

Beta diversity

Species composition varied greatly among grassland types, with seven species shared by all three types (Fig. 8). The overall ß-diversity (ßsor) pattern (Fig. 9A) indicated that the desert steppe differed the most. The meadow sector occupied by Festuca spp. also clustered apart, whereas the other meadow sector, occupied by Stipa spp., clustered with the typical steppe. When the pure turnover (ßsim) is considered (Fig. 9B), the desert was again identified as the ecosystem differing the most. All three sectors of the typical steppe clustered together and were separated by the two meadow sectors. In the analysis of the nestedness component (ßnest) the desert steppe clustered with the sector of typical steppe located at the bottom of the mountain peak and the meadow steppe dominated by Festuca spp., whereas the other meadow sector clustered with the typical steppe (Fig. 9C).