Do different rates of gene flow underlie variation in phenotypic and phenological clines in a montane grasshopper community?

Abstract Species responses to environmental change are likely to depend on existing genetic and phenotypic variation, as well as evolutionary potential. A key challenge is to determine whether gene flow might facilitate or impede genomic divergence among populations responding to environmental change, and if emergent phenotypic variation is dependent on gene flow rates. A general expectation is that patterns of genetic differentiation in a set of codistributed species reflect differences in dispersal ability. In less dispersive species, we predict greater genetic divergence and reduced gene flow. This could lead to covariation in life‐history traits due to local adaptation, although plasticity or drift could mirror these patterns. We compare genome‐wide patterns of genetic structure in four phenotypically variable grasshopper species along a steep elevation gradient near Boulder, Colorado, and test the hypothesis that genomic differentiation is greater in short‐winged grasshopper species, and statistically associated with variation in growth, reproductive, and physiological traits along this gradient. In addition, we estimate rates of gene flow under competing demographic models, as well as potential gene flow through surveys of phenological overlap among populations within a species. All species exhibit genetic structure along the elevation gradient and limited gene flow. The most pronounced genetic divergence appears in short‐winged (less dispersive) species, which also exhibit less phenological overlap among populations. A high‐elevation population of the most widespread species, Melanoplus sanguinipes, appears to be a sink population derived from low elevation populations. While dispersal ability has a clear connection to the genetic structure in different species, genetic distance does not predict growth, reproductive, or physiological trait variation in any species, requiring further investigation to clearly link phenotypic divergence to local adaptation.


