Introduction

Disturbances generate spatially heterogeneous landscape structures that exert an influence on several ecological processes (Lindenmayer et al. 2008). Forest fires, for instance, increase landscape heterogeneity by creating mosaics of unburned and burned patches of different sizes and that feature a wide range of burn severities (Turner et al. 2003). Residual patches constitute habitat islands that can serve as refuges for fire-sensitive species (Gasaway and DuBois 1985; Gandhi et al. 2001) and as sources for the recolonization of burned areas (Keeton 2000; Asselin et al. 2001).

Few studies have examined the relative importance of fire disturbances and residual unburned forest patches on the genetic diversity of natural populations (Banks et al. 2013; Lind-Riehl and Gailing 2015). The post-fire vegetation-free landscape can also facilitate pollen movement in wind-pollinated tree species (Bacles et al. 2005; Shohami and Nathan 2014). A high level of gene flow between remnant populations contributes to homogenize genetic structures in naturally regenerating stands of Nothofagus (Premoli and Kitzberger 2005). Conversely, stand-replacing fire disturbances may decrease the presence of remaining reproductive trees for an extended period of time leading to genetic erosion in the progeny (Premoli and Steinke 2008). This can be particularly true in terms of allelic richness (Widmer and Lexer 2001), and an increase in selfing in post-fire tree cohorts (Shohami and Nathan 2014).

Ecosystem-based forest management in North America now largely uses natural disturbance dynamics to guide decision making (Gauthier et al. 2009). Therefore, it is important to understand the relationship between landscape mosaics and genetic diversity to determine how much habitat is required to meet conservation objectives. Large forest cover is generally considered crucial for conservation due to the relation between the population size of individual species and species diversity (Lindenmayer et al. 2008). In addition, connectivity between remnant stands is most likely to occur when the amount of particular types of land cover are high (Hannon and Schmiegelow 2002). While large areas are important, small and medium-sized remnant stands often have considerable ecological value (Turner et al. 1997). Small fire refuges found in the boreal mixed-woods are characterized by long-term persistence (e.g. ecological continuity) (Ouarmim et al. 2014). These stands can escape consecutive fires for several centuries, up to a few millennia, which is a much longer period than the fire frequency recorded in the surrounding landscape. Their persistence in the landscape can protect populations against climatic variation and disturbance and promote population persistence and demographic stability, as well as increased genetic diversity and stability (i.e. reduced frequency in temporal change of alleles) (Davies et al. 2016). The level of biodiversity of fire refuges in the boreal mixed-wood has yet to be documented and such stands have been largely ignored in conservation policies in North American boreal forests (Gauthier et al. 2009).

In the eastern Canadian boreal landscape, late successional forests and small forest remnants (fire skips) are often dominated by eastern white cedar (EWC, Thuya occidentalis L.), a long-lived, fire-sensitive conifer tree species, that becomes dominant 200–250 years after disturbance (Bergeron 2000). Remnant tree stands represent the sole potential seed sources available for regeneration after forest fires for this species (Asselin et al. 2001). The natural forest remnants of EWC function as refuges in a landscape that is affected by recurrent large-scale wildfires (Ouarmim et al. 2014).

Forest remnants have been identified as historical or potential reservoirs of genetic diversity in several tree species in fire-prone landscapes (Uchiyama et al. 2006.; Mosseler et al. 2003; Premoli and Kitzberger 2005). We took advantage of the co-existence of three types of landscapes characterized by different fire regime and geographical context to understand the effects of their features on the genetic structure of EWC. The analysis included; (i) non-fragmented forests originating from large stand replacing fires on the mainland, (ii) naturally fragmented EWC sites located in an open landscape with fires of mixed severities on islands, and (iii) fire refuges that have remained fire-free for very long periods during the Holocene and constitute small fragments nested in a surrounding landscape matrix of early successional forest. We address the specific question of how does forest spatio-temporal continuity affect the genetic diversity of EWC?

Materials and methods

Study area

