Consequences of extensive habitat fragmentation in landscape-level patterns of genetic diversity and structure in the Mediterranean esparto grasshopper

Anthropogenic habitat fragmentation has altered the distribution and population sizes in many organisms worldwide. For this reason, understanding the demographic and genetic consequences of this process is necessary to predict the fate of populations and establish management practices aimed to ensure their viability. In this study, we analyse whether the spatial configuration of remnant semi-natural habitat patches within a chronically fragmented landscape has shaped the patterns of genetic diversity and structure in the habitat-specialist esparto grasshopper (Ramburiella hispanica). In particular, we predict that agricultural lands constitute barriers to gene flow and hypothesize that fragmentation has restricted interpopulation dispersal and reduced local levels of genetic diversity. Our results confirmed the expectation that isolation and habitat fragmentation have reduced the genetic diversity of local populations. Landscape genetic analyses based on circuit theory showed that agricultural land offers ∽1000 times more resistance to gene flow than semi-natural habitats, indicating that patterns of dispersal are constrained by the spatial configuration of remnant patches of suitable habitat. Overall, this study shows that semi-natural habitat patches act as corridors for interpopulation gene flow and should be preserved due to the disproportionately large ecological function that they provide considering their insignificant area within these human-modified landscapes.


Introduction
Anthropogenic habitat fragmentation has altered the distribution and population sizes in many organisms worldwide and can be considered one of the major threats to biodiversity (Noss and Csuti 1994;Fahrig 2002;Lindenmayer and Fischer 2006). As a result of this process, the genetic connectivity and diversity of populations has been often negatively impacted (Frankham 1996;DiBattista 2008). Fragmentation of formerly continuous habitats can modify the dispersal routes of organisms in such a way that only suitable remnant habitats act as corridors for genetic exchange among populations (Fahrig 2007;e.g. Pavlacky et al. 2009;Jha and Kremen 2013). In turn, population subdivision and disruption of gene flow can reduce levels of genetic diversity (Frankham 1996;e.g. Levy et al. 2013;M endez et al. 2014) and, ultimately, compromise the viability of populations and lead to local extinctions (Frankham 2005; e.g. Saccheri et al. 1998). For this reason, understanding the genetic consequences of habitat fragmentation is useful to predict the fate of populations and establish management practices aimed to preserve their evolutionary potential and ensure their long-term persistence (DiBattista 2008).
The integration of genetic and spatial data can help to determine whether postfragmentation habitat configuration has modified dispersal patterns and identify the landscape elements that most contribute to genetic connectivity (Manel et al. 2003;Storfer et al. 2007Storfer et al. , 2010. In the absence of natural barriers to dispersal, the movement of animals is expected to be mostly determined by geographical distance in prefragmentation continuous landscapes but constrained by the spatial distribution of corridors of suitable habitat patches after fragmentation (Zellmer and Knowles 2009;Jha and Kremen 2013). However, the responses of species to landscape fragmentation are difficult to predict and highly dependent on their dispersal capacity and propensity to cross unsuitable habitat patches (Blanchet et al. 2010;DiLeo et al. 2010;Lange et al. 2010). Accordingly, some studies have found that fragmentation results in population connectivity is constrained by corridors of suitable habitat embedded within a hostile habitat matrix (e.g. Pavlacky et al. 2009;Zellmer and Knowles 2009;Jha and Kremen 2013;Ruiz-Gonz alez et al. 2014), whereas others have revealed that human-driven landscape fragmentation has no effect (Qu em er e et al. 2010) or even facilitates dispersal and gene flow (Bacles et al. 2005;Pavlova et al. 2012). Landscape genetics can help to identify discontinuities in gene flow and determine the relative resistance to movement imposed by different landscape elements, offering a powerful tool to understand the impacts of habitat fragmentation and inform on ground conservation practices aimed to maintain or promote population connectivity (Segelbacher et al. 2010).
The Mediterranean region has been modified over centuries of logging and land clearing for grazing and agriculture, constituting one of the areas of the world historically most altered by humans (Blondel and Aronson 1999). As a result of this process, large areas of formerly forested regions have been transformed into agricultural lands and nonfarmed areas are now often reduced to relict semi-natural habitat patches (Blondel and Aronson 1999). This is the case of La Mancha, a region from central Spain that constitutes the largest plain from the country and where land is currently devoted in its great majority to extensive cultures of wheat, barley, vineyards and olive trees. Semi-natural habitats are mostly reduced to a few scattered small hills and canyons, rock outcrops, and saline low grounds not suitable for agriculture. In this study, we analyse whether the configuration of these remnant semi-natural habitat patches within this chronically fragmented landscape has shaped the patterns of genetic diversity and structure in the habitat-specialist esparto grasshopper (Ramburiella hispanica) (Rambur, 1838). This Mediterranean orthoptera exclusively inhabits esparto grass formations (e.g. Stipa tenacissima and Lygeum spartum), which in La Mancha region are ubiquitous plant species in any remnant semi-natural nonagricultural land. We first characterize the habitat of the species in the study area and use circuit theory to generate different isolation-by-resistance (IBR) scenarios and test the relative importance of geographical distance and the distribution of suitable habitats on contemporary patterns of genetic differentiation (McRae 2006;McRae and Beier 2007;McRae et al. 2008). Second, we examine whether habitat fragmentation has negatively impacted local levels of genetic diversity (Frankham 1996;e.g. M endez et al. 2014;Levy et al. 2013). In particular, we (i) predict that agricultural lands constitute barriers to gene flow and (ii) hypothesize that fragmentation has restricted interpopulation dispersal and reduced local effective population sizes, resulting in populations located in more fragmented and isolated habitat patches have lower levels of genetic diversity.

Study sites and sampling
During 2010, we sampled 352 individuals from 18 populations of esparto grasshoppers in La Mancha region, central Spain (~2500 km 2 ; Fig. 1). We aimed to sample 20 specimens and an equal number of males and females in each locality, but samples sizes are often male-biased due to the difficulties to capture females in some sites (Table 1). Specimens were preserved whole in 1500 lL of 96% ethanol and stored at À20°C until needed for genetic analyses. Population code descriptions and further information on sampling localities are given in Table 1 and Fig. 1.

Microsatellite genotyping and basic genetic statistics
We used a salt extraction protocol to purify genomic DNA from a hind leg of each individual (Aljanabi and Martinez 1997). We used 12 highly polymorphic microsatellite markers to genotype each sampled individual (Aguirre  Figure 1 Geographical location of the studied populations of esparto grasshopper (Ramburiella hispanica) (red dots) and land-cover types classified as 'natural vegetation' (noncultivated areas adequate for the species) and 'agricultural land' (mostly cultivated areas and, in a much lesser extent, lagoons, reservoirs, villages and other developed areas). Population codes are described in Table 1.  Table S1). Amplifications were conducted in 10-lL reaction volumes containing approximately 20 ng of template DNA, 19 reaction buffer (67 mM Tris-HCL, pH 8.3, 16 mM (NH 4 ) 2 SO 4 , 0.01% Tween-20, EcoStart Reaction Buffer; Ecogen, Madrid, Spain), 2 mM MgCl 2 , 0.2 mM of each dNTP, 0.15 lM of each dye-labelled primer (FAM, PET, VIC or NED) and 0.1 U of Taq DNA EcoStart Polymerase (Ecogen). The PCR programme used was 9 min denaturing at 95°C followed by 40 cycles of 30 s at 94°C, 45 s at the annealing temperature (Table S1) and 45 s at 72°C, ending with a 10-min final elongation stage at 72°C. Amplification products were electrophoresed using an ABI 310 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA), and genotypes were scored using GENEMAPPER 3.7 (Applied Biosystems).
Microsatellite loci were tested for departure from Hardy-Weinberg equilibrium within each sampling population using an exact test (Guo and Thompson 1992) based on 900 000 Markov chain iterations as implemented in the program ARLEQUIN 3.1 (Excoffier et al. 2005). We also used ARLEQUIN 3.1 to test for linkage disequilibrium between each pair of loci for each sampling population using a likelihood-ratio statistic, whose distribution was obtained by a permutation procedure (Excoffier et al. 2005). We applied sequential Bonferroni corrections to account for multiple comparisons (Rice 1989).

