Resilience of the boreal forest in response to Holocene fire-frequency changes assessed by pollen diversity and population dynamics

Additional keywords: climate, fire, Holocene, lake sediments, numerical analysis, pollen diversity.

Fire is a global process (Flannigan et al. 2009) and a key ecological factor impacting on many of the world’s ecosystems (Bond et al. 2005), including the boreal forests (Johnson 1992). However, the function of fire-frequency changes on plant diversity and community structure remains to be understood (Krawchuk et al. 2006) on multi-century scales by using high-resolution paleofire and high-vegetation analyses. Indeed, plant community diversity can be altered by changes in fire frequency, whereas the transformed vegetation can provide a different fuel composition that subsequently feeds new fire regimes. Furthermore, the increase in fire frequency expected with global warming (Flannigan et al. 2009) might move the boreal ecosystems outside Holocene conditions (Girardin and Mudelsee 2008). Indeed, the current global warming influences natural disturbances such as fire (Bergeron et al. 2010), and thus plant dynamics through fire–climate interactions (Overpeck et al. 1990; de Groot et al. 2003).

Over the recent millennia, fire frequency has changed several times in the North American boreal forests (e.g. Clark 1990; Lynch et al. 2004), as well as in other boreal regions (Clark et al. 1998; Carcaillet et al. 2007; Katamura et al. 2009). Because these fire events can cover a range of situations that are not encountered during the period of direct observation (20th–21st centuries), it is necessary to examine the long-term plant responses to past fire frequencies. Ultimately, this approach can provide potential scenarios for future transformations of boreal forests (hereafter ‘BF’) and help to preserve the ecological integrity of these ecosystems (Landres et al. 1999).

Today, the resilience of BF to changes in fire frequency is poorly known and a greater understanding is required because these forests are to be sustainably managed (Pastor et al. 1998; Chapin et al. 2004). Indeed, the resilience of boreal ecosystems remains a major societal question in light of global current and future changes in fire frequency, let alone for harvesting practices, as this biome stores a third of global terrestrial carbon (Apps et al. 1993). Although most studies focussed on short-term community dynamics (e.g. Haeussler et al. 2004; Rydgren et al. 2004), very few studies have attempted to examine the BF community–disturbance complex over centuries or millennia. Ecosystems can react to modification of the fire regime by continuous or discontinuous changes in plant composition. However, the magnitude of fire-frequency changes is crucial for determining whether communities will change abruptly or progressively. The discontinuous transformation of communities is linked to the notion of thresholds (Hooper et al. 2005). Although BF are generally resilient to a single fire after 150 to 200 years of tree dynamics (e.g. Bergeron 2000), the abrupt change in mean fire interval (MFI) from e.g. 100 to 500 years, or vice versa, could have consequences on the composition of communities by suppressing or facilitating tree growth, or transforming one ecosystem into another (Le Goff and Sirois 2004; Jasinski and Payette 2005).

In the Mixedwood (southern) BF of Québec, the MFI changed abruptly ∼2200–2000 years ago (Carcaillet et al. 2001; Cyr et al. 2009). A long MFI occurred between 6800 and 2200 calibrated years before present (hereafter ‘cal BP’), but MFI decreased after 2200–2000 cal BP in the Mixedwood, likely resulting from more frequent intrusions of dry air masses from the Pacific (Carcaillet et al. 2001). In the Coniferous (northern) BF, the fire regime was characterised by longer to shorter MFI until ∼3000 cal BP, followed by longer intervals until ∼1000 cal BP, and the longest MFI since 800 years ago (Ali et al. 2009a). This increasing MFI results from rising available moisture in the region, ultimately linked to reduction in summer insolation (millennial-scale orbital forcing; Hély et al. 2010).

Here, we re-examine the links between changes in MFI and vegetation inferred from two published pollen datasets. The pollen-frequency diagram did not demonstrate any relationship between pollen-inferred vegetation and fire frequency (Carcaillet et al. 2001), whereas studies in western North American BF indicate that MFI and vegetation composition are linked (Hu et al. 2006; Higuera et al. 2008). Because the former use of pollen data was intuitive, we stress the need for statistics and ecological indices to pinpoint whether vegetation and fire are related, especially when pollen-assemblage changes are subtle (Genries et al. 2009). Indeed, several diversity indices can help to detect the pollen-inferred vegetation transformation. We also consider plant-population dynamics by the use of a standardised pollen influx combined with Principal Component Analysis (PCA). PCA provides a temporal summary of species behaviour, and the standardisation of the pollen influxes provides a quantification of population density changes (Fuller 1998) or productivity (Seppä et al. 2009). Finally, the pollen trajectories are analysed to determine the timing of vegetation convergence. The similarity of vegetation history between the sites compared with the timing of fire frequency allows investigation of the long-term effects of fire-frequency on plant cover, along with the vegetation’s resilience towards major local changes. We hypothesise that long MFI during the mid-Holocene would have promoted late successional species such as Abies balsamea (L.) Mill., Thuja occidentalis L. and Picea glauca (Moench) Voss. Similarly, the short MFI that occurred during the late Holocene should have favoured early successional plant species such as Populus tremuloides Michx., Betula papyrifera Marsh., Alnus viridis ssp. crispa (Ait.) Turril, and Pinus banksiana Lambert.

