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

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.


| INTRODUC TI ON
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 . 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és-Bravo, 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-markrecapture approaches or standard genetic methods (Bekkevold, Gross, Arula, Helyar, & Ojaveer, 2016;Chapuis et al., 2008Chapuis et al., , 2011Ibrahim, 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 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).
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 -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 ge-  (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).

| 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/m 2 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 m 2 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).

| 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 F I G U R E 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 (N e ) 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.

| Genomic data analyses
We used the different programs distributed as part of the StackS  (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).

| 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).
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 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.

| 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

| 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).

| 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 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.

| 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.

| Geographical and environmental drivers of genetic differentiation
We tested for the presence of isolation-by-distance (

| 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 (H O ) and expected (H E ) heterozygosity, major allele frequency (P), and the Wright's inbreeding coefficient (F IS ) (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 (N e ) over time.
These analyses were restricted to populations with eight genotyped individuals (see Figure 2 and Table 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.

| Genetic structure
The PCA based on all populations showed a clear separation of  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).

| Geographical and environmental drivers of genetic differentiation
Multiple matrix regression with randomization (MMRR) analyses showed that genetic differentiation [F ST /(1 − F ST )] among Iberian populations was not significantly associated with geographical or environmental (PC1, PC2, and PC3) distances (Table 2).

| Genetic diversity and past demographic history
Population genetic statistics (P, H O , H E , π, F IS ) calculated in StackS for all positions (polymorphic and nonpolymorphic) are presented in Table 1. Population genetic diversity (π) and standardized multilocus heterozygosity (sMLH) for each individual are shown in Figure 5a.
All estimates of population genetic diversity were significantly in-  a subsequent expansion at the onset of the Holocene to reach contemporary effective population sizes of 220,000-350,000 diploid individuals (Figure 2c; Figure A7).

| D ISCUSS I ON
We employed a large single-nucleotide polymorphism dataset (~50,000 loci) to infer historical demographic trends and understand the neutral and selective processes shaping spatial patterns of genomic variation in the Moroccan locust, a grasshopper that has traditionally emerged as an devastating pest in overgrazed areas and caused extensive agricultural damage in Spain (e.g., Buj-Buj, 1992) and many other areas across its ample distribution range (Latchininsky,  Table 1. 1998). Our analyses, focused on the Iberian Peninsula and remote populations from the Canary Islands, indicated that the species is characterized by widespread gene flow and low levels of genetic differentiation at regional scales. They also showed little differences in past demography and levels of genetic diversity and differentiation between outbreak and nonoutbreak populations. Finally, genome scans revealed that a small fraction (~0.2%-1.0%) of the sequenced loci are putatively under divergent selection and might be involved in local adaptation processes in response to the different ecological and environmental gradients experienced by the species.  Table 1) and Structure analyses showed that genetic drift after divergence for the cluster corresponding to the Canary Islands was much higher than that estimated for Iberian populations (Pritchard et al., 2000). As expected, levels of genetic differentiation among populations of the Moroccan locust were similar, albeit in some cases sensibly higher, than microsatellite-based estimates obtained for other locusts and pest grasshoppers sampled at a similar spatial scale (Table   A3). However, it should be noted that the different kinds of markers employed (i.e., SNPs versus microsatellites) make difficult the direct comparison of our estimates of genetic differentiation with those obtained in previous studies and, thus, small differences in absolute values must be interpreted with extreme caution (DeFaveri, Viitaniemi, Leder, & Merila, 2013;Ryynanen, Tonteri, Vasemagi, & Primmer, 2007). Overall, our results are in agreement with previous studies on other locust (e.g., Chapuis et al., 2008;Chapuis et al., 2011) and pelagic marine species (e.g., Als et al., 2011;Hoarau, Rijnsdorp, Veer, Stam, & Olsen, 2002) finding very low levels of genetic differentiation and indicate occasional genetic drift and differentiation in nonoutbreaking, solitarious or phase transition populations persisting at low densities in isolated pockets of suitable habitats (Chapuis et al., 2014;Ibrahim et al., 2000).

