Phenotypic and genetic divergence in a cold‐adapted grasshopper may lead to lineage‐specific responses to rapid climate change

Species responses to global warming will depend on intraspecific diversity, yet studies of factors governing biogeographic patterns of variability are scarce. Here, we investigate the evolutionary processes underlying genetic and phenotypic diversity in the flightless and cold‐adapted grasshopper Sigaus piliferus, and project its suitable space in time.


| INTRODUC TI ON
Understanding the forces shaping intraspecific variability is a keystone in evolutionary biogeography and conservation, providing important insights on the factors determining biotic responses to environmental changes (Seaborn et al., 2021).Spatial patterns of intraspecific variation reflect the influence of past environmental fluctuations leading to divergent evolutionary trajectories (Gillespie & Roderick, 2014;Grant & Grant, 2017;Hewitt, 2004).Quaternary climates, and in particular late glacial stages, were major drivers of biogeographic reorganisation around the world (Davis & Shaw, 2001;Hewitt, 2004;Stewart et al., 2010), and it is apparent that the effects of Pleistocene climate change on biota covaried with latitude, topography, and intrinsic species traits such as dispersal behaviour and ecological preferences (Byrne, 2008;Davis & Shaw, 2001;Đurović et al., 2021;Hewitt, 2004;Stewart et al., 2010).Polar ice cap extension and retreat in the Northern Hemisphere is inferred to be a key driver of population range shifts in many species, whereas Southern Hemisphere landscapes were little effected in this way (Byrne, 2008;Carmelet-Rescan et al., 2021;Meza-Joya et al., 2023;Trewick et al., 2000).Conditions during glacial phases appear to have promoted range expansion of cold-adapted biota, especially in regions less affected by polar ice sheet such as mountains in the tropics (Hewitt, 2004), the Japanese Alps (Ikeda & Setoguchi, 2007), and Southern Alps in New Zealand (Trewick et al., 2000).These range shifts shaped distinct genetic and phenotypic variants upon which deterministic (e.g.selection) and stochastic forces (e.g.genetic drift) acted to mould modern patterns of intraspecific diversity (Gillespie & Roderick, 2014;Grant & Grant, 2017;Hewitt, 2004).
The New Zealand Southern Alps (Kā Tiritiri o te Moana) dominate the South Island landscape, climate and to a large extent biogeography, but North Island, New Zealand (Te Ika-a-Māui, Aotearoa) has a geologically younger and more volcanic landscape.Here the study of biogeographic partitioning has to consider the extensive geophysical remodelling that was coincident with Plio-Pleistocene climate shifts (Figure 1).Most of southern North Island emerged as land since the late Pliocene (3 Mya) and late Pleistocene land bridging allowed exchange of terrestrial biota between the main islands (Trewick & Bland, 2012).Rapid uplift of the southern North Island axial ranges during the Pleistocene (~1 Mya) generated heterogeneous environments (Trewick & Bland, 2012) that appear to have stimulated lineage diversification (Heenan & McGlone, 2013;Marshall et al., 2012;Prebble et al., 2021).More recent sea straits (~450 kyr) may have isolated terrestrial populations but gene flow was possible during lowered sea levels of the later glaciations (Trewick et al., 2017;Trewick & Bland, 2012).Repeated and lengthy Pleistocene glacial phases (Rother et al., 2014), drove vegetation turnover that may have repeatedly reshaped animal biogeography (e.g.Trewick et al., 2000).In particular a shift between forest during warmer stages and alpine-like scrub grassland during cooler, drier phases (McGlone, 1985).Volcanism in the centre of the island caused local extinction (Trewick et al., 2011), and genetic data indicate intraspecific secondary contacts in many plant and animal species (Shepherd et al., 2022), and instances of hybridisation have been detected between cytogenetically divergent populations of the tree wētā Hemideina thoracica (Morgan-Richards et al., 2000).
Here we examine drivers of intraspecific diversity in Sigaus piliferus Hutton (Orthoptera: Acrididae), the sole North Island representative of an endemic New Zealand radiation of flightless, cold-adapted grasshoppers (Koot et al., 2022).According to fossil calibrated phylogenetic analysis this species diverged from a common ancestor with the larger high-alpine Sigaus villosus Salmon ~14 Mya (17-11 Mya) during the Middle Miocene (Koot et al., 2020), so its current patchy and relatively low latitude distribution could reflect responses to late Neogene geoclimatic conditions.Like its South Island counterparts, S. piliferus is common in alpine areas, but unlike them, it also naturally occurs in sub-alpine habitats, where it seems to be rare and subject to local population extinction over the last 80 years (Trewick & Morris, 2008).As expected of a small, flightless endemic insect, this species displays apparent phenotypic (Bigelow, 1967) and genotypic (Trewick & Morris, 2008)  Integration of genetic, phenotypic and climatic data provides a robust framework to investigate possible drivers of local adaptation and to identify populations that require conservation resources to avoid loss of evolutionary potential (Dowle et al., 2015;Merilä & Hendry, 2014;Mullen et al., 2009).Genetic analyses using mitochondrial sequences provide information about the influence of palaeo-geophysical changes on population connectivity and demography, and correlation with morphological variation.
Ecological niche models allow us to infer persistence and vulnerability of suitable environmental space through time based on climate proxies.Observations from other cold-adapted New Zealand grasshoppers (e.g.Carmelet-Rescan et al., 2021;Meza-Joya et al., 2023) suggest glacial climates could have allowed S. piliferus range to expand compared to interglacials.By exploring correlations between phenotypic variation and environmental and genetic (mtDNA) divergence we provide insights into the processes underlying population structure and differentiation across this species range.The active geophysical history of North Island must have been influential for the genetic (Trewick & Morris, 2008) and morphological (Bigelow, 1967) adaptation, conservation, intraspecific variation, Quaternary climates and volcanism, range shifts populations of this species, so we would expect both deterministic and stochastic evolutionary processes to explain current patterns of intraspecific divergence.If as expected, past climates influenced the spatial distribution of intraspecific diversity in this grasshopper, then anthropogenic warming during this century will also influence surviving diversity.