The study was conducted in western Quebec, ∼300 km south of James Bay. Lac Francis and Lac à la Pessière are headwater lakes, located on eskers that are partly embedded in clayey lake sediments. The ‘Northern Clay Belt’ is the region of the Proglacial Lake Ojibway, which drained abruptly when the residual ice-sheet collapsed, allowing a northward outflow to Hudson Bay at 8470 cal BP (Barber et al. 1999). Subsequently, organic sediment accumulated in the small lakes that were left behind. Despite the location of the study sites on eskers, which are atypical for the Clay Belt, present-day statistical similarity in the length of stand-replacing fire cycle on esker, rock, and clay sites (Bergeron et al. 2004) allows us to extend our results to the entire region.

Lac Francis (hereafter L. Francis, 48°31′35″N 78°28′20″W; 305 m above sea level (ASL)) is in the Mixedwood BF, and Lac à la Pessière (L. Pessière, 49°30′30″N 79°14′25″W; 280 m ASL) is in the Coniferous BF. The lakes are small, strongly coloured, dimictic closed basins. The maximum water depth is 6 m in L. Francis, and 16 m in L. Pessière. The water columns develop a strong thermal stratification at 1.5-m depth in L. Francis and 4 m in L. Pessière, related to the size of the lakes: 0.8 ha for L. Francis and 4.0 ha for L. Pessière; anoxia is reached at ∼3 and 7 m respectively. Anoxic conditions promote sediment stratification, and hence the organic sediments are laminated, except for those deposited during the last millennium in L. Francis.

The vegetation in the Mixedwood part of the Clay Belt is a mixed closed-canopy forest, characterised by Abies balsamea and Betula papyrifera. Around L. Francis, the canopy is dominated by Picea glauca and Abies balsamea on the mesic soils. Species at their northern limit around L. Francis include Pinus strobus L. and Pinus resinosa Ait. (Bergeron et al. 1997). Stands of certain species (e.g. P. strobus, Thuja occidentalis, A. balsamea, and Acer spicatum Lam.) are significantly fragmented just a few tens of kilometres north of L. Francis, at ∼49°N (e.g. Messaoud et al. 2007). This latitudinal change of population structure corresponds to the transition between the Coniferous and Mixedwood BF and has not been in equilibrium with the climate since at least the Little Ice Age (Hofgaard et al. 1999). Furthermore, fire frequencies have not been different in the Coniferous and Mixedwood BF for at least the last 300 years, suggesting that current fire regimes cannot explain the current vegetation pattern (Bergeron et al. 2004).

Coniferous BF is dominated by Picea mariana on wet and mesic soils and by Pinus banksiana with Cladina lichen on dry-sandy soils. L. Pessière is surrounded by old-growth Picea mariana forest on wet soils and by clearcuts on mesic soils, where Betula papyrifera, Populus tremuloides, Picea mariana and Abies balsamea regenerate. In the Coniferous BF, A. balsamea populations are located in areas that are not exposed to frequent fires, such as islands or lake shores (Sirois 1997), and have been sensitive to increasing MFI or fire severity during the Holocene (Ali et al. 2008).

Coring techniques, and laboratory treatment for pollen and charcoal preparation as well as the chronologies obtained are detailed in Carcaillet et al. (2001). The following numerical assessment of the pollen data results in a summarisation of the timing and magnitude of palynological changes at each site, and a comparison of the differences between the two sites through time. Nevertheless, since Carcaillet et al. (2001), the numerical methods for revealing paleofires have markedly improved (Gavin et al. 2006; Higuera et al. 2008; Ali et al. 2009b), although the main principles are conserved, i.e. decomposing the sedimentary charcoal records into low-frequency and peak components, and based on statistics, selecting the most significant charcoal peaks to reconstruct the fire events. Consequently, we were able to reconstruct the local fire history from the previous charcoal raw series using the program CharAnalysis.

The peak component is assumed to be composed of two subpopulations, namely noise that represents variability in sediment mixing, sampling and analytical and naturally occurring noise, and fire events, representing charcoal input from local fires. For each sample, we used a Gaussian mixture model to identify the noise distribution. A 99th percentile of the noise distribution was used as a threshold for separating samples into fire and non-fire events. The low-frequency components are estimated with a locally weighted regression using a 400- or 500-year window. For each record, we chose the window width that maximised a signal-to-noise index and the goodness-of-fit between the empirical and modelled peak distributions. Fire history was described by quantifying the variation of fire return intervals (FRI; years between two consecutive fire events) over time, smoothed using a LOWESS with a 1000-year window.

The rate of change (ROC) is a dissimilarity measurement based on a Euclidian distance that quantifies the magnitude of palynological changes per time unit between two adjacent pollen spectra (Jacobson and Grimm 1986). Here, the dissimilarity is based on squared-chord distance (Birks 1977; Gavin et al. 2003) computed on the square root of pollen frequencies. ROC is a powerful method to detect rapid vegetation change (Jacobson and Grimm 1986; Shuman et al. 2005) such as those expected with changes in fire frequency. Different estimations of the ROC were carried out, first with a taxa list that took into account the full pollen dataset (i.e. n_taxa = 37) and, second, with a selected taxa list that included only those with a mean frequency greater than 1% (i.e. n_taxa = 18).

The pollen diversity varies according to the vegetation composition and structure surrounding the study site (Jackson and Lyford 1999; Broström et al. 2005). If the vegetation changes in terms of plant density or landscape structure, the processes of pollen dispersal and accumulation will be altered, thus affecting the pollen diversity. Changes in landscape or community structure may be due to changes in the disturbance regimes or in climate. In addition, changes in pollen richness may provide a means of quantifying changes in vegetation diversity through time, which can be assessed by ‘rarefaction analysis’ (Birks and Line 1992), an index based on the calculation of the expected number of taxa (E(T_n), Eqn 1) in a random sample of n individuals from a larger population of N individuals containing T taxa (Heck et al. 1975).

