History or demography? Determining the drivers of genetic variation in North American plants.

Understanding the impact of historical and demographic processes on genetic variation is essential for devising conservation strategies and predicting responses to climate change. Recolonization after Pleistocene glaciations is expected to leave distinct genetic signatures, characterised by lower genetic diversity in previously glaciated regions. Populations' positions within species ranges also shape genetic variation, following the central-marginal paradigm dictating that peripheral populations are depauperate, sparse and isolated. However, the general applicability of these patterns and relative importance of historical and demographic factors remains unknown. Here, we analysed the distribution of genetic variation in 91 native species of North American plants by coupling microsatellite data and species distribution modelling. We tested the contributions of historical climatic shifts and the central-marginal hypothesis on genetic diversity and structure on the whole data set and across subsets based on taxonomic groups and growth forms. Decreased diversity was found with increased distance from potential glacial refugia, coinciding with the expected make-up of postglacially colonised localities. At the range periphery, lower genetic diversity, higher inbreeding levels and genetic differentiation were reported, following the assumptions of the central-marginal hypothesis. History and demography were found to have approximately equal importance in shaping genetic variation.

INTRODUCTION

Disentangling the determinants shaping species’ distributions and genetic diversity is key to understand and conserve biodiversity and ultimately predict responses to ongoing global change. The identification of the spatial distribution of genetic variation can elucidate the mechanisms of evolution and speciation, shed light on the processes maintaining geographical ranges, improve climate change forecasting, anticipate the spread of invasive species, and pinpoint conservation‐priority populations (Guo, 2014; Hampe & Petit, 2005; Howes & Lougheed, 2008). Phylogeographic analyses can be employed to study the role of ecological factors and mechanistic processes, such as past climatic shifts and demographic fluctuations, on the genetic structure of populations. Yet, despite a large body of research on the population genetics of individual species, few studies have looked at general patterns of genetic structure and variation across multiple species’ ranges at large spatial scales (Gaston, 2009).

Historical climate‐driven changes in species range limits are well known to still affect present‐day genetic diversity (Durka, 1999; Hampe & Petit, 2005; Hewitt, 2000). The Pleistocene ice ages led to southwards range contractions, brought by the extinction of northern populations when temperatures plummeted, with multiple areas across continents serving as glacial refugia for species during the ice ages (Beatty & Provan, 2011; Hewitt, 2000). The end of the Pleistocene led to a subsequent northward expansion of many species in the wake of deglaciation (Hewitt, 1996). The historical climate changes over the past 20,000 years have led to shifting, fragmentation, and reconnection of species’ ranges, which have resulted in a legacy of genetic consequences for contemporary populations (Alvarez et al., 2009; Comes & Kadereit, 1998; Taberlet et al., 1998). One major consequence is that previously glaciated areas are expected to show reduced genetic diversity as a result of sequential founder effects during the post‐glacial recolonisation from refugia (Comes & Kadereit, 1998; Schonswetter et al., 2005; Taberlet et al., 1998). The effect of the ice ages was especially strong in temperate areas. In particular, the northern part of North America was covered with two immense ice sheets (Cordilleran and Laurentide; Beatty & Provan, 2010; Taberlet et al., 1998) at the Last Glacial Maximum (LGM). The changing conditions following the melting of the glacial ice make North America an ideal natural laboratory to study patterns of post‐glacial colonisation and geographic variation (Hewitt, 2000; Pulgarín‐R & Burg, 2012).

The distribution of genetic variation within a species is expected to also be shaped by demographic and evolutionary processes at range margins. Declining environmental suitability towards the periphery is predicted to result in decreasing population density (Kirkpatrick & Barton, 1997), thus reducing population size, gene flow and connectivity. Consequently, marginal populations tend to exhibit low genetic diversity and high genetic differentiation (Brussard, 1984; Pfeifer et al., 2009; Sagarin & Gaines, 2002). However, the broad applicability of this biogeographic paradigm – the central‐marginal hypothesis– is debated since patterns of population genetic variation across large spatial scales are highly variable and usually species‐specific (Sagarin & Gaines, 2002; Yakimowski & Eckert, 2007). For example, Eckert et al. (2008) criticised studies for not including a proper quantitative measure of centrality/peripherality or estimates of population sizes. Additionally, the central‐marginal hypothesis assumes concordance between the geographic and environmental spaces, but this assumption might not always hold, since ecological marginality does not always imply spatial peripherality, and vice versa (Pironon et al., 2015; Soule, 1973).