F I G U R E 1
North Island, New Zealand palaeogeography against the modern coastline, indicating major geophysical features since late Pliocene time.Tectonic uplift (Trewick & Bland, 2012), lowered sea level (~28 kyr; Rother et al., 2014) and continuous forests contraction at the last glacial phase (Newnham et al., 2013), with area affected by widespread pyroclastic flow from Oruanui eruption in the Taupō Volcanic Zone (TVZ) at ~25.5 kyr (Wilson, 2001).Sampling of Sigaus piliferus location names are available in Table S1.1 in Appendix S1.Map projection: NZGD2000.

| Sample collection and DNA extraction
We collected Sigaus piliferus from 33 locations spanning its documented range (Figure 1; Table S1.1 in Appendix S1).Most specimens came from elevations above 1000 m, but lower elevation populations were also identified and sampled.Grasshoppers were collected during summer when they were active (December to March 2004-2022) and frozen before being preserved in 95% ethanol.Taxonomic identification, sex, and maturity were determined based on morphological traits (Bigelow, 1967).Whole genomic DNA was extracted from leg muscle of 186 specimens (Table S1.1 in Appendix S1) using a solvent-free Proteinase K and salting-out method (Trewick & Morgan-Richards, 2005).Extracted DNA was amplified by PCR for the mitochondrial protein-coding NADH-dehydrogenase 2 (ND2) gene under standard conditions, using the primers HopND2_147F and HopND2_1286R (Carmelet-Rescan et al., 2021).This marker has a higher capacity to accumulate haplotype diversity in these grasshoppers than the commonly analysed COI locus (Carmelet-Rescan et al., 2021).Sequencing reactions used BigDye Terminator v.3.1 (Life Technologies) with signal capture on an ABI-3730xl System (ThermoFisher).DNA sequences were edited and aligned in Geneious vR10 (Kearse et al., 2012).

| Population genetic structure and demography
Genealogical relationships among ND2 haplotypes were inferred using a statistical parsimony network algorithm implemented in PoPArt 1.7 (Leigh & Bryant, 2015).Statistical support for distinct mitochondrial lineages was evaluated with a Bayesian phylogenetic analysis in MrBAyes 3.2.6 (Ronquist et al., 2012) with Sigaus villosus as outgroup taxon (Koot et al., 2020).The analysis used four chains on two runs for 20 6 generations, with sampling frequency of 20 3 generations and a burn-in of 0.10.Pairwise genetic distances between population samples were estimated using MeGA X (Kumar et al., 2018) based on the Kimura 2-parameter (K2P; Kimura, 1980) and Tamura-Nei (TN; Tamura & Nei, 1993) models with 1000 bootstrap replicates.We used spatial Principal Component Analysis (sPCA) in the R (R Core Team, 2022) package 'adegenet' 2.1.3(Jombart, 2008) to examine non-random spatial patterns of genetic variability.Spatial distances were based on the K nearest neighbours method (K = 25), which allowed possible connections at regional scale, and statistical significance was assessed using an Eigenvalue test (Montano & Jombart, 2017) with 10,000 permutations.We tested for genetic differentiation among the identified sPCA genetic clusters using an AMOVA with the R package 'poppr' 2.9.3 (Kamvar et al., 2014) with 10,000 permutations.
Matrilineal DNA variability in population samples (n > 5) was estimated using DnAsP 5.10 (Librado & Rozas, 2009) to compute haplotype (h) and nucleotide diversity (π).Pairwise Φ ST values (Excoffier et al., 1992) were computed to test genetic differentiation among populations with the R package 'haplotypes' 1.1.2(Aktas, 2020) with 10,000 replicates.Historical demography was inferred using mismatch distributions under the constant population size model with DnAsP.The fit of the observed distribution to the null model was evaluated using the Harpending's r raggedness index (Harpending, 1994).We also calculated three neutrality statistics to detect departures from the mutation-drift equilibrium indicative of population size change under the assumption of locus neutrality: Tajima's D (Tajima, 1989), Fu's F S (Fu, 1997) and Ramos-Onsins' R 2 (Ramos-Onsins & Rozas, 2002).Statistical significance (p < .05except for F S where p < .02)was evaluated with 10,000 simulations.

