Influence of northern limit range on genetic diversity and structure in a widespread North American tree, sugar maple (Acer saccharum Marshall)

Abstract Due to climate change, the ranges of many North American tree species are expected to shift northward. Sugar maple (Acer saccharum Marshall) reaches its northern continuous distributional limit in northeastern North America at the transition between boreal mixed‐wood and temperate deciduous forests. We hypothesized that marginal fragmented northern populations from the boreal mixed wood would have a distinct pattern of genetic structure and diversity. We analyzed variation at 18 microsatellite loci from 23 populations distributed along three latitudinal transects (west, central, and east) that encompass the continuous–discontinuous species range. Each transect was divided into two zones, continuous (temperate deciduous) and discontinuous (boreal mixed wood), based on sugar maple stand abundance. Respective positive and negative relationships were found between the distance of each population to the northern limit (D_north), and allelic richness (A R) and population differentiation (F ST). These relations were tested for each transect separately; the pattern (discontinuous–continuous) remained significant only for the western transect. structure analysis revealed the presence of four clusters. The most northern populations of each transect were assigned to a distinct group. Asymmetrical gene flow occurred from the southern into the four northernmost populations. Southern populations in Québec may have originated from two different postglacial migration routes. No evidence was found to validate the hypothesis that northern populations were remnants of a larger population that had migrated further north of the species range after the retreat of the ice sheet. The northernmost sugar maple populations possibly originated from long‐distance dispersal.

This shift in range could be particularly marked in peripheral populations located along their distribution edge (Iverson, Schwartz, & Prasad, 2004). Generally isolated and smaller in size (effective population sizes) than core populations (Vucetich & Waite, 2003), these peripheral populations are expected to have a lower level of genetic diversity and a higher level of genetic differentiation among populations following genetic drift or bottlenecks (Eckert, Samis, & Lougheed, 2008;Waters, Fraser, & Hewitt, 2013). In theory, bottlenecked populations are expected to experience a reduction in allelic richness and a limited decrease in heterozygosity. The consequence of this drift would be a reduction in the number of rare alleles (Nei, Maruyama, & Chakrabordy, 1975). The central-marginal hypothesis predicts that those patterns of genetic diversity and structure are the consequences of ecological marginality in populations found at the periphery of the species range compared to core populations (Eckert et al., 2008). In contrast, "range shift following the last glacial maximum" hypothesis proposes that the patterns of genetic diversity in leading edge (i.e., colonizing front) populations are influenced by past climate-driven range dynamics rather than by demographic and evolutionary processes that are predicted in the central-marginal model (Hampe & Petit, 2005). One problem in distinguishing the relative importance of these processes is that predicted genetic patterns, such as lower genetic diversity and higher differentiation in northern edge populations, may result either from range shifts or from environmental conditions at the edge of a species range. In addition, wind-pollinated tree species may be buffered against the effect of fragmentation on the genetic structure of peripheral populations due to long-distance gene flow that contributes to a decrease in differentiation among populations (Kremer et al., 2012).
In the Northern Hemisphere, northern peripheral populations often exhibit local adaptations to cold stress (cold hardiness and earlier bud set in fall) and can potentially be maladapted to climate warming.
On the one hand, lower chilling requirements can result in early bud flush in spring under milder environments and greater susceptibility to spring frost damage (Howe et al., 2003). On the other hand, climate warming could also improve their survival and increase the success of sexual reproduction (Alberto et al., 2013).
It has major economic value as a source of saw timber and syrup production in Canada and the United States (Godman et al., 1990).
Sugar maple reaches its northern range at the transition between the temperate deciduous and boreal mixed-wood forests (Saucier, Grondin, Robitaille, & Bergeron, 2003). Palynological reconstructions support the idea of a constant presence of sugar maple since its establishment in the boreal mixed-wood forest (Liu, 1990;Richard, 1993). However, two studies have reported the presence of locally larger populations in the past along the shore of the St.
Lawrence River (Jetté & Richard, 1992;Richard, Larouche, & Lortie, 1992). Richard and Grondin (2009) reported the establishment of disjointed sugar maple stands around 8,500 calendar years BP (cal. yr BP) in northern Québec. Like many other meridional tree species in the Northern Hemisphere, sugar maple is predicted to migrate north of its current range with anticipated climate change (Goldblum & Rigg, 2005;Iverson et al., 2008).
Earlier studies of sugar maple genetic variation were conducted with provenance tests. These studies revealed clinal variations in response to temperature and moisture stresses (Kriebel, 1957).
Sugar maple provenances also showed adaptation to altitudinal gradients for respiration, photosynthesis, and specific leaf weight (Ledig & Korbobo, 1983). Analysis of neutral markers (allozymes and RAPD) revealed weak or no genetic differentiation among sugar maple populations in Canada (Perry & Knowles, 1989;Young, Warwick, & Merriam, 1993) and across its range (Gunter, Tuskan, Gunderson, & Norby, 2000). On a regional scale in Canada, Young et al. (1993) reported the presence of moderate levels of genetic structure among populations, while Diochon, Rigg, Goldblum, and Polans (2003) and Perry and Knowles (1989) found no structure among populations.
In this study, highly polymorphic microsatellite markers developed for sugar maple (Graignic, Tremblay, & Bergeron, 2013) were used to study the genetic diversity and structure of sugar maple populations. Sugar maple populations were selected along three transects that encompass the continuous-discontinuous species range at the transition between the temperate deciduous and boreal mixed-wood forests. Sexton, McIntyre, Angert, and Rice (2009) highlighted the importance of the replication of edges to disentangle equilibrium limitations (migration-drift balance) from demographic nonequilibrium characteristics (population expansion, various scenarios of colonization) of range margins. We hypothesized that marginal fragmented northern populations from the boreal mixed-wood forest would have (1) a lower level of genetic diversity and higher genetic differentiation than populations from the temperate deciduous forest; (2) lower correlations between genetic and geographic distances would be observed in marginal populations; and (3) asymmetric gene flow would occur from core to edge populations under both models. Replication of the edges makes it possible to determine whether the spatial distribution of genetic diversity in sugar maple populations reflects lasting signatures of postglacial range expansion ("range shift following the last glacial maximum"; Hampe & Petit, 2005), or if they are characteristic of peripheral populations per se (central-marginal hypothesis; Eckert et al., 2008). Graignic, Tremblay, and Bergeron (2016) previously reported the presence of higher genetic diversity and a lower inbreeding coefficient in mature trees than in younger regeneration cohorts in three natural sugar maple populations. Mature tree and sapling cohorts were sampled to explore whether the distribution of genetic diversity was different between cohorts in a larger area.

