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

AIM: 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. LOCATION: South‐eastern Iberia. METHODS: We quantified genomic variation in the two focal taxa and coupled ecological niche models and a spatiotemporally explicit simulation approach based on coalescent theory to determine the relative statistical support of a suite of competing demographic scenarios representing contemporary population isolation (i.e. a predominant role of genetic drift) versus historical connectivity and post‐glacial colonization of sky islands (i.e. pulses of gene flow and genetic drift linked to Pleistocene glacial cycles). RESULTS: Inference of spatial patterns of genetic structure, environmental niche modelling and statistical evaluation of alternative species‐specific demographic models within an approximate Bayesian computation framework collectively supported genetic admixture during glacial periods and post‐glacial colonization of sky islands, rather than long‐term population isolation, as the scenario best explaining the current distribution of genomic variation in the two focal taxa. Moreover, our analyses revealed that isolation in sky islands has also led to extraordinary genetic fragmentation and contributed to reduce local levels of genetic diversity. MAIN CONCLUSIONS: This study exemplifies the potential of integrating genomic and environmental niche modelling data across biological and spatial replicates to determine whether organisms with similar habitat requirements have experienced 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.