The study area is located in northwestern Quebec, in the Lake Duparquet Research and Teaching Forest (79º10′W–48º30′N), which covers an approximately 80 km2 boreal landscape and more than 170 islands. The forest is dominated by balsam fir (Abies balsamea [L.] Mill.), EWC, white spruce (Picea glauca [Moench] Voss), and black spruce (Picea mariana [Mill.]BSP). Pioneer species, such as paper birch (Betula papyrifera Marsh), trembling aspen (Populus tremuloides Michx.), and jack pine (Pinus banksiana Lamb.), occupy large areas following a disturbance (Bergeron and Dubuc 1989).The periodic occurrence of fires since regional deglaciation (cal. 8000 BP) has been well documented (Carcaillet et al. 2001). Eight major fires burned in this forest between 1717 and 1964 (1964, 1944, 1916, 1870, 1847, 1823, 1797, 1760) (Bergeron 1991, 2000) whereas several fires of mixed severity occurred on the islands (Bergeron et al. 1991; Drobyshev et al. 2010). The study forest was selected because it has been minimally affected by human intervention (Bergeron 2000) .

The species

Eastern white cedar (EWC; Thuya occidentalis L.), also known as northern white cedar, has a large natural range in North America. It is a monoecious and wind-pollinated conifer species. Its flowers usually develop in late April or May; cones mature in August but may remain attached to the tree until the following spring. Seeds are mainly disseminated by wind from September to November. Eastern white cedar is ill adapted to fire and surviving seed trees are needed to reinvade burned areas (Asselin et al. 2001; Bergeron et al. 2004). Relatively small groups of mature trees disperse their seeds by wind to colonise a nearby area affected by disturbance and this process is strongly influenced by the distance from seed sources (Asselin et al. 2001). Subsequent expansion of EWC occurs some decades after the break up of the post-fire cohort. The pattern of recruitment is attributed to prolific seeding and layering on a seedbed devoid of deciduous litter (Bergeron 2000).

Sampling

Nine sites (Fig. 1; Table 1) covering three different types of landscape characterized by different fire regimes were selected (Bergeron 1991; Ouarmim et al. 2014). Three small fire refuges (< 3 ha) located in an area that last burned in 1944 were selected to represent sites fragmented in a forest landscape dominated by early successional broadleaf tree species. The long-term fire history has been reconstructed for the three sites [f441 (Georges), f442 (Venteux), and f443 (Cadeau)] (Ouarmim et al. 2014). The three sites have remained fire-free for very long periods during the Holocene (520–4095 years) and only burned during what appear to have been exceptionally severe fires. The site f443 remained fire-free for at least 4095 years during the mid-Holocene (6275–2180 cal. BP). Long fire-free periods were recorded at site f442 (520 years, between 715 and 195 cal. BP) and f441 (at least 1470 years) during the late Holocene. The resulting total number of fire events detected was 9 at f441, 3 at f442 and 18 at f443. The mean fire-free interval was estimated to be 455 years at f441, 222 years at f442 and 291 years at f443. The number of EWC adult trees ranged from 20 (f442), 25 (f441) to > 100 (f443) (Table 1). Three islands in the lake (> 100 EWC adult trees)—is39, is42, and is134—last burned in1889, 1825, and after 1825 (precise year unknown), respectively, representing the lake-island landscape, were selected. Three old sites (f60, f97, f23), which had 500 to more than 1000 EWC adult trees per hectare and were last burned in 1760, 1797, and 1823, were selected to represent non-fragmented sites on the mainland in a forest dominated by late successional coniferous tree species (Bergeron 2000). The fire history on the mainland and islands was reconstructed by using fire scars and dendrochronological analysis (Bergeron 2000). The mean fire-free intervals on the islands and the mainland are similar and are estimated to be 74 and 63 years respectively.

Fig. 1
figure 1

Map of the study area showing the locations of the nine Eastern white cedar sites from three landscape types in Quebec, Canada. The four clusters (K = 4) from the STRUCTURE analysis (Pritchard et al. 2000) are displayed in blue, yellow, orange, and green

Full size image
Table 1 Location and genetic diversity parameters in nine Eastern white cedar (Thuya occidentalis L.) sites from three landscape types in Quebec, Canada
Full size table

Only EWC trees with a diameter at breast height (d.b.h.) of more than 25 cm were sampled for foliar or cambial root tissues. The distance between selected trees was greater than 20 m whenever possible to reduce the probability of selecting closely related individuals. For each site, 20–30 trees (Table 1) were selected and used for microsatellite genotyping. All the sampled materials were kept in a freezer (− 20 °C) before the genetic analysis.

DNA extraction and genotyping