| INTRODUC TI ON
There is increasing recognition that species' responses to environmental change depend on their potential for plastic and adaptive responses (Hoffmann & Sgrö, 2011), and yet the evolutionary potential of populations remains difficult to determine due to the complex genetic architecture and quantitative nature of key life-history and physiological traits (Lande, 1982). Steep environmental gradients are a promising context in which to study the evolution of these phenotypic traits, as environmental selection is strong and genetic differences among populations often reflect the forces of selection rather than long periods of genetic isolation (Savolainen, Lascoux, & Merilä, 2013). However, it is well known that life-history and physiological phenotypes manifest both from plasticity in response to environmental conditions (Angilletta, Steury, & Sears, 2004;Buckley, Nufio, Kirk, & Kingsolver, 2015;Nylin & Gotthard, 1998) and local adaptation associated with genetic differentiation (Dahlhoff et al., 2008;Natarajan et al., 2015). Understanding how gene flow facilitates or impedes genomic divergence among populations, and whether emergent phenotypic variation in life-history and physiological traits is limited by gene flow rates, remains an important challenge for predicting how populations respond to environmental change.
One classic example of environmentally driven phenotypic clines involves insect species found along elevational gradients, where higher elevation populations exhibit shorter wings, concomitant with a reduction in flight capacity and dispersal distance (Wagner & Liebherr, 1992). Temperature and the length of the active season decrease rapidly with elevation, and these two factors are particularly important for insects, in which rates of growth and development, digestive efficiency, and mortality are tightly linked to the thermal environment (Deutsch et al., 2008;Frazier, Huey, & Berrigan, 2006;Scriber & Lederhouse, 1992). Higher elevation insect populations are often larger than their lower elevation counterparts due to slower development in cooler environments, especially when wings are reduced in size (Angilletta et al., 2004;McCulloch & Waters, 2018). Alternatively, some species (especially fully winged species) reverse this pattern, becoming smaller in order to rapidly complete development during a limited active season (Dearn, 1977;Klok & Harrison, 2013;Shelomi, 2012). A decline in fecundity with elevation is common and is attributed to a decline in body size (for some species) and/or departure from physiological optima (Hodkinson, 2005). Thus, there often appear to be fitness trade-offs, yet it remains unclear whether these trait differences are due to adaptive phenotypic plasticity or local adaptation. Despite extensive work on a variety of species along elevation gradients, there has been limited investigation of underlying gene flow patterns. Broad surveys of plant and animal studies have shown frequent genetic and phenotypic divergence with altitude (Keller, Alexander, Holderegger, & Edwards, 2013;Ohsawa & Ide, 2008). For example, several studies have shown strong differences in allele frequencies of metabolic enzymes in insect populations (Orr, Porter, Mousseau, & Dingle, 1994;Rank & Dahlhoff, 2002). However, more effort is needed to integrate explicit measurements of gene flow with these well-studied elevational patterns, in order to understand directionality and rates of gene flow.
The genetic basis of life-history and physiological divergence has important implications for climate change responses, as phenotypic plasticity allows individuals to track local environmental optima and variation over short timescales, but may ultimately slow evolutionary responses to more pronounced environmental change (Ghalambor et al., 2015;Ghalambor, McKay, Carroll, & Reznick, 2007;Price, Qvarnström, & Irwin, 2003). Local adaptation to different environments could speed up evolutionary responses to climate change and potentially allow greater phenotypic change, but is unlikely to have evolved if populations have high gene flow (Bridle & Vines, 2007;Kremer et al., 2012;Lenormand, 2002;Moritz et al., 2012). At present, it remains difficult and costly to determine genetic adaptation from large insect genomes. Fortunately, patterns of genetic structure and gene flow at neutral markers can be informative about underlying evolutionary processes (Crispo, 2008). As a general expectation, if phenotypic variation along an elevational gradient is driven by local adaptation, as opposed to phenotypic plasticity, one would expect to see evidence of genetic differentiation and limited gene flow, aligned with patterns of phenotypic variation (Kawecki & Ebert, 2004;Sultan & Spencer, 2002). While phenotypic variation might exist simply as a result of genetic drift, it is not expected to covary systematically in relation to an environmental gradient. Reducing gene flow over strong environmental gradients is often required for local adaptation to emerge (Cassel-Lundhagen, Kaňuch, Low, & Berggren, 2011). This assumption about gene flow generally holds, except when phenotypic traits are controlled by a small number of genes of large effect (Pfeifer et al., 2018;Tigano & Friesen, 2016), as these phenotypes could be maintained by divergent selection despite strong gene flow.
Here, we test the hypothesis that clinal differences in growth and reproduction "(life-history traits)" and physiological traits of four ecologically divergent grasshopper species Buckley, Nufio, & Kingsolver, 2014;Levy & Nufio, 2015) are associated with genetic differentiation and reduced demographic estimates of gene flow (local adaptation). We examine a large set of genome-wide single nucleotide polymorphisms for each species to estimate genetic parameters. Our study focuses on a steep environmental gradient near Boulder, Colorado, spanning 2,000 m elevation over a horizontal distance of ~35 km, where Melanoplus boulderen sis and Aeropedellus clavatus, short-winged, early-season species, can be compared to Camnula pellucida and Melanoplus sanguinipes, long-winged, mid to late-season species. Short-winged and flightless grasshopper species should have reduced dispersal potential, show greater genetic differentiation among geographical populations and a greater potential for local adaptation than long-winged species (Knowles, 2000;Ortego, García-Navas, Noguerales, & Cordero, 2015;Tinnert & Forsman, 2017). The two long-winged species are considered more dispersive, because in addition to occurring at a broader range in elevation, they are occasionally collected as "accidentals" (nonresident migrants) at high-elevation sites along the gradient (Alexander, 1964). Phenotypic clines in morphological, physiological, and reproductive traits have been well documented among these species (Table 1, Figure 1), with notable differences among our four focal taxa corresponding to ecological divergence Levy & Nufio, 2015).
While phenotypic variation is evident in all four grasshopper species, it is unclear how that variation is related to geographical, environmental, and genetic factors. We first compare population genetic structure and test for an association between genetic divergence and observed phenotypic clines. Ideally, such a comparison would be based on the additive genetic variation in phenotypic traits, but lacking such data for each population, we focus on correlations between phenotypic traits and geographical, environmental and genetic variables. Because genetic differentiation can emerge even with substantial gene flow (such as in a source-sink model or with strongly asymmetrical gene flow: Gaggiotti, 1996;Lowe & Allendorf, 2010), we next develop models varying the demographic processes that underlie population genetic patterns. We employ model selection approaches to compare these competing spatial models of gene flow for each grasshopper species. We then extend this analysis to investigate whether an ecological determinant of gene flow, phenological overlap among populations, is consistent with our observed demographic estimates. Phenological mismatch has been observed to directly reduce gene flow across elevation gradients, particularly in the downslope direction (Peterson, 1997). Limited temporal overlap among populations along an elevational gradient would reduce gene flow and increase genetic structure. Our previous research has shown that phenology in these grasshopper species is temperature dependent and occurs earlier in warmer years (Nufio & Buckley, 2019). Our hypothesis is that the short-winged species will exhibit genetic differentiation, limited phenological overlap, and gene flow, and have a strong association between genetic divergence and phenotypic differentiation along the elevational cline. Conversely, long-winged species might exhibit limited genetic differentiation, a source-sink pattern of gene flow, high phenological overlap, and have no association of genetic divergence with phenotypic differentiation.