| Morphological variation
We applied metric and geometric morphometrics to identify phenotypic variation and test for correlation with other factors (Figure S2.1 in Appendix S2).Four non-independent body dimensions (hind femur length, FL; hind femur width, FW; pronotum length, PL; pronotum width, PW) and two ratios (FL/FW and PW/PL) were recorded from 156 adult grasshoppers (n = 100♀, 56♂; Table S2.1 in Appendix S2) using an Olympus SZX7 stereomicroscope with Olympus SC100 image capture and Olympus cellSens Dimension v1.6 software (Olympus Corporation, Japan).These metrics were chosen as they are not arbitrarily affected by preservation (Bigelow, 1967), correlate with body mass (Meza-Joya et al., 2022), and are informative of shape and size variation in S. piliferus (Bigelow, 1967).Geometric analysis of pronotum shape used 14 digitalised landmarks on the dorsal perimeter of the pronotum in images from 161 adult grasshoppers (n = 101♀, 60♂; Table S2.2 in Appendix S2) using tPsDiG2 2.29 (Rohlf, 2015).This structure is informative for analysing size-independent shape variation in New Zealand grasshoppers (Carmelet-Rescan et al., 2021;Morgan-Richards et al., 2022).Images were taken using a Canon EOS 600D with EF100 mm f2.8 USM macro lens (Canon Inc., Japan) mounted on a vertical stand (Kaiser Fototechnik, Germany).A generalised Procrustes analysis (GPA) was applied to remove non-shape variation from the dataset using the R package 'geomorph' 4.0.0 (Adams et al., 2020).Repeatability analyses indicated any technical measurement error was biologically negligible (Figure S2.2 in Appendix S2).
Preliminary analyses showed the studied traits are sexually dimorphic in this grasshopper (Figure S2.3 in Appendix S2) and thus, only adult females were analysed hereafter due to better sample size and more homogeneous geographic sampling compared to males.
We modelled the most likely number of morphological clusters within our datasets (size and shape) without prior information on specimen origin using Gaussian mixture modelling with the R package ' Mclust' 5.0.2 (Scrucca et al., 2016).For the traditional approach, analyses considered the total dataset (i.e.four metrics and two ratios), while geometric data clustering was computed using meaningful components from a Principal components analyses (PCA) on size-uncorrected Procrustes coordinates (i.e.PC1).Exploratory multivariate models with PERMANOVA design with 10,000 permutations (shape ∼ log(centroid size)) indicated size influence on shape was minimal (R 2 = .041,p = .0004)and classification performance increases when using shape and size at the same time rather than only shape (Mitteroecker & Schaefer, 2022).Canonical variate analyses (CVA) with 10,000 permutations were conducted in the R package 'Morpho' 2.9 (Schlager, 2016) to investigate the degree of dissimilarity between model-based groups using the Mahalanobis distance (Md) with Bonferroni adjusted p value (Klingenberg, 2013).

