Testing the role of ancient and contemporary landscapes on structuring genetic variation in a specialist grasshopper

3110 | Ecology and Evolution. 2017;7:3110–3122.www.ecolevol.org Received: 4 October 2016 | Revised: 31 December 2016 | Accepted: 24 January 2017 DOI: 10.1002/ece3.2810 ORIGINAL RESEARCH Testing the role of ancient and contemporary landscapes on structuring genetic variation in a specialist grasshopper Víctor Noguerales1 | Pedro J. Cordero1 | Joaquín Ortego2 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2017 The Authors. Ecology and Evolution published by John Wiley & Sons Ltd. 1Grupo de Investigación de la Biodiversidad Genética y Cultural, Instituto de Investigación en Recursos Cinegéticos - IREC (CSIC, UCLM, JCCM), Ciudad Real, Spain 2Department of Integrative Ecology, Estación Biológica de Doñana (EBD-CSIC), Seville, Spain Correspondence Víctor Noguerales, Grupo de Investigación de la Biodiversidad Genética y Cultural, Instituto de Investigación en Recursos Cinegéticos - IREC (CSIC, UCLM, JCCM), Ciudad Real, Spain. Email: victor.noguer[email protected] Funding information Ministerio de Economía y Competitividad, Grant/Award Number: BES-2012-053741; Severo Ochoa, Grant/Award Number: SEV2012-0262; Ramón y Cajal, Grant/Award Number: RYC-2013-12501; Ministerio de Economía y Competitividad and European Social Fund, Grant/Award Number: CGL201125053 and CGL2014-54671-P; Junta de Comunidades de Castilla-La Mancha and European Social Fund, Grant/Award Number: POII10-0197-0167 and PEII-2014023-P; European Regional Development Fund, Grant/ Award Number: UNCM08-1E-018 Abstract Understanding the processes underlying spatial patterns of genetic diversity and structure of natural populations is a central topic in evolutionary biogeography. In this study, we combine data on ancient and contemporary landscape composition to get a comprehensive view of the factors shaping genetic variation across the populations of the scrublegume grasshopper (Chorthippus binotatus binotatus) from the biogeographically complex region of southeast Iberia. First, we examined geographical patterns of genetic structure and employed an approximate Bayesian computation (ABC) approach to compare different plausible scenarios of population divergence. Second, we used a landscape genetic framework to test for the effects of (1) Late Miocene paleogeography, (2) Pleistocene climate fluctuations, and (3) contemporary topographic complexity on the spatial patterns of population genetic differentiation. Genetic structure and ABC analyses supported the presence of three genetic clusters and a sequential westtoeast splitting model that predated the last glacial maximum (LGM, c. 21 Kya). Landscape genetic analyses revealed that population genetic differentiation was primarily shaped by contemporary topographic complexity, but was not explained by any paleogeographic scenario or resistance distances based on climate suitability in the present or during the LGM. Overall, this study emphasizes the need of integrating information on ancient and contemporary landscape composition to get a comprehensive view of their relative importance to explain spatial patterns of genetic variation in organisms inhabiting regions with complex biogeographical histories. KEYWORDS Bayesian inference, climate niche modeling, genetic diversity, genetic structure, isolation by resistance, topographic complexity 1 | INTRODUCTION Understanding the mechanisms that shape spatial patterns of genetic diversity and structure is a central topic in evolutionary biogeography (Habel et al., 2015; Peterman, Connette, Semlitsch, & Eggert, 2014; Yannic et al., 2014). Presentday landscape configuration and the geographical distribution of suitable habitats, jointly with speciesspecific ecological characteristics, define contemporary interpopulation dispersal and realized gene flow (Castillo, Epps, Davis, & Cushman, 2014; Edwards, Keogh, & Knowles, 2012). However, past climate changes and ancient geological events have also greatly altered the spatial configuration of corridors and barriers to gene flow (He, Edwards, & Knowles, 2013; Pepper, Doughty, Arculus, & Keogh, 2008). Such temporal shifts in landscape structure and dispersal routes have often | 3111 NOGUERALES Et AL. left genetic signatures in contemporary populations that are useful to track back in time their past demographic trajectories (He et al., 2013; Lanier, Massatti, He, Olson, & Knowles, 2015). Thus, the study of how presentday and past landscape composition have impacted gene flow is necessary to get a comprehensive view of the processes underlying spatial patterns of genetic diversity and structure of natural populations, which can ultimately help to predict their responses to ongoing or future environmental changes (Fordham, Brook, Moritz, & NoguésBravo, 2014; Yannic et al., 2014). Quaternary climatic fluctuations, characterized by cold glacial stages alternated with warm interglacial periods, have strongly influenced the demography of many organisms during the past 2 million years (2–0.04 Mya) (Hewitt, 2000). During glacial periods, the distribution ranges of most species from temperate zones contracted and their populations persisted in refugia located at lower elevations or latitudes (Homburg et al., 2013; Qu et al., 2014). Conversely, the populations from cooladapted species expanded during glacial periods and shrank during interglacials (Canestrelli & Nascetti, 2008). Under any scenario, populations from regions subjected to major climate changes experience fluctuating demographic dynamics that, ultimately, are expected to reduce their effective population sizes, increase genetic drift, and erode local levels of genetic diversity (Brown & Knowles, 2012; Carnaval, Hickerson, Haddad, Rodrigues, & Moritz, 2009; Yannic et al., 2014). However, populations from climatically unstable areas can recurrently go extinct and be recolonized by immigrants from multiple source populations, which can increase local levels of genetic diversity via admixture (Ortego, Gugger, & Sork, 2015; Petit et al., 2003). Thus, the stability of climatically suitable habitats can impact patterns of genetic diversity and admixture in opposite directions, a possibility that has been generally overlooked (Ortego, Gugger, et al., 2015). Beyond Quaternary climatic fluctuations and contemporary landscape features, much older paleogeological events such as uplifting of mountain ranges or the emergence of islands and sea corridors are also considered important factors responsible of geographical patterns of genetic differentiation in many taxa (Ceccarelli et al., 2016; MastrettaYanes, MorenoLetelier, Pinero, Jorgensen, & Emerson, 2015; Papadopoulou, Anastasiou, Keskin, & Vogler, 2009). Although ancient geological changes are known to underlie the spatial patterns of genetic divergence found in several organisms (Abellán, Arribas, & Svenning, 2012; Cheng et al., 2016; Opell, Helweg, & Kiser, 2016; Ortego, Bonal, Cordero, & Aparicio, 2009), in many other cases the genetic signals left by paleogeological events are expected to have been totally or partially eroded as a result of gene flow promoted by subsequent landscape changes (Graham, Hendrixson, Hamilton, & Bond, 2015; Pepper et al., 2008). Thus, examining landscape configuration at different time periods can help to better understand the mechanisms by which intraspecific genetic diversity and differentiation arise and are maintained (Reilly, Corl, & Wake, 2015). The mountainous area of southeast Iberia has undergone remarkable geological changes that have shaped the complex biogeographical history of the region. Geological reconstructions based on stratigraphic and sedimentary data show that the emergence of mountain chains in the Tortonian (c. 12 Mya) configured a mosaic of islands (hereafter Betic Islands) at the confluence of European and African continental platforms. The rotation of the Betic Islands toward the Iberian Peninsula, in combination with sedimentation processes, resulted in their fusion to the continent and the configuration of a continuous emerged landscape that is currently conformed by the Prebetic, Penibetic, and Subbetic mountain ranges of southeast Iberia (Braga, Martín, & Quesada, 2003; Braga et al., 2010; Martín, Braga, Aguirre, & Puga-Bernabéu, 2009). Furthermore, this area is also the southernmost limit of the influence of Quaternary glaciations (c. 2–0.04 Mya) in Europe, during which vast portions of land were free of permanent ice at elevations below 2,500 m.a.s.l. (Hughes & Woodward, 2008) and constituted an important refugium for biota from temperate habitats (Hewitt, 2000). The magnitude and complexity of these paleogeological and climate events are considered the most important engines of diversification and genetic structuring of many taxa in the region (Andújar, GómezZurita, Rasplus, & Serrano, 2012; Faille, Andujar, Fadrique, & Ribera, 2014; Fromhage, Vences, & Veith, 2004). For all these reasons, southeast Iberia is an ideal template for testing the combined effects of ancient and more contemporary climate and landscape changes on spatial patterns of genetic diversity and structure of local populations (Faille et al., 2014). The scrublegume grasshopper (Chorthippus binotatus binotatus Charpentier, 1825) (Orthoptera: Acrididae) is a winged Orthoptera with a 1year generation time (Figure 1; Defaut, 2011). This species is primarily distributed in montane regions from southwest Europe, including France and the Iberian Peninsula (Defaut, 2011). The scrublegume grasshopper is an oligophagous species that exclusively feeds on some scrublegume taxa from the tribe Genisteae (Defaut, 2011). In southeast Iberia, the host plants (primarily Erinacea anthyllis and, more occasionally, Echinospartum boissieri, Genista versicolor, and Ulex parviflorus) form scattered vegetation patches located at moderate to high elevations (>1,200 m.a.s.l.). This fact restricts the distribution of the scrublegume grasshopper to the different mountain ranges of the region (Prebetic, Penibetic and Subbetic systems) (Defaut, 2011; FIGURE1 Scrublegume grasshopper (Chorthippus binotatus binotatus), the study organism. The photography shows a male specimen on a legume host plant of the genus Ulex (tribe Genisteae). Photography by Víctor Noguerales 3112 | NOGUERALES Et AL. Table S1). Thus, the narrow ecological requirements of the scrublegume grasshopper, the patchy distribution of its host plants, and the limited dispersal abilities of the species (low flying capacity; V.N., P.J.C and J.O., pers. obs.) have resulted in most of its populations from southeast Iberia being currently highly fragmented and separated by extensive lowlands of unsuitable habitats (Defaut, 2011). Here, we used the scrublegume grasshopper as model system to analyze the contribution of contemporary (presentday topography and distribution of climatically suitable habitats) and historical (paleoclimatebased distribution of suitable habitats and Late Miocene paleogeography) factors on shaping spatial patterns of genetic diversity and structure across the populations of the species from southeast Iberia. In particular, we first (1) examined geographical patterns of genetic structure and employed an approximate Bayesian computation (ABC) framework to compare different plausible scenarios of population divergence (Beaumont, 2010; Cornuet et al., 2014). Second, we (2) applied circuit theory to test whether observed patterns of genetic differentiation are explained by a comprehensive suite of isolationbyresistance (IBR) scenarios (McRae, 2006; McRae & Beier, 2007), including paleogeography at different time periods since Late Miocene (c. 12.0–7.0 Mya; Martín et al., 2009), current and last glacial maximum (LGM, c. 21 Kya) climate suitability and stability, and contemporary topographic complexity (TC). Finally, we (3) tested the hypothesis predicting more genetic diversity in populations from areas with high past and present climate suitability and stability since the LGM. 2 | MATERIAL AND METHODS 2.1 | Population sampling In 2012 and 2013, we collected 354 individuals from 19 populations of scrublegume grasshopper from southeast Iberia (~80,000 km2) (Table S1; Figures 2–4). Our sampling included populations from all mountain ranges in the region (Prebetic, Penibetic, and Subbetic ranges) and covered the entire elevation range of the scrublegume grasshopper in the study area (958–2,314 m.a.s.l.; Table S1). This allowed us to sample populations from different habitats such as alpine and Mediterranean scrublegume formations. Specimens were collected using a butterfly net, and the whole body was preserved in 2ml vials with 96% ethanol and stored at –20°C until needed for DNA extraction. Our sampling was performed under licenses from the “Junta de Comunidades de CastillaLa Mancha,” “Junta de Andalucía,” and “Gobierno de la Región de Murcia.” Population codes and more information on sampling sites are presented in Table S1. 2.2 | Microsatellite genotyping and basic genetic statistics We extracted genomic DNA from a hind leg of each individual using a salt extraction protocol (Aljanabi & Martinez, 1997). Each individual was genotyped at 18 speciesspecific microsatellites markers (Basiita et al., 2016). All microsatellite markers were polymorphic in all populations, and the most common alleles were shared across all populations. We performed PCRs and genotyping following the procedure described in Ortego, Aguirre, Noguerales, and Cordero (2015) and Basiita et al. (2016). We tested for deviations from Hardy– Weinberg equilibrium, linkage disequilibrium (LD), and the presence of null alleles as described in Noguerales, Cordero, and Ortego (2016). Two loci (Cbin16 and Cbin36) were discarded from all downstream FIGURE2 Paleogeographic maps showing the spatial configuration of emerged lands in the study area during the (a) Early Tortonian (c. 12.0–11.6 Mya), (b) Late Tortonian (c. 8.0–7.3 Mya), and (c) Earliest Messinian (c. 7.2–7.0 Mya) according to Martín et al. (2009). Yellow dots indicate the location of sampled populations (number codes as in Table S1). Dashed lines represent continental limits in the present. Inset map from panel (a) shows the location of our study area within the Iberian Peninsula (a) (b) (c) North-Betic strait Guadalquivir basin Mediterranean Sea Atlantic Ocean Guadalhorce gateway 1 2 3 45 6 78 9 10 11 12 13 14 15 16 17 19 18 1 2 3 45 6 7 8 9 10 11 12 13 14 15 16 17 19 18 1 2 3 45 6 78 9 10 11 12 13 14 15 16 17 19 18 | 3113 NOGUERALES Et AL. analyses because of HW disequilibrium in all populations and the presence of null alleles. We did not find evidence for LD between any pair of loci in any sampling population after sequential Bonferroni corrections (Rice, 1989). 2.3 | Analyses of genetic structure We estimated population genetic differentiation calculating FST values between all pairs of sampling populations. Significance of genetic differentiation between all pairs of populations was tested with Fisher’s exact tests after 10,000 permutations using ARLEqUiN 3.5 (Excoffier & Lischer, 2010). pValues were corrected using a sequential Bonferroni adjustment (Rice, 1989). Due to the frequent presence of null alleles in Orthoptera (Keller, Holderegger, & van Strien, 2013), we also calculated pairwise FST values corrected for null alleles (FSTNA) using the socalled ENA method implemented in the program FREENA (Chapuis & Estoup, 2007). We inferred genetic structure using Bayesian clustering analyses in StRUctURE 2.3.3 (Falush, Stephens, & Pritchard, 2003; Pritchard, Stephens, & Donnelly, 2000). We considered correlated allele frequencies and an admixture model without prior information on population origin. We performed 10 independent runs for each value of assumed number of genetic clusters (K = 1–12) with a burnin period of 200,000 steps and a run length of 1,000,000 Markov chain Monte Carlo cycles. The number of genetic clusters (K) best fitting the data set was defined using log probabilities [Pr(X|K)] (Pritchard et al., 2000) and the ΔK method (Evanno, Regnaut, & Goudet, 2005). We used the Greedy algorithm in the program cLUmpp 1.1.2 (Jakobsson & Rosenberg, 2007) to align replicated runs and average individual assignment probabilities for the most likely K values. Finally, we used DiStRUct 1.1 (Rosenberg, 2004) to produce bar plots displaying probabilities of individual membership to each inferred genetic cluster. We also examined the spatial genetic structure considering geographical coordinates of sampling sites as a priori information in the Bayesian clustering method implemented in tESS 2.3.1 (Chen, Durand, Forbes, & Francois, 2007; Durand, Jay, Gaggiotti, & François, 2009). We used the conditional autoregressive (CAR) Gaussian model of admixture with a linear trend surface, updating the spatial interaction parameter (ψ), initially set to the default value 0.99. The variance term (initially set to 1) permitted to update during the course of runs. CAR model was chosen in order to avoid overestimation of the most likely K in the presence of genetic clines (François & Durand, 2010; Guillot, 2009). We ran 20 independent replicates for each value of K = 2–12 using 50,000 sweeps of which 10,000 were used as burnin period. The best supported number of genetic clusters (K) was estimated using the deviance information criterion (DIC) values and stabilization of the Q matrix of posterior probabilities (Chen et al., 2007; Gao, Bryc, & Bustamante, 2011). For each KMAXvalue considered, we conducted 180 additional replicate runs up to a total of 200 replicates. We used the 10 runs with the lowest DIC values to align and average individual assignment probabilities with cLUmpp before being represented using DiStRUct as indicated above for StRUctURE analyses. Complementarily, we constructed a phylogenetic tree to visualize the genetic relationships between all populations. We used the program pOpULAtiONS 1.2.31 (Langella, 1999) to obtain a neighborjoining tree based on pairwise CavalliSforza and Edwards (Dc) genetic distances (CavalliSforza & Edwards, 1967). Finally, we carried out analyses of molecular variance (AmOvAs) to examine the partitioning of FIGURE3 Climate niche modeling for the scrublegume grasshopper in southeast Iberia for (a) the present and (b) the last glacial maximum (LGM, c. 21 Kya). Panel (c) shows climate stability estimated as the sum of pixel values of current and LGM climate suitability maps. The LGM maps represent the average climate suitability index of the projections obtained from CCSM and MIROC climate models. Gray scales refer to climate suitability (range: 0–1) and climate stability (range: 0–2), with increasingly darker shades of gray indicating increasing climate suitability and stability. Inset map from panel (a) shows the location of our study area within the Iberian Peninsula (a) (b) (c) Prebetic system Penibetic system Subbetic system 3114 | NOGUERALES Et AL. the genetic variation among and within regions and populations as defined by five population grouping hypotheses. Populations were pooled according to their historical location in the three different paleogeographical Late Miocene scenarios (see Figure 2 and section “Landscape genetic analyses”) and their current location in the main mountain ranges of the region (Prebetic, Penibetic, and Subbetic systems; see Table S1 and Figure 3a). Additionally, we tested the grouping scheme used for ABC analyses (see next section). AmOvAs were performed in ARLEqUiN 3.5 (Excoffier & Lischer, 2010), and the significance of the variance components was tested using 10,000 permutations of the original data. 2.4 | Approximate Bayesian computation In order to infer the evolutionary and demographic history of the scrublegume grasshopper in the region, we compared four plausible scenarios of population divergence using an ABC approach (Beaumont, 2010). To simplify the analyses, we defined three main groups (groups A, B, and C) of populations by pooling sampling sites according to their geographical location and the results from AmOvAs (Table S2) and Bayesian clustering analyses (StRUctURE and tESS) (e.g., Tsuda, Nakao, Ide, & Tsumura, 2015). Note that although K = 2 was the most supported clustering solution for both StRUctURE and tESS analyses, K = 3 revealed further hierarchical genetic substructure with geographical coherence (see the “Results” section and Figure S1). In group A, we included populations 1–6 (western populations); group B, 7–12 (southeastern populations); and group C, 13–19 (northeastern populations) (see Figures 4 and 5). The topology of each scenario was designed considering the connectivity of populations according to Bayesian clustering analyses (Figure 4 and Figure S1c). The scenarios tested were the following: (1) Scenario I, null model: The three groups diverged simultaneously; (2) Scenario II, sequential splitting model from west to east: Group A split from group B and C at t2, and these two groups subsequently split at t1; (3) Scenario III, sequential splitting model from east to west: Group C split from groups A and B at t2, and these two groups split at t1; (4) Scenario IV, splitting model from FIGURE4 Sampling sites of scrublegume grasshoppers and genetic structure based on Bayesian clustering analyses. Pie charts on the map represent the genetic assignments for each sampling population according to StRUctURE analyses. For each population, left and right pie charts represent the admixture proportions considering K = 2 and K = 3, respectively. Circle size is proportional to the number of genotyped individuals in each population. Code numbers are described in Table S1. On the bottom, barplots represent the assignment of individuals to each genetic group according to tESS analyses considering K = 2 (top) and K = 3 (bottom). Each individual corresponds to a vertical bar, which is partitioned into Kcolored segments that represent the individual’s probability of belonging to the cluster with that color. Vertical black lines separate individuals from different populations. On the right, neighborjoining tree based on CavalliSforza and Edwards chord distances. Colors are according to StRUctURE analyses based on K = 3. Inset map shows the location of our study area within the Iberian Peninsula 1 19 17 16 18 15 ´ ñ 14 5 12 11 10 9 8 7 2 6 3 4 13 16.Cabras 19.Espuña 14.Poyotello 15.Pinar 12.María 17.Carrasca 18.Almerara 13.Cazorla 7.Otero 8.Ragua 9.Gador 11.Filabres 10.Baza 5.Mágina 4.Pandera 2.Tejeda 6.Arana 1.Ronda 3.Parapanda ´ | 3115 NOGUERALES Et AL. central to peripheral populations: Group B split from groups A and C at t2, and these two groups subsequently split at t1 (Figure 5). We conducted all the computations using DiyAbc 2.0.4 (Cornuet et al., 2014). We generated 3 millions of simulated data sets per scenario considering a generalized mutation model and no singlenucleotide indels (Table S3). The summary statistics (SS) used in ABC analyses are described in Table S3. We performed preevaluation of scenarios and prior distributions in DiyAbc to adjust the priors of Ne and t to their most appropriate values (see Table S3), assuming a uniform prior probability distribution for them. To avoid biases in parameter estimates, we selected the subset of seven microsatellites markers with lower frequency of null alleles as estimated in the program FREENA. Selection of the most probable scenario, confidence in scenario choice (type I and II errors), model checking, and estimation of the posterior distribution of all parameters under the best supported model were performed as described in Ortego, Noguerales, Gugger, and Sork (2015). 2.5 | Landscape genetic analyses We applied circuit theory (McRae, 2006; McRae & Beier, 2007) and a multiple matrix regression with randomization (MMRR) approach (Wang, 2013) to examine the relative contribution of a suite of IBR scenarios to explain patterns of genetic differentiation in our study populations. Specifically, we tested nine different hypothetical scenarios of population connectivity, which included (1) three paleogeographic scenarios defined by the spatial configuration of emerged lands at different time periods (Early Tortonian, Late Tortonian, and Earliest Messinian); (2) three scenarios based on the distribution of climatically suitable habitats since the LGM (current climate suitability, LGM climate suitability, and climate suitability stability since the LGM); (3) a scenario of population connectivity defined by contemporary TC; (4) an isolationbydistance (IBD) scenario representing the geographical distance between each pair of populations. Below we describe in detail the methods followed to generate these scenarios and test their relative contribution to contemporary patterns of genetic differentiation. 2.5.1 | Paleogeographic scenarios To test the possible effect of the complex geological history of the study region on contemporary patterns of genetic differentiation, we considered three paleogeographic scenarios: Early Tortonian (12.0– 11.6 Mya), Late Tortonian (8.0–7.3 Mya), and Earliest Messinian (7.2–7.0 Mya) (see Figure 2). We used ARcGiS 10.0 (ESRI, Redlands, CA, USA) to create vector layers for emerged lands based on geological models from Martín et al. (2009). Then, we transformed vectors layers for each scenario into raster maps with a 30 arcsec (c. 1 km) resolution that were finally used as inputs in ciRcUitScApE (McRae, 2006; McRae & Beier, 2007) (see below for details on ciRcUitScApE analyses). 2.5.2 | Climatic suitability scenarios We modeled the potential climate distribution of scrublegume grasshopper at different time periods to investigate whether the spatial distribution of climatically suitability habitats are relevant factors shaping observed patterns of genetic differentiation in the study populations. For this purpose, we built a climate niche model (CNM) using the maximum entropy presenceonly algorithm implemented in mAxENt 3.3.3 (Phillips, Anderson, & Schapire, 2006; Phillips & Dudik, 2008) based on current climate. We used a total of 85 occurrence points obtained from the Global Biodiversity Information Facility, the literature (Defaut, 2011), and our own sampling. To construct the models, we used the 19 bioclimatic variables available in WorldClim and downloaded at 30 arcsec (c. 1 km) resolution (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005). Variables retained in the final models were selected following several complementary criteria (Vega et al., 2010). At first, we used ENMtOOLS (Warren, Glor, & Turelli, 2010) to examine colinearity among variables, in order to retain a single layer among those with a high Pearson correlation coefficient (r > .85). Then, we used the Jackknife of regularized training gain procedure implemented in mAxENt to retain the variables with the maximum contribution to the model. We discarded the worst and highest correlated predictors among the whole set of variables, conducted a new model with the remaining variables, and repeated this backward process until the final model only retained the best explanatory and less correlated variables (Vega et al., 2010). Model evaluation statistics were produced from 10 crossvalidation replicate model runs. To obtain the distribution of scrublegume grasshopper during the LGM (LGM, c. 21 Kya), we projected contemporary species–climate relationships to the LGM using two atmospheric circulation models: the Community Climate System Model (CCSM3; Collins et al., 2006) and the Model for Interdisciplinary Research on Climate (MIROC 3.2; Hasumi & Emori, 2004) from the Paleoclimate Modelling Intercomparison Project Phase II (PMIP2; Braconnot et al., 2007). LGM layers were downloaded from WorldClim at 2.5 arcmin and interpolated to 30 arcsec resolution. To reduce the level of uncertainty FIGURE5 Scenarios compared using an approximate Bayesian computation (ABC) approach (t# represents time in number of generations; N# represents effective population sizes during each time period) Group C (northeastern) Group A (western) Group B (southeastern) 3116 | NOGUERALES Et AL. arising from different past projections, we averaged climate suitability scores from projections based on CCSM and MIROC models to obtain a consensus LGM map of climatically suitable areas. In addition, we summed current and LGM climate suitability layers to generate a map of climate suitability stability, with pixel values ranging from 0 (minimum climate suitability in both periods) to 2 (maximum climate suitability in both periods). All GIS calculations were conducted in ARcGiS10.0. Finally, current, LGM, and stability climate suitability raster maps were used as inputs in ciRcUitScApE to calculate IBR distance matrices (see below for details). 2.5.3 | Contemporary topographic complexity scenario We investigated the role of contemporary TC as a potential factor shaping patterns of genetic differentiation in our study populations. We calculated the surface ratio index for each cell from a presentday digital elevation model using “DEM SURFAcE tOOLS” (Jenness, 2013) in ARcGiS 10.0. Surface ratio is an index of TC, with values close to one indicating flat areas and values higher than one indicating a more abrupt relief with deeper slopes (Jenness, 2004). Calculations were conducted on a 90m resolution digital elevation model from NASA Shuttle Radar Topographic Mission (SRTM Digital Elevation Data). Although no information is available on the dispersal distance and home range of the study species, the high resolution of the digital elevation model is expected to capture well the TC relevant for a mediumsize grasshopper with a suspected low dispersal ability. The final raster map was transformed to 30 arcsec (c. 1 km) resolution and used as input in ciRcUitScApE (see below for details). 2.5.4 | CirCuitsCape analyses We used ciRcUitScApE 4.0 (McRae, 2006; McRae & Beier, 2007) to calculate resistance distance matrices between all pairs of populations considering an eightneighbor cell connection scheme. The raster layers generated for the nine different scenarios of population connectivity were used as inputs in ciRcUitScApE. For the three paleogeographic scenarios, the raster layers included two element classes: “emerged land” and “sea water.” We considered that “sea water” was the main landscape feature limiting the dispersal of terrestrial fauna in the study area during these periods. We generated different IBR scenarios assigning different resistance values to “sea water” (10, 50, 100, 500, 1,000, 10,000), an approach that allowed us to identify the optimal ratio of landscape resistance between both landscape elements that best fit our data on genetic differentiation (e.g., Andrew, Ostevik, Ebert, & Rieseberg, 2012; Ortego, Aguirre, et al., 2015). To test the effect of IBD, we calculated pairwise resistance distances on a completely “flat” landscape based on a raster layer in which all cells had an equal value (conductance = 1). This IBD resistance model is expected to yield similar results than a matrix of Euclidean geographical distances, but it is more appropriate for comparison with others competing models also generated with ciRcUitScApE (VeloAntón, Parra, ParraOlea, & Zamudio, 2013). 2.5.5 | Statistical analyses We used a MMRR approach to examine the relative contribution of all IBR and IBD scenarios to explain patterns of genetic differentiation in our study populations (Wang, 2013). We tested the two matrices of genetic differentiation (FST and FSTNA) against all pairwise resistance distance matrices representing the nine different IBR/IBD scenarios. We used a backward procedure to select final models, eliminating nonsignificant 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 (Noguerales et al., 2016; Ortego, Gugger,et al., 2015). 2.6 | Analyses of genetic diversity and admixture Allelic richness (AR) standardized for sample size was calculated for each population using HpRARE (Kalinowski, 2005). We estimated the genetic admixture of populations using a genetic admixture index (GADMIX) obtained from the probabilities of population membership to each genetic cluster inferred by StRUctURE analyses (Ortego, Gugger, et al., 2015). This index was designed to standardize the degree of genetic admixture across populations with different probabilities of membership to different genetic clusters (Ortego, Gugger, et al., 2015), and its advantages and potential caveats are those inherent to StRUctURE analyses (Falush et al., 2003; Pritchard et al., 2000). GADMIX ranges from 0 (indicating no admixture, i.e., genetically pure populations assigned to a single genetic cluster) to 1 (indicating maximum admixture, i.e., genetically admixed populations with an equal probability of membership to each inferred genetic cluster). We used generalized linear models (GLMs, using a Gaussian error distribution and an identity link function) and an informationtheoretic model selection approach to analyze AR and GADMIX (Burnham & Anderson, 2002). Models for AR included as independent variables current climate suitability (HSCUR), LGM climate suitability (HSLGM), and climate suitability stability (HSSTA). Models for GADMIX included HSSTA as independent variable. Longitude and latitude were included as additional covariates in models for both AR and GADMIX to take in account possible geographical clines of genetic diversity and admixture (e.g., Guo, 2012). We calculated average HSCUR, HSLGM, and HSSTA with ARcmAp 10.0 at different spatial scales using buffers of 1, 10, and 100 km2 around sampling locations. Given that the precision of AR and GADMIX estimates may differ among populations due to differences in sample sizes, we used a weighted least square method where weight equals the sample size for each studied population. GLMs were built in the R package LmE4 (Bates, Maechler, Bolker, & Walker, 2015; R Core Team, 2015) and model selection and averaging were performed using the R package mUMiN (Barton, 2015) as detailed in Noguerales, Traba, Mata, and Morales (2015) and Ortego, Aguirre, et al. (2015). 3 | RESULTS 3.1 | Population genetic structure We found that most pairs of populations were genetically differentiated. In particular, 139 of 171 pairwise FST values (~81%) were significantly | 3117 NOGUERALES Et AL. higher than zero after sequential Bonferroni correction (Table S4). Significant pairwise FST values ranged from .025 to .164, whereas pairwise FSTNAvalues were slightly lower and ranged from .008 to .144. Pairwise FST and FSTNAvalues were highly correlated (Mantel r = .993; p < .001). The populations from Ronda, Parapanda, and Tejeda, located at the westernmost portion of the study area, exhibited the highest levels of genetic differentiation with the rest of populations. Analyses in StRUctURE showed a best supported number of clusters for K = 2 according to the ΔK method. The first cluster included the Western populations, whereas the second cluster included the remaining populations located in the east part of the study area. However, log probabilities [Ln Pr (X|K)] steadily increased from K = 2 to K = 5 (Figure S1a). Individual assignment probabilities to a certain genetic cluster were moderately high up to K = 5 and the spatial distribution of genetic variation exhibited geographical consistency, but most populations showed a considerable degree of genetic admixture (Figure 4, Figure S1c). Genetic clustering analyses in tESS resulted in an optimal K = 4 according to the DIC criterion (Figure S1b), but one of the inferred clusters represented a “ghost cluster” with no individual assigned to it (see Chen et al., 2007; Guillot, Estoup, Mortier, & Cosson, 2005). When K = 3 was considered, the first cluster included the Western populations, the second cluster included the southeastern populations, and the third cluster included the northeastern populations of the study area (Figure 4). tESS and StRUctURE analyses yielded similar results for K = 2 and K = 3 (Figure S1c). The result of the neighborjoining tree based on CavalliSforza and Edwards chord distances (Dc) was also congruent with the results from Bayesian clustering analyses (Figure 4). Finally, AmOvA analyses indicated most genetic variance was attributed to differences within populations (>90%, for all grouping hypotheses; Table S2). The population grouping hypothesis that explained the highest percentage of total variation attributed to differences among groups was the one used for ABC analyses (Table S2). 3.2 | Approximate Bayesian computation The scenario considering sequential population divergence from west to east (scenario II) had the highest posterior probability based on both direct and logistic regressionbased estimates, and its 95% confidence interval did not overlap with those obtained for others scenarios that showed much lower support (Table 1). Observed data fell within simulated data (all SS Ps > .2) for scenario II, suggesting good model fit. Type I and II errors were .394 and .372, respectively, and RMAE values were moderate in most cases (Table 1). Considering the 1year generation time of scrublegume grasshopper (Defaut, 2011), the western genetic group (group A) diverged from the southeastern and northeastern groups (groups B and C, respectively) ~215,000 years ago (t2) (95% CI: 67,600–342,000 years ago), whereas these two groups split ~42,000 years ago (t1) (95% CI: 5,540–165,000 years ago) (Table 2). Assuming a constant mutation rate, the posterior estimates of effective populations sizes (Ne) indicated no important demographic changes after the different splitting events (Table 2). 3.3 | Climate niche modeling The variables included in the final CNM were temperature seasonality (BIO4), mean temperature of the driest quarter (BIO9), precipitation of the driest month (BIO14), and precipitation of the warmest quarter (BIO18). This model had a very high value of area under the curve (AUC; 0.982 ± 0.007), indicating overall good performance. The predicted distribution of scrublegume grasshopper in the present (Figure 3a) is consistent with its observed fragmented distribution. The distribution of the species was more extensive and populations were better connected during the LGM than at present time (Figure 3b). However, populations located in the western portion of our study area (corresponding with group A in ABC analyses) have remained highly isolated during both the LGM and the present. Scenario Posterior probability 95% CI Type I error Type II error I.0230 [0.0222–0.0253] II .9618 [0.9595–0.9640] .394 .372 III .0059 [0.0054–0.0064] IV .0085 [0.0078–0.0093] TABLE1 Posterior probability for each of the four tested scenarios and 95% confidence intervals (CI) based on the weighted polychotomous logistic regression approach for approximate Bayesian computation (ABC) analyses. Type I and type II errors for the best supported scenario (in bold) are indicated TABLE2 Posterior parameter estimates (median and 95% confidence intervals) for the best supported scenario (scenario 2, see Figure 5). Estimates are based on 1% of simulated data sets closest to the observed values. Relative median absolute errors (RMAE) based on 500 pseudoobserved data sets are also indicated for each parameter Parameter Median q0.025 q0.975 RMAE N1 520,000 194,000 736,000 .279 N2 407,000 113,000 715,000 .246 N3 600,000 277,000 740,000 .258 N1anc 303,000 39,800 690,000 .394 N2-3 345,000 46,600 695,000 .372 Nx57,000 24,800 541,000 .384 t142,700 5,540 165,000 .404 t2215,000 67,600 342,000 .213 μ7.75 × 10−6 4.24 × 10−6 2.85 × 10−5 .357 N1, effective population size of group A; N2, effective population size of group B; N3, effective population size of group C; N1anc, effective population size of the ancestral group A; N2-3, effective population size of the ancestral groups BC; Nx, effective population size of the most ancestral population, t1, time (in generations = years) to the most recent divergence event; t2, time (in generations = years) to the most ancient divergence event (see scenarios in Figure 5); μ, mean mutation rate. 3118 | NOGUERALES Et AL. 3.4 | Landscape genetic analyses Tejeda and Ronda populations are currently located in areas that were not likely to form emerged lands during Early and Late Tortonian, and for this reason, they were excluded from landscape genetic analyses. Thus, we performed MMRR analyses using the 17 populations presumably located on permanently emerged lands since the late Miocene in order to make our landscape genetic analyses comparable across all tested scenarios and time periods. Considering these 17 populations, only resistance distances based on contemporary TC and IBD were significantly associated with genetic differentiation (Ps < .006) (Table 3). However, only TC was retained into the final model (β = .826, t = 7.56, p = .004) (Figure 6). Analyses based on FSTNA gave similar results (Table 3), but models had slightly lower values of r2. Analyses considering all populations (n = 19) yielded qualitatively analogous results (data not shown). 3.5 | Analyses of genetic diversity and admixture GADMIX [K = i, 5] indexes obtained considering different K values were highly correlated among them (all r > .598, all Ps < .007). Population genetic admixture based on any GADMIX [K = i, 5] index and AR were also correlated (all r > .451, all Ps < .050). Model selection results showed that AR was not significantly associated with longitude, latitude, or HSCUR, HSLGM, or HSSTA at any analyzed spatial scale (all unconditional 95% CIs of the predictors crossed zero; Table S5). Likewise, population genetic admixture (based on any GADMIX [K = i, 5] index) was not significantly associated with longitude, latitude, or HSSTA at any analyzed spatial scale (all unconditional 95% CIs of the predictors crossed zero; Table S6). 4 | DISCUSSION Assuming ecological niche stability through time, our niche model revealed a moderate shift in the distribution of climatically suitable habitats for the scrublegume grasshopper in southeast Iberia during the last 21,000 years (NoguésBravo, 2009). Climate niche modeling indicated that the potential distribution of the species is more fragmented in the present than during the LGM, a pattern congruent with the increased population connectivity during glacial periods inferred for many other montane species from temperate regions (BlancoPastor, FernándezMazuecos, & Vargas, 2013; VeloAntón et al., 2013). Continuous climatically suitable habitats connected the southern foothills of the Prebetic mountain range and the eastern portion of the Penibetic system during the LGM, an area where the species is not present today as verified by our own surveys. Our CNM also outlined the isolation of the Western populations since LGM, which has probably contributed to shape their strong genetic differentiation with the rest of the populations within the study area. The isolation of the Western populations could have occurred during the last interglacial (LIG) period (c. 120,000–140,000 years ago), when the scrublegume grasshopper probably showed a distribution similar to that in Model FST FSTNA r2βt p r2βt p IBD .298 .824 7.543 .004 .273 .811 7.097 .011 TC .299 .826 7.561 .006 .274 .814 7.118 .021 HSCUR .031 .163 2.084 .264 .055 .224 2.815 .086 HSLGM .220 .446 6.16 .051 .194 .429 5.662 .064 HSSTA .210 .432 5.978 .053 .185 .417 5.521 .062 Early Tortonian .115 .320 4.180 .063 .096 .301 3.780 .055 Late Tortonian .101 .296 3.886 .066 .084 .278 3.522 .052 Earliest Messinian .088 .278 3.612 .065 .072 .258 3.237 .062 For paleogeographical models, we considered high resistance values for sea water (=100) and low for emerged lands (=1). Table shows the results based on the 17 populations presumably located on permanently emerged lands since the Late Miocene. TABLE3 Results of univariate matrix regressions with randomization (MMRR) for genetic differentiation [FST and FST corrected for null alleles (FSTNA)] in relation to different isolationbyresistance (IBR) scenarios: geographical distance (IBD), contemporary topographic complexity (TC), current climate suitability (HSCUR), last glacial maximum climate suitability (HSLGM), climate suitability stability (HSSTA), and three paleogeographical models (Early Tortonian, c. 12.0–11.6 Mya; Late Tortonian, c. 8.0–7.3 Mya; Earliest Messinian, c. 7.2–7.0 Mya) FIGURE6 Relationship between genetic differentiation (FST) and resistance distances calculated using ciRcUitScApE on the basis of contemporary topographic complexity (TC) 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00 0.20 0.60 1.00 1.40 1.80 Genetic differentiation (F ST ) TC resistance distance