Geographical isolation versus dispersal: Relictual alpine grasshoppers support a model of interglacial diversification with limited hybridization

Alpine biotas are paradigmatic of the countervailing roles of geographical isolation and dispersal during diversification. In temperate regions, repeated distributional shifts driven by Pleistocene climatic oscillations produced both recurrent pulses of population fragmentation and opportunities for gene flow during range expansions. Here, we test whether a model of divergence in isolation vs. with gene flow is more likely in the diversification of flightless alpine grasshoppers of the genus Podisma from the Iberian Peninsula. The answer to this question can also provide key insights about the pace of evolution. Specifically, if the data fit a divergence in isolation model, this suggests rapid evolution of reproductive isolation. Genomic data confirm a Pleistocene origin of the species complex, and multiple analytical approaches revealed limited asymmetric historical hybridization between two taxa. Genomic‐based demographic reconstructions, spatial patterns of genetic structure and range shifts inferred from palaeodistribution modelling suggest severe range contraction accompanied by declines in effective population sizes during interglacials (i.e., contemporary populations confined to sky islands are relicts) and expansions during the coldest stages of the Pleistocene in each taxon. Although limited hybridization during secondary contact leads to phylogenetic uncertainty if gene flow is not accommodated when estimating evolutionary relationships, all species exhibit strong genetic cohesiveness. Our study lends support to the notion that the accumulation of incipient differences during periods of isolation were sufficient to lead to lineage persistence, but also that the demographic changes, dispersal constraints and spatial distribution of the sky islands themselves mediated species diversification in temperate alpine biotas.