| Ecological niche modelling
To estimate and project the potential niche of S. piliferus we compiled a presence-absence database (geodetic datum WGS84) including data for 14 related, endemic New Zealand grasshopper species, spanning 1967-2016 (Koot et al., 2022), plus recent sample collections and verified, high resolution observations (uncertainty <1 km 2 ) from iNaturalist NZ (https:// inatu ralist.nz).We also included target-group absences (Mateo et al., 2010) (1960-1990), and future (2061-2080).We projected two future scenarios with mean temperature increases of 1.4°C (RCP4.5)and 3.0°C (RCP8.5) in New Zealand by 2090 (Ministry for the Environment, 2018).Additional variables capturing geodiversity and topographic traits were discarded early in the variable selection process (Table S3.1 in Appendix S3).Layer files were at 30 arc-seconds resolution (~1 km 2 ), except LGM which was at 2.5 arc minutes (~5 km 2 ).All layers were cropped to the extent of New Zealand, which represents the potential accessible area for the species during the relevant time period given land connections between major islands of the archipelago during Pleistocene glacial phases (Trewick & Bland, 2012).Layer processing used QGIS 3.16.1 (QGIS Development Team, 2020).To reduce multicollinearity, we selected a subset of five independent and ecologically informative variables (Table S3.1 in Appendix S3) using the R package 'fuzzySim' 3.0 (Barbosa, 2015).
Three modelling methods were used (Generalised Boosting, Random Forest and Maximum Entropy) with settings optimised by testing alternative hyperparameters (Tables S3.2 and S3.3 in Appendix S3) with the R package 'SDMtune' 1.1.4(Vignali et al., 2020).Final models were refitted using the most relevant environmental predictors and the best performing hyperparameters in the R package 'biomod2' 3.4.6(Thuiller et al., 2020), through a spatialblocking cross-validation approach (Figure S3.2 in Appendix S3) with the R package 'blockCV' 2.1.1 (Valavi et al., 2019).Five runs were completed for each modelling method, resulting in 15 individual models.Model performance was evaluated using the area under the receiver operator characteristic curve (AUC) and true skill statistic (TSS).Models with high predictive power (TSS and AUC > 0.9; Figure S3.3 in Appendix S3) were ensembled using a weighted mean approach (EMmw), and variable importance for the EMmw (Figure S3.4 in Appendix S3) estimated (Fletcher et al., 2016).Time projections were generated using the EMmw and binary (presence/ absence) maps were built using the cut-off value that maximised the rate of presences and absences correctly predicted.Extrapolation risk was assessed through a multivariate environmental similarity surface (MESS) analysis (Elith et al., 2010) with the R package 'dismo' 1.3-3 (Hijmans et al., 2020).

| Climate niche analyses
We assessed niche differentiation within S. piliferus using a principal component analysis (PCA) approach (Broennimann et al., 2012), as implemented in the R package 'ecospat' 3.5.1 (Di Cola et al., 2017) with 5000 random points across North Island to summarise the climate space available for each mtDNA lineage.Niche overlap between lineages was calculated using a corrected version of the Schoener's D index (Di Cola et al., 2017), which ranges from 0 (no overlap) to 1 (total overlap) of climate space.We used the niche similarity test (Warren et al., 2008) to assess whether the observed D values were more (i.e.niche conservatism) or less similar (i.e.niche divergence) than expected by chance using 1000 permutations.We quantified niche marginality (M) and specialisation (S) for each lineage with the R package 'CENFA' 1.1.0(Rinnan, 2018).The higher the marginality, the more the niche deviates from the mean climate conditions and the higher the specialisation, the narrower the niche (Hirzel et al., 2002).This analysis used the actual presence data and the climatic background was represented by the subset of climate variables used for niche modelling (Table S3.1 in Appendix S3), as they are not strongly correlated and are important for predicting this grasshopper climate envelope.

| Processes underlying population divergence
We used multiple matrix regression with randomization (MMRR) to explore the potential processes underlying differentiation among populations of S. piliferus.For this, we built pairwise Euclidean distances for body phenotype (FL, FW, PL, and PW) with the R package 'vegan' 2.5-7 (Oksanen et al., 2020), and pronotum shape with the R package 'Geomorph' using a PERMANOVA design with 10,000 iterations (shape ∼ population).Euclidian distances for the subset of variables used for niche modelling (Table S3.1 in Appendix S3) plus mean annual temperature (hereafter temperature) were estimated using the R package 'Stats' 4.0.3(R Core Team, 2022), as they are expected to capture climate-related factors influential for phenotypic variation in orthopterans (e.g.Whitman, 2008).Elevation was excluded as it was highly correlated with temperature (r = .808,p = .001).Pairwise ND2 Φ ST and topographically-corrected geographic distances (km) were used as proxies of genetic and geographic variation (see above).We performed MMRRs with 10,000 permutations in the R package 'PopGenReport' 3.0.7 (Adamack & Gruber, 2014).We fitted separate models to assess the effects of geography, climate, and genetic divergence on body phenotype and then in pronotum shape differentiation.We used a backward elimination procedure to retain the most informative variables to be included in the final models by fitting a starting model with all explanatory variables, and then removing iteratively those with the lowest contribution to the model (i.e.highest p-value) until only significant variables remained (p < .05).Analyses were conducted considering 27 population samples for which we had genetic and phenotypic data.

| Population genetic structure and demography
We obtained a 747 bp alignment of mitochondrial ND2 sequence for 186 individual Sigaus piliferus (Table S1.1 in Appendix S1).The analysed fragment exhibited an A + T bias (A + T content = 0.752, G + C content = 0.248) as is typical for insects (Trewick & Morgan-Richards, 2005).The aligned ND2 sequences rendered an unambiguous region with 47 amino acid substitutions and 155 polymorphic nucleotide sites, comprising 91 haplotypes.
Genealogical relationships based on a parsimony haplotypic network revealed a deep north-south phylogeographic split either side of the Manawatū Basin (Figure 2).The southern lineage included 32 samples from five locations in the Tararua Range, while the northern lineage comprised 154 samples from 25 locations (Table S1.1 in Appendix S1).These deepest lineages were separated by 42 mutational steps and average genetic distances between 9.6% (TN) and 9.7% (K2P).Mean within-lineage distance were 2.0% (K2P and TN) in the north, and 1.1% (K2P and TN) in the south (Figure S1.1 in Appendix S1).The more widespread northern lineage had some geographic structure, the most salient separating most samples from Kaimai Range and other populations by 36 mutational steps (6.8% TN; 6.9% K2P).This general structure was retrieved by our phylogeny (posterior probabilities >.95), yet shallow divisions were not resolved (Figure 2 The examined mtDNA ND2 fragment showed high haplotype diversity (h = 0.965 ± 0.007) and moderate nucleotide diversity (π = 0.038 ± 0.003) within S. piliferus.Haplotype diversity was marginally higher in the southern lineage (0.984 ± 0.012 vs. 0.950 ± 0.009) but nucleotide diversity was almost two times more in the northern lineage (0.020 ± 0.002 vs. 0.011 ± 0.002), likely reflecting the influence of unequal sample sizes (northern = 154 vs. southern = 32 samples).When examining population samples (n > 5), haplotype diversity ranged from 0.011 to 0.861 and higher values were generally found in those from central North Island.
In turn, nucleotide diversity ranged from 0.0003 to 0.0150 and values were generally higher at northernmost (Kaimai Range) and southernmost samples (Tararua Range), with differences reaching two orders of magnitude regardless of sample size (Table S1.2 in Appendix S1).Genetic differentiation among population samples was found to be moderate to strong and significant in most cases Tauhara, Mangaio, Kaweka, Kuripapango and Whanahuia (Figure 2; Table S1.3 in Appendix S1).The star shape of the haplotypic network including these populations is indicative of recent population expansion (Figure 2).Statistical support for demographic expansion was also found for one south population sample (Herepai), others fitted the model of constant population size (Figure 2; Table S1.3 in Appendix S1).

| Morphological variation
Phenotypic clustering (size and shape) comprising two overlapping morphological clusters was broadly consistent with the major genetic partitions (Figure 3).The best-fitted Gaussian model for linear size (ellipsoidal, equal shape; BIC = −879.84)mostly overlay the two major lineages of the molecular analyses (22 misassignments, 88% accuracy; Figure 3a).In general, females in lower latitudes/elevations (north) have larger femora and pronota than those in higher latitudes/elevations (south), although misclassifications were evident in the central area.Shape variation was only loosely associated with major genetic structure (40 misassignments, 60% accuracy;  Considering extrapolation risk (Figure S3.6 in Appendix S3), hindcasts indicate suitable space was more extensive and connected during interglacial (LIG = 21,448 km 2 and MH = 13,776 km 2 ) than glacial periods (LGM = 6868 km 2 ).North Island habitat space at the last interglacial was predicted to be similar to current predictions, whereas much less habitat was predicted at the LGM and this restricted to north (Coromandel-Kaimai ranges), northeast (Raukumara Range) and south (west Tararua Range).Climate warming during the mid-Holocene pushed suitable conditions upslope, resulting in the loss of 36% of suitable space.Similarly, the suitable habitat envelope is expected to move poleward and upward under future warming scenarios (Figure 5; Figure S3.5 in Appendix S3).After accounting for extrapolation risk (Figure S3.6 in Appendix S3), habitat loss was predicted to be more severe under the most pessimistic climatic scenario (RCP 8.5 = 7930 km 2 vs. RCP 4.5 = 10,139 km 2 ), but in all cases involves loss of at least half of the current suitable space (RCP 4.5 = 49% and RCP 8.5 = 60%) and extensive fragmentation.

| Climate niche analyses
The first two axes of the PCA based on five climate variables (Table S3.1 in Appendix S3) across the whole study area, explained 48.93% and 21.82% of the total environmental variance, respectively (Figure 5a).Niche segregation between mtDNA lineages was closely related to the first PCA axis, which had a strong negative relation with mean temperature of wettest and driest quarter (−0.95 and −0.88, respectively).Niche overlap was minimal (D = 0.107), indicating these lineages have non-equivalent niches.Despite this, we found evidence of niche conservatism as the similarity test indicated that niches were significantly more similar than expected by chance in both recipro-

| Processes underlying population divergence
In the backward elimination procedure, mtDNA distance was removed early as having no association with grasshopper phenotype.
The environmental factors considered explain a moderate amount of phenotypic variation (capture by linear measurements of hind legs and pronotum) in S. piliferus, suggesting adaptive divergence but pronotum shape was not informative (Table 1).The best-fitted model for body phenotype (R 2 = .424,p < .0001)indicated divergence was significantly positively affected by temperature (β = .377,p < .001),diurnal range (β = .375,p = .003),and to a minor degree geographic distance (β = .202,p = .034),and negligible yet negative effects were found for precipitation seasonality and temperature of wettest quarter (Table 1; Figure S4.1 in Appendix S4).In contrast, the F I G U R E 4 Predicted suitable niche space for Sigaus piliferus (cut-off = 488.5) in North Island, New Zealand during three past epochs (last interglacial, last glacial maximum, and mid-Holocene), current time, and two future scenarios (RCP 4.5 and RCP 8.5).Map projection: NZGD2000.
best-fitted model for pronotum shape (R 2 = .070,p = .0024)showed only geographic distance had a weak significantly positive effect on shape differentiation (β = .004,p = .002),but the low explanatory power of this model indicates biological irrelevance.

| DISCUSS ION
As in other mid to high-latitude regions of the world (e.g.Đurović et al., 2021;Endo et al., 2015;Ikeda & Setoguchi, 2007;Stewart et al., 2010), Pleistocene glacial phases allowed alpine conditions to flourish in New Zealand (Heenan & McGlone, 2013;McGlone, 1985), and cold-adapted organisms are inferred to have attained larger ranges and population sizes in response (Carmelet-Rescan et al., 2021;Meza-Joya et al., 2023;Trewick et al., 2000Trewick et al., , 2011)).We tested this expectation for the flightless, cold-adapted New Zealand grasshopper Sigaus piliferus but unexpectedly found a F I G U R E 5 Distinct niche and habitat features for grasshopper lineages within Sigaus piliferus.Niche overlap along the first two axes of a principal component analysis on five climatic variables (see Table S3.Newnham et al., 2013).Our models suggest suitable habitat for S. piliferus during the LGM was restricted to the latitudinal limits of their current range (Figure 4), where distinct climate regimes prevailed (Drost et al., 2007).We suggest this hindcast reflects intraspecific niche variation that is consistent with our field observations of northern and southern populations (Figure 5).Climate based predictions for northern populations are more widely represented across the landscape, indicating the lineage expanded its thermal regimen towards increased warm tolerance compared to southern populations of S. piliferus and sibling species in South Island.In contrast, the southern populations show niche specialisation to a rarer habitat that is expected to be more vulnerable to future warming.
These distinct niche features suggest variability in population-level response to climate change (e.g.Maguire et al., 2018) and could reinforce geographic isolation with future climate warming and promote ecological speciation (see Wiens & Graham, 2005).
Ecological responses to climate and volcanic disturbance at These contrasting histories suggests populations at the range edges are critical for this species adaptive capacity.If these patterns also hold for future climate scenarios, the narrowly distributed southern lineage may be especially threatened.
Volcanism has been an important influence on North Island biogeography during the last two million years (Trewick & Bland, 2012).
Stratigraphic features (Manville et al., 2009;Segschneider et al., 2002) and palynology (Wilmshurst & McGlone, 1996) indicate pyroclastic flow from the most recent Taupō super-eruption (~1.8 kyr, Wilson & Walker, 1985) mantled the landscape, destroyed vegetation and disrupted populations across the centre of the island (Morgan-Richards et al., 2000;Shepherd et al., 2022;Trewick et al., 2011).Together with continuous activity on the central volcanic plateau (Ngāuruhoe, Tongariro and Ruapehu volcanos) since ~0.28 Ma (Cronin et al., 1996;Trewick & Bland, 2012), vegetation loss and recovery (Wilmshurst & McGlone, 1996) would have been critical for herbivorous insect population dynamics.Volcanic destruction of forest in the central volcanic zone led to replacement by pioneer ferns and grasses (Wilmshurst & McGlone, 1996), and the resulting mosaic of open habitats could provide ecological opportunities for grasshopper populations to flourish.
Population recovery after repeated volcanism is typically characterised by recolonisation from diverse source populations outside the volcanic zone (Baranzelli et al., 2020;Muir et al., 2000;Osuna et al., 2020) Joya et al., 2023;Trewick, 2008;Trewick et al., 2000) and coldadapted insects elsewhere (e.g.Borer et al., 2010;Endo et al., 2015;Martinez-Sañudo et al., 2022).However, the elevated topography in North Island is mostly younger (<1 Mya) and this suggests that even divergent S. piliferus lineages (e.g.currently associated with the Kaimai range ~6.9%) were formerly associated with habitat that differs from most South Island grasshoppers.Beyond the central volcanic zone haplotype sharing was limited to two pairs of adjacent populations (Pinnacles-Te Aroha and Pirongia-Karioi).
In keeping with mtDNA signal for two main lineages (northern and southern), we found phenotypic variation roughly comprised two main types ('north' and 'south'), but mtDNA and phenotypic differentiation were not precisely correlated.We found most morphometric variance among grasshoppers in the central part of the species range (Figure 3) consistent with high variation noted by Bigelow (1967).The co-occurrence of disparate ('northern' and 'southern') phenotypes in the central region most likely reflects mixing of populations with distinct nuclear genetics not detected by the single locus mtDNA marker.We suggest that mixed phenotype indicate a non-equilibrium population driven by volcanic remodelling of the habitat that led to extirpation and reinvasion from multiple sources.Future research using multilocus nuclear markers would illuminate the evolutionary history and future trajectories in these populations, and this may be especially useful for conservation planning in these lineages.
The latitudinal gradient is expected to result in warmer conditions to the north compared to the south (further from the equator), and temperature is widely recognised as influencing body size differentiation in Orthoptera (Whitman, 2008).A tendency towards larger size at warmer, lower latitudes and elevations, is the most frequently reported pattern in grasshoppers and other orthopterans (García-Navas et al., 2017a;Mousseau, 1997;Shelomi, 2012), a presumable adaptive response (Walters & Hassall, 2006).We found body size in S. piliferus, inferred from leg dimensions, was associated with temperature gradients, yet body shape, including relative leg dimensions, might also represent adaptation to other environmental attributes such as vegetation type that is itself linked to temperature and seasonality (Wardle, 1985).We found that S. piliferus living in northern sites with montane woody vegetation (e.g.Coromandel- Kaimai ranges) had more elongated bodies and longer hind legs than those in southern open alpine grasslands (e.g.Tararua Range) that had stout bodies and robust femora (Figure 5).This same association of  et al., 2017b;Laiolo & Obeso, 2015;Santos et al., 1994), some of the observed variation can reflect plastic effects, likely in concert with genotypic variation (Laiolo & Obeso, 2015).We propose two S. piliferus ecotypes be recognised ('shrubland' and 'alpine',Figure S4.2 in Appendix S4) on the basis of readily identifiable morphological differences and habitat (Gregor, 1944;Nosil, 2007), as these may have distinct evolutionary trajectories and respond differently to rapid environmental change.
Overall, if niche stability and saturation operate, suitable space for S. piliferus is predicted to diminish during the intensified 'anthropogenic' interglacial now developing around the globe.In keeping with this, loss of suitable space was also predicted at mid-Holocene, which was warmer than today.The extent of habitat loss and fragmentation predicted in North Island (49%-60%) will likely lead to increased habitat patchiness and local extirpation (see Meza-Joya et al., 2023) with different effects on the ecotypes.Climate change is expected to reduce alpine vegetation on mountains and thus grasshoppers living there, but could increase opportunity for the northern shrubland ecotype.While the climatic model strongly predicts the realised niche of S. piliferus, variables used for modelling are proxies for drivers of this species range, and neglect the ecological factors setting range limits and responses to environmental shifts (Wiens, 2011).The existence of relict northern populations on the margins of exotic pine plantations suggest survival under alternative conditions is possible, but such silvicultural systems are in ecological disequilibrium, subject to invasive exotic pests (Crous et al., 2017), and may create population sinks or ecological traps for native biota (Pawson et al., 2010), threatening persistence on these grasshoppers.
Predicting species' fates under anthropogenic climate change remains a challenge due to uncertainties on how to characterise the adaptive capacity of biotic units across the population-species continuum and determine robust metrics for prioritising populations for conservation (Coates et al., 2018;Seaborn et al., 2021).The predicted upward shift of suitable space for southern populations (Tararua Ranges) during climate warming will further reinforce geographic isolation, and this might be an instance where anthropogenic climate change drives evolutionary divergence (see Wiens & Graham, 2005).Contrary to expectations of warming driving selection for smaller body size (Gardner et al., 2011), warming conditions seems to select for larger sizes in parts of this grasshopper system.This highlights the need to incorporate phenotypic data and functional genetics in studies of species response to climate change to gauge the rate and likely outcomes.Despite current levels of intraspecific diversity in S. piliferus indicating adaptive capacity, the rate of habitat loss predicted over the next 50 years of warming will lead to small and fragmented populations (Koot et al., 2022) with limited capability to adapt via contemporary evolution to the novel selective forces arising from rapid environmental pressures (Kardos et al., 2021;Radchuk et al., 2019;Willi et al., 2022).