Here, rarefaction analysis is carried out on the complete set of terrestrial pollen. Studies of the use of rarefaction analysis on pollen data have concluded that the results are controlled by the size of the pollen catchment (Seppä 1998), and by the vegetation composition or structure if the catchment size has not changed (Birks and Line 1992; Seppä 1997). Here, we assume that the size of the pollen catchments has been constant since the onset of organic matter accumulation in the lakes, following the drainage of the residual proglacial waters. Thus, all changes identified in E(T_n) should result from changes in vegetation composition or structure (open or closed, ±fragmented, etc.). The E(T_n) and the ROC were carried out with 3Pbase (© Centre National de la Recherche Scientifique; Guiot and Goeury 1996).

Two other biodiversity indices are used here to estimate variations in the pollen-inferred vegetation composition: Simpson’s (D, Eqn 2) and Shannon’s (H, Eqn 3) indices consider the species number and the abundance. The species abundance is assumed to be correlated with the pollen abundance (e.g. Bradshaw and Webb 1985), and if not, the relationship between pollen and plant abundance remains broadly stable through time. The indices help to compare the structure of those pollen assemblages in which the record of each species is comparable through time, and are calculated here for all the terrestrial taxa. D is computed by determining the proportional contribution of each taxon to the total number of pollen in the whole spectrum:

where p_i is the pollen frequency for the ith taxon. Shannon’s index (H) is calculated as:

Compared with Simpson’s index, Shannon’s index has the advantage of reducing the importance of the abundant taxa through the log-normal transformation of frequencies. Consequently, both indices must be used together to reveal the role of abundant and rare taxa in the diversity scores. The temporal variability in pollen diversity (H and D), was assessed by the Kruskal–Wallis test (KW-test) after averaging the pollen diversity scores per millennia, e.g. 5000 cal BP corresponds to 5500–4500 cal BP.

PCA was used to assess changes in biological assemblages. The best discrimination was obtained with the square roots of 17 terrestrial taxa with a mean frequency above 1% over the entire sequence. The 17 terrestrial taxa summed to 100%.

Two series of PCAs conducted on a correlation matrix were carried out, on each pollen series separately, then on the two series together. The first allows determination of pollen trends at each site independently. The second permits a comparison of the temporal trajectories of the two series on the same basis and at the same time scale (Genries et al. 2009). Only for the joint PCA are the two series transformed into equal time intervals of 100 years. Correlations coefficients (r) between assemblages were computed at the same time to estimate the similarities or dissimilarities between sites.

According to experiments, only small lake basins (∼1 ha) register any local forest disturbances when the pollen signal is expressed in terms of pollen accumulation rate or influx (Koff et al. 2000). To decipher the changes in the main plant populations, pollen influxes were standardised to allow a comparison of the magnitude of changes in pollen accumulation rates (PAR) independently of the pollen productivity differences among taxa and the taphonomical properties of each taxon (Fuller 1998). Each taxon is displayed individually at a scale that is the same for all, allowing an estimate to be made of the magnitude of pollen population changes within and among pollen taxa. However, contrary to Fuller (1998), the rescaled pollen influx (PAR′_i) is here based on a minimax transformation (from 0 to 1), such as:

where PAR′_i is the standardised value of the ith sample PAR_i, and PAR_min and PAR_max are the minimum and maximum of the PAR′_i. After transformation, the data were smoothed.

The mean pollen diameter of Betula is related to the species (Richard 1970) and allows changes in Betula populations to be interpreted (Marcoux and Richard 1993; Wolfe et al. 2000; Fréchette et al. 2006). Any changes in birch populations are interesting in discussing forest community transformation because of the ecological diversity in the Betula species. The main regional birches are: B. alleghaniensis, which is a mature forest species from temperate Mixedwood plains to Mixedwood BF; B. papyrifera, which is a typical species from Mixedwood and Coniferous BF; and shrubby birch (e.g. B. glandulosa), which abounds in mires and in subarctic and subalpine areas. The equatorial diameter of 100 pollen grains of Betula was measured in each spectrum. Betula pollen size was determined as the diameter from the pore tip to the outer surface of the opposite wall.

The trends in the new reconstructions of fire intervals (Fig. 1c, d) are not different from the old reconstructions (Fig. 1a, b), despite differences in the methods to decompose fire events. The L. Francis record shows (Fig. 1c) three periods of fire intervals, with the first before 6000 cal BP when the MFI was short (230 ± 149 years). However, only two fire intervals are recorded, which is poor for comparisons. The main transition occurred at 2200 cal BP, splitting a 4000-year-long period of longer MFI (340 ± 158 years, n = 12), followed by a 2000-year-long period of shorter MFI (192 ± 124 years, n = 12). The 2200-cal BP transition was identified both in the old (Fig. 1a) and in the new reconstruction (Fig. 1c).

Fig. 1.  Mean fire intervals (mean ± s.d., n) inferred from local fire history at Lac Francis, Lac à la Pessière, based on high-resolution analyses of charcoal influx. The old reconstruction of fire intervals (a, b), from Carcaillet et al. (2001), is compared with the new reconstruction based on the up-to-date methods (c, d). The vertical lines indicate changes in local fire frequency. The two new reconstructions (c, d) are rescaled to allow intersite comparison (e) independently from local mean of fire intervals related to vegetation, topography or landscape features.