| Past demographic history
Demographic reconstructions for the Moroccan locust using  (Table 1) and compatible with the less severe impact of Pleistocene climatic oscillations at lower latitudes (Fernández-Palacios et al., 2016;Hewitt, 1999;Snyder, 2016 (Meco et al., 2010(Meco et al., , 2011Snyder, 2016). The overall good correspondence between peaks of population size and warm periods inferred from both fossil records (Meco et al., 2010(Meco et al., , 2011 and genomic data (present study) can be explained by the thermo-

| Putative loci under selection and potential for local adaptation
Genome scans based on F ST outlier tests revealed that a small portion of the sampled genome (between 0.2% and 1.0%, depending on the dataset and method) is putatively under selection, whereas environmental association analyses (EAA) identified a much larger proportion of SNPs (12%-15%) with allele frequencies correlated with one or more environmental gradients. Beyond these differences in numbers, the specific loci identified by the two methods showed little overlap (Figure 3), which is in agreement with the results obtained by previous RADseq-based studies (e.g., Dudaniec et al., 2018;Guo et al., 2016 Jeffery et al., 2018). This might be explained by the very scarce genomic resources available for grasshoppers, with only one draft genome sequenced so far for a species (Locusta migratoria) (Wang et al., 2014) belonging to a different subfamily than the Moroccan locust (Cigliano et al., 2019). Also, the putative signals of selection in certain loci could be a consequence of genetic hitchhiking resulted from linkage disequilibrium between the identified outlier SNPs and nearby nonsequenced genes actually subjected to diversifying selection (Smith & Haigh, 1974). This is particularly relevant considering the extraordinarily large size charactering the genome of grasshoppers, which might have resulted in we have only sampled a relatively small representation of it (Camacho et al., 2015;Wang et al., 2014).

| Conclusions and prospects
Overall, our study shows for the first time that populations of the economically important Moroccan locust present a shallow genetic differentiation and little differences in past demography, putative signatures of selection, and contemporary levels of genetic diversity and structure between outbreak and nonoutbreak populations.
Outbreaks of Moroccan locust are usually linked to considerable cattle densities and overgrazing (Latchininsky, 1998;Louveaux et al., 1996), and it has been suggested that their frequency might increase in the future favored by global warming (Aragón et al., 2013).
Although genetically nondifferentiated populations have often been interpreted as a single panmictic unit (e.g., Hoarau et al., 2002;Schrey, Schrey, Heist, & Reeve, 2008), the homogenizing effects of gene flow do not necessarily indicate demographic dependence or synchrony among populations (Chapuis et al., 2011) and evidence of selective processes inferred in this study suggests that some populations might experience idiosyncratic evolutionary dynamics (e.g., Guo et al., 2016;Pujolar et al., 2014). Future longitudinal and functional genomic studies could help to identify the proximate factors and specific genes underlying the observed putative signatures of selection (e.g., Bakkali & Martín-Blázquez, 2018;Wang et al., 2014) and determine whether these are temporally stable or if, on the contrary, they are ephemeral and restricted to one or a few generations due to the homogenizing effects of gene flow (Babin et al., 2017;Laporte et al., 2016;Pujolar et al., 2014). The low levels of genetic differentiation among Iberian populations of the Moroccan locust revealed by our analyses indicate that even the application of thousands of SNP markers will not be able to identify the source of nascent outbreak populations at regional scales (see also Chapuis et al., 2009). However, these markers would be useful to determine the origin of eventual waves of long-distance migrants from remote populations with genotypic differences. The high connectivity among populations of the Moroccan locust within the Iberian Peninsula also points out that the demographic dynamics of the species largely exceed regional and national boundaries. For this reason, and in order to avoid the implementation of local control actions with limited success and high environmental and economic costs, the development of future monitoring and sustainable management programs will require the coordination across the different administration levels and government authorities involved. Overall, the results of this study provide baseline information about some unknown aspects of the biology of the Moroccan locust with implications to develop integrated management practices aimed at reducing the negative impacts of this species that is expected to experience in the near future range expansions and increased outbreak intensity in response to ongoing climate warming and land use alterations (Abrol, 2014;Aragón et al., 2013;Benfekih, Chara, & Doumandji-Mitche, 2002;Lankau et al., 2011).

ACK N OWLED G EM ENTS
We thank the unconditional support of Mercedes París and Vicenta

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R S ' CO NTR I B UTI O N S
MGS, PJC, and JO conceived the study and collected the samples.
MGS and JO designed the study and analyses. MGS performed the laboratory work and analyzed the data guided by JO. MGS and JO led manuscript preparation with input from PJC.