ACK N OWLED G EM ENTS
We are grateful to the New Zealand Department of Conservation

CO N FLI C T O F I NTE R E S T S TATE M E NT
We declare we have no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
The mtDNA sequences supporting the findings of this study are spatial structure, that suggest a combination of restricted gene flow and local selective pressures.The absence of near relatives of S. piliferus in North Island and ability of this species to inhabit alpine and nonalpine ecosystems is suggestive of ecological flexibility and/or local adaption.
consisting of sites where sampling effort has produced other Orthoptera (Gryllidae and Tettigoniidae), but not S. piliferus.Duplicate records were removed ensuring one datum per pixel (i.e.~1 km 2 habitat area).The resulting database consisted of 1188 site records, including 1053 absences and 135 presences, representing the realised niche of S. piliferus (Figure S3.1 in Appendix S3).As potential niche proxies we used 19 bioclimatic layers from WorldClim v1.4 for five periods (MIROC-ESM and NCAR-CCSM climate models): last interglacial at ∼120-140 kyr (LIG), last glacial maximum at ~28 kyr (LGM), mid-Holocene at ∼6 kyr (MH), current , Figure S1.2 in Appendix S1).Most haplotypes were private (95.5%), and the most common haplotypes (30 and 31), separated by two mutational steps, were found in population samples from central North Island (maximum geographic distance ~135 km).The sPCA Eigenvalue test supported global (p = .0001)but no local struc-ture (p = .9999),consistent with a weak signature of isolation by distance over the entire species range (r = .385,p = .0001)and no regional patterns: north (r = −.006,p = .5426)and south (r = −.454,p = .8728).Interpolation of the most informative sPCA axis (17.9% of variance explained) revealed two major regions with distinct genetic diversity (north and south) in agreement with the structure recovered by our haplotype network and phylogenetic tree (Figure2).AMOVA indicated this partition accounts for 79.92% of total genetic variability, with significant differentiation between north and south lineages (p = .001).

