Multiple cryptic refugia of forest grass Bromus benekenii in Europe as revealed by ISSR fingerprinting and species distribution modelling

Despite not having been fully recognized, the cryptic northern refugia of temperate forest vegetation in Central and Western Europe are one of the most important in the Holocene history of the vegetation on the subcontinent. We have studied a forest grass Bromus benekenii in 39 populations in Central, Western and Southern Europe with the use of PCR-ISSR fingerprinting. The indices of genetic population diversity, multivariate, and Bayesian analyses, supplemented with species distribution modelling have enabled at least three putative cryptic northern refugial areas to be recognized: in Western Europe—the Central and Rhenish Massifs, in Central Europe—the Bohemia–Moravia region and in the Eastern/Western Carpathians. Central Poland is the regional genetic melting-pot where several migratory routes might have met. Southern Poland had a different postglacial history and was under the influence of an Eastern/Western Carpathian cryptic refugium. More forest species should be checked in a west–east gradient in Europe to corroborate the hypothesis on the Western European glacial refugia.


Introduction
A key to understanding the origin of contemporary floras are glacial refugia, i.e. spots where some thermophilous elements could have survived. Certainly, they were located far south of the glacial front, therefore not affected by harsh climate. It is assumed that these refugia, which enabled the recolonization of Central Europe, were mainly in the Mediterranean area. It is widely accepted that the Balkan Peninsula played a significant role in the reconstruction of Central European flora (Taberlet et al. 1998), supplied continuous, at least from the Eemian interglacial, pollen records of trees such as Carpinus sp., Tilia sp., Ulmus sp. and Quercus sp. (Tzedakis et al. 2002).
In recent years, this view was supplemented by the concept of the northern cryptic refugia (Bhgawat and Willis 2008;Provan and Bennett 2008). This hypothesis stated that not only temperate-boreal, but also thermophilous temperate tree species, might have survived the last glacial maximum (LGM, pleniglacial, *21,000 year BP) in refugial areas, ''oases'' in Central Europe, where sufficient warmth and humidity existed in small micro-environmental pockets (Willis et al. 2000). For example, Magri et al. (2006) and Magri (2008) postulated that Fagus sylvatica could survived the last pleniglacial in cryptic northern refugia in the open forests of Central Europe, including the southern Moravia and southern Bohemia and most probably the Eastern Carpathians. Moreover, the macroscopic charcoal records indicate the existence of thermophilous trees such as Carpinus betulus, Quercus sp., Corylus sp., Ulmus sp., and Tilia sp. in the Hungarian landscape during the LGM (Willis et al. 2000). Brewer et al. (2002) recognised primary, full-glacial refugia, and secondary, temporary refugia for Quercus sp., which supported populations of the termophilous trees during the short, climatically unfavourable, late-glacial Younger Dryas stadial. This picture, based on the fossil pollen data, is complemented by chloroplast DNA data (Bordács et al. 2002).
The retracing of the postglacial migration routes of taxa from refugia into the area left by the outgoing glacier is very important for understanding the history of the European Holocene flora. This has been confirmed by many phytogeographical studies on herbaceous plants and animals (Tyler 2002a, b;Kotlík et al. 2006;Stachurska-Swakoń et al. 2012). The palaeobotanical data (studies of macrofossils and pollen analysis) are particularly important in the localization of refugia and in determining the postglacial migration routes of trees and shrubs (Ralska-Jasiewiczowa et al. 2003, 2004Daneck et al. 2011). Unfortunately, there are few comparable studies for forest herbaceous plants because pollen grains are often taxonomically identified only to the genus level.
Some forest species were examined by means of enzyme electrophoresis, including Carex digitata, Melica nutans (Tyler 2002a, b) and Lathyrus vernus (Schiemann et al. 2000), and chloroplast DNA, including Alnus glutinosa , Hedera spp. (Grivet and Petit 2002), Carpinus betulus , and Fraxinus excelsior (Heuertz et al. 2004). All these studies show the importance of the glacial forest refugia in the mountains of southern (the Apennines and Balkans) and western (the Pyrenees) Europe in postglacial recolonization of central and northern Europe.
The molecular analysis of highly variable, non-coding regions of DNA (mini-and microsatellites) has contributed to the identification of past refugia of many European taxa, mainly alpine and subalpine. In the refugial populations an increased genetic divergence, accompanied by a high number of the rare alleles, could be expected from theoretical and empirical studies (Paun et al. 2008;Ronikier 2011), including genetic drift, different selection pressures and mutation (Tyler 2002a). Also, it is expected that the highest population genetic diversity of the forest plant species can be found in places where divergent lineages from separate refugia met in the ''melting pots'' at intermediate latitudes in Central Europe . However, repeated founder events and immigrations from different sources resulting from repeated long-distance dispersal may blur geographical variation patterns (Tyler 2002a).
In this study, we analysed the forest grass Bromus benekenii (Lange) Trimen in 39 populations from Central, Western and Southern Europe to reveal geographic pattern of genetic variation using the PCR-ISSR protocol.
We selected the species for its strict habitat preferences (Kožuharov et al. 1981) and the wide distribution typical to European, temperate forest herb species. It occurs mainly in Central and Eastern Europe, as well as in western Asia (Fig. 1c). Scarce populations of the species occur also in Western Europe (France, Germany, Switzerland, Italy, Austria), southern Scandinavia, the Caucasus and in the mountains of central Asia (Zając and Zając 2009). B. benekenii in Europe is mainly a lowland species, although it may occur also in foothills and the lower mountain areas (Mirek and Piękoś-Mirkowa 2002). It is closely related to fertile beech forest and montane alderwood communities, and it prefers hummocky medium-moist habitats or medium light-exposed forest edges (Balcerkiewicz 2002). It can be also found in dry mesic oak forests on loess, dealpine mountain calcareous grassland, and in the synanthropic vegetation with a high requirement for nutrients or on lime (Š effer et al. 2002;Roleček 2005).
The phylogeographic results were supported by the species distribution modelling (SDM, Guisan and Zimmermann 2000). The SDM techniques are based on statistical or mechanistic approaches to assess the relationship between species distribution and potential determinants with the use of a representative sample of occurrence data from its current geographical range, climatic data and Community Climate System Model (Collins et al. 2004).
In the present study, we performed the genetic analysis of B. benekenii for the following purposes: (1) to test whether the genetic structure of B. benekenii reflects the long-term isolation in putative ''northern'' refugia during the Quaternary; (2) to check a hypothesis on the regional ''melting pots'' where met several migration routs from various refugia; (3) to compare the inferred glacial refugia of B. benekenii with the putative refugial areas of some other forest plant species; (4) to reconstruct the potential distribution of B. benekenii in the LGM using SDM based on the palaeoclimatic scenario and current climatic data to support the main hypotheses.

