Cold-adapted species escape climate warming by latitudinal and/or altitudinal range shifts, and currently occur in Southern Europe in isolated mountain ranges within ‘sky islands.

Here we studied the genetic structure of the Apollo butterfly in five such alpine islands (above 1000 m) in France, and infer its demographic history since the last interglacial, using single nucleotide polymorphisms (ddRADseq SNPs). The Auvergne and Alps populations show strong genetic differentiation but not alpine massifs, although separated by deep valleys. Combining three complementary demographic inference methods and species distribution models (SDMs) we show that the LIG period was highly defavorable for Apollo that probably survived in small population in the highest summits of Auvergne. The population shifted downslope and expanded eastward between LIG and LGM throughout the large climatically suitable Rhône valley between the glaciated summits of Auvergne and Alps. The Auvergne and Alps populations started diverging before the LGM but remained largely connected till the mid-Holocene. Population decline in Auvergne was more gradual but started before (~7 kya versus 800 ya), and was much stronger with current population size ten times lower than in the Alps. In the Alps, the low genetic structure and limited evidence for isolation by distance suggest a non-equilibrium metapopulation functioning. The core Apollo population experienced cycles of contraction-expansion with climate fluctuations with largely inter-connected populations over time according to a ‘metapopulation-pulsar’ functioning. This study demonstrates the power of combining demographic inferences and SDMs to determine past and future evolutionary trajectories of an endangered species at a regional scale.