Glacial connectivity and current population fragmentation in sky islands explain the contemporary distribution of genomic variation in two narrow‐endemic montane grasshoppers from a biodiversity hotspot

Cold‐adapted biotas from mid‐latitudes often show small population sizes, harbour low levels of local genetic diversity and are highly vulnerable to extinction due to ongoing climate warming and the progressive shrinking of montane and alpine ecosystems. In this study, we use a suite of analytical approaches to infer the demographic processes that have shaped contemporary patterns of genomic variation in Omocestus bolivari and Omocestus femoralis, two narrow‐endemic and red‐listed Iberian grasshoppers forming highly fragmented populations in the sky island archipelago of the Baetic System.


| INTRODUC TI ON
Distributional shifts in response to Quaternary climatic oscillations had a dramatic impact on biogeographical patterns of species diversity, abundance and local endemism (Hewitt, 1996;Sandel et al., 2011). These cycles, in particular the high-amplitude climatic changes characterizing the mid and late Pleistocene (<0.8 Ma ;Jouzel et al., 2007), have also shaped the distribution and spatial patterns of genetic variation in many organisms (Hewitt, 2000). However, multiple studies have documented considerable heterogeneity across regions and taxa in the demographic consequences of Pleistocene glacial cycles (Hewitt, 2000). On the one hand, the impact of Pleistocene glaciations strongly depended on latitude and regional topography, with extinction-recolonization dynamics at higher latitudes (i.e. "southern richness to northern purity" paradigm; Hewitt, 1996) and elevational shifts and more complex processes of population fragmentation and connectivity at lower latitudes such as the tropics or temperate regions (e.g. "refugia within refugia" concept; Gómez & Lunt, 2007). On the other hand, the way organisms respond to climate changes strongly depends on species-specific niche requirements and life-history traits, which define favourable/ unfavourable climatic periods (glacials or interglacials; Bennett & Provan, 2008) and their capacity to deal with population fragmentation (e.g. microhabitat preferences, dispersal capacity; Massatti Papadopoulou & Knowles, 2016). These aspects determined the location and extension of Pleistocene refugia, which have played a predominant role on species' persistence during unfavourable climatic periods and acted as source populations from which species expanded their ranges at the onset of more favourable conditions (Bennett & Provan, 2008). Temperate species currently inhabiting low-elevation areas generally restricted their distributions to southern refugia during glacial phases (i.e. glacial refugia) and expanded during interglacial periods whereas cold-adapted species, nowadays presenting fragmented populations at high elevations/latitudes (i.e. interglacial refugia), had much more widespread distributions in glacial stages (Hewitt, 2000).
Mid-and low-latitude alpine and montane taxa represent smallscale replicates of cold-adapted species living at high latitudes (Hewitt, 1996(Hewitt, , 2000. Sky islands are a paradigmatic example of interglacial refugia in which alpine and montane organisms currently form highly isolated populations in mountain ranges separated from each other by intervening valleys with unsuitable environmental conditions (e.g. DeChaine & Martin, 2005; Knowles & Massatti, 2017;McCormack et al., 2008;Mouret et al., 2011). Changing climatic conditions causes suitable habitats to expand or shrink along elevational gradients, leading to alternating periods of population isolation and connectivity (DeChaine & Martin, 2005;Knowles, 2000). Periods of isolation on mountain tops are expected to lead to demographic bottlenecks, genetic drift and divergence among populations, while periods of connectivity allow for dispersal and gene flow among formerly isolated populations (DeChaine & Martin, 2005). Species restricted to sky islands often show unique patterns of population genetic structure resulting from climate-induced distributional shifts (Mouret et al., 2011), processes that in some cases have led to lineage diversification and speciation (Knowles, 2000;Knowles & Massatti, 2017). Therefore, montane habitats are often important hotspots of intraspecific genetic diversity, species richness and local endemism, providing an excellent study system for understanding how climatic changes associated with Pleistocene glacial cycles impacted species distributions, demography and genetic diversification (Mairal et al., 2017). Their study takes particular relevance if we consider that species restricted to sky islands often show small population sizes, harbour reduced levels of local genetic diversity and are highly vulnerable to extinction due to ongoing climate warming and the progressive shrink of alpine ecosystems (Rubidge et al., 2012).
South-eastern Iberia constitutes an important biodiversity hotspot due to the interplay among a complex geological history (Meulenkamp & Sissingh, 2003), an extraordinarily complex topography (Braga et al., 2003) and a limited impact of Quaternary glaciations (Hughes & Woodward, 2017). The presence of vast areas free of permanent ice during the coldest stages of the Pleistocene made this region an important glacial refugium for warm-temperate taxa (Gómez & Lunt, 2007;Hewitt, 2011). At the same time, the mountain ranges of the region (Baetic system; >3,000 m) constitute important interglacial refugia for several cold-adapted organisms that currently persist in severely fragmented populations at high elevation patches of suitable habitat and contribute disproportionally to the extraordinary rates of local endemism of the area (Mota et al., 2002). This offers an ideal biogeographical setting to jointly study the genetic consequences of long-term population fragmentation in this sky island system and assess the time-scale at which past distributional shifts have shaped contemporary patterns of genetic variation in alpine organisms from mid-latitudes.
Here, we use genomic data obtained via restriction siteassociated DNA sequencing (ddRAD-seq; Peterson et al., 2012) and a suite of analytical approaches to shed light on the demographic concerted/idiosyncratic responses to Quaternary climatic oscillations, which can ultimately help to reach more general conclusions about the vulnerability of mountain biodiversity hotspots to ongoing climate warming.

K E Y W O R D S
approximate Bayesian computation, ddRAD-seq, demographic history, environmental niche modelling, landscape genetics, Pleistocene glaciations processes underlying contemporary patterns of genetic variation in the grasshoppers Omocestus bolivari Chopard, 1939, andOmocestus femoralis Bolívar, 1908, two narrow-endemic Iberian species inhabiting high elevation areas in the sky island archipelago of the Baetic System (Presa et al., 2016a(Presa et al., , 2016b. The two taxa are closely related non-sister species that present adjacent but not overlapping distributions (Tonzo et al., 2021; see Figure 1), both occupying similar open habitats of thorny shrub formations (e.g. Erinacea sp., Festuca sp., Juniperus sp.) located at the upper (>1,500 m; see Table S1) thermoclimatic belts of the Mediterranean region (Clemente et al., 1991;Presa et al., 2016aPresa et al., , 2016b. They are graminivorous, univoltine (i.e. a single generation per year) and present similar phenologies, with adult population peaks in July (Clemente et al., 1991). The distribution range of the two species includes sky islands of varying size, ranging from large ones in the main massifs of the region (O. bolivari: Sierra Nevada and Sierra de Baza-Filabres; O. femoralis: Sierra de Cazorla) to tiny patches of suitable habitat located in isolated mountains with maximum elevations close to the lower altitudinal limit of the species (Figure 1). By using these two species of uniform ecological and lifehistory traits as independent replicates, we exemplify how the integration of genomic and environmental niche modelling (ENM) data within a comparative framework can help to elucidate the responses of species with similar habitat requirements to processes of population fragmentation/connectivity driven by Pleistocene climatic oscillations. This can ultimately help to understand the temporal scale at which such processes have shaped contemporary patterns of genetic variation in the focal species and other co-distributed organisms with analogous environmental niches (Knowles & Massatti, 2017).
In a first step, (a) we tested whether observed patterns of genomic variation in the two biological replicates reflect the signals of historical population connectivity (i.e. gene flow during glacial periods) and subsequent colonization and population isolation in sky islands (i.e. genetic drift during interglacial periods) or, rather, if post-glacial population fragmentation and genetic drift after the Pleistocene-Holocene transition have entirely monopolized the genetic makeup of contemporary populations. Specifically, we quantified genetic structure and coupled ecological niche models (ENMs) and a spatiotemporally explicit simulation approach based on coalescent theory to mimic complex demographic processes under a suite of scenarios representing historical gene flow and post-glacial colonization in sky islands versus contemporary population isolation (Knowles & Massatti, 2017;Ray et al., 2010). The expectations under these alternative scenarios were tested and validated against observed genomic data within an approximate Bayesian computation (ABC) framework (e.g. He et al., 2013). In a second step, (b) we analysed the impact of sky island size and the spatial location of contemporary populations on their demographic fate and levels of genetic variation. To this end, we tested the hypothesis that genetic diversity is positively associated with the geographic centrality of populations and habitat patch size and used demographic reconstructions based on genomic data to determine whether changes in effective population size (N e ) over time differed between sky islands of contrasting size and degree of peripheral location (Lira-Noriega & Manthey, 2014).  Table S1). The two species present allopatric distributions and can be easily identified based on diagnostic morphological traits (i.e. shape and relative size of forewings and hindwings; Clemente et al., 1991). According to our extensive surveys in south-eastern Iberia and available records in the literature (Presa et al., 2016a(Presa et al., , 2016b, the sampled populations cover the entire known distribution range of the two species (Figure 1). We stored specimens in 2-ml vials with 96% ethanol and preserved them at −20°C until needed for genomic analyses. All specimens were preserved as vouchers in our Orthoptera collections, with labels for genotyped individuals available at the NCBI Sequence Read Archive

| Genomic library preparation and processing
We used Nucleo Spin Tissue kits (Macherey-Nagel, Düren, Germany) to extract and purify total genomic DNA from a hind leg of each individual. Genomic DNA was individually barcoded and processed in-house using the double-digest restriction site-associated DNA procedure (ddRAD-seq) described in Peterson et al. (2012). Briefly, DNA was doubly digested with EcoRI and MseI restriction enzymes, followed by the ligation of Illumina adaptor sequences and unique 7base pair barcodes. Ligation products were pooled into four libraries, size-selected for fragments between 475 and 580 bp using a Pippin prep machine (Sage Science) and amplified by iProof™ High-Fidelity DNA Polymerase (Bio-Rad) with 12 cycles. Single-read 151-bp sequencing was performed on an Illumina HiSeq2500 platform at The Centre for Applied Genomics (Hospital for Sick Children). We used the different programs distributed as part of the stacks v. 1.35 pipeline (Catchen et al., 2013) to filter and assemble our sequences into de novo loci and call genotypes. All details on sequence assembling and data filtering are presented in Methods S1.

| Population genetic structure
We analysed population genetic structure of the two species using the Bayesian Markov chain Monte Carlo clustering method implemented in the program structure v. 2.3.3 (Pritchard et al., 2000). We LGM ran structure assuming correlated allele frequencies and admixture without using prior population information. We conducted 15 independent runs with 200,000 MCMC cycles, following a burn-in step of 100,000 iterations, for each of the different possible K genetic clusters. We retained the ten runs having the highest likelihood for each value of K and inferred the number of populations best fitting the dataset using log probabilities [Pr(X|K)] (Pritchard et al., 2000) and the ΔK method (Evanno et al., 2005). Complementary to structure analyses, we performed principal component analyses (PCA) as implemented in the r package adegenet (Jombart, 2008).

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

| Testing alternative demographic models
We used the integrative distributional, demographic and coalescent (iDDC) approach (He et al., 2013) and an approximate Bayesian computation (ABC) framework (Beaumont et al., 2002) to test alternative scenarios representing different hypotheses about how landscape heterogeneity (or its lack thereof) and colonization from glacial refugia versus contemporary population isolation in sky islands explain the spatial distribution of genomic variation in our two focal species. This approach is described in He et al. (2013) and consists of three main steps: (a) constructing alternative demographic models representing different hypotheses about the processes shaping spatial patterns of genetic structure and diversity;

| Constructing demographic models
We generated two main sets of models that differ in the hypothetical demographic processes that have shaped spatial patterns of contemporary genetic variation (Table 1) ). Finally, we generated two F I G U R E 1 (a, b) Spatial patterns of genetic diversity and (c-h) environmental niche models (ENM) for Omocestus bolivari and O. femoralis. Panels a-b show genetic diversity (π) of the studied populations (triangles, population codes as in Table S1) interpolated across the current distribution of each species as predicted by their respective ENMs. Panels C-H show the projection of ENMs for (c-d) present and Last Glacial Maximum (LGM) conditions under the (e-f) CCSM4 and (g-h) MIROC-ESM general atmospheric circulation models. Maps for current distributions show records (crosses) used to build the ENMs and the dark yellow colour indicates areas predicted as suitable out of the known limits of species' ranges (i.e. overpredictions, not shown in panels a-b). All maps show suitable areas above the maximum training sensitivity plus specificity (MTSS) logistic threshold of maxent. Grey background represents elevation, with darker areas corresponding to higher altitudes. Latitude and longitude in decimal degrees. The temporal scale of "isolation" (static) and "colonization" (dynamic) demographic models is

| Model choice and parameter estimation
We used an approximate Bayesian computation (ABC) framework to perform model selection and parameter estimation, as implemented in abctoolbox programs (transformer and abcestimator) and r scripts (findPLS) . We used the r package pls v.2.6-0 (Mevik & Wehrens, 2007) and the findPLS script to extract partial least-squares (PLS) components with Box-Cox transformation from the summary statistics of the first 10,000 simulations for each model . The first five PLSs extracted from the summary statistics were used for ABC analyses, as the root-mean-squared error (RMSE) of the three demographic parameters employed (K MAX , m, N ANC ; see Methods S3 for details) for the two species did not decrease significantly with additional PLSs ( Figure S1). We used the linear combinations of summary statistics obtained from the first 10,000 simulations for each model to transform all datasets (observed empirical and simulated datasets) with the program transformer . For each model, we retained the 1,000 simulations (0.5%) closest to observed empirical data and used them to approximate marginal densities and posterior distributions of the parameters with a post-sampling regression adjustment using the ABC-GLM (general linear model) procedure detailed in Leuenberger and Wegmann (2010) and implemented in abcestimator. We used Bayes factors (BF) for model selection (Kass & Raftery, 1995). Note: A higher marginal density corresponds to a higher model support and a high Wegmann's p-value indicates that the model is able to generate data in agreement with empirical data. Bayes factors represent the degree of relative support for the model with the highest marginal density over the other models. R 2 is the coefficient of determination from a regression between each demographic parameter (K MAX , m, N ANC ) and the five partial least squares (PLS) extracted from all summary statistics.
K MAX , carrying capacity of the deme with highest suitability (100%).
m, migration rate per deme per generation.
N ANC , ancestral population size.

| Model validation
To evaluate the ability of each model to generate the empirical data, we calculated Wegmann's p-value from the 1,000 retained simulations . We also assessed the potential for a parameter to be correctly estimated by computing the proportion of parameter variance that was explained (i.e. the coefficient of determination, R 2 ) by the retained PLSs. For the most supported model for each species, we determined the accuracy of parameter estimation using a total of 1,000 pseudo-observation datasets (PODs) generated from prior distributions of the parameters. If the estimation of the parameters is unbiased, posterior quantiles of the parameters obtained from PODs should be uniformly distributed . As with the empirical data, we calculated the posterior quantiles of true parameters for each pseudo run based on the posterior distribution of the regression-adjusted 1,000 simulations closest to each pseudo-observation.

| Genetic diversity and historical changes in effective population size
To further explore the processes determining the spatial distribution of genetic variation in the two species, we (a) tested the association between genetic diversity and the geographic centrality and habitat patch size of the studied populations and (b) used genomic data to reconstruct changes in effective population size (N e ) over time.
We employed linear regressions in spss v. 26.0 to analyse whether genetic diversity (nucleotide diversity, π) is associated with the distance of each population to the species' distribution centroid and contemporary suitable habitat patch size. Given that the precision of genetic diversity estimates may differ among populations due to differences in sample sizes, we used a weighted least-square (WLS) method where weight equals the number of genotyped individuals for each population (Table S1). We calculated nucleotide diversity for each population using the program populations from stacks. The centroid of species distribution was calculated in arcmap v.10.3 on the basis of polygons (patches) of suitable habitat defined by grid cells with a probability of presence of the focal species above the maximum training sensitivity plus specificity (MTSS) logistic threshold for occurrence from maxent (Liu et al., 2005). Similarly, patch size for each study population was calculated considering the area (in km 2 ) of continuous suitable habitat defined as above on the basis of raster maps obtained from projecting ENMs to contemporary bioclimatic conditions.
We inferred changes in effective population sizes (N e ) over time for each studied population using stairway plot (Liu & Fu, 2015 (Figures 1b and 2b). Principal component analyses for both species yielded analogous results to those obtained with structure, separating along PC1 north-western from south-eastern populations and grouping populations according to their geographical location ( Figure S5).

| Environmental niche modelling
We generated final ENMs using the feature class (FC) combination

| Testing alternative demographic models
Based on marginal densities calculated from the 0.5% of retained simulations, the best ranked models in the two species were those incorporating the colonization process from hypothetical glacial refugia (Models A, B, D and E; Table 1). In contrast, models of isolation in contemporary sky islands had considerably lower marginal densities and a difference in Bayes factors > 10 6 with the most supported model (Table 1), indicating strong relative support for colonization models in the two species (Kass & Raftery, 1995 (Kass & Raftery, 1995). In O. femoralis, both dynamic and static (i.e. flat landscape) colonization models based on CCSM projections (Models B and E) were the most supported by our empirical data and statistically indistinguishable (BF < 2; Table 1). The dynamic colonization model based on MIROC bioclimatic conditions also had a low Bayes factor (<20), but its low Wegmann's p-value (<.03) indicates that it was not capable of generating data compatible with our empirical data (Table 1).
Posterior distributions of parameters under the most probable models were considerably distinct from the prior, indicating that the simulated data contained information relevant to estimating the parameters ( Figure S6). Comparison of the posterior distributions before and after the ABC-GLM also showed the improvement that this procedure had on parameter estimates ( Figure S6). The posterior distribution of K MAX in O. femoralis was flat, indicating uncertainty in the estimation of this parameter ( Figure S6b). Accordingly, the coefficients of determination (R 2 ) from a multiple regression between each demographic parameter and the five retained PLS indicated that the employed summary statistics had a moderately high potential to correctly estimate all the parameters except K MAX in O. femoralis (Table 1).
Posterior quantiles of m in O. bolivari, K MAX in O. femoralis and N ANC in both taxa significantly deviated from a uniform distribution, indicating a potential bias in the estimation of these parameters ( Figure S7).

| Genetic diversity and historical changes in effective population sizes
Visual inspection of the spatial distribution of genetic variation in the two species evidenced that peripheral populations and those located in smaller sky islands tended to present lower levels of genetic diversity (Figure 1a,b). Accordingly, genetic diversity (π) in populations of O. bolivari was negatively associated with population centrality (r = −.913, t = −6.72, p < .001; Figure 3a) and showed a positive relationship with patch size (r = 0.733, t = 3.23, p = .010; Figure 3b).
In the case of O. femoralis, genetic diversity was positively associated with patch size (r = 0.878, t = 3.17, p = .050; Figure 3d) but the negative relationship with population centrality did not reach statistical significance (r = 0.743, t = −1.92, p = .150; Figure 3c). stairway plot analyses for O. bolivari revealed similar results across most of the populations, with an abrupt demographic expansion led by the advance of the glaciation till its maximum level around 21 ka BP (Figure 4a,b). Afterwards, populations experienced a decline in N e until recent times (Figure 4a,b). One exception was the highly isolated populations PAND, MAGI and LUJA, which presented a general demographic stability during the last 100 ka (Figure 4b).
In the case of O. femoralis, most populations showed idiosyncratic demographic histories, which ranged from expansions during the last glacial period followed by reductions in N e at the onset of the Holocene (ALME and MARI) to marked bottlenecks around the LGM (ESPU and CAZO) (Figure 4c,d).

| D ISCUSS I ON
Our study contributed to a better understanding of the demographic processes underlying spatial patterns of genetic variation in coldadapted biotas currently forming severely fragmented populations in sky island archipelagos from temperate regions (Hewitt, 1996;Knowles & Massatti, 2017;Perrigo et al., 2020). Specifically, the combination of population genomic data, niche modelling and spatiotemporally explicit coalescent-based simulations supported population connectivity during glacial periods followed by colonization and population isolation in sky islands as the most supported sce- Ultimately, this can lead to more general conclusions about the way organisms and whole communities have responded to Quaternary climatic oscillations. These aspects are of pivotal importance to understand the dynamics of montane and alpine ecosystems from midand low-latitude regions, which not only harbour high levels of local endemism but have been also identified among the most vulnerable to extinction due to ongoing climate warming and the progressive shrink of these habitats (Perrigo et al., 2020).

| Marked genetic structure in sky islands
The GeneƟc diversity (π) (d) transcend their contemporary distributions in specific mountain ranges, which warns against oversimplified isolation models often assumed for alpine and montane organisms (see also Knowles & Alvarado-Serrano, 2010;Knowles & Massatti, 2017). Thus, the idiosyncratic nature of isolation of sky island biotas cannot be fully understood without a thoughtful assessment of the spatial F I G U R E 4 Demographic history of each studied population of (a, b) Omocestus bolivari and (c, d) Omocestus femoralis inferred using stairway plot. Lines show the median estimate of effective population size (N e ) over time (ka) for populations located in (a, c) large and (b, d) small sky islands. Population codes shown in the legend of each panel are described in Table S1  and temporal dimensions of population fragmentation and connectivity (i.e. colonization history, multiple source populations and historical gene flow), which present certain analogies but are also expected to substantially differ from the evolutionary pathways and demographic dynamics dominating true island systems (Flantua et al., 2020).

| Post-glacial colonization of sky islands
Model selection under an ABC framework revealed that the spatially explicit demographic scenarios most supported by empirical data were those in which the spatial distribution of contemporary genomic variation was shaped by processes of post-glacial colonization of sky islands from source populations defined by the spatial location of environmentally suitable areas during the last ice age (Knowles & Massatti, 2017;Mouret et al., 2011). In contrast, models of isolation in sky islands since the onset of the Holocene had a very low statistical support, indicating that genetic drift after post-glacial fragmentation cannot solely explain the genetic makeup of contemporary populations. Although scenarios based on colonization of sky islands were the only ones able to simulate genomic data compatible with observed data, our approach was not able to distinguish between dynamic models considering changes in landscape heterogeneity since the LGM and static models considering a diffusion scenario of expansion from glacial-source populations (see also Knowles & Massatti, 2017). These results are surprising considering that the two focal taxa have strict habitat preferences for high-elevation areas and unsuitable lowlands are expected to act as impassable barriers to dispersal, as evidenced by distribution models  (Slatkin, 1993), which is similar to the expected outcome of migration-genetic drift equilibrium under a static, flat landscape scenario (He et al., 2013;Knowles & Massatti, 2017). Thus, admixture among nearby populations is likely to have reduced the power of our analyses to discriminate between dynamic (Models A-B) and static (Models D-E) scenarios of postglacial colonization, particularly in O. femoralis for which only five populations were available for analyses (Table S1).
The two general circulation palaeoclimatic models supported the expansion of the two focal taxa towards lowland areas during the LGM, but their respective populations were predicted to be much more connected under the MIROC-ESM than under the CCSM model ( Figure 1). Although differences between palaeoclimatic models in predicting past species distributions are a frequent outcome of most studies (e.g. Ramírez-Barahona & Eguiarte, 2014), the relative support for contrasting predictions is rarely validated using independent sources of information (Nogués-Bravo, 2009). Such differences seemed to be captured by our model-based approach, which revealed that glacial-source populations identified for the two focal taxa on the basis of projections of ENMs to bioclimatic conditions under the CCSM general circulation model tended to provide a better fit to our genomic data than those obtained based on MIROC projections (BF > 2 in all cases; Table 1). Although differences in the relative likelihood for scenarios based on the two palaeoclimatic models are small and must be interpreted with extreme caution, our results highlight the potential of spatiotemporally explicit coalescent-based simulations to refine hypotheses about the location of glacial refugial (Bemmels et al., 2019) and validate distributional shifts inferred from ENMs that cannot usually be contrasted with independent sources of information (e.g. fossil records; Fordham et al., 2014;Metcalf et al., 2014).

| Genetic diversity and demographic reconstructions
The spatial distribution of genetic variation in the two species evidenced that peripheral populations and those located in smaller sky islands tend to present comparatively lower levels of genetic diversity ( Figure 1). These findings are consistent with the centre-periphery model, which predicts that populations located on the species range edges become progressively smaller, are spatially more isolated, and harbour lower levels of genetic diversity than populations situated at the core of the distribution (e.g. Lira-Noriega & Manthey, 2014).
In contrast, genomic-based demographic reconstructions indicate that, in general terms, most populations presented larger effective population sizes during the last glacial period and experienced demographic declines at the onset of the Holocene (Figure 4), in agreement with palaeoclimatic-based reconstructions of species distributions ( Figure 1). Inference of changes of N e through time did not reveal obvious differences between populations located in large versus small sky islands, suggesting that contemporary populations, and their ancestral sources, reacted in a similar fashion to Quaternary climatic oscillations independently of their current geographical location. This is not surprising given the small range of the two focal taxa (Presa et al., 2016a(Presa et al., , 2016b, with most distant populations separated by <150 km, and the likely similar regional impact of Pleistocene glacial cycles across the entire distribution of the species. Also, it must be considered that, according to our own field observations, small habitat patches can often sustain high local densities of the species. This might contribute to maintain large effective population sizes and avoid strong genetic drift even in populations currently confined to tiny sky islands. These results indicate that although the genetic makeup of contemporary populations has been shaped in a great extent by historical processes of genetic admixture (glacial periods) and drift during the retreat to interglacial refugia (e.g. Knowles & Massatti, 2017), contemporary isolation has also contributed to erode the levels of genetic diversity of populations persisting in peripheral sky islands of small size (Lira-Noriega & Manthey, 2014; Rubidge et al., 2012). Collectively, our results are in agreement with the "interglacial refugia model" documented for other montane and alpine organisms from temperate and tropical regions in which populations experienced expansions during glacial periods followed by severe demographic contractions and disruption of gene flow during interglacials (Bennett & Provan, 2008;Camacho-Sánchez et al., 2018). It must be noted that this scenario is very different to the dynamics of alpine organisms from temperate regions that were more extensively glaciated (e.g. North America) and for which genetic-based inferences have often supported strong isolation and divergence in glacial refugia and net demographic expansions during warm interglacials ("expanding alpine archipelago model"; e.g. Maier et al., 2019;Schoville et al., 2011). This emphasizes the importance of inferring the demographic fate of montane and alpine biotas under different biogeographical contexts, as these will largely determine, in interaction with taxon-specific traits (e.g. dispersal capacity and habitat specialization), patterns of local endemism, the extent of population fragmentation, and the magnitude of distributional shifts in response to past and future climate changes (Flantua et al., 2020;Massatti & Knowles, 2016).

| CON CLUS IONS
By considering two biological replicates from the sky island archipelago of the Baetic System, our study highlights the potential of integrating different sources of information to infer the evolutionary and demographic processes shaping spatial patterns of genetic variation in cold-adapted faunas from temperate regions. Our spatiotemporally explicit simulations and testing of alternative demographic scenarios supported population expansion during cold periods and subsequent post-glacial colonization and isolation in sky islands as the most likely explanation for the current distribution of genetic variation in alpine and montane biotas from the Mediterranean biodiversity hotspot. Global warming is expected to threaten the persistence of alpine and cold-adapted organisms from mid-and low-latitude regions by extraordinarily reducing the extent of suitable habitats already confined to mountain tops (Rubidge et al., 2012). The similar demographic responses inferred for the two focal taxa to past climate changes suggest that they and other co-distributed organisms with similar life-history traits will probably present concerted responses to ongoing global change. Thus, our results can help to establish unified conservation strategies aimed at preserving the extraordinary rates of local endemism of this and other mountain biodiversity hotspots, from a community-level perspective (Perrigo et al., 2020).

ACK N OWLED G EM ENTS
We are much indebted to Anna Papadopoulou for her valuable help in study design and useful comments, suggestions and corrections on a first draft of the manuscript. We also wish to thank Pedro J. Cordero and Víctor Noguerales for their help during fieldwork and providing samples from some populations, Amparo Hidalgo for support during laboratory work, Sergio Pereira (The Centre for Applied Genomics) for Illumina sequencing and three anonymous referees for their constructive and valuable comments on an earlier version of the manuscript. Logistical support was provided by Laboratorio de Ecología Molecular (LEM-EBD) and Laboratorio de Sistemas de Información Geográfica y Teledetección (LAST-EBD) from Estación Biológica de Doñana. We also thank Centro de Supercomputación de Galicia (CESGA) and Doñana's Singular Scientific-Technical Infrastructure (ICTS-RBD) for access to computer resources. This study was funded by the Spanish Ministry of Economy and Competitiveness and the European Regional Development Fund (ERDF) (CGL2014-54671-P and CGL2017-83433-P). VT was supported by an FPI predoctoral fellowship (BES-2015-73159) from the Spanish Ministry of Economy and Competitiveness.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13306.