| Study area, species, and phenotypic variation
We examine grasshopper populations spanning an elevation gradient of 1,500-3,500 m along the 40th N parallel in Boulder County, Colorado, USA (Table 1, Figure 2). The sites are all grassy meadows, with somewhat denser vegetation at the lower elevation sites. The forest matrix at the montane sites is less continuous and has been more impacted by forest fires than the subalpine site. The alpine site is an open ridgeline, but separated from the subalpine site by continuous forest vegetation (aspens and pines). All four of the species occur in montane foothill habitat up to subalpine habitat, while M. boulderensis (part of the Melanoplus dodgei species complex, (Otte, 2012), A. clavatus and occasionally M. sanguinipes occur in alpine habitats. Melanoplus boulderensis is a neoendemic of the Rocky Mountains, likely having evolved in response to glacial cycles over the last million years (Knowles & Otte, 2000), while the other three species are widespread in North America. The grasshoppers eat a variety of forbs and grasses with the exception of the grass specialist A. clavatus. The species overwinter in an egg diapause and are univoltine (annual) with the exception of A. clavatus, which takes 2 years to develop at higher elevations (Alexander & Hilliard, 1964).  while these species also exhibit higher thermal limits . All species increase feeding and digestion rates with temperature, but the temperature dependence of performance (hopping distance) varies among species . Metabolic rates increase with elevation for M. boulderensis and M. sanguinipes (in regressions accounting for mass and temperature), but do not vary with elevation for the other species (Table 1). However, laboratory-based rearing experiments have suggested that some physiological traits in M. sanguinipes vary as a result of phenotypic plasticity Buckley et al., 2015). Furthermore, phenological patterns are known to be related to developmental plasticity, as the early-season species M. boulderensis and C. pellucida slow their development rate (in terms of physiological time) at the high-elevation site in response to warming (Buckley et al., 2015).