Foliar or cambial root tissues were ground in liquid nitrogen, and genomic DNA was extracted and amplified over 18 microsatellite loci (Appendix 1) using previously described methods (Xu et al. 2013). Briefly, each reaction mixture contained 1 µl of DNA extract, 5 µl of 2X QIAGEN Multiplex PCR Master Mix (QIAGEN), and a final concentration of 0.2 µM for each forward and reverse primer. The PCR program consisted of an initial heat-activation step at 95 °C for 15 min, 36 cycles of 3-step cycling: denaturation at 94 °C for 30 s, annealing at 54 °C for 90 s, extension at 72 °C for 60 s, and a final extension at 60 °C for 30 min. Amplified DNA was sent to Genome Quebec (Montreal, Canada) for fragment analysis on an ABI 3730 genetic analyzer (Applied Biosystems, Carlsbad, California, USA). The results were analyzed using GeneMapper 3.7 software (Applied Biosystems).

Microsatellite Marker analysis

The presence of null alleles and scoring errors were tested in Micro-Checker (Van Oosterhout et al. 2004). The basic genetic diversity estimates within stands, including the number of alleles per locus (Na), observed heterozygosity (Ho), expected heterozygosity (He), and private alleles (PA), were calculated in GenAlex 6.4 (Peakall and Smouse 2006). Allelic richness (AR) was estimated in FSTAT 2.9.3.2 (Goudet 2001). We also tested for linkage disequilibrium (LD) over 18 loci (153 pairs) for all pooled samples. The Bonferroni correction was applied to get the adjusted p-values of the significance test. GENEPOP (version 4.2) (Rousset 2008) was used to test deviations from the Hardy–Weinberg equilibrium (HW) and to estimate P-values using the Markov Chain method with 1000 iterations.

Population genetic structure and differentiation

The samples were grouped for each landscape type, and genetic diversity estimates including AR, Ho, He, inbreeding coefficient (Fis), population differentiation (Fst) and population relatedness (Rel) among the groups (1000 permutations) were compared. The AR was computed for each locus and population using the rarefaction method (based on a minimum sample size of 20 diploid individuals). The pair-wise genetic differentiation among the stands and the associated significance were tested with a G-test and 36,000 permutations. All the analyses were performed in FSTAT 2.9.3.2 (Goudet 2001).

The model-based Bayesian clustering approach implemented in STRUCTURE 2.3.4 (Pritchard et al. 2000) was used to detect any regional gene pool subdivision predating stand formation. It was used to assign individuals to genetic clusters (K) and to estimate admixture proportions (Q) for each individual. Our analyses were based on 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 structures unidentified by standard settings. Ten independent runs were performed for each value of K (1–9) 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 et al. 2005). STRUCTURE HARVESTER version 0.6.93 was used to extract the results and created a graphic of the ∆K criterion (Earl and Vonholdt 2012). The results were visualized for the best K, with DISTRUCTversion 1.1. (Rosenberg 2004).

A principal coordinates analysis (PCoA) was also performed to visualize potential groupings of populations according to landscape features using the covariance-standardized method based on the population genotypic genetic distance in GenAlex 6.4 (Peakall and Smouse 2006).

Bottleneck analysis

We ran a bottleneck analysis to determine whether EWC populations had undergone significant reductions in Ne (i.e., genetic bottlenecks due to landscape features in islands, refuges, mainland). Evidence of bottlenecks at each site was evaluated using two tests: heterozygote excess and allele frequency mode shift. For the first test, heterozygote excess, we used BOTTELENECK 1.2.02 and each site was tested under two mutation models, single-step mutation model (SMM) and two-phase mutation model (TPM). Under infinite allele model (IAM), microsatellite loci have been shown to exhibit heterozygosity excess in stable populations (Luikart and Cornuet 1998). For the TPM, we set the multistep mutation events to 5% and variance to 30. The significance of the test under all models was performed using one- and two-tailed Wilcoxon sign-rank tests with 1000 permutations. The second approach, the allele frequency mode shift, was also conducted in each population using BOTTLENECK. Under recent bottlenecks, rare alleles may occur at lower frequencies leading to a mode-shift distortion of the typical L-shaped allele frequency distribution (Luikart and Cornuet 1998).

Gene flow