| INTRODUC TI ON
The opportunities for divergence in isolation, but also the counteracting effects of gene flow during periods of secondary contact, are quintessential processes of Pleistocene speciation in alpine and montane biotas from temperate regions (Hewitt, 2000). Isolation of populations in glacial  or interglacial (Bennett & Provan, 2008) refugia during the climatic oscillations of the Pleistocene is likely to have exposed them to different selective regimes and increased genetic drift, which ultimately are hypothesized to have promoted divergence and speciation (Hewitt, 1996;Stewart et al., 2010). Conversely, latitudinal displacements to or from glacial refugia and downslope movements towards lower elevation areas during ice ages may contribute to geographical contact of gene pools that had remained isolated over extended periods of time (Knowles & Massatti, 2017;Maier et al., 2019;Tonzo & Ortego, 2021). If incipient speciation during geographical isolation is not accompanied by effective reproductive isolation, secondary contact will lead to post-divergence gene flow (Hewitt, 2000). Depending on the permeability to gene exchange between previously allopatric lineages, the consequences of such process will range from speciation reversal (e.g., Maier et al., 2019) to different levels of introgressive hybridization (e.g., Melo-Ferreira et al., 2005). For these reasons, Pleistocene glacial-interglacial cycles are hypothesized to have acted both as "species pumps" and as "melting pots," creating opportunity for divergence and gene exchange across different stages along the continuum of speciation (April et al., 2013;Ebdon et al., 2021;Haffer, 1969;Hewitt, 1996Hewitt, , 2000Knowles, 2001;Petit et al., 2003). Climate oscillations during the Quaternary are thus expected to have promoted reticulate speciation in many organism groups, rather than a strictly bifurcating evolutionary history of species divergence (Nevado et al., 2018;Thom et al., 2018).
An accurate reconstruction of the history of species divergence is a prerequisite for inferring the tempo and mode of speciation and testing alternative biogeographical and macro-evolutionary hypotheses regarding the processes underlying observed patterns of biological diversity (Nylander et al., 2008;Rangel et al., 2015;Tariel et al., 2016). However, the phylogenetic relationships of recently diverged species can be obscured by unresolved nodes (i.e., polytomies) (e.g., Kutschera et al., 2014;Takahashi et al., 2014). Phylogenetic uncertainties are frequently a consequence of incomplete lineage sorting (ILS) of ancestral polymorphism (Maddison & Knowles, 2006) and/ or deviations from strictly bifurcating lineages due to horizontal gene transfer, hybrid speciation or introgression (Mallet et al., 2016;McBreen & Lockhart, 2006). Identifying the causes of phylogenetic conflict (i.e., reticulation vs. ILS) is essential for distinguishing among alternative evolutionary pathways (e.g., de Manuel et al., 2016;Schrago & Seuanez, 2019;Thom et al., 2018), which can ultimately provide key insights about the pace of speciation (Rosindell et al., 2010;Sukumaran & Knowles, 2017). However, this task is more daunting in recent Pleistocene radiations in which species may have weak reproductive barriers and short interspeciation times that often co-occur with secondary introgression (i.e., post-divergence gene flow) (Nevado et al., 2018;Wen et al., 2016). In the last decade, increased capacity to generate large genomic data sets in nonmodel organisms has been critical to overcoming statistical uncertainties contributed by limited genetic information. This has also driven the development of numerous analytical approaches aimed at resolving gene tree discordances and detecting admixture (reviewed in Payseur & Rieseberg, 2016). Thanks to these analytical advances, we can now test speciation hypotheses that depart from models of divergence in strict isolation, which is key to considering whether introgressive hybridization is a component of the evolutionary portrait of diversification (e.g., de Manuel et al., 2016;Thom et al., 2018).
However, given the assumptions and limitations inherent to each available approach, corroboration across multiple lines of evidence for complex histories of diversification is recommended (for details, see Payseur & Rieseberg, 2016).
Here, we use genomic data to evaluate the countervailing effects of dispersal and isolation on the speciation process and determine whether population isolation driven by Pleistocene glacial cycles triggered the necessary mechanisms for long-lasting genetic cohesion of lineages or if, on the contrary, extensive gene flow during periods of secondary contact have impaired their persistence and the formation of new species (Dynesius & Jansson, 2014). Specifically, we apply an integration of multiple approaches to unravel a Pleistocene diversification history of alpine grasshoppers of the genus Podisma from the Iberian Peninsula (Orthoptera: Acrididae) (Morales-Agacino, 1951). As the southernmost distributional limit of the genus, the region hosts three species distributed in allopatry across different mountain ranges (Cigliano et al., 2021; Figure 1).
The three taxa are distributed at elevations >1200 m, restricted to montane and alpine open habitats dominated by low grasslands and dwarf shrub formations (e.g., Juniperus sp., Vaccinium sp. and Rhododendron sp.), which are interspersed with patches of bare ground and rocks (Presa et al., 2016a(Presa et al., , 2016bZuna-Kratky et al., 2016). As such, their contemporary populations are extremely fragmented across sky islands of suitable habitat embedded in an inhospitable matrix characteristic of the Mediterranean climate (Cigliano et al., 2021;Presa et al., 2016aPresa et al., , 2016bZuna-Kratky et al., 2016).
There are no clear phenological differences among the taxa; the three species are univoltine (i.e., a single generation per year) with adult populations peaking in July-August (Morales-Agacino, 1951;Zuna-Kratky et al., 2016;J. Ortego, personal observation). They are very similar in external appearance and are flightless, with Podisma pedestris and P. carpetana being micropterus and P. cantabricae apterous, although rare macropterous (i.e., fully winged) forms have occasionally been described in P. pedestris (Lemonnier-Darcemont & Darcemont, 2014;Morales-Agacino, 1951). Capture-mark-releaserecapture studies on P. pedestris indicate that these taxa exhibit very low dispersal capacities and a marked philopatric behaviour (Barton & Hewitt, 1981, 1982Mason et al., 1995). Due to this limited dispersal ability and strict habitat requirements, we hypothesize that recurrent pulses of population expansions and contractions during Pleistocene glacial cycles have contributed to genetic isolation and speciation, but that the shifting distributions also generated repeated opportunities for post-divergence gene flow (e.g., Barton, 1980;Keller et al., 2008). To accommodate an evolutionary history that may depart from assumptions of divergence in isolation and to gain insights into the processes underlying speciation that includes the possibility of post-divergence gene flow, we integrate a comprehensive suite of phylogenomic and population genomic approaches with palaeoclimate-based reconstructions of species distributions. Specifically, we apply the multispecies coalescent (MSC) model to infer phylogenetic relationships among taxa and identify nodes with potential conflict that might be indicative of either ILS or reticulation. Then, we perform phylogenetic tests to distinguish ILS from introgression and use a model-based approach to evaluate alternative scenarios of post-divergence gene flow or lack thereof.
Using environmental niche modelling and palaeoclimate-based reconstructions of species distributions, we infer range shifts during glacial/interglacial periods in each species and use this framework to determine which expectations in terms of population fragmentation and secondary contact are most probable given species divergence, past demography and introgressive hybridization estimated based on the genomic data.

| Population sampling
Occurrence records from the literature were used to design sampling and guide collection of specimens from populations repre-  Table S1). Seven specimens of Cophopodisma pyrenaea (Fischer, 1853) (tribe Podismini ;Cigliano et al., 2021) were used as an outgroup in phylogenomic analyses and ABBA/BABA tests (Table S1). Spatial coordinates were recorded using a Global Positioning System (GPS) and whole specimens were preserved at −20°C in 1500 µl of 96% ethanol until needed for genomic analyses. F I G U R E 1 (a) Results of genetic assignments for Podisma pedestris, P. carpetana and P. cantabricae based on the Bayesian method implemented in the program structure. Hierarchical structure analyses were run for all species and independently for populations of P. pedestris and P. carpetana. Analyses are based on a random subset of 10,000 SNPs. Thin vertical lines separate individuals and thick lines demarcate sampling sites, with each individual partitioned into K coloured segments proportional to the individual's estimated ancestry proportions; population codes are described in Table S1