| Genetic sampling, DNA extraction, and genotype-by-sequencing
Grasshoppers were collected from the field in June and August 2016, from three to five different locations, depending on species (Table 1). Twenty adults of each species were collected from each location, where possible, with approximately even numbers of males and females. Specimens were frozen then transferred to 100% ethanol for storage prior to DNA extraction. DNA was extracted from hind femur F I G U R E 1 Species vary in the strength and direction of clines for morphological, physiological, and reproductive phenotypes along the elevation gradient (top). Femur length is a proxy for body size and critical thermal minima and maxima (CT min and CT max , respectively) indicate thermal constraints on performance. We approximate reproductive potential as the number of functional ovarioles (egg-producing tubules). We additionally depict egg and clutch (eggs within a pod) mass. Species and populations (color) additionally vary in the temperature dependence of hopping performance (mean distance, bottom). The data suggest that species are (left to right) cool adapted, warm-adapted, and thermal generalists. We depict means and standard error for field-collected grasshoppers. Data are previously published Levy & Nufio, 2015) muscle tissue using a Qiagen DNeasy ® 96 Blood & Tissue Kit (Qiagen).
We followed the manufacturer's protocol with one modification: 4 μl  In order to infer missing genotype calls, the vcf file for each SNP dataset was processed in FaStphaSe v1.4 (Scheet & Stephens, 2006), which was developed to impute genotypes given the observed allele frequencies and linkage among SNPs in each subpopulation. The number of random starts of the expectation-maximization algorithm was set to 10, and the program was run without haplotype estimation. Missing genotypes were estimated for each sample site within each species, and those genotypes with uncertainty (<70% probability) were left as missing data. SNPs were then filtered by removing loci that had an overall minor allele frequency of <2.0%. After this step, no loci or individuals contained missing data.

| Population genetic analyses
For each SNP we calculated the following quantities: minor allele frequency (MAF), observed (H O ) and expected (H E ) heterozygosity, within-population gene diversity (H S ), and Wright's inbreeding coefficient (F IS ). Each statistic was averaged across loci within populations, and also calculated overall for each species. Additionally, for each species, we calculated overall gene diversity (multilocus heterozygosity; H T ) and D est (Jost, 2008). These measures are useful for assessing the overall genetic diversity found within and among species. To provide quantitative estimates of genetic differentiation among populations, we estimated the fixation index (F ST ) in a pairwise manner using Nei's D (Nei, 1978) and a linearized version based on Latter's equation (Latter, 1972). All summary statistics were calculated using the "basic.stats" function in the hierF-Stat package (Goudet, 2005;Goudet & Jombart, 2015)  As a first step in exploring population structure, we used a principal components analysis (PCA) in R, on centered and scaled SNP datasets for each species. Second, we examined genetic admixture using SnmF (Frichot, Mathieu, Trouillon, Bouchard, & François, 2014), implemented using the "snmf" function in LEA 1.6.0 package (Frichot & François, 2015) in R. For each species, we conducted 10 independent runs at K = 2 to K = N, where K is the number of ancestral populations and N is the total number of populations sampled, and combined the results from each K using CLUMPP (Jakobsson & Rosenberg, 2007), via the "clumppExport" function in pophelper v2.0.0 (Francis, 2017) in R. CLUMPP averages the individual admixture coefficients from each run and these results were visualized using pophelper. As an alternative clustering method, we used the program admixture 1.3.0 (Alexander, Novembre, & Lange, 2009), which implements a maximum likelihood approach to estimating ancestry. We used a fivefold cross-validation procedure to identify the value of K (the number of ancestral populations) for which the model has best predictive accuracy, testing K = 1 to K = N where N is the total number of populations sampled. We conducted 10 independent runs at each value of K and compared the average cross-validation error. Data conversions were done using custom R scripts and plink 1.9 (Chang et al., 2015).

| Modelling rates and direction of gene flow
We set out to test alternative models of gene flow among populations and estimate gene flow rates. We employed an Approximate Bayesian computation (ABC) approach (Beaumont, Zhang, & Balding, 2002) to fit our data to a set of four spatially explicit models for each species ( Figure  size drawn from a log-uniform prior distribution (200-10,000). For each model, we simulated two million datasets each with 500,000 independent sequences of 10 base pairs, and drew SNPs in equal number to those available in the empirical dataset for each species. We then used arlSumStat in arlequin v3.5 (Excoffier & Lischer, 2010) to calculate population genetic summary statistics for each simulated dataset, as well as the SNP datasets for each species.
These statistics included population-and species-level estimates of diversity (total heterozygosity and segregating substitutions, segregating substitutions per population) and differentiation (number of pairwise differences per populations and pairwise F ST ).
ABC analysis was completed using the abc package (Csilléry et al., 2012) in R. First, we used a leave-one-out cross-validation to determine whether the applied ABC approach should, in principle, be able to distinguish between the proposed models. We used a cross-validation sample of 100 and tolerance levels of .001, .005, and .01 (accepting the 2,000, 10,000, or 20,000 of closest estimates for each model, respectively), with posterior model probabilities estimated using both a neural network ) and a rejection method (Beaumont et al., 2002). Second, we replicates in FaStSimcoal2 and R.

