Introduction
Grasshopper plagues, one of the main disasters in the Inner Mongolia grassland, China, have devastating consequences for grassland ecosystems, including desertification and degradation. It is difficult to restore extremely destroyed grasslands within 20 years. Grasshopper plagues seriously affect grazing production and the lives of local people, as well as posing a serious threat to the nation's ecological engineering projects, such as the sandstorm source control project around Beijing and Tianjin that started in 2000 and the natural grassland vegetation restoration and conservation project in Inner Mongolia that started in 2001. Therefore, it is necessary and important to precisely predict the location and scope of grasshopper infestations and to develop effective preventive and control measures.
Spatial heterogeneity of grasshopper occurrence mainly depends on habitat conditions if weather conditions are similar; thus, habitat suitability is an important indicator of grasshopper occurrence. Relevant habitat conditions include soil properties (such as soil type, texture, humidity, temperature, pH, salinity, and inorganic matter content), vegetation characteristics (such as plant community composition, plant abundance, and plant nutrient content), and topographic elements. These factors affect grasshoppers’ selection of oviposition sites, incubation and mortality rates of grasshopper eggs, development of nymphs, and reproduction of adults (Liu et al., Reference Liu, Xi and Chen1984; Kang et al., Reference Kang, Li and Chen1989; Ni et al., Reference Ni, Jiang, Wang, Gong, Wang and Voss2000; Torrusio et al., Reference Torrusio, Cigliano and de Wysiecki2002; Wang & Ni, Reference Wang and Ni2003; Chen, Reference Chen2007; Sirin et al., Reference Sirin, Eren and Çıplak2010; Ebeling et al., Reference Ebeling, Allan, Heimann, Köhler, Scherer-Lorenzen, Vogel, Weigelt and Weisser2013). In addition, human activities, such as heavy livestock grazing, afforestation, and burning plants may strongly affect grasshopper occurrence by changing the habitats (Kang, Reference Kang1997; O'Neill et al., Reference O'Neill, Olson, Rolston, Wallander, Larson and Seibert2003; Bazelet & Samways, Reference Bazelet and Samways2011a , Reference Bazelet and Samways b ; Cease et al., Reference Cease, Elser, Ford, Hao, Kang and Harrison2012). Therefore, in order to monitor and predict grasshopper occurrences, we must explore the relationships between habitat conditions and grasshopper abundance and develop quantitative methods to evaluate habitat suitability.
A number of qualitative studies have explored these relationships and grasshopper habitat suitability since the 1970s. In quantitative studies, statistical methods are used to determine whether a habitat factor or interactions among factors contribute significantly to grasshopper abundance, such as correlation analysis, significance test, and multivariate ordination techniques (Torrusio et al., Reference Torrusio, Cigliano and de Wysiecki2002; Sirin et al., Reference Sirin, Eren and Çıplak2010; Bazelet & Samways, Reference Bazelet and Samways2011a , Reference Bazelet and Samways b ; Cease et al., Reference Cease, Elser, Ford, Hao, Kang and Harrison2012; Ebeling et al., Reference Ebeling, Allan, Heimann, Köhler, Scherer-Lorenzen, Vogel, Weigelt and Weisser2013; Crous et al., Reference Crous, Samways, Pryke, Stewart and Bezemer2014). Zhou et al. (Reference Zhou, Wang, Zhao and Zhang2012) used geostatistics and Kriging interpolation to obtain the spatial pattern of grasshoppers and vegetation characteristics and their spatial correlations. However, statistical methods themselves cannot explicitly determine how habitat factors contribute to grasshoppers distribution, and thus, cannot be used to predict grasshopper infestations. Modelling plays an indispensable role in this aspect. Statistical models, such as the generalized linear model and multivariate binary logistic regression models, are often used to detect grasshoppers’ responses to environmental variables and to predict the distribution and abundance of grasshoppers (Bazelet & Samways, Reference Bazelet and Samways2011a ; Buse & Griebeler, Reference Buse and Griebeler2011; Badenhausser & Cordeau, Reference Badenhausser and Cordeau2012). Hernández-Zul et al. (Reference Hernández-Zul, Quijano-Carranza, Yáñez-López, Ocampo-Velázquez, Torres-Pacheco, Guevara-González and Castro-Ramírez2013) built a population dynamic model of grasshopper growth and development as a function of environmental conditions (including temperature, rainfall, grass biomass, and physical soil properties). Ni (Reference Ni2002) built a fuzzy evaluation model to explore the relationships among topography, soil, and vegetation factors and grasshopper occurrence, and obtained the spatial distribution of grasshoppers in Qinghai Lake region, China.
The Inner Mongolian grasslands are the largest in China and represent typical Eurasian semiarid ecosystems. Nevertheless, to the best of our knowledge and with the exception of our previous study (Zhang et al., Reference Zhang, Zhang and Chen2012), only a few studies have used models to examine the distribution of grasshoppers in Inner Mongolia grasslands, although many qualitative studies have been conducted since the 1980s. In our previous study, we surveyed grasshopper density and topography, soil, vegetation, and land use factors, which might closely relate to grasshopper occurrence in Xianghuangqi County, southwest of Xilingol League in Inner Mongolia. Using these field data, we built a fuzzy evaluation model combined with 3S (GIS, RS, and GPS) technology to estimate the habitat suitability for the potential occurrence of grasshoppers (POG). The results showed that the spatial heterogeneity of POG was mainly related to habitat factors, particularly elevation and soil sand content. Some limitations were present in the study: first, it was assumed that there was a linear relationship between each habitat factor and grasshopper density; secondly, the effect of slope on grasshopper density was eliminated using multivariate analysis of variance; thirdly, only the independent effects of the individual factors were considered in the model, and the two-factor interactions were simplified as a random error.
The present study aimed to overcome these drawbacks to improve the model's ability to estimate habitat suitability and indicate the POG. We used the measured data and the data processing software Excel–GeogDetector to determine the degree to which the spatial pattern of grasshopper density is consistent with that of a habitat factor (Wang et al., Reference Wang, Liao and Liu2010; Wang & Hu, Reference Wang and Hu2012). Excel–GeogDetector is an Excel version of a spatial geographical detector (http://www.sssampling.org/Excel-geodetector/), initially used to analyse risk factors for epidemic diseases and the potential impacts of the interactions of different risk factors on health. It is composed of four detectors: factor detector, interaction detector, risk detector, and ecological detector. The risk detector indicates potential risk areas, and the ecological detector identifies the impact differences of two risk factors. The details on factor detector and interaction detector will be given later. The basic idea of GeogDetector is spatial variance analysis. It can explore the spatial relationship between epidemic diseases and risk factors, while the traditional ordination method of Principal Components Analysis can only reduce the dimension of data to obtain a linear combination of a small number of variables. It also provides a new approach to revealing the hidden factor interactions, with no need for any assumptions or restrictions (Hu et al., Reference Hu, Wang, Li, Ren and Zhu2011). We determined the memberships of grasshopper density ranks for each habitat category and determined the weights of each habitat category by using an analytic hierarchy process. These results were used to construct a model evaluating grasshopper habitat suitability. Finally, the modelled results were validated by comparing them with measured data.
Materials and methods
Study area and organism
The study area is located in Xianghuangqi County, Inner Mongolia (41°56′–42°45′N, 113°32′–114°45′E) (fig. 1). It covers an area of 5.1 × 105 hm2, of which 97.8% is occupied by steppes, with average altitude of 1300 m. The study area experiences a continental middle temperate monsoon semiarid climate. The area receives 3031.6 h of annual mean sunshine, has a mean annual average temperature of 3.1°C, and has mean annual precipitation of 267.9 mm (65.2% of which falls between June and August). The zonal native vegetation cover is temperate bunchgrass steppe, dominated by Stipa krylovii, Cleistogenes squarrosa, and Artemisia frigida. The other vegetation types present include temperate short bunchgrass and short semi-shrub desert steppe, and temperate grass and forb meadow steppe. The zonal soil is chestnut soil, and the intrazonal soil types include aeolian sandy soil, saline meadow soil, saline alkali soil, and lithosol (Yang et al., Reference Yang, Wu, Yao, Wang, Li and Shi2009; Zhang et al., Reference Zhang, Zhang and Chen2012). The great variety of habitats was the main reason for choosing this study area.
Fig. 1. Location of study area and field plots sampled in 2010 and 2011 in Xianghuangqi County, Inner Mongolia, China.
The main grasshoppers occurring in the study area include Oedaleus decorus asiaticus (Bei-Bienko, 1961), Dasyhippus barbipes, Bryodema luctuosum, and Myrmeleotettix palpalis. The present study considered only the occurrence of O. d. asiaticus, the primary pest that influenced the area from 2001 to 2010. O. d. asiaticus is univoltine: the adult density peak occurs around mid-to-late-July, adults oviposit around late July or early August, and the egg pods overwinter until nymphs hatch around early-to-mid-June in the next year. O. d. asiaticus incubates its eggs only at the sites where the adults have oviposited in the previous year and does not disperse until the nymphs become adults after mid-July. O. d. asiaticus is an oligophagous insect, mainly feeding on grass (Chen, Reference Chen2007; Zhang et al., Reference Zhang, Zhang and Chen2012).
Field survey and data sources
Field surveys were conducted every 3–5 km along both sides of the roads and pastoral grazing paths over the entire study area. The 234 and 342 field plots were sampled in early-to-mid-July of 2010 and 2011, respectively (fig. 1), before grasshoppers started to disperse. The 335 and 369 field plots were sampled in late-June-to-early-July of 2012 and 2013, respectively, and they were surveyed again in late-July-to-early-August after grasshoppers might have dispersed.
A plot was a patch with similar habitat conditions. Two or three observers walked side by side from the start locations to the end locations within one plot, in the same direction, to visually count the amount of O. d. asiaticus per m2. In most cases, each observer covered more than 10 m wide and two or three observers could cover the whole area of the plot. We summarized the data recorded by all the observers and estimated the mean density for the plot.
At each plot, some habitat conditions related to grasshoppers were measured. At the centre of each plot, we measured the latitude, longitude, elevation, aspect, and slope. Plant species richness and abundance, and plant phenology and vitality records were qualitative descriptions. Plant species richness was estimated using a direct counting method, but rare species were excluded. Plant species abundance was estimated based on the amount and coverage of plants, recorded in Drude scales: Soc. (sociales), Cop.3 (copiosae3), Cop.2 (copiosae2), Cop.1 (copiosae1), Sp. (sparsae), Sol. (solitariae), and Un. (unicum) (Jaroshenko, Reference Jaroshenko1961). Plant phenology mainly included vegetative period, flowering period, heading period, and fruiting period. Plant vitality was simply classified as strong, medium, and weak growth vigour. Vegetation coverage was visually estimated within three 1 m × 1 m quadrats, where the soil and vegetation conditions indicated the plot. Soil volumetric water content (SVWC) in the top 5 cm of the soil was measured at these quadrats using a soil moisture meter (Aquaterr M–300). Other auxiliary information was also surveyed, such as structure and distribution of gravel and sand, intensity of livestock grazing, and afforestation. We investigated the grassland area and the amount of livestock that each herding family possessed, and then calculated the amount of livestock per unit area to indicate the intensity of livestock grazing. Only a few plots were afforested, and we recorded the afforested sites and tree abundance. The information about preventive or control measures and grasshopper infestations in the previous years were obtained from the local herders and staff at the grassland station.
The 49 and 197 soil samples were taken from the top 5 cm of the soil in 2010 and 2011, respectively, and soil sand contents were measured using the hydrometer method in the laboratory. For the unmeasured plots, soil sand content values were derived from the spatially interpolated data by using the Kriging method in ArcGIS software (spatial resolution of 30 m × 30 m). For each plot, soil type data came from the digital map (cartographic scale of 1:1,000,000) provided by Data Centre for Resources and Environmental Sciences of Chinese Academy of Sciences (http://www.resdc.cn/). In addition, the aspect and slope values for each plot were derived from the digital elevation map (spatial resolution of 30 m × 30 m), after the measured values were found to contain errors; this map was provided by International Scientific Data Service Platform (http://datamirror.csdb.cn/).
Classification of grasshopper density and habitat factors
We classified grasshopper density into four groups: very high (>20 ind. m−2), high (10–20 ind. m−2), low (5–10 ind. m−2), and very low (≤5 ind. m−2).
To derive definitive relationships between grasshopper occurrence and habitat factors, we mainly focused on the relatively stable habitat factors, which were thought to vary little with interannual and seasonal climate change. Given that SVWC and vegetation coverage and their relationships with grasshopper occurrence may have great variations between years and seasons, they were not considered for constructing the model system.
Regarding the continuous habitat factors, elevation was initially divided every 50 m; slope every 2°; soil sand content every 10% (fig. 2). Some subcategories were combined into the same category if the corresponding grasshopper densities were similar or exhibited similar variations (table 1).
Fig. 2. Relationship between grasshopper density and habitat factors in Xianghuangqi County. Aspect: −1, Flat terrain; 1, Sunny slope; 2, Half-sunny slope; 3, Half-shady slope; 4, Shady slope. Vegetation types: 1, Stipa krylovii and Cleistogenes squarrosa steppe; 2, Stipa grandis and Artemisia frigida steppe; 3, S. grandis, Leymus chinensis, and Salsola collina steppe; 4, Allium mongolicum and Carex tristachya steppe; 5, Caragana microphylla Lam. steppe. Soil types: 1, Chestnut soil; 2, Dark chestnut soil; 3, Light chestnut soil; 4, Salinized chestnut soil; 5, Basified chestnut soil. Land cover types: 1, High-coverage grassland; 2, Medium-coverage grassland; 3, Low-coverage grassland; 4, Swampland; 5, Barren rocky land; 6, Bottomland; 7, Dry cropland. Outliers and extreme values denote those that are 1.5 and 3.0 times larger than the interquartile range, respectively. The calculations of the mean grasshopper density excluded the extreme values. The data from the 234 plots sampled in 2010 are shown.
Table 1. Classification of habitat factors affecting grasshopper density in Xianghuangqi County.
Regarding the categorical habitat factors (such as aspect, soil type, vegetation type, and land cover type), the initial classifications were determined based on the general habitat classification information and field surveys (fig. 2). Student's t test was used to determine the differences in grasshopper density that were associated with the initial classifications of a habitat factor. The classifications between which there were significant differences were retained, and the classifications between which there were no significant differences were adjusted or merged until all the classifications were associated with significant differences in grasshopper density (table 1). In addition, the rare categories with very low grasshopper density were combined into the single category ‘Other’. Aspect was divided into five categories: sunny slope, half sunny slope, half shady slope, shady slope, and flat terrain. The slopes on the north side were defined as having a slope of 0°, and the other slopes successively increased from 0° in a clockwise direction. Sunny slope was defined as a slope of 157.5°–247.5°, half sunny slope was 112.5°–157.5° and 247.5°–292.5°, half shady slope was 67.5°–112.5° and 292.5°–337.5°, and shady slope was 0°–67.5° and 337.5°–360°. Flat terrain was expressed as −1.
Evaluation of grasshopper habitat suitability
Effects of individual habitat factors on grasshopper density
The factor detector in the Excel–GeogDetector can be used to investigate the potential risk factors and to test whether a spatial pattern of health risk depends on one or multiple factors.
The health effect H (e.g., heart disease rate, cancer mortality ratio) is taken as a health risk index. The weighted dispersion variance of H of each spatial category of the determinant D (Var d ) is compared with the global variance of H for the same D over the entire study area (σ2 H ). If σ2 H is high and Var d is low, the classified categories of D explain the variations in H well, and D has a strong contribution to H. This effect can be expressed as Power of Determinant (PD):
$$PD = 1 - \displaystyle{{Var_d} \over {{\rm \sigma} _H^2}} = 1 - \displaystyle{{\sum\nolimits_{i = 1}^L {(n_{{\rm D},i} {\rm \sigma} _{H_{{\rm D},i}} ^2} )} \over {n{\rm \sigma} _H^2}} $$
$$n = \sum\limits_{i = 1}^L {n_{{\rm D},i}} $$
where L is the category number of D, n
D,i
is the number of samples associated with the category i of D, and
$\sigma ^2 _{H_{{\rm D},i}} $
is the divisional variance of H for the category i of D. Generally, PD ranges from 0 to 1. The higher the PD value, the greater the effect of D on H. If PD equals 1, the spatial distribution of D controls the spatial pattern of H completely; if PD equals 0, the spatial distribution of D does not affect the spatial pattern of H at all (Wang et al., Reference Wang, Liao and Liu2010; Hu et al., Reference Hu, Wang, Li, Ren and Zhu2011).
In this study, H indicated grasshopper density; D indicated relatively stable habitat factors that might affect grasshopper density. The categories of each habitat factor indicated the corresponding classifications of the factor (table 1). When there were significant differences in grasshopper density among the different categories of a habitat factor, σ2 H was very high, and when there were insignificant differences in grasshopper density for each category of this factor, Var d was very low. When these conditions existed, there was a synergistic relationship between grasshopper density and this habitat factor, and the PD value, which indicates the effect of this factor on grasshopper density, was very high.
Effects of two-factor interactions on grasshopper density
The interaction detector in the Excel–GeogDetector can be applied to test the effects of interactions between two risk factors on disease prevalence. In this study, it was used to analyse the effects of interactions between any two-habitat factors on grasshopper density. One category of a habitat factor was combined with one category of another habitat factor to generate a new classification. The weighted dispersion variance of grasshopper density within all the new combinations for the two factors (Var d ) and the global variance of grasshopper density for the same two factors over the entire study area (σ2 H ) were calculated to obtain the PD value using equation (1). A higher PD value indicated that the interactions between a certain two habitat factors had a greater effect on grasshopper density.
Membership analysis of grasshopper density
Membership analysis was used to express information about the effects of habitat factors on grasshopper occurrence. Here, the membership Q refers to the probability that a habitat category is associated with a given grasshopper density rank, ‘very high,’ ‘high’, ‘low’, or ‘very low’. Given that grasshoppers live in a variety of habitats, and the effects of habitat factors on their distribution are relatively uncertain, the weights of the four ranks of grasshopper density R were set to 1, 0.7, 0.4, and 0, respectively. For each habitat category, Q was statistically calculated using our field surveys and measurements, and then the weighted mean membership for each habitat category was obtained from Q and R.
Suitability weight analysis for each habitat category
The analytical hierarchy process (AHP) was used to analyse the suitability weight for each category of each habitat factor. Building the grasshopper habitat suitability evaluation system was a target layer, the seven habitat factors were criterion layers, and a habitat category was an index layer. For the habitat factor i, the judgment matrix A i was constructed according to the relationships between M p and M q , which denoted the weighted means of the memberships of the four grasshopper density ranks for any category p and q.
$${\bi A}_i = \left[ {\matrix{ {A_{11}} & {A_{12}} & \cdots & {A_{1q}} & \cdots & {A_{1n}} \cr {A_{21}} & {A_{22}} & \cdots & {A_{2q}} & \cdots & {A_{2n}} \cr {A_{31}} & {A_{32}} & \cdots & {A_{3q}} & \cdots & {A_{3n}} \cr \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \cr {A_{\,p1}} & {A_{\,p2}} & \cdots & {A_{\,pq}} & \cdots & {A_{\,pn}} \cr \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \cr {A_{n1}} & {A_{n2}} & \cdots & {A_{nq}} & \cdots & {A_{nn}} \cr}} \right]$$
In this matrix, n is the number of the categories of habitat factor i as well as the order of A i .
$$A_{\,pq} = M_p /M_q $$
A pq > 1 indicates that the category p has greater importance than q.
$${\rm \beta} _p = \left( {\prod\limits_{q = 1}^n {A_{\,pq}}} \right)^{{1 / n}} (\,p = 1,2,3, \ldots n)$$
$${\rm \alpha} _p = \displaystyle{{{\rm \beta} _p} \over {\sum\nolimits_{p = 1}^n {{\rm \beta} _p}}} $$
$${\bi B} = ({\rm \alpha} _1, {\rm \alpha} _2, \ldots, {\rm \alpha} _n )$$
The consistency of the matrix A i was tested by comparing the consistency index (CI) of A i with the same order average random consistency index (RI). If the largest eigenvalue of A i (λmax) is not equal to n, the consistency ratio (CR) may be calculated as
$${\rm CR} = \displaystyle{{{\rm CI}} \over {{\rm RI}}} = \displaystyle{{({\rm \lambda} _{\max} - n)/(n - 1)} \over {{\rm RI}}}$$
where RI is 0.58, 0.90, and 1.12 as n is equal to 3, 4, and 5, respectively. The matrix A i has better consistency only if CR < 0.1, when the vector B may be deemed as the suitability weight set for the n categories of habitat factor i (Li, Reference Li2009; Deng et al., Reference Deng, Li, Zeng, Chen and Zhao2012).
Model to evaluate grasshopper habitat suitability
We improved the equation estimating fuzzy evaluation scores that was initially proposed by Thomas et al. (Reference Thomas, Wotherspoon and Raymond2012) to assess grasshopper habitat suitability by fully considering the effects of individual habitat factors and the interactions between any two factors:
$$S_k = \sum\limits_{i = 1}^7 {P_i} Q_{ki} + \sum\limits_{i = 1}^7 {\sum\limits_{\,j = 1}^7 {P_{ij}}} \log \displaystyle{{P_{ij}} \over {Q_{ki} Q_{kj}}} \quad (i \ne j)$$
where i and j indicate habitat factors; k is the rank of grasshopper density, that is, very high (k = 1), high (k = 2), low (k = 3), and very low (k = 4); S k is the fuzzy evaluation score of habitat suitability on the kth rank of grasshopper density; P i is the PD value indicating the effect of the individual factor i on grasshopper density; P ij is the PD value indicating the effect of the interaction between two factors i and j on grasshopper density; and Q ki and Q kj are the memberships of the kth rank of grasshopper density for the categories of factors i and j, respectively. The first item of the equation indicates the integrated effect of individual habitat factors on grasshopper density; the second item indicates the integrated effect of the two-factor interactions. Q ki and Q kj are set to 0.0001 if they are equal to 0.
$$S = \sum\limits_{k = 1}^4 {S_k R_k} $$
where S is the fuzzy evaluation score of grasshopper habitat suitability, and higher S values mean more suitable habitat. R k is the weight of the kth rank of grasshopper density.
The partial field data were used for constructing the evaluation system, including the data from the 234 plots sampled in 2010 and 282, 275, and 309 plots selected from the total 342, 335, and 369 plots sampled during the first field survey in 2011, 2012, and 2013, respectively. We calculated the P i values for the seven individual habitat factors using the factor detector in the Excel–GeogDetector and their P ij values using the interaction detector based on all of the selected plot data. We then calculated the memberships of the four grasshopper density ranks (Q 1, Q 2, Q 3, and Q 4) for each habitat category in each distinct year from 2010 to 2013. Next, the evaluation scores S 1, S 2, S 3, and S 4 for all these plots were estimated using an equation (9) and the final score S was estimated using another equation (10).
Results
Effects of habitat factors on grasshopper density
The results of using the factor detector in the Excel–GeogDetector showed that the PD values for the seven habitat factors based on the 2010–2013 field data decreased in the following order: elevation (0.440), soil sand content (0.418), aspect (0.328), slope (0.203), vegetation type (0.061), soil type (0.044), and land cover type (0.027).
The results of using the interaction detector in the Excel–GeogDetector showed that the effects of the interactions between any two habitat factors on grasshopper density (P ij ) were greater than the effect of either individual habitat factor (P i or P j ). In all the pairwise combinations, the interaction between elevation and soil sand content most significantly affected grasshopper density. There was a nonlinear and synergistic relationship between elevation and soil sand content; the effect of their interaction was greater than the sum of the independent effects of the two individual factors. The interaction between elevation and aspect exhibited the second greatest effect on grasshopper density, followed by the interaction between elevation and slope. The effects of the interactions between soil type and soil sand content, vegetation type, and land cover type on grasshopper density were greater than the effect of any individual factor, and they were smaller than the sum of the effects of the two individual factors. The interaction between vegetation type and land cover type exhibited a similar relationship (table 2). In short, habitat factors do not independently or linearly affect grasshopper density.
Table 2. Effects of two-factor interactions on grasshopper density in Xianghuangqi County.
Note: The values in the parentheses are the PD values, indicating the effects of individual habitat factors on grasshopper density. P i ↑P j indicates that the effect of the interaction between the habitat factors i and j was greater than the sum of the independent effects of the two individual factors. P i ↑↑P j indicates that the effect of the interaction between the habitat factors i and j was greater than the effect of any individual factor and smaller than the sum of the effects of the two individual factors. The calculations were based on all of the selected plot data from 2010 to 2013.
Memberships of grasshopper density ranks and suitability weights for each habitat category
The membership analysis and habitat suitability weight analysis exhibited similar relationships between grasshopper density and habitat category for the 4 years, although the memberships themselves were different. The results showed that the grasshoppers occurred most often at elevations of 1300–1400 m, followed by elevations of 1200–1300 m. Flat terrain was most suitable for grasshoppers, followed by the sunny, half sunny, and half shady slopes. A slope of 4–6° was most suitable for grasshoppers, followed by a slope of 6–8°. Grasshoppers occurred most often in typical chestnut soil, followed by dark and light chestnut soil. Grasshoppers occurred most often at sites with soil sand content of 70–80%; sites with soil sand content that was lower than 60% or higher than 80% were not suitable for grasshoppers. Temperate bunchgrass steppe dominated by S. krylovii and C. squarrosa was noticeably favoured by grasshoppers. Non-grassland sites were not suitable for grasshoppers. Medium- and low-coverage grasslands were most suitable for grasshoppers, followed by high-coverage grassland in 2010 and 2011, but high-coverage grassland was most suitable for grasshoppers, followed by medium-coverage grassland in 2012 and 2013. In table 3, we have included only the results from 2011 because of space constraints.
Table 3. Membership of grasshopper density rank and suitability weights for each category of habitat factors in Xianghuangqi County, 2011.
1 Vegetation type: 1, Stipa krylovii and Cleistogenes squarrosa steppe; 2, S. grandis and Artemisia frigida steppe; 3, S. grandis, Leymus chinensis, and Salsola collina steppe.
Considering the effects of the two-factor interactions, grasshoppers were most likely to inhabit sites with an elevation of 1300–1400 m and soil sand content of 70–80% when other conditions were similar. Sites that were flat with an elevation of 1300–1400 m, or sites with an elevation of 1300–1400 m and a slope of 4–6°, also showed a high probability of containing more grasshoppers when other conditions were similar.
Evaluation and validation of grasshopper habitat suitability
The results showed that the modelled evaluation scores of grasshopper habitat suitability (S) corresponding to the four grasshopper density ranks (>20, 10–20, 5–10, and <5 ind. m−2) ranged from 466 to 545 (S 1) for 93.3% plots in 2010, 85.7% plots in 2011, and 100% plots in 2012 and 2013, 361–465 (S 2) for 85% plots in 2010, 100% plots in 2011, 88.9% plots in 2012, and 92% plots in 2013, 251–360 (S 3) for 74.1% plots in 2010, 93.1% plots in 2011, 82.0% plots in 2012, and 88.2% in 2013, and 170–250 (S 4) for 90.4% plots in 2010, 90.8% plots in 2011, 91.2% plots in 2012, and 89.8% in 2013, respectively.
We found that few grasshoppers were present at some suitable plots during the first field survey, before they could move to other sites, in 2012 and 2013, whereas there were many grasshoppers at the same plots during the second field survey, after they might have moved. This indicated that the grasshoppers did not oviposit at these plots in the previous year although the plots were suitable, but later some adult grasshoppers moved to these plots. For these plots, we used the values of grasshopper density obtained during the second field survey to calibrate those obtained during the first field survey when assessing the modelled S values. After calibration, the percentage of plots with the evaluation score S 3 and the corresponding grasshopper density rank (5–10 ind. m−2) increased from 82.0 to 83.6% in 2012 and from 88.2 to 90.6% in 2013, and that with S 4 and the corresponding grasshopper density rank (<5 ind. m−2) increased from 91.2 to 93.1% in 2012 and from 89.8 to 93.4% in 2013. In 2010 and 2011, we only conducted a survey around early-to-mid-July; hence, this calibration could not be made.
The modelled S values appeared to provide better indicators of grasshopper density ranks for more than 75.0% plots, except for the modelled S 3 values in 2010. For some plots, the modelled S 3 values were higher, even reaching the range of the modelled S 2 values, suggesting that these plots had more suitable habitat conditions but low grasshopper density. There are several possible reasons for this result: (1) Human efforts to prevent or control grasshopper infestations at some suitable plots reduced the number of grasshoppers living there. This situation could account for about 13.3% of the unexpected results. (2) Some plots with suitable habitat were surrounded by trees that might have harboured birds, which are predators of grasshoppers. This case could account for about 20.0% of the unexpected results. (3) Some suitable plots were adjacent to high hills or main roads and were isolated from other suitable sites containing many grasshoppers. Therefore, these plots might have been prevented from becoming oviposition sites or dispersal targets. This case could account for about 26.7% of the unexpected results. (4) Grasshoppers did not arrive at certain suitable plots and oviposit in the previous year, so that there were few nymphs at these plots before mid-July. Some of these plots had more grasshoppers when individuals immigrated from other sites after mid-July. This has been verified by the data from the second field survey in 2012 and 2013.
The 180 plots sampled in 2011–2013, which were not involved in the creation of the evaluation system, were used to validate the reliability of the finished model. The modelled S values were compared with the measured values of grasshopper density (Gd) for the 180 plots. The results showed that the sites with higher S values also had higher Gd (fig. 3). When the modelled S values were within the range of S 1, 100% plots had a Gd that was higher than 20 ind. m−2. When the modelled S values were within the range of S 2, 92.9% plots had a Gd of 10–20 ind. m−2. When the modelled S values were within the range of S 3, 87.9% plots had a Gd of 5–10 ind. m−2, and when the modelled S values were within the range of S 4, 85.2% plots had a Gd of fewer than 5 ind. m−2.
Fig. 3. Relationship between the modelled fuzzy evaluation scores of grasshopper habitat suitability and the measured grasshopper density in Xianghuangqi County. The data from the 180 plots sampled in 2011–2013 are shown, which were not involved in the creation of the evaluation system.
Discussion
Improvement of existing model system to evaluate grasshopper habitat suitability
The good agreement between the modelled S values and the measured values of grasshopper density demonstrated that the new habitat suitability evaluation system was reliable and that the modelled S values could be applied to indicate the statuses of grasshopper occurrence in different years in the Inner Mongolia grassland. The new model system improved and perfected the previous one that we constructed in our earlier study (Zhang et al., Reference Zhang, Zhang and Chen2012).
Zhang et al. (Reference Zhang, Zhang and Chen2012) used a multi-objective linear weighted function to evaluate the habitat suitability for POG.
$$S = \sum\limits_{k = 1}^4 {\left[ {\sum\limits_{i = 1}^7 {Q_{ki} W_i}} \right]R_k} $$
where W i is the suitability weight for the habitat factor i. The meaning of k, i, S, Q ki , and R k are the same as those in equations (9) and (10).
In the present study, the effects of all the related habitat factors on grasshopper density were evaluated using the factor detector in the Excel–GeogDetector without making any assumptions. After elevation, slope had the fourth strongest effects on grasshopper density, indicating that this factor should not be eliminated from the analysis. The effects of the two-factor interactions on grasshopper density were also quantified using the interaction detector in the Excel–GeogDetector, which provides a new perspective on the study of grasshopper habitat suitability. The results showed that the two-factor interactions significantly enhanced the effects of individual habitat factors on grasshopper density. Thus, it was unreasonable or even incorrect to simplify them as a random error. In short, the improved model was more comprehensive and more systematic than the model built by Zhang et al. (Reference Zhang, Zhang and Chen2012).
Zhang et al. (Reference Zhang, Zhang and Chen2012) used a binary comparison method to obtain the suitability weights (W) for all the habitat factors that were included in the model. The present study used the AHP method to obtain the suitability weights for each category of each habitat factor, rather than the weights of the habitat factors themselves.
Another improvement made by the present study was the reclassification of the seven habitat factors. The classifications of elevation, soil sand content, soil type, and land cover type were similar in Zhang et al. (Reference Zhang, Zhang and Chen2012) and the present study. However, the new classifications of aspect and vegetation type were quite different from the previous ones. The aspect in the present study was divided into five groups ranging from sunny to shady slope, whereas the previous groups were based on nine distinct directions, which might break apart associated measurements. The new vegetation types were based on the surveyed dominant and subdominant species, whereas the previous ones were based on the digitized vegetation type map with a cartographic scale of 1:1,000,000.
However, we note that the model directly evaluates the suitability of habitat, which is closely related to POG, but which may not indicate actual grasshopper occurrence. Taking other factors (such as the actions of humans to prevent or control grasshopper populations, grazing of livestock, the proximity of trees and birds that are predators of grasshoppers, and the isolation of a given site from other suitable sites) into account and calibrating the results using further data in the future are expected to improve the accuracy of the model to indicate the habitat suitability, as well as more closely reflect the actual grasshopper density.
Explanations of the relationship between habitat factors and potential spatial distribution of grasshoppers
The spatial heterogeneity of POG within Xianghuangqi County mainly arose from unevenness in host plants and habitat patterns. The effects of individual habitat factors on grasshopper density varied. The results of the present study demonstrated that grasshopper density was associated most closely with elevation, soil sand content, and aspect, which might be due to the much greater spatial heterogeneity of these habitat factors, relative to the other factors studied. Meanwhile, grasshopper density was associated least closely with land cover type, soil type, and vegetation type. The seemingly paradoxical result might be due to the relatively homogeneous landscape, mainly dominated by temperate bunchgrass steppe and chestnut soil, in which differences in the effects of land cover or vegetation type and soil type on grasshopper dynamics may not be apparent.
Elevation directly affects the temperature, humidity, and wind speed of sites, which can affect soil and vegetation. In Xianghuangqi County, sites with an elevation of more than 1400 m are characterized as having lower air and soil temperatures, higher wind speed, larger pieces of gravel, and weaker plant growth (including the absence of any plants at all). These sites cannot meet the most basic food and ovipositing needs of grasshoppers. By comparison, sites with an elevation of less than 1200 m have higher air and soil temperatures, and when soil water is plentiful they can be densely vegetated, factors that do not favour grasshopper oviposition and incubation (Rourke, Reference Rourke2000; Zhang et al., Reference Zhang, Zhang and Chen2012).
Aspect and slope vary at smaller spatial scales and generally affect the local spatial distribution of grasshoppers. Aspect directly affects environmental conditions related to heat, including solar radiation and temperature. O. d. asiaticus favours warm soil for its oviposition, overwintering, and incubation (Chen, Reference Chen2007). The effects of slope on temperature and water conditions are more complicated. For example, the lower solar radiation and temperatures at steeper slopes are not conducive to O. d. asiaticus; on the other hand, although there are better temperature conditions at some gentler slopes, these sites have poor drainage of surface water from summer rainstorms and spring snowmelt, which reduces their suitability for O. d. asiaticus (Chen et al., Reference Chen, Hao, Pang and Kang2009). Thus, taking into consideration the combined effects of aspect and slope, flat sites appear to be most suitable for O. d. asiaticus, followed by gentler sunny slopes.
O. d. asiaticus occurred only in chestnut soil, whereas meadow soil, lithosol, aeolian sandy soil, and saline-alkali chestnut soil were not suitable for O. d. asiaticus. Meadow soil with higher soil water content and lower soil temperature, stony and hard lithosol, loose aeolian sand soil with poor water holding capacity, and saline-alkali chestnut soil with higher salinity and pH values were unsuitable for the ovipositing, overwintering, incubating, or development needs of most grasshopper species (Edwards & Epp, Reference Edwards and Epp1965; Ni et al., Reference Ni, Wang, Jiang and Zha2007).
Sites with soil sand content of 70–80% were most suitable for O. d. asiaticus. Soil with too high sand content does not hold water well, and does not favour grasshopper oviposition and incubation. Soil with low sand content has poor water and air permeability and low temperature, hardens easily, and does not favour grasshopper oviposition, overwintering, or incubation (Chen, Reference Chen2007; VanDyke et al., Reference VanDyke, Latchininsky and Schell2009).
Plants can provide many resources required by grasshoppers, such as food, appropriate microclimate, habitat sites, and refuge from predators (Torrusio et al., Reference Torrusio, Cigliano and de Wysiecki2002). The spatial distributions of grasses have the greatest effects on grasshopper occurrence (Zhao et al., Reference Zhao, Zhou, Wang, Shi and Gao2011). In Xianghuangqi County, with the exception of the steppes dominated by S. krylovii and C. squarrosa or by Stipa grandis and A. frigida, there were small numbers of O. d. asiaticus in temperate bunchgrass and forb meadow steppe, where the dominant plant species were S. grandis, Leymus chinensis, and Salsola collina. O. d. asiaticus is a grass-feeder, so it is not common at the sites dominated by non-grass vegetation (such as shrubs or forbs) (Chen, Reference Chen2007; Lu et al., Reference Lu, Han and Zhang2008).
These relationships imply that appropriate management practices that change habitat features, such as vegetation type and soil sand content, could reduce the magnitude and frequency of grasshopper infestations, especially at sites with a high POG rank. For example, moderate livestock grazing or rotational grazing decreases the abundance of dominant grasses; sowing green manure crops in rotations improves soil texture and reduces soil sand content below 60 or 50%.
It should be noted that the study selected only relatively stable habitat factors and excluded soil water content and vegetation coverage that exhibit very large responses to climate change in spring and summer, although they may have great effects on grasshopper occurrence. We found that there existed apparent differences in the PD values of soil water content between 2010/2011 (0.346/0.226) and 2012/2013 (0.668/0.704) because of the differences in frequency and magnitude of precipitation within the 20 days before we surveyed. We also found that the relation patterns of weighted mean memberships for different categories of most of the habitat factors were similar in different years, but the patterns varied dramatically for the factors related to vegetation coverage. For example, 30–50% vegetation coverage had the maximum weighted mean membership among the five coverage categories in 2010 and 2011, while it increased to 50–70% in 2012 and >70% in 2013, because vegetation coverage closely associated with precipitation was higher and more plots had high vegetation coverage in 2012 and 2013 than in 2010 and 2011. Therefore, the membership relationship between high-coverage (>50%) grassland and medium-coverage (30–50%) grassland exhibited interannual variations. In the future, we intend to couple soil water content and vegetation coverage with the relatively stable habitat factors when evaluating climate suitability.
In addition, we intend to expand the evaluation to the entire study area, with the aid of spatial maps of habitat factors, and to combine our evaluation with climate suitability evaluation. The availability of data on topography, vegetation, and soil factors and grasshopper density in association with ecological factors will allow the evaluation system to be extended to other regions.
Acknowledgements
This research is financially supported by the National Natural Science Foundation of China (Grants NO. 31270512). We gratefully acknowledge Jinfeng Wang and Chengdong Xu for providing the Excel–GeogDetector software. We are grateful to Ze-Hua Zhang and Hui-Hui Wu from the Institute of Plant Protection, Chinese Academy of Agricultural Sciences for their helps conducting the field survey. We also gratefully acknowledge the editor and anonymous reviewers for their valuable comments and constructive suggestions. Finally, we would like to thank Editage (http://www.editage.cn/index.html) for English language editing.