At L. Pessière, four fire periods are noted (Fig. 1d), as in the old reconstruction (Fig. 1b). However, the first MFI transition is later at 6700 cal BP (Fig. 1d), whereas it was at 7300 cal BP previously (Fig. 1b). The first MFI period (7600–6700 cal BP) shows nine intervals (96 ± 102 years) followed by a 3400-year-long period with longer MFI (166 ± 102 years, n = 20). The second transition at 3300 cal BP matches the previous reconstruction. However, the latest transition at 1000 cal BP is ∼200 years later than in the old reconstruction. Between 3300 and 1000 cal BP, the MFI is shorter (125 ± 98 years, n = 17), followed by the longest historical MFI (182 ± 126 years, n = 6) since 1000 cal BP (Fig. 1d).

The rescaling of the MFI (Fig. 1e) allows an intersite comparison, and reveals a more or less similar pattern of fire history, whatever the local MFI, which are longer around L. Francis compared with L. Pessière. This rescaling indicates that regionally the fire intervals were short before ∼6000 cal BP, as well as between ∼3000 and 1000 cal BP. In contrast, between 6000 and 3000 cal BP and since 1000 years ago, the MFI were longer (Fig. 1e).

The ROC is based on summed data for all terrestrial pollen taxa as well as terrestrial pollen taxa with a frequency greater than 1%. Both approaches produce similar patterns, which are highly correlated for each pollen series (r² = 0.95; P = 0.0001). This indicates that most of the variation between assemblages is due to the most abundant taxa, and that taxa accounting for <1% of the total are insignificant here in explaining the dynamics of pollen-inferred vegetation.

The single elevated ROC at L. Pessière and the first at L. Francis correspond to the afforestation phase and can be dated at 7500 and 6700 years ago respectively (Fig. 2a, b). Before 5500 cal BP, the ROC at L. Francis exhibits fluctuations (instability) resulting from vegetation transformation during the millennia following the afforestation (Fig. 2a). After 5000 cal BP, L. Francis experienced 3000 years of low ROC values with insignificant fluctuations (stability). The main late-Holocene instability occurred between 2000 and 1300 cal BP. This period of higher ROC at L. Francis occurs after the observed decrease in MFI, i.e. 2200 years ago (Fig. 2i). For the past 1200 years, the ROC remained low at L. Francis (Fig. 2a), although the MFI tends to increase up to ∼250 years (Fig. 3i).

Fig. 2.  Indices of vegetation change calculated from the pollen data from Lac Francis and Lac à la Pessière. (a, b) Rate of change (ROC). Three diversity indices: (c, d) Shannon’s index (H); (e, f) Simpson’s index (D); and (g, h) the rarefaction analysis (E(T_n)). The solid horizontal lines show the mean values, and the dashed lines the associated standard deviations. The histograms show the mean values of D, H and E(T_n) computed per millennia; the letters refer to differences based on the Kruskal–Wallis test (different letters mean significant differences). The arrows indicate the observed changes in pollen assemblages, and the grey vertical lines indicate changes in local fire frequency observed at Lac Francis and at Lac à la Pessière.

Fig. 3.  Scores of the principal component analyses (PCA) based on pollen data plotted against the time axis. Axes 1 (a and b) and 2 (c and d) show gradual trends for the two lakes (no abrupt change), suggesting a lack of response to fire frequency change. The grey vertical lines indicate changes in local fire frequency observed at Lac Francis and at Lac à la Pessière. Values correspond to mean fire intervals. Loading ordination of pollen taxa on axes 1 and 2 at (e) Lac Francis, and (f) Lac à la Pessière.

At L. Pessière, the ROC has been very low and stable for the past 7200 years, suggesting that very little vegetation transformation has occurred in the Picea mariana forest (Fig. 2b). Clearly, some small changes are observed during the last 3000 years at L. Pessière, but never of the magnitude of those recorded at L. Francis. The first main pollen modification since afforestation appears between 2700 and 2500 cal BP, and the second between 700 and 300 cal BP. Both are delayed by 400–700 years compared with the recorded changes in fire frequency dated at 3400 and 1000 cal BP (Fig. 2j).

Although the ROC distribution shows some obvious changes, the diversity indices (Fig. 2c–f) and the rarefaction analysis (Fig. 2g, h) have remained quite stable since the afforestation that ended at ∼6500 and ∼7500 cal BP at L. Francis and L. Pessière respectively. At L. Francis, no change in diversity (D, H) was observed (Fig. 2c, e), whereas a little increase in the richness (E(T_n), KW-test P < 0.05) is revealed between 3500 and 1500 cal BP (Fig. 2g), which uncovers the main period of elevated ROC (Fig. 2a). At L. Pessière, two decreases in diversity (D, H) occurred 6500 and 4600 years ago (KW-test, P < 0.05; Fig. 2d, f), and only one in terms of richness (E(T_n), P < 0.05; Fig. 2h).

As there were no major abrupt or long-term changes, according to the rarefaction analysis and diversity indices at L. Francis, there appears to have been no significant modification of pollen input at the lake after the afforestation phase, either in terms of pollen catchment or vegetation structure whatever the MFI (Fig. 1c). These results suggest that the small palynological change recorded by ROC at L. Francis did not affect the basic composition and structure of the vegetation (richness and diversity) and is independent of fire history, but might be based on plant abundance shifts in the forest communities. Conversely, whereas the pollen diagram of L. Pessière was very stable over the Holocene (Carcaillet et al. 2001), the diversity indices (D, H) indicate that the structure of the vegetation was significantly transformed. However, these transformations occurred during periods of stable MFI (Fig. 1d).