| Phenological differences among sites
To provide an ecological measure of gene flow potential between populations, we measured the degree of phenological overlap.
Weekly survey data from 1959 to 1960 were assembled from field notebooks as part of the Gordon Alexander Project (ghopc limate. color ado.edu). Weekly resurveys were conducted 2006-2015 (Nufio et al., 2010) following the original protocol, consisting of 1.5 personhour of sweep netting (divided among 1-3 surveyors) and 0.75 person-hours of searching for adults and juveniles that may have been missed by sweep netting (Alexander & Hilliard, 1969). We assessed phenological overlap between sites for each species using Pianka's symmetric metric of niche overlap, which ranges from 0 to 1 (no to complete overlap) (Fleming & Partridge, 1984). The metric accounts for the proportion of individuals that occur during each sampling period at each site and sums across the survey season. Because sites were sometime sampled on different days within a week, we grouped the sampling period by week. We additionally quantified phenological overlap by calculating the day of year that the 15th and 85th percentile adult was surveyed for each site and species within each season.

| Diversity and population genetic structure
A total of 18,604 to 29,213 SNPs were identified in the four species (  (Table 3) and 84.2 ± 13.8%, assignment to a single ancestral population, respectively).

| Association between phenotypic clines and predictor variables
We examined the correlation of phenotypic differences among populations on the elevational gradient and a set of predictor variables (geographical distance, elevational distance, environmental distance, and genetic distance) using Mantel tests and assuming a linear relationship (Table 4). Notably, only short-winged species had significant phenotypic associations, but not with genetic distance as a predictor variable. Significant (α = .05) tests were found in M. boulderensis for both growth (geographic distance r = .656, p = .04; elevation distance r = .800, p = .04; environmental distance r = .481, p = .04) and nearly significant (α = .10) tests were found for reproductive traits (geographic distance r = .853, p = .08; elevation distance r = .867, p = .08; environmental distance r = .863, p = .08). Significant (α = .05) tests were found in A. clavatus for growth traits (geographic distance r = .642, p = .04; elevation dis tance r = .786, p = .04; environmental distance r = .482, p = .04), and growth traits were nearly significant in association with genetic distance (r = .889, p = .08).

| Modelling gene flow
We compared four competing models of spatial gene flow for each    Note: Models refer to-dn: unidirectional downslope gene flow; up: unidirectional upslope gene flow; ss: bidirectional stepping stone gene flow; am: equal rates of gene flow in an island model; sk: source-sink gene flow. For some species, not all models could be tested using the neural network method because infinite values were produced for lower likelihood models. Therefore, the top two or three models were compared in the neural network, and the model(s) with the highest posterior probability were used for downstream parameter estimation (indicated by *). Goodness-of-fit tests were conducted for the top models, and a significant test indicates a deviation of the model simulated summary statistics from the observed summary statistics.

| Phenology and phenological overlap
TA B L E 6 Median parameter estimates for effective size (N e ) and migration rate (m), with 95% confidence intervals (in parentheses) based on Approximate Bayesian computation using the neural net function

Camnula pellucida
Stepping stone Stepping stone Stepping stone