| Study area and sampling
The study area was located at the northern range limit of sugar maple in Québec (eastern Canada; Figure 1). The presence of sugar maple stands along the latitudinal gradient was estimated from an analysis of large inventory databases (for more details, see Graignic, Tremblay, & Bergeron, 2014). The latitudinal gradient was divided into two zones, based upon the proportion of sugar maple stands in the continuous (temperate deciduous forest) and discontinuous (boreal mixed-wood forest) zones (referred as C and D zones below, re- Sites were distributed along three south-north transects (1, western; 2, center; 3, eastern), with eight sites per transect (four sites per transect per zone; n = 24; we analyzed n = 23 as population 2-D-B needs to be thrown out, see the Markers genetic diversity section below for more details; Figure 1). All the sites were old-growth, uneven-aged stands, selected to be as similar as possible (for more details, see Graignic et al., 2014). Eleven of the sites selected were either old-growth or rare forests that are classified as Exceptional  (Foré, Hickey, Guttman, & Vankat, 1992;Mulcahy, 1975),  Graignic et al. (2014) for more details; population 2-D-B was thrown out due to amplification errors), and proportional membership of each sugar maple sample in clusters for four genetic groups (K = 4) inferred by structure analysis and analyzed. Samples were frozen (−80°C) or dried (using silica gel) until needed for genetic analyses.

| Molecular methods
DNA was extracted using Extract-N-Amp™ Plant PCR Kits (Sigma-Aldrich, Oakville, ON, Canada). All samples were genotyped for 18 variable microsatellite loci using PCR and genotyping protocols as previously described by Graignic et al. (2013). Five different multiplex PCR sets were used, and 37 cycles in the PCR amplification were processed (see Table S1).