The two first axes of the PCA account for almost the same amount of variation (% of variance) for L. Francis and L. Pessière (Fig. 3a–d, Table 1), whereas the visual observation of the two pollen diagrams (in Carcaillet et al. 2001) indicates larger temporal modifications at L. Francis compared with L. Pessière. The first axis for each lake, accounting for 33 and 29% of the variation respectively shows a strong contrast between the initial vegetation at the time of the retreat of residual proglacial waters and the late Holocene (Fig. 3a, b). No abrupt transition is observed, suggesting that a multimillennial, long-term and gradual process dominated the vegetation dynamics around the two lakes. The changes in MFI (Fig. 3e, f) did not significantly affect either the short- or long-term dynamics, although they coincide with the nil values on axis 1 at L. Francis.

Table 1.  Variance (%) and eigenvalues of the three principal component analyses (PCA)

The second PCA axis, which accounts for 17% of the variation at both sites, indicates that the late-Holocene assemblages tended to resemble those from the early Holocene (Fig. 3c, d). At L. Francis, the pollen trend highlighted by the second axis appears associated with the change in MFI dated from 2200 cal BP (Fig. 3c). At L. Pessière, the change in MFI at ∼3400 cal BP is independent of the vegetation trends, but those at ∼1000 cal BP precede the observed modification in vegetation composition (Fig. 3d). Since ∼800–600 years cal BP, the two sites have experienced a synchronous transformation, highlighted on the second PCA axis, which suggests a broad-scale rather than a local-scale process linked to fire regime.

The third PCA axes of the two pollen series, accounting for ∼10% of the variation (Table 1), identify a short period, covered by two to three samples, different from the rest of the pollen series. There is no association here between the vegetation trends and the fire-frequency changes, or the long-term postglacial dynamics. Furthermore, their eigenvalues are less than the mean eigenvalues for each PCA (Table 1), suggesting that these axes are less reliable than the higher-order axes (Legendre and Legendre 1984).

Fig. 4 shows the main results of a PCA performed on the matrix of the two combined pollen series, with time intervals standardised to 100 years. This PCA allows a comparison of the main variation (axes 1 and 2) of the two pollen trajectories, on the same pollen dataset and at the same periodicity. The L. Pessière pollen assemblage was similar to L. Francis 6800 years ago (Fig. 4). After the early-mid Holocene, the two sites have distinct assemblages and diverge during the rest of the mid-Holocene. The L. Francis pollen data are located on the negative sides of axes 1 and 2, whereas the L. Pessière data are on the positive sides of these axes. The late Holocene is characterised by a transformation of the L. Francis assemblage, towards a pollen composition similar to L. Pessière. The observed pollen trajectories, based on the two first axes, cannot explain the full extent of the similarity or difference, because only 45% of the total variation is accounted for by these two axes (Table 1).

Fig. 4.  Principal component analyses (PCA) of the pollen series transformed to equal time intervals of 100 years (above), and the species contribution (below). The variance of axes 1 and 2 is 26 and 20% respectively. The assemblages from the early–middle Holocene are clustered in the two series, and the late Holocene displays similar trends over the last 1200 years.

The correlation coefficient r between the pollen spectra from L. Pessière and L. Francis at the same age is a more comprehensive index of similarity between the two series than PCA axes 1 and 2 together, because it encapsulates 100% of the total variation. The r was calculated using the same list of 17 taxa in both pollen series, expressed as percentages (same datasets as for the PCA). The r scores plotted against time clearly indicate that three main sequences can be recognised (Fig. 5). Between 6800 and 6300 cal BP, the two series were relatively similar; between 6000 and 1200 cal BP, the pollen rain around the two sites was distinct (P = 0.01), whereas the two pollen series display the highest level of similarity during the last 1200 years (r > 0.6, P = 0.01; Fig. 5).

Fig. 5.  Correlation coefficients between synchronised pollen assemblages. The series have been transformed to equal time intervals of 100 years before calculation. For the last 1200 years, the pollen series are highly correlated between the sites, whereas during the middle Holocene, the forests were different.

The positive side of the PCA axis 1 is characterised by the trees currently dominant in the BF, i.e. Picea type (t.) mariana and Abies, and the negative side by Populus t. tremuloides, Cyperaceae, Taxus, Pinus strobus and ‘other trees’ (Fig. 4). Axis 1 divides the assemblages of the afforestation phase (early Holocene on the negative side) from those of the late Holocene (positive side). The second axis places mid-Holocene dominant trees and shrubs, such as Thuja/Juniperus, Alnus t. incana, Pinus strobus, Betula and Taxus on the negative side, against a group of plants on the positive side that occurred during the early Holocene and once again during the last 1200 years, e.g. the trees Alnus t. viridis, Pinus t. banksiana and Picea t. mariana, and other herbs (Fig. 4).

Plotting the standardised pollen influxes against time illustrates the plant-population changes. In Fig. 6, the taxa are grouped by affinities based on the PCA (Fig. 4), which illustrates the similarity of the individual species behaviour to environmental changes.

Fig. 6.  Standardised pollen influxes. The taxa are grouped per sites according to their mutual correlation, resulting from the different principal component analyses (PCA). The grey vertical lines indicate changes in local fire frequency observed at Lac Francis and at Lac à la Pessière. Values correspond to mean fire intervals.

