Phenology dictates the impact of climate change on geographic distributions of six co‐occurring North American grasshoppers

Abstract Throughout the last century, climate change has altered the geographic distributions of many species. Insects, in particular, vary in their ability to track changing climates, and it is likely that phenology is an important determinant of how well insects can either expand or shift their geographic distributions in response to climate change. Grasshoppers are an ideal group to test the hypothesis that phenology correlates with range expansion, given that co‐occurring confamilial, and even congeneric, species can differ in phenology. Here, I tested the hypothesis that early‐ and late‐season species should possess different range expansion potentials, as estimated by habitat suitability from ecological niche models. I used nine different modeling techniques to estimate habitat suitability of six grasshopper species of varying phenology under two climate scenarios for the year 2050. My results suggest that, of the six species examined here, early‐season species were more sensitive to climate change than late‐season species. The three early‐season species examined here might shift northward during the spring, while the modeled geographic distributions of the three late‐season species were generally constant under climate change, likely because they were pre‐adapted to hot and dry conditions. Phenology might therefore be a good predictor of how insect distributions might change in the future, but this hypothesis remains to be tested at a broader scale.

latitudinal shifts are quite variable among species (Beckmann et al., 2015;Chen et al., 2011). There is, as yet, no consistent pattern that explains which insect species exhibit range shifts and which do not.
Life history strategy is often invoked as a determinant of potential for range shifts (Estrada et al., 2015), but there are few phylogenetically controlled studies that contrast different life history strategies within a single clade. Orthopterans (grasshoppers, crickets, and katydids) provide an opportunity to compare sensitivity to climate change among life history strategies, given that co-occurring grasshopper species possess a remarkable functional diversity McClenaghan et al., 2015). In the United Kingdom, for example, warm-adapted, generalist grasshoppers with high dispersal ability are the only species to have undergone range expansion (Beckmann et al., 2015). In the Great Plains of North America, grasshoppers can be broadly partitioned into two life history groups. Early-emerging species, such as Arphia conspersa, Eritettix simplex, and Xanthippus corallipes, overwinter as nymphs and emerge as adults in the spring (Capinera & Sechrist, 1982). These three species reach peak abundance in April or May, several months before most other grasshopper species (Buckley et al., 2021). Lateemerging species, such as Arphia pseudonietana, Opeia obscura, and Phoetaliotes nebrascensis, overwinter as eggs, hatch in early summer, and reach the adult stage by mid-to-late summer in July or August (Branson, 2016;Capinera & Sechrist, 1982). Given their different climatological niches, that is, cold wet spring vs. dry hot summer, and the fact that two of these species are congeners, these six species provide a phylogenetically controlled experiment for how life history might impact how species respond to climate change.
It is possible to compare how species of different life history strategies might respond to climate change using ecological niche models (ENMs). ENMs correlate occurrence records with climate and are often used to predict range expansions. For example, ENMs can identify areas at risk of invasion under future climates (Gong et al., 2020;Kistner-Thomas, 2019) or identify high-priority conservation targets (Garzon et al., 2021), which is critically important as the ranges of many threatened species might collapse in the near future (Lemoine, 2015). One shortcoming is that ENMs rarely account for phenology (Ingenloff & Peterson, 2021); many simply use mean annual temperature or precipitation (Booth et al., 2014;Title & Bemmels, 2018). The use of models that incorporate annual trends might miss spatiotemporal shifts in distribution because habitat suitability can change within a given year (Martinez-Meyer et al., 2004).
For example, an annual model predicted that suitable habitat for the mosquito Aedes aegypti should cover most of Mexico, while a temporally explicit model revealed distinct seasonal shifts in habitat suitability of A. aegypti (Peterson et al., 2005). This is also true for insects in seasonal temperate environments. The shortgrass steppe of Colorado, on average, is 8.7°C and receives 395 mm of rainfall per year. Yet early-season grasshoppers that emerge in May experience an environment that is 22°C and receives 61 mm of rain. Late-season grasshoppers, in contrast, emerge into an arid environment of 29°C and 40 mm of rainfall. Thus, accurate predictions in ENMs require that climatological data match life history data as closely as possible (Ingenloff & Peterson, 2021). Using mean annual temperature or precipitation might over-or underestimate the sensitivity of species to climate change by mischaracterizing their environmental niches.
Here, I modeled how climate change might affect intra-annual spatiotemporal patterns of habitat suitability for six North American grasshoppers that differ in life history strategy. Specifically, I predicted that three early-emerging species, A. conspersa, E. simplex, and X. corallippes, would favor cool, wetter temperatures. Thus, both the southern and northern boundaries of suitable habitat conditions for these three species should move northward (i.e., total range shift) and occur earlier in the year, which would predict an advancing phenology. In contrast, A. pseudonietana, O. obscura, and P. nebrascensis all emerge as adults in July and August, and therefore should have suitable habitat expand northward while maintaining the current southern boundary (i.e., range expansion), and suitable habitat should extend later into the year. To test these hypotheses, I constructed ENMs using nine separate machine learning classification techniques and predicted suitable habitat into the future for two different climate scenarios, with four general circulation models for each climate scenario used to produce an ensemble prediction.