| Local adaptation versus phenotypic plasticity in life-history traits
We tested the hypothesis that dispersal ability would predict patterns of genomic and phenotypic differentiation in four ecologically divergent grasshopper species along an elevation transect in the Rocky Mountains. Relative to the long-winged species, we expected the short-winged species to exhibit stronger genetic differentiation, limited phenological overlap and gene flow, and a strong association between genetic divergence and phenotypic differentiation along the elevational cline. While our results for genomic divergence, gene flow, and phenology were supported, we found more equivocal evidence for associations between genetic and phenotypic variation that might be consistent with local adaptation. . We quantify phenological overlap using Pianka's index, which accounts for the abundance of individuals and ranges from 0 to 1 (no to complete overlap). Lines indicate the phenological overlap (y-axis value) for the populations (elevations) connected by the line in a single survey year. We additionally examine phenological overlap as the shift in the phenology of the 15th and 85th seasonal percentile individual as a function of elevation (right) adaptation, as the strength of selection need not be exceptionally strong to maintain adaptive alleles (Kawecki & Ebert, 2004).
However, we do not see a clear statistical association between genetic distance and phenotypic divergence. In the short-winged species, statistically well-supported differences in life-history and physiological phenotypes along this elevational gradient Levy & Nufio, 2015)   .
Other studies with grasshoppers have found strong genetic structure at a small spatial scale while phenotypic patterns are decoupled from such genetic variation (Ortego et al., 2015;Tinnert & Forsman, 2017). Future efforts to understand phenotypic patterns in our focal species will need to employ reciprocal transplant experiments to directly test predictions about local adaptation and plasticity and to determine if fitness trade-offs persist despite gene flow across these small spatial scales (Hereford, 2009). We currently lack relative fitness data to confirm that the phenotypic clines are adaptive. However, Levy and Nufio (2015) documented differential clines in both reproductive potential and realized reproduction for these four species.

| Patterns of genetic divergence and modelling gene flow along elevational gradients
Across environmental gradients, the amount and direction of gene flow between populations influence the evolutionary processes controlling phenotypic variation. Gene flow is known to facilitate adaptation by increasing the genetic diversity upon which selection can act, but it is more likely to reduce divergent selection on an ecological gradient as it allows an influx of maladaptive alleles to populations along a cline (Lenormand, 2002 The neural network approach is thought to perform better when the number of summary statistics is large and to provide more accurate results Csilléry, Blum, Gaggiotti, & François, 2010). However, the result of the modelling exercise suggests considerable uncertainty remains, particularly in parameter estimates. It is clear that gene flow among populations is present, although limited along the elevational transect, but the direction and rate of gene flow may not be accurate. We acknowledge that our models relied on simplistic assumptions (e.g., a single rate of gene flow across all populations) that likely impact these parameter estimates. However, there is inherent difficulty in estimating gene flow (Marko & Hart, 2011), especially using approaches reliant on summary statistics.
We, therefore, conducted additional simulations to estimate gene flow from the joint allele frequency spectrum of each species, but bootstrapping suggested parameter estimates remained imprecise.
Uncertainty in gene flow estimates may be a consequence of insufficient sampling in our study system, or simply a consequence of the limits of the above methods in inferring demographic parameters in some evolutionary scenarios (Cayuela et al., 2018).

