Evidence for genetic isolation and local adaptation in the field cricket Gryllus campestris

Understanding how species can thrive in a range of environments is a central challenge for evolutionary ecology. There is strong evidence for local adaptation along large‐scale ecological clines in insects. However, potential adaptation among neighbouring populations differing in their environment has been studied much less. We used RAD sequencing to quantify genetic divergence and clustering of ten populations of the field cricket Gryllus campestris in the Cantabrian Mountains of northern Spain, and an outgroup on the inland plain. Our populations were chosen to represent replicate high and low altitude habitats. We identified genetic clusters that include both high and low altitude populations indicating that the two habitat types do not hold ancestrally distinct lineages. Using common‐garden rearing experiments to remove environmental effects, we found evidence for differences between high and low altitude populations in physiological and life‐history traits. As predicted by the local adaptation hypothesis, crickets with parents from cooler (high altitude) populations recovered from periods of extreme cooling more rapidly than those with parents from warmer (low altitude) populations. Growth rates also differed between offspring from high and low altitude populations. However, contrary to our prediction that crickets from high altitudes would grow faster, the most striking difference was that at high temperatures, growth was fastest in individuals from low altitudes. Our findings reveal that populations a few tens of kilometres apart have independently evolved adaptations to their environment. This suggests that local adaptation in a range of traits may be commonplace even in mobile invertebrates at scales of a small fraction of species' distributions.


| INTRODUC TI ON
Organisms can inhabit a broad geographic range, either because they remain competitive despite not having the optimal phenotype in parts of their range, or because they change as their environment changes. Generalists are able to compete effectively for resources in a variety of environments while expressing a single phenotype.
However, better adaptation to particular regions can be achieved through the expression of different phenotypes in different environments. Environmentally adapted phenotypes can occur either through adaptive phenotypic plasticity, where the same genotype produces different phenotypes dependent on environmental conditions (Acasuso-Rivero et al., 2019), or through the evolution of geographically adapted phenotypes. Local adaptation (Williams, 1966) has been studied both for the opportunities it provides for the study of evolution in nature (Kawecki & Ebert, 2004), and more urgently, because of its relevance for predicting biotic responses to climate change (Franks & Hoffmann, 2012;Schilthuizen & Kellermann, 2014).
Understanding whether species with broad geographical distributions remain competitive across a range of environments through a generalist phenotype or through specialization derived from phenotypic plasticity or local adaptation is central to understanding present and future species distributions. In this context, local adaptation is particularly important because the expected ability to cope with climate change of species with large geographical ranges (Thomas et al., 2004) requires that species-level tolerances for variation in environment are also present within individual populations (Kelly et al., 2011).
Local adaptation has been thoroughly studied in plants because of their stationary habit. The overall picture in these studies is that local adaptation is prevalent, although not ubiquitous (Leimu & Fischer, 2008;Savolainen et al., 2007). Like plants, small ectotherms such as insects are acutely exposed to climatic variables and do not typically use metabolism to maintain their temperature. However, unlike plants, they have the potential to use behaviour to modify their exposure to the environment by moving (Kearney et al., 2009;May, 1979;Ørskov et al., 2019). Most studies of local adaptation examine changes across large-scale environmental clines. There is abundant evidence for changes in adaptations to climate across these clines in animals (Keller et al., 2013), including examples in insects revealed by common garden rearing or reciprocal transplant experiments (e.g. Bradford & Roff, 1995;Demont & Blanckenhorn, 2008;Ottenheim et al., 2008). Studies examining heritable variation in lifehistory traits across altitudinal gradients, even at scales of tens of km, have also found evidence for local adaptation (Hodkinson, 2005) (for instance in development rate in grasshoppers (Berner et al., 2004;Laiolo & Obeso, 2015)).
Crickets have a long history in the study of local adaptation.
Masaki's studies of the field cricket Teleogryllus emma revealed heritable differences in traits including the temperature dependency of diapause, development rate and adult body size, between northern and southern populations (Masaki, 1963). This variation was shown to form a cline across ~2,000 km of their range in Japan (Masaki, 1967). Subsequently, similar clines have been found in other crickets (Bradford & Roff, 1995;Matsuda & Numata, 2019;Nystrand et al., 2011;Shimizu & Masaki, 1997). More recently, the discovery of local adaptation in field crickets to acoustically orienting parasitoids, Teleogryllus oceanicus on Hawaii (Zuk et al., 2006) provides an example of very rapid local adaptation in an isolated population (Pascoal et al., 2020). Studies to investigate the size of the fitness benefits that local adaptation provides in nature will be valuable in the future (Keller et al., 2013), but measuring individual fitness in the wild is very demanding. However, previous work on the European field cricket Gryllus campestris has demonstrated that they can be individually phenotyped and genotyped in their natural habitat and fitness measures extracted (Bretman et al., 2011;Rodríguez-Muñoz et al., 2010;. Furthermore, they can be translocated from one burrow to another without any detectable effects on their behaviour (Niemelä & Dingemanse, 2017, pers obs). Gryllus campestris is flightless and its rough grassland habitat is highly fragmented by anthropogenic (intensive agriculture, towns, roads etc.) and natural (rivers, forests, other unsuitable habitat) barriers. Local populations of these crickets are found from sea level to over 1,300 m (Table 1), but whether the species is locally adapted to high and low altitude habitats is unknown. Our aim was to investigate the potential for local adaptation in this model species on a scale of a few tens of km. To do so, we examined whether populations close to the north coast of Spain form genetic clusters according to whether they are at high or low altitude, or whether there is the potential for independent evolution of adaptations to altitude in replicate populations situated either close to sea level or at higher altitude. We then used common-garden rearing experiments to test for genetically determined differences in physiological and life-history traits.
We focus on chill coma recovery time (CCRT), as a proxy for cold tolerance and juvenile growth rate as a key life-history parameter.
CCRT is a widely used metric of cold tolerance in small ectotherms. It is measured as the time it takes to be capable of muscular activity after the reversible loss of motor capacity due to cooling. This elevated-temperature-induced recovery reflects a regaining of ion homeostasis after cooling (Findsen et al., 2014). The capacity to recover rapidly from cold-induced coma is likely to be particularly relevant in environments where overnight temperatures (or other short-term fluctuations) are cold enough to prevent activity. In the winter ant, Prenolepis impairs, CCRT is faster in high altitude populations: (Tonione et al., 2020). In Drosophila species, CCRT has been related to thermal adaptation as it differentiates between species of tropical and temperate origin (Gibert et al., 2001). CCRT responds K E Y W O R D S chill coma, climate adaptation, growth rate, phenotypic plasticity, RAD-seq to selection in D. melanogaster (Gerken et al., 2016), making it a candidate locally adapted aspect of the cold tolerance phenotype of insects.
Our predictions are that the loss of flight ability in this species and the fragmented nature of their contemporary distribution will mean that there are sufficient restrictions on gene flow to allow local adaptation among populations. There are substantial differences in both the biotic and abiotic environments between high and low altitudes. We expect that populations from higher altitude will show faster recovery from cold-induced coma since they are more likely to experience short-term periods when temperatures are too low for activity. We also expect that high altitude populations will develop faster because growth occupies the majority of the crickets' nondiapause life, and the period when it is warm enough for growth is shorter at higher altitude.