| Genomic library preparation and 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 restrictionsite associated DNA sequencing procedure (ddRADseq) described in Peterson et al. (2012). In brief, we digested DNA with the restriction HiSeq2500 platform at The Centre for Applied Genomics (Toronto, ON, Canada). Raw sequences were demultiplexed and preprocessed using stacks version 1.35 (Catchen et al., 2013) and assembled using pyrad version 3.0.66 (Eaton, 2014); see Methods S1 for details on sequence assembling and data filtering. The choice of different filtering and assembling thresholds had a little impact on the obtained inferences (see Section 3; e.g., Eaton, 2014;Ortego et al., 2018). For this reason, unless otherwise indicated, all downstream analyses were performed using data sets of unlinked single nucleotide polymorphisms (SNPs) (i.e., using a single SNP per RAD locus) obtained with pyrad considering a clustering threshold of sequence similarity of 0.85 (W clust =0.85) and excluding loci that were not present in at least 20 individuals (minCov = 20). We used the option related-ness2 in vcftools to calculate the relatedness among all pairs of genotyped individuals and to exclude the possibility that we had sampled close relatives within each study population (Danecek et al., 2011;Manichaikul et al., 2010).

| Quantifying genetic structure
We analysed population genetic structure and admixture using the Bayesian Markov chain Monte Carlo (MCMC) clustering method implemented in the program structure version 2.3.3 (Pritchard et al., 2000). We conducted structure analyses hierarchically, initially analysing data from all populations and species jointly and, subsequently, running independent analyses for subsets of populations assigned to the same genetic cluster in the previous hierarchical-level analysis (Janes et al., 2017;Pritchard et al., 2000). We ran structure using a random subset of 10,000 SNPs with 200,000 MCMC cycles after a burn-in step of 100,000 iterations, and assuming correlated allele frequencies and admixture (Pritchard et al., 2000). We performed 15 independent runs for each value of K genetic clusters, where K ranged from 1 to n + 1 for each data set of n populations, to estimate the most probable number of clusters. We retained the 10 runs with the highest likelihood for each K-value. As recommended by Gilbert et al. (2012) and Janes et al. (2017), we used two statistics to interpret the number of genetic clusters (K) that best describes our data: log probabilities of Pr(X|K) (Pritchard et al., 2000) and ΔK (Evanno et al., 2005). These statistics were calculated as implemented in structure harvester (Earl & vonHoldt, 2012). We used clumpp version 1.1.2 and the Greedy algorithm to align multiple runs of structure for the same K-value (Jakobsson & Rosenberg, 2007) and distruct version 1.1 (Rosenberg, 2004) (Jombart, 2008). Before running the PCAs, we replaced missing data by the mean frequency of the corresponding allele estimated across all samples (Jombart, 2008).

| snapp
For the snapp analyses, we ran two independent replicates of >2 mil-