| Environmental data
To construct climatic niches for each species, I downloaded WorldClim2 climate data (Fick & Hijmans, 2017), which is an interpolated climate dataset covering the years 1970-2000. As I was specifically examining phenological patterns, I used monthly data at a 5 arc-minute resolution. The use of monthly data restricted the environmental variables to average monthly precipitation, average monthly temperature, minimum monthly temperature, and maximum monthly temperature, as other WorldClim2 variables are either seasonal aggregates or unavailable at monthly time steps. Given the extremely high correlation among temperature variables (r > .90 for all temperature combinations, Table S1), I used only mean monthly precipitation (mPPT) and mean minimum monthly temperature (mT min ) for all subsequent analyses.

| Data cleaning, filtering, and pseudoabsences
I cleaned GBIF records following a standard pipeline (Feng et al., 2019;Zurell et al., 2020). First, I dropped any records with null values for latitude, longitude, month, or year. Next, I removed records with a "0" for latitude or longitude. I then dropped any observations that had coordinates identical to those of a US state capital city to within 0.01 decimal degrees, and also dropped any duplicate geographic coordinates except for those observations in different months and years. Once this pre-screening was complete, I visually I then filtered data to remove pseudoreplicates in environmental space. Although many studies advocate spatial filtering, I instead filtered observations on the basis of environmental similarity. Such environmental filtering has shown to be more robust, less biased, and more accurate than spatial filtering (Varela et al., 2014). For the en- Due to the temporal aspect of the hypotheses tested here, I used a phenological approach to generating pseudoabsences (Ingenloff & Peterson, 2021). Briefly, for each species, I calculated the number of observations falling within each month. I then generated the same number of pseudoabsences from the mPPT and mT min for that month. The end product was the same number of observations and pseudoabsences for each species within each month. I chose to use equal numbers of pseudoabsences because a 1:1 ratio of observations:pseudoabsences performs the best for many classification models (Barbet-Massin et al., 2012). I used a simple random pattern, rather than a gridded or weighted approach, because multiple studies demonstrated that simple random pseudoabsences perform at least as well as weighted or stratified pseudoabsences, especially for some of the classification methods used here (Barbet-Massin et al., 2012;Hanberry et al., 2012).

| Ecological niche models
ENMs use correlative approaches to summarize the climatic niche of a species. There is a large degree of uncertainty in ENMs, including uncertainty due to presence-only sampling, spatial biases, and in climate models. Perhaps the largest source of uncertainty is among modeling techniques (Araújo et al., 2005). Different methods make which did not include an interaction between mPPT and mT min .
In this model, y is presence/pseudoabsence (1/0) and z is the log odds of occurrence (i.e., logit transformation).