There are relatively few studies that have tested both the central‐marginal hypothesis and historical influences in a phylogeographical framework, or have attempted to distinguish historical effects on genetic diversity from patterns caused by contemporary geographical variation in population demography and dispersion (Eckert et al., 2008; Gaston, 2009). This is problematic since the patterns in genetic diversity resulting from these two processes can resemble each other. If populations at the northern margin show reduced genetic diversity, is this due to founder effects of postglacial recolonization or due to demographic effects related to the central‐marginal hypothesis? Moreover, almost all evidence regarding these fields comes from studies on single species, which makes it difficult to draw any general conclusions.

Here, we use a novel phylogeographic framework to test the contributions of both historical climatic shifts and the central‐marginal hypothesis on population genetic diversity across the ranges of 91 vascular plant species. We do this by coupling genetic data sourced from the literature and species distribution modelling to analyse the spatial structuring of population genetic variation across the North American continent (Figure 1). Species distribution models were hereby employed as macroecological tools, used to estimate population suitability, act as surrogates of abundance, and develop proxies for colonisation history and population demography. We perform a continental‐scale analysis to identify concordant patterns of population genetic diversity on a large number of unrelated taxa, empirically overcoming the drawbacks commonly associated with demographic and colonisation history studies.

FIGURE 1

Overview of the methodological approach, which couples genetic data and species distribution modelling to test the contributions of historical climatic shifts and the central‐marginal hypothesis on the spatial distribution of genetic variation. The distribution data shown and modelling output are for Tsuga canadensis (see Tables S1 and S2)

MATERIALS AND METHODS

Genetic data

We compiled a genetic database consisting of microsatellite data for 91 native diploid species of North American angiosperms and gymnosperms, with each taxon containing at least three sampling sites including a minimum of five individuals per location. We searched the literature for studies employing microsatellite markers on native vascular plants from North America, which resulted in genetic data published in 67 peer‐reviewed studies (see Table S1). Two species were represented by two studies. We also obtained four unpublished data sets by direct communication with the authors. We were unable to find microsatellite data sets of Arctic species, which are expected to display the strongest genetic consequences from Pleistocene climate shifts given the longer postglacial colonisation routes (Hewitt, 2000). For all studies we tried to obtain the raw allelic data by downloading them from data repositories or by contacting the authors. In total, we obtained the raw data for 52 species. The included species display a wide range of ecological characteristics, abundances, range sizes and life history traits, thus aiming to serve as a representative sample of the continent's vascular flora. The varying range areas included aims to reduce potential biases in the emerging patterns (Supporting Information 1). The sampling and genotyping protocols vary across species, being detailed in the corresponding publications.

Genetic summary statistics

We used the microsatellite data to calculate a set of four summary statistics per population, indicative of population genetic diversity and differentiation: the expected and observed heterozygosity (H S and H O, respectively), the population inbreeding coefficient (F IS), and β (Nei, 1987; Weir & Cockerham, 1984). The expected within‐population diversity H S is calculated only based on the allele frequencies within a population, whereas H O considers the observed frequency of heterozygotes in the population. F IS can be defined by comparing the above‐mentioned heterozygosity measures: (H S–H O)/H S, with high levels indicating an excess of homozygous genotypes compared to Hardy‐Weinberg expectations. β is a population‐specific estimator of the genetic differentiation statistic F ST (Nei, 1987; Weir & Cockerham, 1984). The summary statistics calculations were performed using the function basic.stats() from the hierfstat R package (Goudet, 2005), with every sampling location being treated as a separate population. Clone correction was performed for species with apomictic or vegetative reproduction, by removing all duplicated genotypes to ensure each genotype is represented by a single individual (one ramet per genet) per population. For the species for which the allelic data was not available, we used the published estimates of the summary statistics from the original studies, which were mostly limited to H S, H O and F IS.

