BACKGROUND
understanding the factors structuring genetic variation in natural populations is a paramount topic in evolutionary biology  <cit> .
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 .
the isolation-by-distance  model predicts that genetic differentiation increases with euclidean geographic distance because of limited dispersal and genetic drift  <cit> .
even though this classic model explains spatial patterns of genetic differentiation in a wide range of organisms , it does not consider more sophisticated information than straight-line geographic distances and assumes landscape homogeneity, an unrealistic scenario for most natural systems  <cit> .
recently, the emergence of landscape genetics has explicitly incorporated landscape complexity into the study of evolutionary processes  <cit> , 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 .
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  <cit> .
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 .
recent research has revealed the ubiquity of ibe patterns, highlighting the importance of ecological factors in shaping genetic structure of natural populations .
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  <cit> .
indeed, different isolating mechanisms are not mutually exclusive and gene flow is often influenced by a combination of geographical and ecological factors  <cit> .
given that geography and environment are often highly correlated, disentangling their relative impact on population genetic differentiation harbors inherent analytical difficulties  that have promoted the progressive development of more robust and accurate statistical methods .
however, only a few studies have jointly considered the relative effects of geography, landscape composition and environmental heterogeneity on either landscape-level  <cit>  or range-wide patterns of genetic structure  <cit> .
the pyrenean morales grasshopper   is a narrow-endemic subspecies belonging to the chorthippus binotatus group, exclusively distributed in central aragón and catalonia pyrenean mountains .
it is a winged grasshopper primarily feeding on gramineous herbs  <cit>  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  <dig> and 2400 m.a.s.l.
<cit> .
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  <cit> .
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  <cit> .
these mountains experienced dramatic climate fluctuations during the pleistocene  <cit> , which are expected to have strongly influenced the demographic history and altered the distribution of many organisms of the region  <cit> .
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  <cit> .
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  <cit> .
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  in shaping patterns of genetic differentiation throughout the entire distribution range of a narrowly distributed taxon  <cit> .
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 .
in particular, we used multiple matrix regressions with randomization , distance-based redundancy analyses  and geostatistical modelling based on bayesian inference  <cit>  to test whether the spatial pattern of genetic differentiation in the pyrenean morales grasshopper is explained by i) geographic distances , ii) resistance distances based on topographic complexity and current and past  climate suitability ; and iii) altitudinal and environmental dissimilarity between populations  .
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  <cit> .
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  <cit> .fig.
<dig> workflow summarizing the methodological approach employed in this study to analyze the relative contribution of isolation by distance , isolation by resistance , and isolation by environment  in structuring genetic variation in the pyrenean morales grasshopper.
the response variables and predictors considered for each analytical approach  are indicated.
hs: habitat suitability; lgm: last glacial maximum; lig: last interglacial; mmrr: multiple matrix regression with randomization; dbrda: distance-based redundancy analysis
methods
population sampling
between  <dig> and  <dig>  we collected  <dig> individuals from  <dig> populations of the pyrenean morales grasshopper.
we aimed to sample populations throughout the entire distribution range of the species  based on occurrence-data available in the literature  <cit>  and our own prospection of areas with potentially suitable habitats.
the morales grasshopper is primarily distributed in the south side of the pyrenees  and a small area located in the northeastern part of these mountains .
the latter portion of the species distribution was intensively prospected during our field work but we only found a single population in the area .
overall, we were able to collect individuals from populations located across almost the entire climatic and elevation gradient occupied by the species , including populations from a wide range of habitats  and differing in up to ~ 900 m of elevation .
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  <dig>  specimens were collected in public lands under license from ‘gobierno de aragón’, ‘generalitat de catalunya’ and ‘ordesa y monte perdido national park’.fig.
<dig> population genetic structure in the pyrenean morales grasshopper.
panel a shows sampling sites of the species and phylogenetic relationships among the  <dig> populations inferred from a neighbor-joining  tree based on cavalli-sforza and edwards chord distances .
the tree was plotted on a topographic map of the pyrenees using the software gengis [
103
].
we downloaded topographic data from nasa shuttle radar topographic mission  as 90 m resolution digital elevation model and subsequently transformed to  <dig> arc-sec  resolution for representation.
panel b represents the results of genetic assignment of  <dig> individuals of the pyrenean morales grasshopper based on the bayesian method implemented in structure.
we performed hierarchical analyses for subsets of populations considering the most probable k-value inferred at the previous hierarchical level .
each individual corresponds to a vertical bar partitioned into k-colored segments that represent the individual’s probability of belonging to the cluster with that color.
black lines separate individuals from different populations.
population codes as in table 1
n
microsatellite genotyping and basic genetic statistics
we extracted genomic dna from a hind-leg of each individual using a salt extraction protocol  <cit> .
we amplified and genotyped each individual using the  <dig> microsatellites markers described in  <cit> .
we performed pcr amplifications following the procedure described in  <cit> , run pcr products on an abi  <dig> genetic analyzer  and scored genotypes using genemapper  <dig>  .
we tested for deviations from hardy-weinberg equilibrium  using exact tests  <cit>  based on  <dig>  <dig> markov chain iterations and  <dig>  <dig> dememorization steps as implemented in the program arlequin  <dig>   <cit> .
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  <cit> .
we also used arlequin to test for linkage disequilibrium using a likelihood-ratio statistic, whose distribution was generated with  <dig>  <dig> permutations.
we did not find evidence for linkage disequilibrium between any pair of loci in any sampling population after sequential bonferroni corrections  <cit> .
genetic structure analyses
we estimated population genetic differentiation calculating fst-values between all pairs of sampling populations and testing their significance with fisher’s exact test after  <dig>  <dig> permutations using arlequin  <dig> .
we corrected p-values using a sequential bonferroni adjustment  <cit> .
due to the high frequency of null alleles in orthoptera , we also calculated pairwise fst-values corrected for null alleles  using the so-called ena-method implemented in the program freena  <cit> .
we inferred genetic structure using bayesian clustering analyses in structure  <dig> .
<dig>  <cit> .
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 .
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  with a burn-in period of  <dig>  <dig> steps and a run length of  <dig>  <dig>  <dig> markov chain monte carlo  cycles.
we estimated the best-supported number of genetic clusters with the log probability of the data   <cit>  and the Δk method  <cit> .
we used the ‘full search’ algorithm in the program clumpp  <dig> .
<dig>  <cit>  to align replicated runs and average individual assignment probabilities for the most likely k-value.
finally, we used distruct  <dig>   <cit>  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  to identify clusters of genetically related individuals  <cit> .
although both clustering methods exhibit similarities , several important differences exist in their analytical approaches.
structure suffers from the assumption of hwe and gametic disequilibrium  <cit>  and typically fails to detect isolation-by-distance   <cit>  and hierarchical patterns of spatial genetic structure  <cit> .
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  <cit> .
dapc is a methodological approach that requires data transformation using a principal component analysis  as a prior step to a discriminant analysis .
da partitions genetic variation maximizing differences between clusters while minimizing within-cluster variation.
we implemented dapc analysis in r  <dig> .
<dig> using the package adegenet  <cit> .
at first, we used the ‘find.cluster’ function using all available principal components  to determine the best-supported number of genetic clusters using the bayesian information criterion .
the ‘find cluster’ function runs successive k-means clustering with increasing number of clusters  and provides a bic value for each simulated k-value .
secondarily, we determined the optimal number of pcs for the dapc by cross-validation using the ‘xvaldapc’ function with  <dig> replicates.
we selected the number of pcs associated with the lowest ‘root mean squared error’  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  <dig> .
<dig>  <cit>  to obtain a neighbor-joining  tree based on pairwise cavalli-sforza and edwards  genetic distances  <cit> .
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  <cit> .
climate niche modelling
we used climate niche models  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  <dig> .
<dig>  <cit>  based on current climate layers and using  <dig> cross-validation replicate model runs.
we used a total of  <dig> occurrence points obtained from the literature  <cit>  and our own sampling.
to construct the models, we used the  <dig> present-day bioclimatic variables available in worldclim  <cit>  and downloaded at  <dig> arc-sec  resolution  <cit> .
to obtain the distribution of the pyrenean morales grasshopper in the last glacial maximum  and the last interglacial , we projected contemporary species-climate relationships to these periods.
the lgm layers were based on the community climate system model  from the palaeoclimate modelling intercomparison project phase ii .
we downloaded lgm layers from worldclim at  <dig>  arc-min and interpolated to  <dig> arc-sec resolutions.
lig layers were based on  <cit>  and downloaded from worldclim at  <dig> arc-sec resolution.
according to the suggestions from  <cit> , 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  calculation to address the problems derived from projecting the current distribution into novel climates   <cit> .
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 .
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  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  and we compared its lgm and lig projections with those obtained by mess analyses .
all the output maps from the models were visualized using threshold values based on maximum training sensitivity plus specificity .
topographic complexity
in order to investigate the role of topographic complexity  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’  <cit>  in arcgis  <dig>  .
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  <cit> .
we made calculations on a 90 m resolution digital elevation model from nasa shuttle radar topographic mission  and the final layer was transformed to  <dig> arc-sec  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 .
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  with ‘varimax’ rotation applied to the values of the  <dig> present-day bioclimatic variables from worldclim extracted from sampling sites, occurrence points used in maxent and  <dig>  <dig> 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  <dig> ,  <dig>  and  <dig>  % respectively of the environmental variance .
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  <dig>  .
isolation by resistance matrices
we applied circuit theory to model gene flow between populations and test the effects of different landscape resistance scenarios  on observed patterns of genetic differentiation  <cit> .
we used circuitscape  <dig>   <cit>  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  by calculating pairwise resistance distances on a completely ‘flat’ landscape based on a raster layer in which all cells had an equal value .
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  <cit> .
multiple matrix regression with randomization
all ibr matrices were tested against genetic distance matrices  using multiple matrix regressions with randomization  as implemented in r. in these models, we also included elevation  and climate dissimilarity  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 .
distance-based redundancy analysis
complementary to mmrr analyses, we tested the relationship between genetic differentiation, geography and environment using distance-based redundancy analyses .
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  <cit>  as implemented in r. the genetic distance matrix  was tested against the following variables: i) geographic distances , ii) elevation and iii) population’s pc scores of the first three axes from the pca performed on the environmental data .
euclidean geographical distances between sampled populations were calculated using geographic distance matrix generator  <dig> .
<dig>  <cit> .
geographic distances were tested after transforming the euclidean geographical distance matrix to a continuous rectangular vector by principal coordinates analyses  using the ‘pcnm’ function in the package vegan.
significance of the predictors was assessed using multivariate f-statistics with  <dig> 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  and then we performed a partial dbrda  for each variable while controlling for the influence of geography .
geostatistical simulations and bayesian inference
we quantified the relative effects of geography and environment on genetic differentiation using sunder  <cit> , a novel geostatistical method modelling covariance in allele frequencies between populations as a decreasing function of geographical and ecological distances .
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  best fits to the data.
using the multinomial model, we ran sunder with  <dig> million of iterations for each data set, sampling every  <dig>  <dig> iterations.
we set to update in the mcmc iterations all parameters of the vector θ .
we also set their initial state and upper bounds of their dirichlet prior distributions following suggestions in  <cit> .
we visually checked trace plots for parameters to assure convergence.
we used the 10 % of our data set  as validation set during the cross-validation procedure.
we performed sunder analyses using as environmental matrices the elevation dissimilarity matrix  and the climate dissimilarity distance matrix , which we separately tested against ibr distance matrix based on a completely ‘flat’ landscape .
before analyses, we standardized all distance matrices by their respective standard deviations.
RESULTS
climate niche modelling
we constructed past and present-day final climate niche models using six bioclimatic variables: annual mean temperature , mean temperature of the coldest quarter , annual precipitation , precipitation of the driest month , precipitation seasonality , and precipitation of the warmest quarter .
this model had a very high value of area under the receiving operator characteristics curve , 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 .
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  <dig> 800– <dig> 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 .
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.
<dig> climate niche modelling of the pyrenean morales grasshopper for  the present,  the last glacial maximum  and  the last interglacial .
climatically suitable areas defined using the maximum training sensitivity plus specificity threshold  is in red.
the topography is in the background map, in which whiter colors indicate higher elevations and darker colors indicate lower elevations .
yellow dots represent the eleven sampling sites
population genetic structure
pairwise fst-values ranged from  <dig>  to  <dig>  and  <dig> of  <dig> comparisons were significant after sequential bonferroni correction .
pairwise fstna-values were lower than fst-values and ranged from  <dig>  to  <dig> .
pairwise fst-values were highly correlated with pairwise fstna-values .
structure analyses on all populations showed that the best-supported number of clusters was k =  <dig> according to the Δk method .
these initial analyses detected a strong correspondence between the inferred genetic clusters and their geographic location, even when a broader range of k-values  was evaluated.
subsequent hierarchical analyses on different subsets of populations detected further genetic structuring .
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 .
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– <dig> .
considering k =  <dig> , dapc partitioned all individuals in western , central-west , central-east  and eastern  groups .
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  <cit>  .
the low overlapping of the genetic clusters on the ordination plot indicated high degree of differentiation between them .
when k =  <dig> and k =  <dig> were considered and compared with k =  <dig>  dapc revealed a hierarchical distribution of genetic variation similar to that identified by structure .fig.
<dig> summary of the results of discriminant analysis of principal components .
panel a represents the bayesian information criterion  for each value of k. the minimum value of bic before the first increase or stabilization indicates the best-supported number of genetic clusters .
panel b represents an ordination plot for the first two discriminant axes.
each dot represents one individual and colors and inertia ellipses indicate their assignment to one of the four genetic clusters inferred by dapc.
the up-right graph inset displays the variance explained by the principal component axes used for dapc .
the bottom-right inset displays in relative magnitude the variance explained by the two discriminant axes plotted .
panel c represents whether the individuals  were correctly assigned  to the genetic cluster where they were included a priori  by k-means analyses used to infer the best-supported clustering solution.
colors represent membership probabilities to each genetic cluster  and blue crosses indicate the cluster where the individuals were originally assigned by k-means analyses.
generally, the dapc classification of individuals is consistent with their assignment to the clusters originally identified by k-means analyses .
panel d shows the number of individuals from each population  assigned to each of the four inferred genetic clusters .
the size of black squares is proportional to the number of individuals assigned to the different clusters .
population codes are described in table 1
fig.
<dig> results of clustering analyses for  <dig> individuals of the pyrenean morales grasshopper for different numbers of genetic clusters  based on a discriminant analysis of principal components .
each individual corresponds to a vertical bar partitioned into k-colored segments that represent the individual’s probability of belonging to the cluster with that color.
black lines separate individuals from different populations.
cluster colors in k =  <dig> barplot correspond to those of fig.
4
the result of the nj tree based on dc genetic distances showed that populations were included into monophyletic groups geographically clustered according to the hierarchical structure inferred by structure and dapc analyses .
multiple matrix regression with randomization
genetic differentiation  was positively associated with geographic distance , topographic complexity  and current , lgm  and lig  habitat resistance distances when these variables were included alone into different models .
indeed, mantel tests showed that all these variables were highly inter-correlated .
climatic dissimilarity  was also correlated with all other variables, but with comparatively lower mantel r values, whereas elevation dissimilarity  was only significantly correlated with climdis .
univariate models including ibd and ibrtc provided the highest and most similar coefficients of determination  .
however, only ibrtc was included into the final model after the backward selection procedure .
the rest of variables did not remain significant when they were tested against topographic complexity .
analyses based on fstna yielded similar results, but models had generally lower values of coefficient of determination  .
our results remained similar after sequential bonferroni correction for multiple testing  <cit> .table  <dig> results of univariate matrix regressions with randomization for genetic differentiation  among eleven populations of the pyrenean morales grasshopper in relation with elevation  and climatic  dissimilarity and five isolation by resistance  scenarios: ibd, isolation by distance ; ibrtc, topographic complexity; ibrcurrent, current habitat suitability; ibrlgm, last glacial maximum habitat suitability and ibrlig, last interglacial habitat suitability
f
st
f
r
2
t
p
r
2
t
p
ibd
ibr
tc
ibr
current
ibr
lgm
ibr
lig
elevdis
climdis
predictors with p <  <dig>  in bold
fig.
<dig> relationship between genetic differentiation  and topographic complexity  resistance distances   in eleven populations of pyrenean morales grasshopper.
regression line is shown
distance-based redundancy analysis
marginal tests showed a significant association between genetic differentiation  and geography and one environmental predictor  that explained  <dig> – <dig>  % and  <dig> – <dig>  % of genetic variation, respectively .
climate pc <dig> was explained by a pool of bioclimatic variables related to annual temperature variation .
however, the environmental predictor  remained not significant after accounting for the influence of geography in the conditional test .
our results remained similar after sequential bonferroni correction  <cit> .
analyses based on fstna-values  or performed with pca considering only the six bioclimatic layers employed to build the climate niche model in maxent provided similar results .table  <dig> results of distance-based redundancy analyses  testing the effects of geography, climate and elevation on genetic differentiation among eleven populations of the pyrenean morales grasshopper
f
p
f
p
f
st
geography
climate pc2
f
geography
climate pc2
geography is tested after transforming the euclidean geographical distance matrix to a continuous rectangular vector by principal coordinates analyses .
in marginal tests, we tested each predictor separately, while in conditional  tests geography was always included as covariate.
the proportion of the multivariate genetic variation explained  by a given predictor or set of predictors is indicated.
predictors with p <  <dig>  in bold
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 .
the bayesian posterior estimates of the parameter βg  were smaller than βe , indicating a primordial effect of geographical distances and an absence of isolation by climate or elevation on genetic differentiation .table  <dig> results of bayesian inference and model selection in sunder testing the relative effects of geographical and environmental variables on genetic differentiation among eleven populations of the pyrenean morales grasshopper
− <dig>
elevdis
− <dig>
we separately tested the environmental variables  against an ibd resistance distance matrix .
‘g’ corresponds to models only considering geography, ‘e’ corresponds to models only considering the environmental variable, and ‘g + e’ corresponds to models considering both of them.
for all runs, we show the likelihood of each model based on the validation dataset and the values of β
i parameter.
the β
i parameter quantifies the magnitude of the effect of each variable on genetic covariance .
the most likely model for each comparison is in bold
discussion
despite its small distribution range , 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  and much higher than that shown for ramburiella hispanica at a similar spatial scale , two species presenting a patchy distribution restricted to highly isolated and fragmented habitats  <cit> .
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  <cit> .
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  <dig> 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  <dig> kya according to our niche models  could be linked to historical events  predating the temporal scale addressed in our study  <cit> .
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  <cit> .
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  <cit> .
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  <cit> .
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 .
the short generation time of the studied species  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  <cit> .
we found no support for ecological divergence and local adaptation processes have contributed to population genetic differentiation and the three different employed approaches  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  <cit> .
our results contrasts with other studies that have documented an important role of environment on genetic differentiation in many taxa  <cit> , including some insects such as grasshoppers  <cit> , walking sticks  <cit>  or beetles  <cit> .
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 .
the meta-analysis by  <cit>  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  <cit> .
hence, the lack of effects of environment on gene flow could be due to different reasons, including low environmental variation among sampling sites  <cit>  or local adaptation driven by other ecological factors not considered in this study .
CONCLUSIONS
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 .
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  <cit> .
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  <cit> .
ethics
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 .
all other data supporting the results of this article are included within the article and its additional files.
consent to publish
not applicable.
additional file
additional file 1: table s <dig>   genetic differentiation between eleven populations of the pyrenean morales grasshopper .
table s <dig>  results of mantel tests analyzing the relationship between the different distance matrices  used to evaluate the factors associated with population genetic differentiation in the pyrenean morales grasshopper .
table s <dig>  results of principal component analysis  applied to the values of the  <dig> present day bioclimatic variables obtained from the worldclim dataset.
figure s <dig>  the position in the environmental space  of all known populations of the pyrenean morales grasshopper .
figure s <dig>  results of bayesian clustering analyses in structure to determine the best-supported number of clusters in hierarchical analyses.
abbreviations
aucarea under the curve.
bicbayesian information criterion
bpbefore present
cnmclimate niche models
dapcdiscriminant analysis of principal components
dbrdadistance-based redundancy analyses
hwehardy-weinberg equilibrium
ibdisolation-by-distance
ibeisolation-by-environment
ibrisolation-by-resistance
kyathousand years ago
lgmlast glacial maximum
liglast interglacial
m.a.s.l.meters above sea level
mcmcmarkov chain monte carlo
mmrrmultiple matrix regressions with randomization
mtssmaximum training sensitivity plus specificity threshold
njneighbor-joining
pcaprincipal component analysis
pcoaprincipal coordinates analyses
rmseroot mean squared error
competing interests
the authors declare that they have no competing interests.
authors’ contributions
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.

 understanding the factors structuring genetic variation in natural populations is a paramount topic in evolutionary biology. 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. the isolation-by-distance model predicts that genetic differentiation increases with euclidean geographic distance because of limited dispersal and genetic drift. even though this classic model explains spatial patterns of genetic differentiation in a wide range of organisms, it does not consider more sophisticated information than straight-line geographic distances and assumes landscape homogeneity, an unrealistic scenario for most natural systems. recently, the emergence of landscape genetics has explicitly incorporated landscape complexity into the study of evolutionary processes, 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. 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. 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. recent research has revealed the ubiquity of ibe patterns, highlighting the importance of ecological factors in shaping genetic structure of natural populations. 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. indeed, different isolating mechanisms are not mutually exclusive and gene flow is often influenced by a combination of geographical and ecological factors. given that geography and environment are often highly correlated, disentangling their relative impact on population genetic differentiation harbors inherent analytical difficulties that have promoted the progressive development of more robust and accurate statistical methods. however, only a few studies have jointly considered the relative effects of geography, landscape composition and environmental heterogeneity on either landscape-level or range-wide patterns of genetic structure. the pyrenean morales grasshopper is a narrow-endemic subspecies belonging to the chorthippus binotatus group, exclusively distributed in central aragón and catalonia pyrenean mountains. 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 and 2400 m.a.s.l.. 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. 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. 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 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. in particular, we used multiple matrix regressions with randomization, distance-based redundancy analyses 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, ii) resistance distances based on topographic complexity and current and past climate suitability; and iii) altitudinal and environmental dissimilarity between populations. 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.fig. workflow summarizing the methodological approach employed in this study to analyze the relative contribution of isolation by distance, isolation by resistance, and isolation by environment in structuring genetic variation in the pyrenean morales grasshopper. the response variables and predictors considered for each analytical approach are indicated. hs: habitat suitability; lgm: last glacial maximum; lig: last interglacial; mmrr: multiple matrix regression with randomization; dbrda: distance-based redundancy analysis methods population sampling between and we collected individuals from populations of the pyrenean morales grasshopper. we aimed to sample populations throughout the entire distribution range of the species based on occurrence-data available in the literature and our own prospection of areas with potentially suitable habitats. the morales grasshopper is primarily distributed in the south side of the pyrenees and a small area located in the northeastern part of these mountains. the latter portion of the species distribution was intensively prospected during our field work but we only found a single population in the area. overall, we were able to collect individuals from populations located across almost the entire climatic and elevation gradient occupied by the species, including populations from a wide range of habitats and differing in up to ~ 900 m of elevation. 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  specimens were collected in public lands under license from ‘gobierno de aragón’, ‘generalitat de catalunya’ and ‘ordesa y monte perdido national park’.fig. population genetic structure in the pyrenean morales grasshopper. panel a shows sampling sites of the species and phylogenetic relationships among the populations inferred from a neighbor-joining tree based on cavalli-sforza and edwards chord distances. the tree was plotted on a topographic map of the pyrenees using the software gengis [ 103 ]. we downloaded topographic data from nasa shuttle radar topographic mission as 90 m resolution digital elevation model and subsequently transformed to arc-sec resolution for representation. panel b represents the results of genetic assignment of individuals of the pyrenean morales grasshopper based on the bayesian method implemented in structure. we performed hierarchical analyses for subsets of populations considering the most probable k-value inferred at the previous hierarchical level. each individual corresponds to a vertical bar partitioned into k-colored segments that represent the individual’s probability of belonging to the cluster with that color. black lines separate individuals from different populations. population codes as in table 1 n 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 microsatellites markers described in. we performed pcr amplifications following the procedure described in, run pcr products on an abi genetic analyzer and scored genotypes using genemapper. we tested for deviations from hardy-weinberg equilibrium using exact tests based on markov chain iterations and dememorization steps as implemented in the program arlequin. 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 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 fst-values between all pairs of sampling populations and testing their significance with fisher’s exact test after permutations using arlequin. we corrected p-values using a sequential bonferroni adjustment. due to the high frequency of null alleles in orthoptera, we also calculated pairwise fst-values corrected for null alleles using the so-called ena-method implemented in the program freena. we inferred genetic structure using bayesian clustering analyses in structure.. 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. 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 with a burn-in period of steps and a run length of markov chain monte carlo cycles. we estimated the best-supported number of genetic clusters with the log probability of the data and the Δk method. we used the ‘full search’ algorithm in the program clumpp. to align replicated runs and average individual assignment probabilities for the most likely k-value. finally, we used distruct 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 to identify clusters of genetically related individuals. although both clustering methods exhibit similarities, 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 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 as a prior step to a discriminant analysis. da partitions genetic variation maximizing differences between clusters while minimizing within-cluster variation. we implemented dapc analysis in r. using the package adegenet. at first, we used the ‘find.cluster’ function using all available principal components to determine the best-supported number of genetic clusters using the bayesian information criterion. the ‘find cluster’ function runs successive k-means clustering with increasing number of clusters and provides a bic value for each simulated k-value. secondarily, we determined the optimal number of pcs for the dapc by cross-validation using the ‘xvaldapc’ function with replicates. we selected the number of pcs associated with the lowest ‘root mean squared error’ 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. to obtain a neighbor-joining tree based on pairwise cavalli-sforza and edwards 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 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. based on current climate layers and using cross-validation replicate model runs. we used a total of occurrence points obtained from the literature and our own sampling. to construct the models, we used the present-day bioclimatic variables available in worldclim and downloaded at arc-sec resolution. to obtain the distribution of the pyrenean morales grasshopper in the last glacial maximum and the last interglacial, we projected contemporary species-climate relationships to these periods. the lgm layers were based on the community climate system model from the palaeoclimate modelling intercomparison project phase ii. we downloaded lgm layers from worldclim at arc-min and interpolated to arc-sec resolutions. lig layers were based on and downloaded from worldclim at 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 calculation to address the problems derived from projecting the current distribution into novel climates. 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. 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 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 and we compared its lgm and lig projections with those obtained by mess analyses. all the output maps from the models were visualized using threshold values based on maximum training sensitivity plus specificity. topographic complexity in order to investigate the role of topographic complexity 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. 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 and the final layer was transformed to arc-sec 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. 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 with ‘varimax’ rotation applied to the values of the present-day bioclimatic variables from worldclim extracted from sampling sites, occurrence points used in maxent and 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, and  % respectively of the environmental variance. 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. isolation by resistance matrices we applied circuit theory to model gene flow between populations and test the effects of different landscape resistance scenarios on observed patterns of genetic differentiation. we used circuitscape 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 by calculating pairwise resistance distances on a completely ‘flat’ landscape based on a raster layer in which all cells had an equal value. 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. multiple matrix regression with randomization all ibr matrices were tested against genetic distance matrices using multiple matrix regressions with randomization as implemented in r. in these models, we also included elevation and climate dissimilarity 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. distance-based redundancy analysis complementary to mmrr analyses, we tested the relationship between genetic differentiation, geography and environment using distance-based redundancy analyses. 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 was tested against the following variables: i) geographic distances, ii) elevation and iii) population’s pc scores of the first three axes from the pca performed on the environmental data. euclidean geographical distances between sampled populations were calculated using geographic distance matrix generator.. geographic distances were tested after transforming the euclidean geographical distance matrix to a continuous rectangular vector by principal coordinates analyses using the ‘pcnm’ function in the package vegan. significance of the predictors was assessed using multivariate f-statistics with 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 and then we performed a partial dbrda for each variable while controlling for the influence of geography. 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. 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 best fits to the data. using the multinomial model, we ran sunder with million of iterations for each data set, sampling every iterations. we set to update in the mcmc iterations all parameters of the vector θ. 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 as validation set during the cross-validation procedure. we performed sunder analyses using as environmental matrices the elevation dissimilarity matrix and the climate dissimilarity distance matrix, which we separately tested against ibr distance matrix based on a completely ‘flat’ landscape. 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, mean temperature of the coldest quarter, annual precipitation, precipitation of the driest month, precipitation seasonality, and precipitation of the warmest quarter. this model had a very high value of area under the receiving operator characteristics curve, 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. 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 800– 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. 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. climate niche modelling of the pyrenean morales grasshopper for the present, the last glacial maximum and the last interglacial. climatically suitable areas defined using the maximum training sensitivity plus specificity threshold is in red. the topography is in the background map, in which whiter colors indicate higher elevations and darker colors indicate lower elevations. yellow dots represent the eleven sampling sites population genetic structure pairwise fst-values ranged from to and of comparisons were significant after sequential bonferroni correction. pairwise fstna-values were lower than fst-values and ranged from to. pairwise fst-values were highly correlated with pairwise fstna-values. structure analyses on all populations showed that the best-supported number of clusters was k =  according to the Δk method. these initial analyses detected a strong correspondence between the inferred genetic clusters and their geographic location, even when a broader range of k-values was evaluated. subsequent hierarchical analyses on different subsets of populations detected further genetic structuring. 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. 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–. considering k = , dapc partitioned all individuals in western, central-west, central-east and eastern groups. 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. the low overlapping of the genetic clusters on the ordination plot indicated high degree of differentiation between them. when k =  and k =  were considered and compared with k =  dapc revealed a hierarchical distribution of genetic variation similar to that identified by structure.fig. summary of the results of discriminant analysis of principal components. panel a represents the bayesian information criterion for each value of k. the minimum value of bic before the first increase or stabilization indicates the best-supported number of genetic clusters. panel b represents an ordination plot for the first two discriminant axes. each dot represents one individual and colors and inertia ellipses indicate their assignment to one of the four genetic clusters inferred by dapc. the up-right graph inset displays the variance explained by the principal component axes used for dapc. the bottom-right inset displays in relative magnitude the variance explained by the two discriminant axes plotted. panel c represents whether the individuals were correctly assigned to the genetic cluster where they were included a priori by k-means analyses used to infer the best-supported clustering solution. colors represent membership probabilities to each genetic cluster and blue crosses indicate the cluster where the individuals were originally assigned by k-means analyses. generally, the dapc classification of individuals is consistent with their assignment to the clusters originally identified by k-means analyses. panel d shows the number of individuals from each population assigned to each of the four inferred genetic clusters. the size of black squares is proportional to the number of individuals assigned to the different clusters. population codes are described in table 1 fig. results of clustering analyses for individuals of the pyrenean morales grasshopper for different numbers of genetic clusters based on a discriminant analysis of principal components. each individual corresponds to a vertical bar partitioned into k-colored segments that represent the individual’s probability of belonging to the cluster with that color. black lines separate individuals from different populations. cluster colors in k =  barplot correspond to those of fig. 4 the result of the nj tree based on dc genetic distances showed that populations were included into monophyletic groups geographically clustered according to the hierarchical structure inferred by structure and dapc analyses. multiple matrix regression with randomization genetic differentiation was positively associated with geographic distance, topographic complexity and current, lgm and lig habitat resistance distances when these variables were included alone into different models. indeed, mantel tests showed that all these variables were highly inter-correlated. climatic dissimilarity was also correlated with all other variables, but with comparatively lower mantel r values, whereas elevation dissimilarity was only significantly correlated with climdis. univariate models including ibd and ibrtc provided the highest and most similar coefficients of determination. however, only ibrtc was included into the final model after the backward selection procedure. the rest of variables did not remain significant when they were tested against topographic complexity. analyses based on fstna yielded similar results, but models had generally lower values of coefficient of determination. our results remained similar after sequential bonferroni correction for multiple testing.table results of univariate matrix regressions with randomization for genetic differentiation among eleven populations of the pyrenean morales grasshopper in relation with elevation and climatic dissimilarity and five isolation by resistance scenarios: ibd, isolation by distance; ibrtc, topographic complexity; ibrcurrent, current habitat suitability; ibrlgm, last glacial maximum habitat suitability and ibrlig, last interglacial habitat suitability f st f r 2 t p r 2 t p ibd ibr tc ibr current ibr lgm ibr lig elevdis climdis predictors with p <  in bold fig. relationship between genetic differentiation and topographic complexity resistance distances in eleven populations of pyrenean morales grasshopper. regression line is shown distance-based redundancy analysis marginal tests showed a significant association between genetic differentiation and geography and one environmental predictor that explained –  % and –  % of genetic variation, respectively. climate pc was explained by a pool of bioclimatic variables related to annual temperature variation. however, the environmental predictor remained not significant after accounting for the influence of geography in the conditional test. our results remained similar after sequential bonferroni correction. analyses based on fstna-values or performed with pca considering only the six bioclimatic layers employed to build the climate niche model in maxent provided similar results.table results of distance-based redundancy analyses testing the effects of geography, climate and elevation on genetic differentiation among eleven populations of the pyrenean morales grasshopper f p f p f st   geography   climate pc2 f   geography   climate pc2 geography is tested after transforming the euclidean geographical distance matrix to a continuous rectangular vector by principal coordinates analyses. in marginal tests, we tested each predictor separately, while in conditional tests geography was always included as covariate. the proportion of the multivariate genetic variation explained by a given predictor or set of predictors is indicated. predictors with p <  in bold 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. the bayesian posterior estimates of the parameter βg were smaller than βe, indicating a primordial effect of geographical distances and an absence of isolation by climate or elevation on genetic differentiation.table results of bayesian inference and model selection in sunder testing the relative effects of geographical and environmental variables on genetic differentiation among eleven populations of the pyrenean morales grasshopper − elevdis − we separately tested the environmental variables against an ibd resistance distance matrix. ‘g’ corresponds to models only considering geography, ‘e’ corresponds to models only considering the environmental variable, and ‘g + e’ corresponds to models considering both of them. for all runs, we show the likelihood of each model based on the validation dataset and the values of β i parameter. the β i parameter quantifies the magnitude of the effect of each variable on genetic covariance. the most likely model for each comparison is in bold discussion despite its small distribution range, 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 and much higher than that shown for ramburiella hispanica at a similar spatial scale, two species presenting a patchy distribution restricted to highly isolated and fragmented habitats. 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 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 kya according to our niche models could be linked to historical events 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. 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. 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. the short generation time of the studied species 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. we found no support for ecological divergence and local adaptation processes have contributed to population genetic differentiation and the three different employed approaches 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. our results contrasts with other studies that have documented an important role of environment on genetic differentiation in many taxa, including some insects such as grasshoppers, 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. 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.  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. 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. 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. ethics 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. all other data supporting the results of this article are included within the article and its additional files. consent to publish not applicable. additional file additional file 1: table s genetic differentiation between eleven populations of the pyrenean morales grasshopper. table s results of mantel tests analyzing the relationship between the different distance matrices used to evaluate the factors associated with population genetic differentiation in the pyrenean morales grasshopper. table s results of principal component analysis applied to the values of the present day bioclimatic variables obtained from the worldclim dataset. figure s the position in the environmental space of all known populations of the pyrenean morales grasshopper. figure s results of bayesian clustering analyses in structure to determine the best-supported number of clusters in hierarchical analyses. abbreviations aucarea under the curve. bicbayesian information criterion bpbefore present cnmclimate niche models dapcdiscriminant analysis of principal components dbrdadistance-based redundancy analyses hwehardy-weinberg equilibrium ibdisolation-by-distance ibeisolation-by-environment ibrisolation-by-resistance kyathousand years ago lgmlast glacial maximum liglast interglacial m.a.s.l.meters above sea level mcmcmarkov chain monte carlo mmrrmultiple matrix regressions with randomization mtssmaximum training sensitivity plus specificity threshold njneighbor-joining pcaprincipal component analysis pcoaprincipal coordinates analyses rmseroot mean squared error competing interests the authors declare that they have no competing interests. authors’ contributions 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. 
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. 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. 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. 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. 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 and local adaptation processes in the evolution of spatial genetic structure of the morales grasshopper, 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, topographic complexity and present and past habitat suitability models, and environmental dissimilarity. different clustering analyses revealed a deep genetic structure that was best explained by ibr based on topographic complexity. chorthippus saulcyi moralesiecological niche modellinghierarchical genetic structureisolation by environmentisolation by resistancelandscape geneticspyreneestopographic complexityhttp://dx.doi.org/ /501100003329ministerio de economía y competitividadcgl2011-25053cordero pedro j. junta de comunidades de castilla-la manchapoii10-0197-0167peii-2014-023-pcordero pedro j. european regional development funduncm08-1e-018cordero pedro j. issue-copyright-statement© the author 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 and local adaptation processes in the evolution of spatial genetic structure of the morales grasshopper, 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, topographic complexity and present and past habitat suitability models, and environmental dissimilarity.  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. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords chorthippus saulcyi moralesiecological niche modellinghierarchical genetic structureisolation by environmentisolation by resistancelandscape geneticspyreneestopographic complexityhttp://dx.doi.org/ /501100003329ministerio de economía y competitividadcgl2011-25053cordero pedro j. junta de comunidades de castilla-la manchapoii10-0197-0167peii-2014-023-pcordero pedro j. european regional development funduncm08-1e-018cordero pedro j. issue-copyright-statement© the author 2016 