Sampling
Leaf material of 319 individual plants from 39 populations of B. benekenii was sampled across Western, Central and Southern Europe in Austria, Bosnia and Herzegovina, Czech Republic, France, Germany, Hungary, Italy, Montenegro, Poland, Romania, Slovakia and Slovenia (Fig. 1a, b; Table 1) in July-August 2011. Fragments of leaves were picked out from the individuals distributed in a minimum of 3 m distance from each other and then dried in silica gel. Voucher specimens were deposited in the herbarium of the Institute of Botany, Jagiellonian University in Kraków (KRA).

Species distribution modelling (SDM)
The main source of occurrence data used to models calibration was the GBIF database (http://www.gbif.org) and the Interactive Agricultural Ecological Atlas of Russia and Neighbouring Countries (http://www.agroatlas.ru). Data on the occurrence of species were also supplemented by numerous field studies carried out by the authors across the Europe. After removing duplicate records a total of 627  et al. (1965-1992). The colours of the dots refer to the three ISSR genetic groups resolved using Bayesian analysis (Fig. 2). Size of circle is directly proportional to the rarity index DW (see Table 1 for exact values) Bromus benekenii in Europe 1439  lins et al. 2004) of past climate were used to project potential distribution in the LGM. Bioclimatic variables downscaled to 2 0 30 00 resolution were obtained from the WorldClim dataset (Hijmans et al. 2005, http://www. worldclim.org), together with present-day climate data at the same resolution. To avoid redundancy from a list of all original WorldClim layers we selected variables that were possibly weakly correlated with each other (Pearson's r \ 0.7). Finally, we considered seven environmental variables as potential predictors of the B. benekenii distribution; namely: mean diurnal range of temperature, temperature seasonality, mean temperature of the wettest quarter, mean temperature of the driest quarter, precipitation seasonality, precipitation of the warmest quarter and precipitation of the coldest quarter. All environmental layers were cropped to the same area between 10W-60E longitude and 35N-75N latitude and prepared using GRASS GIS software version 6.4 (http://grass.osgeo.org).
Recent comparative studies showed that different modelling techniques calibrated on the same species can produce different results (Elith et al. 2006), particularly when models are used to project distributions of species into independent climatic scenarios (Thuiller 2004). Therefore, we used nine different algorithms implemented in BIO-MOD software version 1.1-7 (Thuiller et al. 2009  25 -50% 50 -75% >75% Probability of presence: Occurrence data of B. benekeniiused in model calibration and evaluation Coastline in the last glaciation (~21,000 years BP) Current coastline and administrative boundaries The consensus BIOMOD projections of climatic niche of B. benekenii onto current climatic conditions (a) and past climate model CCSM (b), derived by median of models predictions. Darker areas represent higher probability of occurrence. The locations of occurrence records used in calibration and evaluation of models were also marked by white dots in a assigned for training and 30 % for testing each time. As only occurrences data were available, pseudo-absences were generated randomly to fill the absence component of the models. The predictive power of each model was tested using the area under the Receiver-Operating Characteristic (ROC) curve and Cohen's Kappa coefficient (Manel et al. 2001). Models calibrated with collected records of B. benekenii were projected onto current (Fig. 2a) and LGM (*21,000 year BP, Fig. 2b) climatic conditions in Europe. We also performed ensemble forecasting (Thuiller et al. 2009) to generate final consensus models and to identify areas classified as suitable by the majority of the algorithms. This methodological approach allows one to eliminate artifacts generated by individual models and to identify areas classified as suitable by the majority of the SDM algorithms, particularly spatially isolated areas on the northern limits of the B. benekenii geographical range, which could constitute northern cryptic refugia.

DNA extraction and ISSR analysis
The Inter Simple Sequence Repeats (ISSR) method is based on highly polymorphic sequences of satellite DNA, consisting of a number of nucleotide sequences (microsatellite) tandemly repeated in thousands of copies. PCR reaction products are segments of DNA located between microsatellite regions and includes microsatellite sequences (Stepansky et al. 1999). DNA was isolated from fully developed leaves without damage symptoms caused by insects, and mould. DNA was extracted with Genomic Mini AX Plant (A&A Biotechnology). The core sequence of ISSR markers (Table 2) consisted of 2-5 repeats, in total 15-18 nucleotides in whole primer were used in PCR reaction (Stepansky et al. 1999). Amplification was carried out with a 25 ll reaction mixture comprising: a 2.5 ll tenfold concentrated reaction buffer supplied by the Taq DNA polymerase manufacturer (Fermentas), 1.5 mM MgCl 2 , 0.19 mM of each dNTPs (Fermentas), 27 pmol primer, 100 ng template DNA and 1.4 U of Taq polymerase (Taq DNA Polymerase (recombinant), Fermentas). Reactions were conducted with a 2720 thermal cycler (Applied Biosystems). The annealing temperature for primers ISSR2, ISSR4, ISSR7 was 44°C, and for ISSR1, ISSR3, ISSR6 was 47°C. Optimal conditions for the reaction were as follows: initial denaturation: 94°C, 5 min; 42 amplification cycles: denaturation 94°C, 59 s; annealing 44°C (47°C), 59 s; elongation 72°C, 59 s; final elongation 72°C, 7 min.
A negative control reaction without a DNA template was included in each amplification. Products were subjected to electrophoresis in 1.5 % agarose gel stained with ethidium bromide (50 ll/100 ml) at 100 V for about 1.5 h. Bands were observed and archivised with an Imagemaster VDS (Pharmacia, Amersham). Original software Liscap Capture version 1.0 was also applied.
For analysis of band patterns a GelScan version 1.45 (Kucharczyk TE) software was used (http://www.web statsdomain.com/tags/gelscan/). Thanks to the opportunity to create a calibration curve based on the band pattern of markers lengths (GeneRuler TM 100 bp, Fermentas), it was possible to determine the molecular weight of the resulting amplification products. ISSR reproducibility tests (Bonin et al. 2004) included within plate (n = 12) and between plate (n = 9) replicates independently analysed from the DNA extracts.

Data analysis
The amplification products were scored as a presence/ absence matrix of binary data. Percentage of the polymorphic bands (%Pol) and Nei's (1973) gene diversity index (h) were computed in accordance with the Bayesian method based on a non-uniform prior distribution of allele frequencies (Krauss 2000), implemented in an AFLP-surv version 1.0 (Vekemans et al. 2002). Shannon's index of diversity I and total number of fragments (bands) per population FT are sensitive to the number of sampled individuals (Kučera et al. 2008). The rarefaction technique was used to estimate the expected allelic richness at a locus for a fixed sample size. The rarefaction was made using R software (http://www.r-project.org). Rarefied probabilities served as a basis for the calculation of Shannon's index I and Simpson's index S. The Simpson's 1 -S is equivalent to Nei's gene diversity h and represents the probability that two randomly selected alleles in a population are different. We calculated the rarefied allelic richness for two individuals, an unbiased estimator of S index, according to Hurlbert (1971). The mean number of bands per individual in a population FA was calculated and the rarity index DW, corresponding to ''frequency-down-weighted marker values'' per population (Schönswetter and Tribsch 2005) was computed using AFLPdat (Ehrich 2006). The DW marker values were calculated as a number of occurrences of a marker in the population divided by the occurrences of that marker in the entire data set. Finally, these values were summed up and corrected for the total number of markers and individuals (Ehrich 2006). High DW values are expected in long-term isolated populations (Paun et al. 2008).
The correlations between genetic diversity indices and sample size were tested using Pearson's correlation coefficient. The differences in molecular variability among geographic and genetic groups were tested with Kruskal-Wallis H test. The calculations were carried out with STATISTICA 10 (http://www.statsoft.com).
The genetic portioning among 319 individuals in 39 populations of B. benekenii was estimated by means of STRUCTURE, version 2.3.3 (Pritchard et al. 2000), applying a Bayesian model-based clustering algorithm for the use of dominant markers (Falush et al. 2007). The numbers of K = 2-10 groups were tested in five replications per each K. A burn-in period 200,000, followed by 1 million Markov chain Monte Carlo (MCMC) repetitions were used. An admixture model with uncorrelated allele frequencies was applied. The dominant ISSR data were analysed by treating each class of genotypes as being, effectively, haploid alleles, according to the software documentation. The estimation of the optimal number of groups was based on the likelihood of partitions, estimates of posterior probability provided in STRUCTUTRE output, examined as a function of increasing K (Pritchard et al. 2000) and DK values, estimating the change in the likelihood function with respect to K and estimated as an indicator of the most reliable clustering structure (Evanno et al. 2005). A model-based algorithm implemented in the program STRUCTURE was calculated on a mega computer at the Bioportal at the University of Oslo (http://www.bio portal.uio.no). Similarity between runs was estimated using the symmetric similarity coefficient (Nordborg et al. 2005) with Structure-sum R-script (Ehrich 2006).
The genetic structure of populations and variation levels were assessed by an analysis of molecular variance (AM-OVA; Excoffier et al. 1992). Significance levels were determined using 1,023 permutations. The permutations were carried out at three different hierarchic levels: within populations, among populations and among groups of populations (geographical groups of populations based on a priori 10 geographical regions and detected by Bayesian analysis for K = 3), calculated with Arlequin 3.5 (Excoffier and Lischer 2010).
A neighbour-net diagram was based on the population matrix of Nei & Li genetic distances and bootstrapped using 1000 replicates with SPLITSTREE 4.12 software (Huson and Bryant 2006). A principal coordinate analysis (PCoA) was performed using simple matching similarity coefficient (Sokal and Michener 1958) with NTSYSpc version 2.11 (Rohlf 2002). Hallden et al. (1994) considered the simple matching coefficient to be the more appropriate measure of similarity when closely related taxa are considered. The minimum spanning tree among regions was computed with NTSYSpc. The rarefied frequency data were angular transformed according to the formula x 0 = SQRT ARCSIN (x).

Prediction of LGM distributions
The most common method of model evaluation is the threshold-independent Area Under the Curve (AUC) of the receiver operating characteristic (ROC) plot. The AUC values obtained for individual models ranged from 0.86 to 0.95. According to the criteria described by Swets (1988), AUC values between 0.8 and 0.9 are considered to indicate good models, and values between 0.9 and 1-excellent models. Among the models, the RF model achieved the highest mean AUC and Kappa values (0.94 and 0.74, respectively). The AUC scores greater than 0.9 were also achieved by the GAM, GBM, FDA and MARS models. Because of the large number of results, only the final consensus models derived by the median of the models predictions are presented here.
The final consensus projection onto LGM climatic scenario identified the most suitable conditions (probability [0.75) only in the western Alps and Savoy. The areas with a predicted probability of over 0.5 are mainly located in the Rhone valley, Massif Central, in the Pyrenees, Lower Austria, Julian Alps and on the Balkan Peninsula, and along the Dinaric Alps. The majority of the SDM also predicted small, isolated areas with suitable conditions on the northern limits of potential range: in the Rhenish Massif, the Harz Mountains and in the valley of the Danube River, and on the eastern limits, with moderate probability, in several places in the Carpathians (Fig. 2b).
Genetic diversity within-and among-populations PCR-ISSR analysis resulted in 446 polymorphic markers obtained from 319 individuals in 39 populations. The number of bands generated by particular primers per individual varied from 3 to 13, with a mean ranged from 5 to 8 (Table 2). Data quality tests indicated a high repeatability across of the ISSR bands of above 97 %. The individual number of markers ranged from 29 (B34, Slovakia) to 56 (B36, southern Poland) with a mean of 41.7 (5.2 SD). The total (rarefied) number of bands per population Ftr was the lowest (57)  The rarefaction procedure decreased the correlation (Pearson's coefficient) between total number of bands per population and population size from r = 0.85 for raw data to r = 0.52 for rarefied data. The correlation between the Nei's gene diversity h (with AFLPsurv) and Simpson's index of diversity 1 -S equaled r = -0.85. The rarity index DW was not significantly (p [ 0.05) correlated with any of the genetic diversity measures, as well as with the sample size. The rarefied number of bands per population FTr was correlated with 1 -S index of diversity (r = -0.74). The mean number of bands per population was weakly correlated with the rarefied mean number of bands per population (r = 0.52, p [ 0.05). The percentage of polymorphic loci %Pol and Shannon's index I were highly correlated with the rarefied number of bands per population, r = 0.84 and r = 0.89, respectively. All these correlations had a significance level at p \ 0.001.
The Kruskal-Wallis H tests gave significant results for the genetic structuring of 38 (excluding B08) populations among 10 geographic regions for the three indices: individual number of bands FA (p = 0.0223), rarefied number of bands per population FTr (p = 0.0339), and 1 -S (p = 0.0082).

Genetic relationships among regions and STRUCTURE groups
Bayesian groupings of individuals provided a confirmation of the above results. In the STRUCTURE analysis, the optimal DK and the highest similarity coefficient (0.86 ± 0.06 SD) were given only to those runs in which K = 3 (Fig. 1S, Supported Materials). The results showed the existence of the three genetic groups in B. benekenii individuals, irrespective of the geographic population structure (Figs. 1, 3). The genetic Group 1 consisted of the Moravian-Bohemian B01, B02, B08, the German population from Swabian Jura B03, two populations from Slovakia (B33 and B34), and the populations from Hungary (B25, B26, B30) and the Balkans (B31 and B32). The genetic Group 2 consisted of two populations from southern Moravia B06 and B07, the Sudetes (B10, B11, B13), central Poland (B14-B17, B19), and Western Europe (Germany B23, Italy B20 and France B21 and B22). The genetic Group 3 had populations from Germany (B04, B05), central Poland (B15), Romania (B41, B42, B44), southern Poland (B45, B46), and Slovakia (B27).
The Kruskal-Wallis H tests gave significant results for the genetic structuring of 38 (excluding B08) populations among the three genetic groups for the two indices: FA (p = 0.0046) and DW (0.0044). FA index was the highest in the Group 3 (44.7 ± 0.94 SE), and the lowest in the Group 1 (39.4 ± 0.94SE). DW index was the highest in the Group 1 (2.33 ± 0.29 SE), and the lowest in the Group 3 (0.96 ± 0.29 SE).
The results of AMOVA indicated a markedly higher proportion of variation within populations then between populations (F ST = 0.33, Table 3). At the level of ten geographical regions, the variation among groups was 7.90 %, 25.75 % among populations within groups, and 66.35 % within populations, with F ST = 0.33. The three STRUCTRE groups yielded slightly lower variation among groups 6.74 %, 27.85 % within groups, with F ST = 0.34.
The neighbour-net diagram confirmed the three genetic groups among the populations of B. benekenii (Fig. 4). Most of the groups was strongly supported within 87-100 %, and only the two groups were supported in 57 and 58 %. Within  The PCoA ordination analysis at the individual level confirmed the weak structuring of the data set. Individuals from different regions are scattered along Axis 1 and 2, points to their internal heterogeneity (Fig. 5a). Only individuals from three regions: Balkans, Hungary and Romania are relatively genetically homogenous. The first axis delimited mainly the Balkans, Hungary and most of Slovak and Czech populations, and on its right side, southern Poland (B45, B46) and Slovakia (B27). An intermediate position had the Sudetes and a part of the Czech populations. The second and third axes delimited Romania, a part of southern Poland (B36, B38), central Poland (B14) and Germany (B04, B05). An intermediate position along the second axis had most individuals from central Poland and Western Europe. The genetic Group 1 of the STRUCTURE analysis was the most homogenous (Fig. 5b). The genetic Group 3 was divided into the two parts. One of them consisted of southern Poland and Slovakia. The diagram showed a clear genetic STRUCTURE difference within the Czech (CZ) populations that belong to the Genetic groups 1 and 2 (the lower part of the diagram) and similarity between the Slovak (SL) and southern Polish (SPL) populations within the Group 3 (Fig. 5a, b).
A PCoA with a minimum spanning tree calculated at the geographic regional level showed the distinctiveness of the Balkan, Hungarian and Sudetic territory and the existence of two nodes (Fig. 6). The first was formed by Slovakianarea populations connected with the Balkan and southern Poland. The Slovakian node was linked also with the central Polish territory, which formed the second node, here cross links to Romania, Germany, and the Czech Republic. In the same cluster Germany was connected with Western Europe (Italy and Austria).  Table 1, the geographic regions as in Fig. 3. Bootstrap values C50 are given. Group colours correspond to those in Fig. 3 Table 3 Analysis of molecular variance (AMOVA) of Bromus benekenii for 39 populations, ten geographical regions, and three genetic groups

Putative melting pots
Postglacial migrations of European biota are one of the most important issues in the historical biogeography of the Quaternary, as they shaped the present vegetation and geographical distribution of most animal and plant species.
DNA variation, both for tree and non-tree species, is a very informative way to investigate postglacial colonization. In the case of B. benekenii, it should be acknowledge as a wind pollinating forest grass, having dispersed by epizoochory. We found three genetic groups of the populations concordant with the regional geographic level. We hypothesised that they have various phylogeographic status including putative glacial refugia and melting pots, i.e. regions where postglacial migratory routs met (Taberlet et al. 1998). The status of a region could be inferred from the a priori expectations or assumptions. For example, it is expected the generally high genetic richness of the populations in European contact zones . Such a pattern was found for example in Quercus in Scandinavia, where the marginal northern occurrence of the species was independently colonised from two varying directions , and in the present paper in B. benekenii in central Poland (B17, B19) and Germany (B23). Those three ''northern'' populations are scattered among other having low values of genetic richness within 50-52°of latitude. All they were probably outside refugial areas and represent newly established Holocene populations with relatively high values of genetic indices I, FA and DW. The tentative explanation could their origin by hybridization between different genetic lineages that originated from different refugia (Walter and Epperson (2001). The first refugium was located probably in W Europe (Genetic Group 2) and the second in the Carpathians area (Genetic Group 3). Otherwise, it seems a rule that in refugial areas the occurrence of rare, ''endemic'' private alleles is high (Paun et al. 2008), as for example in a subalpine Aconitum bucovinense at the range margin in the Eastern Carpathians (Boroń et al. 2011).

Putative glacial refugia
The numerical analyses of the populations of B. benekenii in Central Europe clearly distinguished the three genetic groups assembled from various regions. They likely represent the different glacial-interglacial migration histories. In addition, the genetic/geographic groups are differentiated in respect of the genetic diversity indices and their distribution fits to the results of climatic reconstructions. It makes their interpretation more plausible and enables the generation of relevant phylogeographic hypotheses.
The most conspicuous findings concern the genetic Group 1, presumably refugial, consisted of the Bohemian-Moravian populations B01, B02, B08 and the German population B03, sharing seven bands. The group is genetically homogenous, supported in 100 %, and possesses one of the highest values of the rarity DW index (2.597-4.479). A high value of rarity index DW is often related to the refugial status of a population. The genetic Group 1, including the three Bohemian-Moravian populations plus German population B03, has the highest (and statistically significantly different) value of the DW index, in comparison to the remaining groups. It suggests the existence of a massive cryptic northern forest grass refugium in the Moravian Highlands.
One of the earliest record of the Holocene beech expansion in Central Europe comes from sites located in the Moravian Highlands (southern Moravia) and southern Bohemia, where beech was sporadically present c. 10,200 year BP and attained 2 % isopollen a 1,000 years later (Magri et al. 2006). Then, a rapid expansion of the species towards the western border of the central Bohemian basin and towards the Western Carpathians 6-5 kyr 14 C bp was noted (Magri 2008). According to Latałowa et al. (2004) Fagus could have reached southwestern Poland via slow migration from the Czech and Moravia lands into the Sudetes and to the Western Carpathians. Thus, most probably B. benekenii could have survived in the southern Moravian-southern Bohemian refugium and then expanded, through the Moravian Gate, to the Sudetes and central Poland. The postulated migration and colonization of new areas was accompanied by repeated bottlenecks that typically lead to reduced genetic diversity northward, a ''leading edge'' model (Hewitt 1999). Such a phenomenon was found for Fagus in Poland (Sułkowska et al. 2012). Also, the glacial Moravian refugium (Č eskomoravské mezihoří Hills, Č eskomoravska vrchovina Hills) was postulated for a forest shrub Lonicera nigra (Daneck et al. 2011). Also, the German population B03 from the Swabian Jura seems also to be refugial, however, more populations from the region should be analysed to make a statement more conclusive.
The next subgroup within the genetic Group 1 was formed by the populations from Slovakia (B34), Pannonia/ Carpathians (B25, B26, and B30, Hungary), and the Balkan countries (B31 and B32), with no common bands. It is a very uniform genetic group that seems to represent also refugial populations. It is characterized by high values of the genetic diversity indices, especially FTr index. Generally, the Balkans, Carpathians and adjoining Carpathian Basin were shelter for Quercus sp. (Bordács et al. 2002). According to Hendrych and Hendrychová (1979) the migrations from the southern areas were one of the most important ways by which the forest flora of Slovakia was enriched in the postglacial period. However, we cannot rule out that the forest grass could have survived the LGM in the Slovak Carpathians, or more south, in the Pannonian Basin.
The genetic Group 2 group consists of the two subgroups. The first, supported in 87 %, (PCoA) is formed by the populations: B07 (S Moravia), B10, B11, B13 (the Sudetes), and B14 (central Poland), sharing nine bands. These populations, opposite to the previous group, had low values of genetic diversity. Their phylogeographic status is unclear. Most probably they represent secondary populations originated from an unrecognized forest refugium.
The second subgroup consists of two populations B17 and B19 from central Poland, four populations B18, B20, B21, B22 from Western Europe and one population from the Rhenish Massif (B23, Germany). The populations from central Poland, the Rhenish Massif and central France shared five bands. We tend to speculate that genetic richness of the German B23, and possibly the French B21, as well as B18 from Austria and B20 from N Italy, is the effect of their refugial status (FA 38.7-47.1, DW 1.341-1.834). In our opinion they represent unrecognized forest refugia in Western Europe, including the Rhenisch and Central Massifs (see CCSM model, Fig. 2b). The high genetic richness of the two central Polish populations B17 (FA = 41.2, DW = 2.247) and B19 (45.3, 1.502, respectively) is clearly not the effect of their refugial status but rather a consequence of the admixture of divergent lines (a melting pot). First, B. benekenii could have migrated in the Holocene from the putative refugia north of the Alps, e.g. in Rhenish Massif (B23, Germany) eastwards. This direction roughly fits the postglacial migration for example of Carpinus betulus onto the Central European Lowlands (Ralska-Jasiewiczowa et al. 2003). Secondly, the forest species could have survived the LGM in the Eastern Carpathians and then migrated northwards to the Polish Highlands (see below) and, sporadically, to central-eastern Poland, for example the population B19. Such a southnorth migratory route in Poland was found for Picea abies and Quercus sp. (Dering et al. 2008).
The Massif Central is known as a refugial from other phylogeographic/palaeobotanic studies including Abies (Magri et al. 2006) and the tall-herb Cicerbita alpina (Michl et al. 2010).
The genetic Group 3 was formed by populations from different geographical regions. Here could be found populations from central and southern Poland, Romania, Germany and Western Europe. It seems a transitory group and represents presumably the expansion of the forest grass species from refugial areas or represents wide distribution of the similar genotype that could be the result of a previously continuous presence (Kotlík et al. 2006;Fér et al. 2007).
The first subgroup, highly supported (99 %), consists of two populations from central and southern Germany B04 and B05 and one population from central Poland B15. These populations have one the highest number of bands per individual FA. These index well-characterized regions were repeated long-distance dispersal from different sources (melting pot) could have occurred. This is an interesting group because it turns our attention to the possible migrations of the species in European lowlands north of the Alps. None of the populations could be shown as a source for the remaining ones because all of them have the same status in respect of their genetic diversity. In a study on AFLP genetic structure of Carpinus betulus in Europe a general trend of decreasing gene diversity FTr and %Pol was detected from east to west, i.e. away from the unrecognized Eastern European refugia (Coart et al. 2005).
The existence of these two genetic and geographic subgroups clearly shows that the Polish populations from central and southern parts of the country had different postglacial histories and sources of origin. As the population from the central part of Poland has originated probably from the Moravian and W European regions, forming a melting pot, the Carpathian populations in southern Poland had sources in the Romanian E Carpathians and Slovakian W Carpathians. Such a pattern seems typical, taking into consideration the distinctness of the Carpathians as a refugial area for many alpine and subalpine species , Ronikier 2011), Fagus (Liepelt et al. 2009), Lonicera (Daneck et al. 2011), Polygonatum (Kramp et al. 2009), and forest animals (Provan and Bennett 2008). Jankovská et al. (2002) report the sporadic presence of some climatically more demanding trees such as Corylus, Ulmus, Quercus and Carpinus under a very cold continental climate during the LGM in the Slovakian W Carpathians. Similar results were obtained for other thermophilous organisms, including the forest rodent bank vole Clethrionomys glareolus (Kotlík et al. 2006) and a shrub Rosa pendulina L. (Fér et al. 2007).

Species distribution modelling
The results of species distribution modelling confirmed some results of PCR-ISSR molecular analysis. Postulated LGM refuges (inferred from the presence of rare alleles) of populations in the Balkans and in southern France were predicted with probability [0.5. Some of the results obtained by Svenning et al. (2008) using different LGM climate simulations (LMDZHR and S3P) and SDM methods confirm the high potential diversity of temperate species in these areas. They suggest that the view of the LGM landscape in Europe as largely treeless needs to be revised. Interestingly, most of the SDM performed on the CCSM climate reconstruction also identified suitable conditions for small scattered areas, e.g. in the Rhenish Massif (from which the B23 population comes), in the Harz Mountains and in the Danube valley (Fig. 2b). Putative refuges of B. benekenii in the Romanian Eastern Carpathians and Slovak Western Carpathians were also predicted, albeit with a lower probability (0.25-0.5). The projection onto CCSM climate reconstruction predicted the highest probability of species presence ([0.75) only in small areas in the northern-western Alps. This result is somewhat surprising since there is evidence that this area was under an ice cap (Kelly et al. 2004). It may be due to the fact that the palaeoclimatic scenario used in these studies was downscaled from coarse-grained outputs of the climatic scenarios provided by general circulation model. Statistical downscaling usually combines macroclimatic anomalies between past and current conditions with current highresolution climate data. Unfortunately, there is no general agreement about the best way of downscaling these simulations to properly define the relationship between largeand local-scale climate variables (Varela et al. 2011). Generally, these results suggest a strong relationship of B. benekenii occurrence to the mountainous areas and foothills in LGM. Longitudinal arrangement of the main mountain ranges in Europe may promote the survival of temperate species, particularly on slopes with a southerly exposure; moreover higher precipitation in higher altitudes could support preference of forest species for mountain areas. While comparing the present results with those of Leroy and Arpe (2007, Fig . 6b), it is clear that the distribution of the thermophilous deciduous trees (i.e. for those whose limiting factor is a minimum temperature [2.5°C) is more geographically restricted and confined to a few locations in Southern and Western Europe. These locations include about 10 regions in the Iberian Peninsula, the Riviera, large parts of the Apennines and Ancona, Dalmatia, isolated pots in Greece and Thracia, as well as in southeastern Europe in portions of Turkey, including the east cost of the Black Sea.

Concluding remarks
In conclusion, the forest grass B. benekenii had probably multiple refugia in Western Europe (Massif Central, Rhenish Massif) and in Central Europe, including the Moravian-Bohemian area, E and W Carpathians. The refugia in Southern Europe in the Balkan Mountains were also predicted from the SDM modelling. The putative migratory routes of the forest grass species agree with the well-documented postglacial migrations of Carpinus and Fagus. The relationships between the central Polish and German/France populations seem to be linked with the migration of Carpinus from the western European refugia. However, the reverse migration from east to west along the European Lowlands is also highly probable from the putative forest refugia on the Podolian-Volhynian Highlands of the Ukraine (Szafer and Zarzycki 1972) or from the Eastern Carpathians Magri et al. 2006;Liepelt et al. 2009;Ilnicki et al. 2011), forming typical melting pot. The migration of the grass to central Poland from the south-westerly direction, via the Moravian Gate, is probably linked with the Moravian-Bohemian refugium, well documented for Fagus (Magri et al. 2006;Magri 2008). From here the grass species could have migrated both to Poland as well as to Germany. Nevertheless, existence of a cryptic refugium in Germany cannot be ruled out and it could serve as source for some the central Polish populations.
In Poland the genetic structure of the populations was quite different in its southern vs. central parts. In southern Poland the influence of the Carpathians is obvious, including their both E and W parts. In central Poland postglacial migrations of B. benekenii converged and formed a mosaic of the genetic groups, forming typical melting pot. The results obtained encourages one to seek the forest glacial refugia in Central Europe with the use of other forest herbaceous species, to better formulate the hypotheses regarding their occurrence in the region.