Species distribution modelling

Species distribution models were built for all 91 species to evaluate habitat suitability and develop proxies for colonisation history and population demography. The georeferenced occurrences from the sites included in the genetic data set were complemented by downloading species records from the Global Biodiversity Information Facility (GBIF) data portal (GBIF dois in Table S1). Occurrences were restricted to North America, and duplicate locations were removed to avoid pseudoreplication. The records were manually cleaned by removing any suspicious data points in uninhabitable locations for each species. All species were represented by at least 10 unique presence records, which is the recommended minimum number required to calibrate a species distribution model (SDM; Proosdij et al., 2016). Species occurrences totalled 97,074 unique records, ranging from 10 to 9594 occurrences per species.

To model the species distributions, the 19 bioclimatic predictors of the WorldClim v.1.4 data set (http://www.worldclim.org/) were obtained for past and present scenarios at a 2.5 arc‐minute spatial resolution (Hijmans & Elith, 2012). The used paleoclimate data resulted from simulations from a global climate model (GCM) for the LGM (approximately 22,000 years ago), as estimated by the MIROC‐ESM climate model (Watanabe et al., 2011). Current conditions represent interpolations of observed climatic data from 1960–1990. The GTOPO30, a global 30 arc‐second digital elevation model (DEM) was retrieved from the USGS EROS archive (https://www.usgs.gov/centers/eros). The GTOPO30 was aggregated by a cell factor of five to achieve a 2.5 arc‐minute resolution, resulting in a DEM for the present scenario. In order to obtain a DEM for the LGM, when the sea level was 120–135 m lower than at present (Clark & Mix, 2002), the GEBCO_2019 Grid, a global 15 arc‐second bathymetry DEM (GEBCO, 2019) was downloaded. The GEBCO_2019 Grid was rasterised and aggregated by a cell factor of 10 to achieve a 2.5 arc‐minute resolution, followed by clipping it using an LGM bioclimatic layer as a mask. Employing altitude as a variable for modelling is advised against when the SDM aims to project to past climatic conditions (Raes & Aguirre‐Gutiérrez, 2018). Thus, slope and aspect (the direction of the slope) were derived from the DEMs and included as variables instead. All bioclimatic and topographic variables were clipped to the extent of the North American continent, from the southernmost point of Panama to the northernmost point of Canada, excluding Greenland. Data layer manipulations were performed with ArcGIS (ESRI).

To avoid collinearity amongst environmental predictors, which can result in overfitting (Graham, 2003; Peterson et al., 2007), the number of predictor variables was reduced by removing highly correlated parameters, as given by a Spearman's rank correlation test (r s > 0.7). When deciding which of two correlated variables to retain, we aimed to closely capture the key determinants of physiological processes limiting distributions of plants, considering the ample range in ecological preferences displayed by the species included. The retained variables for modelling were mean diurnal temperature range (Bio2), temperature annual range (Bio7), mean temperature of wettest quarter (Bio8), annual precipitation (Bio12), precipitation seasonality (Bio15), aspect, and slope.

Species distribution models were built for each species under present conditions and then projected onto the LGM conditions, employing the sdm R library (Naimi & Araújo, 2016). Three modelling methods were implemented: Domain (Carpenter et al., 1993), Generalised Linear Model (GLM; McCullagh & Nelder, 1989), and Maximum Entropy (MaxEnt; Phillips et al., 2006), each belonging to one of the three main types of modelling algorithms, being “profile”, “regression”, and “machine learning”, respectively. Cross‐validation was performed to validate each model, with 70% of the data being employed for calibration and 30% for evaluation, with 10 bootstrap replications being run per method. SDM accuracy was evaluated using the area under the curve (AUC) of the receiver operating characteristic (ROC) plot (Hanley & McNeil, 1982), a threshold‐independent measure that is relatively insensitive to prevalence (McPherson et al., 2004). AUC values range from zero to one; values close to one indicate maximum fit, whereas values under 0.5 (half of the area under the ROC curve) indicate the model prediction is no better than a random prediction. All models performed better than random expectation, with the lowest AUC value being 0.6. Ensembles of model forecasts were then fitted by combining the three modelling techniques. By employing a consensus, errors (sensitivity to data, lack of absence data, errors in environmental variables) tend to cancel each other out in ensembles, thus producing a more robust and conservative solution (Araújo & New, 2007; Diniz‐Filho et al., 2009). The ensembles were built using weighted averaging based on the AUC statistic to cope with model variability and to improve the reliability of model predictions. The R script used for the species distribution modelling is available in the supplement.

Ecological data

We used the output of the SDM to calculate for every population in each species four measures that quantify ecological suitability, colonisation history, and population demography: (1) habitat suitability under the current conditions, (2) habitat suitability during the LGM, (3) distance to range edge under current conditions, and (4) distance to potential glacial refugium. The habitat suitability was taken directly from the ensemble forecasting produced as output of the different SDMs (Anderson & Martínez‐Meyer, 2004; Diniz‐Filho et al., 2009). Thus, the ecological suitability S of each population under the two modelled time frames was defined as the average value of occurrence provided by each model in the ensemble. For estimating population centrality/peripherality, an innovative quantitative measure was developed. For this, the suitability data was transformed into presence/absence data by setting a threshold, using the Max SSS approach of Liu et al. (2005), which maximises the sum of model sensitivity and specificity. Out of 13 threshold selection methods, Max SSS was found to perform best when only presence data is available (Liu et al., 2013). The distance D to the species’ range edge under present bioclimatic conditions was then computed by calculating the closest distance between each population and the contour of the generated binary presence map. Distances to range edge had negative values if populations were found outside the predicted species range. Finally, distance from a potential glacial refugium was calculated as the closest distance from each population to a suitable area under LGM after creating a binary presence map based on the ensemble prediction for the LGM. The development of this proxy enables the identification of previously unknown glacial refugia, particularly relevant for species with northern refugia that a latitudinal proxy would not account for. Distance to range edge computations were performed in ArcGIS version 10.2 (ESRI) and distance to suitable areas under different scenarios were calculated using the geosphere R package (Hijmans et al., 2019), both employing a geodesic method.

Hypothesis testing

Linear mixed effect models were used to test the relative contribution of distance to edge range, distance to a suitable area under the LGM, and present and LGM suitability in shaping the four population genetic parameters. Additional linear mixed effects models were used to investigate the relationship between the distance to edge range and distance to a suitable area under the LGM with present and LGM suitability. The species name was used as a nested random effect in the analyses to account for phylogenetic nonindependence. In addition to analysing the data set as a whole, we also performed the analysis on subsets of the data, based on either the two main taxonomic groups (Gymnosperms and Angiosperms) or the growth forms of the species (trees, shrubs, and herbs).

The linear mixed effects models were performed using the lme4 package (Bates et al., 2013) in R version 3.5. Model selection was performed in terms of parsimony (based on AIC) and variance explained. The variance explained was calculated using the methods proposed by Nakagawa and Schielzeth (2013) as implemented in the MuMIn package, which provides the total variance explained by fixed and random effects and allows the calculation of variance explained by each fixed effect (Barton, 2011). p‐values were calculated using the Satterthwaite (1946) approximations, using the lmerTest R package (Kuznetsova et al.,2017), standard Bonferroni correction was then applied (Bonferroni, 1936). In addition, linear models were run in order to add regression lines to the plots, which aid the visual interpretation of the results of the linear mixed effect models.

RESULTS

The database we compiled includes 1406 populations across 91 vascular plant species, spanning the whole North American continent, except for the arctic regions. These included 885 angiosperms and 521 gymnosperm populations; with respect to life history, the populations could be classified into 512 herbs, 78 shrubs and 816 trees. For 829 populations, the genetic summary statistics were computed using the raw allelic microsatellite data, whereas for the remaining populations the available published estimates were taken (Table S2). Overall minimum and maximum estimates were 0.00 and 0.92 for the expected heterozygosity (H S), 0.00 and 1.00 for the observed heterozygosity (H O), −1.00 and 1.00 for the inbreeding coefficient (F IS), and −0.27 and 1.00 for the population differentiation statistic β.

The species distribution modelling performed on all 91 species had high predictive power for the relationship between the species’ distribution and the bioclimatic variables (Table S1). Of the 2730 species distribution models (three SDMs × 10 bootstrap replications × 91 species), 88% had AUC values above 0.85. AUC values of the MaxEnt models were generally higher than those of GLM models, which in turn were higher than those of Domain, pointing towards a higher predictive power of MaxEnt over GLM and Domain for most species.

The analyses combining the genetic summary statistics for all species with the output of the SDM revealed a legacy of past glaciation in the genetic data. The genetic diversity –as measured with H S– significantly decreased with increasing distance to suitable areas during the LGM (Figure 2). Furthermore, there was a significant positive relationship between the suitability of populations under LGM conditions and both H S and H O (Figure 3). No significant patterns were found for F IS and β, (Figures 2 and 3) even though F IS showed a slight decrease with increasing distance to suitable areas under the LGM. Additionally, populations displayed significantly lower ecological suitability when further away from suitable areas under the LGM (Figure S1).

FIGURE 2

Relationship between the distance from populations to a suitable area under LGM conditions and the four genetic parameters. Significance was tested using linear mixed effects models (Table 1). Regression lines of were added to the figures to aid visual inspection, with the linear model of H e and F IS showing p < .001

FIGURE 3

Relationship between the suitability of populations under LGM conditions and the four genetic parameters. Significance was tested using linear mixed effects models (Table 1). Regression lines were added to the figures to aid visual inspection, with the linear model of H e showing p < .001

TABLE 1

Output of linear mixed model analysis: degrees of freedom (d.f.), percentages of variance explained (variance) and p‐values are given for the global models of each genetic summary statistic including the distance to suitable area under the LGM, the suitability under the LGM, the distance to range edge, and present suitability. The genetic statistics are the expected and observed heterozygosity (H S and H O, respectively), the population inbreeding coefficient (F IS), and beta (β). The total variance explained by the global and minimal models of each genetic summary statistic are reported in the bottom rows. The total variance explained by the historical and demographic proxies is given in the rightmost column, calculated as the sum of the standardised variances explained for the four genetic summary statistics. For each summary statistic, the variables included in the minimum model, according to AIC, are given in bold; variables excluded from the minimum model are given in black

		H S	H O	F IS	β	Total variance of parameters (%)
Distance to suitable area under LGM	d.f.	1310	1198	1110	822	22.1
	Variance (%)	18	1	1.1	2
	p‐value	0.0026	0.3690	1.0000	1.0000
Suitability under LGM	d.f.	1287	1185	1090	817	27.5
	Variance (%)	10.8	2.2	0	14.5
	p‐value	0.0030	0.0464	0.8088	0.3846
Distance to range edge	d.f.	1323	1217	1126	814	22.6
	Variance (%)	14.8	2.1	1.9	3.8
	p‐value	0.0007	1.0000	1.0000	1.0000
Present suitability	d.f.	1327	1209	1114	819	19.1
	Variance (%)	1.3	4.5	3.6	9.7
	p‐value	0.0149	0.3472	1.0000	0.1691
Total variance of global model (%)		79.4	88.2	85.3	62.2
Total variance of minimal model (%)		79.4	88.2	–	63.4

When the full data set was analysed, it also showed clear support for the central‐marginal hypothesis. Marginal populations indeed had a lower ecological suitability as evidenced by a significant positive relationship between distance to range edge and suitability under the current climatic conditions (Figure S2). This was accompanied by a significant increase in genetic diversity (H S) as distance to the range edge incremented (Figure 4). Furthermore, the value of H S showed a significant negative correlation with the present suitability of the populations (Figure 5). No significant patterns were observed for H O, β and F IS (Figure 5). While 13 species displayed negative distance values, these did not bias the observed patterns as the tests remained significant when negative values were removed.

FIGURE 4

Relationship between the distance from populations to the range edge and the four genetic parameters. Significance was tested using linear mixed effects models (Table 1). Regression lines were added to the figures to aid visual inspection, with the linear models of H s and e showing p < .001

FIGURE 5

Relationship between population present ecological suitability and the four genetic parameters. Significance was tested using linear mixed effects models (Table 1). Regression lines were added to the figures to aid visual inspection, with the linear models of H s, H e and β showing p < .01

The standardised variance explained in the global models of the genetic parameters totalled 22.1% by the distance to suitable area under the LGM, 27.5% by LGM suitability, 22.6% by the distance to edge range, and 19.1% by present suitability (Table 1). For the expected heterozygosity, the minimal model included both effects of the ice ages and effects of the central‐marginal hypothesis; this minimal model was composed of the distance to edge range, the distance to a suitable area under the LGM, and the present suitability, with fixed and random effects explaining a total variance of 79.4%. The observed heterozygosity's minimal model was given by the present suitability, accounting for 88.2% of the data variation. Model selection was not performed for the inbreeding coefficient, as no variables were deemed significant. Finally, the minimal model for genetic differentiation was given by present suitability, explaining 63.4% of the variance. For all models, the reported values for the percentage of variance explained also includes the nested random effect of species name, showing that there are strong differences between species in genetic diversity and differentiation.

When the data was divided into subsets based on taxonomic groups or life forms, the results were more complex but still showed support for both the influence of the ice ages and central‐marginal processes (Table S3). When the gymnosperms were analysed separately, both H S and H O decreased with increasing distance to suitable area during the LGM, and increased with habitat suitability during the LGM, indicating influence of the ice ages. The present range edge and suitability were not included in the minimum model according to AIC, though the distance to the present range edge did show a significant P‐value. For angiosperms, H S was similarly correlated with distance to suitable area and habitat suitability during the LGM, but H O was influenced more by the distance to the present range edge and the present suitability, showing the influence of central‐marginal processes. For the growth forms, trees mostly showed correlations with the present range edge (H S and β), while herbs showed correlations only with the present suitability (H S, H O, and β). Shrubs did not show any significant correlations, probably due to the very small number of shrub species in the data set.

DISCUSSION

Here, we used a phylogeographical framework combining species distribution modelling and population genetics to assess the influences of Pleistocene climate shifts and central‐marginal processes on contemporary genetic patterns. Our analysis of genetic data of 91 vascular plant species clearly shows the effects of both Pleistocene climatic events and central‐marginal processes in spatially structuring genetic variation in North America. Of the more than 1400 included populations, those located at species’ range margins and most distant from glacial refugia were found to have significantly reduced genetic diversity. Our multispecies approach grants increased statistical power compared to studies exploring these processes at the single‐species level, allowing us to draw general conclusions at a large geographical scale.

The effects of the Pleistocene ice ages are evident in the decrease in genetic diversity, as measured using H S, with increased distance from suitable areas under the LGM, i.e. glacial refugia. This reduction in diversity is similar to what is observed in simulation models on the genetic make‐up of postglacially colonised regions (Hewitt, 1999, 2000; Pironon et al., 2015). The observed genetic patterns are most likely due to the larger population sizes and more stable population dynamics at refugial localities, as well as repeated bottlenecks at the advancing (generally northern) range edge during postglacial colonisation (Comes & Kadereit, 1998; Hewitt, 1996). Our results therefore indicate that the shifting and fragmentation of species’ geographical ranges in the past 20,000 years has played an important role in moulding genetic variation across large taxonomic and spatial scales (King & Ferris, 1998; Parks et al., 1994; Sewell et al., 1996; Soltis et al., 1997). Nonetheless, we must be cautious when interpreting diversity patterns, as patterns do not solely emerge from simple models of postglacial colonisation (Petit et al., 2003), and the relationship with the LGM variables was no longer significant when the data was divided by growth form.

The effects of central‐marginal dynamics are evident in the negative correlation of genetic diversity with both the distance from the range edge and ecological suitability. These effects were visible both in the full data set and for the angiosperms, but interestingly not for the gymnosperms. The central‐marginal hypothesis predicts range limits arise because peripheral populations occur in marginal habitats and cannot adapt to conditions beyond the range edges (Haag & Ebert, 2004; Kawecki, 2008). This demographic instability can induce low effective population sizes and frequent bottlenecks, leading to reduced genetic diversity in margins, such as observed in our data (Micheletti & Storfer, 2015; Sagarin & Gaines, 2002; Sexton et al., 2009). Although we do not have population size estimates, we assessed the ecological suitability of populations within their ranges. Diniz‐Filho et al. (2009) propose employing this variable as a macroecological surrogate of abundance, which ensures the logical application of the central‐marginal model to complex spatial abundance patterns as it considers species’ ecological, and not geographical, ranges. The applicability of this approach in our data was supported by a clear central‐peripheral pattern that was observed in the ecological suitability for the analysed populations as this was negatively correlated with the distance to the range margins.

Though we found strong patterns in genetic diversity (H S and to some extent H O) related to the proxies, no significant associations were found for the inbreeding coefficient F IS; the minimum model in our linear mixed model analysis of F IS did not include any of the explanatory variables. Theoretical models have suggested that self‐fertilisation may be favoured at the range margins, especially at postglacially expanding range fronts, given the possible advantages of reproductive assurance or local adaptation (Arnaud‐Haond et al., 2006; Hargreaves & Eckert, 2014). Indeed, several studies on single species have shown selfing to be more frequent in marginal populations, due to low population sizes and environmental stress at range margins, leading to elevated inbreeding levels (Aldrich & Hamrick, 1998; Barrett, 2002; Schoen et al., 1996). Our results across a taxonomically wide range of species suggests that this may not be a general phenomenon.

Similar to F IS, no significant relationships were found in the analysis of the full data set for the genetic differentiation (β, a single population estimate of F ST) of populations with the rest of the metapopulation, though the minimum model according to AIC did include the present suitability as an explanatory variable. This lack of a significant relationship contrasts with numerous theoretical models and empirical studies reporting a decrease in genetic differentiation along postglacial expansion routes (Austerlitz et al., 1997; Excoffier et al., 2009). On the other hand, marginal populations are expected to be more isolated from the central ones, so it is possible that the effects of central‐marginal processes and the Pleistocene recolonisation are counteracting each other. It is also possible that major population structure lies elsewhere, such as among different genetic clusters associated with recolonisation from multiple isolated refugia (Ursenbacher et al., 2015).

The influence of the ice ages and the central‐marginal effects varied depending on the life history and taxonomy of the species, as shown by the analyses of the data subsets. Long‐lived perennials such as trees are known to generally show higher diversity and weaker population structure than short‐lived perennials and annuals (Hamrick & Godt, 1996; Nybom, 2004), generally indicating larger population sizes and higher connectivity. This is consistent with our results, as a significant relationship between β and the present suitability was found for herbs, but not for shrubs or trees. Additionally, this can also explain why this same correlation was not significant for gymnosperms –which are mainly trees– but did show significance for angiosperms –which also include herbs. Similarly, the lack of support for central‐marginal dynamics in gymnosperms may be an effect of the higher population connectivity of trees, which may quickly erode any effects of smaller population sizes at the distribution margins. However, it is also important to realise that the subsetting of the data has led to a decrease in statistical power, which may have affected the inferences, especially for the shrubs.

We employed a species distribution modelling approach to develop effective proxies for colonisation history and population demography, and as such the followed approach does have several limitations and makes several assumptions with regard to the ecology of the study species. One assumption is that the species distribution models only base their predictions on the included environmental variables. We used the same set of variables for all species in order to allow for better comparisons among species, but not all variables may be equally applicable to all species. A general concern around species distribution models is that they do not account for biotic interactions, potential for rapid adaptations or time lag (Godsoe & Harmon, 2012; Hällfors et al., 2015). It is also important to acknowledge that both the SDM and the modelled climate data for the LGM represent extreme extrapolations, and therefore can only be taken as a rough approximation of the situation during the ice ages. For a better overview of the locations of glacial refugia, SDM results should be combined with palynological data, though this has also been proven to be challenging due to the coarseness of the palynological record (Birks, 2019). In addition, the approach assumes that all species have reached their equilibrium distribution. However, it is possible that some species have not yet filled their potential ranges in response to rapid postglacial warming, resulting in incomplete SDMs and inaccurate estimations of the distances to glacial refugia. This would result in an underestimation of the emerging genetic patterns, which would be expected to grow stronger during future recolonization dynamics.

Despite the inherent assumptions in the distribution modelling, we believe that our approach is a major step up from many previous studies that have worked with the simple assumption that latitude is an adequate surrogate for either recolonisation history or central‐marginal processes (Eckert et al., 2008). Instead, our SDM‐based approach allowed us to test and separate the effects of both these phenomena. Our results also show that the distribution models, which had high power for predicting the present distribution, can be used as a valuable tool to infer causation of extant genetic variation patterns across a large scale. Notably, the main purpose of our modelling was not to get an accurate prediction of actual habitats, but rather to get an estimate of the overall range of the species, for which these types of models are widely used.

Both Pleistocene historical and demographical processes appear to have equal and nonexclusive importance in shaping genetic variation. Most genetic summary statistics displayed the expected responses to these proxies; however, only a quarter of the comparisons were significant and the variation explained by the fixed effect factors was relatively low. Most variation in the summary statistics comes from within‐species determinants, as a large proportion of the variance was explained by the species random effects. This was not unexpected given the widely different demographic histories of the species studied, and it highlights the role of potential additional variables acting at an intraspecific level. Numerous studies have emphasized how the genetics of populations is shaped by ecological factors, topography, anthropogenic factors –such as land use and habitat fragmentation–, biotic interactions, and life‐history traits –including the aforementioned breeding strategy and dispersal capacities (Alvarez et al., 2009; Kuittinen et al., 1997; Meirmans et al., 2011; Yeaman & Jarvis, 2006). Additionally, each study included in our meta‐analysis used different genotyping and sampling strategies and covered varying proportions of the range, as well as employing different microsatellite markers, characterised by their own mutation rates. All combined, these factors have introduced noise into our analyses and potentially blurred the overall pattern, probably resulting in an underestimation of the emerging patterns. Despite the large influence of this random factor on our analysis, we were able to detect the effects of two vital processes on species ranges and their intraspecific genetic variation.

Bridging the effects of historical range shifts and contemporary demographic processes on species’ genetic constitutions over large scales has many implications to evolutionary biology. The potential for adaptation to new environments largely depends on the amount of genetic variation available, commonly sourced from standing genetic variation or arising through gene flow among populations. Central populations exhibit higher diversity and increased connectivity that can favour gene flow and adaptive introgression, conferring them an advantage towards changing conditions. Understanding this interplay between environmental conditions, genetic variation, connectivity and gene flow is crucial to address conservation needs in the face of ongoing global change. This study provides the foundation for further phylogeographical analyses to continue integrating the processes dictating genetic variation across species ranges.

CONFLICT OF INTERESTS

The authors declare no competing interests.

AUTHOR CONTRIBUTIONS

P.G.M. and J.L.‐D. designed the study and compiled the database. J.L.‐D. built the species distribution models and performed the data analysis, under supervision of P.G.M. J.L.‐D. wrote the manuscript, with input from P.G.M.

Figures S1‐S2

Click here for additional data file.

Supplementary Material S1

Click here for additional data file.

Supplementary Material S2

Click here for additional data file.

Supplementary Material S3

Click here for additional data file.