| Markers genetic diversity
The 2-D-B population had to be removed from the analysis due the low number of successful amplification, which was likely caused by poor sample conservation. Accordingly, a total of 913 individuals were selected for the subsequent analysis. For each locus, the total number of alleles (A T ), mean number of alleles (A), mean observed (H O ) and expected (H E ) heterozygosity, and inbreeding coefficient (F IS ) were estimated using fstat 2.9.3.2 (Goudet, 2001). Departure from Hardy-Weinberg equilibrium (HWE) per locus in each population and linkage equilibrium between all pairs of loci in each population were calculated using an exact test implemented in genepop 4.2.1 (Rousset, 2008).
Markov chain parameters for HWE were 10,000 dememorizations, followed by 500 batches of 5,000 iterations per batch. We corrected for multiple comparisons using a sequential Bonferroni adjustment of p-values to a predetermined experiment-wise error rate of .05 (Rice, 1989). Null allele frequencies were estimated using freena (10,000 replicates; Chapuis & Estoup, 2007). This software was chosen because it uses the Dempster, Laird, and Rubin (1977) algorithm, which provided the most accurate estimate of several algorithms tested in Chapuis and Estoup (2007). The presence of null alleles overestimated fixation index (F ST ) and genetic distances (Chapuis & Estoup, 2007

| Genetic diversity and differentiation between populations and groups
The mean number of alleles per locus (A), mean allelic richness (A R ), H O , H E , pairwise F ST , mean pairwise F ST , and F IS were estimated using fstat. The A R was calculated using rarefaction method, based upon the size of the smallest number of samples. The indices were estimated for each population and on mature trees and saplings for each population separately. To assess whether population diversity parameters differed between zones and transects, A R , H O , H E , F ST , and F IS were calculated by population (all individuals pooled, mature trees, and saplings) within each zone (discontinuous and continuous), transects (1, 2, and 3), and zone in each transect separately using fstat (tested for significance using 1,000 permutations).  Graignic et al., 2014). The population size of sugar maple was calculated separately for mature tree and sapling using their densities within the quadrat and stand surface areas. For stand surface areas, we used data from EFE and ecoforestry maps that were obtained from MFFPQ (MFFPQ 2016a,b). Both central-marginal and range shift models predicted that the closer the populations are to the northern limit, the more they are differentiated and the less they are diversified. We fixed the northern limit between balsam fir-white birch (A. balsamea-B. papyrifera) and balsam fir-yellow birch (A. balsamea-B. alleghaniensis) bioclimatic domains, and the distance of each sampled population to this limit was estimated with arcmap™ 10.0 (Environmental Systems Research Institute, ESRI, CA, USA); negative values were assigned to sites located to the north of the northern limit, that is, in the balsam fir-white birch (A. balsamea-B. papyrifera) bioclimatic domain. Correlations between variables were examined, and assumptions of normality and homoscedasticity were graphically tested on global models. Where needed, the data were square-root-transformed to meet assumptions. Fits for the 19 models were compared using the Akaike information criterion (AIC; aictab in the AIC c modavg library; Mazerolle, 2011), corrected for small sample sizes (AIC C ). We computed AIC C and ΔAIC C weights (ω) to determine the strength of evidence for each model (Burnham & Anderson, 2004). We then performed multimodal inference (modavg function in AIC c modavg library), where required (ΔAIC C ≤ 4), to calculate the parameter estimates and unconditional 95% confidence intervals (CI).

| Bottleneck tests
Evidence of bottlenecks in each population was evaluated using three tests: heterozygote excess, allele frequency mode shift, and M-ratio.
For the first test, heterozygote excess, we used bottleneck 1.2.02 (Cornuet & Luikart, 1996) and each population was tested under three mutation models, namely IAM (infinite allele model), SMM (stepwisemutation model), and TPM (two-phase mutational model; Piry, Luikart, & Cornuet, 1999). The significance of the test under all models was performed using one-tailed Wilcoxon signed-rank test with 1,000 permutations. The allele frequency mode shift was also conducted in each population using bottleneck (Luikart & Cornuet, 1998).
We also calculated the M-ratio, that is, the mean ratio of the number of alleles (k) to the range in allele size (r), using m_p_val software (Garza & Williamson, 2001). This method assumes that rare alleles are lost faster during a bottleneck, reducing the number of observed allelic states (k) faster than the size range of those alleles (r), thereby resulting in a reduced M-ratio (M = k/r). To determine the significance of a ratio, we compared our observed value to critical values (M c ) that were established using the critical_m program (Garza & Williamson, 2001).
When the observed M-value was significantly (p < .05) below M c , we assumed that a population had experienced a significant reduction in size (Garza & Williamson, 2001). To calculate M c and p-values, we used three different ancestral theta values (θ = 4N e μ = 1, 5 and 10; N e , the effective population size; the mutation rate, μ) to account for a wide range of effective population sizes; and the parameters were set as recommended by Garza and Williamson (2001). Each set of simulations consisted of 10,000 iterations.