The first group contains taxa such as Populus (t. tremuloides + t. balsamifera), Pinus strobus, Taxus, and Thuja/Juniperus, which abound a few decades after the drainage of the proglacial waters at each site (Fig. 6a). The high pollen values of Thuja/Juniperus observed here are likely due to T. occidentalis (Liu 1990; Richard 1995). The two sites show the same pattern, with an immediate expansion of Populus, P. strobus and Taxus during the afforestation. The dominance of Populus has never been as high again as it was during the early Holocene. P. strobus and Taxus quickly expanded with Populus, but they maintained elevated pollen influxes through the mid-Holocene, followed by a slight decline up to the present day (Fig. 6a). The smallest P. strobus and Taxus influxes occurred over the past 4600 years at L. Francis and the past ∼1200 years at L. Pessière. Thuja shows a maximum influx between 6100 and 4600 cal BP at L. Francis and between 6600 and 4600 cal BP at L. Pessière (Fig. 6). The diachronic expansion of Thuja during the early Holocene is easily explained by the asynchrony of the afforestation at the two sites. However, the similarity of the timing and pattern of the Thuja decline suggests the effect of a broad-scale environmental process (Fig. 6). In both sites, the first change in MFI seems associated with the increase in MFI, with a progressive decrease of Populus, P. strobus and Taxus before the rise in MFI, and an increase in Thuja afterwards (Fig. 6a). Conversely, no similar pattern occurred during the last 1000 years, while the MFI was increasing. This feature indicates that the change in MFI is not a property linked to the dynamics of this grouping of tree species. Furthermore, the empirical model of vegetation dynamics (Bergeron and Dubuc 1988; DeGrandpré et al. 1993) indicates that Taxus should abound when the MFI is minimal, and vice versa, whereas the pollen influx of Taxus is maximal with short MFI and decreases when MFI are long. No obvious relationship is characterised between this first group of taxa and the changes in MFI (Fig. 6a).

The second group of taxa includes Betula, Alnus t. incana, other broad-leaved trees, and herbs, including Cyperaceae (Fig. 6b), i.e. taxa generally abundant in early successional communities. The dynamics of these taxa in relation to fire history is rather complex. When MFI is short (before 6000 and after 2200 cal BP at L. Francis and, before 6700 and after 3300 cal BP at L. Pessière), high influxes are achieved gradually over several centuries (Fig. 6b). The relationship between the dynamics of this second group of taxa appears partially explained by the MFI pattern (Fig. 6d).

The third group of taxa includes Alnus t. viridis and conifers typical of eastern BF, i.e. Pinus t. banksiana, Picea (type mariana plus glauca), and Abies (Fig. 6c). Although they have been present since the early afforestation, their influxes increased only from 1600 years ago at L. Francis, i.e. 600 years after the observed fire-frequency increase. At L. Pessière, P. t. banksiana and A. viridis were important up to 6800 years ago, then they declined slightly during the period of long fire intervals, finally expanding again between 3000 and 700 cal BP in association with Picea and Abies during the period of short MFI (Fig. 6c). The pollen abundance of this group at both sites was always greatest during the periods of short MFI. Although the exact timing of the decline and expansion of pollen influxes of these taxa does not perfectly match the timing of the fire-frequency changes, it is likely that a relationship does exist, but with a delayed response. The broad-scale history of Picea, Abies, P. t. banksiana, and A. t. viridis is variable, suggesting that local-scale factors affect these taxa.

The birch populations display several changes, clearly distinguished by the ROC calculated on the basis of Betula pollen diameter (Fig. 7a), and by the function describing mean Betula pollen diameter per assemblage (Fig. 7b). A PCA, performed on the 18 size-classes of Betula diameter (19–36 µm) indicate three stages at L. Francis. The first, between 6800 and 4600 cal BP, is characterised by pollen of both small (18–19 µm) and large (29–37 µm) diameter. A diameter of 18–19 µm corresponds to dwarf Betula (i.e. B. pumila L. or B. glandulosa Michx.) and 29–37 µm to both B. cf. papyrifera (pro parte) and cf. alleghaniensis (Richard 1970). The second stage is characterised by a larger mean pollen diameter, with grains of 29 to 37 µm (B. cf. alleghaniensis). The third stage, since 1700 cal BP, is associated with a smaller mean diameter, with pollen mostly between 20 and 28 µm (dwarf Betula and B. cf. papyrifera) (Fig. 7b). B. cf. papyrifera was likely not absent following the afforestation phase, but the presence of this taxon is only statistically significant over the last 1700 years.

Fig. 7.  (a) The rate of Betula change (ROC_Betula) against time. ROC_Betula is calculated from changes in the frequency of pollen grains in each diameter class. (b) Mean Betula pollen diameter plotted against time. The solid horizontal lines show the mean values over the Post-glacial, and the dashed lines the associated standard deviation. The thin vertical lines show the size-class limits deduced from the principal component analyses (PCA). The grey vertical lines indicate changes in local fire frequency observed at Lac Francis and at Lac à la Pessière.

At L. Pessière, only two stages were recognised, the first before 1800 cal BP and the second after this date. The first stage is characterised by large pollen grains (27–36 µm) corresponding to B. cf. papyrifera (pro parte) and to B. cf. alleghaniensis. The second stage is associated with diameter-size classes between 19 and 25 µm, i.e. B. cf. papyrifera (pro parte) and dwarf Betula. The second Betula stage is also marked by a high variation in pollen diameter compared with the earlier stage (Fig. 7b).

