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

Although thermophilous and arid‐dwelling 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 Maghrebian‐Levantine 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 arid‐dwelling faunas.


| INTRODUC TI ON
Narrowly distributed thermophilous, arid-adapted, 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 & Blasco-Zumeta, 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 arid-dwelling 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 & Blasco-Zumeta, 1998).
Miocene-Pliocene 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 Mediterranean-Atlantic 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 post-Messinian colonization and considerable genetic cohesiveness between thermophilous and arid-adapted 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 human-mediated 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 & Blasco-Zumeta, 1998). This end is supported by the presence in Iberian semiarid habitats of multiple relict species shared with Maghrebian, Saharo-Arabian and Irano-Turanian regions and whose distributions in the region likely predate anthropic deforestation (Le Driant & Carlon, 2020;Ribera & Blasco-Zumeta, 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ánchez-Gómez et al., 2013), arthropods (Bolívar, 1897;Pascual & Aguirre, 1996), and vertebrates (Barrientos et al., 2009;. Two main hypotheses have been postulated to explain the distribution of these faunas and floras in arid habitats from southern Europe: (i) long-term persistence of relictual biotas that presented a much wider distribution during the Miocene-Pliocene 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 sea-level oscillations Ribera & Blasco-Zumeta, 1998;Sanmartín, 2003).
In the present study, we integrate genome-wide nuclear data with eustatic sea-level 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:  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 post-Messinian biogeographic connections between the two continents despite the barrier effect of the Mediterranean Sea.

K E Y W O R D S
Dericorys carthagonovae, Dericorys millierei, genetic fragmentation, palaeodistribution modelling, phylogenomic inference, Pleistocene glacial cycles, transmarine dispersal 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ínez-Solano et al., 2004;Ribera & Blasco-Zumeta, 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 stepping-stone 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 red-listed D.
carthagonovae is compatible with human-induced habitat fragmentation or, alternatively, a consequence of ancient processes predat- F I G U R E 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 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;

