Spatiotemporal diversification of forest understorey species reveals the existence of multiple Pleistocene forest refugia in Central Europe

Spatiotemporal diversification of forest understorey species reveals the existence of multiple Pleistocene forest refugia in Central Europe | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Spatiotemporal diversification of forest understorey species reveals the existence of multiple Pleistocene forest refugia in Central Europe Camille Voisin , View ORCID Profile Philipp Kirschner , Eliška Záveská , Božo Frajman , Karl Hülber , Johannes Wessely , Wolfgang Willner , Peter Schönswetter , Pau Carnicero doi: https://doi.org/10.1101/2025.05.13.653645 Camille Voisin 1 Department of Botany, University of Innsbruck , Innsbruck, Austria 2 Laboratoire d’Écologie Alpine, Université Grenoble Alpes, Université Savoie Mont Blanc , CNRS, Grenoble, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Philipp Kirschner 1 Department of Botany, University of Innsbruck , Innsbruck, Austria Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Philipp Kirschner For correspondence: philipp.kirschner{at}gmail.com Eliška Záveská 3 Institute of Botany of the Czech Academy of Sciences , Průhonice, Czechia Find this author on Google Scholar Find this author on PubMed Search for this author on this site Božo Frajman 1 Department of Botany, University of Innsbruck , Innsbruck, Austria Find this author on Google Scholar Find this author on PubMed Search for this author on this site Karl Hülber 4 Department of Botany and Biodiversity Research, University of Vienna , Vienna, Austria Find this author on Google Scholar Find this author on PubMed Search for this author on this site Johannes Wessely 4 Department of Botany and Biodiversity Research, University of Vienna , Vienna, Austria Find this author on Google Scholar Find this author on PubMed Search for this author on this site Wolfgang Willner 4 Department of Botany and Biodiversity Research, University of Vienna , Vienna, Austria 5 Vienna Institute for Nature Conservation & Analyses (VINCA) , Vienna, Austria Find this author on Google Scholar Find this author on PubMed Search for this author on this site Peter Schönswetter 1 Department of Botany, University of Innsbruck , Innsbruck, Austria Find this author on Google Scholar Find this author on PubMed Search for this author on this site Pau Carnicero 1 Department of Botany, University of Innsbruck , Innsbruck, Austria 6 Department of Animal Biology , Plant Biology and Ecology, Autonomous University of Barcelona , Bellaterra, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract During Pleistocene cold stages, European temperate forests were not only restricted to refugia in the southern European peninsulas. Rather, there is increasing evidence for survival of trees also in isolated patches further north, termed “northern refugia”. While their existence is undisputed, based on what is known from a handful of tree species, there is very limited knowledge about forest understoreys. Here, we fill this gap by examining the evolutionary histories of three Central European forest understorey species (FUS; Aposeris foetida, Cardamine trifolia, Hacquetia epipactis). To do so, we use a set of exploratory and explicit analyses utilizing genomic data and ecological niche models, and interpret these data following an a priori defined framework. We identify the northwestern Balkan Peninsula as the primary diversification center for the three species but found additional northern refugia in the Alps, the Carpathians, and the Apennines. Divergence times indicated pre-Last Glacial Maximum (LGM) diversification within each FUS, suggesting persistence of forest islands in Central Europe during Pleistocene cold stages rather than exclusively post-LGM colonization. We conclude that FUS thrived in scattered northern refugia. This refines our understanding of past forest dynamics and further supports widespread long-term persistence of forest patches in Central Europe during Pleistocene cold stages. Introduction The geographic range of temperate forests, which are currently the dominant zonal vegetation of large parts of Europe ( Leuschner & Ellenberg, 2017 ), has been strongly reshuffled by Quaternary climatic fluctuations ( Birks & Tinner, 2016 ; Magri et al., 2006 ; Petit et al., 2003 ). Phases of drastic range reductions to small and disjunct refugia during cold stages alternated with phases of expansions leading to continuous distributions during warm stages such as the Holocene ( Comes & Kadereit, 1998 ). Last Glacial Maximum (LGM, ∼ 21 ka) refugia for temperate deciduous trees were initially proposed only in the southern European peninsulas ( Hewitt, 1999 ; Taberlet et al., 1998 ). However, subsequent studies suggested additional refugia north of these peninsulas, that are often also named ‘northern refugia’ or ‘cryptic refugia’ ( Stewart & Lister, 2001 ; Postolache et al ., 2017 ; conceptual discussion in Rull, 2010 ; hereafter only the term ‘northern refugia’ will be used). Several northern refugia for temperate forests were proposed for Central Europe, the northern Apennines, the southwesternmost Alps and the southern Carpathians ( Fig. 1 ). Specifically, for European beech ( Fagus sylvatica ), one of the dominant trees of European temperate forests, macrofossil and palynological evidence alongside genetic data suggested northern refugia in the northwestern Balkan Peninsula, the southwestern fringe of the Western Alps, the Western Carpathians and the Apuseni Mountains ( Fig. 1 ; Magri et al ., 2006 ; Brus, 2010 ). Genomic data further supported a northern beech refugium in the Apennines ( Marchesini et al., 2023 ). The existence of an isolated refugium of broad-leaved trees in the Euganean Hills (Colli Euganei) at the southern margin of the Alps ( Gubler et al., 2018 ; Kaltenrieder et al., 2009 ) and in a former hydrothermal field in Moravia ( Hošek et al., 2024 ) were further suggested by pollen data as well as retrospective evaluation of current microclimatic variability. For coniferous trees, additional refugia were discovered at the southern and south-eastern edge of the Alps and in the adjacent Po Plains ( Gugerli et al., 2023 ; Monegato et al., 2015 ; van der Knaap et al., 2005 ), the forelands of the easternmost Central Alps (e.g. Huber et al ., 2010 ) and even the easterly adjacent Pannonian Plains ( Fig. 1 ; Willis & van Andel, 2004 ; Gugerli et al ., 2023 ). Download figure Open in new tab Fig. 1. Potential refugia of the studied forest understorey species in Central Europe and neighboring areas during the Last Glacial Maximum (LGM, ∼21 ka). Green areas refer to the current distribution of European beech (taken from( Willner et al., 2023 )) and reflect the potential natural distribution of temperate forest. Light blue areas indicate the maximum extent of the LGM glaciation ( Ehlers et al., 2011 ). 1 Western Alps ( Magri et al., 2006 ); 2 northern Apennines ( Guido et al., 2020 ; Marchesini et al., 2023 ; Willner et al., 2009 ); 3 Southern Limestone Alps ( Brus, 2010 ; Daneck et al., 2016 ; Ravazzi et al., 2012 ; Slovák et al., 2012 ; Záveská et al., 2021 ); 4 Northeastern Limestone Alps ( Záveská et al., 2021 ); 5 northwestern Balkan Peninsula ( Birks & Willis, 2008 ; Brus, 2010 ; Kirschner et al., 2023 ; Magri et al., 2006 ; Monegato et al., 2015 ; Willis & van Andel, 2004 ; Willner et al., 2009 ); 6 Western Carpathians ( Birks & Willis, 2008 ; Daneck et al., 2016 ; Magri et al., 2006 ; Magyari, 2002 ; Urbaniak et al., 2018 ; Willis & van Andel, 2004 ); 7 Southern Carpathians & Apuseni Mountains ( Kirschner et al., 2023 ; Magri et al., 2006 ); 8 central Balkan Peninsula ( Daneck et al., 2016 ; Kirschner et al., 2023 ; Magri et al., 2006 ; Slovák et al., 2012 ); Most dominant European tree species are wind-pollinated, which facilitates extensive pollen-mediated gene flow over long distances, and promotes homogenisation of gene pools ( Buiteveld et al., 2007 ; Comps et al., 2001 ; Wessinger, 2021 ). While this ultimately confers blurred genetic signals, especially on smaller spatial scales, the application of a large number of genomic markers has recently refined our knowledge of genetic structure in some wind-pollinated tree species on a continental scale ( Milesi et al., 2024 ). In contrast to trees, most temperate forest understorey species (hereafter FUS) are insect-pollinated herbs with limited pollen and seed dispersal capacity, which implies reduced rates of gene flow and more conserved patterns of genetic variation as compared to trees, especially on smaller spatial scales. Phylogeographic studies on FUS thus refined our understanding of the Pleistocene history of temperate forests in Europe in unparalleled detail ( Kirschner et al., 2023 ; Slovák et al., 2012 ; Urbaniak et al., 2018 ; Záveská et al., 2021 ). Specifically, alongside corroborating refugia inferred for trees, they pointed at the existence of additional forest refugia along the southern margin of the Southern Limestone Alps ( Fig. 1 ; Cyclamen purpurascens , Slovák et al., 2012 ; Helleborus niger , Záveská et al., 2021 ), and in the Northeastern Limestone Alps ( H. niger; Záveská et al., 2021 ). Beyond fossil records and phylogeographic evidence, ecological niche models (ENMs) have emerged as the third major source of evidence to infer areas with climatic stability throughout the glacial-interglacial climate fluctuations of the Quaternary ( Gavin et al., 2014 ). When projected to past climatic conditions, such models allow the integration of ecological data and provide detailed information on environmentally suitable areas under different climatic conditions ( Kirchheimer et al., 2018 ). Specifically, ENMs supported the presence of European beech in northern refugia, such as the Carpathians ( Saltré et al., 2013 ), and also their absence from northern Europe ( Svenning, Normand, & Skov, 2008 ). In contrast to trees, ENMs were rarely used to reconstruct refugia of FUS. Willner et al . 2023 modelled dispersal processes to explain the highly idiosyncratic distribution patterns of a set of species including the ones used in this study. They found that the small (relative to those of trees) and fragmented ranges of many FUS are best explained by a combination of broad ecological niches and rare medium- and long-distance dispersal events. Stochasticity is thus an important determinant of current species ranges, explaining the idiosyncratic distribution patterns of the study species despite strong similarities in refugia, ecological tolerances and dispersal abilities. Approaches combining genomic data and ENMs have been extensively used to provide large-scale inferences about past range shifts and refugial dynamics of species and entire biomes ( Carnicero et al., 2022 ; Kirschner et al., 2020 , 2023 ; Pan et al., 2020 ; Rota et al., 2024 ). Genetic and ecological cues commonly interpreted as evidence of refugial survival of populations or lineages are presented in Fig. 2 . In essence, alongside pre-LGM divergence followed by strong isolation, such cues rely on the spatial aggregation of private genetic variation and the geographic distribution of distinct genetic lineages ( Widmer & Lexer, 2001 ; Tribsch & Schönswetter, 2003 ; Paun et al ., 2008 ; Šrámková-Fuxová et al ., 2017 ). However, a direct translation of such patterns into hypotheses about the number and location of cold-stage refugia would be overly simplistic, as demographic processes such as bottlenecks could ultimately produce very similar patterns of genetic variation in a much shorter time. As a possible solution to this dilemma, demographic modelling based on genomic data allows discriminating between evolutionary scenarios, such as long-term refugial survival and post-LGM foundation of a population, in a rigorous statistical framework ( Gutenkunst et al., 2009 ), providing at the same time deeper insights into these processes ( Kirschner et al., 2023 ; Rota et al., 2024 ; Záveská et al., 2021 ). Download figure Open in new tab Fig. 2. Criteria applied for distinguishing between evolutionary hypotheses involving multiple refugia versus a single refugium for the studied forest understorey species. Black dots refer to the set of sampled populations. Light blue areas indicate the maximum extent of the LGM glaciation (∼21 ka; Ehlers et al ., 2011 ). A . Distribution of populations with elevated private allelic richness (reflected by dot size) as an indicator for long-term stability of a population and as indication for a refugial area ( Widmer & Lexer, 2001 ). B . Number and geographic distribution of genetic clusters (multiple geographically structured clusters vs. single uniform cluster). C . Branching of phylogenetic networks (distinct branching reflecting geography vs. shallow or no branching). D . Age of divergence between distinct groups or lineages (pre-LGM splits vs. post LGM splits). E . Number and location of areas providing suitable climatic niches under LGM climatic conditions (multiple, and separated suitable areas vs. a single suitable area). Here, we combine information from genomic data and ENMs, all obtained from a range-wide population sampling of three FUS ( Aposeris foetida, Cardamine trifolia and Hacquetia epipactis ) with a strong preference for forests dominated by European beech. Their ranges, although overlapping with previously proposed forest refugia, are discontinuous and much smaller than those of beech, suggesting limited dispersal capabilities ( Willner et al., 2009 ). After establishing that the study species constitute a suitable proxy for temperate forests using climate and vegetation data, we asked the following questions. (1) Do the disjunct current distributions of FUS result from long-term survival in multiple glacial refugia followed by limited dispersal to adjacent areas or rather from post-LGM colonisation from a single refugium? (2) Do the identified refugia overlap across FUS in spite of their divergent ranges, and do they match areas that have previously been inferred as temperate forest refugia? Materials and Methods Study species Three entomophilic species with strong preference for deciduous forests over carbonate bedrock usually dominated by European beech were studied. Aposeris foetida (Asteraceae; distribution area: Fig. 3A ) is a hemicryptophyte rosette plant; clonal spread is negligible and there is no adaptation to wind dispersal ( Abs, 1994 ; Chytry et al., 2021 ). The species is distributed from the colline to the subalpine belt (40–2,080 m a.s.l; Chytrý et al ., 2016 ). Cardamine trifolia (Brassicaceae; distribution area: Fig. 3A ) is a geophyte with pronounced vegetative reproduction via rhizomatous growth ( Chytry et al., 2021 ) distributed from the colline to the high montane belt (55–1,860 m a.s.l; Chytrý et al ., 2016 ). Hacquetia epipactis (Apiaceae; distribution area: Fig. 3A ) is a geophyte with a short rhizome and thus limited vegetative reproduction ( Chytry et al., 2021 ). It thrives from the colline to the montane belt (50–1,600 m a.s.l.; Chytrý et al ., 2016 ). Download figure Open in new tab Fig. 3. Structure of genetic variation in Aposeris foetida , Cardamine trifolia and Hacquetia epipactis based on RADseq data. A . Geographic distribution of private allelic richness. The size of circles reflects the number of private alleles in each population (absolute number of private alleles in Table A1; map including population ID in Supplementary Fig. 1). Green polygons indicate the current distribution of each species ( Willner et al., 2023 ). B . Geographic distribution of genetic clusters. Pie charts depict the proportional representation of genetic clusters within each population, as identified by the optimal Bayesian clustering solution ( Aposeris foetida , K=6; Cardamine trifolia , K=6; Hacquetia epipactis , K=2). The extent of glaciation at the Last Glacial Maximum (∼21 ka) is indicated in light blue ( Ehlers et al., 2011 ). Colored outlines and labels correspond to geo-genetic groups used in genetic demographic modeling; groups not used are outlined and labelled in light gray. C . Phylogenetic networks; colored lines indicate geo-genetic groups as in B. Pie charts show the proportion of genetic clusters in the populations belonging to the respective branches, as identified by genetic clustering. Collection of plant material Samples of A. foetida , C. trifolia and H. epipactis were collected between 2016 and 2019 from 123, 92 and 50 sites, respectively, aiming at covering their entire distribution ranges ( Fig. 3A , Supplementary Table 1). Sampling localities (here termed ‘populations’) were numbered from west to east, preceded by species acronyms (Af, Ct and He; Supplementary Fig. 1). From each population, leaves were collected into silica gel from three individuals spaced at least 5 m to avoid sampling of clones; a herbarium voucher was made and deposited in the herbarium of the University of Innsbruck (IB). Sampling was carried out in accordance with the Nagoya Protocol and collection permits were obtained if required (details can be obtained from P. Schönswetter on request). DNA extraction Total genomic DNA was extracted from silica-gel dried material with two different methods. For A. foetida , we used a modified CTAB-protocol ( Tel-Zur et al., 1999 ). Prior to extraction with a CTAB buffer the ground tissue was washed three times with a buffer containing sorbitol to remove polysaccharides. DNA extracts were purified using the NucleoSpin gDNA clean-up kit (Macherey-Nagel, Düren, Germany). DNA of C. trifolia and H. epipactis was extracted using the innuPREP Plant DNA Kit II (Analytik Jena, Jena, Germany) as described by ( Kirschner et al., 2022 ). The DNA concentration of all extracts was estimated using a Qubit 4 fluorometer (ThermoFisher Scientific, Waltham, MA, USA). RADseq protocol, identification of RADseq loci and SNP calling Single-digest RADseq libraries were prepared using the restriction enzyme PstI (New England Biolabs; Ipswich, Masschusetts, USA) and a protocol adapted from ( Paun et al., 2016 ). Briefly, we started with 110 ng DNA per individual and ligated 100 mM P1 adapters to the restricted samples. Shearing by sonication was performed with a M220 Focused-ultrasonicator (Covaris; Woburn, Massachusetts, USA) with settings targeting a size range of 200–800 bp and a mode at 400 bp (peak in power: 50, duty factor 10%, 200 cycles per burst and treatment time 90 s at 20°C). Libraries were sequenced on Illumina HiSeq (Illumina; San Diego, California, USA) at VBCF NGS Unit as 100 bp single-end reads ( http://www.vbcf.ac.at/ngs/ ). The raw reads were quality filtered and demultiplexed based on individual-specific barcodes using Picard BamIndexDecoder included in the Picard Illumina2bam package (available from https://github.com/wtsinpg/illumina2bam ) and the program process_radtags.pl implemented in Stacks 2.3( Catchen et al., 2013 ; Rochette et al., 2019 ). The RADseq loci were further assembled, and SNPs were called using the ‘denovo_ map.pl’ pipeline also implemented in Stacks. Denovo_map.pl was run on subsets of the raw data to infer the parameters for an optimal loci yield following( Paris et al., 2017 ). Consequently, a dataset used for subsequent phylogenetic reconstruction was built using a minimum coverage to identify a stack of 5× (-m 5); a maximum number of differences between two stacks in a locus in each sample of three / six (-M 3 / -M 6) and a maximum number of differences among loci to be considered as orthologous across multiple samples of three / six (-n 3 / -n 6) for A. foetida and C. trifolia /H. epipactis . The program populations implemented in Stacks 2.3 ( Catchen et al., 2013 ; Rochette et al., 2019 ) was used to export a set of RADseq loci for the different analyses with the following filtering options: minimum percentage of individuals across populations required to process a locus of 80% / 70% (-R 0.8 / -R 0.7) for A. foetida / C. trifolia and H. epipactis ; and maximum observed heterozygosity required to process a nucleotide site at a locus of 65% (-- max-obs-het 0.65). In addition, the --write-single-snp flag was used to select only a single (first) SNP per fragment to avoid linked loci for STRUCTURE ( Pritchard et al., 2000 ) and δaδi ( Gutenkunst et al., 2009 ) analyses. These SNPs were exported in structure format (-- structure) for STRUCTURE ( Pritchard et al., 2000 ) analysis and vcf format (--vcf) for all other analyses. The datasets were converted in nexus format for SplitsTree with vcf2phylip( Ortiz, 2019 ). Finally, for δaδi( Gutenkunst et al., 2009 ) analyses, SNPs (without filtering for minimum percentage of individuals across populations) were filtered using VCFtools 0.1.16( Danecek et al., 2011 ) by keeping only sites with depth values over all individuals (DP) between 10 and 50 (--minDP 10; --maxDP 50) and genotypes with a quality above 20 (--minGQ 20). Exploratory analyses of SNP data To estimate the number of private alleles and the nucleotide diversity (π) per population, we used the program populations in Stacks 2.3 ( Catchen et al., 2013 ; Rochette et al., 2019 ). In addition, Neighbour-Nets were produced with SplitsTree 4.17.1 ( Huson & Bryant, 2006 ). To determine the number of genetic clusters and the extent of admixture, we used STRUCTURE 2.3.4 ( Pritchard et al., 2000 ) with ten replicate runs for each K, the number of groups, ranging from 1 to 10, carried out using a burn-in of 20,000 MCMC iterations followed by 200,000 additional iterations. The dataset included 340 individuals and 13,738 loci for A. foetida , 232 individuals and 8,438 loci for C. trifolia and 127 individuals and 10,808 loci for H. epipactis . We summarised the results with CLUMPAK ( Kopelman et al., 2015 ). The optimal K was identified as the K, where the increase in likelihood started to flatten out, the results of replicate runs were similar, and the clusters were non-empty. Additionally, the deltaK criterion was employed, reflecting an abrupt change in likelihood of runs at different K ( Evanno et al., 2005 ). Demographic modelling based on genomic data First, intraspecific genetic groups were defined within each of the three species based on a combination of STRUCTURE clustering, Neighbour-Net networks and geographic distributions (further on termed as ‘geo-genetic groups’). In doing so, only populations for which major cluster assignment in STRUCTURE was > 50%, and assignment to second most frequent cluser was < 35% were included, as the underlying demographic signal in admixed populations may be hidden by gene flow ( Chikhi et al., 2010 ). These criteria were e.g. not fulfilled for the central group of A. foetida (see Results). Intraspecific geo-genetic groups were further defined as a set of populations for which at least two of the three following criteria applied: i) the populations belonged to a distinct genetic cluster (from STRUCTURE), ii) the populations formed a divergent branch in the phylogenetic network, and/or iii) the populations were geographically very isolated. For each species, alternative divergence scenarios were tested between pairs of geo-genetic groups (see Supplementary Table 2 for a list of all pairwise comparisons). This was done using the dadi pipeline 3.1.1 ( Portik et al., 2017 ) that is based on the diffusion approximation method of δaδi that analyses 2D-JSFS( Gutenkunst et al., 2009 ). Model selection included nine predefined models capturing different divergence scenarios ( Charles et al., 2018 ; Portik et al., 2017 ; Záveská et al., 2021 ); graphical representation in Supplementary Fig. 2). The selected models capture changes in population size, migration, and divergence times between groups over one or two epochs. Most importantly, the selected models also allow identification of the mode of divergence, e.g. if it was due to a vicariance event (three models) or a founder event (defined as post-divergence exponential growth of the derived population ( Charles et al., 2018 ). To evaluate directionality of divergence events, the latter models were tested in both directions, meaning that each model was tested twice placing each geo-genetic group as source and founding populations. The 2D-JSFS were prepared with easySFS ( https://github.com/isaacovercast/easySFS ) and with the Find-Best-Projections template of ( Portik et al., 2017 ). Down-projecting allows maximising the number of segregating sites by reducing the number of individuals (alleles) in the 2D-JSFS. The numbers of alleles for each down-projection for each group are shown in Supplementary Table 2. Analyses in the dadi pipeline ( Portik et al., 2017 ) were done using the following settings. All models were initially optimised with 3-fold perturbed random starting parameters followed by two rounds of 2-fold and one round of 1-fold perturbation to estimate the log-likelihood of the SFS given the model (folds = 3,2,2,1). For each round, multiple replicates were run (reps = 60,70,70,80) and parameter estimates were used from the best scoring replicate (highest log-likelihood) to seed searches in the following round (maxiters = 10,10,10,15). Across all analyses, optimised parameter sets of each replicate and a multinomial approach were used to estimate the log-likelihood of the 2D-JSFS given the model. Retaining only a single SNP per RADseq locus allows assuming that loci are unlinked, thus the returned log-likelihood values represent the true likelihoods according to ( Portik et al., 2017 ). Model selection was done using the Akaike information criterion (AIC) and Akaike weights (ωi; Burnham & Anderson, 2002 ). In addition, goodness-of-fit tests were done following ( Barratt et al., 2018 ). This involved running 100 simulations to obtain a distribution of log-likelihood scores. Pearson’s chi-squared test statistic was used to subsequently check whether the empirical values from the original model selection were inside these distributions. Assuming that SNPs were unlinked, these simulations represent parametric bootstraps ( Barratt et al., 2018 ). Simulations were consequently used to obtain confidence intervals for the parameters ‘fraction of the ancestral population’ (‘s’) and ‘divergence time’ (‘T’, ‘T1’ and ‘T2’), which were then transformed into biological estimates following ( Gutenkunst et al., 2009 ). This was done using the mutation rate of Arabidopsis thaliana (7e-9 base substitutions per site per generation; Ossowski et al ., 2010 ) and generation times (10 years for weakly clonal A. foetida and H. epipactis ; 20 years for strongly clonal C. trifolia ). To ascertain whether the simulated times were significantly older than the LGM, a one-sided Wilcoxon test was conducted using the R package ‘stats’ ( R Core Team, 2023 ) to compare the distribution of the simulations with a null hypothesis (i.e. lower than LGM) set to 20 ka, which allows obtaining a confidence interval and a pseudomedian value. Plastid DNA sequencing Five plastid DNA regions totalling 6,500, that is atpB-rbcL ( Piwczyński et al., 2018 ), ndhF – trnL , ndhJ – trnT , trnH–psbA and trnQ – trnK ( Shaw et al., 2005 , 2007 ), were inspected for variability in A. foetida . For C. trifolia , the five before-mentioned plastid regions plus rpo B- trn C and trn G -trn S ( Shaw et al., 2005 , 2007 ), totalling 8,200 bp, were screened. The trnH – psbA spacer was most variable and was thus amplified and sequenced for one individual per population from 111 and 83 populations in A. foetida and C. trifolia (Supplementary Table 1). Amplification was performed in a reaction mix (total volume 20 µl) containing 9 µl REDTaq ReadyMix PCR Reaction Mix with MgCl 2 (Sigma-Aldrich, St. Louis, MO), 1 µl of the primers (10 µM) trnH and psbA ( Shaw et al., 2005 ) and 1 µl DNA template. Cycling conditions were 5 min at 95°C, 30 cycles of 30 sec at 95°C, 60 sec at 48°C and 60 sec at 72°C, followed by 10 min at 72°C. For H. epipactis , four plastid DNA regions totalling 7,500 bp, that is ndhJ – trnT , ndhF – trnL , trnG – trnS and trnQ – trnK ( Shaw et al., 2005 , 2007 ) were inspected for variability. The trnQ – trnK region was the only variable and its terminal parts (parts of trnQ – rps16 spacer and rps16 – trnK spacer) were thus amplified and sequenced for one individual per population from 40 and 39 populations, respectively (Supplementary Table 1). Despite several attempts we failed to amplify both the rps16–trnK and trnQ–rps1 6 spacers for individuals of populations He_2, He_19, He_20, He_34 and He_39. Amplification was performed in a reaction mix (total volume 20 µl) containing 7.7 µl REDTaq ReadyMix PCR Reaction Mix with MgCl 2 (Sigma-Aldrich, St. Louis, MO), 1 µl of both primers (10 µM) trnQ and trnK ( Shaw et al., 2007 ), 1 µl BSA (1 mg/ml) and 1 µl DNA template. Cycling conditions were 5 min at 95°C, 35 cycles of 30 sec at 95°C, 30 sec at 60°C and 4 min at 65°C, followed by 10 min at 65°C. PCR programs were run on Eppendorf 5331 thermocyclers (PE Applied Biosystems, Foster City, CA). PCR products were purified with E. coli Exonuclease I and FastAP (Thermosensitive Alkaline Phosphatase; Fermentas, St. Leon-Rot, Germany) following the manufacturer’s instructions. Sequencing was carried out at Eurofins Genomics (Ebersberg, Germany) using the primer psbA for A. foetida and C. trifolia and the primers trnK and trnQ for H. epipactis . Contigs were assembled, edited and sequences aligned using Geneious Pro 5.5.9 ( Kearse et al., 2012 ). We constructed statistical parsimony networks using PopART ( Leigh & Bryant, 2015 ) applying the TSC network option ( Clement et al., 2000 ). Gaps were treated as fifth character state and indels longer than 1 bp were reduced to single base pair columns allowing those structural mutations to be counted as single base pair mutations only. Modelling environmental suitability Current climatic conditions and those at the LGM were retrieved from the Chelsa Climate database ( Karger et al., 2017 , 2023 ). We selected the LGM data based on the Community Climate System Model version 4 (CCSM4; Gent et al ., 2011 ) as this model outperforms other global circulation models in Europe ( Fordham et al., 2017 ). We used three bioclimatic variables highly relevant for the occurrence of the modelled species: mean annual temperature (bio1), annual precipitation sum (bio12) and precipitation seasonality (bio15). These climatic variables, which are available with a spatial resolution of 30″, were projected to a grid of 1 × 1 km cell size using the nearest neighbour method and checked to avoid high correlation (|Pearson’s r| < 0.6). As additional environmental variable we used topographical roughness as the standard deviation of the elevations of 100 × 100 m cells (derived from https://www.eea.europa.eu/data-and-maps/data/copernicus-land-monitoring-service-eu-dem ) within each 1 × 1 km cell. We used ENMs relating species occurrences (presences/absences) to the four environmental variables described above to identify suitable areas for each of the three study species under current and LGM climatic conditions. Locations of sampled populations were used as presences, while 100 pseudo-absences were drawn using the ‘sre’ method within the BIOMOD2 package ( Thuiller et al., 2024 ). Thereby, pseudo-absences were selected outside a species’ range envelope constructed by the default setting of BIOMOD2 (i.e., trimming the occurrences along each environmental axis to 95%). The selection of pseudo-absences was repeated ten times. ENMs were parameterized within the BIOMOD2 framework by means of four modelling techniques using their default settings: Generalised Linear Models (GLM), Generalised Additive Models (GAM), Boosted Regression Trees (GBM) and Random Forests (RF). Model runs were replicated three times for each set of pseudo-absences using 80% of the occurrence data for model parameterization and the remaining 20% for model evaluation using the True Skill Statistic score (TSS; Allouche et al ., 2006 ). Based on the resulting 120 parameterized models per species (3 replicates × 10 pseudo-absence datasets × 4 modelling techniques), we calculated the environmental suitability as ensemble projections of potential species ranges under current and LGM climatic conditions as weighted (by the TSS value of single models) mean of the projected occurrence probabilities of the single models. Niche overlap Niche overlap between the European beech and each study species was computed following ( Broennimann et al., 2012 ). This comprises the calculation of two-dimensional density distributions (dds) of occurrences for each species and the determination of the niche overlap in this environmental space. First, dds were derived from a PCA based on values of the four environmental variables (standardised to zero mean and unit variance) used for the ENM (including cells all over Europe). Maximum and minimum scores of the first and second axis of the PCA were used to define boundaries of the environmental space, which was divided into 100 bins along each axis resulting in a grid of 10,000 cells each, representing a unique combination of environmental variables. The dd of occurrence points of each lineage was calculated using the function ecospat.grid.clim.dyn applying the adehabit kernel smoother included in the R package ecospat ( Di Cola et al., 2017 ). As a metric for the niche overlap between European beech and each study species, Schöneŕs D ( Schoener, 1970 ) was calculated using the function ecospat.niche.overlap from the same package ( Di Cola et al., 2017 ). Vegetation analysis Vegetation data from the European vegetation archive (EVA; Chytrý et al ., 2016 ) was used to estimate the amount of variability in the species composition of communities where at least one of the three studied species or European beech were recorded. In total, the data set comprised 19,084 vegetation surveys (i.e., plot records) of 100–200 m 2 , of which 1,417 contained A. foetida , 498 C. trifolia and 221 H. epipactis . A total of 18,051 plots contained European beech. All plot records used were limited to the study area (EVA plot identifiers are provided in Supplementary Data 1). An ordination analysis (detrended correspondence analysis, DCA) was done using the R-package ‘vegan’ ( Oksanen et al., 2022 ). Results Illumina reads and population structure Restriction site-associated DNA sequencing (RADseq) datasets resulted in 296,281, 141,559 and 219,077 loci for A. foetida , C. trifolia and H. epipactis after quality filtering, with a mean coverage of 16.0× (±3.2×), 12.3× (±1.0×) and 13.6× (±2.9×), respectively. All resulting raw RADseq data are available in the NCBI Short Archive as BioProject PRJNA824654 (accession numbers SAMN2773504–SAMN27735362, SAMN27409840–SAMN27410080, SAMN28084126–SAMN28084261 for A. foetida , C. trifolia and H. epipactis , respectively; details are given in Supplementary Table 1). The number of private alleles per population ranged from 37 to 991 in A. foetida , from 12 to 296 in C. trifolia and from 14 to 518 in H. epipactis (Supplementary Table 1; Fig. 3A ). Nucleotide diversity (π) amounted to 0.014–0.043 in A. foetida , 0.054–0.167 in C. trifolia and 0.033–0.104 in H. epipactis (Supplementary Table 1; Supplementary Fig. 3). Results from Bayesian clustering analyses for each species are summarised for the respective optimal solution in Fig. 3 , and in detail in Supplementary Fig. 4. For A. foetida , the genetic clustering analysis resulted in an optimal separation into six clusters (best K = 6) based on 13,738 unlinked single nucleotide polymorphisms (SNPs; Fig. 3B ). Six geo-genetic groups were defined following the elaborated rationale, and after excluding few strongly admixed populations ( Fig. 3B,C ; Supplementary Fig. 4). Despite meeting two criteria, that is geographic coherence and genetic separation of populations, the Center group was not divided into distinct groups due to the extensive admixture cline. This was done to circumvent complications in downstream analyses. For C. trifolia , the genetic clustering analysis resulted in an optimal separation into six clusters (best K = 6; Fig. 3B ) based on 8,438 unlinked SNPs. After excluding strongly admixed populations, eight geo-genetic groups were defined using the elaborated criteria ( Fig. 3B,C , Supplementary Fig. 4). For H. epipactis , the genetic clustering resulted in an optimal separation into two genetic groups (best K=2; Fig. 3B ) based on 10,808 unlinked SNPs. Following the exclusion of admixed populations, three geo-genetic groups were defined ( Fig. 3B,C , Supplementary Fig. 4). Demographic history Quality assessment and parameter estimates of demographic models are provided in Supplementary Table 2. Parameter estimates from goodness-of-fit tests for the best models are given in Supplementary Table 3. In the following, all parameter values correspond to median values inferred by goodness-of-fit tests. In all cases, the increase of the log likelihood values and decrease of variance across replicates, for each model, indicated convergence into the best-ranked model. For A. foetida , founder event models with migration fitted the data best ( Fig. 4 ). When modelling the split between the pairs of geo-genetic groups Af_nAlps/Af_Center, Af_nApen/Af_Center and Af_wAlps/Af_Center, the size of the founding fraction ‘s’ was comparably large (s > 0.2; Supplementary Table 3). Initial analyses did not allow the identification of any ancestral relationships between the two subclusters of Af_Center. To reduce complexity in the modelling, this entity was not divided into further groups. For the colonisation of Af_Carp from Af_cBalk, we found strong support for a founder event involving a small founder population (s = 0.05). Finally, Af_Center was colonised from Af_cBalk via a founder event with considerable population size (s = 0.14). Download figure Open in new tab Fig. 4. Results of demographic modeling based on genomic data for Aposeris foetida , Hacquetia epipactis and Cardamine trifolia . The panel to the left shows the three evolutionary scenarios that have been evaluated; each scenario is also represented by a symbol (square, circle, triangle) which corresponds with the symbol indicating the best scenario for each species in the right panel. Maps show the geo-genetic groups that were tested against each other; numbers indicate the pairs included in the respective modeling run. The white trapezoids between geo-genetic groups indicate the directionality of the separation whereas the wider side of the trapezoid part faces towards the respective source group. The aspect ratio of sides facing each group proportionally indicates the fraction of the ancestral group involved in establishing the derived group in each modeling run (parameter s; absolute values in Supplementary Table 3). Barplots show the distribution of split age estimates (n = 100) between two groups; numbers correspond to the pairs on the map and the red lines indicate the Last Glacial Maximum (LGM; 21 ka). Results of C . trifolia are presented in two maps to ensure readability. All eight pairwise comparisons in C. trifolia identified divergence via vicariance followed by secondary contact as the most likely model ( Fig. 4 ). The second-best models received low support in all cases except for the divergence between Ct_nAlps and Ct_Carp (ωi = 0.19). This model indicated a founder event from Ct_Carp to Ct_nAlps involving a small founder population (s = 0.03), followed by secondary contact. For H. epipactis the divergence between He_Carp and He_nwBalk was modelled as vicariance followed by secondary contact. In both best models (ωi = 0.81 and ωi = 0.17, respectively), the split between He_seAlps and He_nBalk was modelled as a pre-LGM vicariance event followed by secondary contact. For each best-performing model, visualisations of the two-dimensional joint site frequency spectra (2D-JSFS) for empirical data and model data, and corresponding residuals are shown in Supplementary Fig. 5. The results of the goodness-of-fit tests indicated a good fit of all models, and resulted in median divergence times predating the LGM (Supplementary Fig. 5, Supplementary Table 3). One-sided Wilcoxon tests confirmed pre-LGM divergence for all studied pairs of geo-genetic groups (Supplementary Table 3). In some cases, however, the simulated divergence times yielded values post-dating the LGM, e.g. in the pair Af_cBalk/Af_Carp, where 12% of the simulated values indicated post-LGM divergence ( Fig. 4 , Supplementary Table 2, Supplementary Table 3). In this particular case, the dating result should be interpreted with caution. Plastid DNA sequencing The trnH – psbA sequences were between 433 (Af_112) and 473 (Af_43) bp long for A. foetida , and 313 bp long for C. trifolia . For H. epipactis , the trnQ–rps16 and rps16–trnK spacer sequences were 843 and 694 bp long, respectively. GenBank accession numbers are given in Supplementary Table 1). The haplotype networks for A. foetida , C. trifolia and H. epipactis identified a central haplotype in each species, that was widespread throughout the distribution range in each species (Supplementary Fig. 6) and present in 75%/83%/81% of the investigated populations, respectively. Divergent haplotypes differed only by one or two mutations, and were distributed across the species’ range as well, with the more common one, aggregated along the Southern and south-eastern Alps in all three species (see Supplementary Fig. 6 for more detail). Environmental niche models, niche overlap and vegetation analysis Regions predicted to be highly suitable by the ENMs were largely consistent across the study species for LGM conditions ( Fig. 5 ) as well as for current day conditions (Supplementary Fig. 7). They comprised large parts of the South and the Southeast of the study area during the LGM, mostly outside of the species current distribution ranges. In contrast, present-day predictions covered well the actual distributions of the species (including all identified refugia) and, in addition, largely matched the current distribution of European beech. All ENMs showed high evaluation scores (TSS = 0.91, 0.93 and 0.93 for A. foetida, C. trifolia and H. epipactis , respectively). Download figure Open in new tab Fig. 5. Predicted environmental suitability (indicated by the color shade ramp) for Aposeris foetida (A), Cardamine trifolia (B) and Hacquetia epipactis (C) at the Last Glacial Maximum (∼21 ka; LGM) obtained from environmental niche models (ENM). Coastlines reflect the lower sea level at the LGM (−120 m; Becker et al ., 2015 ); LGM glaciation ( Ehlers et al., 2011 ) is indicated in light blue. The niches of A. foetida and C. trifolia covered that of European beech to similar proportions, while H. epipactis showed a slightly narrower niche ( Fig. 6 ). Each study species rarely occurred outside the niche of European beech. Aposeris foetida and H. epipactis extended towards higher and lower annual precipitation sums represented by bio12. Cardamine trifolia extended towards higher precipitation seasonality and topographic roughness (bio15 and topo). Hacquetia epipactis showed two ecological centres of distribution. The empirical niche overlap measured as Schöneŕs D between European beech on the one hand and A. foetida , C. trifolia and H. epipactis on the other hand is D=0.37, D=0.32 and D=0.40, respectively. Download figure Open in new tab Fig. 6. Niche of each forest understorey species (FUS; (shaded area) compared to that of European Beech ( Fagus sylvatica ; contour line) in environmental space (i.e. based on the first and second axis of a PCA using four environmental variables). Shading shows the density of the occurrences of the FUS. A proportion of 37.3% and 30.4% of inertia are explained by the first and second axis of the PCA, respectively. The insert illustrates the contribution of the environmental variables: bio1, annual mean temperature; bio12, annual precipitation; bio15, precipitation seasonality; topo, standard deviation of elevations. The first and second DCA axis of the vegetation plot records explained 33.2% and 25.9% of variance, respectively. The plot records containing any of the three study species were placed within the variability of records with European beech. Aposeris foetida covered most of the ordination space occupied by records with European beech, whereas C. trifolia and H. epipactis filled a considerably smaller proportion (Supplementary Fig. 8). Discussion The classic hypothesis of Pleistocene cold-stage refugia of European temperate tree species being located in the South of the European peninsulas ( Hewitt, 1999 ; Taberlet et al., 1998 ) has been challenged by work demonstrating that trees also survived further north in “northern refugia” ( Fig. 1 ). We show that this is also the case in three FUS and outline evolutionary scenarios that involve multiple northern forest refugia in Central Europe. The inferred FUS refugia showed a large spatial overlap between the three species, coinciding well with some of the previously proposed temperate forest refugia ( Fig. 1 ), and also suggesting the presence of forests in areas where this was previously unclear such as the northeasternmost Alps. Among all the identified northern refugia, the northwestern Balkan Peninsula played the most eminent role, further highlighting the status of this region as a diversification centre of European forest biota in general. We show that the three FUS Aposeris foetida , Cardamine trifolia and Hacquetia epipactis are well-suited proxies for reconstructing the late Quaternary history of Central European temperate forests because their climatic niches are almost entirely restricted to the climatic niche of European beech ( Fig. 6 ) and they occur largely in the same plant communities as European beech (Supplementary Fig. 8). However, some FUS might have been resilient against past changes in the tree composition of temperate forests ( Záveská et al., 2021 ). As a consequence, we suggest that the refugia inferred in this study might indicate the existence of refugia of broadly defined – deciduous, coniferous or mixed – temperate forests. The northwestern Balkan Peninsula acted as primary diversification centre The northwestern Balkan Peninsula and adjacent areas were the primary centre of diversification and acted as the main refugial area for the investigated FUS. As expected for refugia, populations with the highest numbers of private alleles, a proxy for long-term stability ( Fig. 3A ) and elevated nucleotide diversity (Supplementary Fig. 3) were found in this area. Most of the genetic clusters occurring there exhibit ranges extending beyond the northwestern Balkan Peninsula, which either underlines the region’s role as source for the post-LGM colonisation of other areas or, alternatively, reflects vicariance resulting from past wider continuous distributions ( Španiel & Rešetnik, 2022 ; see following sections). Accordingly, we conclude that this area has acted as a major refugium for temperate forests, which is in line with evidence from recent studies on other FUS ( Kirschner et al., 2023 ; Slovák et al., 2012 ; Záveská et al., 2021 ) and on European beech ( Fig. 1 ; Magri et al ., 2006 ; Willner et al ., 2009 ). The northwestern Balkan Peninsula was, however, not a single, homogeneous refugium harbouring a panmictic group of refugial populations; rather, it was subdivided and comprised distinct, weakly isolated groups of populations ( Fig. 3 ; Slovák et al ., 2012 ). We therefore suggest a scenario where post-LGM expansion from these isolated refugia and subsequent gene flow have caused the weak, but still extant substructure and admixture within the area ( Fig. 3B,C ). The persistence of genetic structure in the presence of gene flow during post-LGM recolonisation contrasts the assumption that gene flow erases any genetic signal for refugia ( Hewitt, 2000 ; Keppel et al., 2012 ), but is a prevalent pattern given sufficient resolution of markers ( Carnicero et al., 2022 ; Kolář et al., 2015 ). Survival in multiple northern refugia beyond the northwestern Balkan Peninsula Beyond the northwestern Balkan Peninsula, northern forest refugia were previously suggested in the Western and Southern Carpathians, the western, northeastern and southern margin of the Alps, the central Balkan Peninsula and/or the northern Apennines (referenced summary in Fig. 1 ; details in the Supplementary Information). We find refugia linked to almost all disjunct parts of the distribution areas of the studied FUS ( Figs. 3 , 7 ). This suggests relative distributional stasis (i.e. no large-scale expansion beyond the refugial range) since the LGM, a phenomenon often observed in species with limited dispersal abilities ( Carnicero et al., 2022 ; Dullinger et al., 2012 ; Kropf et al., 2002 ; Schneeweiss & Schönswetter, 2010 ; Svenning, Normand, & Kageyama, 2008 ). Indeed, simulations of the postglacial range formation of our study species showed that in the absence of local refugia their disjunct distribution patterns can only be explained by rare long-distance dispersal events ( Willner et al., 2023 ). In addition, demographic models based on genomic data found substantial migration between the tested groups for all three study species ( Fig. 4 , Supplementary Table 3). However, geo-genetic groups originating from different refugia maintained their genetic integrity over time, even if barriers must have been permeable at different time horizons. This was supported by the best-performing models ( Fig. 4 , Supplementary Tables 2, 3), and corresponds to evolutionary scenarios outlined for European temperate tree species which similarly exhibited strong genetic structuring despite recurring gene flow ( Milesi et al., 2024 ). In case of C. trifolia and H. epipactis , almost all models estimated gene flow after an epoch of isolation (Supplementary Fig. 5), which can be interpreted as secondary contact of populations in the course of post-LGM range expansion, which may also support the absence of haplotype structure (Supplementary Fig. 6). We speculated that this indicates ample distributions of the investigated FUS during a warmer and more humid past such as the Holocene climatic optimum ( Doláková et al., 2010 ), when temperate forest thrived also in lowlands such as the Pannonian Plains (the lowlands between the Alps and the Carpathians). Intraspecific diversification predating the LGM provides strong support for survival in multiple northern refugia. Divergence time estimates derived from the demographic models ( Fig. 4 , Supplementary Fig. 5, Supplementary Table 3), suggested splits significantly before the LGM (median divergence times: A. foetida 40–79 ka, C. trifolia 556–897 ka, H. epipactis 554–828 ka; Supplementary Table 3; full range of estimates in Fig. 4 ), providing direct evidence against post-LGM colonisations. In A. foetida , however, a small fraction of simulations also resulted in post-LGM splits (12% in nAlps/center; 4% in cBalk/Carp; Fig. 4 , Supplementary Table 3, Supplementary Fig. 5), acknowledging that the involved areas (nAlps, cBalk) might not have been refugia. In addition, no explicit demographic test was done between the Southern Limestone Alps and the northwestern Balkan Peninsula in A. foetida due to strong admixture of these groups. It thus remains open whether this species also survived in refugia along the southern margin of the Alps. Divergence events in C. trifolia and H. epipactis were best represented by demographic models capturing vicariance whereas founder event models represented divergence events in A. foetida best. As argued above, this suggests ample, though probably discontinuous distributions of C. trifolia and H. epipactis during interglacials, which were probably similar to the situation in the Holocene. Such phases were followed by divergence driven by isolation in refugia due to unfavourable climatic conditions during cold stages ( Fig. 4 , details provided in Supplementary Information). Instead, founder events in A. foetida may point to an important role of stochastic dispersal in attaining its present distribution (i.e., long-distance dispersal; Willner et al ., 2023 ). However, we emphasise that in the demographic modelling approach used, founder events were defined as an exponential increase of population size of the founded population after divergence ( Charles et al., 2018 ), and consequently, exponential increase of the founding population size may always be resolved as a founder event, irrespective of the initial size of the founding population (s parameter; Supplementary Fig. 2). This contrasts the strict definition of a founder event, which per definition should involve only a small number of individuals ( Mayr, 1954 ; Templeton, 2008 ). Indeed, in A. foetida three out of five models estimated large proportions of the source populations’ effective population size in the founding population (s > 0.2; Supplementary Table 3), hinting at vicariance rather than strict founder events. Partially asynchronous divergence and ubiquitous dispersal limitation Intraspecific splits within C. trifolia and H. epipactis occurred between ca. 550 ka and 900 ka, significantly predating the LGM ( Fig. 4 , Supplementary Table 3). As such, these diversification ages are comparable with what has been observed for groups and lineages of other European FUS ( Kirschner et al., 2023 ; Valtueña et al., 2012 ; Záveská et al., 2021 ), between Carpathian endemic forest plants and insects, and their sister species ( Mráz & Ronikier, 2016 ), and, recently, also for temperate European tree species (> 600 ka; Marchesini et al ., 2023 ; Milesi et al ., 2024 ). Notably, in many of the abovementioned instances (including C . trifolia and H . epipactis ), periods of diversification roughly coincide with the Mid-Pleistocene Transition (0.7–1.25 Ma), when cold stages became significantly longer (ca. 100 ka instead of ca. 40 ka), and more intense ( Herbert, 2023 ), and forests were recurrently replaced by the expanding steppes ( Kirschner et al., 2022 ). In contrast, and as outlined above, divergence within A. foetida occurred more recently ( Fig. 4 , Supplementary Table 3; Supplementary Fig. 6), and was likely shaped by stochastic dispersal events. These younger dispersal events might, however, have erased the signal of older differentiation as has been observed when high numbers of propagules were involved ( Todesco et al., 2016 ). In any event, the lack of synchronicity of divergence events in the study species highlights that codistribution does not necessarily imply congruent evolutionary trajectories ( Dawson, 2014 ). We conclude that the intraspecific diversifications of FUS, although not being synchronous, are much older than the last glacial period, and suggest that instead the Mid-Pleistocene Transition may have significantly shaped the evolution of European temperate forest biota. Only some of the LGM refugia inferred from genomic data were predicted to be suitable under LGM climatic conditions by ENMs ( Fig. 5 ), while under current climatic conditions, large areas outside of the study species’ actual distribution ranges were found to be suitable (Supplementary Fig. 7). This likely suggests that the fundamental niches of the study species are considerably larger than the currently realised niches ( Willner et al., 2023 ); only the latter are captured by the ENMs. Such discrepancies can be caused by strongly restricted dispersal abilities, preventing the species from reaching areas with climatic conditions that would actually be within their fundamental niche. Indeed, strong dispersal limitation of the studied FUS was suggested by simulations ( Willner et al., 2023 ). It is, for example, very unlikely that the studied FUS inhabited the large areas in Southern Europe predicted as suitable under LGM conditions ( Fig. 5 ), and went completely extinct there after the last cold stage. Most of these regions still harbour deciduous forests similar to those, where the study species occur today ( Willner et al., 2009 ), rendering it likely that the species have actually never reached these areas. Moreover, species might have adapted to the changing climate since the LGM ( Lavergne et al., 2010 ; Luqman et al., 2023 ; Wessely et al., 2022 ) which may have contributed to the observed incongruence of highly suitable regions and suggested refugia ( Figs. 5 , 7 ). Temperate forest refugia in Europe – a synthesis The refugia suggested for the individual study species ( Fig. 7 ) overlap and mostly match refugia previously suggested for trees based on palynological records and macrofossils ( Fig. 1 ; for an in depth discussion of refugial dynamics for each FUS and each region see the Supplementary Information). Specifically, the northwestern Balkan Peninsula, an unequivocal source of colonisation of Central Europe by European beech ( Magri et al., 2006 ), emerged as a strongly supported refugium based on genomic data from our three study species and additional FUS ( Fig. 7 ; Záveská et al ., 2021 ; Kirschner et al ., 2023 ). The southwestern fringe of the Western Alps, another well-established refugium ( Magri et al., 2006 ), was supported for A. foetida . Refugia at the southern margin of the Alps ( Fig. 1 ; Magri et al ., 2006 ; Kaltenrieder et al ., 2009 ; Gubler et al ., 2018 ) were indicated for C. trifolia and Helleborus niger , with inconclusive evidence for A. foetida and Cyclamen purpurascens ( Figs. 4 , 7 ; Slovák et al ., 2012 ; Záveská et al ., 2021 ). The same applies to the northern Apennines, which were identified as a refugium for C. trifolia and possibly also for A. foetida ( Figs. 4 , 7 ). Download figure Open in new tab Fig. 7. Refugia inferred for understorey species (FUS) of Central European temperate forests based on evidence from genome-wide SNPs ( Aposeris foetida: ‘ A . f .’ , Cardamine trifolia : ‘C . t .’, Hacquetia epipactis : ‘H . e .’, this study; E . carniolica : ‘ E . c .’, Kirschner et al ., 2023 ; H . niger : ‘ H . n .’, Záveská et al ., 2021 ). Photo credits: Cäcilia Lechner-Pagitz ( A . f ., C . t ., H . e ., H . n .); Božo Frajman ( E . c .). The coloring of the outer tiles and the inner tiles of each circle indicate if a species occurs in the area and if there is genomic evidence for refugial survival of the species, respectively. Genomic data further support a cold-stage FUS refugium in the Northern Limestone Alps for A. foetida and possibly also for C. trifolia ( Figs. 4 , 7 ). This is remarkable, given that the existence of temperate forests during the LGM has only recently been suggested by genomic data in this region, and so far only in its northeasternmost periphery ( Záveská et al., 2021 ). For the eastern part of the study area, our data support previous findings of glacial refugia ( Kirschner et al., 2023 ; Magri et al., 2006 ; Marchesini et al., 2023 ). These refugia have been identified in the Western Carpathians ( H. epipactis and possibly C. trifolia ), the Southern Carpathians and the Apuseni Mountains ( A. foetida , Figs. 1 , 7 ) as well as the central Balkan Peninsula ( A. foetida , C. trifolia ). Taken together, our results fully support that the Northwestern Balkan Peninsula was the preeminent refugium for the investigated FUS; with increasing distance from there, the number of FUS occurring in each potential refugium decreases ( Fig. 7 ). This may suggest that smaller refugia within the area may have limited capacity to preserve a major fraction of temperate forest diversity, similar to the limited carrying capacity of small islands in island biogeography ( MacArthur & Wilson, 1967 ). In any event we emphasise that – disregarding three uncertain cases highlighted in Fig. 7 – all of the disjunct FUS occurrences derive from populations locally present before the LGM, suggesting widespread occurrence of temperate forests through glacial periods. The existence of widespread northern refugia implies shorter migration distances in recolonizing Europe and explains how temperate forests could have colonised Europe rapidly after the Pleistocene. Based on evidence from FUS, we demonstrate the widespread occurrence of such refugia, and thereby contribute to resolving the long-standing enigma of Reid’s paradox (J. S. Clark et al., 1998 ). However, in the same line as argued by Záveská et al . (2021) , our findings hint at the existence of patches of temperate forests in cold-stage conditions rather than providing direct evidence for large-scale deciduous forests in such periods. Some of the outlined refugia may have been dominated by conifers, while still offering suitable conditions for FUS. Indeed, evidence for the existence of different coniferous trees during cold-stages exists for most of the suggested northern refugia or close areas ( Fig. 1 ; Willis & van Andel, 2004 ; van der Knaap et al ., 2005 ; Huber et al ., 2010 ; Monegato et al ., 2015 ; Gugerli et al ., 2023 ). Even under present day conditions A. foetida and C. trifolia are known to occasionally occur in subalpine dwarf mountain pine ( Pinus mugo ) communities. Finally, we conclude that although FUS refugia do not directly indicate the presence of continuous cold stage forests, they do indicate the presence of trees and the existence of closed canopy cover in all the respective refugial areas outlined in this study. Competing interests The authors declare no competing interests. Author Contributions P.C., P.K., and P.S. designed the study. E.Z. generated the genomic data. C.V., P.C. and P.K. analyzed genomic data. J.W., K.H., and W.W. were responsible for ecological modelling. C.V., K.H., P.C., P.K., P.S., and W.W. co-wrote the manuscript. All authors contributed to the development of the manuscript and improved earlier drafts of the paper. Data availability Short read sequence data generated in this study is available from NCBI short read archive BioProjects PRJNA824654 (accession numbers SAMN2773504–SAMN27735362, SAMN27409840–SAMN27410080, and SAMN28084126–SAMN28084261). Data from chloroplast DNA sequencing was uploaded to NCBI genbank and accession numbers are provided in Supplementary Table 1. Source data for vegetation analyses were obtained from European Vegetation Archive (EVA) and plot identifiers are provided in Supplementary Data 1). Acknowledgments The paper is dedicated to Christoph Dobeš, who contributed significantly to the research project, from which this paper emerged. He died in a mountain accident in March, 2020. This work was financed by the Austrian Science Fund (FWF, project P29413 ‘Range formation of beech forest understory herbs’ to PS). Marianne Magauer and Daniela Pirkebner performed most of the lab work. Faruk Bogunić, Christoph Dobeš, Julian Haider, Judita Kochjarová, Bohumil Mandák, Ivana Rešetnik, Sabahudin Solaković, Simon Stifter and Ernesto Venturi sampled most of the studied populations. Christoph Dobeš and Iwona Dullinger obtained sampling permits in accordance with national and international laws, and we thank the respective authorities for granting these permits (permits are available upon request). We also thank W. Thuiller and his team (LECA, France) for providing advice in statistics. Analyses were partly conducted using the LEO HPC infrastructure of the University of Innsbruck. Finally, we want to thank the custodians of the European Vegetation Archive (EVA) and its member databases for the permission to use plot data for our analysis. Funder Information Declared FWF Austrian Science Fund, https://ror.org/013tf3c58 , P29413 References ↵ Abs , C . ( 1994 ). Populationsökologie als geobotanische Kausalanalyse . Berichte Der Reinhold-Tüxen-Gesellschaft , 6 , 149 – 163 . OpenUrl ↵ Allouche , O. , Tsoar , A. , & Kadmon , R . ( 2006 ). Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS) . Journal of Applied Ecology , 43 ( 6 ), 1223 – 1232 . doi: 10.1111/j.1365-2664.2006.01214.x OpenUrl CrossRef Web of Science ↵ Barratt , C. D. , Bwong , B. A. , Jehle , R. , Liedtke , H. C. , Nagel , P. , Onstein , R. E. , Portik , D. M. , Streicher , J. W. , & Loader , S. P . ( 2018 ). Vanishing refuge? Testing the forest refuge hypothesis in coastal East Africa using genome-wide sequence data for seven amphibians . Molecular Ecology , 27 ( 21 ), 4289 – 4308 . doi: 10.1111/mec.14862 OpenUrl CrossRef ↵ Becker , D. , Verheul , J. , Zickel , M. , & Willmes , C . ( 2015 ). LGM paleoenvironment of Europe—Map [Dataset] . CRC806-Database . doi: 10.5880/SFB806.15 OpenUrl CrossRef ↵ Birks , H. J. B. , & Tinner , W. ( 2016 ). Past forests of Europe . In J. San-Miguel-Ayanz , G. Caudullo , D. de Rigo , T. Houston Durrant , & A. Mauri (Eds.), European Atlas of Forest Tree Species (1st ed., p. e01493b+) . Publications Office of the European Union . doi: 10.2760/233115 OpenUrl CrossRef ↵ Birks , H. J. B. , & Willis , K. J . ( 2008 ). Alpines, trees, and refugia in Europe . Plant Ecology & Diversity , 1 ( 2 ), 147 – 160 . doi: 10.1080/17550870802349146 OpenUrl CrossRef Web of Science ↵ Broennimann , O. , Fitzpatrick , M. C. , Pearman , P. B. , Petitpierre , B. , Pellissier , L. , Yoccoz , N. G. , Thuiller , W. , Fortin , M.-J. , Randin , C. , Zimmermann , N. E. , Graham , C. H. , & Guisan , A . ( 2012 ). Measuring ecological niche overlap from occurrence and spatial environmental data . Global Ecology and Biogeography , 21 ( 4 ), 481 – 497 . doi: 10.1111/j.1466-8238.2011.00698.x OpenUrl CrossRef ↵ Brus , R . ( 2010 ). Growing evidence for the existence of glacial refugia of European beech ( Fagus sylvatica L.) in the south-eastern Alps and north-western Dinaric Alps . Periodicum Biologorum , 112 ( 3 ), 239 – 246 . OpenUrl ↵ Buiteveld , J. , Vendramin , G. G. , Leonardi , S. , Kamer , K. , & Geburek , T . ( 2007 ). Genetic diversity and differentiation in European beech ( Fagus sylvatica L.) stands varying in management history . Forest Ecology and Management , 247 ( 1 ), 98 – 106 . doi: 10.1016/j.foreco.2007.04.018 OpenUrl CrossRef ↵ Burnham , K. P. , & Anderson , D. R . ( 2002 ). Model Selection and Multimodel Inference ( 2nd ed .). Springer . doi: 10.1007/b97636 OpenUrl CrossRef PubMed ↵ Carnicero , P. , Wessely , J. , Moser , D. , Font , X. , Dullinger , S. , & Schönswetter , P . ( 2022 ). Postglacial range expansion of high-elevation plants is restricted by dispersal ability and habitat specialization . Journal of Biogeography , 49 ( 10 ), 1739 – 1752 . doi: 10.1111/jbi.14390 OpenUrl CrossRef PubMed ↵ Catchen , J. , Hohenlohe , P. , Bassham , S. , Amores , A. , & Cresko , W . ( 2013 ). Stacks: An analysis tool set for population genomics . Molecular Ecology , 22 ( 11 ), 3124 – 3140 . doi: 10.1111/mec.12354 OpenUrl CrossRef PubMed Web of Science ↵ Charles , K. L. , Bell , R. C. , Blackburn , D. C. , Burger , M. , Fujita , M. K. , Gvoždík , V. , Jongsma , G. F. M. , Kouete , M. T. , Leaché , A. D. , & Portik , D. M . ( 2018 ). Sky, sea, and forest islands: Diversification in the African leaf-folding frog Afrixalus paradorsalis (Anura: Hyperoliidae) of the Lower Guineo-Congolian rain forest . Journal of Biogeography , 45 ( 8 ), 1781 – 1794 . doi: 10.1111/jbi.13365 OpenUrl CrossRef ↵ Chikhi , L. , Sousa , V. C. , Luisi , P. , Goossens , B. , & Beaumont , M. A . ( 2010 ). The confounding effects of population structure, genetic diversity and the sampling scheme on the detection and quantification of population size changes . Genetics , 186 ( 3 ), 983 – 995 . doi: 10.1534/genetics.110.118661 OpenUrl Abstract / FREE Full Text ↵ Chytry , M. , Danihelka , J. , Kaplan , Z. , Wild , J. , Holubová , D. , Novotný , P. , Řezníčková , M. , Rohn , M. , Dřevojan , P. , Grulich , V. , Klimešová , J. , Leps , J. , Lososová , Z. , Pergl , J. , Sádlo , J. , Smarda , P. , Štěpánková , P. , Tichý , L. , Axmanová , I. , & Pyšek , P. ( 2021 ). Pladias database of the Czech flora and vegetation . Preslia , 93 , 1 – 87 . doi: 10.23855/preslia.2021.001 OpenUrl CrossRef ↵ Chytrý , M. , Hennekens , S. M. , Jiménez-Alfaro , B. , Knollová , I. , Dengler , J. , Jansen , F. , Landucci , F. , Schaminée , J. H. J. , Aćić , S. , Agrillo , E. , Ambarlı , D. , Angelini , P. , Apostolova , I. , Attorre , F. , Berg , C. , Bergmeier , E. , Biurrun , I. , Botta-Dukát , Z. , Brisse , H. , … Yamalov , S . ( 2016 ). European Vegetation Archive (EVA): An integrated database of European vegetation plots . Applied Vegetation Science , 19 ( 1 ), 173 – 180 . doi: 10.1111/avsc.12191 OpenUrl CrossRef ↵ Clark , J. S. , Fastie , C. , Hurtt , G. , Jackson , S. T. , Johnson , C. , King , G. A. , Lewis , M. , Lynch , J. , Pacala , S. , Prentice , C. , Schupp , E. W. , Webb , T. , III , & Wyckoff , P. ( 1998 ). Reid’s paradox of rapid plant migration: Dispersal theory and interpretation of paleoecological records . BioScience , 48 ( 1 ), 13 – 24 . doi: 10.2307/1313224 OpenUrl CrossRef Web of Science Clark , P. U. , Archer , D. , Pollard , D. , Blum , J. D. , Rial , J. A. , Brovkin , V. , Mix , A. C. , Pisias , N. G. , & Roy , M . ( 2006 ). The middle Pleistocene transition: Characteristics, mechanisms, and implications for long-term changes in atmospheric pCO2 . Quaternary Science Reviews , 25 ( 23 ), 3150 – 3184 . doi: 10.1016/j.quascirev.2006.07.008 OpenUrl CrossRef GeoRef PubMed Web of Science ↵ Clement , M. , Posada , D. , & Crandall , K. A . ( 2000 ). TCS: A computer program to estimate gene genealogies . Molecular Ecology , 9 ( 10 ), 1657 – 1659 . doi: 10.1046/j.1365-294x.2000.01020.x OpenUrl CrossRef PubMed Web of Science ↵ Comes , H. P. , & Kadereit , J. W . ( 1998 ). The effect of Quaternary climatic changes on plant distribution and evolution . Trends in Plant Science , 3 ( 11 ), 432 – 438 . doi: 10.1016/S1360-1385(98)01327-2 OpenUrl CrossRef Web of Science ↵ Comps , B. , Gömöry , D. , Letouzey , J. , Thiébaut , B. , & Petit , R. J . ( 2001 ). Diverging trends between heterozygosity and allelic richness during postglacial colonization in the European beech . Genetics , 157 ( 1 ), 389 – 397 . doi: 10.1093/genetics/157.1.389 OpenUrl Abstract / FREE Full Text ↵ Danecek , P. , Auton , A. , Abecasis , G. , Albers , C. A. , Banks , E. , DePristo , M. A. , Handsaker , R. E. , Lunter , G. , Marth , G. T. , Sherry , S. T. , McVean , G. , & Durbin , R . ( 2011 ). The variant call format and VCFtools . Bioinformatics , 27 ( 15 ), 2156 – 2158 . doi: 10.1093/bioinformatics/btr330 OpenUrl CrossRef PubMed Web of Science ↵ Daneck , H. , Fér , T. , & Marhold , K . ( 2016 ). Glacial survival in northern refugia? Phylogeography of the temperate shrub Rosa pendulina L. (Rosaceae): AFLP vs. chloroplast DNA variation . Biological Journal of the Linnean Society , 119 ( 3 ), 704 – 718 . doi: 10.1111/bij.12619 OpenUrl CrossRef ↵ Dawson , M. N . ( 2014 ). Natural experiments and meta-analyses in comparative phylogeography . Journal of Biogeography , 41 ( 1 ), 52 – 65 . doi: 10.1111/jbi.12190 OpenUrl CrossRef ↵ Di Cola , V. , Broennimann , O. , Petitpierre , B. , Breiner , F. T. , D’Amen , M. , Randin , C. , Engler , R. , Pottier , J. , Pio , D. , Dubuis , A. , Pellissier , L. , Mateo , R. G. , Hordijk , W. , Salamin , N. , & Guisan , A. ( 2017 ). ecospat: An R package to support spatial analyses and modeling of species niches and distributions . Ecography , 40 ( 6 ), 774 – 787 . doi: 10.1111/ecog.02671 OpenUrl CrossRef ↵ Doláková , N. , Roszková , A. , & Přichystal , A . ( 2010 ). Palynology and natural environment in the Pannonian to Holocene sediments of the Early Medieval centre Pohansko near Břeclav (Czech Republic) . Journal of Archaeological Science , 37 ( 10 ), 2538 – 2550 . doi: 10.1016/j.jas.2010.05.014 OpenUrl CrossRef Web of Science ↵ Dullinger , S. , Willner , W. , Plutzar , C. , Englisch , T. , Schratt-Ehrendorfer , L. , Moser , D. , Ertl , S. , Essl , F. , & Niklfeld , H . ( 2012 ). Post-glacial migration lag restricts range filling of plants in the European Alps . Global Ecology and Biogeography , 21 ( 8 ), 829 – 840 . doi: 10.1111/j.1466-8238.2011.00732.x OpenUrl CrossRef Web of Science ↵ Ehlers , J. , Gibbard , P. L. , & Hughes , P. D . ( 2011 ). Quaternary Glaciations – Extent and Chronology: A Closer Look ( 1st ed. , Vol. 15 ). Elsevier . ↵ Evanno , G. , Regnaut , S. , & Goudet , J . ( 2005 ). Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study . Molecular Ecology , 14 ( 8 ), 2611 – 2620 . doi: 10.1111/j.1365-294X.2005.02553.x OpenUrl CrossRef PubMed Web of Science ↵ Fordham , D. A. , Saltré , F. , Haythorne , S. , Wigley , T. M. L. , Otto-Bliesner , B. L. , Chan , K. C. , & Brook , B. W . ( 2017 ). PaleoView: A tool for generating continuous climate projections spanning the last 21 000 years at regional and global scales . Ecography , 40 ( 11 ), 1348 – 1358 . doi: 10.1111/ecog.03031 OpenUrl CrossRef ↵ Gavin , D. G. , Fitzpatrick , M. C. , Gugger , P. F. , Heath , K. D. , Rodríguez-Sánchez , F. , Dobrowski , S. Z. , Hampe , A. , Hu , F. S. , Ashcroft , M. B. , Bartlein , P. J. , Blois , J. L. , Carstens , B. C. , Davis , E. B. , de Lafontaine , G. , Edwards , M. E. , Fernandez , M. , Henne , P. D. , Herring , E. M. , Holden , Z. A. , … Williams , J. W. ( 2014 ). Climate refugia: Joint inference from fossil records, species distribution models and phylogeography . New Phytologist , 204 ( 1 ), 37 – 54 . doi: 10.1111/nph.12929 OpenUrl CrossRef PubMed Web of Science ↵ Gent , P. R. , Danabasoglu , G. , Donner , L. J. , Holland , M. M. , Hunke , E. C. , Jayne , S. R. , Lawrence , D. M. , Neale , R. B. , Rasch , P. J. , Vertenstein , M. , Worley , P. H. , Yang , Z.- L. , & Zhang , M . ( 2011 ). The community climate system model version 4 . Journal of Climate , 24 ( 19 ), 4973 – 4991 . doi: 10.1175/2011JCLI4083.1 OpenUrl CrossRef Web of Science ↵ Gubler , M. , Henne , P. D. , Schwörer , C. , Boltshauser-Kaltenrieder , P. , Lotter , A. F. , Brönnimann , S. , & Tinner , W . ( 2018 ). Microclimatic gradients provide evidence for a glacial refugium for temperate trees in a sheltered hilly landscape of Northern Italy . Journal of Biogeography , 45 ( 11 ), 2564 – 2575 . doi: 10.1111/jbi.13426 OpenUrl CrossRef ↵ Gugerli , F. , Brodbeck , S. , Lendvay , B. , Dauphin , B. , Bagnoli , F. , van der Knaap , W. O. , Tinner , W. , Höhn , M. , Vendramin , G. G. , Morales-Molino , C. , & Schwörer , C. ( 2023 ). A range-wide postglacial history of Swiss stone pine based on molecular markers and palaeoecological evidence . Journal of Biogeography , 50 ( 6 ), 1049 – 1062 . doi: 10.1111/jbi.14586 OpenUrl CrossRef ↵ Guido , M. A. , Molinari , C. , Moneta , V. , Branch , N. , Black , S. , Simmonds , M. , Stastney , P. , & Montanari , C . ( 2020 ). Climate and vegetation dynamics of the Northern Apennines (Italy) during the Late Pleistocene and Holocene . Quaternary Science Reviews , 231 , 106206 . doi: 10.1016/j.quascirev.2020.106206 OpenUrl CrossRef ↵ Gutenkunst , R. N. , Hernandez , R. D. , Williamson , S. H. , & Bustamante , C. D . ( 2009 ). Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data . PLoS Genetics , 5 ( 10 ), e1000695 . OpenUrl CrossRef ↵ Herbert , T. D . ( 2023 ). The Mid-Pleistocene Climate Transition . Annual Review of Earth and Planetary Sciences , 51 (Volume 51, 2023), 389 – 418 . doi: 10.1146/annurev-earth-032320-104209 OpenUrl CrossRef ↵ Hewitt , G. M . ( 1999 ). Post-glacial re-colonization of European biota . Biological Journal of the Linnean Society , 68 ( 1–2 ), 87 – 112 . doi: 10.1111/j.1095-8312.1999.tb01160.x OpenUrl CrossRef GeoRef Web of Science ↵ Hewitt , G. M . ( 2000 ). The genetic legacy of the Quaternary ice ages . Nature , 405 ( 6789 ), 907 – 913 . doi: 10.1038/35016000 OpenUrl CrossRef PubMed Web of Science ↵ Hošek , J. , Pokorný , P. , Storch , D. , Kvaček , J. , Havig , J. , Novák , J. , Hájková , P. , Jamrichová , E. , Brengman , L. , Radoměřský , T. , Křížek , M. , Magna , T. , Rapprich , V. , Laufek , F. , Hamilton , T. , Pack , A. , Di Rocco , T. , & Horáček , I. ( 2024 ). Hot spring oases in the periglacial desert as the Last Glacial Maximum refugia for temperate trees in Central Europe . Science Advances , 10 ( 22 ), eado6611 . doi: 10.1126/sciadv.ado6611 OpenUrl CrossRef PubMed ↵ Huber , K. , Weckström , K. , Drescher-Schneider , R. , Knoll , J. , Schmidt , J. , & Schmidt , R . ( 2010 ). Climate changes during the last glacial termination inferred from diatom-based temperatures and pollen in a sediment core from Längsee (Austria) . Journal of Paleolimnology , 43 ( 1 ), 131 – 147 . doi: 10.1007/s10933-009-9322-y OpenUrl CrossRef GeoRef Web of Science ↵ Huson , D. H. , & Bryant , D . ( 2006 ). Application of phylogenetic networks in evolutionary studies . Molecular Biology and Evolution , 23 ( 2 ), 254 – 267 . doi: 10.1093/molbev/msj030 OpenUrl CrossRef PubMed Web of Science ↵ Kaltenrieder , P. , Belis , C. A. , Hofstetter , S. , Ammann , B. , Ravazzi , C. , & Tinner , W . ( 2009 ). Environmental and climatic conditions at a potential glacial refugial site of tree species near the Southern Alpine glaciers. New insights from multiproxy sedimentary studies at Lago della Costa (Euganean Hills, Northeastern Italy) . Quaternary Science Reviews , 28 ( 25 ), 2647 – 2662 . doi: 10.1016/j.quascirev.2009.05.025 OpenUrl CrossRef GeoRef ↵ Karger , D. N. , Conrad , O. , Böhner , J. , Kawohl , T. , Kreft , H. , Soria-Auza , R. W. , Zimmermann , N. E. , Linder , H. P. , & Kessler , M . ( 2017 ). Climatologies at high resolution for the earth’s land surface areas . Scientific Data , 4 ( 1 ), 170122 . doi: 10.1038/sdata.2017.122 OpenUrl CrossRef PubMed ↵ Karger , D. N. , Nobis , M. P. , Normand , S. , Graham , C. H. , & Zimmermann , N. E . ( 2023 ). CHELSA-TraCE21k – high-resolution (1 km) downscaled transient temperature and precipitation data since the Last Glacial Maximum . Climate of the Past , 19 ( 2 ), 439 – 456 . doi: 10.5194/cp-19-439-2023 OpenUrl CrossRef ↵ Kearse , M. , Moir , R. , Wilson , A. , Stones-Havas , S. , Cheung , M. , Sturrock , S. , Buxton , S. , Cooper , A. , Markowitz , S. , Duran , C. , Thierer , T. , Ashton , B. , Meintjes , P. , & Drummond , A . ( 2012 ). Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data . Bioinformatics , 28 ( 12 ), 1647 – 1649 . doi: 10.1093/bioinformatics/bts199 OpenUrl CrossRef PubMed Web of Science ↵ Keppel , G. , Van Niel , K. P. , Wardell-Johnson , G. W. , Yates , C. J. , Byrne , M. , Mucina , L. , Schut , A. G. T. , Hopper , S. D. , & Franklin , S. E. ( 2012 ). Refugia: Identifying and understanding safe havens for biodiversity under climate change . Global Ecology and Biogeography , 21 ( 4 ), 393 – 404 . doi: 10.1111/j.1466-8238.2011.00686.x OpenUrl CrossRef Web of Science ↵ Kirchheimer , B. , Wessely , J. , Gattringer , A. , Hülber , K. , Moser , D. , Schinkel , C. C. F. , Appelhans , M. , Klatt , S. , Caccianiga , M. , Dellinger , A. , Guisan , A. , Kuttner , M. , Lenoir , J. , Maiorano , L. , Nieto-Lugilde , D. , Plutzar , C. , Svenning , J.-C. , Willner , W. , Hörandl , E. , & Dullinger , S . ( 2018 ). Reconstructing geographical parthenogenesis: Effects of niche differentiation and reproductive mode on Holocene range expansion of an alpine plant . Ecology Letters , 21 ( 3 ), 392 – 401 . doi: 10.1111/ele.12908 OpenUrl CrossRef PubMed ↵ Kirschner , P. , Perez , M. F. , Záveská , E. , Sanmartín , I. , Marquer , L. , Schlick-Steiner , B. C. , Alvarez , N. , Steiner , F. M. , & Schönswetter , P . ( 2022 ). Congruent evolutionary responses of European steppe biota to late Quaternary climate change . Nature Communications , 13 ( 1 ), Article 1. doi: 10.1038/s41467-022-29267-8 OpenUrl CrossRef PubMed ↵ Kirschner , P. , Záveská , E. , Gamisch , A. , Hilpold , A. , Trucchi , E. , Paun , O. , Sanmartín , I. , Schlick-Steiner , B. C. , Frajman , B. , Arthofer , W. , Steiner , F. M. , & Schönswetter , P . ( 2020 ). Long-term isolation of European steppe outposts boosts the biome’s conservation value . Nature Communications , 11 ( 1 ), 1968 . doi: 10.1038/s41467-020-15620-2 OpenUrl CrossRef PubMed ↵ Kirschner , P. , Záveská , E. , Hülber , K. , Wessely , J. , Willner , W. , Schönswetter , P. , & Frajman , B . ( 2023 ). Evolutionary dynamics of Euphorbia carniolica suggest a complex Plio–Pleistocene history of understorey species of deciduous forest in southeastern Europe . Molecular Ecology , 32 ( 19 ), 5350 – 5368 . doi: 10.1111/mec.17102 OpenUrl CrossRef ↵ Kolář , F. , Píšová , S. , Záveská , E. , Fér , T. , Weiser , M. , Ehrendorfer , F. , & Suda , J . ( 2015 ). The origin of unique diversity in deglaciated areas: Traces of Pleistocene processes in north-European endemics from the Galium pusillum polyploid complex (Rubiaceae) . Molecular Ecology , 24 ( 6 ), 1311 – 1334 . doi: 10.1111/mec.13110 OpenUrl CrossRef ↵ Kopelman , N. M. , Mayzel , J. , Jakobsson , M. , Rosenberg , N. A. , & Mayrose , I . ( 2015 ). Clumpak: A program for identifying clustering modes and packaging population structure inferences across K . Molecular Ecology Resources , 15 ( 5 ), 1179 – 1191 . doi: 10.1111/1755-0998.12387 OpenUrl CrossRef PubMed ↵ Kropf , M. , Kadereit , J. W. , & Comes , H. P . ( 2002 ). Late Quaternary distributional stasis in the submediterranean mountain plant Anthyllis montana L . ( Fabaceae) inferred from ITS sequences and amplified fragment length polymorphism markers. Molecular Ecology , 11 ( 3 ), 447 – 463 . doi: 10.1046/j.1365-294x.2002.01446.x OpenUrl CrossRef PubMed ↵ Lavergne , S. , Mouquet , N. , Thuiller , W. , & Ronce , O . ( 2010 ). Biodiversity and climate change: Integrating evolutionary and ecological responses of species and communities . Annual Review of Ecology, Evolution, and Systematics , 41 (Volume 41, 2010), 321 – 350 . doi: 10.1146/annurev-ecolsys-102209-144628 OpenUrl CrossRef Web of Science ↵ Leigh , J. W. , & Bryant , D . ( 2015 ). popart: Full-feature software for haplotype network construction . Methods in Ecology and Evolution , 6 ( 9 ), 1110 – 1116 . doi: 10.1111/2041-210X.12410 OpenUrl CrossRef PubMed ↵ Leuschner , C. , & Ellenberg , H . ( 2017 ). Ecology of Central European Non-Forest Vegetation: Coastal to Alpine , Natural to Man-Made Habitats ( 1s t ed.). Springer . doi: 10.1007/978-3-319-43048-5 OpenUrl CrossRef ↵ Luqman , H. , Wegmann , D. , Fior , S. , & Widmer , A . ( 2023 ). Climate-induced range shifts drive adaptive response via spatio-temporal sieving of alleles . Nature Communications , 14 ( 1 ), 1080 . doi: 10.1038/s41467-023-36631-9 OpenUrl CrossRef PubMed ↵ MacArthur , R. H. , & Wilson , E. O. ( 1967 ). The Theory of Island Biogeography (REV-Revised) . Princeton University Press ; JSTOR. http://www.jstor.org/stable/j.ctt19cc1t2 ↵ Magri , D. , Vendramin , G. G. , Comps , B. , Dupanloup , I. , Geburek , T. , Gömöry , D. , Latałowa , M. , Litt , T. , Paule , L. , Roure , J. M. , Tantau , I. , Van Der Knaap , W. O. , Petit , R. J. , & De Beaulieu , J.-L. ( 2006 ). A new scenario for the Quaternary history of European beech populations: Palaeobotanical evidence and genetic consequences . New Phytologist , 171 ( 1 ), 199 – 221 . doi: 10.1111/j.1469-8137.2006.01740.x OpenUrl CrossRef PubMed Web of Science ↵ Magyari , E . ( 2002 ). Holocene biogeography of Fagus sylvatica L. and Carpinus betulus L. in the Carpathian-Alpine region . Folia Historico-Naturalia Musei Matraensis , 26 , 15 – 35 . OpenUrl ↵ Marchesini , A. , Silverj , A. , Torre , S. , Rota-Stabelli , O. , Girardi , M. , Passeri , I. , Fracasso , I. , Sebastiani , F. , & Vernesi , C . ( 2023 ). First genome-wide data from Italian European beech ( Fagus sylvatica L.): Strong and ancient differentiation between Alps and Apennines . PLoS ONE , 18 ( 7 ), e0288986 . doi: 10.1371/journal.pone.0288986 OpenUrl CrossRef PubMed ↵ Mayr , E. ( 1954 ). Change of genetic environment and evolution . In J. Huxley , A. C. Hardy , & E. B. Ford (Eds.), Evolution as a Process . Allen and Unwin . https://www.blackwellpublishing.com/ridley/classictexts/mayr.pdf ↵ Milesi , P. , Kastally , C. , Dauphin , B. , Cervantes , S. , Bagnoli , F. , Budde , K. B. , Cavers , S. , Fady , B. , Faivre-Rampant , P. , González-Martínez , S. C. , Grivet , D. , Gugerli , F. , Jorge , V. , Lesur Kupin , I. , Ojeda , D. I. , Olsson , S. , Opgenoorth , L. , Pinosio , S. , Plomion , C. , … Pyhäjärvi , T . ( 2024 ). Resilience of genetic diversity in forest trees over the Quaternary . Nature Communications , 15 ( 1 ), 8538 . doi: 10.1038/s41467-024-52612-y OpenUrl CrossRef PubMed ↵ Monegato , G. , Ravazzi , C. , Culiberg , M. , Pini , R. , Bavec , M. , Calderoni , G. , Jež , J. , & Perego , R . ( 2015 ). Sedimentary evolution and persistence of open forests between the south-eastern Alpine fringe and the Northern Dinarides during the Last Glacial Maximum . Palaeogeography, Palaeoclimatology, Palaeoecology , 436 , 23 – 40 . doi: 10.1016/j.palaeo.2015.06.025 OpenUrl CrossRef GeoRef ↵ Mráz , P. , & Ronikier , M . ( 2016 ). Biogeography of the Carpathians: Evolutionary and spatial facets of biodiversity . Biological Journal of the Linnean Society , 119 ( 3 ), 528 – 559 . doi: 10.1111/bij.12918 OpenUrl CrossRef ↵ Oksanen , J. , Simpson , G. L. , Blanchet , F. G. , Kindt , R. , Legendre , P. , Minchin , P. R. , O’Hara , R. B. , Solymos , P. , Stevens , M. H. H. , Szoecs , E. , Wagner , H. , Barbour , M. , Bedward , M. , Bolker , B. , Borcard , D. , Carvalho , G. , Chirico , M. , Caceres , M. D. , Durand , S. , … Weedon , J. ( 2022 ). Vegan: Community ecology package . R package version 2.6-4 [R package]. https://CRAN.R-project.org/package=vegan ↵ Ortiz , E. M. ( 2019 ). vcf2phylip v2.0: Convert a VCF matrix into several matrix formats for phylogenetic analysis . 10.5281/zenodo.2540861 (Version 2.0) [Python]. doi: 10.5281/zenodo.2540861 OpenUrl CrossRef ↵ Ossowski , S. , Schneeberger , K. , Lucas-Lledó , J. I. , Warthmann , N. , Clark , R. M. , Shaw , R. G. , Weigel , D. , & Lynch , M . ( 2010 ). The rate and molecular spectrum of spontaneous mutations in Arabidopsis thaliana . Science , 327 ( 5961 ), 92 – 94 . doi: 10.1126/science.1180677 OpenUrl Abstract / FREE Full Text ↵ Pan , D. , Hülber , K. , Willner , W. , & Schneeweiss , G. M . ( 2020 ). An explicit test of Pleistocene survival in peripheral versus nunatak refugia in two high mountain plant species . Molecular Ecology , 29 ( 1 ), 172 – 183 . doi: 10.1111/mec.15316 OpenUrl CrossRef ↵ Paris , J. R. , Stevens , J. R. , & Catchen , J. M . ( 2017 ). Lost in parameter space: A road map for stacks . Methods in Ecology and Evolution , 8 ( 10 ), 1360 – 1373 . doi: 10.1111/2041-210X.12775 OpenUrl CrossRef ↵ Paun , O. , Schönswetter , P. , Winkler , M. , Consortium , I. , & Tribsch , A . ( 2008 ). Historical divergence vs. contemporary gene flow: Evolutionary history of the calcicole Ranunculus alpestris group (Ranunculaceae) in the European Alps and the Carpathians . Molecular Ecology , 17 ( 19 ), 4263 – 4275 . doi: 10.1111/j.1365-294X.2008.03908.x OpenUrl CrossRef PubMed ↵ Paun , O. , Turner , B. , Trucchi , E. , Munzinger , J. , Chase , M. W. , & Samuel , R . ( 2016 ). Processes driving the adaptive radiation of a tropical tree ( Diospyros , Ebenaceae) in New Caledonia, a biodiversity hotspot. Systematic Biology , 65 ( 2 ), 212 – 227 . doi: 10.1093/sysbio/syv076 OpenUrl CrossRef PubMed ↵ Petit , R. J. , Aguinagalde , I. , De Beaulieu , J. L. , Bittkau , C. , Brewer , S. , Cheddadi , R. , Ennos , R. , Fineschi , S. , Grivet , D. , Lascoux , M. , Mohanty , A. , Müller-Starck , G. , Demesure-Musch , B. , Palmé , A. , Martín , J. P. , Rendell , S. , & Vendramin , G. G. ( 2003 ). Glacial refugia: Hotspots but not melting pots of genetic diversity . Science , 300 ( 5625 ), 1563 – 1565 . doi: 10.1126/science.1083264 OpenUrl Abstract / FREE Full Text ↵ Piwczyński , M. , Wyborska , D. , Gołębiewska , J. , & Puchałka , R . ( 2018 ). Phylogenetic positions of seven poorly known species of Ferula (Apiaceae) with remarks on the phylogenetic utility of the plastid trnH-psbA , trnS-trnG , and atpB-rbcL intergenic spacers . Systematics and Biodiversity , 16 ( 5 ), 428 – 440 . doi: 10.1080/14772000.2018.1442374 OpenUrl CrossRef ↵ Portik , D. M. , Leaché , A. D. , Rivera , D. , Barej , M. F. , Burger , M. , Hirschfeld , M. , Rödel , M.-O. , Blackburn , D. C. , & Fujita , M. K . ( 2017 ). Evaluating mechanisms of diversification in a Guineo-Congolian tropical forest frog using demographic model selection . Molecular Ecology , 26 ( 19 ), 5245 – 5263 . doi: 10.1111/mec.14266 OpenUrl CrossRef ↵ Postolache , D. , Popescu , F. , Paule , L. , Ballian , D. , Zhelev , P. , Fărcaş , S. , Paule , J. , & Badea , O . ( 2017 ). Unique postglacial evolution of the hornbeam ( Carpinus betulus L.) in the Carpathians and the Balkan Peninsula revealed by chloroplast DNA . Science of The Total Environment , 599 – 600 , 1493–1502. doi: 10.1016/j.scitotenv.2017.05.062 OpenUrl CrossRef ↵ Pritchard , J. K. , Stephens , M. , & Donnelly , P . ( 2000 ). Inference of population structure using multilocus genotype data . Genetics , 155 ( 2 ), 945 – 959 . doi: 10.1111/j.1471-8286.2007.01758.x OpenUrl Abstract / FREE Full Text ↵ R Core Team . ( 2023 ). R: A language and environment for statistical computing . Vienna, Austria : R foundation for Statistical Computing . Http://www.R-project.org [Computer software]. R Foundation for Statistical Computing . https://www.R-project.org/ ↵ Ravazzi , C. , Badino , F. , Marsetti , D. , Patera , G. , & Reimer , P. J . ( 2012 ). Glacial to paraglacial history and forest recovery in the Oglio glacier system (Italian Alps) between 26 and 15 ka cal BP . Quaternary Science Reviews , 58 , 146 – 161 . doi: 10.1016/j.quascirev.2012.10.017 OpenUrl CrossRef GeoRef ↵ Rochette , N. C. , Rivera-Colón , A. G. , & Catchen , J. M . ( 2019 ). Stacks 2: Analytical methods for paired-end sequencing improve RADseq-based population genomics . Molecular Ecology , 28 ( 21 ), 4737 – 4754 . doi: 10.1111/mec.15253 OpenUrl CrossRef PubMed ↵ Rota , F. , Carnicero , P. , Casazza , G. , Nascimbene , J. , Schönswetter , P. , & Wellstein , C . ( 2024 ). Survival in nunatak and peripheral glacial refugia of three alpine plant species is partly predicted by altitudinal segregation . Molecular Ecology , 33 ( 9 ), e17343 . doi: 10.1111/mec.17343 OpenUrl CrossRef ↵ Rull , V . ( 2010 ). On microrefugia and cryptic refugia . Journal of Biogeography , 37 ( 8 ), 1623 – 1625 . doi: 10.1111/j.1365-2699.2010.02340.x OpenUrl CrossRef ↵ Saltré , F. , Saint-Amant , R. , Gritti , E. S. , Brewer , S. , Gaucherel , C. , Davis , B. A. S. , & Chuine , I . ( 2013 ). Climate or migration: What limited European beech post-glacial colonization? Global Ecology and Biogeography , 22 ( 11 ), 1217 – 1227 . doi: 10.1111/geb.12085 OpenUrl CrossRef ↵ Schneeweiss , G. M. , & Schönswetter , P . ( 2010 ). The wide but disjunct range of the European mountain plant Androsace lactea L. (Primulaceae) reflects Late Pleistocene range fragmentation and post-glacial distributional stasis . Journal of Biogeography , 37 ( 10 ), 2016 – 2025 . doi: 10.1111/j.1365-2699.2010.02350.x OpenUrl CrossRef ↵ Schoener , T. W. ( 1970 ). Nonsynchronous spatial overlap of lizards in patchy habitats . Ecology , 51 ( 3 ), 408 – 418 . doi: 10.2307/1935376 OpenUrl CrossRef Web of Science ↵ Shaw , J. , Lickey , E. B. , Beck , J. T. , Farmer , S. B. , Liu , W. , Miller , J. , Siripun , K. C. , Winder , C. T. , Schilling , E. E. , & Small , R. L . ( 2005 ). The tortoise and the hare II: Relative utility of 21 noncoding chloroplast DNA sequences for phylogenetic analysis . American Journal of Botany , 92 ( 1 ), 142 – 166 . doi: 10.3732/ajb.92.1.142 OpenUrl Abstract / FREE Full Text ↵ Shaw , J. , Lickey , E. B. , Schilling , E. E. , & Small , R. L . ( 2007 ). Comparison of whole chloroplast genome sequences to choose noncoding regions for phylogenetic studies in angiosperms: The tortoise and the hare III . American Journal of Botany , 94 ( 3 ), 275 – 288 . doi: 10.3732/ajb.94.3.275 OpenUrl Abstract / FREE Full Text ↵ Slovák , M. , Kučera , J. , Turis , P. , & Zozomová-Lihová , J . ( 2012 ). Multiple glacial refugia and postglacial colonization routes inferred for a woodland geophyte, Cyclamen purpurascens : Patterns concordant with the Pleistocene history of broadleaved and coniferous tree species . Biological Journal of the Linnean Society , 105 ( 4 ), 741 – 760 . doi: 10.1111/j.1095-8312.2011.01826.x OpenUrl CrossRef Web of Science ↵ Španiel , S. , & Rešetnik , I . ( 2022 ). Plant phylogeography of the Balkan Peninsula: Spatiotemporal patterns and processes . Plant Systematics and Evolution , 308 ( 5 ), 38 . doi: 10.1007/s00606-022-01831-1 OpenUrl CrossRef ↵ Šrámková-Fuxová , G. , Záveská , E. , Kolář , F. , Lučanová , M. , Španiel , S. , & Marhold , K . ( 2017 ). Range-wide genetic structure of Arabidopsis halleri (Brassicaceae): Glacial persistence in multiple refugia and origin of the Northern Hemisphere disjunction . Botanical Journal of the Linnean Society , 185 ( 3 ), 321 – 342 . doi: 10.1093/botlinnean/box064 OpenUrl CrossRef ↵ Stewart , J. R. , & Lister , A. M . ( 2001 ). Cryptic northern refugia and the origins of the modern biota . Trends in Ecology & Evolution , 16 ( 11 ), 608 – 613 . doi: 10.1016/S0169-5347(01)02338-2 OpenUrl CrossRef Web of Science ↵ Svenning , J.-C. , Normand , S. , & Kageyama , M . ( 2008 ). Glacial refugia of temperate trees in Europe: Insights from species distribution modelling . Journal of Ecology , 96 ( 6 ), 1117 – 1127 . doi: 10.1111/j.1365-2745.2008.01422.x OpenUrl CrossRef Web of Science ↵ Svenning , J.-C. , Normand , S. , & Skov , F . ( 2008 ). Postglacial dispersal limitation of widespread forest plant species in nemoral Europe . Ecography , 31 ( 3 ), 316 – 326 . doi: 10.1111/j.0906-7590.2008.05206.x OpenUrl CrossRef Web of Science ↵ Taberlet , P. , Fumagalli , L. , Wust-Saucy , A.-G. , & Cosson , J.-F . ( 1998 ). Comparative phylogeography and postglacial colonization routes in Europe . Molecular Ecology , 7 ( 4 ), 453 – 464 . doi: 10.1046/j.1365-294x.1998.00289.x OpenUrl CrossRef PubMed Web of Science ↵ Tel-Zur , N. , Abbo , S. , Myslabodski , D. , & Mizrahi , Y . ( 1999 ). Modified CTAB procedure for DNA isolation from epiphytic cacti of the genera Hylocereus and Selenicereus (Cactaceae) . Plant Molecular Biology Reporter , 17 , 249 – 254 . doi: 10.1023/A:1007656315275 OpenUrl CrossRef Web of Science ↵ Templeton , A. R . ( 2008 ). The reality and importance of founder speciation in evolution . BioEssays , 30 ( 5 ), 470 – 479 . doi: 10.1002/bies.20745 OpenUrl CrossRef PubMed Web of Science ↵ Thuiller , W. D. , Georgers , D. , Gueguen , M. , Engl , R. , Breiner , F. , Lafourcade , B. , Patin , R. , & Blanchetea , H. ( 2024 ). Biomod2: Ensemble platform for species distribution modeling . R package . Https://biomodhub.github.io/biomod2/ (Version R package version 4.2-6-2) [R]. https://biomodhub.github.io/biomod2/ ↵ Todesco , M. , Pascual , M. A. , Owens , G. L. , Ostevik , K. L. , Moyers , B. T. , Hübner , S. , Heredia , S. M. , Hahn , M. A. , Caseys , C. , Bock , D. G. , & Rieseberg , L. H . ( 2016 ). Hybridization and extinction . Evolutionary Applications , 9 ( 7 ), 892 – 908 . doi: 10.1111/eva.12367 OpenUrl CrossRef PubMed ↵ Tribsch , A. , & Schönswetter , P . ( 2003 ). Patterns of endemism and comparative phylogeography confirm palaeo-environmental evidence for Pleistocene refugia in the Eastern Alps . Taxon , 52 ( 3 ). doi: 10.2307/3647447 OpenUrl CrossRef ↵ Urbaniak , J. , Kwiatkowski , P. , & Ronikier , M . ( 2018 ). Postglacial history and current population genetic diversity of a central-European forest plant Hacquetia epipactis . Preslia , 90 ( 1 ), 39 – 57 . doi: 10.23855/preslia.2018.039 OpenUrl CrossRef ↵ Valtueña , F. J. , Preston , C. D. , & Kadereit , J. W . ( 2012 ). Phylogeography of a Tertiary relict plant, Meconopsis cambrica (Papaveraceae), implies the existence of northern refugia for a temperate herb . Molecular Ecology , 21 ( 6 ), 1423 – 1437 . doi: 10.1111/j.1365-294X.2012.05473.x OpenUrl CrossRef PubMed Web of Science ↵ van der Knaap , W. O. , van Leeuwen , J. F. N. , Finsinger , W. , Gobet , E. , Pini , R. , Schweizer , A. , Valsecchi , V. , & Ammann , B. ( 2005 ). Migration and population expansion of Abies , Fagus , Picea , and Quercus since 15000 years in and across the Alps, based on pollen-percentage threshold values . Quaternary Science Reviews , 24 ( 5 ), 645 – 680 . doi: 10.1016/j.quascirev.2004.06.013 OpenUrl CrossRef GeoRef Web of Science ↵ Wessely , J. , Gattringer , A. , Guillaume , F. , Hülber , K. , Klonner , G. , Moser , D. , & Dullinger , S . ( 2022 ). Climate warming may increase the frequency of cold-adapted haplotypes in alpine plants . Nature Climate Change , 12 ( 1 ), 77 – 82 . doi: 10.1038/s41558-021-01255-8 OpenUrl CrossRef ↵ Wessinger , C. A . ( 2021 ). From pollen dispersal to plant diversification: Genetic consequences of pollination mode . New Phytologist , 229 ( 6 ), 3125 – 3132 . doi: 10.1111/nph.17073 OpenUrl CrossRef PubMed ↵ Widmer , A. , & Lexer , C . ( 2001 ). Glacial refugia: Sanctuaries for allelic richness, but not for gene diversity . Trends in Ecology & Evolution , 16 ( 6 ), 267 – 269 . doi: 10.1016/S0169-5347(01)02163-2 OpenUrl CrossRef PubMed ↵ Willis , K. J. , & van Andel , T. H. ( 2004 ). Trees or no trees? The environments of central and eastern Europe during the Last Glaciation . Quaternary Science Reviews , 23 ( 23 ), 2369 – 2387 . doi: 10.1016/j.quascirev.2004.06.002 OpenUrl CrossRef GeoRef Web of Science ↵ Willner , W. , Di Pietro , R. , & Bergmeier , E. ( 2009 ). Phytogeographical evidence for post-glacial dispersal limitation of European beech forest species . Ecography , 32 ( 6 ), 1011 – 1018 . doi: 10.1111/j.1600-0587.2009.05957.x OpenUrl CrossRef Web of Science ↵ Willner , W. , Wessely , J. , Gattringer , A. , Moser , D. , Záveská , E. , Dullinger , S. , Schönswetter , P. , & Hülber , K . ( 2023 ). Postglacial range formation of temperate forest understorey herbs – insights from a spatiotemporally explicit modelling approach . Global Ecology and Biogeography , 32 , 1046 – 1058 . doi: 10.1111/geb.13677 OpenUrl CrossRef ↵ Záveská , E. , Kirschner , P. , Frajman , B. , Wessely , J. , Willner , W. , Gattringer , A. , Hülber , K. , Lazić , D. , Dobeš , C. , & Schönswetter , P . ( 2021 ). Evidence for glacial refugia of the forest understorey species Helleborus niger (Ranunculaceae) in the Southern as well as in the Northern Limestone Alps . Frontiers in Plant Science , 12 , 839 . doi: 10.3389/fpls.2021.683043 OpenUrl CrossRef View the discussion thread. Back to top Previous Next Posted May 15, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Spatiotemporal diversification of forest understorey species reveals the existence of multiple Pleistocene forest refugia in Central Europe Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Spatiotemporal diversification of forest understorey species reveals the existence of multiple Pleistocene forest refugia in Central Europe Camille Voisin , Philipp Kirschner , Eliška Záveská , Božo Frajman , Karl Hülber , Johannes Wessely , Wolfgang Willner , Peter Schönswetter , Pau Carnicero bioRxiv 2025.05.13.653645; doi: https://doi.org/10.1101/2025.05.13.653645 Share This Article: Copy Citation Tools Spatiotemporal diversification of forest understorey species reveals the existence of multiple Pleistocene forest refugia in Central Europe Camille Voisin , Philipp Kirschner , Eliška Záveská , Božo Frajman , Karl Hülber , Johannes Wessely , Wolfgang Willner , Peter Schönswetter , Pau Carnicero bioRxiv 2025.05.13.653645; doi: https://doi.org/10.1101/2025.05.13.653645 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Evolutionary Biology Subject Areas All Articles Animal Behavior and Cognition (7633) Biochemistry (17680) Bioengineering (13889) Bioinformatics (41927) Biophysics (21445) Cancer Biology (18585) Cell Biology (25491) Clinical Trials (138) Developmental Biology (13373) Ecology (19897) Epidemiology (2067) Evolutionary Biology (24308) Genetics (15606) Genomics (22494) Immunology (17736) Microbiology (40385) Molecular Biology (17175) Neuroscience (88583) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4822) Physiology (7641) Plant Biology (15149) Scientific Communication and Education (2045) Synthetic Biology (4293) Systems Biology (9822) Zoology (2271)