The correlation between the geographic and genetic distances (Fst) (Weir and Clark Cockerham 1984) among the populations was tested with a Mantel test (999 permutations) in GenAlex 6.4 (Peakall and Smouse 2006). A Bayesian multilocus genotyping procedure implemented in BAYESASS v. 1.3 using MCMC methods (Wilson and Rannala 2003) was used to estimate the direction and rate of recent immigration. This approach was used to determine the relative contributions of gene flow versus drift on EWC population genetic structure and the effect of landscape features on contemporary gene flow. We also examined the pattern of directionality in migration patterns. The number of generations was estimated to be lower than five generations (or < 150 years) given the estimated generation time of 20–30 years. The number of iterations for the chain was 3,000,000, with 999,999 being burn-in, and the sampling frequency was 2000. An increase in the iterations did not affect the output. We repeated the analysis three times to verify its stability.

Environmental association analyses

Our aim was to establish how the mean fire intervals and the number of fire events were associated with the genetic structure among stands. A Multiple Regression Model (MRM) approach implemented in the program GESTE ver. 2.0 (Foll and Gaggiotti 2006) was used to evaluate the effect of these factors on the genetic structure (Fst) of EWC populations. GESTE is a Bayesian inference method that relates Fst values for each local population to environmental factors using a generalised linear model. Data from microsatellites were run against the two uncorrelated stand-specific environmental variables. Ten pilot runs of length = 5000 were used to adjust the acceptance rates for each parameter of the MCMC chain to the recommended range of 0.25–0.45. Thereafter, 50,000 iterations were performed per run, and a thinning interval of 20 was used.

Results

Microsatellite Marker analysis

All microsatellite loci examined were polymorphic, having a mean number of alleles per locus of 9.7. (Appendix 1). Their exclusion power (PI) was very high, and the probabilities of a match between the profiles of two unrelated individuals in a randomly mating EWC population ranged between 1.0 × 10−14 and 7.0 × 10−15. We did not detect any null alleles or scoring errors. The total number of alleles per locus ranged from 3 (TO605, TO418) to 19 (TP11) and the effective number of alleles ranged from 1.731 (TO512) to 7.197 (TO791) (Appendix 1). Among the 18 loci, 11 had a significant excess of heterozygotes (P < 0.00278, adjusted 5% nominal level: 0.00278), and three had a significant deficit of heterozygotes (P < 0.00278) (data not shown). Among 153 loci pairs, 145 were not in LD (P > 0.0007, adjusted 5%: 0.000327), and only eight were weakly linked (P = 0.000330) and could be considered to be in LD (data not shown).

Genetic diversity parameters within sites are listed in Table 1. The Na ranged from 4.9 (f442) to 6.4 (f23), with an average value of 5.8. The AR averaged 5.4, ranging from 4.9 (f442) to 5.8 (f23). The mean Ho was high (0.718), ranging from 0.674 (f443) to 0.774 (is39), and the mean He was 0.619 (Table 1). The inbreeding coefficients (FIS) were generally negative, and a significant deviation from the HW equilibrium with an excess of heterozygotes was detected in the global test (FIS = − 0.171; PHW = 0.0000). All sites contained private alleles, with the mainland having the highest number of private alleles (15), followed by the islands (8) and small fire refuges (5).

Population genetic structure and differentiation

There were significant differences in AR, Fst and Rel among the three different landscape groups (Table 2). Small fire refuges had the lowest AR (5.061; P = 0.02) and the highest Fst (0.052; P = 0.05). Average pair-wise relatedness estimates (r) for each site were always less than 0.125 (third-order relationship) and varied from 0.007 (is134) to 0.056 (f442). No significant difference was observed for the other indices (Ho, He, Fis). Most of the pair-wise Fst values (34 out of 36) were in the lower range (Fst < 0.06, which is typical of long-lived outbreeding, wind-pollinated tree species) but were significantly different from zero, indicating the presence of genetic differentiation (Table 3). Within each group, the pair-wise Fst values were the highest between small fire refuges (0.0473–0.0588), followed by the mainland old stands (0.0188–0.0287) and the islands (0.009–0.0242).

Table 2 Comparison of genetic diversity parameters among nine Eastern white cedar (Thuya occidentalis L.) sites from three landscape types in Quebec, Canada
Full size table
Table 3 Population differentiation (pair-wise Fst) among nine Eastern white cedar (Thuya occidentalis L.) sites from three landscape types in Quebec, Canada
Full size table