K-Neighbors Classifier (KNC):
A KNC uses a simple "vote-counting" method to assign a point to a class. Essentially, an unknown point (test data) is mapped into environmental space with training data. The algorithm counts the n nearest neighbors and assigns the test point to the class with the majority or plurality of neighbors. The output can be converted into a probability by counting the fraction of n points belonging to a given class. For the model here, I used n = 5 equally weighted neighboring points, and the distances between training points and the test points in environmental space were determined via Euclidean distance.

Gaussian Process Classifier (GPC):
Gaussian process models treat data as arriving from a multivariate distribution, generated by an unknown function: is the function describing the variability of x in space, m(x) is the mean function, and K(x, x′) is the kernel/covariance function. Because the kernels allow for covariance among observations that varies with the distance of observations, continuous Gaussian process models are popular for time series and spatial modeling, where they are known as "kriging" (Brahim-Belhouari & Bermak, 2004;Roberts et al., 2013). GPCs extend Gaussian process models to a binomial response using latent variables, much like logistic regression: where z(x) is a latent variable achieved by the logistic transformation of pseudoabsence (0) and presence (1) data. In practice, we often assume a constant mean: such that the kernel choice dictates the shape of the function.
Researchers have advocated GPCs for ENMs because they are often more accurate than other classification methods, such as boosted regression trees, generalized additive models, and generalized linear models (Golding & Purse, 2016). Here, I constructed ENMs from GPCs using the radial basis function: where α is a scaling parameter determining the magnitude of process noise and l is a length parameter that determines the smoothness of the function.

Decision Tree Classifier (DTC):
DTCs are nonparametric, supervised machine learning techniques that construct decision trees using if/then rules from training data in order to infer the class of the test points. Essentially, decision trees split the data into groups then conduct logistic regressions to classify the training data. The split with the highest predictive ability is taken as the first decision criteria to generate two new groupings within the next level of the tree. The procedure proceeds iteratively within each grouping until a maximum tree depth is achieved.
These models are simple, fast, and nonlinear, but can be prone to overfitting, particularly if a tree is too deep. For the model here, I used the Gini criteria to evaluate the quality of a given split, with a maximum tree depth of five levels. I required each group to have a minimum of two samples.

Random Forest Classifier (RFC):
An RFC is a "meta"-classifier that constructs a number of DTCs from random subsamples of the training data and averages the outputs. For the RFC here, I generated 100 random DTCs using the same Gini criteria.

Artificial Neural Network (ANN):
ANNs with multilevel perceptrons approximate the way human brains process information by allowing computing nodes, called neurons, to process and share information to inform an output. ANNs consist of three layers of nodes: input nodes, hidden process nodes, and output nodes.
The input layer contains nodes for each feature (i.e., explanatory variable), hidden process nodes combine features with a weighted linear function, and an output function uses a nonlinear function to transform the hidden process nodes into a binary or continuous response. Several authors have advocated using ANNs for ENMs (Maravelias et al., 2003), in particular because they outperform many other methods for constructing ENMs, such as classification trees, generalized linear models, generalized additive models, and spatial interpolators (Segurado & Araújo, 2004). I trained the linear weights using a stochastic gradient optimizer, and the nodes were translated into a real output using the rectified linear unit function max(0, x). The ANN here had one hidden layer with 100 nodes, and a regularization parameter α = 1.

Ada Boost Classifier (ABC):
The ABC is similar to RFC, in that it relies on multiple DTCs. However, whereas RFCs generate 100 random DTCs and then average the outputs, ABCs proceed iteratively, repeatedly fitting the same DTC on the training data but with the weights of incorrect cases adjusted so the classifier focuses on more difficult cases. I used the SAMME.R algorithm, stopping at a maximum of 50 iterations. The SAMME.R algorithm is an updating algorithm that uses the probabilities of belonging to each class for each point as weights in an exponential loss function used to assess model fit (Pedregosa et al., 2011).

