Geographic Variation in Migratory Grasshopper Recruitment under Projected Climate Change

Abstract

Climate change is expected to alter prevailing temperature, precipitation, cloud cover, and humidity this century, thereby modifying insect demographic processes and possibly increasing the frequency and intensity of rangeland and crop impacts by pest insects. We leveraged ten years of migratory grasshopper (Melanoplus sanguinipes) field surveys to assess the response of nymph recruitment to projected climate conditions through the year 2040. Melanoplus sanguinipes is the foremost pest of grain, oilseed, pulse, and rangeland forage crops in the western United States. To assess nymph recruitment, we developed a multi-level, joint modeling framework that individually assessed nymph and adult life stages while concurrently incorporating density-dependence and accounting for observation bias connected to preferential sampling. Our results indicated that nymph recruitment rates will exhibit strong geographic variation under projected climate change, with population sizes at many locations being comparable to those historically observed, but other locations experiencing increased insect abundances. Our findings suggest that alterations to prevailing temperature and precipitation regimes as instigated by climate change will amplify recruitment, thereby enlarging population sizes and potentially intensifying agricultural pest impacts by 2040.

1. Introduction

Climate is the dominant force driving biogeographic patterns at large spatial scales [1,2] and exerts an overriding influence on insect distributions, abundances, and demographic processes [3,4,5,6,7]. Through a combination of anthropogenic forcing and natural variability, climate change is expected to alter prevailing temperature and precipitation patterns this century, propelling departure from historic conditions and potentially increasing the regularity and severity of extreme weather events in the western United States (US) [8,9,10,11]. Graduated climate change occurring over multiple years may alter soil nitrogen deposition rates, decrease near-surface carbon stores, and limit soil moisture availability as well as alter plant nutrient content integral to phytophagous insect ecology [12,13]. Over shorter time frames, weather extremes have the capacity to affect insect demographic structure and instigate population booms or crashes within or between seasons [14,15,16]. Climate change has wide ranging implications for insect populations [17,18,19,20], particularly for insect groups like grasshoppers (Orthoptera: Acrididae) that hold central roles in ecosystem functioning and agricultural pest management.

Grasshoppers are simultaneously the principal grassland invertebrate herbivores [21] and the most impactful rangeland insect pests in the US [22]. Although vital to nutrient processing and ecosystem function [23], grasshopper population densities often show dramatic increases that beget substantial economic harm to cultivated crops and rangeland forage [22,24]. Despite their potential as focal species for climate change research, enormous uncertainty exists in foreseeing how climate change may affect grasshopper populations. As examples, climate change may diminish critical grasshopper contributions to ecosystem function, amplify the extent and intensity of grasshopper outbreaks, or otherwise modify grasshopper population dynamics to impact food security and agricultural economic viability [13,20,25,26]. As an initial step in understanding how climate change may affect grasshopper populations in the US, process models were developed to forecast future reproductive rates for the migratory grasshopper (Melanoplus sanguinipes).

M. sanguinipes (Msang) is found throughout North America and feeds on a wide variety of plants and crops. As a pest species, Msang is frequently recognized as the foremost pest of grain, oilseed, pulse, and rangeland forage crops in the US, but also impacts wheat, barley, oats, vegetables, and ornamentals when grasshoppers are at high density [7,27,28]. Under typical conditions, adult Msang grasshoppers (gh) reach densities ranging from 0 to 3 gh/m${}^2$ and nymphs (immature juveniles) may attain densities between 0 and 7 gh/m${}^2$; however, during population surges and outbreaks, total densities may be well over 50 gh/m${}^2$ [29,30]. Msang population dynamics are shaped by a three-stage life cycle [29], in which nymphs hatch in the spring, progress through five instar developmental phases over the period of about a month, and then emerge as adults during summer. Adult females are fecund soon after emergence (2–3 weeks) and deposit pods (egg clutches) in the soil, where soil temperature and moisture largely govern the rate of embryonic development before winter diapause [31,32]. Msang is typically univoltine in Canada and the US; however, multivoltine populations have been identified in southern portions of its range [33]. In addition to mediating embryonic development, temperature and moisture availability control hatch timing, food plant availability, and behavioral aspects of Msang feeding, migration, and reproductive ecology [16,34,35,36,37]. Modifications to prevailing temperature and precipitation regimes as prompted by climate change are clearly of major consequence to Msang population dynamics and ultimately to US and Canada agricultural interests [38].

