Introduction

In times of a rapidly changing climate, understanding how contemporary climate change affects the European tree flora is of high importance (Lenoir & Svenning 2015). Rapid climate change might cause a spatial homogenization of phylogenetic diversity (Saladin et al. 2020). While the prediction of future tree species distributions is still prone to high uncertainty (Pecchi et al. 2019), understanding the effects of past climatic changes during the Pleistocene (2.6 million to 11,700 years ago) on the spatial variation of phylogenetic diversity will help to support the prediction of biodiversity dynamics. Throughout the Pleistocene, repeated glacial cycles forced temperate tree species to move to glacial refugia, mostly in Mediterranean regions south of the Alps from where they migrated northward during warmer periods. As a consequence of the geographical barriers of the Mediterranean Sea and the Alps, the European tree flora is species-poor compared to the nemoral forests in Northern America or East-Asia (Schröeder 1998), and mainly cold- and drought-tolerant taxa occur (Svenning 2003; Leuschner & Ellenberg 2017).

Hence, the population structure of the present European tree flora is strongly affected by the post glacial colonisation with clear signatures of Mediterranean glacial refugia (Saladin et al. 2020). During the inter-glacial periods, trees expanded poleward from the refugia into suitable habitat, and such recolonisation events resulted in range shifts and expansions. These postglacial demographic expansions leave genetic signatures which are displayed in many temperate species, with postglacial migration routes being established on variation in genetic markers. Individuals colonising the north from adjacent refugia leads to an east–west pattern of genetic structure, which has been found in a number of studies on trees such as oaks, ash, crabapple and grey alder (Petit et al. 2002; Heuertz et al. 2004; Cornille et al. 2013; Mandak et al. 2016). This process of recolonisation also produces another gradient due to post-glacial founder effects and a steady increase in genetic drift along the axis of expansion, resulting in decreased allelic richness in the more recently colonised higher latitude northern portion of the range relative to the refugia. Indeed, recent paleoclimate modelling conducted in the genera Pinus and Quercus strongly supports a Mediterranean refugia hypothesis with a suitable Last Glacial Maximum (LGM) habitat in the Iberian Peninsula, south France, Italy and the Balkans (Di Pasquale et al. 2020). However, in eastern Europe, where there has been persistent ambiguity and uncertainty regarding the location of refugia in the Balkan peninsula, findings from a study of Carpinus betulus supports a multiple refugia-within-refugia scenario and the possible presence of ‘suture zones’, all revealing a complex evolutionary history (Postolache et al. (2017). On the other hand, an increasing body of studies argue for additional northern refugia (Willis et al. 2000; Svenning et al. 2008; Provan & Bennett 2008; Bhagwat & Willis 2008; Parducci et al. 2012; Ali et al. 2020; Saladin et al. 2020). According to this northern refugia hypothesis, given boreal species such as Picea abies L. might have survived the LGM in southern Central Europe refugia, for example such as the Carpathians or even ice-free refugia in Scandinavia, while the majority of nemoral tree species migrated to the Mediterranean region (Birks 2015).

In order to shed further light on the postglacial colonization hypotheses, we selected field maple (Acer campestre L.) as a suited candidate woodland tree species due to its widespread distribution across Eurasia in a transition between the Mediterranean and Euro-Siberian ecoregions including southern parts of Scandinavia. Acer campestre is a medium-sized, subdominant and mesophile nemoral tree species, which has a wide ecological range, although it is more common in mesophile stands, especially deciduous oak forests and has its main distribution area in Central and Southern Europe (Leuschner & Ellenberg 2017; Roloff et al. 2016; Zecchin et al. 2016). Its high social capabilities make A. campestre a characteristic mixed species of the Central and East European broad-leaved forests as it does not form pure stands across its natural range. Usually found in the lower canopy, it has an important role in formation of the vertical structure of xerothermic oak forests and Mediterranean shrub vegetation. In the arid lowlands, A. campestre is often found in a codominant position with oaks; while in humid areas or higher elevations, it has only a limited competitiveness contrary to beech and hornbeam (Nagy 2004). Acer campestre has not been planted widely outside its natural range, except as an ornamental tree. In northern Europe, it has been introduced as a garden plant in hedges and as solitary trees in gardens and parks. In Sweden, the distribution of A. campestre is dominated by garden escapes and foreign gene pools swamp the indigenous population (Wahlsteen 2021a, 2021b). However, given its low commercial importance, A. campestre is not normally silviculturally managed, and it often grows in spontaneously established and semi-natural populations (Zecchin 2016). Therefore, its distribution range and ecogeographical genetic patterns are supposedly nature states (Nagy & Ducci 2004; Chybicki et al. 2014). However, considering the human impact of the European landscape during the Holocene, translocations of plants cannot be ruled out.

Grimm and Denk (2014) outline putative evolution pathways and subsequent dispersal for A. campestre using internal transcribed spacer (ITS) sequences. Major ITS variants were formed in eastern Turkey, Thrace, Bulgaria and Central Europe, with an origin in Turkey. Roughly, the ITS variants can be divided into three main population subdivisions, a north-western, a central and an eastern. The most likely scenario for the dispersal of A. campestre is an origin in Iran, where the species share a common ancestor with Acer cappadocicum Gled. (Li et al. 2019). Acer campestre would probably disperse from the Middle East to west along the northern Mediterranean area. After ice melting, the area of distribution should expand north and reach the UK and Scandinavia at last. Although ITS data frames early dispersals (mainly Neogene), we hypothesise that the population of A. campestre has a south-eastern and a north-western population structure connected by a gradient as a result of pre-glacial dispersal and population expansion after the Last Glacial Maximum.

For A. campestre, we employed approximate Bayesian computation (ABC) to test for different population divergence scenarios in order to propose a hypothesis on glacial refugia impacting the present population structure. Microsatellites have been proved useful in distinguishing factors which influence the demography of populations and can help to reveal phylogeographic signals, often exhibiting genetic patterns which are congruent on those obtained from plastid markers (e.g. Zinck and Rajora 2016). In Fraxinus excelsior and Malus sylvestris, nuclear microsatellites were employed alone to establish hypotheses on post-glacial migration based on approximate Bayesian computation (Heuertz et al. 2004; Cornille et al. 2013).

Besides glacial refugia, we assume that also allele frequency clines along the east to west expansion axis of the species can heavily impact the genetic structure of the A. campestre population. With the expansion, allele frequency clines (AFCs) can be observed, and there is strong evidence theoretically, and more recently empirically, that genetic drift at the expansion edge can significantly influence genetic patterns, producing alleles which can stochastically ‘surf’ on the expansion front, increasing to high frequencies and in some cases reaching large spatial distributions (Excoffier & Ray 2008; Edmonds et al. 2004).

In order to cover major regions of the distributional area of A. campestre, we sampled 61 populations from East Sussex in England to southern Sweden, the greater Caucasus in Azerbaijan and Babor mountains in Algeria, covering an area of 3.3 million km2 in total. The aim and scope of this work is (i) to reveal the genetic structures of A. campestre on a continental level and to put the species in a context of other, earlier investigated forest tree species, (ii) to investigate population diversity clines over the entire distributional range, and (iii) to synthesize the results into a postglacial colonization hypothesis.

Materials and method

Sampling

For the present study, individuals that were present at the same location were regarded as a population. Sixty-one natural populations were sampled, representing 606 individual trees with 10 as the mean number of sampled individuals per population (Table S1). Leaves were collected during summer 2020 from mature trees with some variation in age to avoid sampling of siblings.

DNA extraction and microsatellite genotyping

Twelve microsatellite loci previously described for other closely related Acer species were cross-amplified in A. campestre: MAP 9, MAP 40, MAP 46 (Pandey et al., 2004), Aop116, Aop132, Aop943 (Segarra-Moragues et al. 2008), Am116, AM118, SM14, SM29, and SM60 (Kvesic et al. 2020), and AsT (Harmon et al. 2017). All were PCR-amplified using DFS-“Hot” Taq DNA Polymerase (GeneOn Bioscience, Ludwigshafen am Rhein, Germany) in 1X PCR buffer (20mM Tris-Cl, pH 8.3, 10mM (NH4)2SO4, 35mM KCl, 1.8 mM MgCl2, 0.1mg/ml BSA).

The PCR conditions were as follows: Multiplex 1 (MAP 9, MAP 46, Aop132, AsT) was amplified using a modified touchdown PCR program (according to Chybicki et al. 2014) with initial denaturation at 95 °C for 2 min, five touchdown cycles (94 °C for 30 s, 61 °C (−1 °C/cycle) for 30 s and 72 °C for 1 min) and 28 regular cycles (94 °C for 30 s, 55 °C for 90 s and 72 °C for 1 min; Multiplex 2 (Aop116, Aop943, and MAP 40) was amplified using conventional PCR profile comprising initial denaturation at 95 °C for 2 min and 29 cycles of 94 °C for 30 s, 56.5 °C for 90 s, and 72 °C for 1 min; Multiplex 3 (Am116, AM118, SM14) was amplified using conventional PCR profile comprising initial denaturation at 94 °C for 30 s and 32 cycles of 59 °C for 30 s and 72 °C for 1 min. SM29 and SM60 were amplified individually (36 cycles comprising 94 °C for 30 s, 56 °C for 30 s, and 72 °C for 1 min for SM29, and 36 cycles comprising 94 °C for 30 s, 59.5 °C for 30 s, and 72 °C for 1 min for SM60). The published reverse primers were modified by addition of 5’ GTTTCTT tails to control Taq polymerase terminal transferase adenylation activity (Brownstein et al., 1996). PCR products were diluted 1:15 in deionized formamide and sized via a 3730XL DNA analyser using Genescan 400HDRox size standard (Applied Biosystems, UK). Exported data was analysed using Genotyper 2.0 software (Applied Biosystems, UK).

Genetic diversity and spatiality

Population diversity (number of different alleles, observed and expected heterozygosity, Fst and private alleles) was calculated in GenAlex (Peakall & Smouse 2006), allelic richness with a minimum of 10 individuals in Popgenereport (Adamack & Gruber 2014), Fis and significance of deviation from Fis in Arlequin 3.5.2.2 (Excoffier et al. 2005). Number of alleles, allele size range, observed and expected heterozygosity, Hardy-Weinberg equilibrium and polymorphic information content was calculated in Cervus 3.0.7 and frequency of null alleles was calculated in freeNA (Chapuis & Estoup 2007). To decide whether to keep all loci or remove those with higher null allele frequencies, we compared He and Fst between the two datasets with and without such loci. The allele frequencies were plotted in QGIS (version 3.16.2) to reveal geographical distributions or frequency clines over the distribution range. We then carried out correlation between genetic diversity of populations and latitude to test for declining genetic diversity trends from south to north, which is an expected signature of founder effects along colonisation routes.

For the six alleles exhibiting a geographic cline, a Mantel tests were conducted to reveal any correlation to the five chosen climate variables (see IBE below). Geo-positioned observed allele frequencies from all 61 populations were imported to QGIS; the values were interpolated by the inverse distance weighting method to create a Geo tiff-file. On this data file, a geographic 1° grid was superimposed from which 717 data points were sampled by a sliding window approach similar to that proposed in Pereira et al. (2018). This method compensates for potential stochastic differences between populations, allowing the identification of global patterns of variation throughout the colonization axis (Pereira et al. 2018). The set of interpolated allele frequency data was tested against longitude and latitude in linear regressions to reveal any geographical clines along the expansion axis.

Population genetic structure

To estimate the number of population genetic clusters (K), we used a number of complementary methods: the two model-based approaches structure (2.3.4; Pritchard et al. 2000) and Geneland (4.0.9; Guillot et al. 2005); and also the multivariate non-model based discriminant analysis of principal components (DAPC; Jombart et al. 2010). Structure is a non-spatial method which clusters on genetic information alone, and Geneland is a spatial method which integrates genetic and geographic information, thus allowing the identification of genetic clusters and the influence that spatial clustering has on the population structure. Discriminant analysis of principal components was used to further explore the population structure.

With structure, the posterior distribution was approximated with 1,000,000 iterations after 100,000 samples that were disregarded for burn-in. The admixture model was combined with both the independent and the correlated allele frequencies models. The use of the independent model was justified by several geographically isolated populations with probably low rate of admixture. According to Porras-Hurtado et al. (2013), the correlative model can give structure more power at detecting subtle structures and the authors suggest that it is prudent to use it, as it would identify undetected correlation but without affecting the result if the correlation does not exist. The analysis was performed assuming the number of clusters K = 1 –20, with 20 repetitions for each K. The optimal K was determined based on the Delta K statistic (Evanno et al. 2005) as well as the statistics MedMedK, MedMeanK, MaxMedK and MaxMeanK based on the Peuchmaille method (Puechmaille 2016), which better accounts for uneven datasets, in thresholds varying from 0.5 to 0.8 in steps of 0.1. Both methods were conducted via the web interface StructureSelector (Li & Liu 2018). Based on the structure results, proportion of admixture related to given membership was calculated for all populations (i.e. the remaining proportions of the memberships not inferred to the major group).

For Geneland, the model was implemented in a MCMC scheme, using 100,000 iterations, every 100th iteration saved, K ranging from 1 to 20, and the correlated allele frequency model was used, which allows the detection of structure in the presence of low genetic differentiation. However, Guillot et al. (2012) caution that although the correlative model can be more powerful, it can be prone to algorithm instability and sensitivity to departures from model assumptions (e.g. isolation by distance) and always recommend comparing results with the uncorrelated model.

For DAPC, the chord distance (Cavalli-Sforza & Edwards 1967) between populations was computed using the freeNA software (Chapuis & Estoup 2007). From the distance matrix, a principal components analysis was conducted resulting in 60 PC:s which were kept for K-mean clustering. To confirm the optimal K solutions proposed by the Delta K and the Peuchmaille approaches for the structure result, K was set to 3, 4, and 9 also for the DAPC analysis. Discriminant analysis was subsequently conducted for all 61 populations and the three optimal solutions. All computations were conducted in the software Past 4.03 (Hammer et al. 2001).

Population divergence history

For probabilistic analysis among alternative hypotheses for the history of gene pool divergence, the approximate Bayesian computation procedure was employed as it is implemented in the software diyABC (Cornuet et al. 2014). However, ABC exhibits caveats in low ability to recover true demographic models and to estimate accurate time for population divergence (Cabrera & Palsbøll 2017). Considering this, for the highly admixed population of A. campestre, we employ the method with caution and will interpret the results as general insight into the past demography of the species. From the K3 model derived from the structure analysis with a western, central, and an eastern gene pool, a threshold of 80% affinity to one of these gene pools permitted inclusion in the ABC analysis. In this way, sampled localities with high admixture or doubtful affinity were eliminated to achieve a more reliable result. In total, the eastern gene pool was represented by 81 individuals and the western and the central by 92 each.

For the three gene pools, several initial computations were conducted, all representing three main glacial refugium hypotheses (Figure S1):

Three refugia hypothesis: the three gene pools are derived from three independent glacial refugia, a western, a central and an eastern (scenario 1).

Two refugia hypothesis: Western and eastern refugia are ancestral for two gene pools, and the central gene pool is derived from one of those (scenarios 6 and 7). A western and an eastern refugia are ancestral for two gene pools, and the eastern gene pool is derived from one of those (scenarios 8 and 9). An eastern and a central refugia are ancestral for two gene pools, and the western gene pool is derived from one of those (scenarios 10 and 11). The central gene pool is derived from admixture by the western and the eastern gene pool (scenario 12).

One refugium hypothesis: One western refugium is ancestral for all three gene pools (scenarios 2 and 3). One eastern refugium is ancestral for all three gene pools (scenarios 4 and 5).

Initially, all eleven scenarios were tested simultaneously, and thereafter, a hierarchical approach was employed reducing the eleven scenarios in pairs until a final scenario was revealed. Reducing the number of scenarios simultaneously tested to three, significantly improves the accurate recovery rate (Cabrera & Palsbøll 2017). In diyABC, the genetic model was set to default with the generalized stepwise mutation model, but locus AsT was set to 4 bp motif, summary statistics was set to mean number of alleles, mean genic diversity and mean size variance for one sample and Fst, classification index and (dμ)2 distance for two sample summary statistics. One million datasets were simulated for each scenario in the reference table.

Isolation by distance and environment (IBD, IBE)

The hypothesis of isolation by distance was tested by a Mantel test to test the occurrence of a positive correlation between a genetic matrix (chord distances) and geographic distances. This was conducted on the global dataset and on the optimal K solutions from the Bayesian cluster methods. Mantel test was conducted in PAST 4.03 based on chord distances for genetic data and Euclidean distances for environmental data. Chord distances were computed in freeNA software.

Whether the sampled populations exhibit a genetic connection to climate parameters was tested by Mantel tests. This approach exhibits caveats but has been suggested appropriate for continuous populations or populations interconnected by high rates of gene flow (Frichot et al. 2015). Climate data was downloaded from worldclim.org and the four variables — temperature seasonality (BIO4), mean temperature of the wettest quarter (BIO8), mean temperature of the driest quarter (BIO9) and precipitation of driest month (BIO14) — were used for the subsequent analyses. As shown by several authors (i.e. Di Pasquale et al. 2020 and Pradhan 2016), these four variables have the lowest correlation, and to use only the most differing variables reduces the number of statistical tests and the collinearity among the data.

Climate variable BIO4 (temperature seasonality) is correlated with temperature over the coldest periods of the year (BIO6, BIO11 and BIO7). This climate variable is a measure in percent of temperature change over the course of the year, the larger the percentage, the greater the variability of temperature. Climate variable BIO8 (mean temperature of wettest quarter) is not correlated with any other variables. BIO9 (mean temperature of driest quarter) is related to precipitation and temperature over the year (BIO1, BIO6 and BIO 11). BIO14 (precipitation of driest month) is correlated with a precipitation variable of the driest time of the year (BIO17). The chosen variables may impact the geographical distribution and the adaptation of the drought resistant and heat loving species of A. campestre.

Genetic bottlenecks in populations

The software Bottleneck 1.2.0.2 (Piry et al. 1999) was used to determine genetic bottleneck events for all included populations. The coalescent process was simulated under the two-phase model (TPM) with default settings. The significance of heterozygosity excess or deficiency was assessed by sign test and Wilcoxon test because these are the most powerful and robust when used with less than 20 loci. The test for mode-shift was not considered since it required more than 30 individuals in the sampled population. The software also conducts a standardized differences test with a minimum of 20 polymorphic loci. Since our data matrix included maximum 12 loci, this approach was disregarded.

Results

Genetic diversity and spatiality

Twelve polymorphic loci were used for all subsequent analyses (Table 1). The total number of alleles per locus varied from three (AsT) to 37 (MAP 40) with a total of 179 alleles and an average of 14.9. 12 alleles had a restricted geographic eastern range (Aop132: 158, Aop116: 102, SM29: 285, SM14: 80, 84, 90, 96, AM116: 243, 247, 251, 263, 267; see Figure S2), no alleles were restricted only to the western range, and most alleles had a general distribution over the total distribution.

Table 1 The genetic diversity of 12 loci in Acer campestre
Full size table

On a continental level four alleles (2%) exhibited a pattern similar to a geographical cline and increased with western distribution (longitude; SM14: 100 (Fig. 1) and 112, Aop116: 103 and 116 (Figure S3–5)). On a regional level three alleles showed a similar pattern. In the western range SM14: 114 increased with northern distribution (latitude) and in the eastern range MAP 40: 112 also increased with northern distribution (Figure S6–7).

Fig. 1
figure 1

Increase of allele frequency of SM14: 100 by western distribution suggests allele surfing following the range expansion. Arrows indicate dispersal direction from a west Asian origin (Grimm et al. 2012). The graph (inset) describes the exponential correlation between sample distribution in longitude versus increase in allele frequency, based on 717 samples of interpolated data

Full size image

The average frequency of null alleles across all loci was 6.9% with two loci exceeding 10% (Aop132, 12% and MAP 9R, 15%, Table 1). Presence of null alleles pushed mean population He from 0.55 to 0.57 and Fst from 0.078 to 0.080. Since the number of markers used can be considered to be at the lower range and the impact of null alleles rather low, we decided to keep all loci in downstream analyses. The mean expected heterozygosity (0.67) exceeded observed heterozygosity (0.50) and all loci (except AsT) deviated from Hardy-Weinberg equilibrium (Table 1).

The Mantel tests revealed that the frequency of three alleles (Aop116-116, SM14-100 and SM14-112) correlated with temperature seasonality (BIO4) and one of them (SM14-100) also with mean temperature of wettest quarter (BIO8) and with mean temperature of driest quarter (BIO9). For precipitation of driest month (BIO14), there was no correlation.

Population diversity

The mean number of alleles within populations ranged from 3 (Poland) to 6.25 (Serbia) with a mean of 4.63. Allelic richness measured after rarefaction ranged from 2.18 (Sweden) to 3.30 (Serbia) with a mean of 2.69 (Table S1). Mean Fst among populations was 0.08. The inbreeding coefficient (Fis) ranged from −0.33 (Sweden) to +0.39 (Serbia). However, none of the populations with a negative value deviated significantly from zero when the calculation was permuted, thus, strong inbreeding seems not likely to affect the populations (Table S1). Nevertheless, Fis is sensitive to sample size, and the results are probably affected by low population sample size and Wahlund effects, so presented results in Table S1 should be interpreted with caution. Expected heterozygosity ranged from 0.42 (Czech) to 0.70 (Serbia) (mean 0.57) (Table S1) with highest values in southern Europe (Serbia, Bulgaria, Croatia, Slovenia, Italy, southern France and Spain) and lower values in central Europe. Geographically, allelic richness (after rarefaction) correlated overall with expected heterozygosity with higher variation in southern Europe and lower variation in northern, western and central Europe (Fig. 2A, C). For the three main divisions obtained from the Bayesian cluster analyses (see results below), the western and the eastern divisions showed increasing population diversity loss with increasing north distribution. Private alleles were found in 16 of 61 populations with the highest number in Greece (5) and Serbia (3) (Fig. 2D).

Fig. 2
figure 2

Maps of overall population genetic variation for 61 populations of Acer campestre. A Genetic variation measured as expected heterozygosity increases in southern Europe with its center in Serbia and Hungary. B Inbreeding coefficient (Fis) for all populations. C Genetic variation measured as allelic richness increases in southern Europe with its center in Serbia and Hungary. D Number of private alleles per population. Populations without private alleles are omitted

Full size image

Population genetic structure

The two allele frequency models in the structure analysis gave almost identical results when lower K numbers (K=3) were compared. The difference between the two models in membership coefficient (Q) varied from 0.02 to 0.09 with mean 0.06. When higher K numbers were compared (K=9, derived from the MedK model result, see below), the correlation model tended to split some populations into higher Q values (i.e. higher admixture rate). The Q value varied from 0.02 to 0.05 with a mean of 0.03 (Figure S8). Since the difference between the two models are low, and the correlation model provides a more sensitive and conservative interpretation, all further analysis is based on this model.

The assumption of three subpopulations was found to be the optimal according to the delta K model (Fig. 3). The delta K graph revealed one peak at K3 and no secondary peaks. The populations from France, Denmark and England exhibited a high affinity to a western population subdivision. Average proportions of membership over populations were 66% with 34% admixture from the other population subdivisions. Populations from Poland, Hungary and Serbia revealed a high affinity to an eastern population subdivision with average proportions of membership over populations 63% with 37% admixture from the other population subdivisions. A third, broken subdivision is composed of populations from Sweden, Germany, Hungary, Algeria and Azerbaijan with 63% membership proportion and 37% admixture from other population subdivisions. Populations with the highest membership proportions to the broken division were Sweden (92%), Germany Eibelstadt (85%) and Algeria (87%) with five populations showing membership proportions higher than 80%. For the western division the highest proportions were found in France Migné-Auxances (94%), Denmark Vester Skerninge (86%) and France Eyliac (87%), and six populations showed membership proportions higher than 80%. For eastern division, the highest proportions were Czech Hnanice (90%), Poland Czerniejewo (90%) and Hungary Vásárosmiske (84%), and five populations showed membership proportions higher than 80%.

Fig. 3
figure 3

Structure bar plot (below) and geographical projection of the optimal K according to the Delta K statistic with membership proportions shown as pie charts. The analysis reveals a broken population subdivision (purple), an eastern subdivision (blue) and a western subdivision (orange). The DeltaK graph indicates K3 as optimal with no secondary peaks (upper inset). The most probable population divergence scenario is depicted in lower inset with corresponding coloured arrows in the main figure for the west (W) and the east (E) gene pools, indicating suggested post glacial expansion axes

Full size image

All four estimators of the MedK model were computed for four thresholds, but only MedMeaK and MaxMeak were considered for further analysis. The two latter have been shown superior for populations dispersed via a stepping stone model (Puechmaille 2016), a model corresponding to the IBD results shown below. For threshold, 0.5 MedMeaK = 11 and MaxMeak = 12 were suggested, for threshold 0.6–0.7 MedMeaK = 9–8 and MaxMeaK = 9 and 9 were suggested and for threshold 0.8 MedMeaK = 4 and MaxMeaK = 5 were suggested. Puechmaille (2016) explains that similar results would suggest a clear signal, why we choose K = 9 (from threshold 0.6–0.7) as a contrasting solution to the delta K statistic. This model reveals a high degree of admixture among populations (Figure S9), but a western versus an eastern fraction is still distinguishable. A central genepool extends southwards to Algeria (purple) and several populations of Sicily share a southern genepool (green). The central, purple, genepool has the lowest admixture rate of 58% followed by the western, orange of 65% and the western, pink of 69%.

The Geneland analysis using the uncorrelated model revealed the maximum a posteriori estimate of K = 2 with two clusters consisting of an eastern and a western division (Fig. 4), very similar to the structure K = 3 solution. The simultaneous analysis with the correlated model revealed a broken spatial population structure with related populations split over the continent.

Fig. 4
figure 4

Geneland results for K = 2 using the spatial model with uncorrelated allele frequencies. Pie charts illustrate estimated posterior probability of population membership by posterior mode. Inset, posterior probability isocline, which indicates extensions of the genetic populations

Full size image

The principal components analysis generated 60 PC:s where only the first five components exhibited significant influence on the result according to the broken stick model. However, for further analysis all PC:s were kept. K was set to 3 and 9 to represent the subdivisions suggested by the delta K statistic and the MedK model for the structure analysis. For both suggestions of optimal K, the DAPC analyses overall agree with the results of the structure analysis, dividing the population into two main subdivisions (eastern and western) and several intermediates (Figure S10–11).

Population divergence history

When all twelve scenarios were tested simultaneously, both the direct and the logistic regression suggested the ‘three refugia hypothesis’ (scenario 1) to be the most likely. The posterior distribution of parameter t (divergence time) for scenario 1, was estimated to 1450 generations ago (95% CI = 529–3060). The effective population size of the common ancestral population was estimated to 1090, the eastern to 8250, the western to 3160 and the central to 7590. Also the pairwise reduction ended in scenario 1 with estimated posterior distribution of divergence time to 1340 generations ago (95% CI = 469–2930, Table 2). The effective population size of the ancestral population was 876, the eastern 7920, the western 3030 and the central 7400. Consequently, regardless of sampling strategy, the assumed ancestral population experienced a considerable population expansion at the divergence time. When the generation turnover time for A. campestre is estimated to 20–60 years depending on the habitat (open fields versus forest), the population expansion causing the split can be dated to 29,000 to 87,000 years before present, just before the last glacial maximum (LGM) some 26 500 to 15,000 years BP (Clark et al. 2009).

Table 2 Demographic parameters for the three refugia hypothesis derived from the pairwise scenario set
Full size table

Isolation by distance and environment (IBD, IBE)

Test for isolation by distance (IBD) was conducted at the continental level for all 61 populations and revealed a significant correlation between geographic and genetic distances (p = 0.015, R = 0.37).

When pairwise chord values for the three main divisions obtained from the structure analysis (western, central and eastern, see results below) were subjected to Mantel test against spatial coordinates, the correlation was significant for the eastern division (from Bulgaria versus latitude R=0.54, p=0.009) and the central division (from Germany Randersacker and Gerschheim versus longitude R= 0.80, p=0.043) but not for the western. For the Geneland divisions, a significant correlation was revealed within the western division between Germany and England (R=−0.68, p=0.013), in the Mediterranean division between Sicily and Serbia (R=0.58, p=0.037) and between Sicily and the northern populations (R=0.68, p=0.015). For the eastern division, a south-north gradient was evident between Hungary and Poland (R=0.58, p=0.014).

The Mantel tests between genetic distance and temperature seasonality (isolation by environment, IBE) revealed a weak correlation (R=0.19, p=0.01) as the mean temperature of the wettest quarter (R=0.12, p=0.03). The other three climate variables did not exhibit any correlations.

Genetic bottlenecks in populations

The Wilcoxon test for heterozygosity excess averaged across loci revealed significant results indicating that eight populations have undergone bottleneck events and for the sign, test 6 populations revealed significant results. If both methods are combined, the following nine populations have undergone a bottleneck event: France (Migné-Auxances) Czech Republic (Hnanice), Hungary (three populations) and Republic of Serbia (four populations) (Figure S12).

Discussion

Continental level genetic structure and population diversity clines

We identified three A. campestre gene pools without competing peaks but with high degree of admixture from the others according to the delta K model of the structure results. Although the work of Kvesic et al. (2020) and Chybicki et al. (2014) was employed on a more restricted geography, the high degree of admixture found by these authors is consistent with our results, using the same model of optimal K. When the number of gene pools (K) was increased, in general also the membership proportions increased, indicating that more gene pools did not contribute to better reveal population structures. This is consistent with expectations for population structure for populations experiencing isolation by distance and for species with dispersal models explained by a stepping stone model. The population of A. campestre shows a west to east climate change (Fig. 3) with no isolated gene pools possible for the structure algorithm to detect. This division is consistent with the results from ITS variants presented by Grimm and Denk (2014) revealing a south-eastern and a north-western population structure connected by a gradient. The K=3 model is the most intuitive solution supposed to represent three gene pools derived from three glacial refugia. Considering limitations in the quality and number of markers in combination with limited sampling size, we propose the lower K solution to explain the population geographical structure. Geneland, on the other hand, evaluates the presence of population structure in georeferenced genetic data by inferring and explicitly identifying genetic discontinuities along the landscape. Thus, Geneland reveals genetic divisions neatly united in the landscape (Fig. 4) contrary to the structure result with closely related populations scattered over the distribution range. However, the main patterns of separated west and east gene pools are still evident also in the Geneland result, and we conclude that both algorithms generate similar interpretations of the genetic data.

For the continental level, our results show a significant positive correlation between genetic differentiation and geographic distance indicating that populations are spatially genetically structured, with isolation by distance (IBD) playing an important role. For the genetic subdivisions obtained from the Bayesian methods, the structure divisions (west and east) revealed a clear signal of IBD whence the corresponding divisions for the Geneland result showed no such signal. This may be explained by a mathematical artefact due to geographically larger and more heterogeneous divisions from the structure result or that the results actually represent a biological signal explained by range expansions from different glacial refugia.

When climate factors were tested against genetic differentiation the significance of the results were low and the correlation weaker than for geographic distance. This is expected when using neutral markers, but also indicates a stronger signal for genetic structure impacted by geographic distance rather than varying climate conditions. For Fraxinus excelsior, Heuertz et al. (2004) interpret the significant IBD patterns within southeastern Europe as a result of colonisations from different glacial refugia. For Fagus sylvatica, Jump et al. (2006) concluded that in continuous forests, no spatial genetic structure was evident from isolation by distance. The authors further found population differentiation to be very low, indicating that at the spatial scale of the study, the continuous forest area is subject to panmictic breeding. In contrast, the authors detected significant isolation by distance between forest fragments. The influence of gene flow between fragments therefore declines relative to the effects of genetic drift as the geographical distance between them increases. This pattern indicates that fragmentation of the beech forest in this region has led to the breakup of the panmictic breeding system observed in the neighbouring area of continuous forest. Acer campestre exhibits a fragmented distribution as a secondary forest species and the tree never forms pure, continuous stands. The spatial genetic structure identified as isolation by distance can probably be addressed both by a history of colonisations from different glacial refugia, by the present fragmented distribution due to modern land management and by a natural niche as a secondary forest species.

When drastic differences in allele frequencies between geographic regions often have been interpreted as signatures of positive selection, later simulation studies have shown that neutral surfing during spatial expansion could indeed lead to very similar patterns just by chance (Excoffier and Ray 2008). Our results suggest that the increase of frequency in several alleles are correlated with the direction of a post glacial population expansion or an early dispersal of the species from an ancient west Asian origin (Fig. 1), rather than changing climate conditions. The observation of surfing alleles is congruent with the hypothesised dispersal of the species outlined by Grimm and Denk (2014) based on ITS variants, but it also strengthens the assumption of a western and an eastern glacial refugium (Figure S6–7). We conclude that allele surfing has occurred during population expansion in A. campestre and that this phenomenon has a considerable impact on the genetic structure of the species (Fig. 1, S3–5). The biological meaning of allele surfing can thus be explained as a relict genetic signature present in species experienced a rapid range expansion. Since this signature is found in neutral alleles, it represents geographical rather than adaptational traits.

Postglacial colonization

The geographical range of the three gene pools and the decreasing genetic diversity by increasing latitudes suggest that A. campestre contracted its range to a few glacial refugia, which probably were in the south of France and in the Balkans. We show that the population diversity increased in the south-eastern part of the distribution area (especially in Serbia and Hungary), according to our geographical projection (Fig. 2). This pattern is congruent with the frequently observed south-north gradient in genetic diversity (“southern richness to northern poverty”) on the northern hemisphere is thought to reflect the range expansions induced by the glacial cycles (Hewitt 1996). This trend is often observed at species’ range margins where smaller effective population sizes and stronger geographical isolation make the populations more susceptible to the loss of genetic variation (Tollefsrud 2009). For A. campestre, the differences in population diversity can be addressed by the expansion of the ice during the last glacial period with refugial remnants in southern Europe. During recolonisation, it is likely that genetic diversity got lost due to environmental adaptation and multiple founding events (Le Corre & Kremer 1998; Wegmann et al. 2006 and Dlugosch & Parker 2008). Another signature may be bottlenecks that may occur when a population experience a sudden reduction of size, expected in range rear populations. Our tests for bottleneck events revealed that ten of the 61 populations have experienced a recent sudden reduction in population size. High bottleneck statistics in recently recolonized areas may not only reflect signatures of recent population bottlenecks, but also population admixture. For instance, when two or more populations with identical allele composition but different allele frequencies are mixed, the number of alleles remains the same, but the distribution of allele frequencies becomes more even, producing a relative increase of gene diversity (Heuertz 2004). Our results show that the Swedish population and one of the French (Migné-Auxances) have both particularly low admixture rates and have experienced recent bottlenecks. This signature of bottlenecks may probably be explained by a typical founder event where a limited gene pool of a few individuals expands into a new population. For several of the east European populations, the admixture rate is high, and the signature may be explained by population admixture. However, our results do not confirm that edge populations should be more susceptible to bottlenecks. The results reveal that the populations detected for bottlenecks are scattered over the range of the species with most signatures in the eastern and central regions.

Considering the main findings of the Bayesian cluster analyses and the approximate Bayesian computation (ABC) demography reconstruction, we suggest a post-Pleistocene expansion in A. campestre from three main areas (or gene pools). The estimated divergence time for the three gene pools can be dated to 29,000 to 87,000 years ago, just before the last glacial maximum (LGM). The ABC results suggest a considerable population migration at the time of the end of the Weichselian glaciation when virgin habitats became available. An eastern gene pool was evident from the structure results analysed by both the delta K and the MedK method and from the Geneland analysis. The eastern gene pool exhibits a clear signal for isolation by distance (IBD), and the genetic difference and population diversity loss increases by the distance from the Balkan region. We interpret these results as the Balkan may have acted as a glacial refugia for A. campestre. However, nuclear microsatellites rather express more recent population genetic events, so the results may be interpreted with caution until further work based on plastid data confirms this hypothesis.

Also, a western gene pool was revealed from the both Bayesian analyses with a clear signal of isolation by distance and population diversity loss from a putative refugium area in the Pyrenees-southern France. For the western gene pool, the level of allelic richness was lower than that in the eastern populations (Fig. 2). These differences in genetic variability may be due to much more severe climatic episodes (i.e. arid and cold) during the Quaternary Period in this region than in other parts of Europe (Petit et al. 2002). Tree populations that survived successive ice ages in southwest Europe were restricted to a few small suitable areas and were thus smaller than those in other parts of Europe.

The origin of the third gene pool (the broken subdivision) present in the scattered populations in Algeria, Greece, Azerbaijan and Sweden is more demanding to explain. Translocation of plants or seeds by humans may be an explanation of the interrupted distribution; conscious introduction of field maple into Algeria by colonizers may be likely and the population in Greece may origin from introduced garden plants. The ABC analysis suggested a simultaneous divergence for the three gene pools, opening for a solution explained by a third glacial refugium. Among the sampled localities dominated by the third gene pool most are located in Germany and Hungary, assuming a northern of Alp refugium. The northern refugia hypothesis postulating that several tree species grew north of the Alps in micro-environmentally favourable sites, originally developed by Willis et al. (2000), has resulted in some controversies. For example, while Parducci et al. (2012) provide evidence that pine and spruce trees survived even in an ice-free refugium in northwestern Norway, Tzedakis et al. (2013) claim that temperate trees were absent north of 45°N during the late Pleniglacial, and questions the evidence of northern refugia for Fagus sylvatica. Just recently, however, Alsos et al. (2020) showed that the Andøya lake on the northwest coast of Norway indeed must have been ice free for more than 23,400 years BP, and the authors do not preclude local growth of Picea abies. As a consequence, more recent studies assume the presence of various glacial refugia in Europe besides the main potential Mediterranean refugia (Saladin et al. 2020). In our study, the K3 solution from the structure analysis suggest a strong influence of a third gene pool in Germany and Hungary (Fig. 3, purple), but there is no evidence of IBD or any spatial population diversity loss among these populations and the Geneland analysis included these populations in larger divisions (eastern and western). However, weak statistical signals may be addressed to low power caused by limited sampling in these geographical areas.

The palynological literature on A. campestre is highly limited and reported findings are rather young. Thus, is it difficult to strengthen a northern refugia hypothesis for our species. However, records more than 7200 years old are found in German Brandenburg, west of Berlin (Wolters 1999 & 2002) and records more than 7800 years in the southern North Sea (Wolters et al. 2010). This can be compared to the oldest reported pollen from Bulgaria from 8200 to 9300 years BP (Connor et al. 2013). Our sampling is too limited to suggest an exact location for a north of Alps refugium, and further plastid works is necessary to build a convincing hypothesis. Thus, our present work cannot reliably explain the origin of the third gene pool, but leaves that question open for future work.

Conclusions

Because of its widespread distribution across Eurasia, A. campestre is a suited candidate species to address the postglacial colonization hypotheses. In this species, we identified hierarchical isolation by distance, clines in allele frequencies, gradients of genetic diversity and a population structure, indicating a phylogeographic signal. Altogether, this supports a ‘serial founder events’ model characterised by drift and extreme drift in the form of allele surfing along multiple postglacial recolonisation axes. Acer campestre shows a high degree of admixture among populations and typical signatures of isolation by distance with no naturally delimited subpopulations. The population structure is highly impacted by geographical, rather than climatological factors such as allele frequency clines and alleles strongly limited to geographical areas. Our data also suggest that the population structure still today harbours signatures of post-glacial migrations from several glacial refugia.