| Study system
We conducted the study over 4 years using G. campestris collected from 18 locations in Northern Spain (Table 1, Figure 1). These crickets have a single annual generation, and in this region, which is towards the southerly limit of their range, they live from sea level up to around 1,500 m. Winter snowfalls can occur at low altitudes and high altitude sites are covered in snow for several months in winter, during which time the crickets have an obligatory diapause as lateinstar nymphs. Summer temperatures can reach the high 30°C at low altitudes. The mean height (± SD) of low altitude populations is 108 ± 54 m and of high altitude populations is 1,160 ± 108 m. We do not have precise climatological data for collection sites, but the lower air density is expected to reduce the air temperature by ~7°C at high versus low altitude sites. Mean daily temperatures in 2018 at a site at 127 m (Asturias airport (N 43 33.74, W6 02.00)) and one at 1,480 m (Valgrande ski station (N 42 58.89, 5 46.54)) were 14.2 and 7.2°C respectively (datosclima.es, 2020) which suggests a temperature difference between our high and low altitude sites in the region of 5-6°C. The first adults appear in late April or early May at low altitudes and a few weeks later at high altitudes. The longest lived adults live until late June.
In 2016, we collected crickets from 11 wild populations including the meadow, where a long-term study ("WildCrickets") measuring individual fitness in the wild has been carried out for over a decade (Rodríguez-Muñoz, Boonekamp, et al., 2019;Rodríguez-Muñoz et al., 2010). Ten crickets from each population were used for a RADseq population genetics study. In light of the population structure, we identified in this study (below), in 2018 and 2019 we collected crickets from 3 of the same populations (1/2 black/white circles) and from 7 new populations (black circles) (Figure 1).
There is widespread evidence for local adaptation at large spatial scales, particularly along large-scale ecological clines (Keller et al., 2013). Habitat variation frequently occurs at much smaller scales, but studies of local adaptation within patchworks of Note: León was included as an outgroup, being located on the inland plain on the southern side of the Cantabrian range whereas the other populations are all on the northern side of the range.
TA B L E 1 Sampling locations (see Figure 1). 2016 samples were used for the RAD-seq analysis, 2018 for the chill coma recovery experiments and 2019 for the growth rate experiment neighbouring habitats are scarce. Our collection sites were chosen with the aim of identifying populations at either high or low altitudes that are nevertheless in close proximity to one another. To avoid potentially closer ancestry within altitude treatments than between them, we chose 5 low and 5 high altitude populations each of which is isolated from other populations in the same altitude treatment by at least one major river, but is in the same water catchment as a population in the opposite altitude treatment. There are numerous barriers to dispersal for these flightless insects, but large rivers are likely to be a very difficult barrier for crickets to cross. Due to an issue with access, site H9 included 2 meadows 1 km apart (see Table 1); however, as there is a lot of suitable habitat between the two H9 collections, we expect them to represent a single deme.
All study animals were collected in late April and early May using burrow traps (Meadows, 2013) and taken to our facility at site L15.
There is a correlation between altitude and latitude due to the East-West orientation of the Cantabrian Mountains ( Figure 1), but it is modest, and the climatic effects of the change in altitude are very much larger than any potential latitudinal effects.