The objective of the current study was to leverage ten years of Msang field survey data to assess the response of nymph recruitment under projected climate conditions through the year 2040. Nymph recruitment is a more variable, context-dependent, and sensitive measure of population change than is total abundance [39]. To achieve this objective, a multi-level, joint modeling framework was developed that individually assessed nymph and adult life stages while concurrently accounting for observation bias connected to preferential sampling. To ensure that density-dependence between Msang life stages was realistically captured, the model incorporated shared spatiotemporal effects that permitted adult–nymph interactions to inform demographic estimates. After detailing model construction, implications for future Msang population dynamics are discussed.

2. Materials and Methods

2.1. Study Domain and Observation Data

The State of Wyoming in the western United States (US) served as the study domain for analysis (Figure 1). Wyoming is geographically located between 41${}^{\circ }$ and 45${}^{\circ }$ North Latitude and −104${}^{\circ }$ and −111${}^{\circ }$ West Longitude. Wyoming has an areal extent of approximately 253,596 km${}^2$ and exhibits an elevation range from 945 m at the Belle Fourche River to 4209 m at Gannett Peak. Major ecoregions in the study area include the Northwestern Great Plains, Western High Plains, the middle Rocky Mountains, and intermontane basin [40].

The study incorporated a total of 1772 point-level Msang observations collected during field surveys between 2011 and 2020. Data were collected from approximately 1100 unique Wyoming locations by staff from the US Department of Agriculture, Animal and Plant Health Inspection Service, Plant Protection and Quarantine (PPQ) program. Data attributes provided with the data indicated the sample location (longitude and latitude) and grasshopper abundance counts for nymph (874 locations) and adult (898 locations) life stages. Grasshopper sample locations are illustrated in Figure 2.

2.2. Climate and Environmental Data

Due to high levels of uncertainty in climate forecasts [41], a consensus approach [42,43] was applied using three different climate models developed in conjunction with the sixth phase Coupled Model Intercomparison Project (CMIP6) [44]. Employing the model independence criteria developed by Sanderson et al. [45], Sanderson et al. [46], three terrestrial Global climate models (GCM) were adopted for analysis, these included the IPSL-CM6A-LR [47], CanESM5 [48], and the MIROC6 [49]. The IPSL-CM6A-LR, CanESM5, and MIROC6 GCMs were developed by different research groups and were found to exhibit elevated inter-model pairwise distances in multidimensional space (above the mean as assessed across 38 different models, see Sanderson et al. [45]), suggesting model independence. Nineteen bioclimatic variables [50] were derived at 2.5 min spatial resolution from each GCM projection based on four Shared Socioeconomic Pathways (SSPs) SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 for the time period horizon 2021–2040. Brief descriptions for all bioclimatic variables are provided in the Supplementary Information. The SSPs are future climate scenarios under various land-use, population growth, and energy consumption assumptions. As described by [51], the SSPs represent outlooks that range from sustainable growth (SSP1-2.6) or a middle of the road path with reduced climate vulnerability (SSP1-2.6), to future settings in which climate mitigation is a low international priority (SSP3-7.0) and fossil fuel use is intensive (SSP5-8.5). Climate forecasts were down-scaled and bias corrected using WorldClim (v2.1) as a baseline, which corresponds to 30-year average climate conditions [50].