There is no apparent relationship between the timing of fire changes and that of Betula types. The increase in the ROC_Betula occurred at approximately the same time at the two lakes (1900 cal BP at L. Francis; 1800 cal BP at L. Pessière; Fig. 7a), suggesting a broad-scale process. At present, B. alleghaniensis is rare in the region of L. Francis (isolated individuals), and does not occur near L. Pessière, whereas B. papyrifera is abundant in both mixed and coniferous BF. Around L. Pessière, there are currently many large areas of swamps and swamp-forests where B. glandulosa is common.

The new paleofire reconstructions (Fig. 1), following the up-to-date protocol of Higuera et al. (2009), reveal greater fire frequency than the earlier analysis employed (Carcaillet et al. 2001). The new interpretation is much more rational in terms of statistical analyses to distinguish charcoal peaks from the noise, which is illustrated by a Gaussian mixture. More fires are observed during periods of long MFI, thus increasing the frequency, whereas in contrast fewer fires are noted during periods of short fire intervals, resulting in a decrease in frequency. Finally, the new reconstructions provide results that better match the modern range of MFI in the area (Bergeron et al. 2004).

In the long term, the vegetation has changed slightly, but less than expected on the basis of observed local fire history (Fig. 1) and of the empirical fire–vegetation model deduced from stand structure analyses in relation to the time since the last fire (e.g. Bergeron 2000; Gauthier et al. 2000). Delays in the response of fire-sensitive species (Pinus banksiana and Picea spp.) were recorded at L. Francis, along with corresponding changes in the abundance of fire-susceptible species such as Abies (Fig. 6c). Furthermore, vegetation indices and PCA do not highlight any systematic changes associated with changes in MFI (Fig. 3a–d). The shortest MFI between 2200 and 1000 cal BP at L. Francis and 3300 and 1300 cal BP at L. Pessière might have influenced the spread of Betula and herbs to the disadvantage of conifers. Finally, the dynamics of pollen taxa are often synchronous at the two sites, following the same general patterns independently of MFI. Examples of species illustrating this similarity include Pinus strobus, Populus, Taxus, Thuja, and ‘other broad-leaved’ species (Fig. 6), as well as Betula types (Fig. 7). Together, these observations suggest that the magnitude of the MFI changes was not sufficiently strong to transform the plant cover over the long term. Species might have been influenced, but with a delay, e.g. Picea and P. t. banksiana (Fig. 6c).

The eastern BF appears resilient to changes in fire intervals. The diversity of these forests could be an argument to explain the stability of the ecosystems facing changes in the fire regime. Conversely, the postglacial Alaskan BF, which is less diverse in arboreal taxa, experienced marked transformation of its forest communities in relation to Holocene fluctuations in MFI (Lynch et al. 2004; Hu et al. 2006; Higuera et al. 2009).

The resilience of the vegetation facing changes in MFI suggests that the current increase in fire intervals in eastern Canada, associated with global warming since the mid-19th century, would not result in significant changes in the vegetation composition, contrary to previous assertions (e.g. Weber and Flannigan 1997). However, logging could create new historical pathways (Payette and Delwaide 2003). This increase in MFI results from a reduction in the frequency of summer droughts (Bergeron and Archambault 1993) and a change in air mass circulation favouring elevated summer moisture in the Mixedwood forest (Girardin et al. 2009).

Although the forests since the time of afforestation can be described as ‘boreal’, their composition was different from that observed today, as is demonstrated by the trajectories in the PCA (Figs 4, 5). The forests have never returned to their initial composition, despite similarities between the early-mid- and the late Holocene at each site. However, this convergent trend represents at most only 17% of the total variation in the pollen dynamics (Fig. 3c, d; Table 1). The difference between the early postglacial and present-day forests at L. Pessière is smaller than that at L. Francis, according to the trajectories recorded by PCA (Fig. 4).

The original woodland differed from that at present in that Picea, Abies (Fig. 4) and Betula papyrifera (Fig. 7b) were less common, and Pinus strobus, Taxus, Thuja (Fig. 4), Betula cf. alleghaniensis (Fig. 7b) and ‘broad-leaved trees’ were more common. This applied to both L. Francis and L. Pessière until ∼6000 cal BP (Fig. 6). Because P. strobus pollen is aerially transported over long distances, some doubt as to its true local abundance might be raised. However, the higher influx of P. strobus pollen during the mid-Holocene (Fig. 6a) is associated with higher values for Taxus, which is locally dispersed (Richard 1976). Moreover, macrofossil analysis carried out for radiocarbon dating of the L. Francis sediments revealed male cone scales and staminate scales of P. strobus spanning 6800–3700 cal BP (Carcaillet et al. 2001). Finally, macro-remains of P. strobus dating back to 7000 cal BP have been excavated from peat near Lake Bissonnette at ∼49°N, located between L. Francis and L. Pessière (Terasmae and Anderson 1970). Hence, there is clear evidence in the form of pollen and macrofossils of a northward expansion of P. strobus during the mid-Holocene, in association with present-day taxa that are characteristic of mixed forests south of 48°N (Taxus canadensis, Thuja occidentalis, B. alleghaniensis).

BF were quite stable in terms of pollen diversity since the afforestation (Fig. 2c–h), especially for L. Francis. However, moderate changes affected the forest several times according to the ROC (Fig. 2a, b). At L. Pessière, where subtle changes in ROC are recorded (Fig. 2b) and the basic forest composition has remained apparently unchanged since the afforestation, we observed the highest magnitude of diversity decrease, suggesting a loss of community structure in favour of the main dominant species (spruce, jack pine, birch). The forest appears to have changed by gradual and slight rearrangement of species within the community. Indeed, no profound and rapid community transformation, from one type to another, was identified over the last 7000 years.