| RAD tag sequencing
We sampled 10 adult crickets from each of the 11 populations sam- Raw read processing was done as in Settepani et al. (2017) to obtain data for further analyses. In brief, raw reads (forward and reverse) were concatenated or overlapped before quality filtering and demultiplexing. If paired reads were discarded in the quality filtering step, it was often due to poor quality of the reverse read. In those cases, we used the forward read for further analyses if it passed quality filtering. See Settepani et al. (2017) for further details.
We used pyRAD v2.15 (steps 3-7) (Eaton, 2014) to analyse RAD tag data using the following parameter settings: 'clustering thresh- We used Python scripts, also adapted from Settepani et al. (2017) to further process the output of the pyRAD analyses to estimate nucleotide diversity (π) (Tajima, 1983) and genetic structure. To perform population genetic analyses, it was necessary to remove any bias arising from the difference in RAD loci length. We therefore concatenated accepted loci into one alignment in a random order before analysis. For all diversity estimates, we split each sample into 20 subsets of equal size and estimated the population genetic parameter on each subset. This allowed us to obtain estimates of certainty on point estimates; here confidence intervals. All confidence intervals were obtained by bootstrapping over the 20 separate estimates and were used for statistical comparisons. We used Mega6 (Tamura et al., 2013) to estimate nucleotide diversity (π).
F I G U R E 1 Locations of cricket populations near the north coast of Spain.
Populations prefixed with an 'L' are at low altitudes, populations prefixed with an 'H' are at high altitudes (see table 1).
Populations sampled for our RAD-seq study in 2016 are shown as ○ or ◐ or as * for the WildCrickets meadow (Rodríguez-Muñoz et al., 2010). Populations sampled for the common-garden rearing studies of phenotypic traits in 2018 and 2019 are shown as • or ◐. Grey shaded areas are land above 1,600 m, and lines south of the coast are major rivers Genetic differentiation (F ST ) was estimated by 1−π P /π S with π P being the average nucleotide diversity of the populations under consideration and π S the total nucleotide diversity of the species. We tested for isolation by distance (correlation between genetic distance and geographical distance) using Mantel tests with 999 permutations (Oksanen et al., 2019), after excluding the outgroup Léon populations. Population structure analyses were performed in fastSTRUC-TURE (Raj et al., 2014), with K from 1-5. Graphical representation of fastSTRUCTURE results were obtained with a custom script (available here: https://github.com/sheng lin-liu/plotF astSt ructure).

| Chill coma recovery assays
In spring 2018, we collected ten adult or last-instar crickets of each sex from each of populations L1-L5 and H6-H10 and took them back to the lab at the University of Exeter's Penryn campus, UK.
Crickets were maintained in individual 7 × 7 cm boxes provided with ad-lib food (standard rodent diet) and continuous access to water through a 5 ml glass vial with a cotton wool stopper in each box.
Boxes were kept in a single temperature-controlled room at 25°C with a 7 a.m.-10 p.m. light cycle reflecting the contemporaneous day length in Asturias. Once adult, females were each mated to a male from their own population, by placing a male in the female's box and observing them until mating occurred. The position of the boxes was rearranged at least weekly to ensure that potential temperature gradients within the room could not systematically affect populations.
Once females from all populations had been adult for at least 3 days, and at least 5 females per population had mated successfully, all females were provided with wet sand for egg laying. It is possible that a small proportion of females had mated in the field prior to their laboratory mating. Because females store sperm, this means that some families that we subsequently assume to be full-siblings will include some half-siblings; however, this is conservative in relation to all our analyses. Females were left for a period of 48 hr before eggs were collected and transferred to a petri dish containing cottonwool saturated with water and maintained at 25°C in an incubator (1 petri dish per female). Eggs were spread out across the petri dish to reduce the potential for mould to spread between them. As nymphs hatched they were transferred to boxes of the same design as those We used a thermometer to monitor the water temperature, adding ice to the water bath if required to maintain the temperature between 0 and 1°C. After 1 hr the vials were removed from the ice bath and crickets tipped onto their backs into individual petri dishes in a 25°C CT room with typical indoor laboratory/office lighting levels. If any individual(s) landed the right way up, a small paintbrush was used to flip them over (if any nymphs needed to be flipped over, all nymphs in that trial were also brushed). The time taken for each individual to return to an upright position was measured by two observers watching 5 crickets each and recording the time to the nearest second using a voice recorder. Observers were blind to the altitude of the populations that were tested. The process was repeated over 5 days measuring 448 nymphs between 4 and 14 days old. To increase our sample size, the procedure was repeated with 300 nymphs of 21 to 35 days old and with 205 nymphs of 109 to 120 days old using unique nymphs for all assays. The third assay was conducted in our lab at the University of Aarhus (because of the availability of personnel), using identical procedures to those in the UK apart from slightly greater variation in the temperature of the room where the assays were carried out, which ranged between 25-28°C in Aarhus. Although the experiment was divided into 3, this was in order to increase our sample size, and not because we had any specific hypotheses about developmental stage and CCRT. Because we weighed crickets at these different points in their development, we acquired some data about growth rates. Mass and age are highly correlated, and hence, it does not make sense to include both of them in an analysis. We chose to use mass rather than age in our analysis of effects of altitude on chill coma recovery times (CCRT) because mass directly affects rates of heat loss and gain. However, because crickets for these assays were not kept individually and the density in each box changed frequently due to mortality we cannot use these data to examine growth rate reliably. Therefore, in the following year we carried out the specifically designed growth rate assays described below.

| Growth rate assays
In spring 2019, we collected 440 wild field crickets Gryllus campestris (198 males and 242 females) from the same populations sampled in 2018, L1-L5 and H6-H10. Crickets were flown back to our Penryn laboratory and the same husbandry and breeding protocol was followed as described above for our chill coma recovery assays. However, whereas in the CCRT study we were unable to follow the parentage of each individual (because crickets were kept in groups), in this study we kept every nymph in an individual cage, allowing us to retain information about parentage. Nymphs began to hatch on the 22 May 2019, and we continued to collect them until we reached our target of 1,000 nymphs with a maximum of 10 hatchlings from any individual female. Each nymph was placed into an individual 7 × 7 cm plastic box with a perforated lid for aeration, and a piece of cardboard egg carton as a substrate. They were provided with a water vial (as described above) and an excess of a mixture of two ground up foods ('Pets at Home Rat Nuggets' (protein 16%, crude fibre 4%, crude fat 5%, crude ash 5%) and 'Burgess Excel Nature's Blend Rabbit Nuggets' (protein 12.6%, crude fibre 19%, crude fat 3.6%, crude ash 6.5%)). Each nymph box was labelled with only a number to ensure the study was conducted blind and was allocated to one of two temperatures 23°C or 28°C. These temperatures were chosen because our pilot studies showed very low growth rates below 21°C and very high mortality above 30°C. Also, the 5 degree temperature difference is in the region of the 7 degree difference that we expect to find between high and low altitudes due to the difference in air density. Boxes were distributed among 2 incubators at each temperature with L:

| Growth rate statistical methods
We analysed the effects of age and temperature on body mass using a bivariate character state model approach (Via & Lande, 1985;Wilson et al., 2010 (see tutorial 2)). This model is mathematically identical to a univariate model with a random slope for temperature or age, but the bivariate model estimates the residual variance for each of the two states individually, whereas the univariate model estimates a common residual variance. Because the variance of body mass can be expected to be highly heterogeneous among ages for growing animals, it would be problematic to assume a symmetrical common residual variance among age groups. We therefore preferred using a bivariate model. We included a pedigree random effect based on assumed full-sib relationships to partition variance into heritable (genetic and maternal) and nonheritable (environmental) effects. We treated mass as a bivariate response variable at the two ages of day 0 (mass on the day of hatching) and day 39. The intercepts of these two variables estimate the average growth between day 0 and day 39 (we also show the measurements made on days 13 and 26 in Figure 6, but these were not included in the analysis). The pedigree random effect estimates the variance due to the combined genetic and nongenetic parental effects. We assume that parental effects are relatively small given the small size of eggs and lack of any parental care, and the relatively small spermatophore, which females do not consume. However, we cannot rule out paternal and maternal effects, and previous studies on field crickets indicate that not only parental genetic but also epigenetic factors can be transmitted to the germ cells affecting their biochemical composition (Marshall et al., 2009). It is currently unknown whether such nongenetic inheritance affects offspring development and there are fiercely divergent views on the implications for quantitative genetic analyses. Our study design lacks the pedigree complexity that is required for estimating epigenetic variance components (Nadeau, 2009;Tal et al., 2010). Given these limitations, we assume that the additive effects of nongenetic heritable factors are relatively small compared to the effects of additive genetic variance, but there is a possibility that this assumption is wrong. Likewise, dominance and epistatic genetic effects have been generally assumed to be proportionally small (Falconer & Mackay, 1996) and we therefore assume that the genetic variance estimated by the pedigree primarily reflects additive genetic variance. Hence, we use the additive genetic (V A ) and the residual (V R ) variance components to estimate the narrow sense heritability of mass at t(0) and t(39), i.e. h 2 (t) = V A(t) /(V A(t) +V R(t) ). However, we acknowledge that our estimate of additive genetic variance could be confounded by the exclusion of dominance, epistatic and nonheritable factors from the model assumptions, potentially resulting in an overestimation of the additive genetic variance and narrow sense heritability.
We also estimated the genetic correlation among the two age , where cov G is the additive genetic covariance between V A(0) and V A(39) . Genotype by environment (G × E) or Genotype by age (G × Age) interactions are observed when the genetic correlation is significantly below 1 and / or when the two genetic variances significantly differ among the two states (as above, assuming maternal effects are small). We used the coefficient of genetic variation (CV A ) for comparing genetic variances among character states to scale the genetic variance to the mean of character states, essentially controlling for growth. We calculated the evolvability by taking the square of CV A as described by Houle (1992) and Garcia-Gonzalez et al. (2012). Against this model background we tested the effects of altitude of the source population and rearing temperature and their potential interaction by including these variables as fixed factors. In addition, we tested whether there was a G × E interaction using a similar bivariate model but with same age nymphs (all day 39) reared at the two different temperature conditions. Note that the G × E model structure is similar to the G × Age model described above. We used ASReml v.4 in R for all analyses.

| Nucleotide diversity
We sequenced ~15 million bases per individual of which >110 K sites were polymorphic. There was no significant difference in nucleotide diversity between lowland and highland populations (F (1,8) = 3.3, p = .11); lowland populations on average have slightly higher population nucleotide diversity (π) (Figure 2). Two populations, H8 and H10 have lower nucleotide diversity than the other populations, whereas our outgroup, the remote southerly population (León), has the highest population nucleotide diversity (0.0039).

| Population structure
Our fastSTRUCTURE analysis shows that the most likely number of genetic lineages is three (K = 3; model selection based on deviance information criteria) (Figure 3 and Figure S1). H17 is in its own genetic lineage (orange), L11, L12 and H16 form a second genetic lineage ('Western' (blue)), and populations L13, L14, L15, H8, H9 and H10 form a third genetic lineage ('Eastern' (green)). It  (Tajima, 1983) per population estimated from RAD-seq data. Error bars represent 95% confidence intervals. Populations are grouped by altitude (L/H). Dotted horizontal lines represent averages for lowland and highland populations F I G U R E 3 Genetic structure: fastSTRUCTURE (Raj et al., 2014) was used to cluster individuals based on genetic information from RADseq data. The number of genetic clusters representing the lowest model complexity above which prediction errors do not vary significantly was 3 (K = 3). In the graphical representation, individuals are grouped by population (10 individuals per population). The colours represent each individual's membership of the three genetic lineages. See Fig. S1 for full results of K1-5 F I G U R E 4 Isolation by distance: pairwise genetic differentiation (F ST ) between populations as a function of the physical distance between them. The dotted line represents a linear regression of all pairwise populations excluding León (see results). Green dots represent divergence between populations both of which are in the green cluster (as identified in our structure analysis (Figure 3)) (15 pairs); blue dots are divergences within the blue cluster (3 pairs). '+'s indicate divergence between a population in the green cluster and one in the blue (19 pairs), and black circles indicate divergence between the orange population (H17) and blue or green populations (9 pairs). Squares indicate divergence between the León population and each of the other 10 populations structure within these lineages: K = 1 is the best model in both cases ( Figure S2). This pattern is consistent with three independent colonization events from a more southerly source population expanding towards the north coast, or that individuals from the three separate northern lineages migrated south. The León population, which is much further south on the landward side of the Cantabrian mountain range appears to be a mix of the three coastal lineages; however, there is a significant likelihood that this is an artefact of the fastSTRUCTURE analysis due to unbalanced sample sizes and/or because this population in genetically divergent (Lawson et al., 2018). We note that one individual from the H17 population seems to originate from the blue genetic lineage.
We cannot exclude that this is due to contamination, but it may also reflect a recent migration event. We retained this individual in all analyses.
If the southerly 'outgroup' population (León) is excluded, there is evidence of isolation by distance among the populations; this is illustrated by the distribution of data and linear regression line in Figure 4. However, it is more appropriate to analyse these data with a Mantel test (Séré et al., 2017). This confirms the pattern: (Mantel statistic based on Pearson's product-moment correlation (Oksanen et al., 2019) r: .574, p = .001). It is apparent that the orange lineage is the most divergent with average F ST to the green and blue lineages of 0.21, whereas average F ST between green and blue is 0.14. Average F ST within the green and blue lineages is 0.10 and 0.09 respectively (see Table S1 for all between population F ST values).
Chill coma recovery times (CCRT) were faster for nymphs from highland populations ( Figure 5). Our analysis used altitude and mass as fixed effects; we used mass rather than age because these two are highly correlated and mass directly affects rates of heat loss and There were differences in growth rate among populations, revealed by analysis of variance in individual mass at the times of CCRT assays, but this was not significantly related to differences between altitudes (Altitude: F (1,929) = 2.8, p = .09, Population within Altitude: F (8,929) = 5.5, p < .0001).

| Growth rate
After 39 days of growth, the expected large effect of temperature was evident ( Figure 6); there was also a difference in mass between crickets with parents from high and low altitude, which is apparent in crickets reared at the higher temperature ( Figure 6). Hatchling mass was significantly heritable (h 2 = 60%). Despite this high heritability, the coefficient of genetic variation (CV A ) was relatively low (CV A = 0.098) resulting in an evolvability of 0.96%. Hatchling mass was indistinguishable among the treatments, showing that there is no accidental initial size bias with respect to treatment. Source population altitude affected hatchling mass in the direction such that nymphs with parents from low altitude were heavier (Table 2). This effect of altitude was estimated over and above variation partitioned into the pedigree random effect, and hence, there is no evidence that the observed heritability is due to genetic differentiation with respect to altitude background. An additional experiment would be needed to rule out the possibility that there are maternal effects on hatchling size.
As expected, temperature treatment strongly affected growth rate. After 39 days, nymphs in the 28°C group were on average 4 times heavier than nymphs reared at 23°C. This strong effect of treatment presumably swamped the genetic-and altitude effects observable at day 0, because after 39 days of treatment the heritability was reduced to 28.5% and the main effect of altitude disappeared (Table 2; the main effect of altitude is reflected by the L23 estimate). Despite the lower heritability, the CV A was relatively high at 0.359 resulting in an evolvability of 12.9%. Despite the strong effect of temperature treatment, mass at day 0 showed a weak phenotypic correlation with mass at day 39 (r P = .23) ( Table 2). Note that the phenotypic correlation was weaker than the genetic correlation (r G = .56), which is a common phenomenon and is caused by the fact that the heritability of mass is substantially less than 100%; that is, when the heritabilities are exactly 1 then the phenotypic correlation approximates the genetic correlation, which is rarely the case (Cheverud, 1988). Source population altitude remained an important factor, because the effect of temperature treatment depended on altitude: Growth acceleration by temperature was more than 1.5fold greater for nymphs from low altitude populations compared to those with parents from high altitude populations (Table 2; Figure 6).
Growth rate (from day 0 to day 39) showed significant genetic variation as indicated by the substantially higher coefficient of genetic variation at day 39 (Table 2; Figure 6). If there were no genetic variation in growth then coefficients of variation would be of similar size among age groups. The higher coefficient of genetic variation at day 39 is over and above the direct effects of temperature and altitude. Hence, these results show that there is significant Genotype F I G U R E 5 Chill coma recovery time in five cricket populations from low and high altitude (means and standard errors). Horizontal bars represent the grand mean for high and low altitudes. Nymphs from each population were assayed at three developmental stages corresponding to 9 ± 5, 28 ± 7 and 113 ± 7 days post-hatching, data from these stages are pooled for this figure by Age variation in addition to effects of temperature and altitude on growth. High temperature during rearing and having parents from a low altitude population increase body mass, but source population altitude only affected body mass at day 0 and temperature only affected body mass at day 39. The effect of temperature depends on source population altitude as indicated by a significant interaction, resulting in the highest body mass in the low altitude at high temperature group (see Figure 6).
Lastly, we investigated whether in addition to the observed genetic variation in growth rate (G × Age), there was also evidence for genetic variation in the response to temperature (G × E) as such effects were reported in a similar experiment . However, the bivariate model (with body mass at the two temperature conditions after 39 days of treatment as the response variables) yielded unreliable results due to singularities and/ or estimates that were highly dependent on initial starting values of the ASReml optimization procedures. We attribute these errors to an insufficient statistical power to detect the G x E interaction that we previously observed, because (a) the previously estimated genetic variation in response to the same temperature treatment was only moderate despite being estimated after a longer period of temperature exposure (65 days) and (b) the sample size of our present study is substantially smaller.

| D ISCUSS I ON
Our aim was to investigate potential local adaptation in a typical flightless insect in populations separated by only a few tens of km.
Theoretical models predict that there will be selection for adaptation to altitude in life-history traits (Abrams et al., 1996), and it is clear that existing variation of this type has the potential to be TA B L E 2 Character state bivariate model estimating the effects of altitude (low vs. high) and temperature (23°C vs. 28°C) on mass at hatching and after 39 days of treatment Note: V A and V R denote the variance components of the pedigree and residual random effects respectively. CV A denotes the coefficient of additive genetic variation. h 2 is the estimated narrow sense heritability. r G is the genetic correlation among day 0 versus day 39 treatment groups. r P is the phenotypic correlation among the two ages. 'Estimate' reflects the calculated intercepts for all the four temperature and altitude groups, derived from the model estimates, and corresponds with the means plotted in Figure 6.

F I G U R E 6 Growth trajectories.
(a) Body mass in relation to age, altitude ancestry and rearing temperature (95% CI vertical bars represent 1.96 × SE estimated by the model shown in Table 2). (b) Day 0 mass of hatchlings plotted against their mass after 39 days. The phenotypic correlation between day 0 and day 39 mass is shown in Table 2 informative in relation to how ecosystems respond to climate change (Hodkinson, 2005). Our study was carried out on field crickets, not least because of their potential for future transplant experiments in which behaviour and life histories are monitored in detail at an individual level (Makai et al., 2020;Niemelä & Dingemanse, 2017;Rodríguez-Muñoz, Boonekamp, et al., 2019;Rodríguez-Muñoz et al., 2010). Crickets are also a typical temperate species with a single annual generation and a summer breeding season. Their dispersal is more restricted than that of flying insects and much larger vertebrates, but nevertheless, observations of local adaptation to altitude in G. campestris would provide a strong indication that similar adaptation may be widespread. In this scenario, rather than viewing local adaptation as a pattern occurring on scales of hundreds or thousands of km, it might be that many species have a patchwork of genetic adaptation mapping environmental heterogeneity.
If there is substantial gene flow between populations, then independent evolution will be less likely, even if the ~1,000 m difference in altitude between high and low altitude populations creates large environmental differences. Also, if geographically distinct high altitude populations represent a single deme, then our replicate sampled populations may not represent independent samples (and the same goes for our low altitude populations). Our RAD-seq study reveals that there is significant genetic differentiation among populations but not between high and low altitudes. Populations vary in their nucleotide diversity (Figure 2) with an average level of nucleotide diversity towards the lower end of values that have been reported for other arthropod species (Leffler et al., 2012;Lozier, 2014;Romiguier et al., 2014). Such low diversity is likely to at least partly be explained by inevitable underestimation caused by the RAD sequencing procedure (Arnold et al., 2013). However, it may additionally be the con- studies of Anuran populations (Bachmann et al., 2020), which concluded that gene flow was a significant inhibitor of local adaptation on steep environmental gradients.
Because G. campestris has an obligatory winter diapause, it is very challenging to rear multiple generations in a common environment to remove the potential for parents to transfer environmental effects to offspring. Our approach was to catch males and females just before, or just after they became adult and to keep them in a common environment for a period of from a few days to a few weeks before we began collecting eggs. Our chill coma recovery experiments provide convincing evidence that nymphs hatching from eggs laid by adults from high altitude recover from chill coma more rapidly than those from low altitude. The nymphs used in these assays are of a developmental phase that we would expect to find in the wild in the late summer. At low altitudes, the temperature will never fall below zero at this time of year, but overnight sub-zero temperatures do occur at high altitude, and become increasingly common later in the year. As discussed above, we cannot determine whether this is a genetic difference or a parental effect. The mechanistic basis for a parental effect on chill coma recovery time is unclear and we cannot find any published examples thereof, but further experiments would be needed to rule them out. The obvious interpretation of our findings is that crickets at high altitude have evolved physiological adaptations that allow them to recover from chill coma more rapidly. This matches the general pattern seen in inter-species comparisons in Drosophila, where species of tropical origin have longer CCRT than species of temperate origin (Gibert et al., 2001). CCRT seem to be an independent part of the cold tolerance phenotype (Gerken et al., 2016); for instance, Teets and Hahn (2018) found minimal overlap in genes associated with cold shock response and chill coma recovery in Drosophila. This suggests that the differences we observe between high and low altitudes are a direct effect of thermal adaptation, reflecting fitness benefits of being able to rapidly regain motion after extreme cooling events (Findsen et al., 2014).
Previous work  identified significant heritability of growth rate, again using nymphs reared in sibling groups. In our growth rate study, nymphs were reared individually to eliminate environmental effects resulting from competition (including faster growing individuals potentially having to compete with faster growing siblings). We found that newly hatched nymphs from parents from low altitudes were significantly heavier than those with parents from high altitudes. In contrast with inheritance of chill coma recovery (discussed above), variation in egg size provides an obvious potential mechanism for a maternal effect on hatchling mass. This could create an association between maternal environment and offspring mass at hatching Relationships between female body size and egg size have been identified in some studies of insects and have been found to be absent in others (references in Sturm (2016)). The only study on Orthopterans that we are aware of, found that in the cricket Teleogryllus commodus there is a strong relationship between female mass and fecundity but very little variation in egg size (Sturm, 2016). A study on the meadow grasshopper, Chorthippus parallelus, (Köhler, 1983) cited in (Berner et al., 2004) found negligible maternal influence on offspring embryonic development. These findings suggest that the differences we observe are more likely to be genetic than parental effects, but either remains possible. Our most striking result in relation to body mass and growth can be seen in the mass of nymphs at the end of our growth rate measurement period in Figure 6. The effect of temperature on nymph growth rate differed between the offspring of parents from high and low altitudes (as shown in our analysis).
At the lower rearing temperature, there was relatively little difference between high and low altitude crickets, but at the higher temperature, crickets from low altitude grew much faster than those from high altitude. This result is in the opposite direction of our expectation that high altitude populations would develop faster because of selection for rapid development as a result of the shorter breeding season at higher altitude. Optimisation models (Abrams et al., 1996) show that although the optimal response to a shorter period available for growth is typically faster growth, slower growth can also be favoured. Whether or not this occurs depends on factors such as the mortality rate and the relationship between adult size and fitness. Both of these are parameters that might vary substantially between high and low altitude sites, for instance because of differences in predator communities in the two types of environment. Another possible explanation for the faster growth of low altitude crickets at higher temperatures is that this reflects an evolved capacity of these populations to exploit high temperatures to grow more rapidly. This capacity may be weaker in high altitude populations because they do not experience such consistently high temperatures. The nature of the trade-offs that might drive the type of adaptation remain to be investigated. In general, constraints on growth rate are poorly understood (Dmitriew, 2011), with a likelihood that strongly environment-dependent factors such as predation risk are likely to be important in understanding why animals do not grow faster (Dmitriew, 2011). Overall, our RADseq study provides clear evidence of the potential for local adaptation on a scale of tens of km in a flightless, but otherwise mobile and widespread temperate insect. Our investigations of chill coma recovery time and growth rate in nymphs reveal what are likely to be heritable differences in these traits between populations at high and low altitudes, suggesting that local adaptation in these traits has occurred. Our findings suggest that local adaptation may be an important aspect of population viability in temperate insects. This has obvious consequences for how insect populations are likely to be impacted by climate change.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/jeb.13911.

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https://datad ryad.org/stash/ datas et/doi:10.5061/dryad.j9kd5 1ccx.