(
Φ ST > 0.459, p < .05; Figure S1.3 in Appendix S1), consistent with limited gene flow.Mismatch distributions for the whole dataset and the two main lineages were multimodal, a signature of demographic equilibrium.Tajima's D and Ramos-Onsins & Rozas' R 2 values agreed with the inference of population stability, however, Fu's F S test was significantly negative for the whole dataset (−18.143,p = .026)and both lineages (northern = −20.120,p = .011;southern = −11.631,p = .003),implying an excess of alleles as expected from recent population expansion.This inference was supported by non-significant Harpending's r values in all cases.When analysing individual population samples (n > 5), compelling statistical support for expansion was found for five populations in central North Island and nearby areas:

Figure 3b )
Figure 3b) as depicted by the best-fitted Gaussian model (univariate, equal variance model; BIC = −302.89).Overall, females in the north tend to have more triangular and waisted pronota than those in the Evaluation scores for individual models indicated good to excellent performance (Figure S3.3 in Appendix S3), and the ensemble model had high predictability, with a mean AUC score of 0.996.Sensitivity F I G U R E 2 Relief map of North Island, New Zealand with sampling locations for genetic analyses of Sigaus piliferus (coloured by geographic region).Phylogenetic (left) and genealogical (right top) relationships among mtDNA ND2 haplotypes revealed a deep north-south phylogeographic split and phylogeographic structure among northern populations.This major structure was supported by interpolation of the first global axis of a spatial PCA (right bottom).A star-like pattern of haplotype diversity (black arrowhead) was apparent in the Central Plateau area.Collapsed phylogenetic clades (posterior probability ≥ .95) are shown as triangles with colours indicating the proportion of terminals per geographic region (complete information on relationships and node support are available in Figure S1.2 in Appendix S1).Populations with statistical support for demographic expansion are shown.Map projection: NZGD2000.(i.e.fraction of presences correctly predicted) and specificity (i.e.fraction of absences correctly predicted) were also high, 100% and 99.6%, respectively.The cut-off threshold that maximised the proportion of presences and absences correctly predicted by the model was 488.5.Precipitation and temperature related variables had the largest influence on the suitable envelope estimated for S. piliferus.Precipitation seasonality, precipitation of driest month, and mean temperature of driest quarter were the most important predictor variables (importance score ≥25% in all cases; Figure S3.4 in Appendix S3).Current predictions closely matched the known range of S. piliferus throughout most elevated areas in North Island (19,723 km 2 ), except those in Taranaki and Northland districts, where the species has not been recorded (Figure 4; Figure S3.5 in Appendix S3).
cal directions: southern-northern (p = .045)and northern-southern comparisons (p = .048).No evidence of niche divergence was detected in any direction from the climate model (p ≥ .885).Niche specialisation and marginality were higher for the southern (S = 13.895,M = 3.993) than for the northern lineage (S = 1.673,M = 3.538), indicating that populations in the south have narrower climate niches that greatly differ from the average climate conditions in the island, hence, more specialised climate requirements.