The PCoA analysis revealed a clustering of sites, with the first and second axes accounting for 33.86 and 24.15% of the molecular variance, respectively (Fig. 2). A clear gradient appeared along the first axis, with the samples from the three fire refuges sites (f441, f442, f443) separated from the mainland and the island sites. The second axis mostly reflects the variance among sites, indicating a higher differentiation among the three fire refuges (f441, f442, f443). By contrast, sites from the mainland and from the islands were clumped together, except for f60 and is39 (which were closest to the mainland) (Fig. 2).

Fig. 2
figure 2

Principal coordinates analysis (PCoA) of genetic distances among nine sites of Eastern white cedar from three landscape types. Axis 1 explained 33.86% of the variation, and axis 2 explained 24.15%; Lake-islands landscape is39, is42, is134; Non-fragmented forest landscape (mainland sites): f60, f97, f23; Fragmented forest landscape (fire refuges): f441, f442, f443

Full size image

All loci were used in the Bayesian analysis. The most likely number of clusters (K = 4) was inferred by plotting L (K) and confirmed by plotting ∆K (Fig. S1 in Appendix 2). The four clusters were displayed in blue, yellow, orange, and green (Fig. 1). The patterns were consistent with the results of the PCoA and pair-wise Fst analyses. The geographic structuring of the clusters was weak, and no individuals formed a unique cluster (Fig. S2 in Appendix 2). The nine sites contained four genetic groups that were highly admixed. Individuals from fire refuges were more homogeneous (Fig. 1). There was a different pattern of the membership of clusters (Q) among the three small fire refuges (i.e., > 50% of Q: f441, orange 58.9%; f442, orange 48.3% and yellow 30.4%; f443, green 63.5%). Similar patterns of Q were observed in two of the islands (is42, blue 50.7% and yellow 21.2%; is134, blue 42.3% and yellow 26.8%) and two of the mainland sites (f97 blue 38% and yellow 29.5%; f60 blue 24.3% and yellow 39%).

Bottlenecks

According to BOTTLENECK and after corrections for multiple tests, there was no significant evidence for demographic changes in all sites under SMM. An excess in heterozygosity was detected in three sites (f442, f97, f23) under TPM. All sampled sites maintained alleles in frequencies expected for stable populations in mutation-drift equilibrium and showed a typical L-shaped allele frequency distribution.

Gene flow in the landscape

No significant correlation between geographical and genetic distances (Fst) was detected (R = 0.123, P = 0.252). The recent immigration analysis showed that two of the mainland sites (f97, f23) served as the main source of migrants (Appendix 3). Sites f441 and f442 received a large proportion of migrants from f23 (0.225 and 0.271, respectively). Sites receiving migrants from f97 were as follows: f443 (0.273), is42 (0.291), is134 (0.291), and f60 (0.294). The population is39 was a receiver from both f97 (0.136) and f23 (0.155). The mainland site f97 had the highest proportion of non-migrants (0.988).

Environmental association analysis

We tested whether mean fire interval and number of fire events were associated with EWC genetic structure. Regression approaches employed in GESTE found that runs with just a single environmental factor produce higher probability models when the mean fire interval was included in the model compared to when the factor was excluded (Table 4). Two-factor tests did not show the highest probability models to be those that included both factors and their interaction term in all cases.

Table 4 Results from GESTE analysis
Full size table

Discussion

The differences in genetic diversity among the three landscape types

Landscape features may either act as barriers to isolate populations or enhance dispersal and gene flow (Antolin et al. 2006; Manel et al. 2003). Negative values of the fixation index (F) and an excess of heterozygotes were found at most loci in EWC populations. A slightly negative Fis, indicating a slight excess of heterozygotes, is commonly reported in natural populations of conifer species. However, there was no evident landscape influences on this parameter. The most plausible explanation is that natural selection favours heterozygous individuals, a pattern commonly reported in many tree species (Hoebee et al. 2006; Oddou-Muratorio et al. 2004). The level of heterozygosity tends to increase in mature age classes, probably because of selection against homozygotes (Belletti et al. 2007).

