Genomic inferences in a thermophilous grasshopper provide insights into the biogeographic connections between northern African and southern European arid-dwelling faunas

Journal of Biogeography. 2021;00:1–15. | 1wileyonlinelibrary.com/journal/jbi Received: 15 March 2021 | Revised: 28 August 2021 | Accepted: 31 August 2021 DOI: 10.1111/jbi.14267 RESEARCH ARTICLE Genomic inferences in a thermophilous grasshopper provide insights into the biogeographic connections between northern African and southern European ariddwelling faunas Joaquín Ortego1 | María José GonzálezSerna2 | Víctor Noguerales3,4 | Pedro J. Cordero2,5 This is an open access article under the terms of the Creat ive Commo ns Attri bution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. Journal of Biogeography published by John Wiley & Sons Ltd. Handling Editor: Greer Dolby 1Department of Integrative Ecology, Estación Biológica de Doñana (EBDCSIC), Seville, Spain 2Grupo de Investigación de la Biodiversidad Genética y Cultural, Instituto de Investigación en Recursos Cinegéticos (IREC, CSICUCLMJCCM), Ciudad Real, Spain 3Department of Biological Sciences, University of Cyprus, Nicosia, Cyprus 4Instituto de Productos Naturales y Agrobiología (IPNACSIC), La Laguna, Spain 5Escuela Técnica Superior de Ingenieros Agrónomos (ETSIA) de Ciudad Real, Universidad de CastillaLa Mancha (UCLM), Ciudad Real, Spain Correspondence Joaquín Ortego, Department of Integrative Ecology, Estación Biológica de Doñana, EBDCSIC, Avda. Américo Vespucio 26, E41092 Seville, Spain. Email: joaquin.or[email protected] Funding information European Social Fund, Grant/Award Number: CGL201125053, CGL201454671P, CGL201680742R and CGL201783433P; Spanish Ministry of Economy, Industry and Competitiveness Abstract Aim: Although thermophilous and ariddwelling relict biotas constitute a singular component of European biodiversity of high conservation value, we still largely ignore their biogeographic history. In this study, we investigate the geographical diversification of the MaghrebianLevantine crested grasshopper and its colonization of semiarid habitats of southeastern Iberia to gain insights into the historical processes underlying the biogeographic connections between northern African and southern European ariddwelling faunas. Location: Mediterranean region. Taxon: Crested grasshoppers Dericorys millierei and Dericorys carthagonovae (Orthoptera: Dericorythidae). Methods: We used genomic data (ddRADseq) to quantify the genetic structure of populations, infer the phylogenetic relationships among them, estimate divergence times, and elucidate the demographic processes accompanying the colonization of southeastern Iberia. Genomicbased inferences were interpreted in the light of eustatic sealevel reconstructions and species’ range dynamics derived from palaeodistribution modelling at fine temporal resolution. Results: Clustering analyses showed a strong genetic structure and phylogenomic inference revealed that Iberian populations are nested within a Maghrebian clade. Molecular dating analyses indicated that all lineages diverged during the Pleistocene (<1.6 Ma), with point estimates coinciding with glacial periods and the accompanying sea level drops. According to palaeodistribution modelling, the species experienced severe range contractions during the coldest stages of the Pleistocene. Main conclusions: Our results indicate that the colonization of the Iberian Peninsula likely took place during the marked sea level drops characterizing the highamplitude climatic oscillations of the late Quaternary (<0.5 Ma), which considerably reduced overseas distances between northern African and southern European landmasses and 2 | ORTEGO ET al. 1 | INTRODUCTION Narrowly distributed thermophilous, aridadapted, and steppedwelling biotas constitute a singular component of European biodiversity, often including relict species that are the only living representatives in the continent at different taxonomic ranks (Husemann et al., 2014; Kajtoch et al., 2016; Ribera & BlascoZumeta, 1998). In some instances, these taxa show remarkable genetic distinctiveness (i.e., vicariant lineages, subspecies or even species) with respect to core distributions in Central Asia or North Africa, representing a unique evolutionary legacy of high conservation value (Husemann et al., 2014; Kajtoch et al., 2016; Kirschner et al., 2020). Despite recent advances, the temporal and geographic origin of thermophilous and ariddwelling European species is not yet well understood, which in part might be due to the extraordinarily dynamic geological history of the region and the difficulty to distinguish among alternative biogeographical scenarios (Noguerales et al., 2021; Ribera & BlascoZumeta, 1998). MiocenePliocene movements of African and Asian continental plates led to the permanent closing of the eastern end of the current Mediterranean Sea (c. 23– 14 Ma) and the temporal closure of the MediterraneanAtlantic seaways during the Messinian Salinity Crisis (c. 5.96– 5.33 Ma; Bialik et al., 2019; Meulenkamp & Sissingh, 2003), which contributed to faunal and floral exchanges between Africa, Europe and Asia (e.g., Faille et al., 2014; Manafzadeh et al., 2014; Sanmartín, 2003). However, several studies have also found postMessinian colonization and considerable genetic cohesiveness between thermophilous and aridadapted biotas of Europe and those from North Africa (Husemann et al., 2014) and Central Asia steppes (Kirschner et al., 2020), indicating that their current disjunct distributions are most likely a consequence of range expansionfragmentation dynamics linked to Pleistocene climatic oscillations (e.g., Habel et al., 2010; Noguerales et al., 2021). In some other cases, the role of humanmediated dispersal or historical introductions in the distribution of some thermophilous organisms in southern Europe cannot be discarded (Husemann et al., 2014). A paradigmatic case of European thermophilous relict biotas is exemplified in the biogeographic connections between the semideserts characterizing southeastern Iberia and arid regions from North Africa and the Middle East (Le Driant & Carlon, 2020). Although it has been long speculated about the anthropic origin of semidesert areas from the Iberian Peninsula, mounting biogeographical evidence points to the persistence of at least some naturally deforested enclaves through the Pleistocene linked to arid spots with gypsum and saline soils (Ribera & BlascoZumeta, 1998). This end is supported by the presence in Iberian semiarid habitats of multiple relict species shared with Maghrebian, SaharoArabian and IranoTuranian regions and whose distributions in the region likely predate anthropic deforestation (Le Driant & Carlon, 2020; Ribera & BlascoZumeta, 1998). In the specific case of southeastern Iberia, these taxa include strictly thermophilous, xerophytic and deserticolous plants (Cabello et al., 2003; Le Driant & Carlon, 2020; SánchezGómez et al., 2013), arthropods (Bolívar, 1897; Pascual & Aguirre, 1996), and vertebrates (Barrientos et al., 2009; Graciá, Giménez, et al., 2013). Two main hypotheses have been postulated to explain the distribution of these faunas and floras in arid habitats from southern Europe: (i) longterm persistence of relictual biotas that presented a much wider distribution during the MiocenePliocene and expanded during the partial desiccation of the Mediterranean Sea in the Messinian Salinity Crisis; (ii) Quaternary colonization linked to recurrent expansionscontractions of suitable habitats and alterations in the proximity between northern African and southern European landmasses fuelled by Pleistocene climatic and eustatic sealevel oscillations (Graciá, Giménez, et al., 2013; Ribera & BlascoZumeta, 1998; Sanmartín, 2003). In the present study, we integrate genomewide nuclear data with eustatic sealevel and palaeodistribution reconstructions at fine temporal resolutions to shed light on the historical processes underlying the geographical diversification and colonization history of thermophilous faunas shared between northern Africa and southern European arid and semiarid habitats. Specifically, we focus on two closely related taxa of thermophilous crested grasshoppers: Dericorys millierei Bonnet & Finot, 1884 and Dericorys carthagonovae Bolívar, 1897 (Orthoptera: Dericorythidae). Dericorys millierei presents a wide transMediterranean distribution, with a continuous range across the Maghrebian region (Morocco, Algeria, Tunisia and Libya) and disjunct populations in a small area of the Middle East (Israel, Palestine and Jordan; Figure 1). In contrast, D. carthagonovae is a narrowly distributed taxon exclusively present in semiarid areas of southeastern Iberia (Figure 1), where it forms highly fragmented populations linked to vegetation growing in salty and brackish grounds (Verdú et al., 2011). The narrow distribution of D. carthagonovae and the continuous decline of its populations due to extensive destruction of suitable habitats for agricultural and urban might have eased transmarine exchanges of terrestrial faunas. These findings emphasize the high relevance of the Maghreb region as a source of European thermophilous biotas and corroborate postMessinian biogeographic connections between the two continents despite the barrier effect of the Mediterranean Sea. KEYWORDS Dericorys carthagonovae, Dericorys millierei, genetic fragmentation, palaeodistribution modelling, phylogenomic inference, Pleistocene glacial cycles, transmarine dispersal | 3 ORTEGO ET al. development has led to the inclusion of the species in the IUCN Red List of Threatened Species with the category “Endangered” (Hochkirch et al., 2016; Verdú et al., 2011). Although D. carthagonovae is a singular species of high conservation concern, being the only representative of the genus in Europe, its taxonomic status is controversial. The species was first recorded in southeastern Iberia by Bolívar (1897), who described it as a “variety” of D. millierei. The taxon was subsequently upgraded to species rank without any justification in the synonymic catalogue by Kirby (1910), a status that has been accepted and used since then (Cigliano et al., 2021). We first tested the contrasting hypotheses that the current ranges of the two taxa are a consequence of Miocene persistence followed by vicariance events after the Messinian Salinity Crisis (>5.3 Ma; e.g., MartínezSolano et al., 2004; Ribera & BlascoZumeta, 1998) or if, instead, their distributions resulted from more recent pulses of range expansion and fragmentation linked to Pleistocene climatic oscillations (<2.6 Ma; e.g., Fritz et al., 2009; Noguerales et al., 2021; Stöck et al., 2008). Second, we tested the hypothesis that the colonization of the Iberian Peninsula took place coinciding with sealevel lowering during glacial periods, which might have increased the chance of successful passive dispersal by rafting and steppingstone dispersal (Houle, 1998; Husemann et al., 2014). Finally, we quantified spatial patterns of genetic structure and tested whether the timing of genetic subdivision among populations of the redlisted D. carthagonovae is compatible with humaninduced habitat fragmentation or, alternatively, a consequence of ancient processes predating the impacts of anthropogenic activities (GonzálezSerna et al., 2019; Zellmer & Knowles, 2009). FIGURE 1 (a) Phylogenetic relationships among populations as inferred by snapp (3115 SNPs), (b) genetic assignments of populations based on the Bayesian method implemented in the program structure, and (c) map showing the approximate distribution ranges (shaded areas) and geographical location of sampling localities (dots) for Dericorys carthagonovae (in blue) and D. millierei (in red) across the Mediterranean region. snapp tree shows the first (blue) and second (red) most supported topologies and Bayesian posterior probabilities are indicated on the nodes (* = 1). structure analyses were run for all populations (10,000 SNPs) and independently for populations of D. carthagonovae (10,000 SNPs, pie charts on right). Map in EPSG:4326 (WGS84) projection and population codes as described in Table 1. Inset image shows a male of D. carthagonovae (picture by Francisco Rodríguez) MARR BOUL HOCE KAIR N 400 km K= 2 K= 3 * * * * K= 4 K= 5 K= 6 AGUI GATA MUJIPOLA * * >2,000 m <100 m (a) (c) (b) HOCE BOUL KAIR MUJI MARR POLA GATA AGUI 4 | ORTEGO ET al. 2 | MATERIALS AND METHODS 2.1 | Population sampling We sampled populations of Dericorys carthagonovae in the Iberian Peninsula (Spain, n = 3 populations) and Dericorys millierei in the Maghreb (Morocco and Tunisia, n = 4 populations) and the Middle East (Jordan, n = 1 population; Table 1; Figure 1). We used occurrence records available in the literature to design sampling and collect specimens from populations covering the entire distribution ranges of the two taxa (Figure 1). We obtained genomic data for 35 individuals of the two taxa, with an average of four individuals per locality (range = 2– 5; Table 1). Samples of Dericorys lobata lobata (Brullé, 1840) (10 individuals), Dericorys lobata luteipes Uvarov, 1938 (four individuals), and Dericorys minutus Chopard, 1954 (one individual) collected from the Canary Islands (Table 1) were used as outgroups in phylogenomic analyses. We registered spatial coordinates using a Global Positioning System (GPS) and preserved whole specimens at −20°C in 1500 μl ethanol 96% until needed for genomic analyses. Further details on sampling locations are provided in Table 1. 2.2 | Genomic library preparation and genomic data processing We used NucleoSpin Tissue (MachereyNagel) kits to extract and purify DNA from a hind leg of each individual. We processed genomic DNA into one genomic library using the doubledigestion restrictionsite associated DNA sequencing procedure (ddRADseq) described in Peterson et al. (2012). In brief, we digested DNA with the restriction enzymes MseI and EcoRI (New England Biolabs) and ligated Illumina adaptors including unique 7bp barcodes to the digested fragments of each individual. We pooled ligation products and sizeselected them between 475 and 580 bp with a Pippin Prep machine (Sage Science). We amplified the fragments by PCR with 12 cycles using the iProofTM HighFidelity DNA Polymerase (BIORAD) and sequenced the library in a singleread 151bp lane on an Illumina HiSeq2500 platform at The Centre for Applied Genomics (Toronto, ON, Canada). Raw sequences were demultiplexed and preprocessed using stacks v. 1.35 (Catchen et al., 2011, 2013) and assembled using pyrad v. 3.0.66 (Eaton, 2014; e.g., Ortego et al., 2018). Methods S1 provides all details on sequence assembling and data filtering. 2.3 | Genetic structure analyses We analysed population genetic structure and admixture using structure v. 2.3.3 (Pritchard et al., 2000). We ran two independent structure analyses, one including all populations of D. carthagonovae and D. millierei and another focused on the three populations of D. carthagonovae. In both cases, we ran structure using a random subset of 10,000 SNPs, with 200,000 MCMC cycles after a burnin step TABLE 1 Locality, country, code, latitude, longitude, number of genotyped individuals (n), and genetic diversity statistics (π, nucleotide diversity; Hd, haplotype – gene – diversity; θ, population size parameter) for each sampled species and population Species Locality Country Code Latitude Longitude nπHd θ Dericorys carthagonovae Santa Pola Spain POLA 38.205615 −0.613583 50.0008 0.0552 0.0097 Dericorys carthagonovae Águilas Spain AGUI 37.432109 −1.524628 50.0006 0.0426 0.0052 Dericorys carthagonovae Cabo de Gata Spain GATA 36.780914 −2.231885 40.0003 0.0252 0.0026 Dericorys millierei Marrakesh Morocco MARR 31.649482 −7.926335 50.0012 0.0877 0.0092 Dericorys millierei Boulaajoul Morocco BOUL 32.896938 −4.970910 40.0031 0.1396 0.0146 Dericorys millierei Al Hoceima Morocco HOCE 35.193799 −3.864230 50.0020 0.0934 0.0125 Dericorys millierei Kairouan Tunisia KAIR 35.682503 10.224624 50.0086 0.3440 0.0554 Dericorys millierei Wadi Al Mujib Jordan MUJI 31.446110 35.796170 20.0036 0.1751 0.0226 Dericorys lobata lobata Lanzarote Spain LANZ 28.862974 −13.855643 10 — — — Dericorys lobata luteipes Fuerteventura Spain FUER 28.390268 −13.861184 4 — — — Dericorys minutus Gran Canaria Spain GCAN 28.147701 −15.694875 1 — — — | 5 ORTEGO ET al. of 100,000 iterations, and assuming correlated allele frequencies and admixture. We conducted 15 independent runs for each value of Kclusters, where K ranged from 1 to n + 1 for each dataset of n sampled populations. We retained the ten runs having the highest likelihood for each value of K and evaluated the number of genetic clusters that best describes our data according to log probabilities of the data (LnPr(X|K; Pritchard et al., 2000) and the ΔK method (Evanno et al., 2005), as implemented in structure harvester (Earl & vonHoldt, 2012). 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 population membership. Complementary to Bayesian clustering analyses, we performed a principal component analysis (PCA) as implemented in the r v. 4.0.3 (R Core Team, 2021) package ‘adegenet’ (Jombart, 2008). Before running the PCA, we replaced missing data by the mean frequency of the corresponding allele estimated across all samples (Jombart, 2008). 2.4 | Phylogenomic inference First, we reconstructed the phylogenetic relationships among populations of D. carthagonovae and D. millierei using three independent analytical approaches: snapp v. 1.3 (Bryant et al., 2012), bpp v. 4.1 (Flouri et al., 2018), and svdquartets (Chifman & Kubatko, 2014). Second, we used phylonetworks (SolísLemus et al., 2017) and treemix v. 1.12 (Pickrell & Pritchard, 2012) to assess the potential presence and direction of gene flow between nonsister lineages that might result in conflicting phylogenetic relationships and distort tree topology. Methods S2 provides all details on the specific settings used to perform phylogenomic analyses. 2.5 | Estimation of divergence time We used analysis A00 in bpp to estimate the posterior distribution of divergence times (τ; Flouri et al., 2018; Rannala & Yang, 2003). We ran the analyses using the same dataset and settings considered for tree inference analyses in bpp described in Methods S2. We estimated divergence times using the equation τ = 2μt, where τ is the divergence in substitutions per site estimated by bpp, μ is the per site mutation rate per generation, and t is the absolute divergence time in years (Huang et al., 2020; Walsh, 2001). We considered the mutation rate per site per generation of 2.8 × 10−9 estimated for Drosophila melanogaster (Keightley et al., 2014; e.g., Tonzo et al., 2020). Finally, we used paleo sealevel reconstructions to test whether the colonization of the Iberian Peninsula took place coinciding with the lowering of sea levels during the coldest stages of the Pleistocene (Miller et al., 2011). Specifically, we considered the sealevels estimated by Miller et al. (2011) at each time period contained within the high posterior density (HPD) intervals of divergence times and used onesample t tests to determine whether they significantly differ from sea level at present time (i.e., 0 m a.s.l.). 2.6 | Demographic analyses We used the compositelikelihood simulationbased approach implemented in fastsimcoal2 (Excoffier et al., 2013) to estimate the timing of colonization of southeastern Iberia, which is expected to coincide with a demographic bottleneck (i.e., a founder event) predating in situ geographical diversification (e.g., Graciá, Giménez, et al., 2013). We considered that northernmost populations POLA and AGUI share a most recent common ancestor, as supported by phylogenomic analyses (see Section 3) and the comparatively much lower composite likelihood of pilot runs for alternative topological relationships. We calculated a folded joint site frequency spectrum (SFS) considering a single SNP per locus to avoid the effects of linkage disequilibrium. To remove all missing data for the calculation of the joint SFS, minimize errors with allele frequency estimates and maximize the number of variable SNPs retained, each population group was downsampled to n1 of individuals (i.e., four individuals for POLA and AGUI and three individuals for GATA; Table 1) using the easySFS.py script (I. Overcast, https://github.com/isaac overc ast/ easySFS). The SFS contained 5637 variable SNPs. Because invariable sites were excluded from likelihood calculations (‘removeZeroSFS’ option in fastsimcoal2), we fixed the effective population size for one of the demes (POLA) to enable the estimation of other parameters (Excoffier et al., 2013). The effective population size fixed in the model was calculated from the level of nucleotide diversity (π) and estimates of mutation rate per site per generation (μ; 2.8 × 10−9; Keightley et al., 2014). Nucleotide diversity (π) was estimated from polymorphic and nonpolymorphic loci using dnasp v. 6.12.03 (Rozas et al., 2017). The model was run 100 replicated times considering 100,000– 250,000 simulations for the calculation of the composite likelihood, 10– 40 expectationconditional maximization (ECM) cycles, and a stopping criterion of 0.001 (Excoffier et al., 2013). Point estimates for the different demographic parameters were selected from the replicate with the highest maximum composite likelihood. Finally, we calculated confidence intervals of parameter estimates from 100 parametric bootstrap replicates by simulating SFS from the maximum composite likelihood estimates and reestimating parameters each time (Excoffier et al., 2013). 2.7 | Population genetic diversity We calculated levels of haplotype (gene) diversity (Hd) and nucleotide diversity (π) of the different populations using dnasp and tested whether they differ between taxa (oneway ANOVAs) and are explained by geography (i.e., latitude and longitude; linear regressions). Additionally, we calculated contemporary population size parameters (θ) and their respective 95% high posterior density (HPD) intervals in snapp as detailed for phylogenomic analyses in Methods S2. 6 | ORTEGO ET al. 2.8 | Environmental niche modelling We built an environmental niche model (ENM) to predict the geographic distribution of climatically suitable habitats for D. millierei and D. carthagonovae from the last glacial maximum (LGM, 22 ka) to present. To build the ENM, we used the maximum entropy algorithm implemented in maxent v.3.3.3 (Phillips et al., 2006; Phillips & Dudik, 2008) and the 19 bioclimatic variables from the CHELSA database (as described at http://chels aclima te.org/biocl im/) interpolated to 30arcsec resolution (Karger et al., 2017a, 2017b). To estimate environmental suitability from the LGM to present, we projected the ENM to bioclimatic conditions during the last 22,000 years at 100year time intervals (i.e., from 1990 CE to the LGM). Bioclimatic layers at these temporal snapshots are based on a variant of the CHELSA v. 1.2 algorithm (Karger et al., 2017) on the TraCE21 ka data (Liu et al., 2009) and are available at a high resolution (30arcsec) from the CHELSA database (https://chels aclima te.org/; Karger et al., 2017; Yannic et al., 2020). As several lines of evidence indicate that D. millierei and D. carthagonovae should be synonymized (see Section 4), we built a single ENM based on records available for the two currently recognized taxa. Further details on ENM are presented in Methods S3. 3 | RESULTS 3.1 | Genomic data The average number of reads retained per individual after the different quality filtering steps was 2,021,895 (range = 886,652– 3,682,528 reads; Figure S1). On average, this represented 84% (range = 75%– 86%) of the total number of reads recovered for each individual (Figure S1). Final datasets obtained considering a clustering threshold of sequence similarity of 0.85 (Wclust = 0.85) and discarding loci that were not present in at least 50% individuals (minCov = 50%) contained 18,550 SNPs for the dataset including all populations of D. carthagonovae and D. millierei, and 36,483 SNPs for the dataset only including the three populations of D. carthagonovae. 3.2 | Genetic structure analyses structure analyses including populations of D. carthagonovae and D. millierei 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 = 6 (Figure S2a). For K = 2, the two genetic clusters separated populations of D. carthagonovae and D. millierei (Figure 1; Figure S3). Only the population from Tunisia (KAIR) was admixed, with c. 25% of probability of assignment to the genetic cluster mainly represented in the Iberian D. carthagonovae (Figure 1; Figure S3). Analyses for K = 3– 6 sequentially split the different populations of D. millierei in different genetic clusters, which showed no signatures of genetic admixture among them (Figure 1; Figure S3). Analyses focused on the three populations of D. carthagonovae showed that LnPr(X|K) reached a plateau at K = 3 and ΔK peaked at the same K value (Figure S2b). In these analyses, K = 2 split the southernmost population GATA from POLA and AGUI (Figure 1; Figure S3). For K = 3, the three populations of D. carthagonovae were assigned to different genetic clusters that showed no admixture among them (Figure 1; Figure S3). Principal component analysis (PCA) separated D. millierei from D. carthagonovae along the PC1, whereas populations of D. millierei split along the PC2 in the three main genetic clusters (MARR, BOULHOCE, and KAIRMUJI) identified by structure analyses (Figure S4). 3.3 | Phylogenomic inference Phylogenomic analyses revealed that D. carthagonovae is monophyletic and nested within D. millierei, which is a paraphyletic taxon (Figure 1; Figure S5). The population of D. millierei from Marrakech (MARR) was sister to the remaining populations, including those from the Iberian D. carthagonovae, which shared a most recent common ancestor with the Tunisian population (KAIR) of D. millierei (Figure 1; Figure S5). The only incongruence among the different analyses was the phylogenetic position of the population from Jordan (MUJI). bpp and svdquartets analyses supported that MUJI was sister to the subclade including D. carthagonovae (POLA, AGUI, and GATA) and the Tunisian population (KAIR) of D. millierei (Figure S5). However, snapp analyses supported that MUJI was sister to the clade including D. carthagonovae and the remaining populations of D. millierei (excluding MARR; Figure 1). Although this was the only topology contained in the 95% HPD tree set and all nodes were fully supported, the second most supported topology yielded by snapp was identical to that obtained by bpp and svdquartets (Figure 1). All nodes were also fully supported in bpp analyses (Figure S5). Phylogenetic inference using svdquartets was little affected by different schemes of data filtering and all SNP datasets yielded the same topology (Figure S5; e.g., Noguerales et al., 2018; Takahashi et al., 2014). However, the phylogenetic relationships among populations within the clade including D. carthagonovae (POLA, AGUI, and GATA) and the Tunisian (KAIR) and Jordanian (MUJI) populations of D. millierei were not always well resolved by svdquartets, particularly in those analyses based on matrices retaining a lower number of SNPs (i.e., minCov = 25% and 50%; Figure S5). phylonetworks and treemix analyses showed that models considering a strictly bifurcating tree with no introgression edges (i.e., m = 0) are statistically indistinguishable (ΔAIC < 1.4) or more supported than models with one or more migration events, indicating no evidence for postdivergence gene flow between nonsister lineages (Table S1). phylonetworks and treemix retrieved the same topology than bpp and svdquartets (Figure S6). 3.4 | Divergence time estimation bpp analyses (A00 model) estimated that all populations diverged from a common ancestor during the early Pleistocene (c. 1.6 Ma; | 7 ORTEGO ET al. Calabrian age; Figure 2). All populations of the Maghrebian D. millierei split during the early and middle Pleistocene, starting with the divergence of MARR from the rest of the populations (c. 1.6 Ma; Calabrian age) and ending with the split of HOCE and BOUL (c. 0.25 Ma; Chibanian age; Figure 2). Finally, Iberian populations of D. carthagonovae diverged among them (c. 0.27– 0.14 Ma) and from their sister lineage of D. millierei (KAIR; c. 0.44 Ma) during the middle Pleistocene (Chibanian age; Figure 2). The split of the different lineages took place during glacial periods, when sea levels were significantly below current shoreline (onesample t tests, t < −11.70, p < 0.001 for all nodes; Figure S7). The colonization of the Iberian Peninsula was estimated to take place coinciding with the Mindel glaciation (Figure 2), when the sea level dropped to the minimum value of the entire Pleistocene (−123 m; Figure S7). 3.5 | Demographic analyses fastsimcoal2 analyses showed that the most recent common ancestor of southeastern Iberian populations experienced a demographic bottleneck during the Mindel glaciation, which resulted in a reduction of ancestral effective population sizes by c. 39% (Table 2; Figure 3). Remarkably, the point estimate for the timing of the demographic bottleneck (437 ka) is very similar to the stem age (440 ka) calculated in bpp for the divergence between Tunisian and the most recent common ancestor of southeastern Iberian populations (Figure 2). It must be noted, however, that there is considerable uncertainty around the estimation of parameters for more ancient demographic events (especially θANC and TBOT; Table 2), which can in part be explained by the very small samples sizes available for each population (4– 5 individuals/population; Table 1). In situ geographical diversification of D. carthagonovae was estimated to take place during the last interglacialglacial transition (RissWürm), with an initial split of GATA from the rest of the populations (c. 130 ka) followed shortly after by the divergence between POLA and AGUI (c. 124 ka; Table 2; Figure 3). 3.6 | Population genetic diversity Levels of nucleotide diversity (π) and haplotype diversity (Hd), and estimates of the population size parameter (θ) did not significantly differ between D. millierei and D. carthagonovae (oneway ANOVAs, π: F1,6 = 3.27, p = 0.121; Hd: F1,6 = 4.09, p = 0.090; θ: F1,6 = 2.25, p = 0.184; Table 1; Figure 4). Estimates of genetic diversity were not correlated with latitude (π: r = 0.21, F1,6 = 0.29, p = 0.611; Hd: r = 0.29, F1,6 = 0.54, p = 0.489; θ: r = 0.13, F1,6 = 0.11, p = 0.752) nor longitude (π: r = 0.43, F1,6 = 1.40, p = 0.282; Hd: r = 0.47, F1,6 = 1.75, p = 0.234; θ: r = 0.45, F1,6 = 1.50, p = 0.266). Population of D. millierei from Tunisia (KAIR) stood out for its high levels of genetic diversity, which were significantly higher than those observed in the remaining study populations (onesample t tests; π: t = −14.34, FIGURE 2 Phylogenetic tree and divergence times estimated using bpp for the analysed populations of Dericorys carthagonovae and D. millierei across the Mediterranean region. Bars on nodes indicate 95% highest posterior densities (HPD) of divergence times estimated considering a genomic mutation rate of 2.8 × 10−9 per site per generation and a oneyear generation time. Background colours indicate geological divisions of the Quaternary (red: early Pleistocene; yellow: middle Pleistocene; blue: upper Pleistocene; green: Holocene). Sea level estimates based on Miller et al. (2011). Population codes as described in Table 1. Inset image shows a male of D. carthagonovae (picture by Francisco Rodríguez) Time (Ma) Sea level (m) 1.75 1.00 0.75 0.50 0.25 0.001.50 1.25 MARR BOUL HOCE MUJI AGUI GATA KAIR D. milliereiD. carthagonovae -0 -40 -80 -120 POLA A B C D E G F TABLE 2 Parameters inferred from coalescent simulations with fastsimcoal2 under a model of colonization of southeastern Iberia (i.e., hypothetically coinciding with a demographic bottleneck resulted from a founder event) followed by in situ geographical diversification of Dericorys carthagonovae (see Figure 3 for details) Parameter Point estimate Lower bound Upper bound TBOT 437,434 230,154 959,203 TDIV1 130,165 129,963 166,359 TDIV2 124,281 116,718 134,751 θANC 297,960 55,935 466,391 θBOT 119,098 57,953 132,580 θPOLAAGUI 6991 10,397 42,264 θPOLA 140,746 — — θAGUI 57,150 52,521 62,973 θGATA 41,499 39,356 47,646 Note: Table shows point estimates and lower and upper 95% confidence intervals for each parameter, which include the timing of population size change (TBOT) and divergence (TDIV1 and TDIV2), and mutationscaled ancestral (θANC, θBOT and θPOLAAGUI) and contemporary (θPOLA, θAGUI and θGATA) effective population sizes. Estimates of time are given in units of generations (or years, with 1 generation per year). Note that contemporary effective population size for the population POLA (θPOLA) was calculated from its levels of nucleotide diversity (π) and fixed in fastsimcoal2 analyses to enable the estimation of all other demographic parameters (see Section 2.6 for further details). 8 | ORTEGO ET al. p < 0.001; Hd: t = 12.58, p < 0.001; θ: t = −17.96, p < 0.001; Table 1; Figure 4). 3.7 | Environmental niche modelling A threshold (T) feature class and a regularization multiplier of 2 minimized AICc across the set of tested models. After removing highly correlated variables (r > 0.9) and those with a zero percent contribution, the model retained seven bioclimatic variables (sorted by percent contribution, BIO12: 37.4%; BIO3: 34.5%; BIO15: 9.1%; BIO2: 9.0%; BIO6: 6.6%; BIO5: 3.2%; BIO8: 0.1%). Projections of the ENM to bioclimatic conditions from the LGM to present at 100year intervals revealed that the extent of suitable areas, as identified using the maximum training sensitivity plus specificity threshold for species presence (Liu et al., 2005), sharply increased from the LGM to the onset of the Holocene (c. 12 ka), reached a maximum during the Holocene Climate Optimum (c. 9000 to 5000 years ago), and gradually declined since then (Figure 5). In the same line, environmentally suitable areas for the species were considerably fragmented during the LGM and these became much more connected along the Mediterranean coast of North Africa during the Holocene (Figure 5). 4 | DISCUSSION Our results support the Pleistocene connectivity between northern African and southern European arid habitats and the strong genetic cohesiveness of thermophilous terrestrial faunas shared between the two continents (Husemann et al., 2014; e.g., Graciá, Giménez, et al., 2013; Habel et al., 2010; Noguerales et al., 2021). Divergence time estimates indicate that all transMediterranean populations of D. millierei diverged during the Quaternary (<1.6 Ma) and phylogenomic and demographic analyses place the colonization of southeastern Iberia in the midto late Pleistocene (<0.5 Ma), supporting a transmarine migration event that might have taken place coinciding with sealevel lowering during glacial maxima (Graciá, Giménez, et al., 2013; Noguerales et al., 2021). Although these grasshoppers are good flyers and have a high intrinsic dispersal capacity, the FIGURE 3 Schematic of the demographic model used in fastsimcoal2 analyses to estimate the timing of colonization (i.e., hypothetically coinciding with a demographic bottleneck resulted from a founder event) of southeastern Iberia and in situ geographical diversification of Dericorys carthagonovae. Parameters include the timing of population size change (TBOT) and divergence (TDIV1 and TDIV2), and mutationscaled ancestral (θANC, θBOT and θPOLAAGUI) and contemporary (θPOLA, θAGUI and θGATA) effective population sizes. Point estimates yielded by fastsimcoal2 were used to scale the different time events (T) (left axis) and effective population sizes (θ, proportional to box width). Contemporary effective population size for the population POLA (θPOLA, hatched box) was calculated from its levels of nucleotide diversity (π) and fixed in fastsimcoal2 analyses to enable the estimation of all other demographic parameters (see Section 2.6 for further details). Demographic parameter values and confidence intervals are detailed in Table 2 and population codes are described in Table 1 TDI V1 θBOT TBOT TDI V2 0.5 0.2 0.0 0.4 0.3 0.1 θANC Time (Ma) North South | 9 ORTEGO ET al. marked genetic structure of their contemporary populations suggests that strong dependency on severely fragmented habitats has resulted in ancient disruptions of gene flow even among geographically close populations (Figure 1). 4.1 | Pleistocene dispersal and fragmentation Our phylogenomic and dating analyses supported a Pleistocene divergence (<1.6 Ma) of all populations of D. millierei, indicating that the transMediterranean distribution and marked genetic structure of contemporary populations of the species have been most probably shaped by pulses of population expansionfragmentation linked to the highamplitude climatic oscillations of the late Quaternary (Hewitt, 2004; Figure 5). This adds to the accumulating phylogeographic evidence supporting dispersal across North Africa at different time periods (e.g., Beddek et al., 2018; Escudero et al., 2010; Noguerales et al., 2021; PérezCollazos et al., 2009; Veríssimo et al., 2016), which indicates that this region has served as an important migration corridor for numerous terrestrial organisms and played a major role on faunal and floral exchanges between Asia, Africa and Europe (Husemann et al., 2014; Sanmartín, 2003). The main eastwest split separating Moroccan from Tunisian populations is also congruent with findings from numerous previous studies in which the Moulouya River valley (AlgeriaMorocco border) and the Kabylia region (central Algeria) have been identified as the main phylogeographic breaks across numerous organismal groups (e.g., land snails: Guiller & Madec, 2010; reptiles and amphibians: Beddek et al., 2018; plants: SánchezGómez et al., 2013; Taib et al., 2020). Formal testing of concordance in divergence times across codistributed taxa would help to distinguish among alternative biogeographical hypotheses and understand whether spatially similar phylogeographic structures correspond to contrasting evolutionary processes or to analogous responses to the past geological and/or climatic dynamics of the region (e.g., Papadopoulou & Knowles, 2015; Wan et al., 2021). The limited realized dispersal of D. millierei, as evidenced by the pronounced genetic structure of its populations, points to range fragmentation, rather than longdistance dispersal, as the most likely explanation for the current distribution gap of the species in Egypt (Figure 1; Cigliano et al., 2021). Accordingly, our palaeodistribution reconstructions at fine temporal resolution supported the presence of corridors of suitable habitats across the Mediterranean coast of Africa that connected the central Maghreb region and the Middle East during the warmer stages of the Pleistocene (i.e., the Holocene Climate Optimum; Figure 5d). This distribution gap in northeastern Africa is analogous to that reported for two other ariddwelling taxa with transMediterranean disjunct distributions of Pleistocene origin: the saltmarsh grasshopper (Mioscirtus wagneri; Noguerales et al., 2021) and the spurthighed tortoise (Testudo graeca; Fritz et al., 2009). Collectively, these results suggest that the contraction of suitable habitats in northeastern Africa is the most parsimonious explanation for the contemporary disjunct distributions observed in some thermophilous organisms that probably presented wider distributions across the Mediterranean region during the warmest stages of the Pleistocene (Noguerales et al., 2021; Ribera & BlascoZumeta, 1998). 4.2 | Transmarine colonization of southeastern Iberia Our phylogenomic analyses revealed that all populations of the Iberian D. carthagonovae are monophyletic and embedded within the MaghrebianLevantine D. millierei, which is a paraphyletic taxon (Figure 2). This finding is in agreement with previous studies showing that southern European lineages of numerous thermophilous taxa are nested within North African clades (reviewed in Husemann et al., 2014). Estimates of divergence time indicate that the Iberian Peninsula was colonized from the Maghreb region during the Pleistocene (<0.5 Ma), supporting a postMessinian transmarine dispersal event. The split of the different lineages of both D. millierei and D. carthagonovae probably took place during the coldest stages of the Pleistocene (i.e., glacial periods), when sea levels were below the current shoreline (Figure 2) and populations became highly fragmented according to palaeodistribution modelling (Figure 5f). Remarkably, phylogenomic and demographic analyses suggest that southeastern Iberia was colonized from the Maghreb region coinciding with the severe Mindel glaciation (c. 440 ka; Figures 2 and 3), a period that has been estimated to present the lowest sea level of the entire Pleistocene (−123 m; Figure S7; Miller et al., 2011). During this time, overseas distances between northern African and FIGURE 4 Estimates of the population size parameter (θ) (median ± 95% high posterior density intervals) inferred by snapp for the analysed populations of Dericorys carthagonovae and D. millierei across the Mediterranean region. Population codes as described in Table 1 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 D. milliereiD. carthagonovae MARR BOUL HOCE MUJI AGUI GATA KAIR POLA Theta (θ)