| svdquartets
We ran svdquartets to estimate the evolutionary relationships of populations from each species (i.e., a population/species tree;  by evaluating 10,000 random quartets from the data set; uncertainty in relationships was quantified using 100 bootstrapping replicates. Given the low computational burden of svdquartets in comparison with snapp, we analysed six SNP matrices obtained by setting different values of clustering thresholds (W clust = 0.85 and 0.90) and minimum taxon coverage for a locus (minCov = 10, 20 and 30) (see Methods S1). This allowed us to assess the impact of different proportions of missing data and number of loci on the estimated topology and patterns of branch support (Huang & Knowles, 2016;Noguerales et al., 2018;Takahashi et al., 2014).

| Analyses of introgression
Although phylogenetic analyses tended to support P. carpetana and P. cantabricae as sister taxa, the relationships among the three species often remained unresolved (see Section 3.3). To determine the role of ILS versus introgression in explaining such conflicting phylogenetic relationships, we tested the possibility of post-divergence gene flow using a comprehensive suite of approaches detailed below.
Note that these analyses (with the exception of phylonetworks) were carried out sequentially using one representative population per species, and then testing all population combinations, because intraspecific population structure (see Figure 1) can confound analyses that assume panmixia within species.

| Phylogenetic networks
We used Species Networks applying Quartets (snaq) implemented informative (Bernardes et al., 2007). We used these gene trees and phylonetworks to estimate quartet concordance factors (CFs), defined as the proportion of genes that support each possible relationship between each set of four taxa. We used the topology obtained with snapp as a starting tree and estimated the best phylogenetic network testing a varying number of reticulation events (h from 0 to 5), each optimized with 10 independent runs. The optimal number of reticulation events was chosen using a heuristic approach by plotting negative pseudo-likelihood scores against h-values, as recommended by the authors (Solís-Lemus et al., 2017).

| D-statistics
We used four-taxon ABBA/BABA tests based on the D-statistic to test for introgression as an explanation for conflicting phylogenetic relationships (Durand et al., 2011). Briefly, for the sister species P1 and P2 (i.e., P. cantabricae and P. carpetana, respectively), which diverged from a common ancestor with P3 (i.e., P. pedestris), and the outgroup O (i.e., C. pyrenaea), the D-statistic is used to test the null hypothesis of no introgression (D = 0) between P3 and P1 or P2. D-values significantly different from zero indicate gene flow between P1 and P3 (D < 0) or between P2 and P3 (D > 0). We performed ABBA/BABA tests in pyrad and used 1000 bootstrap replicates to obtain the standard deviation of the D-statistic and significance levels (Eaton & Ree, 2013). We ran ABBA/BABA tests sequentially for each of the six different speciespopulation combinations (i.e., using each population as a representative for a species). Only populations with six or more genotyped individuals were considered for these analyses (see Table S1). We ran these analyses using six different genetic data sets obtained by setting different clustering thresholds (W clust =0.85 and 0.90) and minimum taxon coverage for a given locus (minCov =10, 20 and 30).

| Population graphs
We analysed the potential presence of introgression and determined the direction of gene flow using treemix version 1.12 (Pickrell & Pritchard, 2012). treemix fits a population graph based on population allele frequencies and a Gaussian approximation to genetic drift, inferring patterns of splits and admixtures. We ran treemix analyses considering the same six species-population combinations used for ABBA/BABA tests, assuming independence of all SNPs with a window size of one SNP (k = 1). Using an estimated ML tree rooted with the outgroup C. pyrenaea, we tested a range of migration events (m from 0 to 4) and determined the best fit model for the data by plotting Ln(likelihood) scores against m-values.

| Models of interspecific gene flow
We used fastsimcoal2 (Excoffier et al., 2013) to evaluate the fit of the data to 10 alternative divergence models that considered different scenarios of interspecific gene flow (see Figure S1); the timing of gene flow was modelled as a time interval, with an estimate for the time gene flow was initiated (T INTROG1 ) and the time that it ended (T INTROG2 ; Figure S1). We estimated the composite likelihood of the observed data (analysing one SNP per locus) given a specified model using the site frequency spectrum (SFS) and the simulation-based approach implemented in fastsimcoal2 (Excoffier et al., 2013). Separate analyses of one population per species were performed considering six different species-population combinations, as done for ABBA/BABA tests and treemix analyses. Because invariable sites were not included in the SFS, we fixed the effective population size for one species (P. cantabricae) to enable the estimation of other parameters in fastsimcoal2 (Excoffier et al., 2013); the fixed effective population size was calculated from the level of nucleotide diversity (π = 0.0005) and the mutation rate per site per generation (2.8 × 10 −9 ) estimated for Drosophila melanogaster (Keightley et al., 2014), which is similar to the spontaneous mutation rate estimated for the butterfly Heliconius melpomene (2.9 × 10 −9 ; Keightley et al., 2015). To remove all missing data for the calculation of the joint SFS (as required), each population group was downsampled to five individuals using the easySFS.py script (I. Overcast, https://github.com/isaac overc ast/easySFS).
Each 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). We used an information-theoretic model selection approach based on Akaike's information criterion (AIC) to determine the probability of each model given the observed data (Burnham & Anderson, 2002).

| Inference of past demographic history
We reconstructed the past demographic history from the SFS using the program stairway plot version 2.1, which does not require wholegenome sequence data or reference genome information (Liu & Fu, 2015. We computed the SFS for each population as described in the previous section and ran stairway plot fitting a flexible multi-epoch demographic model, considering 1 generation per year (Barton & Hewitt, 1981), assuming a mutation rate of 2.8 × 10 −9 per site per generation (Keightley et al., 2014), and performing 200 bootstrap replicates to estimate 95% confidence intervals.

| Environmental niche modelling
We estimated environmental niche models (ENM) to (i) predict the geographical distribution of climatically suitable areas for the three species both in the present and during the last glacial maxi- have been proven to be enough to develop ENMs with a good predictive power using maxent (e.g., Papes & Gaubert, 2007;van Proosdij et al., 2016;Wisz et al., 2008). We used the r package ENMeval (Muscarella et al., 2014) to conduct species-specific parameter tuning and determine the optimal feature class (FC) and Climatically suitable areas for each species and time period were identified by converting the logistic outputs from maxent into binary maps using the maximum training sensitivity plus specificity (MTSS) threshold value for occurrence (Liu et al., 2005).

| Quantifying genetic structure
For the structure analyses, LnPr(X|K) plateaued at K = 3 and ΔK peaked at the same K-value ( Figure S3a), which corresponds to the three taxa, with no sign of genetic admixture among them (i.e., individual and population probabilities of membership = 1; Figure 1a). structure analyses performed separately on P. pedestris and P. carpetana revealed a strong population genetic structure within each species ( Figure 1a). Two genetic clusters inferred for P. pedestris ( Figure S3b) group individuals by the two analysed populations for this species (AUL and AIG), with no signs of genetic admixture (Figure 1a). For P.
carpetana, the most likely number of clusters was K = 2 according to the ΔK criterion, but LnPr(X|K) steadily increases up to K = 4 ( Figure   S3c). These analyses reveal a north-south hierarchical genetic structure, with signs of admixture restricted to some nearby populations from the Iberian System (DEM-URB and URB-MON; see Figure 1).
PCAs separate well the three taxa and most populations within taxa, supporting the results from structure ( Figure S4).

| Phylogenomic inference
The monophyly of all taxa and the same species relationships were estimated by both snapp and svdquartets (Figures 2 and S5).
Phylogenetic relationships among species are well supported with snapp (posterior probability >0.98; Figure 2a), but not with svdquartets (Figures 2b and S5), although the estimates from svdquartets are robust to different schemes of data filtering and assembling ( Figure   S5). The phylogenetic relationships among geographically proximate populations of P. carpetana were not well resolved by either snapp or svdquartets (Figure 2; see also Figure S5). In snapp, the three topologies contained in the 95% HPD tree set differed only in the population relationships inferred for P. carpetana (Table S2). These unresolved population relationships within P. carpetana are not unexpected given evidence of gene flow among nearby populations located in the same mountain range from structure ( Figure 1a) and PCAs ( Figure S4).  Table S1. Picture shows a male of Podisma pedestris, which is morphologically similar to the other two species

| Phylogenetic networks
phylonetwork analyses revealed that all models involving reticulation events (h > 0) fit our data better than models considering strict bifurcating trees (h = 0) ( Figure S6). Negative pseudo-likelihood scores decrease sharply from h = 0 to h = 2 and remain unaltered or with a very small improvement for h > 2 ( Figure S6), suggesting that the best-fitting phylogenetic model includes two introgression events. One inferred introgression event (γ A ) is from P. pedestris into P. carpetana, with ~11% of gene copies in the ancestor of P. carpetana traced to the ancestor of the two populations of P. pedestris

| D-statistics
A statistically significant excess of ABBA patterns (D > 0) supports post-divergence gene flow between P. pedestris (P3) and P.
carpetana (P2) ( Table 1). This result holds irrespective of which population-species combinations were analysed, or the data filtering and assembling scheme used in generating the data set (Table S3).

| Population graphs
treemix analyses consistently support a single migration event ( Figure  S7) of directional gene flow from P. pedestris to P. carpetana ( Figures   3 and S8).

| Models of interspecific gene flow
fastsimcoal2 analyses performed for all population-species combinations consistently show that the best supported scenario is one with asymmetrical gene flow from P. pedestris to P. carpetana (Figure 4; Model B in Table S4). Considering the 1-year generation time of these species (Barton & Hewitt, 1981), the split between P. pedestris and the two other taxa is estimated to have taken place ~638,000-992,000 years ago, during the early-middle Pleistocene (Figure 4; Table S5). The split between P. carpetana and P. cantabricae is estimated as ~131,000-155,000 years ago, during the middle Pleistocene ( Figure 4; Table S5).
Gene flow from P. pedestris to P. carpetana is inferred to have taken place during the middle-late Pleistocene, between ~108,000-147,000 and 87,000-120,000 years ago (Figure 4; Table S5). It should be noted that estimates for the 95% confidence intervals for the oldest demographic parameters (θ ANC , θ CAR-CAN and T DIV1 ) are much wider than those for more recent events (θ PED , θ CAR , T DIV2 , T INTROG1 and T INTROG2 ) ( Figure 4; Table S5), which is consistent with the lower accuracy of fastsimcoal2 to estimate more ancient events, such as those involving species formation (Excoffier et al., 2013).

| Inference of past demographic history
stairway plot analyses suggest the three species experienced parallel demographic responses to climate warming since the end of the last glacial period ( Figure 5). More specifically, all analysed populations from the three species show demographic declines that generally follow the LGM and reduced their effective population sizes (N e ) by >95% ( Figure 5). We note that these population size estimates differ from those of the parameterized divergence model (Figure 4), but that the divergence model did not include population size change parameters because of the complexity it would have added to the alternative tested models (Knowles, 2009).

| Environmental niche modelling
The low OR MTP (OR MTP <0.01) and high AUC TEST (AUC TEST >0.99) for the ENM of each species indicate their high discriminatory power and low degree of overfitting, respectively (for details on model performance and parameters, see Table S6). Climatically suitable areas predicted by ENMs yield distribution patterns highly congruent with the present-day observed distributions for the three species ( Figure 6).
Only very small areas in mountain ranges far from the current distribution of each species are (over-)predicted as suitable ( Figure 6).

| Determinants of species pump or meltingpot processes
In areas with temperate climates, such as the Mediterranean region, cold-adapted species with narrow climatic niches are currently limited to small and isolated patches of high-elevation habitat (i.e., sky islands; Flantua et al., 2020;Knowles & Massatti, 2017). conditions into what is now unsuitable habitat is also predicted to have driven expansion of cold-adapted species (Hewitt, 2000). LGM. This is consistent with the current distribution of climatically suitable habitats and the low dispersal capability of the species that nowadays persist in small and highly fragmented interglacial refugia (Bennett & Provan, 2008;Stewart et al., 2010). Thus, both extreme isolation and severe demographic bottlenecks after the LGM created the perfect scenario for genetic differentiation via strong genetic drift and a fragmented population structure.
In the different sky island archipelagos of the Iberian Peninsula, like other montane regions across the globe, the repeated climateinduced distributional shifts and associated demographic conditions (e.g., bottlenecked and fragmented populations) experienced by their biotas represent the quintessential setting for a species pump diversification process (Flantua et al., 2020;Haffer, 1969;Papadopoulou & Knowles, 2015a;Wallis et al., 2016). However, this dynamic can transition into a melting-pot scenario-that is, the repeated cycles of distributional shifts result in the loss of incipient divergences due to gene flow during range expansions (Klicka & Zink, 1997;Maier et al., 2019). This tipping point between distributional shifts promoting vs. inhibiting speciation is expected to vary among species and geographical settings. The inherent dispersal constraints of flightlessness, coupled with climatic adaptation TA B L E 1 Analyses of introgression using four-taxon D-statistic (ABBA/BABA) tests Note: Analyses were performed for each of the six species-population combinations separately using Cophopodisma pyrenaea as an outgroup. All tests were highly significant (q < .001) after a false discovery rate (FDR) adjustment (5%) to control for multiple tests. Population codes are described in Table S1. n, number of retained SNPs; D (± SD), D-statistic and corresponding standard deviation; z, z-statistic; q, p-values adjusted at a FDR of 5%.

F I G U R E 3
Maximum-likelihood tree inferred with treemix (4569 snps) showing the most likely migration event (m = 1). The direction of gene flow is represented with an arrow coloured according to the percentage of alleles (weight) originating from the source. Codes of the specific populations of each species included in the analysis are described in Table S1. Note that analogous treemix analyses were run considering all other species-population combinations (see Figure  The geographical context of climatic-induced distributional shifts of habitats and their constituent inhabitants are also important. That is, the dynamics of colonization itself may be instrumental in determining the distribution of genetic variation and differentiation (see Knowles & Massatti, 2017).
With the sky islands in the Iberian Peninsula embedded in a landscape of unsuitable habitat, the likelihood of dispersal may certainly be reduced, but not impossible, especially with the projections of a more expansive distribution of suitable habitat for temperate taxa in the past ( Figure 6). Yet, the presence of past corridors of suitable habitat among the isolated contemporary populations does not necessarily mean that they were utilized or that the three Podisma species actually came into secondary contact. Low dispersal ability and a marked philopatric behaviour of these flightless species (Barton & Hewitt, 1981, 1982Mason et al., 1995), ecological constraints (e.g., biotic interactions; Hampe, 2004;Ortego & Knowles, 2020), and the strong fragmentation and small sizes of severely bottlenecked populations (Figures 1 and 5) might have hampered their capacity to colonize remote areas that were climatically suitable during glacial periods (Kearney & Porter, 2009;Wiens et al., 2009). For instance, it is likely that the current tiny distribution of P. cantabricae (<25-50 km 2 ; Presa et al., 2016b) in a topographically highly complex region severely limited its capacity to reach predicted environmentally suitable, but  Table S1. Parameter estimates for analogous fastsimcoal2 analyses run considering all other species-population combinations are presented in Table S5. *Note the effective population size of Podisma cantabricae (θ CAN ) was calculated from the level of nucleotide diversity (π) and fixed in fastsimcoal2 analyses (see the Section 2.5.4 for further details) distant areas (i.e., located > 1000 km away from its current range) during the LGM (Figure 6). In sum, incipient divergences may be lost among some sky island populations, but not others, and similarly, the opportunity for past hybridization among currently allopatric Podisma taxa may depend on the contemporary and past geographical configuration of suitable habitats within the dynamic ranges of each species (Knowles & Massatti, 2017;Tonzo & Ortego, 2021).
When the speciation process is viewed through the lens of divergence with gene flow, rather than divergence in isolation, it invites a shift in perspective about the controls on speciation (Aguilée et al., 2018;Harvey et al., 2019). For example, in addition to the traditional focus on the rate of evolution of reproductive isolation as a control on diversification (in which potential gene flow associated with cycles of climate-induced distributional shifts is thwarted), given that divergence in Podisma does not fit a divergence in isolation model (Figures 3 and 4), other factors may be involved in the maintenance of incipient divergence. Both the differences in the relative timing of divergence, as well as differences in inferred gene flow among Podisma species and population lineages (Figures 2-4), point to varying degrees in the permeability of lineage boundaries. This suggests that the differing consequences of hybridization (i.e., the varied possibilities and effects of gene flow) across the landscape may control diversification dynamics. That is, opportunities and extent of gene flow across the landscape may determine whether repeated distributional shifts act as a species-pump vs. a melting pot, not just the rate of reproductive isolation (Aguilée et al., 2018;Dynesius & Jansson, 2014;Harvey et al., 2019).
From this perspective of divergence with gene flow (as opposed to divergence in isolation), in which the rate of reproductive isolation is not the only factor controlling diversification, the timing of divergence can take on new meaning and provide new insights into the speciation process. It is notable that species boundaries are maintained even though some have remained semipermeable for extended periods of time (i.e., introgression between nonsister species P. pedestris and P. carpetana; Figure 4). In other words, reproductive isolation may be viewed as having been more or less effective in reducing the loss of incipient divergences, especially given that projections of species distributions predict some overlap among all Podisma species during glacial periods. In addition, because these past distributions were not contiguous, but rather some dispersal corridors were more limited in geographical scope than others ( Figure 6), it suggests that the opportunities for gene flow (or, conversely, the extent of geographical isolation) may have had a prominent effect on the permeability of species boundaries, in addition to any role the rate of reproductive isolation might have played in speciation. In fact, while we cannot exclude the possibility that the rate of evolution of reproductive isolation as a key determinant of the likelihood of speciation (as well as the genetic cohesion of all the species), we can make an argument that Pleistocene speciation in Podisma grasshoppers probably had other contributing factors. For example, if speciation was indeed promoted by the fragmentation and isolation of populations during the relatively short interglacial periods, rather than the geologically longer glacial periods when the grasshopper distributions are projected to have been more widespread ( Figure 6), it would imply the development of reproductive isolating mechanisms correspondingly much more rapid than in classical models where displacements into isolated glacial refugia promoted speciation (e.g., Hewitt, 1996;Knowles, 2001; also see Klicka & Zink, 1997;Ebdon et al., 2021 for arguments against Pleistocene speciation because of the rapidity of glacial cycles).
F I G U R E 6 Current and Last Glacial Maximum (LGM) distributions for each species as predicted by environmental niche models (ENMs). Colours indicate areas predicted to be occupied by each species according to the maximum training sensitivity plus specificity (MTSS) logistic threshold of their respective ENM (Table S6). Predicted distributions for the LGM are based on the MIROC-ESM and the CCSM4 general atmospheric circulation models. Elevation is shown by grey shading, with darker areas corresponding to higher elevations. Yellow colour in current distribution maps indicate small areas (barely visible) predicted as suitable by ENMs but located outside the known distribution ranges of each species (i.e., over-predictions). Current distribution maps show sampling localities (black dots with white rings) for each species; the small map inset shows the position of the species distribution on the Iberian Peninsula

| Determinants of permeable species boundaries
Although ENMs predicted that the distributions of the three taxa largely overlapped during the LGM, genomic data only supported gene flow from P. pedestris and P. carpetana. Different reasons could explain this specific history of limited introgressive hybridization.
Lack of historical hybridization between some species pairs might be a consequence of limited past connectivity and opportunity for gene flow due to geographical isolation (i.e., the geographically distant contemporary ranges of P. pedestris and P. cantabricae; Figure 1).
Alternatively, some of the studied species might have quickly evolved pre-or postzygotic reproductive isolation mechanisms (i.e., speciation in strict isolation) as a by-product of allopatric divergence (Coyne & Orr, 2004) or via reinforcement after secondary contact (Pfennig, 2016;Servedio & Noor, 2003). For instance, the two sister species P. carpetana and P. cantabricae currently have geographically adjacent distributions in northwestern Iberia (Figure 1), and multiple instances of secondary contact during the estimated distributional expansions ( Figure 6) might have accelerated the rapid evolution of reproductive isolation (Coyne & Orr, 2004;Hoskin et al., 2005).
Accordingly, previous studies have found that hybrid dysfunction, genetic incompatibilities and the evolution of pre-and postzygotic reproductive isolation mechanisms is frequent in contact zones between species (Bailey et al., 2004), subspecies (Virdee & Hewitt, 1994) and chromosomal races (Barton & Hewitt, 1981) of montane and alpine grasshoppers.
An intriguing finding is that historical gene flow between P. pedestris and P. carpetana was asymmetric (Figures 2-4). Unidirectional introgression might have resulted from extensive hybridization during periods of secondary contact followed by repeated backcrossing between hybrids and only one parental species (e.g., Field et al., 2011;Kirschel et al., 2020). Asymmetric gene flow could also be explained by a higher capacity of the donor species to disperse into the range of the recipient one (Jacquemyn et al., 2012;Ortego et al., 2021). In the absence of reproductive barriers, the two species will interbreed and first-generation hybrids will more often mate with the most abundant local species, resulting in introgressive hybridization. Although P. pedestris is generally a micropterus flightless species, long-winged individuals have been frequently described in the literature (Lemonnier-Darcemont & Darcemont, 2014 and references therein) and this polymorphism has been suggested to favour population connectivity and contribute to the colonization of suitable habitats (Zuna-Kratky et al., 2016). In contrast, the two other Podisma species from the Iberian Peninsula are either apterous (P. cantabricae) or micropterus (P. carpetana) and long-winged forms have never been reported (Morales-Agacino, 1951;Presa et al., 2016aPresa et al., , 2016b. In Orthoptera species presenting wing polymorphism, macropterous forms seem to be occasional and occur at low frequencies within populations. However, these forms have been found to be integral for range expansions at extremely short spatiotemporal scales (Hochkirch & Damerau, 2009), which might be particularly exacerbated under changing environmental conditions such as those imposed by Quaternary climatic oscillations (Simmons & Thomas, 2004). In the context of our study, increased availability of suitable habitats for colonization in the transition from interglacial to glacial periods might have led to selection for macropterous forms in peripheral populations at the expanding margins and favoured dispersal of P. pedestris into the distribution range of P. carpetana (Hochkirch & Damerau, 2009;Noguerales et al., 2016).

| CON CLUS IONS
Our integrative analyses have provided limited evidence of interspecific gene flow during prolonged periods of projected extensive secondary contact and emphasize the genetic cohesiveness of all species within the alpine Podisma complex. These findings support the notion that the interplay among Pleistocene-driven isolation (i.e., confinement in interglacial refugia), landscape composition (i.e., spatial configuration of sky islands), and species' traits (i.e., flightlessness) can trigger the necessary mechanisms for long-lasting genomic diversification and speciation in alpine and montane biotas (Dynesius & Jansson, 2014;Knowles, 2001).
Our comprehensive suite of distributional, demographic and phylogenomic analyses also provide a mechanistic explanation for the uncertain phylogenetic relationships among the studied grasshopper species and collectively highlight the important role of Quaternary climatic oscillations in promoting diversification and genetic fragmentation of relictual alpine organisms from temperate regions that currently persist as highly isolated populations in disparate mountain ranges. Irrespective of which factors control diversification (i.e., the rate of reproductive isolation vs. dispersal, and hence opportunities for gene flow), the time of speciation supports a model of Pleistocene divergence. This, in its own right, means that Podisma grasshoppers of the sky islands from the Iberian Peninsula constitute an ideal system to further investigate some intriguing questions that the different independent data types and analytical procedures raise about the speciation process. Future experimental crosses might reveal the presence of pre-and postzygotic barriers to interspecific gene flow that could clarify whether the absence of genomic evidence for introgressive hybridization among most species-pairs is a consequence of reproductive isolation and the completion of the speciation process or resulted from limited opportunity for hybridization due to population isolation in sky islands and dispersal constraints during range expansions.

ACK N OWLED G EM ENTS
We thank Amparo Hidalgo-Galiana for library preparation, Víctor

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