The level of genetic differentiation (Fst) among EWC sites detected in this study was low and in the range of values reported for coniferous tree species (Kremer et al. 2012; Savolainen et al. 2007). However, the genetic differentiation of sites increased significantly from the lake-islands landscape (Fst = 0.015) and mainland sites (Fst = 0.023) to the fragmented fire refuges (Fst = 0.052). Our study also reveals a significant decrease in allelic richness from the non-fragmented (AR = 5.59) to the lake-islands (AR = 5.42) and the fragmented fire refuges (AR = 5.06). This result suggests that allelic richness is more sensitive to population size than other estimates of genetic diversity (e.g., heterozygosity), or that variation in this parameter is easier to detect.

Patterns of genetic structure, bottlenecks, and gene flow

Increased levels of genetic differentiation (Fst) and lower diversity (AR) across the three landscape types suggest that the forest matrix acts as a barrier to gene flow, especially in the refuges, which were nested in an early successional forest matrix as small fragments. This fragmentation probably creates barriers to the dispersal of pollen since the surrounding broadleaf forest canopy hinders gene movement between stands. This was supported by the association between Fst and the mean fire-free intervals detect with GESTE (Table 4). This association does not imply causality but may arise when habitat differentiation acts as a barrier to gene flow. The mean fire-free intervals are much longer in fire-refuges causing environmental isolation through time and higher genetic differentiation even in stands that are spatially close. In addition, it is plausible that the size of the EWC population in fire refuges has remained small over an extended time period as these stands escaped fires several times over the Holocene (Ouarmin et al. 2014). Rare alleles were lost possibly due to genetic drift, leading to a decrease in genetic diversity over time in terms of allelic richness (Widmer and Lexer 2001). Despite this possibility, no significant deviation from the mutation-drift equilibrium was detected, and no evidence was found that EWC fire-refuges were either severely reduced in their genetic potential or that they might have suffered from strong demographic declines.

Conversely, the absence of obstacles in the lake-islands landscape probably increases the effectiveness of both seed and pollen dispersal. Gauthier et al. (1992) have reported the presence of a similar level of genetic variability between lake-island and mainland jack pine populations. The level of gene flow in jack pine island populations was sufficiently high to counteract the effect of isolation. Boreal tree species are predominantly outcrossing, and pollen is wind-dispersed. Pollination distances that lead to successful mating (as determined by genetic parentage analysis) range between 3 and 100 km (Kremer et al. 2012). The mean seed dispersal distance of EWC seeds is shorter and is estimated at 45 to 60 m, with a maximum dispersal distance between 60 and 115 m (Asselin et al. 2001). However, conifer seeds and cones can glide more easily on the ice of the lake in winter than through densely forested areas. These life cycle and landscape characteristics facilitate the inter-stand gene exchanges among naturally fragmented EWC populations in an open lake-island landscape. A wind-pollinated, wind-dispersed tree within an open landscape is expected to exchange a substantial number of migrants (Bacles et al. 2006).

Although small EWC fire refuges are fragmented, there is still an exchange of genes occurring among the EWC sites as revealed by the absence of significant isolation by distance (IBD) between the genetic (Fst) and geographic distances. The level of gene flow among the EWC sites buffers the effects of fragmentation on the genetic characteristics of the small EWC fire refuges. The life cycle of EWC and the long juvenile phase of tree species may also have played a role. The presence of individuals from several age or size classes with overlapping generations is characteristic of EWC stands. The pattern of establishment of EWC is constant and showed a typical inverse J age structure (Bergeron 2000). Following a colonization event and before the first cohort of trees reaches reproductive age (during the first decades), there is a constant arrival of new migrants (seed flow) and no in-stand reproduction. Therefore, in a newly founded population, there is always a non-negligible part of the space that is already occupied by juveniles from seeds that arrived years earlier. In their study, Lesser and Jackson, (2013) show that alleles accumulated rapidly in ponderosa pine (Pinus ponderosa) populations following initial colonization and that contemporary levels of genetic diversity are formed early in the development of a tree population. At population sizes of approximately 100 individuals, allele accumulation is saturated. High levels of gene flow in the early stages of population growth result in a rapid accumulation of alleles and create relatively homogenous genetic patterns among populations.