GIS analyses
In our study area, the esparto grasshopper distributes in areas with formations of the host plants S. tenacissima and L. spartum, which are ubiquitous plant species present in any noncultivated habitat patch from La Mancha region. The esparto grasshopper has never been recorded in agricultural lands or developed areas (Pardo and G omez 1995 and references therein). For this reason, we classified our study area in two landscape element classes: 'natural vegetation' (noncultivated semi-natural areas optimal for the species) and 'agricultural land' (mostly cultivated areas and, in a much lesser extent, lagoons, reservoirs, villages and other developed areas unsuitable for the species). This information was obtained by digitalizing habitat patches from the most recent aerial pictures available at Centro Nacional de Informaci on Geogr afica from Spain (http:// centrodedescargas.cnig.es/CentroDescargas/). We used ARCMAP 10.0 (ESRI, Redlands, CA, USA) to create a vector layer that was then transformed into a raster grid (pixel size = 50 m) to be used in subsequent landscape genetic analyses (see below).

Genetic diversity
For each population, we calculated allelic richness (A R ) standardized for sample size using the program HP-RARE (Kalinowski 2005) and gene diversity (H E ) using F STAT (Goudet 1995). We used an information-theoretic modelselection approach to analyse which variables contribute to explain patterns of genetic diversity (A R and H E ) across the studied populations of esparto grasshopper. We considered four explanatory variables: (i) the proportion of suitable habitat (i.e. 'natural vegetation'; see previous section) within a 1 km radius of the sampling site, which we hypothesize to be associated with local effective population sizes (N e ) (see also Levy et al. 2013 for a similar approach); (ii) average genetic differentiation (F ST ) of each population with all other populations, an estimate of population isolation (e.g. Wang et al. 2011;Ortego et al. 2012); (iii) latitude; and (iv) longitude. We used general linear models (GLM) with a normal error structure and identity link function as implemented in the R 3.0.0 package LME4 (R Core Team 2012). The precision of A R estimates may differ among populations due to small differences in sample sizes, and we took this into account using a weighted least square method, where weight equals the sample size for each studied population. We ranked the resulting models following a model-selection approach on the basis of the Akaike's information criterion corrected for small sample size (AIC c ; Burnham and Anderson 1998). AIC c values for each model were rescaled (DAIC c ) calculating the difference between the AIC c value of each model and the minimum AIC c obtained among all competing models (i.e. the best model has DAIC c = 0). Models with DAIC c ≤ 2 were considered equivalent (Burnham and Anderson 1998). In cases where model selection as a function of AIC c did not give a single model, we performed an averaging of equivalent models (i.e. with DAICc ≤ 2; Burnham and Anderson 2002). We calculated the mean of the predictor estimators, their unconditional standard errors (SE) and confidence intervals (CIs), and the relative importance of each variable in the final averaged model (Σ xi, the sum of Akaike weights of models with DAIC c ≤ 2 in which the variable was included). The effect of predictor variables was considered significant if the 95% CI of their estimators did not cross zero (Burnham and Anderson 2002). Model selection and averaging was performed using the R package AICCMODAVG (R Core Team 2012).

Genetic structure
We estimated genetic differentiation between populations calculating pairwise F ST values and testing their significance with Fisher's exact tests after 10 000 permutations as implemented in ARLEQUIN 3.1 (Excoffier et al. 2005). Critical P-values for pairwise tests of allelic differentiation were determined using a sequential Bonferroni adjustment (Rice 1989). Due to the frequent presence of null alleles in grasshoppers (e.g. Chapuis et al. 2008;Blanchet et al. 2012a), we used the program FREENA to estimate null allele frequencies and calculate pairwise F ST values corrected for null alleles using the so-called ENA method (Chapuis and Estoup 2007). We used the randomization method implemented in F STAT 2.9.3 (10 000 permutations) to analyse differences between males and females in interpopulation genetic differentiation (F ST ), deviation from Hardy-Weinberg equilibrium (F IS ), mean assignment index (mAIc) and variance of the assignment index (vAIc) (Goudet et al. 2002), parameters that are informative about sex-biased patterns of dispersal (e.g. Ortego et al. 2011). Finally, we analysed genetic structure using the Bayesian clustering analysis implemented in the program TESS 2.3.1, which incorporates geographical coordinates as a priori information (Chen et al. 2007;Durand et al. 2009). We ran TESS 2.3.1 under the conditional autoregressive (CAR) Gaussian model of admixture with a linear trend surface (Durand et al. 2009), which is expected to be more robust against overestimation of K max in the presence of genetic clines (Guillot 2009;Francois and Durand 2010). We conducted 20 independent replicates for each value of K = 2-12 using 50 000 sweeps and a burn-in period of 10 000 sweeps. The admixture parameter (a) and the interaction parameter (w) were initially set to a = 0.99 and w = 0.6. We used deviance information criterion (DIC) values and stabilization of the Q-matrix of posterior probabilities to determine the optimal number of clusters (i.e. K max ) for the data. Once K max was deduced, 180 additional replicate runs were conducted to yield a total of 200 replicate runs for K max . We used the 10 runs with the lowest DIC values to calculate average individual admixture proportions with CLUMPP 1.1.2 (Jakobsson and Rosenberg 2007), which were visualized as a bar plot using DISTRUCT 1.1 (Rosenberg 2004).

Landscape genetic analyses
We applied circuit theory to model gene flow across a spatially heterogeneous landscape and determine the impact of isolation-by-distance (IBD) and different IBR scenarios on observed patterns of genetic differentiation (McRae 2006; e.g. McRae and Beier 2007). We used CIRCUITSCAPE 3.5.8 to calculate resistance distance matrices between all pairs of populations considering an eight-neighbour cell connection scheme (McRae 2006). We used the raster layer obtained as described in the section 'GIS analyses', which includes the two landscape element classes ('natural vegetation' and 'agricultural land') that a priori are the most likely to determine the distribution and dispersal patterns in our study species. We generated different IBR scenarios assigning a range of resistance values to both habitat classes (Table 2). This allowed us to identify the optimal ratio of landscape resistance between 'natural vegetation' and 'agricultural land' habitat classes that best fit our data of genetic differentiation (e.g. Andrew et al. 2012;Seymour et al. 2013). To test the effect of IBD, we calculated a matrix of Euclidean geographical distances between sampled populations using GEOGRAPHIC DISTANCE MATRIX GENERATOR 1.2.3 (Ersts 2011). We also generated a matrix of resistances in CIRCUITSCAPE considering an entirely 'flat' landscape, that is based on a raster layer in which all cells have equal resistance (resistance = 1). This matrix of flat resistance distances is expected to yield similar results than the matrix of Euclidean geographical distances, but the former has been suggested to be more appropriate for comparison with models of IBR generated with CIRCUITSCAPE (Lee-Yaw et al. 2009; Munshi-South 2012; Jha and Kremen 2013). Geographical distance and IBR matrices were tested against matrices of pairwise F ST values using a multiple matrix regression with randomization (MMRR) approach (Wang 2013). We used the MMRR function script implemented in R 3.0.2 (Wang 2013). Finally, we determined how well data on pairwise genetic differentiation fit the different IBR models (using the coefficient of determination, R 2 ) with varying levels of 'natural vegetation'/'agricultural land' resistance ratios (e.g. Andrew et al. 2012).

Microsatellite data
All microsatellite markers were highly polymorphic across all populations, with 12-40 alleles per locus (Table S1). After applying sequential Bonferroni corrections to compensate for multiple statistical tests, two loci (RhA113 and RhC1) consistently deviated from HWE across all the studied populations and were excluded from further analyses ( Table S1). The frequency of null alleles in the different loci ranged from moderately low (<0.10; RhA108, RhB107, RhC112, RhD2 and RhB2) to high (≥0.20; RhA105, RhA112, RhA2 and RhC2) values (Table S1). We did not find any evidence of genotypic linkage disequilibrium at any pair of loci in any population (exact tests; all Ps > 0.05).

Genetic diversity
A R and H E for each population are indicated in Table 1.
Our most parsimonious models revealed that both A R and H E increased significantly with the proportion of suitable habitat and decreased with isolation estimated as average genetic differentiation (F ST ) with other populations (Table  S2 and Table 3; Fig. 2). All other tested models for A R had a DAIC c value > 2 (Table S2). For H E , the model also including longitude was equivalent to the best ranked model (DAIC c value = 0.22; Table S2). However, unconditional CIs for longitude crossed zero, indicating that this variable had no significant effect (Table 3).

Genetic structure
Pairwise F ST values ranged from 0.0002 to 0.0469, and 73 of the 153 pairwise comparisons were significant after sequen- tial Bonferroni correction (Table S3). The randomization method implemented in F STAT showed that genetic differentiation (F ST ) (P = 0.671), deviation from Hardy-Weinberg equilibrium (F IS ) (P = 0.586), mAIc (P = 0.959) and vAIc (P = 0.348) did not differ between males and females. Individual-based analyses in TESS resulted in a K max = 4, but three genetic clusters were scarcely represented in most populations ( Fig. 3; Table S4).

Landscape genetic analyses
Euclidean geographical distances (IBD) and flat resistance distances were highly correlated (r = 0.99, P < 0.0001) and they were not significantly associated with genetic differentiation, either when they were included alone in the model (Table 2; Fig. 4A) or together with any IBR matrix (all Ps > 0.2). However, genetic differentiation was positively and significantly associated with those IBR matrices considering the minimum resistance for 'natural vegetation' habitat (=1) and high resistance values for 'agricultural land' (>50) ( Table 2; Fig. 4B). Model fit increased with the 'agricultural land'/'natural vegetation' resistance ratio, but the strength of the relationship between genetic differentiation and IBR stabilized beyond a ratio 1000:1 (Table 2; Fig. 5).
Analyses based on F ST corrected for null alleles gave very similar results, but best models reached slightly higher values of R 2 ( Fig. 5B; Table 2). Further analyses based on F ST values calculated only considering the five loci with low fre-quencies of null alleles (RhA108, RhB107, RhC112, RhD2 and RhB2; Table S1) provided analogous results (not shown), indicating that the effects of null alleles on the obtained results are minimal (e.g. Phillipsen et al. 2015).

Discussion
Our results confirmed our original expectations that landscape configuration has impacted gene flow and genetic diversity of esparto grasshopper, indicating that historical fragmentation of natural habitats as a result of land clearing for agriculture has shaped its dispersal patterns and local effective population sizes. Analyses of spatial patterns of genetic structure showed the presence of a shallow genetic differentiation in this species, with a high degree of genetic admixture and low F ST values ( Fig. 3; Table S3). The observed patterns of genetic structure in the esparto grasshopper contrast with those reported in other specialist grasshoppers showing deep genetic structures at similar or much shorter spatial scales (Ortego et al. 2011(Ortego et al. , 2012Keller et al. 2013a) and can be considered comparable to the levels of genetic differentiation found in widespread generalist orthopterans that are likely to be only moderately impacted by habitat fragmentation (Wiesner et al. 2011;Blanchet et al. 2012b;Keller et al. 2013b). The esparto grasshopper is a specialist species that exclusively inhabits S. tenacissima and L. spartum grass formations in the study area. However, these host plant species are ubiquitous in any seminatural habitat patch not devoted to agriculture, which may have contributed to maintain moderately high levels of gene flow within the study area despite a considerable landscape fragmentation. This contrasts with the scenario faced by the co-distributed grasshopper Mioscirtus wagneri, a habitat specialist with low dispersal capacity exclusively depending on the plant Suaeda vera for feeding (Ortego et al. 2010(Ortego et al. , 2012. This plant exclusively grows in saline low grounds, which are relict environments within the study area and submitted to a degree of fragmentation comparatively much higher than that experienced by the habitats occupied by the esparto grasshopper (Ortego et al. 2010). In contrast with other studies on orthopterans, we did not find any genetic signature of sex-biased dispersal in the esparto grasshopper (Bailey et al. 2007;Ortego et al. 2011;Kindler et al. 2012). Thus, the fact that this species has a relatively high flying capacity (J. Ortego and P. J. Cordero, unpublished data) and that both sexes disperse from their natal areas at a similar rates may have contributed to increase gene flow and resulted in subtle genetic structures at the landscape scale here analysed (Ortego et al. 2011).
Despite the high potential for gene flow in this species and the shallow patterns of genetic structure above described, the dispersal routes of the species have been constrained by the spatial configuration of remnant semi-natu- Table 3. General linear models (GLMs) for (a) standardized allelic richness (A R ) and (b) gene diversity (H E ). A single model with DAICc ≤ 2 was obtained for A R . For H E , we performed model averaging of the two best ranked equivalent models (DAICc ≤ 2) to obtain parameter estimates and unconditional standard errors (S.E.) (see Table S2 in Supporting information). The relative importance of each variable is indicated (Σ xi, sum of Akaike weights of models with DAICc ≤ 2 in which the variable was present). Bold type indicates significant variables, that is variables for which their unconditional 95% confidence interval (CI) did not cross zero. ral habitat patches within a matrix of land extensively devoted to agriculture. Landscape genetic analyses based on circuit theory showed that agricultural land offers~1000 times more resistance to gene flow than semi-natural habitats. The model exclusively including geographical distances between populations (i.e. IBD) was not significant and model fit increased with higher ratios of agricultural land/natural habitat resistance (Fig. 5), indicating that the subtle genetic differentiation observed within the study area probably reflects the impact of farming and land clearing on the species' demographic dynamics. However, agricultural land cannot be considered an absolute barrier to dispersal in this species, as many of the semi-natural habitat patches where our populations are located are highly isolated within a matrix of cultures (see Fig. 1) (Coulon et al. 2006). Lack of interpatch dispersal through agricultural lands should have resulted in stronger patterns of population genetic differentiation than we actually found, particularly if we consider the long time elapsed since fragmentation occurred, the short generation time of the species (1 year) and the fact that some of the studied habitat patches are too small and cannot sustain population sizes sufficiently large to buffer the effects of genetic drift (Frankham 1996). Analyses of genetic diversity have shown that the proportion of suitable habitat around sampling sites is positively associated with local levels of genetic diversity, suggesting that populations located in less fragmented habitat patches sustain higher effective population sizes (e.g. Levy et al. 2013;M endez et al. 2014). Further, genetic diversity was negatively correlated with average genetic differentiation with all other populations, indicating that isolation and limited gene flow have also contributed to erode genetic variability in some populations (e.g. Wang et al. 2011). This suggests that effective population sizes of the studied populations are not above a threshold that prevents the loss of genetic diversity and/or that the high potential for gene flow suggested by the subtle patterns of genetic structure observed in this species is not sufficient to counterbalance the effects of genetic drift (Lange et al. 2010).
In a broader context, the results of this study illustrate the importance of analysing the genetic consequences of extensive habitat destruction in ubiquitous species with high potential for gene flow. Understanding the consequences of habitat fragmentation in widespread species can help to determine the 'minimum' impact of this process for the entire community of related species (Lange et al. 2010;Keller et al. 2013b). Our study indicates that even widespread species that can persist in small habitat patches and show considerable gene flow also suffer from genetic drift and loss of genetic diversity due to habitat fragmentation. Thus, ubiquitous species less prone to suffer the effects of habitat fragmentation can inform on the landscape-level demographic processes experienced by other formerly codistributed species that may have already gone locally extinct and provide a late-warning signal of the genetic consequences of historical habitat fragmentation. For these reasons, information on pervasive species, generally disregarded in most conservation genetic studies (but see Lange et al. 2010;Keller et al. 2013b), can greatly contribute to define the ultimate genetic impacts of habitat fragmentation, establish management practices aimed to restore patch connectivity and evaluate the efficiency of conservation actions in many regions of the world that have historically experienced massive land clearing linked to agricultural practices.

Conclusions and implications for conservation
Our study indicates that a few remnant semi-natural habitat patches within a chronically and extensively fragmented landscape act as functional corridors that facilitate interpopulation gene flow and shape local levels of genetic diversity in the esparto grasshopper. These results, in conjunction with those of previous studies, indicate that the impact of land clearing for agriculture on dispersal patterns and gene flow can strongly vary even among related organisms that are expected to similarly respond to habitat fragmentation (Lange et al. 2010;Keller et al. 2013a,b;Levy et al. 2013). The preservation of these semi-natural patches may be particularly important for species with limited dispersal capacity and/or showing preferences for some microhabitats more geographically restricted. The extraordinarily deep genetic structure previously reported for the highly specialist M. wagneri suggests that local extinctions are not likely to be compensated by recurrent recolonizations in this species (Ortego et al. 2010(Ortego et al. , 2012, a pattern that considerably differs from the remarkable gene flow and metapopulation dynamics characterizing the more widespread R. hispanica (present study). Thus, a general recommendation derived from both this and previous studies in the area would be implementing management practices aimed to promote the conservation of organisms that are ecologically dissimilar, but prioritizing those species that are dispersal-limited and more likely to benefit from increasing or maintaining population connectivity (Gaublomme et al. 2011;Keller et al. 2013b). Given that most lands are private properties devoted to agriculture, management should Figure 3 Genetic structure of the studied populations of esparto grasshopper (Ramburiella hispanica). Figure shows the genetic assignment based on the Bayesian method implemented in the program TESS for different numbers of genetic clusters (K). Each individual is represented by a thin vertical line, which is partitioned into K-coloured segments that represent the individual's probability of belonging to the cluster with that colour. Black lines separate individuals from different populations. Population codes are described in Table 1. therefore focus on preserving existing natural habitat patches or enhancing dispersal through riparian corridors (Keller et al. 2013a). Although some habitats have been recently protected or proposed for protection in the study region, these initiatives have been up to now mostly focused on saline/hypersaline lagoons and low grounds of particular importance due to their unique plant and animal communities (e.g. Cirujano-Bracamonte and Medina-Domingo 2002; Cordero and Llorente 2008). However, less attention has been paid to other vulnerable habitats such as esparto grass formations that provide important ecological functions and contribute to maintain biodiversity in a disproportionate way, particularly if we consider the insignifi-cant area that they represent within these extensively human-modified landscapes (Manning et al. 2006). These habitats, often perceived as unproductive lands with no economic value, are submitted to different sources of human disturbance, including indiscriminate ploughing and aerial insecticide spraying for pest management, misleading habitat restoration practices (e.g. non-native pine plantations), uncontrolled waste dumping, and recurrent vegetation damage caused by livestock and off-road driving, among others (P. J. Cordero and J. Ortego, personal  Goodness of fit for models of landscape resistance considering different resistance ratios for agricultural land and natural vegetation. Panels show the coefficient of determination (R 2 ) for models analysing genetic differentiation (panel A: F ST ; panel B: F ST corrected for null alleles) in relation to isolation-by-resistance (IBR) distance matrices plotted against resistance ratios for 'agricultural land' and 'natural vegetation' used to calculate resistance-based distances with CIRCUITSCAPE. Resistance ratios for both habitat classes are log-transformed for illustrative purposes. Filled dots indicate significant models. observations). Thus, regional conservation policies aimed to avoid these practices together with environmental education activities to generate local people's awareness for the preservation of these remnants habitats would greatly contribute to their long-term conservation.

Acknowledgements
We wish to thank to Conchi C aliz for help in sample collection and genotyping. Two anonymous referees provided useful discussion and valuable comments on an earlier draft of this manuscript. JO was supported by Juan de la Cierva, Severo Ochoa (EBD) (SEV-2012-0262) and Ram on y Cajal Research Fellowships. MPA was supported by a technician contract from Ministerio de Ciencia e Innovaci on. VN was supported by a FPI predoctoral fellowship from Ministerio de Econom ıa y Competitividad. This work received financial support from grants CGL2011-25053 (Ministerio de Ciencia e Innovaci on and European Social Fund), POII10-0197-0167 (Junta de Comunidades de Castilla-La Mancha and European Social Fund) and UNCM08-1E-018 (European Regional Development Fund).

Supporting Information
Additional Supporting Information may be found in the online version of this article: Table S1. Microsatellite loci used to genotype esparto grasshoppers (Ramburiella hispanica). Table S2. Model selection to assess the association between population genetic diversity and cover of suitable habitat, average genetic differentiation (F ST ) with other populations, latitude and longitude. Table S3. Pair-wise population F ST and F ST corrected for null alleles. Table S4. Proportion of membership of the studied populations to each of the four genetic clusters identified by the Bayesian method implemented in the program TESS.