F
I G U R E 3 Relief maps of North Island, New Zealand showing the spatial distribution of phenotypic assignment inferred by Gaussian clustering for adult female Sigaus piliferus based on body size (a) and pronotum shape (b).Pies represent population samples scaled by the number of individuals in each cluster.Phenotypic mixing is apparent in the species' mid-range (central North Island).Stacked bar plots (bottom) show assignment probabilities to each morphological cluster for each individual arranged by latitude.Canonical variate frequencies (right) supported group clustering based on Gaussian modelling.Map projection: NZGD2000.
1 in Appendix S3) indicates limited niche sharing (a); top and right plots show the density distribution along each axis.Suitable space within the available background across the marginality (x) and specialization (y) axes for the northern (b) and southern lineage (c); the centroids of the used space and the projections of the predictors are shown.Grasshoppers within the northern lineage (d) typically occur in montane woody vegetation (e.g.The Pinnacles, Coromandel Range), while those within the southern lineage (e) are found on open alpine grasslands (e.g.Herepai Peak, Tararua Range).Silhouettes indicate the typical morphology of grasshoppers (to the same scale) found in each type of habitat.TA B L E 1 Final multiple matrix regression with randomization models obtained through backward stepwise selection procedure indicating factors contributing to phenotypic differentiation among populations of Sigaus piliferus.
scenario consistent with reduced and fragmented habitat during glacial phases.Overlying the projected ecological response to changing climate we found evidence of demographic disruption from volcanic activity in the central part of the species' distribution.Contraction of suitable space for S. piliferus at the last glacial maximum challenges expectations for an alpine specialist in this region where palaeoecological evidence indicates unglaciated alpinelike plant communities dominated the landscape(McGlone, 1985; the LGM were largely overlapping, which limits confident separation of their influence.Increased interglacial habitat for S. piliferus would have enabled the range shifting that is apparent from demographic inference of population expansion and weak or no signature of isolation-by-distance across this species' range.Extant populations from areas less affected by volcanism towards the edges of the range indicate contrasting evolutionary histories.Niche models indicate persistence of suitable conditions in the north at the LGM (Coromandel-Kaimai ranges), and gradual contraction since then.This agrees with demographic inferences from the sole population sample (n = 10) from this region (Te Aroha) indicating population stability.In contrast, a fragmented LGM habitat in the south (Tararua Range) appears to have coalesced and we found demographic signatures of expansion in the best sampled (n = 12) population (Herepai) from this region indicating recent expansion from few ancestors.
locomotory morphology and vegetation type has been documented in European grasshoppers(García-Navas et al., 2017b) and suggests a common pattern of adaptive ecotypic divergence.While gene-environment interactions often underlie clinal phenotypic variation, such as ecotypic divergence in body size and shape (García-Navas for granting access and approval to collect from the Conservation Estate in North Island (Permits 37024-FAU, TW-32116-FAU, W/E 31465/FAU, ECHB-15515-RES, TT-15419-FAU).David Carmelet and Eliana Ramos provided valuable advice on shape analyses.We also acknowledge our field partners: Eliana Ramos and Mari Nakano.The R script for Figure5awas kindly provided by Javier Fernández-López and Irene Villa-Machío.Valuable comments and constructive criticism on earlier versions of this manuscript were provided by three anonymous referees.This research was supported by a doctoral scholarship from Massey University and a Massey University Doctoral Conference Grant (awarded to FLMJ).Open access publishing facilitated by Massey University, as part of the Wiley -Massey University agreement via the Council of Australian University Librarians.
. Extant populations of S. piliferus in the central volcanic available in Genbank at https:// www.ncbi.nlm.nih.gov/ genbank, accession numbers PP277795-PP277981.Phenotypic, genetic and distribution data are available at Dryad Data Repository (DOI: 10.5061/dryad.r4xgxd2m1).The climate data for ecological niche modelling and niche analyses are freely available at https:// www.world clim.org/ data/ world clim21.html.