The level of Fst measures historical population connectivity, and Bayesian assignment tests complemented it by revealing contemporary connectivity (over the last several generations) (Holderegger et al. 2010). This analysis reveals a source (f97, f23)—sink (other sites) dynamic that has rarely been identified in landscape genetic studies (Storfer et al. 2010). It supports the hypothesis that mainland EWC forests are the main source of gene flow. The clear source-sink directions are likely because that northwest winds are dominant during the time of pollen and seed dispersal (Denneler et al. 2008). Variation in the wind direction and surface wind speed affect the level of gene flow and ultimately the genetic structure of the populations (Nathan et al. 2011; Wright et al. 2008).

The clustering of EWC sites into four groups, as found in the STRUCTURE results, might indicate four evolutionary lineages recolonized during the last glacial maximum. However, this explanation seems less likely given the high level of genetic admixture of the majority of stands. In addition, from an evolutionary perspective, the EWC sites studied here are relatively young, as they colonized the area approximately 8000 years ago (Carcaillet et al. 2001). Furthermore, the nine sites were separated by only a few kilometres (< 4 km) and probably have the same biogeographic history. Hence, we presume that common ancestry combined with genetic drift in fire-refuges better explains the observed patterns.

Conservation value of fire residuals in fire-prone landscape

Habitat patchiness is recognized as one of the salient features of landscapes by population geneticists, ecologists and conservationists. Human activities and natural disturbances (i.e., fire), largely contribute to population fragmentation at the landscape scale. Forest remnants have been identified as historical or potential reservoirs of genetic diversity in tree species in fire-prone landscapes, for example; Betula maximowicziana (Uchiyama et al. 2006), Picea rubens (Mosseler et al. 2003) and Nothofagus dombeyi (Premoli and Kitzberger 2005). However, the populations need to be maintain at sufficient density and size to avoid the erosion of genetic diversity and the effects of inbreeding.

The mainland EWC stands maintained high genetic diversity (higher AR) and function as the main source of gene flow that ensures the maintenance of genetic variability in the small EWC fire refuges. Thus, the conservation of large mainland old growth EWC should be given careful attention, and should be protected effectively. Small fire refuges are also important because of their lower fire susceptibility and their ability to persist in the boreal forest landscape over centuries. Fire refuges have an intrinsic conservation value in landscapes affected by spatially heterogeneous fires and they are important for population persistence through disturbance events. These fire refuges may be particularly important for the long-term viability of EWC because they serve as a source of recolonization for areas more affected by disturbances.

Conclusions

Our data revealed a substantial level of gene flow in the EWC among different landscapes in the boreal mixed-wood. The lake landscape facilitates gene flow, whereas the fragmented terrestrial forest increases population differentiation in EWC. Both wind and landscape features shape the genetic structure of EWC. It is not possible to dismantle the respective effect of fires and landscapes. In the boreal forest, fire regimes are controlled by a strong interaction with landscape at the local scale (Bergeron 1991). The consequence is that on the islands, the fires of mixed severity and the proximity of the water promotes the survival of few EWC individuals nearby. EWC populations are re-established faster (lower differentiation, higher gene flow). On the mainland the EWC populations take time (often more than 200 years after fire) to recover. EWC can be eliminated by a fire event and have to return over long distances. In refuges, populations are small and can be potentially isolated for a long period (higher differentiation and relatedness). But they can also grow in the absence of fire around and maintain themselves in the very long term in the landscape. Our study provides evidence that the mainland EWC stands within a successional forest matrix possess a higher allelic diversity than the small EWC fire refuges. Our suggestions for forest managers and practitioners are that during harvesting, small-patch retention is good, but it is not enough for conservation of genetic diversity in EWC. Large-scale management plans that lead to boreal landscapes containing (relatively) large tracts of old-growth forest are desirable (Bergeron and Fenton 2012). Fire-refuges must also be given special attention. These sites are particularly important in a context where climate change is predicted to increase the burn rate over the boreal forest by the end of the twenty-first century (Girardin and Mudelsee 2008). The cumulative impact of increasing the fire incidence and the current rates of harvesting, which still prevail in this region, contributes to the erosion of ecological resilience by reducing ecosystem variability in time and space. Fire refuges are small, dispersed in the forest matrix and are likely to be lost. The conventional models of forest protection might not be appropriate, and fire refuges could benefit from protection even if there is no accumulation of genetic diversity or rare alleles because of their continuity. Such areas are important for EWC population persistence under long-term climatic change because of the lower susceptibility of these sites to fire events (Wilkin et al. 2016).