The main overall result of the late-Holocene changes is observed at L. Francis (Figs 2a, 7c). This vegetation transformation, combined with the gradual intracommunity changes at L. Pessière, reveals a convergence of the two sites over the last 1200 years (Figs 5, 6). This convergence suggests a decrease in diversity at the landscape scale resulting from a decline in differences between the sites. Moreover, since the process prompted a greater change in the pollen assemblages at L. Francis than at L. Pessière, there appears to have been a progressive southerly increase in the dominance of boreal species at the expense of species more characteristic of the modern Mixedwood forests.

All regional pollen studies indicate a similar vegetation history since the retreat of the proglacial waters from the areas of the Mixedwood BF (Richard 1980; Liu 1990; present study) and the Coniferous BF (Garralla and Gajewski 1992; Gajewski et al. 1993; present study). Pinus strobus, Thuja occidentalis and Taxus (all studies), and Betula cf. alleghaniensis (present study) abound more around 6000 cal BP in the Mixedwood BF than today. The abundance of these plants characterises the mixed-deciduous forests located at present south of L. Francis around 47°N latitude. P. strobus, Thuja occidentalis, Taxus, and B. cf. alleghaniensis show a long decline until the late Holocene, whereas Picea experiences a progressive increase right up to the present day, starting after afforestation in the Coniferous BF and for the past 1400 years in the Mixedwood BF. Internal processes, e.g. competition, rate of regeneration, pedogenesis and humus accumulation, can be enhanced. The increase in soil acidification and buffering capacity by humus accumulation was proposed as an explanation for long-term lacustrine geochemical processes, as deduced from diatom assemblages (Philibert et al. 2003). This suggests that soil maturation is able to affect the millennial-scale changes of vegetation records. Furthermore, sediment analyses show that paleofires do not remove enough of the humus layer, thus creating a protection for soil that has little chance of being rejuvenated (Carcaillet et al. 2006), which would be unfavourable for ecosystem productivity (Lecomte et al. 2006). However, the broad-scale character of the changes in pollen assemblages supports the hypothesis that climate controlled the long-term plant dynamics. The maximum abundance of mixed-deciduous forest species occurred during the mid-Holocene, which was characterised in eastern Canada by increasing degree-days resulting from higher solar irradiance (Sawada et al. 1999; Kerwin et al. 2004), and a reduction in annual precipitation (Lavoie and Richard 2000; Viau and Gajewski 2009). Subsequently, the late Holocene, with a progressive decrease in seasonality, accompanied by an increase in winter precipitation, particularly over the last 3000 years (Lavoie and Richard 2000; Moos et al. 2009) would have been more favourable for boreal woody taxa. Therefore, it is reasonable to propose long-term climatic change as the main explanation for the observed millennial vegetation trends over the last 7000 years.

The timing of the observed vegetation change (Figs 5, 6) is roughly synchronous with other vegetation changes in north-eastern North America (Anderson et al. 1986; Gajewski 1987; Miller and Livingstone 1993; Marcoux and Richard 1995; Rhodes and Davis 1995; Lavoie and Filion 2001) that occurred between 1800 and 1200 cal BP, depending on the proxies used and the quality of chronologies. In each of these studies, the trend corresponds to a more or less progressive rise in the abundance of boreal species, relative to more thermophilous species. Furthermore, the change in birch pollen diameter, starting 1800 cal BP (Fig. 7), corresponds to the disappearance of Betula alleghaniensis and the relative increase of B. papyrifera or dwarf species such as B. pumila or B. glandulosa. Similar trends in birch-pollen diameter have been observed in Quebec (Richard 1980; Richard et al. 1982) during the late Holocene, and have been interpreted also in terms of the disappearance of B. alleghaniensis (Marcoux and Richard 1995). These Betula records support the hypothesis that broad-scale processes were more responsible for the detected changes than local events. However, the timing is not precise, for a few hundred years separate the different data in the studies, and no geographical pattern can be identified. Simple problems associated with dating (dating numbers, type of dated material, age–depth models, sediment disturbances, etc.) or the lack of availability of such Betula-diameter reconstructions may explain this coarse time correlation with vegetation change over the last 1500 years. Because a north-eastern American trend of vegetation ‘borealisation’ (i.e. increases in boreal taxa) appears in many studies (Richard 1994), it is difficult to support a local explanation for the recent expansion of boreal trees like Picea mariana, Pinus banksiana, Abies, Betula papyrifera and Alnus viridis ssp. crispa. Broad-scale spatial processes, such as those triggered by reduced seasonality in response to orbital forcing, eventually combined with local-scale processes (changes in the disturbances regime, e.g. fire, insect outbreaks) appear a likely explanation.

Changes in fire frequency have not triggered major transformation of the eastern North American BF. An important resilience is observed after changes in MFI. This study underlines the fact that although fire is a major property of BF ecosystems, susceptible thresholds were never reached to promote new communities. However, different pollen trajectories are recognised in the Mixedwood and the Coniferous BF, the causes of which date back to afforestation. The pollen composition has been significantly converging for the last 1200 years owing to an increase in boreal species at the southern site (Mixedwood BF). This process started at least 2000 years ago, when MFI decreased in the southern site. Changes in MFI alone do not explain the historical and current distribution of vegetation but they may have been partly responsible for accelerating the climatic process of ‘borealisation’, likely resulting from orbital forcing.