- Research article
- Open Access
Hierarchical genetic structure shaped by topography in a narrow-endemic montane grasshopper
BMC Evolutionary Biologyvolume 16, Article number: 96 (2016)
Understanding the underlying processes shaping spatial patterns of genetic structure in free-ranging organisms is a central topic in evolutionary biology. Here, we aim to disentangle the relative importance of neutral (i.e. genetic drift) and local adaptation (i.e. ecological divergence) processes in the evolution of spatial genetic structure of the Morales grasshopper (Chorthippus saulcyi moralesi), a narrow-endemic taxon restricted to the Central Pyrenees. More specifically, we analysed range-wide patterns of genetic structure and tested whether they were shaped by geography (isolation-by-distance, IBD), topographic complexity and present and past habitat suitability models (isolation-by-resistance, IBR), and environmental dissimilarity (isolation-by-environment, IBE).
Different clustering analyses revealed a deep genetic structure that was best explained by IBR based on topographic complexity. Our analyses did not reveal a significant role of IBE, a fact that may be due to low environmental variation among populations and/or consequence of other ecological factors not considered in this study are involved in local adaptation processes. IBR scenarios informed by current and past climate distribution models did not show either a significant impact on genetic differentiation after controlling for the effects of topographic complexity, which may indicate that they are not capturing well microhabitat structure in the present or the genetic signal left by dispersal routes defined by habitat corridors in the past.
Overall, these results indicate that spatial patterns of genetic variation in our study system are primarily explained by neutral divergence and migration-drift equilibrium due to limited dispersal across abrupt reliefs, whereas environmental variation or spatial heterogeneity in habitat suitability associated with the complex topography of the region had no significant effect on genetic discontinuities after controlling for geography. Our study highlights the importance of considering a comprehensive suite of potential isolating mechanisms and analytical approaches in order to get robust inferences on the processes promoting genetic divergence of natural populations.
Understanding the factors structuring genetic variation in natural populations is a paramount topic in evolutionary biology [1, 2]. The genetic structure of populations is primarily determined by inter-population dispersal rates and realized gene flow, which in turn are shaped by geography, environment, historical processes and, more frequently, their combined effects [3–5]. The isolation-by-distance (IBD) model predicts that genetic differentiation increases with Euclidean geographic distance because of limited dispersal and genetic drift [6, 7]. Even though this classic model explains spatial patterns of genetic differentiation in a wide range of organisms (, but see review in ), it does not consider more sophisticated information than straight-line geographic distances and assumes landscape homogeneity, an unrealistic scenario for most natural systems [9, 10]. Recently, the emergence of landscape genetics has explicitly incorporated landscape complexity into the study of evolutionary processes [9, 10], bringing new approaches that take into account the ability of organisms to disperse across different landscape features according with the resistance that they offer to movement (i.e. isolation-by-resistance, IBR [11, 12]; e.g., [13, 14]).
Beyond geography and the spatial configuration of connecting corridors and isolating barriers, environment can also play a major role in shaping spatial patterns of genetic differentiation [4, 15]. This occurs when populations inhabiting ecologically dissimilar habitats experience limited gene flow due to the low performance of immigrants arriving to areas where they may not be locally adapted or as consequence of the reluctance of individuals to cross or establish in unfamiliar habitats (i.e. isolation-by-environment, IBE [15–17]). Recent research has revealed the ubiquity of IBE patterns, highlighting the importance of ecological factors in shaping genetic structure of natural populations (see meta-analyses in [4, 18]). The contribution of environment relative to geography in explaining spatial patterns of genetic variation may vary depending on species characteristics and ecological features of the natural systems [18, 19]. Indeed, different isolating mechanisms are not mutually exclusive and gene flow is often influenced by a combination of geographical and ecological factors [5, 20, 21]. Given that geography and environment are often highly correlated, disentangling their relative impact on population genetic differentiation harbors inherent analytical difficulties (see  for a discussion about eco-spatial autocorrelation) that have promoted the progressive development of more robust and accurate statistical methods [23–27]. However, only a few studies have jointly considered the relative effects of geography, landscape composition and environmental heterogeneity on either landscape-level [28, 29] or range-wide patterns of genetic structure .
The Pyrenean Morales grasshopper (Chorthippus saulcyi moralesi) (Orthoptera: Acrididae) is a narrow-endemic subspecies belonging to the Chorthippus binotatus group, exclusively distributed in central Aragón and Catalonia Pyrenean mountains (see , for taxonomic status and detailed description). It is a winged grasshopper primarily feeding on gramineous herbs  and patchily distributed across a gradient of habitats, including submediterranean shrub formations, mesophile grasslands, montane shrubby vegetation and subalpine open grasslands located at elevations ranging between 1100 and 2400 m.a.s.l. [30, 31]. Its distribution range is restricted to a narrow longitudinal axis with an east–west orientation characterized by a gradient of precipitation and temperature, from Atlantic to Mediterranean climate regimes . The Pyrenees constitute the natural northern border of the Iberian Peninsula and present a high topographic and environmental complexity, rich biodiversity and considerable number of endemic species . These mountains experienced dramatic climate fluctuations during the Pleistocene , which are expected to have strongly influenced the demographic history and altered the distribution of many organisms of the region . Despite the wide variety of habitats and altitudinal ranges occupied by the Pyrenean Morales grasshopper, populations at elevations lower than 1400 m and above 2100 m are anecdotal [30, 31]. This suggests that valley bottoms and mountain tops, together with the complex topographic complexity of the region, may act as barriers to dispersal in this species [35, 36]. Thus, our study system has a great potential to examine the relative role of geography, environmental heterogeneity and present and past configuration of suitable habitats (i.e. corridors and barriers to dispersal) in shaping patterns of genetic differentiation throughout the entire distribution range of a narrowly distributed taxon . We first analyzed spatial patterns of genetic structure and then employed a suite of complementary statistical approaches to test three different plausible scenarios of population differentiation (see Fig. 1 for a summary of the workflow employed in this study). In particular, we used multiple matrix regressions with randomization (MMRR, ), distance-based redundancy analyses (dbRDA, ) and geostatistical modelling based on Bayesian inference  to test whether the spatial pattern of genetic differentiation in the Pyrenean Morales grasshopper is explained by i) geographic distances (IBD), ii) resistance distances based on topographic complexity and current and past (last glacial maximum and last interglacial) climate suitability (IBR); and iii) altitudinal and environmental dissimilarity between populations (IBE) (Fig. 1). If geography, topography or corridors/barriers defined by habitat suitability are identified as the major drivers of genetic structure, then migration-drift equilibrium and neutral divergence can be regarded as the main evolutionary force shaping genetic discontinuities . A predominant or independent significant effect of environment on disrupting gene among populations would point to a role of ecological divergence and local adaptation processes in the evolution of spatial genetic structure .
Between 2012 and 2014, we collected 202 individuals from 11 populations of the Pyrenean Morales grasshopper. We aimed to sample populations throughout the entire distribution range of the species (~7 000 km2; Fig. 2a) based on occurrence-data available in the literature [30, 31] and our own prospection of areas with potentially suitable habitats. The Morales grasshopper is primarily distributed in the south side of the Pyrenees (Spanish Pyrenees) and a small area located in the northeastern part of these mountains (French Pyrenees). The latter portion of the species distribution was intensively prospected during our field work but we only found a single population in the area (Err, France; Table 1). Overall, we were able to collect individuals from populations located across almost the entire climatic and elevation gradient occupied by the species (Table 1; Additional file 1: Figure S1), including populations from a wide range of habitats (from submediterranean shrub formations to subalpine grasslands) and differing in up to ~ 900 m of elevation (Table 1). We preserved specimens in 2 ml vials with 96 % ethanol and stored at − 20 ° C until DNA extraction. Population codes and more information on sampling sites are in Table 1. Specimens were collected in public lands under license from ‘Gobierno de Aragón’, ‘Generalitat de Catalunya’ and ‘Ordesa y Monte Perdido National Park’.
Microsatellite genotyping and basic genetic statistics
We extracted genomic DNA from a hind-leg of each individual using a salt extraction protocol . We amplified and genotyped each individual using the 18 microsatellites markers described in . We performed PCR amplifications following the procedure described in , run PCR products on an ABI 310 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA) and scored genotypes using GENEMAPPER 3.7 (Applied Biosystems, Foster City, CA, USA).
We tested for deviations from Hardy-Weinberg Equilibrium (HWE) using exact tests  based on 900 000 Markov chain iterations and 100 000 dememorization steps as implemented in the program ARLEQUIN 3.5 . We discarded two loci from all downstream analyses because of heterozygosity departure from HWE in all populations, probably due to the presence of null alleles according to MICRO-CHECKER analyses . We also used Arlequin to test for linkage disequilibrium using a likelihood-ratio statistic, whose distribution was generated with 10 000 permutations. We did not find evidence for linkage disequilibrium between any pair of loci in any sampling population after sequential Bonferroni corrections .
Genetic structure analyses
We estimated population genetic differentiation calculating F ST-values between all pairs of sampling populations and testing their significance with Fisher’s exact test after 10 000 permutations using ARLEQUIN 3.5. We corrected P-values using a sequential Bonferroni adjustment . Due to the high frequency of null alleles in Orthoptera (e.g., [14, 46]), we also calculated pairwise F ST-values corrected for null alleles (F STNA) using the so-called ENA-method implemented in the program FREENA [ 47 ].
We inferred genetic structure using Bayesian clustering analyses in STRUCTURE 2.3.3 [48, 49]. We performed an iterative approach to assess the hierarchical genetic structure that could underlie broad genetic clustering patterns identified by STRUCTURE analyses including all populations (for a similar approach, see ). After an initial global analysis including all populations, we analyzed subsequent subsets of the data corresponding to the respective genetic clusters identified in the previous hierarchical level. For all analyses, we considered correlated allele frequencies and an admixture model without prior information on population origin. We performed ten independent runs for each value of K (1 to 12 for the complete dataset; and 1 to n + x for reduced datasets of n populations, being x a number to achieve at least three ΔK values) with a burn-in period of 200 000 steps and a run length of 1 000 000 Markov chain Monte Carlo (MCMC) cycles. We estimated the best-supported number of genetic clusters with the log probability of the data [Ln Pr (X|K)]  and the ΔK method . We used the ‘full search’ algorithm in the program CLUMPP 1.1.2 [ 52 ] to align replicated runs and average individual assignment probabilities for the most likely K-value. Finally, we used DISTRUCT 1.1  to produce bar plots displaying probabilities of individual membership to each inferred genetic cluster.
Complementary to the Bayesian clustering method, we used a discriminant analysis of principal components (DAPC) to identify clusters of genetically related individuals . Although both clustering methods exhibit similarities (e.g. they do not require a priori delimitation of populations), several important differences exist in their analytical approaches. STRUCTURE suffers from the assumption of HWE and gametic disequilibrium  and typically fails to detect isolation-by-distance (IBD)  and hierarchical patterns of spatial genetic structure . However, the multivariate analyses implemented in DAPC do not lay on the assumptions of STRUCTURE Bayesian analyses and could be more efficient to detect complex patterns of genetic differentiation . DAPC is a methodological approach that requires data transformation using a principal component analysis (PCA) as a prior step to a discriminant analysis (DA). DA partitions genetic variation maximizing differences between clusters while minimizing within-cluster variation. We implemented DAPC analysis in R 3.1.2 using the package ADEGENET [54, 55]. At first, we used the ‘find.cluster’ function using all available principal components (PCs) to determine the best-supported number of genetic clusters using the Bayesian Information Criterion (BIC). The ‘find cluster’ function runs successive K-means clustering with increasing number of clusters (K) and provides a BIC value for each simulated K-value (i.e. K-value with the lowest BIC is the ‘optimal’ number of clusters). Secondarily, we determined the optimal number of PCs for the DAPC by cross-validation using the ‘xvalDapc’ function with 100 replicates. We selected the number of PCs associated with the lowest ‘root mean squared error’ (RMSE) value. We ran DAPC using all the available discriminant functions and calculated the assignment probability of individuals to each cluster, which were represented with barplots using DISTRUCT.
We constructed a phylogenetic tree to evaluate the genetic relationships among all populations. We used the program POPULATIONS 1.2.31  to obtain a neighbor-joining (NJ) tree based on pairwise Cavalli-Sforza and Edwards (D c ) genetic distances . This genetic distance is the most accurate to yield the correct tree topology for microsatellite markers under a variety of evolutionary scenarios without making assumptions in relation to mutation rates or constant population sizes .
Climate niche modelling
We used climate niche models (CNMs) at different time periods to investigate whether current and past climate suitability are relevant factors shaping observed patterns of genetic differentiation among populations of the Pyrenean Morales grasshopper. For this purpose, we built a CNM using the ‘maximum entropy presence-only’ algorithm implemented in MAXENT 3.3.3 [59, 60] based on current climate layers and using 50 cross-validation replicate model runs. We used a total of 47 occurrence points obtained from the literature [30, 31] and our own sampling. To construct the models, we used the 19 present-day bioclimatic variables available in WorldClim  and downloaded at 30 arc-sec (c. 1 km) resolution . To obtain the distribution of the Pyrenean Morales grasshopper in the Last Glacial Maximum (LGM, c. 21 kya BP) and the Last Interglacial (LIG, c. 120–140 kya BP), we projected contemporary species-climate relationships to these periods. The LGM layers were based on the Community Climate System Model (CCSM3, ) from the Palaeoclimate Modelling Intercomparison Project Phase II (PMIP2, ). We downloaded LGM layers from WorldClim at 2.5 arc-min and interpolated to 30 arc-sec resolutions. LIG layers were based on  and downloaded from WorldClim at 30 arc-sec resolution. According to the suggestions from , we limited the geographic extent of the climate layers to an area approximately 20 % larger than the known distribution range of the species in order to avoid model over-fitting. We used multivariate environmental similarity surfaces (MESS) calculation to address the problems derived from projecting the current distribution into novel climates (i.e. LGM and LIG periods) . We used this approach to identify and discard climate layers with areas where the predictions should be treated with caution, due to the variables are outside the range present in the training data (for more details see ). We carried out MESS analyses iteratively excluding one variable in each step until discarding all out of range LGM and LIG variables compared to present-day variables. Finally, we checked the reliability of our past and current climate models following two approaches. At first, we developed a new current climate model using the variables with greater 5 % importance (as selected by Jackknife of regularized training gain procedure) and we compared its similarity with the current model generated by MESS analyses. Second, we developed a new model including only the most informative variable (Bio 1) and we compared its LGM and LIG projections with those obtained by MESS analyses (for a similar approach, see [69, 70]). All the output maps from the models were visualized using threshold values based on maximum training sensitivity plus specificity (MTSS).
In order to investigate the role of topographic complexity (TC) as a potential factor shaping patterns of genetic differentiation, we calculated the surface ratio index for each cell from a digital elevation model using ‘DEM SURFACE TOOLS’  in ARCGIS 10.0 (ESRI, Redlands, CA, USA). Surface ratio is an index of topographic complexity, with values close to one indicating flat areas and values higher than one indicating an abrupt relief and deep slopes . We made calculations on a 90 m resolution digital elevation model from NASA Shuttle Radar Topographic Mission (SRTM Digital Elevation Data, ) and the final layer was transformed to 30 arc-sec (c. 1 km) resolution for subsequent analyses. Additionally, we used the digital elevation model to calculate a matrix of differences in elevation between each pair of studied populations (i.e. an elevation dissimilarity matrix).
Environmental characterization of populations
In order to analyze the potential role of environment as a driving factor of genetic differentiation, we characterized the environmental space of each population using a principal component analysis (PCA) with ‘varimax’ rotation applied to the values of the 19 present-day bioclimatic variables from WorldClim extracted from sampling sites, occurrence points used in MAXENT and 1 000 randomly distributed points in the study area. This procedure allowed us to capture the environmental variation of the study area and avoid potential bias resulting from just considering environmental conditions from the sampling sites. Then, we obtained for each population the PC scores of the first three PCs, which explained the 73.18, 10.92 and 8.38 % respectively of the environmental variance (Additional file 1: Figure S1). Finally, we calculated environmental dissimilarity between each pair of populations using Euclidean distances for the obtained PC scores using the ‘dist’ function in R. We performed PCA analysis using IBM SPSS 21.0 (IBM Coorp., Armonk, NY, USA).
Isolation by resistance matrices
We applied circuit theory to model gene flow between populations and test the effects of different landscape resistance scenarios (IBR) on observed patterns of genetic differentiation [12, 74]. We used CIRCUITSCAPE 4.0  to calculate resistance distance matrices between all pairs of populations considering an eight-neighbor cell connection scheme. We obtained different IBR distance matrices considering as inputs in CIRCUITSCAPE the following raster layers: current, LGM and LIG climate niche suitability and topographic complexity. We assigned pixel values of climate niche suitability maps as conductance values and pixel values of topographic complexity layers as resistance values. We also used CIRCUITSCAPE to test for the effect of isolation-by-distance (IBD) by calculating pairwise resistance distances on a completely ‘flat’ landscape based on a raster layer in which all cells had an equal value (conductance = 1). This IBD resistance model yields similar results than a matrix of Euclidean geographical distances, but it is more appropriate for comparison with others competing models also generated with CIRCUITSCAPE [ 75 , 76 ].
Multiple matrix regression with randomization
All IBR matrices were tested against genetic distance matrices (pairwise F ST and F STNA-values) using multiple matrix regressions with randomization (MMRR, ) as implemented in R. In these models, we also included elevation (ELEVDIS, see ‘Topographic complexity’ section) and climate dissimilarity (CLIMDIS, see ‘Environmental characterization of populations’ section) distance matrices in order to test a possible pattern of IBE. We used a backward procedure to select final models, removing non-significant variables from an initial full model including all explanatory predictors. We tested the significance of the remaining variables again until no additional term reached significance (e.g., ).
Distance-based redundancy analysis
Complementary to MMRR analyses, we tested the relationship between genetic differentiation, geography and environment using distance-based redundancy analyses (dbRDA, ). This approach is based on a multivariate multiple regression and estimates the percentage of genetic variation explained by a given predictor or set of predictors. We performed dbRDA using the ‘capscale’ function in the package VEGAN  as implemented in R. The genetic distance matrix (pairwise F ST or F STNA-values) was tested against the following variables: i) geographic distances (IBD), ii) elevation and iii) population’s PC scores of the first three axes from the PCA performed on the environmental data (see ‘Environmental characterization of populations’ section). Euclidean geographical distances between sampled populations were calculated using GEOGRAPHIC DISTANCE MATRIX GENERATOR 1.2.3 . Geographic distances were tested after transforming the Euclidean geographical distance matrix to a continuous rectangular vector by principal coordinates analyses (PCoA) using the ‘pcnm’ function in the package VEGAN. Significance of the predictors was assessed using multivariate F-statistics with 9999 permutations using the ‘anova.cca’ function included in the package VEGAN. We first analyzed the relationship between the genetic distance matrix and each variable separately (marginal test) and then we performed a partial dbRDA (conditional test) for each variable while controlling for the influence of geography (fitted as covariate).
Geostatistical simulations and Bayesian inference
We quantified the relative effects of geography and environment on genetic differentiation using SUNDER , a novel geostatistical method modelling covariance in allele frequencies between populations as a decreasing function of geographical and ecological distances (, see also ). SUNDER uses a Bayesian framework with a MCMC algorithm to estimate the magnitude of the effects of these variables and implements a model selection procedure by cross-validation to assess which sub-model (e.g. with or without the effect of environment) best fits to the data. Using the multinomial model, we ran SUNDER with 10 million of iterations for each data set, sampling every 1 000 iterations. We set to update in the MCMC iterations all parameters of the vector θ (α, variance of allele frequencies; β G and β E , magnitude of the effect of geography and environment respectively on genetic covariance; γ, smoothness of spatial variation of allele frequencies; δ, variation in the allele frequency of a population departing from the other populations). We also set their initial state and upper bounds of their Dirichlet prior distributions following suggestions in . We visually checked trace plots for parameters to assure convergence. We used the 10 % of our data set (sites × locus) as validation set during the cross-validation procedure. We performed SUNDER analyses using as environmental matrices the elevation dissimilarity matrix (ELEVDIS) and the climate dissimilarity distance matrix (CLIMDIS), which we separately tested against IBR distance matrix based on a completely ‘flat’ landscape (IBD). Before analyses, we standardized all distance matrices by their respective standard deviations.
Climate niche modelling
We constructed past and present-day final climate niche models using six bioclimatic variables: annual mean temperature (Bio1), mean temperature of the coldest quarter (Bio11), annual precipitation (Bio12), precipitation of the driest month (Bio14), precipitation seasonality (Bio15), and precipitation of the warmest quarter (Bio18). This model had a very high value of area under the receiving operator characteristics curve (AUC = 0.919 ± 0.067 SD), indicating overall good performance. The predicted habitat suitability area in the present was consistent with the current distribution of the species, but some areas in the northern and west side of the Pyrenees where the Morales grasshopper has not been recorded were also predicted to be suitable for the species (Fig. 3a). Projections of the present-day climate niche envelope to the LGM and LIG suggesting that the Pyrenean Morales grasshopper has experienced important distributional shifts during the Pleistocene. During the LGM, areas above 1 800–2 000 m.a.s.l. resulted unsuitable for the species and its potential distribution range expanded to areas of lower altitude across the western and northern side of the Pyrenees (Fig. 3b). Conversely, the potential distribution range of the species during the LIG expanded to higher altitudes but its peripheral geographical limits were similar than in the present (Fig. 3c).
Population genetic structure
Pairwise F ST-values ranged from 0.021 to 0.216 and 53 of 55 comparisons were significant after sequential Bonferroni correction (Additional file 1: Table S1). Pairwise F STNA-values were lower than F ST-values and ranged from 0.014 to 0.148. Pairwise F ST-values were highly correlated with pairwise F STNA-values (Mantel r = 0.941; P < 0.001).
STRUCTURE analyses on all populations showed that the best-supported number of clusters was K = 2 according to the ΔK method (Fig. 2b; Additional file 1: Figure S2). These initial analyses detected a strong correspondence between the inferred genetic clusters and their geographic location, even when a broader range of K-values (K = 2–7) was evaluated. Subsequent hierarchical analyses on different subsets of populations detected further genetic structuring (Fig. 2b). Individual assignment probabilities to a certain genetic cluster were generally high and the distribution of the hierarchical genetic structure exhibited congruence with the geographical location of the studied populations. Consecutive hierarchical analyses revealed that each population constituted a single cluster, although many pairs of populations showed a considerable degree of genetic admixture (Fig. 2b).
DAPC analyses identified also a deep spatial genetic structure in concordance with the geographic location of populations. The minimum BIC values were obtained for K = 3–5 (Fig. 4a). Considering K = 4 (the K-value with the lowest BIC value), DAPC partitioned all individuals in western (TOR-NER-SAR populations), central-west (CHI-ASP), central-east (BOI-PER) and eastern (CAR-ERR-CRE and RAS) groups (Fig. 5). Discriminant functions based on DAPC analyses correctly assigned most individuals to the genetic cluster where they were assigned a priori by K-means analyses used to infer the best-supported clustering solution  (Fig. 4c). The low overlapping of the genetic clusters on the ordination plot indicated high degree of differentiation between them (Fig. 4b). When K = 3 and K = 5 were considered and compared with K = 4, DAPC revealed a hierarchical distribution of genetic variation similar to that identified by STRUCTURE (Fig. 2b and Fig. 5).
The result of the NJ tree based on D c genetic distances showed that populations were included into monophyletic groups geographically clustered according to the hierarchical structure inferred by STRUCTURE and DAPC analyses (Fig. 2a, b and Fig. 5).
Multiple matrix regression with randomization
Genetic differentiation (F ST) was positively associated with geographic distance (i.e. IBD, resistance distances based on a completely ‘flat’ landscape), topographic complexity (IBRTC) and current (IBRCURRENT), LGM (IBRLGM) and LIG (IBRLIG) habitat resistance distances when these variables were included alone into different models (all Ps < 0.012). Indeed, Mantel tests showed that all these variables were highly inter-correlated (Additional file 1: Table S2). Climatic dissimilarity (CLIMDIS) was also correlated with all other variables, but with comparatively lower Mantel r values, whereas elevation dissimilarity (ELEVDIS) was only significantly correlated with CLIMDIS (Additional file 1: Table S2). Univariate models including IBD and IBRTC provided the highest and most similar coefficients of determination (r 2) (Table 2). However, only IBRTC was included into the final model after the backward selection procedure (Fig. 6). The rest of variables did not remain significant when they were tested against topographic complexity (all Ps > 0.13). Analyses based on F STNA yielded similar results, but models had generally lower values of coefficient of determination (r 2) (Table 2). Our results remained similar after sequential Bonferroni correction for multiple testing .
Distance-based redundancy analysis
Marginal tests showed a significant association between genetic differentiation (F ST and F STNA-values) and geography and one environmental predictor (climate PC2) that explained 46.80–53.66 % and 40.12–45.70 % of genetic variation, respectively (Table 3). Climate PC2 was explained by a pool of bioclimatic variables related to annual temperature variation (Bio2, Bio3 and Bio7; Additional file 1: Table S3). However, the environmental predictor (climate PC2) remained not significant after accounting for the influence of geography in the conditional test (Table 3). Our results remained similar after sequential Bonferroni correction . Analyses based on F STNA-values (Table 3) or performed with PCA considering only the six bioclimatic layers employed to build the climate niche model in MAXENT provided similar results (data not shown).
Geostatistical simulations and Bayesian inference
SUNDER analyses indicated that the models that best fit to the genetic data were those exclusively including the geographic component (Table 4). The Bayesian posterior estimates of the parameter β G (representing the magnitude of the effect of geography) were smaller than β E (representing the magnitude of the effect of environment), indicating a primordial effect of geographical distances and an absence of isolation by climate or elevation on genetic differentiation (Table 4).
Despite its small distribution range (<200 km between the most distant populations), the Pyrenean Morales grasshopper exhibits a strong genetic structure congruent with the geographical location of its different populations. Our microsatellite-based clustering and phylogenetic analyses revealed a strong hierarchical structure and a relatively low degree of inter-group genetic admixture. In most cases, the distinct genetic clusters corresponded to single populations, which clustered in the hierarchically superior genetic group according to the main mountain chains of the region. In comparison with other Mediterranean orthoptera, the Pyrenean Morales grasshopper has a population genetic structure as remarkable as that found at a wider spatial scale in Mioscirtus wagneri (global F ST ~ 0.19) and much higher than that shown for Ramburiella hispanica at a similar spatial scale (global F ST ~ 0.02), two species presenting a patchy distribution restricted to highly isolated and fragmented habitats [14, 81, 82]. The genetic structure found in the Pyrenean Morales grasshopper is in accordance with isolation driven by geographical factors, in which the presence of deep barriers disrupting gene flow between populations primarily explained levels of genetic differentiation . Despite climate warming during LIG and present has probably led to upward altitudinal shifts in the species distribution, its populations have presented a continuous distribution and exhibited a similar connectivity during both cold and warm periods characterizing the last 120 kya. Our niche models suggest that only a few contemporary populations probably went extinct during the LGM, but this might have had a little impact on global patterns of genetic structure if nearby populations, with a similar genetic makeup, colonized present-day suitable habitat patches. The absence of the species in areas identified as climatically suitable and stable during the last 120 kya according to our niche models (i.e. large areas in the peripheral northern and western portion of the Pyrenees) could be linked to historical events (such as extinctions) predating the temporal scale addressed in our study . It could also depend on constraints of our climate models not capturing other important abiotic or biotic interactions that may contribute to the distribution of the species in those areas . Thus, our analyses suggest that the strong genetic structure found across the species distribution range has not arisen as consequence of long-term isolation driven by Pleistocene climatic oscillations that shaped range-wide patterns of genetic structure in many other taxa from temperate regions [34, 85, 86].
Our different landscape genetic approaches confirmed that neutral divergence resulted from the isolating effects of topography primarily drove the deep patterns of population genetic differentiation observed in the Morales grasshopper. The resistance model based on topographic complexity was the best fit to our data, indicating that physical features defining the abrupt landscape characterizing the Pyrenees shaped genetic differentiation. In particular, deep valleys with a north-south orientation and slopes that are generally steep and create canyons and ridges on the landscape crisscross the central and eastern portion of the Pyrenees inhabited by Morales grasshopper. Hence, these topographic features could become impassable barriers and restrict gene flow as has been previously documented for other species with limited dispersal capacities and inhabiting regions of remarkable topographical roughness [35, 36, 87, 88]. The remaining analyzed landscape factors, such as resistance-based distances informed by current and past climate niche models, did not show either a significant association with genetic differentiation after controlling for the influence of topographic complexity or geographical distances. The lower explanatory power of CNM-based resistance distances may be related with the fact that they are not capturing well microhabitat structure, which has been found to be highly relevant in determining the distribution and demography of grasshoppers [89–92]. The short generation time of the studied species (=1 year) may have also resulted in contemporary patterns of genetic differentiation are not capturing the genetic signal left by dispersal routes defined by habitat corridors during the past. This fact contrasts with patterns found in species with longer generation times and that are likely to show a time lag in their response to changing climatic conditions [77, 93].
We found no support for ecological divergence and local adaptation processes have contributed to population genetic differentiation and the three different employed approaches (MMRR analysis using dissimilarity matrices, dbRDA analysis using raw variables and Bayesian inference) confirmed the consistency of this result. We did not find support either for altitude as an isolating mechanism despite elevation gradients have been previously found to be a significant driver of genetic and phenotypic variation in grasshoppers and many other organisms [94, 95]. Our results contrasts with other studies that have documented an important role of environment on genetic differentiation in many taxa [4, 18], including some insects such as grasshoppers [96, 97], walking sticks  or beetles . This discrepancy may be in part due to these studies considered species exhibiting narrow feeding preferences and analyzed ecological dissimilarity in terms of host-plant associations, an aspect that may have a higher impact on genetic divergence than the climate or elevation gradients considered in our study [96–99]. The meta-analysis by  showed that isolation-by-ecology is more frequent than IBD in insects, particularly in species with strong patterns of genetic structure. Considering the high degree of genetic differentiation among our study populations, we can discard the hypothesis that a high level of gene flow has counteracted the potential disruption of gene flow driven by local adaptation processes mediated by environmental heterogeneity . Hence, the lack of effects of environment on gene flow could be due to different reasons, including low environmental variation among sampling sites  or local adaptation driven by other ecological factors not considered in this study (e.g. distinct host plants or habitat structure ).
This study emphasizes the importance of examining jointly different scenarios of population isolation to understand their contribution to the spatial distribution of genetic variation across a species range. Our analyses evidence the importance of topographic complexity in determining patterns of genetic differentiation, indicating that limited dispersal and drift, due to scarce population connectivity, is shaping the genetic structure found in our study system (e.g., ). Further research harnessing high-throughput sequencing will provide a better understanding about the potential association between loci under selection and different ecological factors, which may help to identify genomic regions involved in local adaptation processes [15, 100]. Exploring the relationship between environmental features and genetic and phenotypic patterns of variation could also provide insights about the potential interplay of evolutionary and ecological processes in shaping range-wide patterns of genetic differentiation [101, 102].
This study did not require ethical approval. Specimens were collected under license from ‘Gobierno de Aragón’, ‘Generalitat de Catalunya’ and ‘Ordesa y Monte Perdido National Park’. Our sampling procedures did not affect the survival of the studied populations.
Availability of data and materials
Nuclear microsatellite data are available in the LabArchives repository (http://dx.doi.org/10.6070/H4QC01JB). All other data supporting the results of this article are included within the article and its additional files.
Consent to publish
area under the curve.
Bayesian information criterion
climate niche models
discriminant analysis of principal components
distance-based redundancy analyses
thousand years ago
last glacial maximum
meters above sea level
Markov chain Monte Carlo
multiple matrix regressions with randomization
maximum training sensitivity plus specificity threshold
principal component analysis
principal coordinates analyses
root mean squared error
Carnaval AC, Hickerson MJ, Haddad CFB, Rodrigues MT, Moritz C. Stability predicts genetic diversity in the Brazilian Atlantic forest hotspot. Science. 2009;323:785–9.
Yannic G, Pellissier L, Ortego J, Lecomte N, Couturier S, Cuyler C, et al. Genetic diversity in caribou linked to past and future climate change. Nat Clim Change. 2014;4:132–7.
Lee CR, Mitchell-Olds T. Quantifying effects of environmental and geographical factors on patterns of genetic differentiation. Mol Ecol. 2011;20:4631–42.
Shafer ABA, Wolf JBW. Widespread evidence for incipient ecological speciation: a meta-analysis of isolation-by-ecology. Ecol Lett. 2013;16:940–50.
Wang IJ, Glor RE, Losos JB. Quantifying the roles of ecology and geography in spatial genetic divergence. Ecol Lett. 2013;16:175–82.
Wright S. Isolation by distance. Genetics. 1943;28:114–38.
Slatkin M. Isolation by distance in equilibrium and nonequilibrium populations. Evolution. 1993;47:264–79.
Jenkins DG, Carey M, Czerniewska J, Fletcher J, Hether T, Jones A, et al. A meta-analysis of isolation by distance: relic or reference standard for landscape genetics? Ecography. 2010;33:315–20.
Manel S, Holderegger R. Ten years of landscape genetics. Trends Ecol Evol. 2013;28:614–21.
Manel S, Schwartz MK, Luikart G, Taberlet P. Landscape genetics: combining landscape ecology and population genetics. Trends Ecol Evol. 2003;18:189–97.
McRae BH. Isolation by resistance. Evolution. 2006;60:1551–61.
McRae BH, Beier P. Circuit theory predicts gene flow in plant and animal populations. Proc Natl Acad Sci U S A. 2007;104:19885–90.
Ruiz-Gonzalez A, Cushman SA, Madeira MJ, Randi E, Gómez-Moliner BJ. Isolation by distance, resistance and/or clusters? Lessons learned from a forest-dwelling carnivore inhabiting a heterogeneous landscape. Mol Ecol. 2015;24:5110–29.
Ortego J, Aguirre MP, Noguerales V, Cordero PJ. Consequences of extensive habitat fragmentation in landscape-level patterns of genetic diversity and structure in the Mediterranean esparto grasshopper. Evol Appl. 2015;8:621–32.
Wang IJ, Bradburd GS. Isolation by environment. Mol Ecol. 2014;23:5649–62.
Nosil P, Vines TH, Funk DJ. Reproductive isolation caused by natural selection against immigrants from divergent habitats. Evolution. 2005;59:705–19.
Nosil P. Ecological speciation. New York: Oxford University Press; 2012.
Sexton JP, Hangartner SB, Hoffmann AA. Genetic isolation by environment or distance: which pattern of gene flow is most common? Evolution. 2014;68:1–15.
Funk DJ, Nosil P, Etges WJ. Ecological divergence exhibits consistently positive associations with reproductive isolation across disparate taxa. Proc Natl Acad Sci U S A. 2006;103:3209–13.
Crispo E, Bentzen P, Reznick DN, Kinnison MT, Hendry AP. The relative influence of natural selection and geography on gene flow in guppies. Mol Ecol. 2006;15:49–62.
Edwards DL, Keogh JS, Knowles LL. Effects of vicariant barriers, habitat stability, population isolation and environmental features on species divergence in the south-western Australian coastal reptile community. Mol Ecol. 2012;21:3809–22.
Meirmans PG. The trouble with isolation by distance. Mol Ecol. 2012;21:2839–46.
Balkenhol N, Waits LP, Dezzani RJ. Statistical approaches in landscape genetics: an evaluation of methods for linking landscape and genetic data. Ecography. 2009;32:818–30.
Wang IJ. Examining the full effects of landscape heterogeneity on spatial genetic variation: a multiple matrix regression approach for quantifying geographic and ecological isolation. Evolution. 2013;67:3403–11.
Botta F, Eriksen C, Fontaine MC, Guillot G. Enhanced computational methods for quantifying the effect of geographic and environmental isolation on genetic differentiation. Methods Ecol Evol. 2015;6:1270–7.
Kierepka EM, Latch EK. Performance of partial statistics in individual-based landscape genetics. Mol Ecol Resour. 2015;15:512–25.
Forester BR, Jones MR, Joost S, Landguth EL, Lasky JR. Detecting spatial genetic signatures of local adaptation in heterogeneous landscapes. Mol Ecol. 2016;25:104–20.
Ferrer ES, García-Navas V, Bueno-Enciso J, Barrientos R, Serrano-Davies E, Cáliz-Campal C, et al. The influence of landscape configuration and environment on population genetic structure in a sedentary passerine: insights from loci located in different genomic regions. J Evol Biol. 2016;29:205–19.
Munshi-South J, Zolnik CP, Harris S. Population genomics of the Anthropocene: urbanization is negatively associated with genome-wide variation in white-footed mouse populations. Evol Appl. 2016;9:546–64.
Defaut B. Preliminary revision of Chorthippus of the binotatus group (Charpentier, 1825) (Caelifera, Acrididae, Gomphocerinae). Materiaux Orthopteriques et Entomocenotiques. 2011;16:17–54.
Llucia-Pomares D. Revision of the Orthoptera (Insecta) of Catalonia (Spain). Monografias SEA, vol. 7. Zaragoza: Sociedad Entomológica Aragonesa; 2002.
García-Barros E, Gurrea P, Luciañez MJ, Cano JM, Munguira ML, Moreno JC, et al. Parsimony analysis of endemicity and its application to animal and plant geographical distributions in the Ibero-Balearic region (western Mediterranean). J Biogeogr. 2002;29:109–24.
Jiménez-Sánchez M, Rodríguez-Rodríguez L, García-Ruiz JM, Domínguez-Cuesta MJ, Farias P, Valero-Garcés B, et al. A review of glacial geomorphology and chronology in northern Spain: timing and regional variability during the last glacial cycle. Geomorphology. 2013;196:50–64.
Hewitt G. The genetic legacy of the Quaternary ice ages. Nature. 2000;405:907–13.
Wang IJ. Environmental and topographic variables shape genetic structure and effective population sizes in the endangered Yosemite toad. Divers Distrib. 2012;18:1033–41.
Castillo JA, Epps CW, Davis AR, Cushman SA. Landscape effects on gene flow for a climate-sensitive montane species, the American pika. Mol Ecol. 2014;23:843–56.
Sobel JM, Chen GF, Watt LR, Schemske DW. The biology of speciation. Evolution. 2010;64:295–315.
Legendre P, Anderson MJ. Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments. Ecol Monogr. 1999;69:1–24.
Pflüeger FJ, Balkenhol N. A plea for simultaneously considering matrix quality and local environmental conditions when analysing landscape impacts on effective dispersal. Mol Ecol. 2014;23:2146–56.
Aljanabi SM, Martínez I. Universal and rapid salt-extraction of high quality genomic DNA for PCR-based techniques. Nucleic Acids Res. 1997;25:4692–3.
Basiita RK, Bruggemann JH, Cai N, Cáliz-Campal C, Chen C, Chen J, et al. Microsatellite records for volume 7, issue 4. Conserv Genet Resour. 2015;7:917–44.
Guo SW, Thompson EA. A Monte-Carlo method for combined segregation and linkage analysis. Am J Hum Genet. 1992;51:1111–26.
Excoffier L, Lischer HEL. ARLEQUIN suite version 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10:564–7.
Van Oosterhout C, Hutchinson WF, Wills DPM, Shipley P. MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol Ecol Notes. 2004;4:535–8.
Rice WR. Analyzing tables of statistical tests. Evolution. 1989;43:223–5.
Chapuis MP, Lecoq M, Michalakis Y, Loiseau A, Sword GA, Piry S, et al. Do outbreaks affect genetic population structure? A worldwide survey in Locusta migratoria, a pest plagued by microsatellite null alleles. Mol Ecol. 2008;17:3640–53.
Chapuis MP, Estoup A. Microsatellite null alleles and estimation of population differentiation. Mol Biol Evol. 2007;24:621–31.
Falush D, Stephens M, Pritchard JK. Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics. 2003;164:1567–87.
Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.
Papadopoulou A, Knowles LL. Genomic tests of the species-pump hypothesis: recent island connectivity cycles drive population divergence but not speciation in Caribbean crickets across the Virgin Islands. Evolution. 2015;69:1501–17.
Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14:2611–20.
Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23:1801–6.
Rosenberg NA. DISTRUCT: a program for the graphical display of population structure. Mol Ecol Notes. 2004;4:137–8.
Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 2010;11:94.
R Core Team. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2014. http://www.R-project.org/. Accessed 29 Dec 2014.
Langella O. POPULATIONS 1.2.31 software. 1999. http://bioinformatics.org/populations/. Accessed 25 Feb 2015.
Cavalli-Sforza L, Edwards AWF. Phylogenetic analyses: models and estimation procedures. Evolution. 1967;21:550–70.
Takezaki N, Nei M. Genetic distances and reconstruction of phylogenetic trees from microsatellite DNA. Genetics. 1996;144:389–99.
Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Model. 2006;190:231–59.
Phillips SJ, Dudik M. Modeling of species distributions with MAXENT: new extensions and a comprehensive evaluation. Ecography. 2008;31:161–75.
WorldClim: Global Climate Data. http://www.worldclim.org/. Accessed 17 Oct 2015.
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.
Collins WD, Bitz CM, Blackmon ML, Bonan GB, Bretherton CS, Carton JA, et al. The Community Climate System Model version 3 (CCSM3). J Clim. 2006;19:2122–43.
Braconnot P, Otto-Bliesner B, Harrison S, Joussaume S, Peterchmitt JY, Abe-Ouchi A, et al. Results of PMIP2 coupled simulations of the Mid-Holocene and Last Glacial Maximum - part 1: experiments and large-scale features. Clim Past. 2007;3:261–77.
Otto-Bliesner BL, Marsha SJ, Overpeck JT, Miller GH, Hu AX. CAPE Last Interglacial Project Members. Simulating arctic climate warmth and icefield retreat in the last interglaciation. Science. 2006;311:1751–3.
Anderson RP, Raza A. The effect of the extent of the study region on GIS models of species geographic distributions and estimates of niche evolution: preliminary tests with montane rodents (genus Nephelomys) in Venezuela. J Biogeogr. 2010;37:1378–93.
Alvarado-Serrano DF, Knowles LL. Ecological niche models in phylogeographic studies: applications, advances and precautions. Mol Ecol Resour. 2014;14:233–48.
Elith J, Kearney M, Phillips S. The art of modelling range-shifting species. Methods Ecol Evol. 2010;1:330–42.
Lanier HC, Massatti R, He Q, Olson LE, Knowles LL. Colonization from divergent ancestors: glaciation signatures on contemporary patterns of genomic variation in Collared Pikas (Ochotona collaris). Mol Ecol. 2015;24:3688–705.
Massatti R, Knowles LL. Microhabitat differences impact phylogeographic concordance of codistributed species: genomic evidence in montane sedges (Carex L.) from the rocky mountains. Evolution. 2014;68:2833–46.
Jenness JS. DEM SURFACE TOOLS. Jenness Enterprises. 2013. http://www.jennessent.com/arcgis/surface_area.htm. Accessed 15 Jun 2015
Jenness JS. Calculating landscape surface area from digital elevation models. Wildlife Soc B. 2004;32:829–39.
NASA Shuttle Radar Topographic Mission: SRTM Digital Elevation Data. 2015. http://srtm.csi.cgiar.org/. Accessed 15 Jun 2015
McRae BH, Dickson BG, Keitt TH, Shah VB. Using circuit theory to model connectivity in ecology, evolution, and conservation. Ecology. 2008;89:2712–24.
Jha S, Kremen C. Urban land use limits regional bumble bee gene flow. Mol Ecol. 2013;22:2483–95.
Velo-Antón G, Parra JL, Parra-Olea G, Zamudio KR. Tracking climate change in a dispersal-limited species: reduced spatial and genetic connectivity in a montane salamander. Mol Ecol. 2013;22:3261–78.
Ortego J, Gugger PF, Sork VL. Climatically stable landscapes predict patterns of genetic structure and admixture in the Californian canyon live oak. J Biogeogr. 2015;42:328–38.
Oksanen J, Blanchet FG, Kindt R, Legendre P, Minchin PR, O’Hara RB et al. VEGAN: community ecology package. R Package Version 2.3–1. 2015. http://r-forge.r-project.org/projects/vegan. Accessed 22 Nov 2015.
Erst PJ. GEOGRAPHIC DISTANCE MATRIX GENERATOR, version 1.2.3. American Museum of Natural History, Center for Biodiversity and Conservation. http://biodiversityinformatics.amnh.org/open_source/gdmg. Accessed 15 Nov 2015.
Bradburd GS, Ralph PL, Coop GM. Disentangling the effects of geographic and ecological isolation on genetic differentiation. Evolution. 2013;67:3258–73.
Ortego J, Aguirre MP, Cordero PJ. Genetic and morphological divergence at different spatiotemporal scales in the grasshopper Mioscirtus wagneri (Orthoptera: Acrididae). J Insect Conserv. 2012;16:103–10.
Ortego J, Aguirre MP, Cordero PJ. Population genetics of Mioscirtus wagneri, a grasshopper showing a highly fragmented distribution. Mol Ecol. 2010;19:472–83.
Rizzo V, Comas J, Fadrique F, Fresneda J, Ribera I. Early Pliocene range expansion of a clade of subterranean Pyrenean beetles. J Biogeogr. 2013;40:1861–73.
Hampe A. Bioclimate envelope models: what they detect and what they hide. Global Ecol Biogeogr. 2004;13:469–71.
Hewitt GM. Genetic consequences of climatic oscillations in the Quaternary. Philos Trans R Soc B-Biol Sci. 2004;359:183–95.
Knowles LL, Richards CL. Importance of genetic drift during Pleistocene divergence as revealed by analyses of genomic variation. Mol Ecol. 2005;14:4023–32.
Murphy MA, Dezzani R, Pilliod DS, Storfer A. Landscape genetics of high mountain frog metapopulations. Mol Ecol. 2010;19:3634–49.
Benham PM, Witt CC. The dual role of Andean topography in primary divergence: functional and neutral variation among populations of the hummingbird, Metallura tyrianthina. BMC Evol Biol. 2016;16:22.
Reinhardt K, Kohler G, Maas S, Detzel P. Low dispersal ability and habitat specificity promote extinctions in rare but not in widespread species: the Orthoptera of Germany. Ecography. 2005;28:593–602.
San M, Gómez G, Van Dyck H. Ecotypic differentiation between urban and rural populations of the grasshopper Chorthippus brunneus relative to climate and habitat fragmentation. Oecologia. 2012;169:125–33.
Keller D, Holderegger R, van Strien MJ. Spatial scale affects landscape genetic analysis of a wetland grasshopper. Mol Ecol. 2013;22:2467–82.
Gauffre B, Mallez S, Chapuis MP, Leblois R, Litrico I, Delaunay S, et al. Spatial heterogeneity in landscape structure influences dispersal and genetic structure: empirical evidence from a grasshopper in an agricultural landscape. Mol Ecol. 2015;24:1713–28.
He Q, Edwards DL, Knowles LL. Integrative testing of how environments from the past to the present shape genetic structure across landscapes. Evolution. 2013;67:3386–402.
Laiolo P, Illera JC, Obeso JR. Local climate determines intra- and interspecific variation in sexual size dimorphism in mountain grasshopper communities. J Evol Biol. 2013;26:2171–83.
Roff DA, Mousseau T. The evolution of the phenotypic covariance matrix: evidence for selection and drift in Melanoplus. J Evol Biol. 2005;18:1104–14.
Grace T, Wisely SM, Brown SJ, Dowell FE, Joern A. Divergent host plant adaptation drives the evolution of sexual isolation in the grasshopper Hesperotettix viridis (Orthoptera: Acrididae) in the absence of reinforcement. Biol J Linn Soc. 2010;100:866–78.
Hernández-Teixidor D, López H, Nogales M, Emerson BC, Juan C, Oromí P. Genetic, morphological, and dietary changes associated with novel habitat colonisation in the Canary Island endemic grasshopper Acrostira bellamyi. Ecol Entomol. 2014;39:703–15.
Nosil P, Egan SP, Funk DJ. Heterogeneous genomic differentiation between walking-stick ecotypes: “Isolation by adaptation” and multiple roles for divergent selection. Evolution. 2008;62:316–36.
Funk DJ, Egan SP, Nosil P. Isolation by adaptation in Neochlamisus leaf beetles: host-related selection promotes neutral genomic divergence. Mol Ecol. 2011;20:4671–82.
Manthey JD, Moyle RG. Isolation by environment in White-breasted Nuthatches (Sitta carolinensis) of the Madrean archipelago sky islands: a landscape genomics approach. Mol Ecol. 2015;24:3628–38.
Sistrom M, Edwards DL, Donnellan S, Hutchinson M. Morphological differentiation correlates with ecological but not with genetic divergence in a Gehyra gecko. J Evol Biol. 2012;25:647–60.
Wang IJ, Summers K. Genetic structure is correlated with phenotypic divergence rather than geographic isolation in the highly polymorphic strawberry poison-dart frog. Mol Ecol. 2010;19:447–58.
Parks DH, Porter M, Churcher S, Wang S, Blouin C, Whalley J, et al. GENGIS: a geospatial information system for genomic data. Genome Res. 2009;19:1896–904.
We wish to thank Bernard Defaut for providing valuable information about sampling locations. Two anonymous referees provided useful and valuable comments on an earlier draft of this manuscript. We acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).
VN was supported by a FPI pre-doctoral scholarship (BES-2012-053741) from Ministerio de Economía y Competitividad. JO was supported by Severo Ochoa (SEV-2012-0262) and Ramón y Cajal (RYC-2013-12501) research fellowships. This work received financial support from research grants CGL2011-25053 (Ministerio de Ciencia e Innovación and European Social Fund), POII10-0197-0167, PEII-2014-023-P (Junta de Comunidades de Castilla-La Mancha and European Social Fund) and UNCM08-1E-018 (European Regional Development Fund).
The authors declare that they have no competing interests.
VN and JO conceived and designed the study. VN carried out the lab work and analyzed the data guided by JO. VN, PJC and JO collected the samples. VN wrote the manuscript with help of JO and inputs from PJC. All authors read and approved the final manuscript.
Genetic differentiation between eleven populations of the Pyrenean Morales grasshopper (Chorthippus saulcyi moralesi). Table S2. Results of Mantel tests analyzing the relationship between the different distance matrices (predictors) used to evaluate the factors associated with population genetic differentiation in the Pyrenean Morales grasshopper (Chorthippus saulcyi moralesi). Table S3. Results of principal component analysis (PCA) applied to the values of the 19 present day bioclimatic variables obtained from the WorldClim dataset. Figure S1. The position in the environmental space (first two principal components of a PCA based on the 19 bioclimatic variables from the WorldClim dataset) of all known populations of the Pyrenean Morales grasshopper (Chorthippus saulcyi moralesi). Figure S2. Results of Bayesian clustering analyses in STRUCTURE to determine the best-supported number of clusters in hierarchical analyses. (DOCX 871 kb)
About this article
- Chorthippus saulcyi moralesi
- Ecological niche modelling
- Hierarchical genetic structure
- Isolation by environment
- Isolation by resistance
- Landscape genetics
- Topographic complexity