Geographic Variation in Migratory Grasshopper Recruitment under Projected Climate Change

: 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 ) ﬁeld 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 ﬁndings 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.


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

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 • and 45 • North Latitude and −104 • and −111 • 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.

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,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.

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].

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 (Λ a,n st ), was concurrently estimated with adult (Y a st ) and nymph (Y n st ) 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: where the logarithm of intensity (Λ a,n st ) for grasshoppers in the adult (superscript a) and nymph (superscript n) stages at each geographic location s (s = 1, 2, 3, . . . , n) and year t (t = 2011, 2012, 2013, . . . , 2020) were estimated as Gaussian random fields (W a,n st ). Following Lindgren et al. [67] and Simpson et al. [68], the matrices Q(κ, τ) used to define the spatial fields (W a,n st ) were approximated using stochastic partial differential equations (SPDE) that facilitated implicit estimation of spatial range (κ) and scale (τ) 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, where adult point pattern intensity (Λ a st ) was approximated as the exponential of a Gaussian random field (W a st ) and an intercept (β a Λ ) in the first tier, and the second tier estimated adult counts (Y a st ) following a Poisson distribution. The log of mean adult abundance (µ a st ) was then described by an intercept (β a 0 ), the climate, soil, and vegetation linear covariates (β a 1 X) presented in Section 2.3, and a time-indexed, non-separable random field (W a st ) shared with the pattern-level. As shown above, the α 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, where the nymph pattern (Λ n st ) was jointly estimated with nymph counts (Y n st ) as done for adults. However, in addition to the random field shared between the nymph pattern and process levels (W n st ), the adult random field (W a st ) 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 α 2 term enumerated interaction between the pattern and process levels for nymphs, while α 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 (β n 1 X) to facilitate comparison of coefficient effect sizes between life stages.

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.

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: Recruitment = ln(Nymphs t /Adults t−1 ) 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].

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). 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) revealed that spatial autocorrelation among adult sample locations fell to approximately zero at a distance of 16.0 km (W a st Spatial Range). By comparison, the spatial range for nymph sample locations (W n st ) was similar at about 17.7 km. Spatiotemporal autocorrelation between years was very high for both adults and nymphs with estimated correlation values exceeding ρ ≥ 0.98 in both instances. Spatiotemporal effects shared between model pattern and process levels indicated strong, positive influence for both adults (α 1 ) and nymphs (α 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 (α 3 ) was judged important based on credible intervals excluding zero. The negative polarity estimated for α 3 signifies that on the average, as adult abundances increase, nymph abundances decrease between successive years. The α 3 coefficient −1.07[0.03 sd, (−1.14, −1.03) 95% CI] is plotted as functions of abundance and recruitment in Figure 5.  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 α 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 α 3 slope is graphed in Figure 5.    Table 2).

Hyperparameter
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). Figure 6. Model performance summary. The graph depicts the relationship of predicted grasshopper abundance (vertical axis) to observed grasshopper counts collected through field survey (horizontal axis). Both axes have been transformed to the log(n + 1) scale to better illustrate value ranges. Points (open circles) plotted to graph symbolize 20% of data randomly selected for model validation with adult samples (black color) distinguished from those for nymphs (red color). The diagonal line bisecting the plot area from the origin signifies perfect prediction (predicted = observed) with the shaded area below the diagonal representing under-prediction and the area above the diagonal indicating over-prediction. Smooth, curvilinear lines illustrate general trend by adult (black color) and nymph (red color) stages. Note that under-prediction is apparent for both adults and nymphs (curvilinear trend lines below diagonal).
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.    to (right). The top row displays projected recruitment rates for each SSP and the bottom row shows the difference between projected rates and historic mean recruitment as illustrated in Figure 8. Red colors in the bottom row signify locations forecast to exhibit increased recruitment rates, whereas locations in blue indicate locations projected to have decreased recruitment relative to historic means. Projections are based on consensus mean averages from three GCMs. Wyoming county boundaries have been overlain as a geographic reference to facilitate comparison between maps. Note that the darkest tones (dark blue colors) in the second row are spatially coincident with regions of high elevation in Figure 1 (low-quality Msang habitat).

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 • C, and feeding stops at temperatures above approximately 45 • C or below 13 • 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 (Figures 6 and 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.

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.