Byline: Frederik Saltre, Remi Saint-Amant, Emmanuel S. Gritti, Simon Brewer, Cedric Gaucherel, Basil A. S. Davis, Isabelle Chuine Keywords: Europe; Fagus sylvatica; Holocene; migration modelling; migration rate; post-glacial migration dynamics; process-based SDM Abstract Aim Despite the recent improvements made in species distribution models (SDMs), assessing species' ability to migrate fast enough to track their climate optimum remains a challenge. This study achieves this goal and demonstrates the reliability of a process-based SDM to provide accurate projections by simulating the post-glacial colonization of European beech. Location Europe. Methods We simulated the post-glacial colonization of European beech over the last 12,000 years by coupling a process-based SDM (PHENOFIT) and a new migration model based on Gibbs point processes, both parameterized with modern ecological data. Simulations were compared with palaeoarchives and phylogeographic data on European beech. Results Model predictions are consistent with palaeoarchives and phylogeographic data over the Holocene. The results suggest that post-glacial expansion of European beech was limited by climate on its north-eastern leading edge, while limited by its migration abilities on its north-western leading edge. The results show a mean migration rate of beech varying from 270myr .sub.-1 to 280myr.sub.-1 and a maximum migration rate varying from 560myr.sub.-1 to 630myr.sub.-1, when limited and not limited by climate, respectively. They also highlight the relative contribution of known and suspected glacial refugia in present beech distribution and confirm the results of phylogeographic studies. Main conclusions For the first time, we were able to reproduce accurately the colonization dynamics of European beech during the last 12kyr using a process-based SDM and a migration model, both parameterized with modern ecological data. Our methodology has allowed us to identify the different factors that affected European beech migration during its post-glaciation expansion in different parts of its range. This method shows great potential to help palaeobotanists and phylogeographers locate putative glacial refugia, and to provide accurate projections of beech distribution change in the future. Article Note: Editor: Niklaus Zimmermann Supporting information: Additional Supporting Information may be found in the online version of this article CAPTION(S): Figure S1 Potential present distribution of Fagus sylvatica simulated by PHENOFIT using ATEAM climate data (left panel, AUC=0.79) and observed present distribution according to the European Potential Natural Vegetation Map (right panel). Colours correspond to the different vegetation formation in which beech is recorded as present whatever its status and coverage (http://www.floraweb.de/vegetation/dnld_eurovegmap.html). Figure S2 Conceptual scheme of the Gibbs-based model. The initial spatial configuration of offspring (black dots) is progressively reorganized in function of the parents (black squares) until the final spatial configuration characterizing the given species is reached. Figure S3 Method used to obtain the IPF of tree cohorts in a 25km2 forest stand from the IPF of tree individuals in a 4 km2 forest stand. IPF parameter sets for tree individual are (1.5; -1.4), (6.8; -2.4), (17.9; -8.9), (25.8; -3.4) and (200; -1.4). IPF parameter sets for tree cohorts are (75; -60), (250; -30), (1000; 10), (1500; -80) and (3000; -90). The first value is the distance in meter and the second value is the interaction I[bar]. Figure S4 Conceptual scheme of the IPF calibration. The IPF parameters are fitted to the spatial pattern of the forest stand which is described by the pair correlation function g(r) (Pommerening 2002). For each set of IPF parameters drawn during the optimization process, a point pattern is simulated according to a non-homogeneous Gibbs process and a g(r) function is calculated. This g(r) function is compared to the observed g(r) function. The fitting is achieved using a least square criteria and the simulated annealing algorithm of Metropolis etal. (1953). LSS expresses the least sum of square. Figure S5 The cohort model. The 4 km2 beech forest stand was subdivided regularly (on the left panel) to obtain cohorts defined by their barycentre (middle panel) and mean density represented by the size of the circle symbolizing the cohort (right panel) and obtained with the histogram (middle panel). Figure S6 Boxplot of differences of arrival date estimated by observation [using a I14C datating method, see detail in Brewer (2002)] and dates estimated by our simulation as a function of six colonization scenarios which differed by climatic influence, migration abilities of beech and potential refugia number and location. Table S1 Parameters used in PHENOFIT. Note that the flowering date parameters are identical to the leaf unfolding date parameters except the F* parameter. This is because Fagus sylvatica has compound buds containing leaves and flowers, and flowers appear at the apex of the shoot a few days after leaf unfolding. Prov 201, 403, 602, 751 correspond to four French geographical provenances of beech as defined by the French Forest Inventory. Table S2 Result of the regression analysis of variance. Row 1 to row 3 assess the agreement significance of one factor (climate, migration or refugia) on simulated beech distributions compared with observations (pollen and macrofossils), conditionally to the other two factors holding constant. The row 4 assesses the same agreement significance between model data, considering all of the three factors. The significance (P value) of the F-statistic was tested using a permutation (Nb permutations) of the residuals. ' Var' is the variance explained by each of the three factors.