Insights into the neutral and adaptive processes shaping the spatial distribution of genomic variation in the economically important Moroccan locust (Dociostaurus maroccanus)

Ecology and Evolution. 2020;00:1–18. | 1www.ecolevol.org Received: 20 August 2019 | Revised: 12 February 2020 | Accepted: 18 February 2020 DOI: 10.1002/ece3.6165 ORIGINAL RESEARCH Insights into the neutral and adaptive processes shaping the spatial distribution of genomic variation in the economically important Moroccan locust (Dociostaurus maroccanus) María José González-Serna1 | Pedro J. Cordero1,2 | Joaquín Ortego3 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. © 2020 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 2Departamento de Ciencia y Tecnología Agroforestal y Genética, Escuela Técnica Superior de Ingenieros Agrónomos (ETSIA), Universidad de Castilla-La Mancha (UCLM), Ciudad Real, Spain 3Department of Integrative Ecology, Estación Biológica de Doñana – EBD – (CSIC), Seville, Spain Correspondence María José González-Serna, Grupo de Investigación de la Biodiversidad Genética y Cultural, Instituto de Investigación en Recursos Cinegéticos – IREC – (CSIC, UCLM, JCCM), Ronda de Toledo, 12, E-13071 Ciudad Real, Spain. Email: mariajose.gon[email protected] Funding information Junta de Comunidades de Castilla-La Mancha and European Social Fund, Grant/ Award Number: PEII-2014-023-P; Dirección General de Investigación y Gestión del Plan Nacional I+D+i and European Social Fund, Grant/Award Number: CGL2011-25053, CGL2014-54671-P, CGL2016-80742-R and CGL2017-83433-P; European Social Fund Abstract Understanding the processes that shape neutral and adaptive genomic variation is a fundamental step to determine the demographic and evolutionary dynamics of pest species. Here, we use genomic data obtained via restriction site-associated DNA sequencing to investigate the genetic structure of Moroccan locust (Dociostaurus maroccanus) populations from the westernmost portion of the species distribution (Iberian Peninsula and Canary Islands), infer demographic trends, and determine the role of neutral versus selective processes in shaping spatial patterns of genomic variation in this pest species of great economic importance. Our analyses showed that Iberian populations are characterized by high gene flow, whereas the highly isolated Canarian populations have experienced strong genetic drift and loss of genetic diversity. Historical demographic reconstructions revealed that all populations have passed through a substantial genetic bottleneck around the last glacial maximum (~21 ka BP) followed by a sharp demographic expansion at the onset of the Holocene, indicating increased effective population sizes during warm periods as expected from the thermophilic nature of the species. Genome scans and environmental association analyses identified several loci putatively under selection, suggesting that local adaptation processes in certain populations might not be impeded by widespread gene flow. Finally, all analyses showed few differences between outbreak and nonoutbreak populations. Integrated pest management practices should consider high population connectivity and the potential importance of local adaptation processes on population persistence. KEYWORDS ddRadSeq, demographic inference, environmental association analyses, genetic structure, local adaptation, pest species 2 | GONZÁLEZ-SERNA Et AL. 1 | INTRODUCTION The genetic makeup of a species is influenced by complex interactions between neutral and selective forces, life-history characteristics, and contemporary and past environmental conditions that collectively shape the evolutionary and demographic fate of its populations (Bernatchez et al., 2018). From a neutral perspective, recently developed analytical methods allow inferring complex demographic processes from genomic data and obtaining robust estimates of population size changes and gene flow at different time scales (Liu & Fu, 2015; Miles et al., 2017; Sherpa et al., 2018). Historical demographic reconstructions, linking population size changes with past environmental fluctuations, can help to predict future demographic trends of species under hypothetical climate change scenarios (Brown et al., 2016; Espindola et al., 2012; Fordham, Brook, Moritz, & NoguésBravo, 2014), whereas contemporary estimates of population genetic connectivity and spatially explicit landscape genetic analyses are useful to obtain baseline information on dispersal rates (Crossley, Rondon, & Schoville, 2019a) and identify corridors for gene flow (Crossley, Rondon, & Schoville, 2019b; Karsten, Addison, Jansen van Vuuren, & Terblanche, 2016; Qin et al., 2018; Venkatesan & Rasgon, 2010; Zepeda-Paulo et al., 2010). From an adaptive perspective, new molecular tools bring the opportunity to infer selection at specific loci displaying patterns of variation contrasting to those characterizing genomic regions only affected by neutral processes such as mutation, genetic drift, migration, and demographic changes (Berdan, Mazzoni, Waurick, Roehr, & Mayer, 2015; Luikart, England, Tallmon, Jordan, & Taberlet, 2003). Identifying loci that exhibit a significant reduction of within-population genetic variability and higher divergence across populations consistent with disruptive selection (Vasemägi & Primmer, 2005), coupled with environment association analyses (EAA) linking variation in allele frequencies with ecological gradients (Rellstab, Gugerli, Eckert, Hancock, & Holderegger, 2015), can help to determine local adaptation processes in response to specific selective forces such as those imposed by climate (Crossley, Chen, Groves, & Schoville, 2017; Dowle et al., 2017; Dudaniec, Yong, Lancaster, Svensson, & Hansson, 2018; Guo, Li, & Merilä, 2016), pesticide application (Crossley et al., 2017; Gassmann, Onstad, & Pittendrigh, 2009; Leftwich, Bolton, & Chapman, 2016), or host plant use (Gassmann et al., 2009; Simon et al., 2015; Soria-Carrasco et al., 2014). Gene flow is generally accepted to constrain local adaptation whereas strong divergent selection is expected to prevent interpopulation realized gene flow (Lenormand, 2002). Thus, joined inference of both neutral and selective phenomena can provide a more comprehensive understanding on the relative role of dispersal and local adaptation processes on structuring genetic variation (Dudaniec et al., 2018; Guo et al., 2016). Pest species, either invasive or native, are responsible of considerable economic losses worldwide and, as a result, public, private, and nonprofit organizations annually invest vast amounts of resources to prevent or mitigate their negative impacts (Enserink, 2004; Skaf, Popov, Roffey, Scorer, & Hewitt, 1990). In turn, such management practices often have undesirable side effects on wildlife, natural ecosystems, and human health (e.g., Baker & Wilkinson, 1990; Carson, 2002). For these reasons, understanding the population dynamics, dispersal routes, demographic history, and idiosyncratic evolutionary processes of pest species is a fundamental step to predict their future impacts and develop informed, less pernicious, and more targeted management practices (Abrol, 2014; Lankau, Jørgensen, Harris, & Sih, 2011). Population genetic approaches have proven useful to address several of the abovementioned aspects, and their potential has exponentially grown in the last years by the generalization of high-throughput sequencing techniques and the possibility of inferring both neutral and adaptive evolutionary processes at an unprecedented resolution (Crossley et al., 2017; Wang et al., 2014). The application of genomic tools is particularly important if we consider that most pest species often show large effective population sizes, high dispersal potential, shallow genetic differentiation, and fluctuating and complex demographic dynamics that are difficult to study using traditional capture–mark– recapture approaches or standard genetic methods (Bekkevold, Gross, Arula, Helyar, & Ojaveer, 2016; Chapuis et al., 2008, 2011; Ibrahim, Sourrouille, & Hewitt, 2000; Kirk, Dorn, & Mazzi, 2013). Locusts are a paradigmatic case of pest species with cyclical outbreaks that cause considerable agricultural losses and remission periods during which local populations either disappear or persist at very low densities (Chapuis et al., 2009; Enserink, 2004; Latchininsky, 1998; Skaf et al., 1990). The prevalence of outbreak and solitary phases varies geographically, with populations from some areas recurrently becoming agricultural pests while those from nonoutbreaking regions often occur at low numbers forming harmless populations (Chapuis et al., 2009; Latchininsky, 1998). An intriguing example of this demographic stochasticity is the case of the mysterious Rocky Mountain locust Melanoplus spretus, a devastating pest species endemic from North American prairies during the 19th century that went inexplicably extinct within a 30-year interval (Chapco & Litzenberger, 2004; Lockwood, 2004). Extreme demographic oscillations charactering locust populations are expected to have a considerable impact on genetic variation and the potential of species to respond to selection and evolve local adaptations. On the one hand, population crashes in the transition phase from gregarious to solitary forms are likely to leave genetic signatures of demographic bottlenecks that are probably ephemeral and blurred by genetic admixture after population expansions during outbreak periods (Chapuis et al., 2008; Ibrahim, 2001; Ibrahim et al., 2000). On the other hand, local adaptation processes in response to spatially varying selective pressures could be impeded by high gene flow (Babin, Gagnaire, Pavey, & Bernatchez, 2017; Lenormand, 2002; Pujolar et al., 2014) or restricted to isolated populations that are not swamped by gene flow from outbreaking populations (Chapuis et al., 2014). Previous single-locus and microsatellite-based studies on different locust species have found very shallow patterns of genetic structure at local/regional scales (e.g., Chapuis et al., 2011), no genetic differentiation between gregarious and solitary phase populations (Chapuis et al., 2008), higher levels of gene flow among outbreaking populations than among recession or nonoutbreak | 3 GONZÁLEZ-SERNA Et AL. populations (Chapuis et al., 2009; Ibrahim, 2001; Ibrahim et al., 2000), and a little impact of recession periods on genetic diversity (Chapco & Litzenberger, 2004; Chapuis et al., 2014; Ibrahim et al., 2000). However, with the exception of the recent sequencing and annotation of the locust Locusta migratoria genome (Wang et al., 2014) and the identification of some genes linked to physiology, phase change, and dispersal capacity (Bakkali & Martín-Blázquez, 2018; Ernst et al., 2015; Martín-Blázquez, Chen, Kang, & Bakkali, 2017; Wang et al., 2014), high-resolution genomic data have not been yet employed to determine fine-spatial scale patterns of genetic structure, perform robust demographic inferences in outbreak and nonoutbreak populations, and assess the potential role of selective processes on shaping spatial patterns of genetic variation in these organisms of great economic importance (e.g., Crossley et al., 2017). The Moroccan locust, Dociostaurus maroccanus (Thunberg, 1815; Figure 1), is a xerophilous species distributed in most of the Western Palearctic, from the Canary Islands to south Kazakhstan (Cigliano, Braun, Eades, & Otte, 2019; Latchininsky, 1998). The species is characterized by its broad polyphagia, extreme voracity, enormous fecundity, extraordinarily fluctuating populating sizes, and high capability to migrate (del Cañizo & Moreno, 1949; el Ghadraoui, Petit, Picaud, & Yamani, 2002; Latchininsky, 1998; Uvarov, 1977). Its distribution is discontinuous and consists of fragmented nonoutbreaking populations in some areas and permanent foci of outbreak populations that cyclically become devastating agricultural pests (del Cañizo & Moreno, 1949; Latchininsky, 1998). It has been reported that Moroccan locusts can move distances of 70–100 km during their entire lifetime (rarely up to 200 km; Latchininsky, 1998), making possible the exchange of individuals between distant populations during swarming phases (Latchininsky, 1998). The species is considered a major agricultural pest of high economic importance, damaging pastures, and a wide variety of crops during outbreaks, which requires extensive control operations and chemical interventions with a tremendous cost year after year in affected countries (e.g., Arias-Giralda, Jiménez-Viñuelas, & Pérez-Romero, 1997; Guerrero et al., 2019; Latchininsky, 2013). Here, we focus on populations of the Moroccan locust from the Iberian Peninsula and the Canary Islands, which represent the west margin of the species’ distribution (Latchininsky, 2013). In the Iberian Peninsula, there are three main vast regions that are foci of population outbreaks and have suffered considerable damages to pastures and crops for centuries: Monegros (Aragón), La Serena (Extremadura) and Valle de Alcudia (Castilla-La Mancha) (AlberolaRoma, 2012; Arias-Giralda, Morales-Agacino, Cobos-Suárez, SopeñaMañas, & Martín-Bernal, 1993). These regions are characterized by considerable cattle overgrazing, a factor that has been related to irruptive population growth in the Moroccan locust (Latchininsky, 1998; Louveaux, Mouhim, Roux, Gillon, & Barral, 1996). Besides, there are other regions from the Iberian Peninsula that represent historically outbreak areas that nowadays only sustain small populations or where the species has traditionally occurred at very low densities in isolated pockets of suitable habitat (Aragón, Coca-Abia, Llorente, & Lobo, 2013; Latchininsky, 1998). The small size of some formerly outbreak populations has been hypothesized to be related with the expansion of agriculture and certain plowing techniques, the destruction and fragmentation of suitable breeding areas linked to land use changes, and the massive application of pesticides to control locusts’ populations (Latchininsky, 1998). On this respect, several Moroccan locust populations from outbreak areas have been considerably reduced by human intervention to the point that many of them have almost disappeared in the last decades (Aragón et al., 2013; Latchininsky, 1998). However, no study has been performed so far to understand the degree of genetic connectivity among populations of the Moroccan locust at regional scales, infer its past demographic history, or determine the potential role of local adaptation processes, information that might help to shed light on key aspects of the ecology, distribution, and evolutionary dynamics of this economically important species. In this study, we use genomic data obtained via restriction site-associated DNA sequencing (ddRADseq) to investigate the relative role of neutral versus selective processes on shaping genetic variation in the Moroccan locust, determine spatial patterns of genetic diversity and structure in solitary and outbreaking populations from the westernmost portion of the distribution of the species, and, ultimately, infer its past demographic history. First, we performed genome scans and environmental association analyses to identify putative loci under selection and evaluate the potential importance of local adaptation processes on shaping genetic differentiation at non-neutral genomic regions. Second, we calculated different estimates of genetic diversity and performed a comprehensive suite of analyses of genetic structure to test the hypothesis of lower levels of genetic variation and increased genetic differentiation in solitary than outbreak populations FIGURE 1 Adult male (top) and nymph (bottom) of Moroccan locust (Dociostaurus maroccanus) from Valle de Alcudia (Ciudad Real, Spain). Photograph by Pedro J. Cordero 4 | GONZÁLEZ-SERNA Et AL. (Chapuis et al., 2014). Finally, we used genomic data to infer the past demography of the studied populations. Specifically, we predict genomic signatures of recent demographic declines in solitary populations and formerly outbreaking populations that have undergone remarkable retreats during the last decades (Chapuis et al., 2014; Ibrahim et al., 2000) and, given the thermophilous character of the Moroccan locust, we also expect that the species has experienced historical bottlenecks during Pleistocene glacial periods and population expansions in interglacials (Meco et al., 2011). 2 | MATERIALS AND METHODS 2.1 | Population sampling Between May and July 2011–2016, we prospected adequate habitats for the Moroccan locust (Dociostaurus maroccanus) (i.e., grazed grasslands, natural sparse vegetation, arid or semidesert steppes, and abandoned agricultural fields; Latchininsky, 1998; Latchininsky, 2013) in the Iberian Peninsula and the Canary Islands. We sampled a total of 21 localities representative of both outbreak and nonoutbreak populations (Figure 2; Table 1), a status defined according with our own field observations during sampling (i.e., densities of more than ~20 individuals/m2 were considered as outbreak populations) and corroborated with information provided by regional government authorities implementing pest management programs. Information about the outbreak or nonoutbreak status of the different sampling populations is presented in Table 1. Adult individuals were collected via sweep netting in an area not higher than 500 m2 around each sampling locality. Fresh whole adult specimens were placed in vials with 2–5 ml ethanol 96% and stored at −20°C until needed for genomic analyses. For this study, we analyzed a total of 5–8 adult individuals per locality (Table 1). 2.2 | Genomic library preparation We used NucleoSpin Tissue kits (Macherey-Nagel, Durën, Germany) to extract and purify genomic DNA from the hind femur of each individual. Genomic DNA from 141 individuals of the 21 sampling FIGURE 2 Geographic location of the studied populations of Moroccan locust in (a) the Iberian Peninsula and (c) the Canary Islands. Black dots indicate those populations analyzed with Stairway plot (n = 8 individuals) and white squares the rest of the populations (n < 8 individuals). Panels (b) and (c) present the inferred demographic profiles for Iberian and Canarian populations, respectively. Lines show the median estimate of effective population size (Ne) over time, assuming a mutation rate of 2.8 × 10–9 and 1-year generation time. Populations with an asterisk indicate pest outbreaks during the sampling year. Population codes are described in Table 1. | 5 GONZÁLEZ-SERNA Et AL. TABLE 1 Geographical location of the studied populations of Moroccan locust, population codes, sampling year, number of individuals per population (n), and scores for the three environmental principal components used in LFMM analyses (PC1, PC2, and PC3). Populations with an asterisk indicate pest outbreaks during the sampling year. Average values of genetic statistics across neutral loci are presented for major allele frequency (P), observed (HO) and expected (HE) heterozygosity, nucleotide diversity (π), and the Wright's inbreeding coefficient (FIS) for all positions (polymorphic and nonpolymorphic) Locality (Province) Code Year nLatitude Longitude PC1 PC2 PC3 P HOHEπFIS Trabanca (Salamanca) TRAB 2015 841.24275 −6.40222 −0.447 −0.283 1.282 0.9994 0.0007 0.0009 0.0010 0.0007 Sando (Salamanca)* SAND 2015 840.96961 −6.10954 −0.573 −0.279 0.630 0.9994 0.0007 0.0009 0.0010 0.0007 Salamanca (Salamanca) SALA 2015 540.93764 −5.66806 −0.608 −0.335 −0.726 0.9994 0.0007 0.0009 0.0010 0.0006 Alhama de Aragón (Zaragoza)* ALHA 2015 841.34571 −1.92277 −1.613 −0.820 −0.939 0.9994 0.0007 0.0009 0.0010 0.0007 Los Llanos de Cáceres (Cáceres) CACE 2015 839.53310 −6.32455 1.072 0.887 0.094 0.9994 0.0007 0.0009 0.0010 0.0008 Puerto de la Berzocana (Cáceres) BERZ 2012 539.44998 −5.42609 −0.349 0.401 1.006 0.9994 0.0007 0.0008 0.0009 0.0005 Castuera (Cáceres)* CAST 2015 838.74827 −5.54450 1.036 1.186 0.466 0.9994 0.0007 0.0009 0.0010 0.0007 Cañada del Hoyo (Cuenca) HOYO 2015 839.93709 −1.99317 −1.787 −0.467 0.010 0.9994 0.0007 0.0009 0.0010 0.0007 Laguna de Tirez (Toledo) TIRE 2015 739.54259 −3.35907 0.144 0.660 −1.063 0.9994 0.0007 0.0009 0.0010 0.0007 Raña Cornicabra (Ciudad Real)* CORN 2011 539.13134 −4.73814 0.550 0.924 −0.471 0.9994 0.0007 0.0009 0.0010 0.0005 Valle de Alcudia (Ciudad Real)* ALCU 2015 738.58958 −4.35354 0.560 1.059 0.502 0.9994 0.0008 0.0009 0.0010 0.0006 Belalcázar (Córdoba)* BELA 2015 538.67942 −5.16172 1.125 1.311 0.342 0.9994 0.0008 0.0009 0.0010 0.0005 El Bonillo (Albacete) BONI 2015 738.93758 −2.56544 −0.528 0.209 −0.507 0.9994 0.0007 0.0009 0.0010 0.0008 La Felipa (Albacete) FELI 2015 539.03552 −1.64832 −0.299 −0.185 −1.642 0.9994 0.0007 0.0009 0.0010 0.0005 Santa Elena (Jaén) SANT 2015 538.33753 −3.51923 0.594 0.935 0.098 0.9994 0.0008 0.0009 0.0010 0.0005 Santiago de la Espada (Jaén) ESPA 2011 538.18128 −2.66677 −1.238 −0.058 0.920 0.9994 0.0007 0.0009 0.0010 0.0005 Jumilla (Murcia)* JUMI 2015 538.49016 −1.24497 0.228 −0.237 −1.815 0.9994 0.0007 0.0009 0.0010 0.0005 Caravaca de la Cruz (Murcia)* CARA 2015 838.10283 −2.01999 −0.208 −0.053 −0.895 0.9994 0.0007 0.0009 0.0010 0.0007 Alpujarra (Granada)* ALPU 2015 836.96846 −3.21425 −0.986 −0.101 1.811 0.9994 0.0007 0.0009 0.0010 0.0007 Tenerife (Santa Cruz de Tenerife) TENE 2016 828.41569 −16.40991 1.343 −2.257 1.449 0.9995 0.0006 0.0007 0.0008 0.0004 El Hierro (Santa Cruz de Tenerife)* HIER 2016 827.77303 −17.95796 1.984 −2.497 −0.555 0.9995 0.0006 0.0007 0.0008 0.0004 6 | GONZÁLEZ-SERNA Et AL. localities (Table 1) was processed into three genomic libraries (47 individuals per library) using the double-digestion restriction siteassociated DNA sequencing procedure (ddRADseq) described in Peterson, Weber, Kay, Fisher, and Hoekstra (2012). In brief, DNA was doubly digested with the restriction enzymes MseI and EcoR1 (New England Biolabs, Ipswich, MA, USA) and Illumina adaptors including unique 7-bp barcodes were ligated to the digested fragments. Ligation products from individuals assigned to each library were pooled, size-selected between 475 and 580 bp with a Pippin Prep (Sage Science, Beverly, MA, USA) machine, and amplified by PCR with 12 cycles using the iProofTM High-Fidelity DNA Polymerase (BIO-RAD, Hercules, CA, USA). Each library was sequenced in single-read 150-bp lane on an Illumina HiSeq2500 platform at The Centre for Applied Genomics (SickKids, Toronto, ON, Canada). 2.3 | Genomic data analyses We used the different programs distributed as part of the StackS v. 1.35 pipeline (process_radtags, ustacks, cstacks, sstacks, and populations) to assemble our sequences into de novo loci and call genotypes (Catchen, Amores, Hohenlohe, Cresko, & Postlethwait, 2011; Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013a; Hohenlohe et al., 2010). For the different filtering and assembling steps with StackS, we used the default parameters recommended by the authors (Catchen et al., 2011; Catchen, Hohenlohe, et al., 2013; Hohenlohe et al., 2010). Reads were demultiplexed and filtered for overall quality using the program process_radtags, retaining reads with a Phred score > 10 (using a sliding window of 15%), no adaptor contamination, and that had an unambiguous barcode and restriction cut site. Raw reads were screened for quality with FaStqc v. 0.11.5 (http://www.bioin forma tics.babra ham.ac.uk/proje cts/fastqc), and all sequences were trimmed to 129-bp using Seqtk (https://github. com/lh3/seqtk) in order to remove low-quality reads near the 3´ ends. Filtered reads of each individual were assembled de novo into putative loci with the ustacks program. The minimum stack depth (m) was set to three, and we allowed a maximum distance of two nucleotide mismatches (M) to group reads into a “stack.” We used the “removal” (r) and “deleveraging” (d) algorithms to eliminate highly repetitive stacks and resolve over-merged loci, respectively. Single nucleotide polymorphisms (SNPs) were identified at each locus, and genotypes were called using a multinomial-based likelihood model that accounts for sequencing errors, with the upper bound of the error rate (ε) set to 0.2 (Catchen et al., 2011; Catchen, Hohenlohe, et al., 2013; Hohenlohe et al., 2010). A catalogue of loci was built using the cstacks program, with loci recognized as homologous across individuals if the number of nucleotide mismatches between consensus sequences (n) was ≤2. Each individual data were matched against this catalogue using sstacks program, and output files were exported in different formats for subsequent analyses using the program populations. For all downstream analyses, we exported only the first SNP per RAD locus (option write_single_snp) and retained loci with a minimum stack depth ≥ 5 (m = 5), that were sequenced in at least 50% of the individuals of each population (parameter r = 0.5), represented in ~66% of populations, and with a minimum minor allele frequency (MAF) ≥ 0.01 to reduce the number of false polymorphic loci due to sequencing errors. As demonstrated in previous studies, the choice of different filtering thresholds affecting the proportion of missing data (p and r in program populations) had little impact on the obtained inferences (e.g., González-Serna, Cordero, & Ortego, 2018; González-Serna, Cordero, & Ortego, 2019; Ortego, Gugger, & Sork, 2018). The resulting files were used for subsequent analyses or converted into other formats using the program pGDSpiDer v.2.1.0.3 (Lischer & Excoffier, 2012). 2.4 | Outlier loci detection and environmental association analyses In a first step, we screened for loci not conforming to neutral expectations using two outlier detection approaches: the coalescent-based FDIST method from arlequin (Excoffier & Lischer, 2010) and the Bayesian approach implemented in BayeScan v.2.1 (Foll & Gaggiotti, 2008). The FDIST method in arlequin was run in two different ways, considering both the nonhierarchical (Beaumont & Nichols, 1996) and hierarchical (Excoffier, Hofer, & Foll, 2009) island models. The nonhierarchical island model was run using 200,000 simulations, 100 demes, and expected heterozygosity ranging from 0 to 1. The hierarchical island model was run grouping populations according to their geographical origin and the results obtained from genetic structure analyses (i.e., Iberian versus Canarian populations; see Results section), using the same settings as the nonhierarchical model, and considering three simulated groups (i.e., the number of defined population groups plus one, as recommended in Excoffier and Lischer (2010). P-values were corrected for multiple testing using the p.adjust function in r (R Core Team, 2018) and loci significantly outside the neutral distribution at a false discovery rate (FDR) of 5% (i.e., q < 0.05) were considered as outliers. BayeScan analyses were run under default settings (thinning interval size of 10; 20 pilot runs of 5,000 iterations; burn-in length of 50,000 iterations), except for an increase of outputted iterations to 10,000. We used 10 (default) prior odds and adopted the same FDR to identify candidate loci putatively under selection as in arlequin analyses (FDR of 5%, q < 0.05). BayeScan and arlequin analyses were run independently for two different genomic datasets, one considering all populations and another one only considering populations from the Iberian Peninsula (e.g., Guo et al., 2016). In a second step, we performed environmental association analyses (EAA) using Latent Factor Mixed Models (LFMM) implemented in the r v.3.3.3 (R Core Team, 2018) package lea (Frichot & François, 2015). This approach uses a stochastic Monte Carlo Markov Chain algorithm to test for associations between environmental/ecological variables and allele frequencies while simultaneously controlling for background levels of population structure (Frichot, Schoville, Bouchard, & François, 2013). As environmental information, we used the 19 bioclimatic variables from the worlDclim dataset interpolated to 30-arcsec resolution (~1 km2 | 7 GONZÁLEZ-SERNA Et AL. cell size) (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005). These variables summarize information about temperature and precipitation and have been commonly used in exploratory analyses aimed to test the broad hypothesis of environment-driven selection (e.g., François, Martins, Caye, & Schoville, 2016; Frichot & François, 2015; de Kort et al., 2014; Rellstab et al., 2015). We extracted the values for these variables from all adjacent cells around each sampling location (i.e., ~9 km2) using bilinear interpolations as implemented in arcGiS 10.3 (ESRI, Redlands, CA, USA). To summarize and reduce strong redundancy among the 19 bioclimatic variables, we ran a principal component analysis (PCA) on all of them (e.g., François et al., 2016; Frichot & François, 2015; de Kort et al., 2014). The first three principal components (PCs) cumulatively accounted for 94.82% of the variance (PC1: 38.58%; PC2: 33.77%; PC3: 22.47%) and were retained for LFMM analyses (Table A1). The contribution of bioclimatic variables to each axis (i.e., the factor loadings for each PC) is presented in Table 1. The MCMC algorithm was implemented for each of the three PCs (i.e., PC1, PC2, and PC3), using 10,000 iterations, 5,000 as burning period, and 5 independent replicates of the analysis. As indicated above for BayeScan and nonhierarchical arlequin analyses, LFMM analyses were run independently for all populations and only considering Iberian populations. The number of latent factors included in the model as a covariate to control for demographic history were defined on the basis of Structure analyses (Pritchard, Stephens, & Donnelly, 2000; see Results section) and sparse non-negative matrix factorization (snmf) analyses implemented in the r package lea (Frichot, Mathieu, Trouillon, Bouchard, & Francois, 2014). We set the number of latent factors (K) at K = 2 for analyses including all populations and K = 1 for analyses focused on Iberian populations. The z-scores over the five replicates were combined and recalibrated using the genomic inflation factor (λ) (Frichot & François, 2015). Finally, we performed a FDR adjustment to control for multiple tests (FDR of 5%, q < 0.05) and identify putative loci under environmental selection for each of the three PCs summarizing bioclimatic variation. 2.5 | Genetic structure We employed three complementary approaches to exhaustively explore spatial patterns of genetic structure in our study system, including (a) principal component analyses (PCA) (Jombart, 2008); (b) classic Structure analyses considering and not considering prior population information (Hubisz, Falush, Stephens, & Pritchard, 2009; Pritchard et al., 2000); and (c) the recently developed spatial method implemented in the r program conStruct (Bradburd, Coop, & Ralph, 2018). In all cases, genetic structure was analyzed only considering putatively neutral loci. To this end, outlier loci detected by arlequin and BayeScan and SNPs identified by LFMM analyses as being putatively under environment-driven selection (a total of 9,346 loci; see Section 3) were conservatively excluded to create datasets only containing neutral loci (e.g., Brauer, Hammer, & Beheregaray, 2016; Ortego et al., 2018). This yielded neutral datasets of 40,179 SNPs for all populations, 42,114 SNPs for Iberian populations, and 33,998 SNPs for Canarian populations. 2.5.1 | Principal component analyses (PCA) In order to visualize the major axes of population genetic differentiation, we performed individual-based principal component analyses (PCA) using the r package aDeGenet (Jombart, 2008). Before running the PCA, we scaled and centered allele frequencies and replaced missing data with mean allele frequencies using the scaleGen function as recommended by Jombart (2008). PCAs were run using all neutral SNPs for the two main datasets (all populations and Iberian Peninsula). 2.5.2 | Structure analyses We inferred genetic structure at neutral loci using the Bayesian Markov chain Monte Carlo clustering method implemented in the program Structure v.2.3.3 (Falush, Stephens, & Pritchard, 2003; Hubisz et al., 2009; Pritchard et al., 2000). We conducted Structure analyses hierarchically, initially analyzing data from all populations jointly and, subsequently, running independent analyses for subsets of populations assigned to the same genetic cluster in the previous hierarchical level analysis (e.g., Coulon et al., 2008; González-Serna et al., 2019). To make analyses computationally tractable, we ran Structure using a single random subset of 10,000 unlinked neutral SNPs. According to previous studies, this number of loci is an order of magnitude higher than that necessary to obtain robust and reproducible results in Structure (e.g., Catchen, Bassham, et al., 2013:1,000 loci; González-Serna et al., 2019:1,250 loci). We ran Structure using 200,000 MCMC iterations after a burn-in step of 100,000 iterations, assuming correlated allele frequencies and admixture, and both considering and not considering prior population information (Hubisz et al., 2009). We performed 15 independent runs for each value of K. In order to ensure analysis convergence, we only retained the ten runs having the highest likelihood for each value of K (e.g., Yannic et al., 2018) and checked that all retained replicates reached a similar solution in terms of individual's probabilities of assignment to each genetic cluster (qvalues; Gilbert et al., 2012). As recommended, we used two statistics to identify the most likely number of genetic clusters (K) (Gilbert et al., 2012; Janes et al., 2017): log probabilities of Pr(X|K) (Pritchard et al., 2000) and the ΔK method (Evanno, Regnaut, & Goudet, 2005). These statistics were calculated as implemented in Structure HarveSter (Earl & vonHoldt, 2012). Finally, we used clumpp v. 1.1.2 and the Greedy algorithm to align multiple runs of Structure for the same K value (Jakobsson & Rosenberg, 2007) and DiStruct v. 1.1 (Rosenberg, 2004) to visualize as bar plots the individual's probabilities of membership to each inferred genetic cluster. 8 | GONZÁLEZ-SERNA Et AL. 2.5.3 | conStruct analyses We used the spatial model implemented in r package conStruct v. 1.0.3 to infer patterns of genetic structure across the Iberian populations and determine whether genetic differentiation is a consequence of continuous (i.e., isolation-by-distance) or discrete (e.g., separation by geographic barriers, etc.) processes (Bradburd et al., 2018). Given that this approach is highly sensitive to missing data, we analyzed a database with no missing data (2,904 unlinked SNPs) as recommended by the authors (Bradburd et al., 2018). We ran conStruct analyses with 5,000 iterations and visually checked for convergence using trace plots (Bradburd et al., 2018). We used a fivefold cross-validation approach to examine predictive accuracy across the range of tested K values and determine the best-fit number of genetic clusters (Bradburd et al., 2018). As done for Structure, we plotted individual coancestry coefficients for the different K values using DiStruct. 2.6 | Geographical and environmental drivers of genetic differentiation We tested for the presence of isolation-by-distance (IBD) and/or isolation-by-environment (IBE) patterns of genetic structure by analyzing the association between genetic differentiation (FST) and geographic and environmental distances among Iberian populations (Sexton, Hangartner, & Hoffmann, 2014; Shafer & Wolf, 2013; Wang, 2013). Genetic differentiation (FST) between all pairs of populations was calculated for the subset of neutral loci (see previous section) using the program populations from StackS (Table A2). Geographic distance between each pair of Iberian populations was calculated using GeoGrapHic DiStance matrix Generator v.1.2.3 (Ersts, 2018). Environmental distances were calculated for each PC (PC1, PC2, and PC3) obtained from a PCA on the 19 bioclimatic variables (see LFMM analyses above for details on PCA) using the “dist” function in r 3.3.3 (R Core Team, 2018). Genetic differentiation [FST/(1 − FST)] was tested against matrices of geographical (log10 transformed) and environmental distances using multiple matrix regressions with randomization (MMRR) as implemented in the “MMRR” function (Wang, 2013) in r 3.3.3 (R Core Team, 2018). Given that geographical/environmental distances are only expected to have positive effect on the degree of genetic differentiation between populations, we used one-tailed hypothesis tests for making statistical decisions regarding the null hypothesis of no effect of independent variables on genetic differentiation (Ruxton & Neuhäuser, 2010; e.g., Berkman, Nielsen, Roy, & Heist, 2013; Yannic et al., 2018). We selected final models following a backward procedure, initially fitting all explanatory terms and progressively removing nonsignificant variables until all retained variables were significant. The significance of the variables excluded from the model was tested again until no additional term reached significance (e.g., Ortego, Aguirre, Noguerales, & Cordero, 2015a). 2.7 | Genetic diversity and past demographic history We only employed neutral loci for calculating genetic diversity statistics and performing demographic inference analyses (Luikart et al., 2003). We used the program populations from StackS to calculate some genetic statistics, including nucleotide diversity (π), observed (HO) and expected (HE) heterozygosity, major allele frequency (P), and the Wright's inbreeding coefficient (FIS) (Catchen, Hohenlohe, et al., 2013). Standardized multilocus heterozygosity (sMLH) was calculated for each individual using the r package inBreeDr (Stoffel et al., 2016). sMLH is an individual-based metric defined as the total number of heterozygous loci in an individual divided by the sum of average observed heterozygosities in the population, over the subset of loci successfully typed in the focal individual (Coltman, Pilkington, Smith, & Pemberton, 1999). We used Stairway plot (Liu & Fu, 2015) to reconstruct the demographic history of the studied populations, a novel model-flexible method based on the site frequency spectrum (SFS) that does not require whole-genome sequence data or reference genome information to infer changes in effective population size (Ne) over time. These analyses were restricted to populations with eight genotyped individuals (see Figure 2 and Table 1), as the calculation of the SFS requires a downsampling procedure to remove missing data. These populations are representative of the distribution of the species across the study area (see Figure 2a). To compute the SFS for each population, we ran the program populations from StackS (Catchen, Hohenlohe, et al., 2013) in order to export the first SNP per RAD locus and retain loci with a minimum stack depth ≥ 5 (m = 5) and that were represented in at least 50% of the individuals of the focal population (r = 0.5). To remove all missing data for the calculation of the SFS and minimize errors with allele frequency estimates, each population was down-sampled to 6 individuals using a custom Python script written by Qixin He and available on Dryad (Papadopoulou & Knowles, 2015). We ran Stairway plot for each population fitting a flexible multi-epoch demographic model, assuming the mutation rate per site per generation of 2.8 × 10–9 estimated for Drosophila melanogaster (Keightley, Ness, Halligan, & Haddrill, 2014), a one-year generation time (Latchininsky, 1998), four different number of random breakpoints [(nseq-2)/4, (nseq-2)/2, (nseq-2)*3/4, and nseq-2], and 200 bootstrap replicates to estimate 95% confidence intervals. 3 | RESULTS 3.1 | Genomic data analyses Illumina sequencing of ddRAD libraries generated > 358 millions of reads in total after first quality filtering using the program process_ radtags. The number of reads per individual before and after different quality filtering steps is shown in Figure A1. Only one sample from population BONI was excluded for subsequent analyses due to low number of reads (Figure A1). The dataset obtained with StackS for all populations contained a total of 49,373 unlinked SNPs. | 9 GONZÁLEZ-SERNA Et AL. 3.2 | Outlier loci detection and environmental association analyses For the dataset including all populations, we identified 318 (0.64%) outlier loci using arlequin (FDIST method) and 503 (1.02%) using BayeScan (Figures 3a,b). When the analyses were restricted to Iberian populations, we identified 93 (0.18%) outlier loci using arlequin and 270 (0.52%) using Bayescan (Figure 3a,b). Several SNPs were commonly identified as outliers by both arlequin and BayeScan analyses (all populations: 74 SNPs, 0.15%, Figure 3a; Iberian populations: 35 SNPs, 0.07%, Figure 3b). Most outlier loci were identified as being putatively under divergent selection in analyses based on both the datasets including all populations (arlequin: n = 275, 86.47%; BayeScan: n = 467, 92.84%) and the one restricted to Iberian populations (arlequin: n = 76, 81.72%; BayeScan: n = 269, 99.63%). Environmental association analyses in LFMM showed that a high number of loci were significantly associated with environmental variation (Figure 3c,d). For the dataset including all populations, LFMM detected 7,710 unique loci (15.62%) associated with at least one PC of environmental variation (PC1: 2,932 loci; PC2: 4,051 loci; PC3: 2,797 loci; Figure A2) and 196 of them (0.40%) showed significant associations with the three PCs (Figure 3c). Similarly, when the analyses were restricted to Iberian populations, LFMM detected 6,104 unique loci (11.87%) associated with at least one PC of environmental variation (PC1: 2,713 loci; PC2: 3,081 loci; PC3: 2,982 loci; Figure A2) of which 363 (0.71%) were shared across all PCs (Figure 3d). Only 42 loci (0.08%) for analyses based on all populations (Figure 3a) and 23 loci (0.04%) for analyses focused on Iberian populations (Figure 3b) showed associations with environmental variation in LFMM and were also identified as FST outliers by arlequin and BayeScan analyses. 3.3 | Genetic structure The PCA based on all populations showed a clear separation of Iberian and Canarian populations (Figure 4a). The PCA for the subset of Iberian populations showed that HOYO was the most differentiated population and some individuals from two other populations (ALHA and ALPU) also tended to stand out from the rest (Figure 4b). These results are in good agreement with those obtained from Bayesian clustering analyses. Structure analyses for all populations and not considering prior population information identified K = 2 as the most likely clustering solution according to both the ΔK criterion and log probabilities of the data [LnPr (X|K)] (Figure A3a). These two clusters showed a very low degree of genetic admixture and split Canarian and Iberian populations (Figure 5c). Remarkably, genetic drift after divergence (parameter F in Structure; Pritchard et al., 2000) for the cluster corresponding to the Canary Islands (Fvalue = 0.419) was more than threefold the estimated for Iberian populations (F-value = 0.121). Structure analyses for K = 3 revealed further genetic structure and showed that one population from the Iberian Peninsula (BERZ) split from the rest of the mainland populations, albeit with a high degree of genetic admixture with the other Iberian populations (~10%–20%; Figure 5c). Structure analyses restricted to Iberian populations identified that the most likely number of clusters was K = 2 according to the ΔK criterion, but LnPr (X|K) steadily increased up to K = 4 (Figure A3c). For K = 2, all the Iberian populations and individuals presented the same proportion of ancestry to the two inferred clusters (~20/80), indicating that they represent ghost clusters with no biological significance (see Chen, Durand, Forbes, & François, 2007; Guillot, Estoup, Mortier, & Cosson, 2005; Tonzo, Papadopoulou, & Ortego, 2019). However, Structure analyses for K = 3 and K = 4 showed that BERZ and HOYO were assigned to two different genetic groups albeit with some degree of genetic admixture with the rest of Iberian populations (Figure 5c). The genetic clusters corresponding to these two populations had much higher estimates of genetic drift after divergence (BERZ: F-value = 0.263; HOYO: F-value = 0.084) than the one representing the rest of Iberian populations (F-value = 0.021). Finally, Structure analyses restricted to populations from the Canary Islands showed K = 2 as the most likely clustering solution according to both the ΔK criterion and LnPr (X|K) (Figure A3e). These two clusters separated HIER and TENE populations, which showed a very low degree of genetic admixture (Figure 5c). Structure analyses ran considering prior population information (Hubisz et al., 2009) yielded qualitatively similar results, albeit in this case BERZ and HOYO presented a very small probability of assignment to their specific clusters in analyses focused on Iberian populations for K = 3 and K = 4 (see Figures A3 and A4). Finally, spatial analyses in conStruct focused on Iberian populations also showed strong admixture and no clear pattern of genetic structure. The predictive accuracy of conStruct analyses sharply increased from K = 1 to K = 2 (Figure A5a). However, layers (i.e., genetic clusters) beyond K = 1 contributed relatively little to total covariance (Figure A5b). This was particularly remarkable for K = 2, in which the second layer contributed less than 3% to total covariance (Figure A5b). Accordingly, the inferred genetic clusters for K > 1 presented high genetic admixture with little geographical congruence and only BERZ for K = 2–3 and ALHA for K = 3 tended to be assigned to different genetic clusters (Figure A6). 3.4 | Geographical and environmental drivers of genetic differentiation Multiple matrix regression with randomization (MMRR) analyses showed that genetic differentiation [FST/(1 − FST)] among Iberian populations was not significantly associated with geographical or environmental (PC1, PC2, and PC3) distances (Table 2). 3.5 | Genetic diversity and past demographic history Population genetic statistics (P, HO, HE, π, FIS) calculated in StackS for all positions (polymorphic and nonpolymorphic) are presented in 16 | GONZÁLEZ-SERNA Et AL. Foll, M., & Gaggiotti, O. (2008). A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: A Bayesian perspective. Genetics, 180(2), 977–993. https://doi. org/10.1534/genet ics.108.092221 Fordham, D. A., Brook, B. W., Moritz, C., & Nogués-Bravo, D. (2014). Better forecasts of range dynamics using genetic data. Trends in Ecology & Evolution, 29(8), 436–443. https://doi.org/10.1016/j. tree.2014.05.007 François, O., Martins, H., Caye, K., & Schoville, S. D. (2016). Controlling false discoveries in genome scans for selection. Molecular Ecology, 25(2), 454–469. https://doi.org/10.1111/mec.13513 Frichot, E., & François, O. (2015). LEA: An R package for landscape and ecological association studies. Methods in Ecology and Evolution, 6(8), 925–929. https://doi.org/10.1111/2041-210X.12382 Frichot, E., Mathieu, F., Trouillon, T., Bouchard, G., & Francois, O. (2014). Fast and efficient estimation of individual ancestry coefficients. Genetics, 196(4), 973–983. https://doi.org/10.1534/genet ics.113.160572 Frichot, E., Schoville, S. D., Bouchard, G., & François, O. (2013). Testing for associations between loci and environmental gradients using Latent Factor Mixed Models. Molecular Biology and Evolution, 30(7), 1687–1699. https://doi.org/10.1093/molbe v/mst063 Gassmann, A. J., Onstad, D. W., & Pittendrigh, B. R. (2009). Evolutionary analysis of herbivorous insects in natural and agricultural environments. Pest Management Science, 65, 1174–1181. https://doi. org/10.1002/ps.1844 Gilbert, K. J., Andrew, R. L., Bock, D. G., Franklin, M. T., Kane, N. C., Moore, J. S., … Vines, T. H. (2012). Recommendations for utilizing and reporting population genetic analyses: The reproducibility of genetic clustering using the program Structure. Molecular Ecology, 21(20), 4925–4930. https://doi.org/10.1111/j.1365-294X.2012.05754.x González-Serna, M. J., Cordero, P. J., & Ortego, J. (2018). Using high-throughput sequencing to investigate the factors structuring genomic variation of a Mediterranean grasshopper of great conservation concern. Scientific Reports, 8(1), 13436. https://doi. org/10.1038/s4159 8-018-31775 -x González-Serna, M. J., Cordero, P. J., & Ortego, J. (2019). Spatiotemporally explicit demographic modelling supports a joint effect of historical barriers to dispersal and contemporary landscape composition on structuring genomic variation in a red-listed grasshopper. Molecular Ecology, 28(9), 2155–2172. https://doi. org/10.1111/mec.15086 Guerrero, A., Ramos, V. E., López, S., Álvarez, J. M., Domínguez, A., CocaAbia, M. M., … Quero, C. (2019). Enantioselective synthesis and activity of all diastereoisomers of (e)-phytal, a pheromone component of the Moroccan locust, Dociostaurus maroccanus. Journal of Agricultural and Food Chemistry, 67(1), 72–80. https://doi.org/10.1021/acs. jafc.8b06346 Guillot, G., Estoup, A., Mortier, F., & Cosson, J. F. (2005). A spatial statistical model for landscape genetics. Genetics, 170(3), 1261–1280. https://doi.org/10.1534/genet ics.104.033803 Guo, B., Li, Z., & Merilä, J. (2016). Population genomic evidence for adaptive differentiation in the Baltic sea herring. Molecular Ecology, 25(12), 2833–2852. https://doi.org/10.1111/mec.13657 Hewitt, G. M. (1999). Post-glacial re-colonization of European biota. Biological Journal of the Linnean Society, 68(1–2), 87–112. https://doi. org/10.1111/j.1095-8312.1999.tb011 60.x Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., & Jarvis, A. (2005). Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology, 25(15), 1965–1978. https:// doi.org/10.1002/joc.1276 Hoarau, G., Rijnsdorp, A. D., Van Der Veer, H. W., Stam, W. T., & Olsen, J. L. (2002). Population structure of plaice (Pleuronectes platessa L.) in northern Europe: Microsatellites revealed large-scale spatial and temporal homogeneity. Molecular Ecology, 11(7), 1165–1176. https:// doi.org/10.1046/j.1365-294X.2002.01515.x Hohenlohe, P. A., Bassham, S., Etter, P. D., Stiffler, N., Johnson, E. A., & Cresko, W. A. (2010). Population genomics of parallel adaptation in threespine stickleback using sequenced RAD tags. PLoS Genetics, 6(2), e1000862. https://doi.org/10.1371/journ al.pgen.1000862 Hubisz, M. J., Falush, D., Stephens, M., & Pritchard, J. K. (2009). Inferring weak population structure with the assistance of sample group information. Molecular Ecology Resources, 9(5), 1322–1332. https://doi. org/10.1111/j.1755-0998.2009.02591.x Ibrahim, K. M. (2001). Plague dynamics and population genetics of the desert locust: Can turnover during recession maintain population genetic structure? Molecular Ecology, 10(3), 581–591. https://doi. org/10.1046/j.1365-294x.2001.01212.x Ibrahim, K. M., Sourrouille, P., & Hewitt, G. M. (2000). Are recession populations of the desert locust (Schistocerca gregaria) remnants of past swarms? Molecular Ecology, 9(6), 783–791. https://doi. org/10.1046/j.1365-294x.2000.00932.x Jakobsson, M., & Rosenberg, N. A. (2007). clumpp: A cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics, 23(14), 1801–1806. https://doi.org/10.1093/bioin forma tics/ btm233 Janes, J. K., Miller, J. M., Dupuis, J. R., Malenfant, R. M., Gorrell, J. C., Cullingham, C. I., & Andrew, R. L. (2017). The K=2 conundrum. Molecular Ecology, 26(14), 3594–3602. https://doi.org/10.1111/ mec.14187 Jeffery, N. W., Bradbury, I. R., Stanley, R. R. E., Wringe, B. F., Van Wyngaarden, M., Lowen, J. B., … DiBacco, C. (2018). Genomewide evidence of environmentally mediated secondary contact of European green crab (Carcinus maenas) lineages in eastern North America. Evolutionary Applications, 11(6), 869–882. https://doi.org/10.1111/ eva.12601 Jombart, T. (2008). aDeGenet: A R package for the multivariate analysis of genetic markers. Bioinformatics, 24(11), 1403–1405. https://doi. org/10.1093/bioin forma tics/btn129 Karsten, M., Addison, P., Jansen van Vuuren, B., & Terblanche, J. S. (2016). Investigating population differentiation in a major African agricultural pest: Evidence from geometric morphometrics and connectivity suggests high invasion potential. Molecular Ecology, 25(13), 3019–3032. https://doi.org/10.1111/mec.13646 Keightley, P. D., Ness, R. W., Halligan, D. L., & Haddrill, P. R. (2014). Estimation of the spontaneous mutation rate per nucleotide site in a Drosophila melanogaster full-sib family. Genetics, 196(1), 313–320. https://doi.org/10.1534/genet ics.113.158758 Kirk, H., Dorn, S., & Mazzi, D. (2013). Molecular genetics and genomics generate new insights into invertebrate pest invasions. Evolutionary Applications, 6(5), 842–856. https://doi.org/10.1111/eva.12071 Lankau, R., Jørgensen, P. S., Harris, D. J., & Sih, A. (2011). Incorporating evolutionary principles into environmental management and policy. Evolutionary Applications, 4(2), 315–325. https://doi. org/10.1111/j.1752-4571.2010.00171.x Laporte, M., Pavey, S. A., Rougeux, C., Pierron, F., Lauzent, M., Budzinski, H., … Bernatchez, L. (2016). RAD sequencing reveals within-generation polygenic selection in response to anthropogenic organic and metal contamination in north Atlantic eels. Molecular Ecology, 25(1), 219–237. https://doi.org/10.1111/mec.13466 Latchininsky, A. V. (1998). Moroccan locust Dociostaurus maroccanus (Thunberg, 1815): A faunistic rarity or an important economic pest? Journal of Insect Conservation, 2(3), 167–178. https://doi. org/10.1023/a:10096 39628627 Latchininsky, A. V. (2013). Locusts and remote sensing: A review. Journal of Applied Remote Sensing, 7(1), 1–19. https://doi.org/10.1117/1. JRS.7.075099 | 17 GONZÁLEZ-SERNA Et AL. Leftwich, P. T., Bolton, M., & Chapman, T. (2016). Evolutionary biology and genetic techniques for insect control. Evolutionary Applications, 9(1), 212–230. https://doi.org/10.1111/eva.12280 Lenormand, T. (2002). Gene flow and the limits to natural selection. Trends in Ecology & Evolution, 17(4), 183–189. https://doi.org/10.1016/ S 0 1 6 9 - 5 3 4 7 ( 0 2 ) 0 2 4 9 7 - 7 Lischer, H. E., & Excoffier, L. (2012). pGDSpiDer: An automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics, 28(2), 298–299. https://doi.org/10.1093/bioin f o r m a t i c s / b t r 6 4 2 Liu, X., & Fu, Y.-X. (2015). Exploring population size changes using SNP frequency spectra. Nature Genetics, 47(5), 555–559. https://doi. org/10.1038/ng.3254 Lockwood, J. A. (2004). Locust: The devastating rise and mysterious disappearance of the insect that shaped the American frontier (p. 320). New York, NY: Basic Books. Louveaux, A., Mouhim, A., Roux, G., Gillon, Y., & Barral, H. (1996). Effect of pastoral activities upon locust populations in the Siroua Massif (Morocco). Revue D'écologie: La Terre Et La Vie, 51(2), 139–151. Luikart, G., England, P. R., Tallmon, D., Jordan, S., & Taberlet, P. (2003). The power and promise of population genomics: From genotyping to genome typing. Nature Reviews Genetics, 4, 981. https://doi. org/10.1038/nrg1226 Martín-Blázquez, R., Chen, B., Kang, L., & Bakkali, M. (2017). Evolution, expression and association of the chemosensory protein genes with the outbreak phase of the two main pest locusts. Scientific Reports, 7(1), 6653. https://doi.org/10.1038/s4159 8-017-07068 -0 Smith, J., & Haigh, J. (2007). The hitch-hiking effect of a favourable gene. Genetics Research, 89(5-6), 391–403. https://doi.org/10.1017/S0016 67230 8009579 Meco, J., Muhs, D. R., Fontugne, M., Ramos, A. J., Lomoschitz, A., & Patterson, D. (2011). Late Pliocene and quaternary Eurasian locust infestations in the Canary archipelago. Lethaia, 44(4), 440–454. https://doi.org/10.1111/j.1502-3931.2010.00255.x Meco, J., Petit-Maire, N., Ballester, J., Betancort, J. F., & Ramos, A. J. (2010). The Acridian plagues, a new Holocene and Pleistocene palaeoclimatic indicator. Global and Planetary Change, 72(4), 318–320. https://doi.org/10.1016/j.glopl acha.2010.01.007 Miles, A., Harding, N. J., Bottà, G., Clarkson, C. S., Antão, T., Kozak, K., … Kwiatkowski, D. P. (2017). Genetic diversity of the African malaria vector Anopheles gambiae. Nature, 552, 96. https://doi.org/10.1038/ natur e24995 Ortego, J., Aguirre, M. P., Noguerales, V., & Cordero, P. J. (2015a). Consequences of extensive habitat fragmentation in landscape-level patterns of genetic diversity and structure in the Mediterranean esparto grasshopper. Evolutionary Applications, 8(6), 621–632. https:// doi.org/10.1111/eva.12273 Ortego, J., Gugger, P. F., & Sork, V. L. (2015b). Climatically stable landscapes predict patterns of genetic structure and admixture in the Californian canyon live oak. Journal of Biogeography, 42(2), 328–338. https://doi.org/10.1111/jbi.12419 Ortego, J., Gugger, P. F., & Sork, V. L. (2018). Genomic data reveal cryptic lineage diversification and introgression in Californian golden cup oaks (section Protobalanus). New Phytologist, 218(2), 804–818. https://doi.org/10.1111/nph.14951 Papadopoulou, A., & Knowles, L. L. (2015). Species-specific responses to island connectivity cycles: Refined models for testing phylogeographic concordance across a Mediterranean Pleistocene aggregate island complex. Molecular Ecology, 24(16), 4252–4268. https://doi. org/10.1111/mec.13305 Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S., & Hoekstra, H. E. (2012). Double digest RADseq: An inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS ONE, 7(5), e37135. https://doi.org/10.1371/journ al.pone.0037135 Pritchard, J. K., Stephens, M., & Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics, 155(2), 945–959. Pujolar, J. M., Jacobsen, M. W., Als, T. D., Frydenberg, J., Munch, K., Jónsson, B., … Hansen, M. M. (2014). Genome-wide single-generation signatures of local selection in the panmictic European eel. Molecular Ecology, 23(10), 2514–2528. https://doi.org/10.1111/mec.12753 Qin, Y.-J., Krosch, M. N., Schutze, M. K., Zhang, Y., Wang, X.-X., Prabhakar, C. S., … Li, Z.-H. (2018). Population structure of a global agricultural invasive pest, Bactrocera dorsalis (Diptera: Tephritidae). Evolutionary Applications, 11(10), 1990–2003. https://doi.org/10.1111/eva.12701 Quesada-Moraga, E., & Santiago-Álvarez, C. (2000). Temperature related effects on embryonic development of the Mediterranean locust, Dociostaurus Maroccanus. Physiological Entomology, 25(2), 191–195. https://doi.org/10.1046/j.1365-3032.2000.00185.x R Core Team. (2018). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Retrieved from https://www.R-proje ct.org/. Rellstab, C., Gugerli, F., Eckert, A. J., Hancock, A. M., & Holderegger, R. (2015). A practical guide to environmental association analysis in landscape genomics. Molecular Ecology, 24(17), 4348–4370. https:// doi.org/10.1111/mec.13322 Rosenberg, N. A. (2004). DiStruct: A program for the graphical display of population structure. Molecular Ecology Notes, 4(1), 137–138. https:// doi.org/10.1046/j.1471-8286.2003.00566.x Ruxton, G. D., & Neuhauser, M. (2010). When should we use one-tailed hypothesis testing? Methods in Ecology and Evolution, 1(2), 114–117. https://doi.org/10.1111/j.2041-210X.2010.00014.x Ryynanen, H. J., Tonteri, A., Vasemagi, A., & Primmer, C. R. (2007). A comparison of biallelic markers and microsatellites for the estimation of population and conservation genetic parameters in Atlantic salmon (Salmo salar). Journal of Heredity, 98(7), 692–704. https://doi. org/10.1093/jhere d/esm093 Schrey, N. M., Schrey, A. W., Heist, E. J., & Reeve, J. D. (2008). Fine-scale genetic population structure of southern pine beetle (Coleoptera: Curculionidae) in Mississippi forests. Environmental Entomology, 37(1), 271–276. https://doi.org/10.1093/ee/37.1.271 Sexton, J. P., Hangartner, S. B., & Hoffmann, A. A. (2014). Genetic isolation by environment or distance: Which pattern of gene flow is most common? Evolution, 68(1), 1–15. https://doi.org/10.1111/evo.12258 Shafer, A. B., & Wolf, J. B. (2013). Widespread evidence for incipient ecological speciation: A meta-analysis of isolation-by-ecology. Ecology Letters, 16(7), 940–950. https://doi.org/10.1111/ele.12120 Sherpa, S., Rioux, D., Goindin, D., Fouque, F., Francois, O., & Despres, L. (2018). At the origin of a worldwide invasion: Unraveling the genetic makeup of the Caribbean bridgehead populations of the dengue vector Aedes aegypti. Genome Biology and Evolution, 10(1), 56–71. https:// doi.org/10.1093/gbe/evx267 Simon, J.-C., d'Alencon, E., Guy, E., Jacquin-Joly, E., Jaquiery, J., Nouhaud, P., … Streiff, R. (2015). Genomics of adaptation to host-plants in herbivorous insects. Briefings in Functional Genomics, 14(6), 413–423. https://doi.org/10.1093/bfgp/elv015 Skaf, R., Popov, G. B., Roffey, J., Scorer, R. S., & Hewitt, J. (1990). The desert locust: An international challenge [and discussion]. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 328(1251), 525–538. Snyder, C. W. (2016). Evolution of global temperature over the past two million years. Nature, 538, 226–228. https://doi.org/10.1038/natur e19798 Soria-Carrasco, V., Gompert, Z., Comeault, A. A., Farkas, T. E., Parchman, T. L., Johnston, J. S., … Nosil, P. (2014). Stick insect genomes reveal natural selection's role in parallel speciation. Science, 344(6185), 738–742. https://doi.org/10.1126/scien ce.1252136 Stoffel, M. A., Esser, M., Kardos, M., Humble, E., Nichols, H., David, P., & Hoffman, J. I. (2016). inBreeDr: An R package for the analysis of 18 | GONZÁLEZ-SERNA Et AL. inbreeding based on genetic markers. Methods in Ecology and Evolution, 7(11), 1331–1339. https://doi.org/10.1111/2041-210X.12588 Tonzo, V., Papadopoulou, A., & Ortego, J. (2019). Genomic data reveal deep genetic structure but no support for current taxonomic designation in a grasshopper species complex. Molecular Ecology, 28(17), 3869–3886. https://doi.org/10.1111/mec.15189 Uvarov, B. (1977). Grasshoppers and locusts: A handbook of general Acridology (vol. 2): Behaviour, ecology, biogeography, population dynamics (Vol. 2, p. 613). London, UK: Centre for Overseas Pest Research. Vasemägi, A., & Primmer, C. R. (2005). Challenges for identifying functionally important genetic variation: The promise of combining complementary research strategies. Molecular Ecology, 14(12), 3623– 3642. https://doi.org/10.1111/j.1365-294X.2005.02690.x Venkatesan, M., & Rasgon, J. L. (2010). Population genetic data suggest a role for mosquito-mediated dispersal of west Nile virus across the western United States. Molecular Ecology, 19(8), 1573–1584. https:// doi.org/10.1111/j.1365-294X.2010.04577.x Wang, I. J. (2013). Examining the full effects of landscape heterogeneity on spatial genetic variation: A multiple matrix regression approach for quantifying geographic and ecological isolation. Evolution, 67(12), 3403–3411. https://doi.org/10.1111/evo.12134 Wang, X., Fang, X., Yang, P., Jiang, X., Jiang, F., Zhao, D., … Kang, L. E. (2014). The locust genome provides insight into swarm formation and long-distance flight. Nature Communications, 5, 2957. https://doi. org/10.1038/ncomm s3957 Yannic, G., Ortego, J., Pellissier, L., Lecomte, N., Bernatchez, L., & Cote, S. D. (2018). Linking genetic and ecological differentiation in an ungulate with a circumpolar distribution. Ecography, 41(6), 922–937. https://doi.org/10.1111/ecog.02995 Zepeda-paulo, F. A., Simon, J.-C., Ramírez, C. C., Fuentes-contreras, E., Margaritopoulos, J. T., Wilson, A. C. C., … Figueroa, C. C. (2010). The invasion route for an insect pest species: The tobacco aphid in the new world. Molecular Ecology, 19(21), 4738–4752. https://doi. org/10.1111/j.1365-294X.2010.04857.x SUPPORTING INFORMATION Additional supporting information may be found online in the Supporting Information section. How to cite this article: González-Serna MJ, Cordero PJ, Ortego J. Insights into the neutral and adaptive processes shaping the spatial distribution of genomic variation in the economically important Moroccan locust (Dociostaurus maroccanus). Ecol Evol. 2020;00:1–18. https://doi. org/10.1002/ece3.6165