Study area edaphic conditions were characterized using 285 soil variables (250 m resolution) obtained from the SoilGrids database hosted by the International Soil Reference and Information Centre ([52], https://www.isric.org/ (accessed on 15 December 2021)). Soil data quantified a wide variety of soil attributes across a range of horizon depths, including but not limited to bulk density, total nitrogen, organic carbon concentration, pH, cation exchange capacity, texture fractions, and more complex characters like organic carbon stocks in topsoils and subsoils. Brief descriptions for all soil variables are provided in the Supplementary Information.

Vegetation variability was assessed using thirteen habitat heterogeneity metrics that enumerated composition and textual variation within Enhanced Vegetation Index (EVI) imagery collected by the Moderate Resolution Imaging Spectroradiometer (MODIS) [53]. The data were downloaded at 30 arc-second resolution from the Earth Environment website (http://www.earthenv.org/ (accessed on 15 December 2021)). The data layers were derived from 16-day MODIS EVI (MOD13Q1 version 5) composites captured between between 2001 and 2005. The data set quantified land surface greenness during peak growing season based on first- and second-order texture metrics, which respectively describe the frequency distribution of EVI pixel values within a defined neighbourhood and the probability of observing a pair of pixel values with a given inter-pixel distance and orientation. The texture metrics are not direct plant species composition measures; rather, they index vegetation spatial variability and arrangement [53]. Brief descriptions for all vegetation variables are provided in the Supplementary Information.

2.3. Variable Decomposition

Owing to the problematic nature of drawing biologically relevant, mechanistic inference from species distribution models [54,55], a datamining approach was implemented to prioritize model performance over variable-specific interpretation of species–environment relationships. To accomplish this, climate (19 variables), soil (285 variables), and vegetation (13 variables) data sets were decomposed to create three synthetic covariates respectively summarizing total climate, soil, and vegetation variation. Although Principal Component Analysis (PCA) has been successfully used to summarize total variance within a study area as a whole [56,57,58,59,60], this common practice was modified and expanded to instead identify the linear combinations (weighted principal components) that best described climate, soil, and vegetation variance with respect to Msang specifically. Principal component weights were estimated using discriminant analyses [61], which provided geographic context by quantifying the proportion of environmental variance [62] that best distinguished locations with grasshopper occurrence from random background locations (prior groups). To maximize discriminate power and avoid over-fitting, Msang occurrences used for initial discriminant analysis were then compared to repeated and randomized samples generated with the a-score function available in the adegenet package [63].

2.4. Statistical Model

Msang nymph and adult life stages were jointly modeled as log-Gaussian Cox processes (LGCP) under preferential sampling. In our application, the LGCP were modeled as non-homogeneous, point processes using individual point locations and associated, point-specific attributes (“marks”) without any aggregation to grid cells or other areal units. Preferential sampling refers to data that is potentially biased due to having been opportunistically collected [64]. Msang observation data were not collected through randomized sampling across the study area as a whole (e.g., cross-sections of good and bad habitat), rather data were acquired during field surveys at locations suspected in advance to be suitable for grasshoppers; that is, Msang field sites were preferentially sampled. Preferential sampling results in point patterns (non-random sampling locations) that may not be independent of the species biological or ecological processes under investigation, thus becoming statistically problematic [64,65,66]. For example, preferential sampling during field survey may have produced observation data in which areas of high Msang abundance are overrepresented. To help account for non-independence between grasshopper survey locations (pattern) and grasshopper abundance (process), the number of sample points per unit area, or point pattern intensity (${\mathsf{\Lambda }}_{st}^{a,n}$), was concurrently estimated with adult ($Y_{st}^a$) and nymph ($Y_{st}^n$) abundances in a four-part, joint model with shared spatiotemporal effects. LGCP are point-based models that enable point intensities to be estimated across continuous surfaces as Gaussian random fields. Adopting the LGCP approach, point patterns were specified as:

$$\begin{array}{cc}\hfill \log \left({\mathsf{\Lambda }}_{st}^{a,n}\right)& ={\mathrm{W}}_{st}^{a,n}\hfill \\ \hfill {\mathrm{W}}_{st}^{a,n}& \sim \mathrm{N}(0,\mathrm{Q}(\kappa ,\tau \left)\right)\hfill \end{array}$$

where the logarithm of intensity (${\mathsf{\Lambda }}_{st}^{a,n}$) for grasshoppers in the adult (superscript a) and nymph (superscript n) stages at each geographic location s$(\mathit{s}=1,2,3,\dots ,n)$ and year t $(\mathit{t}=2011,2012,2013,\dots ,2020)$ were estimated as Gaussian random fields (${\mathrm{W}}_{st}^{a,n}$). Following Lindgren et al. [67] and Simpson et al. [68], the matrices $\mathrm{Q}(\kappa ,\tau )$ used to define the spatial fields (${\mathrm{W}}_{st}^{a,n}$) were approximated using stochastic partial differential equations (SPDE) that facilitated implicit estimation of spatial range ($\kappa$) and scale ($\tau$) parameters based on Matérn covariance. For a detailed statistical description of SPDE methods see Lindgren et al. [67], Krainski et al. [69], or other applications by the current authors [70,71,72,73,74].

In addition to jointly modeling the point pattern and abundance specific to each the adult and nymph stage, population density-dependence dynamics necessitated that nymph estimates be conditional on those made for adults in the prior year. Density-dependent effects can cause Msang abundances in one year to be affected by population numbers from the preceding year [75]. Due to density-dependence, the model was designed such that information from the first two levels, which evaluated adult patterns and processes, was shared with the fourth level used to estimate nymph abundance. More formally, the first two model tiers were,

$$\begin{array}{c}\hfill \mathrm{Adults}=\left\{\begin{array}{c}\mathrm{Pattern}\left\{\begin{array}{c}\log \left({\mathsf{\Lambda }}_{st}^a\right)={\mathrm{W}}_{st}^a\hfill \\ {\mathsf{\Lambda }}_{st}^a=\exp \{{\beta }_{\mathsf{\Lambda }}^a+W_{st}^a\}\hfill \end{array}\right.\hfill \\ \\ \mathrm{Process}\left\{\begin{array}{c}{\mathrm{Y}}_{st}^a\sim \mathrm{Poisson}\left({\mu }_{st}^a\right)\hfill \\ \log \left({\mu }_{st}^a\right)={\beta }_0^a+{\beta }_1^{a}X+{\alpha }^1{\mathrm{W}}_{st}^a\hfill \end{array}\right.\hfill \end{array}\right.\end{array}$$

where adult point pattern intensity (${\mathsf{\Lambda }}_{st}^a$) was approximated as the exponential of a Gaussian random field (${\mathrm{W}}_{st}^a$) and an intercept (${\beta }_{\mathsf{\Lambda }}^a$) in the first tier, and the second tier estimated adult counts ($Y_{st}^a$) following a Poisson distribution. The log of mean adult abundance (${\mu }_{st}^a$) was then described by an intercept (${\beta }_0^a$), the climate, soil, and vegetation linear covariates (${\beta }_1^{a}X$) presented in Section 2.3, and a time-indexed, non-separable random field (${\mathrm{W}}_{st}^a$) shared with the pattern-level. As shown above, the ${\alpha }^1$ coefficient is an interaction term that modified the magnitude of the shared field to control for scale differences between the pattern (Gaussian) and process (Poisson) tiers.

The third- and fourth-tier model levels were structured similarly to the first and second but used to assess nymphs conditional on time-lagged estimates from the adult levels, such that,

$$\begin{array}{c}\hfill \mathrm{Nymphs}|\mathit{Adults}=\left\{\begin{array}{c}\mathrm{Pattern}\left\{\begin{array}{c}\log \left({\mathsf{\Lambda }}_{st}^n\right)={\mathrm{W}}_{st}^n\hfill \\ {\mathsf{\Lambda }}_{st}^n=\exp \{{\beta }_{\mathsf{\Lambda }}^n+{\mathrm{W}}_{st}^n\}\hfill \end{array}\right.\hfill \\ \\ \mathrm{Process}\left\{\begin{array}{c}{\mathrm{Y}}_{st}^n\sim \mathrm{Poisson}\left({\mu }_{st}^n\right)\hfill \\ \log \left({\mu }_{st}^n\right)={\beta }_0^n+{\beta }_1^{n}X+{\alpha }^2{\mathrm{W}}_{st}^n+{\alpha }^3{\mathrm{W}}_{st}^a\hfill \end{array}\right.\hfill \end{array}\right.\end{array}$$

where the nymph pattern (${\mathsf{\Lambda }}_{st}^n$) was jointly estimated with nymph counts ($Y_{st}^n$) as done for adults. However, in addition to the random field shared between the nymph pattern and process levels (${\mathrm{W}}_{st}^n$), the adult random field (${\mathrm{W}}_{st}^a$) was copied to the nymph process level to quantify the relationship between nymphs and adults estimated for the prior year. As described for adults, the ${\alpha }^2$ term enumerated interaction between the pattern and process levels for nymphs, while ${\alpha }^3$ accounted for spatiotemporal relationships between adults and nymphs. The same climate, soil, and vegetation covariates used in estimating adults were also evaluated with respect to nymphs, but with a set of dedicated coefficients (${\beta }_1^{n}X$) to facilitate comparison of coefficient effect sizes between life stages.

As a joint, spatiotemporal model with high-dimensionality, computation was performed on the USDA SCINet High-Performance Computing System (https://scinet.usda.gov/ (accessed on 20 December 2021)) using Integrated Laplace Approximation [76,77,78] and the PARDISO solver [79,80,81]. Spatiotemporal effects were specified with first-order autoregressive structure and weakly informative Penalizing Complexity priors [82,83]. Climate, soil, and vegetation linear effects were assigned vague zero mean normal priors with a 0.0001 precision.

2.5. Model Selection and Validation

To evaluate joint model performance, three models were constructed and then compared; individual models for each the adult (Model1) and nymph (Model2) life stage, and the full joint model detailed in Section 2.4 (Model3). Metrics used for model comparison included the deviance information criterion (DIC), Watanabe–Akaike information criterion (WAIC), and the log-conditional predictive ordinate (lCPO) [84,85]. The DIC and WAIC generally perform the same; however, the DIC sometimes under penalizes complexity, thus both were used for comparison. The WAIC is a within sample predictive score and is a fully Bayesian criterion [84,86]. The lCPO is a leave-one-out cross validation metric. Lower scores for all three measures suggest improved parsimony. Data were partitioned into training (80%) and testing (20%) sets for validation with equal proportions allocated from the adult and nymph stages.

2.6. Recruitment Rate Prediction, Consensus, and Forecasts

Outputs from the joint model described in Section 2.4 (Model3) were combined to calculate the Msang per-capita, nymph recruitment rate. Recruitment was defined comparably to a population per-capita growth rate [87,88], where recruitment was the natural log of current-year nymph abundance divided by adult abundance from the preceding year:

$$\begin{array}{c}\hfill \mathrm{Recruitment}=ln({\mathrm{Nymphs}}_t/{\mathrm{Adults}}_{t-1})\\ \hfill \end{array}$$

Formal prediction was conducted for all locations in the study domain (surveyed and un-surveyed localities) based on 30-year average climate conditions [50], soil properties, and vegetation before forecasting recruitment for the time period horizon 2021–2040 under the four SSP scenarios SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5. Because there is considerable uncertainty regarding GCM accuracy and climate forecasts broadly, a consensus approach was adopted in which prediction was independently performed based on the IPSL-CM6A-LR, CanESM5, and MIROC6 GCMs followed by (mean census) model averaging [89,90].

3. Results

Variable decomposition described in Section 2.3 resulted in three synthetic covariates summarizing climate, soil, and vegetation variation with respect to Msang occurrence (Figure 3). Component loadings for the top five contributing climate variables indicated that temperature was most important in describing Msang occurrence, with mean temperature of the warmest quarter contributing 33.7% to the synthetic covariate followed by temperature diurnal range (16.3%), temperature seasonality (13.3%), day-to-night temperatures oscillations (isothermality, 8.4%), and average temperature during the coldest quarter of the year (7.9%). Among the 285 evaluated soil attributes, humus rich Kastanozems, which are associated with grassland vegetation [91], showed the greatest contribution (5.8%), followed by Calcisols (4.6%), Ustolls (3.6%), Fluvents (3.6%), and Calcids (3.3%). MODIS derived remote sensing data suggested that pixel diversity calculated using the Simpson (61.6%) and Shannon (26.0%) indices accounted for the vast majority of variation in the synthetic vegetation covariate, with image texture metrics contributing slightly more. Synthetic covariates are mapped in Figure 3.

Comparison metrics indicated that the joint model (Model3) designed to concurrently estimate adult and nymph abundance exhibited improved parsimony over models constructed to individually assess the adult and nymph life stages (

Table 1. Model parsimony metrics. Deviance information criterion (DIC), Watanabe–Akaike information criterion (WAIC), and log-conditional predictive ordinate (lCPO) for all models. Codes listed in the Stage column signify models structured to estimate the adult and nymph life stages using individual (I) models or a joint (J) model. Note that the joint model (Model3) used to estimate adult and nymph abundance concurrently exhibited the best parsimony.

Model	DIC	WAIC	lCPO	Stage
Model1	10,322	15,271	1.28	Adult (I)
Model2	9905	15,757	1.20	Nymph (I)
Model3	10,134	14,977	1.13	Adult (J)
	9487	11,467	1.08	Nymph (J)

). In addition to parsimony measures, models for individual life stages differed in the importance and strength of estimated climate, soil, and vegetation covariates. For example, the vegetation covariate was judged to be important based on 95% Credible Intervals excluding zero in the individual, stage-specific models; however, the same vegetation covariate was found not to be important when life stages were modeled jointly (Figure 4).

Model estimated hyperparameters (

Table 2. Model hyperparameter summary. The table lists model estimated hyperparameters described in Section 2.4. First column provides effect name with columns to right giving the estimated mean, standard deviation (sd), and 95% credible intervals (2.5 Q and 97.5 Q). Note that ${\alpha }^3$ is the model estimated interaction coefficient that reflects spatiotemporal relationships between adults and nymphs. The coefficient value is negative indicating that as adult abundance increased, nymph abundance in the following year decreased (i.e., density-dependence). The ${\alpha }^3$ slope is graphed in Figure 5.

Hyperparameter	Mean	sd	2.5 Q	97.5 Q
${\mathrm{W}}_{st}^a$ Spatial Range	16.02	0.44	15.18	16.90
${\mathrm{W}}_{st}^a$ Stdev	0.43	0.01	0.41	0.44
${\mathrm{W}}_{st}^a$ Group $\rho$	0.98	<0.01	0.98	0.98
${\mathrm{W}}_{st}^n$ Spatial Range	17.67	0.57	16.58	18.83
${\mathrm{W}}_{st}^n$ Stdev	0.53	0.01	0.51	0.55
${\mathrm{W}}_{st}^n$ Group $\rho$	0.99	<0.01	0.99	0.99
${\alpha }^1$ (adults)	6.06	0.05	5.98	6.17
${\alpha }^2$ (nymphs)	10.39	0.19	10.12	10.83
${\alpha }^3$ (adults-nymphs)	−1.07	0.03	−1.14	−1.03

) revealed that spatial autocorrelation among adult sample locations fell to approximately zero at a distance of 16.0 km (${\mathrm{W}}_{st}^a$ Spatial Range). By comparison, the spatial range for nymph sample locations (${\mathrm{W}}_{st}^n$) was similar at about 17.7 km. Spatiotemporal autocorrelation between years was very high for both adults and nymphs with estimated correlation values exceeding $\rho \ge 0.98$ in both instances. Spatiotemporal effects shared between model pattern and process levels indicated strong, positive influence for both adults (${\alpha }^1$) and nymphs (${\alpha }^2$), suggesting that grasshopper abundances estimated for sample locations were considerably higher than numbers estimated for un-surveyed locations. The spatiotemporal effect shared between adult and nymph stages (${\alpha }^3$) was judged important based on credible intervals excluding zero. The negative polarity estimated for ${\alpha }^3$ signifies that on the average, as adult abundances increase, nymph abundances decrease between successive years. The ${\alpha }^3$ coefficient −1.07 [0.03 sd, (−1.14, −1.03) 95% CI] is plotted as functions of abundance and recruitment in Figure 5.

Model validation suggested excellent model performance (mean absolute error ≤ 10%) for predicted locations having ten or fewer surveyed abundances (adults or nymphs); however, under-prediction was apparent at locations with abundances greater than ten grasshoppers (Figure 6). Under-prediction was more pronounced for adult Msang but also noted for the nymph stage. Comparing predicted abundances to observed counts for all survey locations (testing and training data combined) suggested that under-prediction was especially evident for surveys during years with especially low or high counts (Figure 7).

Formal prediction for life stage-specific abundances and recruitment rates under historic climate conditions (30-year averages), soil characteristics, and vegetation are illustrated in Figure 8. Predictions reflect mean abundance and recruitment rates based on observations during the ten-year period 2011–2020. Forecasts for the twenty-year time horizon 2021–2040 are shown in Figure 9 for each of the SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 socioeconomic climate scenarios. Forecasts were made using climate, soil, and vegetation coefficients estimated by the full, joint model (Model3) and assume constant soil and vegetation conditions through the year 2040. Estimates were derived from GCM mean consensus.

4. Discussion

M. sanguinipes nymph recruitment rates are anticipated to exhibit strong geographic variation under projected climate change by the year 2040. In Wyoming, model results generally indicated that recruitment will remain approximately the same or slightly decrease across central and southern portions of the state (Figure 9), whereas grasshopper habitats immediately west and southeast of the Bridger–Tetons were forecast to experience increased recruitment. Recruitment within pastures, isolated grasslands, and alfalfa fields west of the Bridger–Tetons (e.g., Afton, Fairview, and Etna, WY) were revealed to show the greatest overall rate increases, but substantial rises were also projected for the northeast corner of Wyoming, in a region encompassing and to the north of the Thunder Basin National Grasslands. This basic pattern was consistent for all four SSP scenarios; however, SSPs derived for scenarios where climate mitigation is a low international priority (SSP3-7.0) and fossil fuel use is intensive (SSP5-8.5) showed markedly elevated recruitment increases in comparison to scenarios based on sustainable growth (SSP1-2.6) and intermediate action to reduce climate vulnerability (SSP1-2.6). Assuming that future Msang nymph survival rates remain approximately the same as historic averages, our principal conclusion is that alterations to prevailing temperature and precipitation regimes as instigated by climate change will amplify recruitment, thereby enlarging population sizes and intensifying agricultural pest impacts by 2040.

Of the nineteen climate variables that we evaluated, the top five contributing to Msang occurrence were connected to temperature variation, with mean temperature of the warmest quarter identified as the most important factor (Figure 3). Feeding by Msang nymphs is optimized at approximately 40 ${}^{\circ }$C, and feeding stops at temperatures above approximately 45 ${}^{\circ }$C or below 13 ${}^{\circ }$C [35]. In California, a region with a broad range of elevations like Wyoming, Msang populations differ in nymphal development rate such that they can complete development faster at higher elevations [92]. Age of adults at first reproduction is also earlier at higher elevations. Thus, a warming climate is likely to be favorable for development and population growth, particularly at higher elevations. Along comparable lines, more southerly Msang populations, like those found in some regions of Arizona, do not require diapause to hatch and the populations are multivoltine, resulting in two generations per year when and where the growing season is sufficiently long [33]. Should the non-diapause attribute move north with climate change, population growth and recruitment in WY, where populations are currently univoltine, could increase much more than what has been predicted here. This scenario would be consistent with other studies that have reported an increased frequency in upslope insect dispersal and upslope demographic shifts linked to climate change [93,94,95]. Additional research is needed to investigate how abiotic climate change might modify Msang development, population growth, and dispersal, and how biotic factors (e.g., interspecific competition, predation, infectious disease) might amplify or attenuate these effects. For example, although increased temperature might facilitate transition of current univoltine populations to a multivoltine pattern, any increases in population size may be offset by elevated disease incidence, because the Msang nymph cuticle is paler to minimize thermal elevation when nymphs develop at higher temperatures, but the adaptation also makes them more susceptible to an insect killing fungus [96].

Analysis of thirteen vegetation metrics suggested that plant diversity better described Msang occurrence than did other satellite-derived measures. Taken together, Simpson and Shannon diversity contributed more than 87% to the synthetic vegetation covariate used in modeling (Figure 3). The composition of plants at a location has been variably argued as both an important or minor indicator of grasshopper occurrence and abundance. In some instances, plant composition has been shown as an important indicator [97,98,99], but in other cases plant effects on generalist grasshopper species have been weak or mixed [100]. Interestingly, we found that the synthetic vegetation covariate was an important predictor of both nymph and adult abundance when each life stage was assessed independently; however, when the two life stages were modeled jointly, vegetation became insignificant (Figure 4). This may indicate that vegetation aids in predicting mere occurrence (habitat suitability), but is less useful in estimating abundance as the presence of another grasshopper (or life stage) better explains abundance in the target group than does plant diversity. A similar pattern was shown by climate and soil covariates in that the magnitude of climate and soil effects were greater when adults and nymphs were modeled separately, but substantially reduced (but still significant) when the two stages were assessed concurrently. Further research is needed to more fully partition the extrinsic environmental factors considered in this study from the intrinsic demographic factors that drive Msang population dynamics.

Msang abundance at a given time and location is a consequence of demographic processes tied to stage-specific environmental tolerances, interactions, and behaviors. From the analytical perspective, we argue that process-based methods, which allow for inclusion of underlying demographic mechanisms [101,102], offer advantages over standard species modeling techniques when extrapolating pest populations to future climate scenarios. As examples, our approach expanded on standard correlative distribution models to incorporate Msang demographic information, life stage density-dependence between successive years, and to account for preferentially sampled field observations. This modeling framework enabled us to detect that Msang nymphs exhibited stronger responses to climate, soil, and vegetative conditions than did adults, and that nymph–adult co-occurrence better explained abundances than did vegetation as represented by remote sensing data (Figure 4). Importantly, model structure also supported enumeration of density-dependence (Figure 5), a type of biotic interaction that can produce direct mortality or contribute to resource limitation that lowers female Msang fecundity and reproductive rates [16,24,103,104].

Finally, we view our model as a disciplined metaphor for reality [105] and our presented maps as geographic theories [106], both of which are based on assumptions but nonetheless complementary to empirical and experimental studies investigating pest insect dynamics. Beyond the uncertainly associated with remote sensing products [107,108], the presented study can be improved in at least two ways. Firstly, although model performance was reasonable for predicted abundances near the average rates observed at surveyed sites, nymph and adult numbers were underestimated at times and locations with especially high or low Msang population densities (Figure 6 and Figure 7). We interpret this result to imply that mechanisms driving population booms (high Msang densities) and busts (low Msang densities) were not fully captured by the explanatory variables included in our model. We suspect that the synthetic variables used in this analysis may have been too general to completely represent how climate change may impact key vegetation characteristics, like those connected to plant quality, composition, and senescence. Similarly, it may also be the case that the temporal or spatial resolution of the variables selected for model inclusion were insufficient to account for variation significantly above or below the mean abundance rates. For example, vegetation data incorporated into the current study did not vary through time as did climate variables and were too coarse to capture intra-annual or transient seasonal dynamics like those linked to the El Niño/Southern Oscillation (ENSO). Given these issues, studies that incorporate seasonal variability or explicitly model how vegetation composition patterns may change under future climate conditions might better portray between-year variation in grasshopper abundances. Drawing meaningful, mechanistic inference from models is often challenging [54,55] and additional work is needed to identify other environmental characteristics or intrinsic population attributes that contribute to Msang boom and bust cycles. Secondly, though our study was conducted at the landscape scale, the Msang species range includes the majority of North America. The presented model would benefit from being scaled up to create a more comprehensive analysis of Msang population and habitat characteristics across the entirety of its distributional range.

5. Summary and Conclusions

We applied ten years of field surveys to assess the response of Msang nymph recruitment to projected climate conditions at the landscape scale. Central to the study were methodological advances that enabled grasshoppers to be appraised using a multi-level, joint modeling framework that separately quantified nymph and adult life stages while concurrently incorporating density-dependence and accounting for observation bias connected to preferential sampling. As a new method, our hierarchical modeling approach may be readily adapted to other species, particularly when there is a need to account for covariation between different life stages, density-dependence across space and through time, or biotic interactions among different species. Climate change is expected to alter prevailing temperature and precipitation patterns this century and our findings indicated that Msang nymph recruitment rates will exhibit strong geographic variation in concert with this climate change. Results indicated that Msang population sizes at many locations will remain comparable to those historically observed; however, alterations to prevailing temperature and precipitation regimes as instigated by climate change are anticipated to amplify recruitment at several locations, potentially enlarging population sizes and intensifying agricultural pest impacts by 2040.