| Genomic library preparation and genomic data processing
We used NucleoSpin Tissue (Macherey-Nagel) kits to extract and purify DNA from a hind leg of each individual. We processed genomic DNA into one genomic library using the double-digestion restriction-site associated DNA sequencing procedure (ddRAD-seq) 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 7-bp barcodes to the digested fragments of each individual. We pooled ligation products and size-selected 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 High-Fidelity DNA Polymerase (BIO-RAD) and sequenced the library in a single-read 151-bp 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(Catchen et al., , 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.

| 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 burn-in step TA B L E 1 Locality, country, code, latitude, longitude, number of genotyped individuals (n), and genetic diversity statistics (π, nucleotide diversity; Hd, haplotype -gene -diversity;  of 100,000 iterations, and assuming correlated allele frequencies and admixture. We conducted 15 independent runs for each value of K-clusters, 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.

| Phylogenomic inference
Second, we used phylonetworks (Solís-Lemus et al., 2017) and treemix v. 1.12 (Pickrell & Pritchard, 2012) to assess the potential presence and direction of gene flow between non-sister 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.
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 sea-level 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 sea-levels estimated by Miller et al. (2011) at each time period contained within the high posterior density (HPD) intervals of divergence times and used one-sample t tests to determine whether they significantly differ from sea level at present time (i.e., 0 m a.s.l.).

| Demographic analyses
We used the composite-likelihood simulation-based 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., . 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  (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 non-polymorphic 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 expectation-conditional 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 re-estimating parameters each time (Excoffier et al., 2013).

| 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 (one-way 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.

| Environmental niche modelling
We built an environmental niche model (ENM) to predict the geo-

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

| 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, BOUL-HOCE, and KAIR-MUJI) identified by structure analyses ( Figure S4).

| Phylogenomic inference
Phylogenomic analyses revealed that D. carthagonovae is monophyletic and nested within D. millierei, which is a paraphyletic taxon ( Figure 1; Figure S5)  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 post-divergence gene flow between non-sister lineages (Table S1). phylonetworks and treemix retrieved the same topology than bpp and svdquartets ( Figure S6).

| Divergence time estimation
bpp analyses (A00 model) estimated that all populations diverged from a common ancestor during the early Pleistocene (c. 1.6 Ma; Calabrian age; Figure 2). All populations of the Maghrebian D. The split of the different lineages took place during glacial periods, when sea levels were significantly below current shoreline (one-sample 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).

| 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 T BOT ; 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 interglacial-glacial transition (Riss-Wü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).

| 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 (one-way ANOVAs, Note: Table shows point estimates and lower and upper 95% confidence intervals for each parameter, which include the timing of population size change (T BOT ) and divergence (T DIV1 and T DIV2 ), and mutation-scaled ancestral (θ ANC , θ BOT and θ POLA-AGUI ) 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).

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

| DISCUSS ION
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). ing with sea-level lowering during glacial maxima Noguerales et al., 2021). Although these grasshoppers are good flyers and have a high intrinsic dispersal capacity, the F I G U R E 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 (T BOT ) and divergence (T DIV1 and T DIV2 ), and mutation-scaled ancestral (θ ANC , θ BOT and θ POLA-AGUI ) 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 T

| Transmarine colonization of southeastern Iberia
Our phylogenomic analyses revealed that all populations of the During this time, overseas distances between northern African and F I G U R E 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 southern European landmasses were probably the shortest since the Messinian Salinity Crisis and lower sea levels might have resulted in the emergence of stepping-stone islands (e.g., shoals) and increased the chance of successful passive dispersal by rafting (Houle, 1998;Husemann et al., 2014). It must be noted, however, that estimates of divergence time must be interpreted with extreme caution considering uncertainty around genomic mutation rates and limited sample sizes for demographic analyses ( Table 1). Lack of samples from Algeria, which could be the actual source populations from which southeastern Iberia was colonized (e.g., , also add considerable spatial and temporal uncertainty in stem age estimates (García-Verdugo et al., 2019). In this line, it has been recently suggested that colonization times most likely lie  Figure 2). Another non-mutually exclusive possibility to explain the high genetic diversity of Tunisian populations is that this region is a melting pot of western and eastern lineages (e.g., Dinis et al., 2019), which could also explain the uncertain phylogenetic position of Jordanian populations (Figure 1; Figure S5). Given that propagule numbers in transmarine dispersal events are expected to be extraordinarily low, a considerable loss of genetic diversity during the colonization of the Iberian Peninsula is, therefore, the most likely scenario (i.e., founder effect; Carson & Templeton, 1984). In this line, genetic diversity of Iberian populations is pretty similar to that inferred for populations from Morocco, but much lower than those obtained for the putative source populations in Tunisia (Figure 4). For instance, the spur-thighed tortoise (Testudo graeca) likely colonized southeastern Iberia from western Algeria after the last glacial period  and relict populations of the Barbary thuja (Tetraclinis articulata) in southeastern Iberia genetically cluster together with Maltese and Tunisian populations rather than with the geographically closer populations from Algeria and Morocco (Sánchez-Gómez et al., 2013). In the same line, southern European populations of the viperine water snake (Natrix maura) are more closely related with Algerian and Tunisian populations than with Moroccan populations (Beddek et al., 2018).
The fact that some of these species were utilized by humans in historical times (e.g., Testudo graeca: direct human consumption; Tetraclinis articulata: timber) has led to suggest the possibility of human-mediated dispersal, although this hypothesis was not formally tested (Sánchez-Gómez et al., 2013) or remained unresolved . Our genomic-based age estimates confirm a post-Messinian colonization of the Iberian Peninsula by the ancestor of D. carthagonovae that largely predates the possibility of human-mediated transportation, supporting the hypothesis of transmarine natural dispersal ( Figure 2). Importantly, these results indicate that Pleistocene faunal exchanges between southern Europe and northern Africa did not exclusively take place across the strait of Gibraltar and the Sicilian Channel and suggest that overseas dispersal of terrestrial organisms between the two continents might be much more common than previously thought (Delicado et al., 2014;Husemann et al., 2014;Stöck et al., 2008).
The north to south decline of genetic diversity in populations of D. carthagonovae (Figures 3 and 4; Table 1) suggests that the ancestral founder populations might have arrived in the northernmost portion of the current distribution of the species in the Iberian Peninsula followed by serial founder effects and loss of genetic diversity during southward range expansions (see also . Fine-scale population sampling in combination with detailed landscape genetic and demographic analyses might help to determine whether the observed cline of genetic diversity reflects the expansion history  or if, rather, it is a consequence of local demographic dynamics linked to availability of suitable habitats and more contemporary patterns of genetic connectivity among remnant population of the species (González-Serna et al., 2019). Anyway, our dating analyses clearly indicate that the genetic divergence of the studied populations (>120 ka) largely predate the timing of extensive anthropogenic impacts in the region, supporting that contemporary genetic fragmentation is not a consequence of habitat destruction resulted from recent human activities (González-Serna et al., 2019;Zellmer & Knowles, 2009). In line with studies on other halophilic species (Ortego et al., 2010(Ortego et al., , 2015, strong genetic structure and lack of genetic exchange between nearby populations emphasize the limited realized dispersal of the species and the low chance that extirpated populations are re-colonized by natural dispersal (Hochkirch et al., 2016).

| Taxonomic and conservation implications
Genealogical and divergence time inferences revealed that D.
carthagonovae is nested within D. millierei and indicate that its ancestor arrived to the Iberian Peninsula much later (< 0.5 Ma) than the estimated divergences among main lineages of D. millierei (>1 Ma). This suggests that D. carthagonovae, originally described by Bolívar (1897) as a "variety" of D. millierei and later upgraded to full species status by Kirby (1910) without any justification or re-description, must be synonymized with D. millierei (for a list of synonyms, see Cigliano et al., 2021). The conservation status of Iberian populations of the species should be re-evaluated but always taking into account its high singularity (i.e., it is the only representative of the genus in Europe) and the intrinsic value of occupied habitats, which have experienced a dramatic destruction in the last decades due to the extensive agricultural and urban development of the region (Hochkirch et al., 2016;Verdú et al., 2011; see also Peñas et al., 2011).

| CON CLUS IONS
This study illustrates the potential of integrating genomic data, eustatic sea-level reconstructions, and palaeodistribution modelling at fine temporal resolution to shed light on the processes underlying the distribution of thermophilous and arid-dwelling faunas shared between southern Europe and arid and desertic regions from North Africa. Our analyses revealed a post-Messinian geographical diversification of the crested grasshopper, supporting Pleistocene range expansion and overseas dispersal as the most parsimonious explanation for the current trans-Mediterranean distribution of D. millierei and the colonization of the Iberian Peninsula, respectively. Collectively, these findings highlight the high relevance of North Africa as a source of thermophilous European faunas and support the strong genetic affinities between the two continents despite the potential barrier effect of the Mediterranean Sea (Beddek et al., 2018;Husemann et al., 2014;Rodríguez-Sánchez et al., 2008). More detailed genetic sampling could help to define with more precision the colonization history of the Iberian Peninsula (e.g., Beddek et al., 2018) and experimental crossing attempts between Iberian and Maghrebian populations would allow to confirm the reproductive cohesiveness of their respective populations (Coyne & Orr, 1989;Saldamando et al., 2005). Future genomic studies focused on other co-distributed relict taxa from the semideserts of southeastern Iberia (e.g., Le Driant & Carlon, 2020;Pascual & Aguirre, 1996) would allow to test whether their colonization was the result of concerted (or idiosyncratic) responses to Pleistocene climatic/eustatic fluctuations, which might ultimately help to reach more generalizable conclusions about the processes underlying the biogeographic connections between African and southern European arid-dwelling biotas (Oaks et al., 2013;Papadopoulou & Knowles, 2015).

ACK N OWLED G EM ENTS
We wish to thank to Anna Papadopoulou for her valuable advice in data analyses, José Miguel Aparicio and Nabil Amor for their help

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