Naïve Bayesian Classifier (NBC):
NBCs are simple classifiers based on Bayes' rule. Bayes' rule can calculate the probability that a given map pixel should belong to a class k (i.e., present/absent) as: where x is the environmental variable, p(k|x) is the probability that a pixel of a given environment x belongs to class k, p(k) is the prior probability of belonging to class k, and p(x) is the probability of the environmental variable occurring in the model. For example, imagine classifying whether a pixel should be suitable habitat for a bird (k = present), depending on whether it is forested or not. In this case, p(x|k) is the probability that a pixel is forest given that a bird is present, or the proportion of times a bird was observed in forests, p(k) is the proportion of sightings of the bird throughout the entire dataset, and p(x) is the proportion of pixels that are forested. This example has a discrete predictor, but Gaussian NBCs extend classification to continuous predictors, such as temperature, by using the Gaussian density distribution to calculate the likelihood of a given temperature given an observation of present or absent: In this case, the probability of a bird being present at a given temperature is This method can be extended to multiple predictors by: Gaussian NBCs, along with the other methods here, can be used as a classification algorithm to model species niches (Guo & Liu, 2010). A drawback of this method is that it assumes independence of the features, but it has been shown to be an accurate method for constructing ENMs (Guo & Liu, 2010 Prior to analyses, both mPPT and mT min were standardized to N(0, 1) distributions to improve model fitting. Data were then split into training and test groups containing 66% and 33% of the data, respectively.
Data were split in a stratified manner to ensure equal proportions of presences/pseudoabsences in both the training and test data. Models were fit to the training data, and then tested for goodness of fit on the test data using the area under receiver operating characteristic curves (AUC-ROC). AUC-ROC scores for each of the nine models were then averaged to produce an "ensemble AUC-ROC" (Araújo et al., 2005).
For every species, I projected the current distribution throughout every month of the year based on WorldClim2 monthly data for mPPT and mT min at 5 arc-minute resolution (Fick & Hijmans, 2017 (Peters & Hausfather, 2020). However, the RCP 8.5 scenario is still useful as a "worst-case" baseline. When reconstructing these climate niches, modeling algorithms varied in their performance, although models performed similarly within a species (Table 1). That is, models within a species produced similar AUC-ROC scores (SD < 0.05), with the exception of O. obscura, where GPCs, NBCs, and QDAs performed exceptionally well p (k|x) ∝ p (x|k) p (k) p k|x 1 , x 2 , …, x n ∝ p x 1 |k p x 2 |k …p x n |k p (k)    F I G U R E 3 Predicted, current distribution of the early-season species Eritettix simplex in March, April, and May throughout the Great Plains of North America under current conditions, RCP 4.5, and RCP 8.5. Predictions are the ensemble/stacked averages from the nine different classifiers. The color palette was chosen so that regions where absence is more likely than presence (probability of occurrence <0.5) are shaded in blue, while regions where presence is more likely than absence (probability of occurrence >0.5) are shaded in reds.
Regions where presence and absence are equiprobable (probability of occurrence ~0.5) are shaded in whites/greys. Panels for RCP 4.5 and RCP 8.5 show the change in habitat suitability across each month, with orange regions denoting an increase in habitat suitability, and purple regions denoting a decrease in habitat suitability. Raw probabilities for each climate scenario are given in the supplemental figures suitable habitat for O. obscura and P. nebrascensis increased by <20% in most months, while A. pseudonietana showed evidence for range collapse under RCP 4.5 ( Figure 8). As described above, much of the increase in suitable habitat for P. nebrascensis was a longitudinal expansion, rather than a latitudinal shift (Figure 7).

| DISCUSS ION
As climate change alters the fundamental abiotic template of most ecosystems, many species are tracking favorable climates northward or to higher elevations. Yet, species vary in their ability to follow suitable climates (Beckmann et al., 2015;Chen et al., 2011). While life history characteristics like dispersal undoubtedly play a role in the capacity for range expansion (Beckmann et al., 2015), I hypothesized that the three early-season species examined here would shift poleward while late-season species would have relatively stable geographic distributions due to shifts in habitat suitability in future climates. In testing these hypotheses, I was able to partially confirm my hypotheses. The three early-season species exhibited range expansions via a poleward shift of the northern range limit while maintaining southern range limits, while the three late-season species F I G U R E 4 Predicted, current distribution of the early-season species Xanthippus corallipes in March, April, and May throughout the Great Plains of North America under current conditions, RCP 4.5, and RCP 8.5. Predictions are the ensemble/stacked averages from the nine different classifiers. The color palette was chosen so that regions where absence is more likely than presence (probability of occurrence <0.5) are shaded in blue, while regions where presence is more likely than absence (probability of occurrence >0.5) are shaded in reds.
Regions where presence and absence are equiprobable (probability of occurrence ~0.5) are shaded in whites/greys. Panels for RCP 4.5 and RCP 8.5 show the change in habitat suitability across each month, with orange regions denoting an increase in habitat suitability, and purple regions denoting a decrease in habitat suitability. Raw probabilities for each climate scenario are given in the supplemental figures appeared largely unaffected by climate change. Thus, it appears that spatially co-occurring species might exhibit different responses to climate change based on phenology, and my work highlights the need to account for emergence phenology in species distribution modeling.
Phenological shifts are common responses to climate change for both plants and insects. In plants, warming often advances emergence and flowering dates (Price & Waser, 1998;Wolkovich et al., 2012). However, not all plant species advance their phenology with warming; the phenological response to warming appears to largely depend on plant life history. Spring species that flower early often advance their phenology, sometimes by several weeks, while species that flower in fall can delay their phenology (Sherry et al., 2007).
Like plants, insects also advance their emergence dates (Ellwood et al., 2012). Yet, few studies have tested whether phenology (i.e., overwintering state) influences how insects alter geographic distributions in response to climate change (but see Poyry et al., 2009).
Grasshoppers are an ideal system to test for such possibilities because co-occurring species, indeed even co-occurring congeners as in the case of Arphia, possess early and late phenologies (Capinera & Sechrist, 1982), providing the opportunity for phylogenetically controlled tests of range expansion.

F I G U R E 5
Predicted, current distribution of the early-season species Arphia pseudonietana in March, April, and May throughout the Great Plains of North America under current conditions, RCP 4.5, and RCP 8.5. Predictions are the ensemble/stacked averages from the nine different classifiers. The color palette was chosen so that regions where absence is more likely than presence (probability of occurrence <0.5) are shaded in blue, while regions where presence is more likely than absence (probability of occurrence >0.5) are shaded in reds.
Regions where presence and absence are equiprobable (probability of occurrence ~0.5) are shaded in whites/greys. Panels for RCP 4.5 and RCP 8.5 show the change in habitat suitability across each month, with orange regions denoting an increase in habitat suitability, and purple regions denoting a decrease in habitat suitability. Raw probabilities for each climate scenario are given in the supplemental figures My study suggests that, as with plants, insect emergence phenology might be an important predictor of how insects respond to climate change. In this study, three early-season grasshopper species did not advance their phenology across their entire range, but only in the northernmost regions of the Great Plains. Viewed spatially, this pattern amounts to a northern range expansion in early spring, and viewed temporally, it amounts to an advanced phenology in northern areas. However, late-season species that share the same geographic distribution as early-season species might shift their distributions less in response to climate change. This is likely because climate change will make much of North America both warmer and drier (Greve et al., 2014;Sheffield & Wood, 2008), an environment to which late-season grasshoppers are already adapted. Importantly, no species here showed a range collapse; all grasshopper species examined here are predicted to maintain, if not expand, their current range size. This matches predictions of many other insects (Au & Bonebrake, 2019;de la Giroday et al., 2012;Wilson et al., 2021), and suggests that climate change might not directly precipitate the decline in insect abundances.
It is also likely that grasshopper populations vary in their sensitivity to climate change depending on latitude. Latitudinal temperature variation plays a strong role in grasshopper life history, with F I G U R E 6 Predicted, current distribution of the early-season species Opeia obscura in March, April, and May throughout the Great Plains of North America under current conditions, RCP 4.5, and RCP 8.5. Predictions are the ensemble/stacked averages from the nine different classifiers. The color palette was chosen so that regions where absence is more likely than presence (probability of occurrence <0.5) are shaded in blue, while regions where presence is more likely than absence (probability of occurrence >0.5) are shaded in reds. Regions where presence and absence are equiprobable (probability of occurrence ~0.5) are shaded in whites/greys. Panels for RCP 4.5 and RCP 8.5 show the change in habitat suitability across each month, with orange regions denoting an increase in habitat suitability, and purple regions denoting a decrease in habitat suitability. Raw probabilities for each climate scenario are given in the supplemental figures northern populations being often being smaller, developing slower, and having lower egg viability than warmer, southern populations (Telfer & Hassall, 1999). There is little information on how grasshopper phenology varies with latitude, but numerous studies have documented delayed eclosion/emergence dates for grasshoppers in cold, high-elevation populations (Buckley et al., 2021;Nufio & Buckley, 2019). Northern populations might therefore have delayed phenologies due to a slower accumulation of growing degree days, and these populations might be primed to advance their phenologies more strongly than southern populations where climate warming will be weaker and growing degree days already accumulate rapidly (Diffenbaugh & Giorgi, 2012 (Hao & Kang, 2004), whereas other species can survive cold temperatures only in the presence of an insulating snowpack (Riegart, 1967). Still other species require eggs to be cold stratified for hatching success (Fisher et al., 1996). For nymphal wintering species, the termination of diapause can be trigger by either photoperiod (Ingrisch, 1984) or accumulation of growing degree days (Fisher et al., 1996). Thus, the large variation in how grasshopper species will respond to changes in winter temperature cannot be captured by a general, correlative model. The only means for an ENM to capture such variation is to include winter temperature; however, winter temperatures are highly correlated with spring temperatures (r > .90, see Table S2). Thus, winter temperatures are redundant with spring temperatures and add relatively little new information to the model.
Projecting species distributions in future climates remains an important avenue of research. Doing so can inform us of habitat potentially at risk from species invasions (Gong et al., 2020;Kistner-Thomas, 2019), identify species at risk of collapse (Lemoine, 2015), and pinpoint regions of high priority for conservation (Garzon et al., 2021). In doing so, researchers must carefully account for source of uncertainty. In this study, I accounted for model uncertainty by using nine different ENM estimation techniques, for projection uncertainty by using four separate GCMs, and for scenario uncertainty by using RCP 4.5 and 8.5. My results suggest that phenology might be a good predictor of how insect distributions might change in the future. For North American grasshoppers, early-season species from cool environments might expand their northern range extent, while late-season species that are already adapted to hot and dry conditions could experience only modest changes in geographic distribution. My results provide tantalizing evidence that phenology could explain a considerable amount of variation in insect species' ability to respond to climate change.
F I G U R E 8 Percent change in suitable habitat area for each species under two climate scenarios. Suitable habitat area was calculated as the number of grid cells where the probability of occurrence was greater than 50%

ACK N OWLED G EM ENTS
I would like to acknowledge helpful comments by several anonymous reviewers. This manuscript was made possible by NSF DEB 1941390.

CO N FLI C T S O F I NTE R E S T
The author declares no conflicts of interest.