| Isolation by distance
Isolation by distance was investigated using a Mantel test (1,000 permutations; Pearson's correlation; using the mantel function in and geographic distance (km). Pairwise F ST and Nei's genetic distance, D S (Nei, 1972), were calculated, respectively, using fstat and populations 1.2.32 (Langella, 1999). Geographic distances were estimated using the earth.dist function in the fossil library in R (Vavrek, 2012). We performed isolation-by-distance analysis for all populations, populations for each zone, and each transect separately.

| Genetic structure
Neighbor-joining (NJ) trees were constructed based on D S and using 1,000 bootstraps over loci in populations, one tree for all populations and one tree using cohort separate populations in two parts. NJ trees are useful for describing genetic relationships among populations. The tree was visualized using treeview (Page, 1996).
T A B L E 1 Genetic variability estimates of sugar maple (Acer saccharum) populations in Québec for all individuals We used the Bayesian clustering method, implemented in structure 2.3.4 (Pritchard, Stephens, & Donnelly, 2000), to assess the number of populations (K), assign individuals to each K cluster, detect the population structure, and show whether geographic populations are assigned in genetic populations. Our analyses were based upon an admixture ancestral model with correlated allele frequencies, and a priori sampling locations were used as prior information (LOCPRIOR setting). LOCPRIOR was used to detect any further structure unidentified by standard settings (Hubisz, Falush, Stephens, & Pritchard, 2009). Twenty independent runs were performed for each value of K (1-26) with a burn-in of 100,000 followed by 200,000 MCMC iterations. The most likely value of K was determined using the ∆K criterion (Evanno, Regnaut, & Goudet, 2005). structure Harvester version 0.6.93 was used to extract the results and created a graphic of the ∆K criterion (Earl & Vonholdt, 2012). The results were visualized for the best K, with distruct version 1.1. (Rosenberg, 2004).
Genetic variation was examined using the hierarchical analysis of molecular variance (AMOVA) based on Φ PT statistic and conducted in genalex 6.5b3 (Peakall & Smouse, 2006). We examined the significant difference (9999 permutations) between transects (1, 2, and 3), between populations within transect and within populations, and between zones (C and D), between populations within zone and within populations for all populations, and for each transect separately. The results of K-means clustering were used to cluster populations into four groups:

| Gene flow and population sizes
We used the Bayesian approach implemented in migrate-n version 3.6 (Beerli, 2006)   as an exponential window prior set for both parameters (min = 0, mean = 5, max = 50 and ∆ = 5 for θ; min = 0, mean = 5, max = 50 and ∆ = 5 for M). We replicated two long chains of 10,000,000 genealogies, which were recorded every 10 steps, after a burn-in period of 10,000. The static heating scheme was set to four chains with temperatures of 100,000.0, 3.0, 1.5, and 1.0, and the swapping interval set at one. The program was run several times adding more sampling schemes and replications. The starting parameters and the resulting estimates were compared until the results were congruent.  Term abbreviations: m_ers_BA, mature sugar maple basal area (m 2 /ha); s_ers_BA, sugar maple sapling basal area (m 2 /ha); m_ers_d, mature sugar maple density (stems/ha); s_ers_d, sugar maple sapling density (stems/ha); m_ers_PS, mature sugar maple population size (stems); s_ers_PS, sugar maple sapling population size (stems); D_north, distance of each site to the northern limit (km); SE, standard error.

| Nuclear microsatellite diversity statistics
The total and mean number of alleles per locus ranged between 4 (SM60) and 31 (SM21A) and from 2.8 (SM60) to 14.8 (SM21A), respectively (Table S1). Significant departure from HWE was observed in 37 (9%) of the 414 possible population-locus combinations (Table S2). Four of the 18 loci (SM22, SM27, SM47, and SM56) failed to meet HWE in ≥6 (26%) populations. For those loci, deviations from HWE were due to heterozygote deficiencies (F IS ≥ 0.227; Table S1). Significant linkage disequilibrium between pairs of loci (within populations) was detected in only one of 3,519 tests after Bonferroni correction (SM14 × SM34 in population 3-C-A). The loci were therefore all considered genetically independent.
Markers that showed a departure from HWE also had high frequencies of null alleles (≥0.10) in most populations (Table S3). Pairwise F ST estimates obtained before and after corrections for null alleles were similar (r = .98, p < .001). Consequently, all loci (18) were retained in the analysis.

| Genetic diversity and differentiation
Population 1-D-C (Rémigny) from the discontinuous zone of transect

| Bottlenecks
Inferences to be drawn from the heterozygosity excess tests  (Table S9). No significant bottlenecks were found for the other θ values (Table S9).

| Genetic structure
The neighbor-joining tree separated the populations of the transect 1 (except for 1-D-C, Rémigny) from the other two transects ( Figure S1).

| Gene flow and population sizes
The largest mean scaled effective population size (θ) was in the eastern area of Québec (transect 3; mean θ = 3.68), and the smallest was in the center (transect 2; mean θ = 1.22; Table 5, Figure S3).
Northern populations also had a low mean scaled effective population size (θ = 1.43). Remarkable levels of gene flow were observed among populations from the three transects, with a slightly more scaled immigration rate (M) from the east (transect 3) than from the other two transects (Table 5, Figure S3). Gene flow was asymmetrical, with higher migration rates from the south (transects 2 and 3) toward the northern populations, while low and equal migration rates were detected between the west (transect 1) and northern populations (Table 5, Figure S3).

ST1→i
2.77 (1.7-3.80) 6.6 (4.6-7.5) -15.5 (14.6-16.4) 15.6 (14.5-16.9) level. How fragmentation affects evolutionary processes is unlikely to be universal and to depend upon several factors, such as the effective population size, the extent of gene flow, and the level of connectivity between populations at the range edge. This is the first time that the genetic diversity and structure of sugar maple populations have been studied over a wide range at its northern range. Our results showed evidence of the impact of marginality on sugar maple genetic structure. Lower allelic richness (A R ) and higher genetic differentiation (F ST ) between populations were detected in the discontinuous zone where populations are at the periphery of the sugar maple range in Québec. Furthermore, gene flow between zones was asymmetric, mostly directed outward from the populations of the temperate deciduous forest. A high level of admixture and a weak population structure were also detected. Most genetic variation was found within sugar maple populations, a frequent pattern for long-lived perennials, outcrossing, and wind-dispersed seeds plants (Nybom, 2004). It may also introduce alleles preadapted to a warmer climate (nonlocal genes; Kremer et al., 2012), which may facilitate adaptation in the boreal mixed-wood sugar maple populations in the context of climate change.
A genetic difference between our western transect and others transects led us to think that migration influence this structure, that is the range shift model rather than the ecological marginality in itself was the best explaining model for sugar maple.

| Genetic diversity levels in sugar maple populations
The eighteen microsatellite markers developed by Graignic et al.
(2013) revealed a high level of genetic diversity in sugar maple populations (see Table 1). Only two studies using microsatellite markers on natural population of sugar maple (Graignic et al., 2016;Khodwekar, Staton, Coggeshall, Carlson, & Gailing, 2015) were available for comparison with our results. Our results were in the range reported by Khodwekar et al. (2015) (H E = 0.822), by Graignic et al.

| Genetic variation and latitudinal gradient
We found a lower allelic richness (A R ) and higher genetic differentiation (F ST ) between populations in the discontinuous zone that represents the periphery of sugar maple range in Québec ( Table 2).
The relationship between the proximity to the northern limit and genetic diversity was negative, whereas it was positive with genetic differentiation (Table S8,  3-D-A) were assigned to a distinct group with structure (Figures 1 and   S3). These results demonstrated the effect of a latitudinal gradient on sugar maple genetic diversity.
However, when these relationships were tested for each transect separately, the pattern (discontinuous/continuous) remained significant only for the western transect (transect 1). The reduction in A R and the increase in F st along the western transect are most likely the consequence of successive founder events that followed the retreat of the ice sheet and are more consistent with the range shift hypothesis (Hampe & Petit, 2005). These relationships were not due to isolation by distance because no significant correlation between genetic and geographic distances was detected within transect (  (Graignic et al., 2014). The association between lower sexual recruitment and lower genetic diversity in peripheral populations was also suggested for Thuja occidentalis L. (Paul, Bergeron, & Tremblay, 2014;Xu, Tremblay, Bergeron, Paul, & Chen, 2012) and for the dwarf thistle Cirsium acaule (L.) A.A. Weber ex Wigg (Jump, Woodward, & Burke, 2003) at their northern limits. The northwestern region is included within a larger physiographic region that was created by lacustrine deposits from the proglacial lakes Barlow and Ojibway (Veillette, 1994). Clay and organic deposits increase from south to north, while rock and till deposits, which constitute the most favorable habitats for maple establishment, tend to decrease. These conditions contribute to further limiting the distribution of the species (Tremblay, Bergeron, Lalonde, & Mauffette, 2002). Lower sexual regeneration, combined with the reduction in the number of suitable sites for recolonization following disturbances, may have limited the long-term maintenance of the species in local stands.
In the eastern and central regions, A r and F st did not change from core to edge populations, indicating no direct effect of fragmentation on genetic diversity of sugar maple populations. Paleoecological records support the constant presence of sugar maple in mixed-wood or boreal forests since its establishment (Liu, 1990;Richard, 1993).
We previously reported constant sugar maple recruitment over time in these northern populations (Graignic et al., 2014). Therefore, it is plausible that the number, size, and connectivity of suitable habitat patches were sufficient to maintain stable levels of gene flow and genetic diversity toward the range edges in populations from the central and eastern transects (Sexton et al., 2009).
No recent bottleneck was detected under the SMM, TPM, and mode shift indicator tests, while most populations were in bottleneck under the IAM (Table S9). Under IAM, microsatellite loci have been shown to exhibit heterozygosity excess in stable populations (Luikart & Cornuet,1998). This result suggests there have been no significant reductions in effective population size. The M-ratio analysis proved to be rather sensitive to the choice of θ. The results showed that only three populations may have undergone an historical bottleneck when the mutation-scaled effective population size was estimated at θ = 1, but no bottleneck was detected at θ = 5 or higher (Table S9). Among these populations, only individuals from population 1-D-C (Rémigny) were assigned to a distinct cluster (defined by structure) and had probably undergone a bottleneck in a distant past (Table 5 and Figure S3; sampled when we compare to the whole sugar maple range, we were able to determine that the range shift model rather than the ecological marginality in itself was the best explaining geographical variation in the genetic structure of sugar maple.

| Special case of 1-D-C, Rémigny
The population from Rémigny (1-D-C) had the lowest level of diversity (A, A R , H E , and F IS ), and the highest F ST (Table 1), and a bottleneck was detected (Table S9; M-ratio at θ = 1). Both neighbor-joining trees and structure analysis clearly disjointed this population from all the others ( Figure S1, 1). These differences were probably due to the recent establishment of sugar maple at this site. A study conducted by Pilon and Payette (2015) at a nearby sugar maple stand (about 7 km away) dated sugar maple's arrival after a fire, which occurred 160-220 yrs ago. According to Pilon and Payette (2015), the presence of sugar maple at Rémigny was sporadic or scarce during the last millennia.
Probably, its presence was too recent and the time lapses too short to reach a level of genetic diversity similar to the other sugar populations. In their work on Pinus ponderosa, Lesser, Parchman, and Jackson (2013) reported that more than 230 yrs following stand establishment was required to reach allele saturation.

| Mature tree and sapling cohorts
In many coniferous tree species, the level of heterozygosity increases with the age of the trees (seeds vs. adults, seeds vs. seedlings, different age classes), possibly due to higher fitness of heterozygous individuals (Bush & Smouse, 1992;Nijensohn, Schaberg, Hawley, & Dehayes, 2005). Ballal (1994) reported that sugar maple embryos had lower heterozygosity and higher F IS than seedlings (1-yr to ≤1 cm basal area) and mature trees (≥30 cm dbh). Graignic et al. (2016) reported the same difference between younger regeneration cohorts (seedlings and saplings) and mature trees. Foré et al. (1992) found no difference between sugar maple embryos, 1-yr seedlings, and three other cohorts (dbh ≤2 cm, 15-25 cm and ≥40 cm). We found a significant positive relationship between mature sugar maple basal areas and H O (Table S8, 3 and Figure 2), but no relationship with F IS . Overall, these results did not support the hypothesis of higher heterozygosity in older cohort.

| Genetic signature of migration routes
We observed a genetic difference between the eastern and the western transects (blue and yellow clusters, see Figure 1). The largest population size was found in the eastern populations followed by the western and the central populations, respectively (Table 5; Figure S3).
This result may possibly represent two distinct migration routes (first from the east, then from the west) that separated a long time ago, or originated from different glacial refugia.
Earlier sugar maple migrations in Québec started about 9,900 cal. yr BP in southeastern Québec (step 1 in Figure S4; Lavoie & Richard, 2000; see all steps in Figure S4), and sugar maple pollen was recorded 6,300 cal.
yr BP in the southwest of Québec (step 5; Vincent, 1973). Therefore, the arrival of sugar maple in western Québec occurred around 3,000 yr later than it did in the east, which corroborated evidence for two distinct migration routes in Québec. Our hypothesis also matches Braun's hypothesis, which was based on pollen data analysis that suggests beech-maple associations moved northward via two routes following ice retreat (for more details, see Figure S4; Braun, 1950). Jackson et al. (2000) identified, from pollen and macrofossils data, two maple glacial refugia during the LGM (21,000 cal. yr BP; Figure S4). Interestingly, a genetic signature of these same refugia was identified for red maple (Acer rubrum L.; McLachlan et al., 2005) and multiple refugia for sugar maple (Vargas-Rodriguez, Platt, Urbatsch, & Foltz, 2015).
Southern populations in central Québec (transect 2) had a lower population size (Table 5; Figure S3) and were clustered with western and eastern populations (transects 1 and 3; Figure 1). We hypothesized that these populations originated from an admixture of westward and eastward Québec migrations (step 7 in Figure S4), with greater contributions from the eastward advance ( Figure S3). Similar routes for northernmost populations could also be drawn, but with considerably more contributions from south-central-eastern populations (transects 2 and 3; step 7 in Figure S4) because (1) very low contributions from the southwestern populations to the northernmost populations were found (Table 5, Figure S3), and (2) western populations arrived later than did the eastern populations in Québec (Lavoie & Richard, 2000;Vincent, 1973).

| CONCLUSIONS
The inclusion of among-edge population comparisons provided in this study provides new insights into the question of whether the "ecology marginality at the periphery" (central-marginal hypothesis; Eckert et al., 2008) or "range shift following the last glacial maximum' (Hampe & Petit, 2005) models best describes the pattern of genetic variation in sugar maple populations. The pattern of variation in genetic diversity along the western transect is best supported by the range shift model. No evidence was found to validate the hypothesis that northern populations were remnants of a larger population that had migrated further north of the species range after retreat of the ice sheet. It is possible that northernmost populations originated from long dispersal distances (Cain, Milligan, & Strand, 2000). Additional paleoecological data would be essential to accurately establish the arrival of the northernmost populations and their migration routes. Moreover, this study was performed on a rather small portion of the sugar maple range. Range-wide sampling might provide different insights.

ACKNOWLEDGMENTS
We thank Normand Villeneuve, Claude Poulin, and Renaud Dostie of the MFFPQ who found many potential sites. Geneviève Trudeau,

CONFLICT OF INTEREST
None declared.

AUTHOR CONTRIBUTIONS
N.G., F.T., and Y.B. conceived the ideas. N.G. collected samples, performed the laboratory work, analyzed the data, and led the writing of this manuscript. F.T. performed some statistical analysis and revised the manuscript.

DATA ACCESSIBILITY
GenBank access numbers of microsatellites markers are listed in Table   S1. Microsatellite data are in DRYAD entry doi: 10.5061/dryad.36634.