| Ecological determinants of gene flow
The amount and direction of gene flow are determined not only by dispersal capacity and population size, but also the degree of nonrandom mating. Variation in, for example, activity and reproductive timing could limit gene flow between populations in different environments (Sexton, Hangartner, & Hoffmann, 2014).
Later phenology at higher elevations may contribute to the prevalence of upslope migration. Consistent will limited phenological overlap reducing gene flow, the alpine populations of A. clavatus and M. boulderensis were found to be particularly genetically isolated and also to exhibit phenotypic differentiation. In contrast, the steepness of the elevational cline in phenology varies among years for M. sanguinipes, which may allow gene flow to override genetic differentiation. The short-winged species, A. clavatus and M. boulderensis, are active early in the season and are also less dispersive among the grasshopper species in our study. Thus, low phenological overlap may act to reduce gene flow and contribute to genetic differentiation. Early seasonal species are often adapted for rapid development, which can lead them to capitalize on warming and exhibit the most pronounced phenological shifts (Buckley et al., 2015). Phenological timing may be further constrained by snowmelt timing at high-elevation sites, leading to low phenological overlap in heavy-snow years. Phenological differences are frequently invoked as an agent of reduced gene flow in plants (Franks & Weis, 2009), but the genetic implications of phenological differences have been less studied for insects (Mopper, 2005). In one example, upslope migration was found to result from butterflies emigrating away from senescing host plants at low elevation (Peterson, 1997). Other studies have found reduced genetic variation and emerging genetic structure in high-elevation insect populations (Jackson et al., 2018;Polato et al., 2017), although it is unclear if this pattern is due to phenological mismatch.
Gene flow patterns could also be driven by how wind direction facilitates dispersal. In particular, the genetic clustering of M. san guinipes populations from the elevational extremes suggests an important role of long distance, passive dispersal by the wind.
During summer in the mountains, the prevailing direction of airflow near the ground is upward during the day and downward at night. This air movement, due to radiation and the gradient of air density, flows counter to the prevailing winds in our study region (Alexander, 1964). Reduced grasshopper activity at night leads to most movement occurring predominately upslope during the day.
Recent survey records suggest the greatest incidence of transient individuals occurs during periods of high temperature and low wind speeds in the upslope direction (C. Nufio, unpublished data; during high winds, grasshoppers are less likely to take flight).
As such, more migrants of M. sanguinipes may reach the alpine from the lowest elevation site than the mid-elevation sites due to higher population densities at low elevation. For M. sanguinipes, surveys found only adults (no juveniles) at upper elevation sites, suggesting a lack of reproduction or limited recruitment in these populations. This is reinforced by observed declines in reproductive potential (proportion of ovarioles that become functional) with elevation for M. sanguinipes (Levy & Nufio, 2015) suggesting that low elevation migrants are poorly suited to the alpine.
Populations composed of transient individuals may therefore represent a sink population, rather than a stable local population (Gaggiotti, 1996).
Climate change could quickly alter genetic structure in these species. Dispersal may increase as warm days increase the upward flow of air near the ground. Population differentiation, as observed in the alpine populations of A. clavatus and M. boulderen sis, could subsequently erode and upslope flow of warm-adapted alleles could facilitate adaptive responses to climate change.
However, climate change could also decrease phenological overlap among populations and exacerbate population differentiation.
High-elevation populations might then contract and experience greater fragmentation (Rubidge et al., 2012). An example of environmental conditions influencing phenological overlap is provided by plant responses to drought: drought alters phenology more for wet-adapted plants than dry-adapted plants, increasing phenological overlap (Franks & Weis, 2008). Because grasshopper populations differ in the temperature dependence of their developmental rates (Buckley et al., 2015), and previous documentation of phenological advancements show an association with climate warming (Nufio & Buckley, 2019), we predict that phenological timing will move earlier in the year and will be exacerbated by changes to snowmelt and plant phenology patterns along the gradient. Interspecific phenological mismatch within sites has been extensively studied as a driver of altered species interactions (Visser & Both, 2005), but intraspecific phenological mismatch may be equally important to climate change outcomes as it reduces gene flow among sites.

ACK N OWLED G M ENTS
We would like to thank Jack Cook for assistance with DNA extrac-

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
All authors assisted in paper writing. L.B.B., R.A.S., and S.D.S. conceived the study and L.B.B., C.R.N, and R.A.S. collected specimens.
R.A.S. and S.D.S. conducted the genetic data collection and analyses, C.R.N. conducted the phenological surveys.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw Illumina reads were deposited in the National Center for Biotechnology Information Short Read Archive (accession nos.