Three large fire years threaten resilience of closed crown black spruce forests in eastern Canada

An emerging paradigm regarding vegetation response to climate warming is that the interaction of weather extremes and disturbance will trigger abrupt changes in ecosystem types by overcoming resilience of dominant species. Black spruce (Picea mariana (Mill.)) ecosystems are widespread across the North American boreal forest, because of ecophysiological adaptations that allowed these communities to thrive in fire-prone areas. We investigated resilience of spruce-moss forests to weather-disturbance interaction after a 3-year period (2005 to 2007) of major fire activity caused by extreme fire weather in eastern Canada. Pre- and post-fire conifer densities and environmental parameters related to seed rain, post-fire seedbeds, microclimate, and post-fire weather were measured in 133 burned stands throughout the closed-crown forest of Quebec. Critically low black spruce (BS) regeneration was observed in almost all of the stands, leading to a decrease in stand density and a shift of species dominance from BS to jack pine (Pinus banksiana (Lamb.), JP). The studied sites were characterized by thick residual organic matter, resulting in a predominance of charred duff, a seedbed associated with low water retention and high variation in temperature. While high levels of JP seedling establishment were reported on this seedbed, it was unfavorable to BS germination and survival in the context of warm and dry weather that prevailed in post-fire summers. In these ecosystems, early vegetation establishment patterns are generally reliable predictors of future stand composition and the exclusion of BS will presumably be maintained through succession. During large fire years, high proportions of the landscape are subjected to the interaction of fire regime and weather that creates unsuitable conditions for BS regeneration. Hence, vegetation change is susceptible to happen at a broad scale. Therefore, the frequency of major fire years could have a decisive influence on the rate of vegetation response to climate change in this biome.


INTRODUCTION
Vegetation response to climate change has provoked numerous hypotheses and requires close monitoring in the next decades (Johnstone et al. 2010b, Walther 2010. One hypothesis suggests that vegetation will respond gradually to climate change and that species will migrate northward as temperature increases (Davis and Botkin 1985, Overpeck et al. 1991, Parmesan and Yohe 2003. Palynological data provide evidence to support this paradigm for North American vegetation dynamics during the last 18,000 years (Prentice and Bartlein 1991). The current distribution of vegetation in boreal biomes is also fundamentally tied to climate (Woodward and Williams 1987). However, at a time scale relevant to policy and forest management (decades to millennia, Chapin et al. 2004), and within species ranges, these ecosystems have shown to be resilient to climate change , Pollock and Payette 2010, Lloret et al. 2012). This resilience is attributable to species longevity, ability to reproduce vegetatively and biological legacies that reduce opportunities for regeneration of alternate species , Barrett et al. 2011). An emerging hypothesis is that climate extremes could provoke faster and more ample ecosystem responses than mean trends (Gutschick andBassiriRad 2003, Jentsch et al. 2007). Extremes in climate may alter disturbance frequency, size and severity which could cause abrupt changes in vegetation if resilience thresholds of the dominant species are exceeded (Scheffer andCarpenter 2003, Chapin et al. 2004). In the boreal forest of North America, fire activity fluctuates greatly from year to year, with episodic major fire years characterized by a long fire season and numerous large fires (Stocks et al. 2003, Kasischke and Turetsky 2006, Kasischke et al. 2010. Large fire years are related to extreme fire weather conditions, more specifically to extended periods with low daily precipitation, low relative humidity and high temperatures (Flannigan and Harrington 1988). When these warm and dry conditions persist, they can also affect post-disturbance regeneration. The few-year time window during which dominant species establish after fire is a critical phase for future stand development (Sirois and Payette 1989). Even if mature trees may cope with environmental stress because of an established network of resource capture, the cohort that is germinating in the post-disturbance phase may be more vulnerable Bliss 1980, Johnstone et al. 2010b). Interactions of climate and disturbance are thus likely to provoke the most important and long lasting effects on ecosystem structure and function Beierkuhnlein 2003, Chapin et al. 2004). In the region of the forest-tundra border in eastern Canada, tree depletion in spruce stands after fire lead to expansion of tundra inside of the northerly boreal forest (Sirois and Payette 1991) or to major shifts in species dominance from black spruce (Picea mariana (Mill.), abbreviated BS) to jack pine (Pinus banksiana (Lamb.), abbreviated JP) (Lavoie and Sirois 1998). Forests in this ecotone of high latitudes are particularly responsive to climate-disturbance interaction, because of the decreased reproductive capacity of trees exposed to cold temperatures (Sirois andPayette 1991, Sirois 2000). Further south, in the commercial forest, the closed spruce-moss stands are expected to be more resilient (Sirois and Payette 1991) with fire initiating direct replacement by the same species assemblage Bliss 1978, Viereck 1983). However, Girard et al. (2008Girard et al. ( , 2009 revealed that, over the last 50 years, 9% of the area occupied by closed spruce-moss forests in Quebec has converted to open sprucelichen woodlands after low severity fires or repeated disturbances. If fire can lead to important declines in spruce density within the BS feathermoss bioclimatic domain, then shifts in species dominance from BS to JP may occur in sites where spruce fails to regenerate. JP, like BS, possesses serotinous cones that allow the species to thrive in fire prone areas (Cayford and McRae 1983). It is a shade-intolerant moderately longlived species that reaches an optimum in seed production earlier than BS (Cayford and McRae 1983). JP is better adapted to areas where the fire cycle is short and pine stands are eventually replaced in the extended absence of fire (Cayford andMcRae 1983, Le Goff andSirois 2004) while BS is a shade tolerant long-lived species that can reproduce vegetatively by layering 50 to 70 years after stand initiation Bliss 1980, Viereck 1983). Fire return interval is thought to play a major role in the respective spatial abundance of the two species in the landscape (Le Goff and Sirois 2004). Lavoie and Sirois (1998) hypothesized that fire severity and postfire weather may be two additional drivers of the shift in species dominance, while Johnstone et al. (2010a) found that post-fire establishment of BS was higher in moist sites characterized by low fire severity in Alaska.
Large fire years represent timely opportunities to study species resilience to climate-disturbance interactions: they correspond to anomalies in both climate and disturbance regimes, and the numerous large fires provide potential study sites across a wide range of landscape and fire conditions (Johnstone et al. 2010a). In the province of Quebec, the years 2005 to 2007 were characterized by major fire activity. Compared to the 1981-2010 period, spring and summer weather of these three years was warmer than normal, with locally important deficits of June to August precipitation (Environment Canada 2010). Large fires occurred early in the season: 1.2 M ha burned over these three years, which represents 12% of the total area burned over the past 36 years (Ministè re des Ressources Naturelles du Québec, MRN 2010).
The objective of this research was to investigate conifer regeneration after this 3-year episode in the commercial forest of Quebec, assuming that tree species establishment in this early successional phase would indicate both future stand composition and potential loss of resilience. We hypothesized that weather-disturbance interaction could cause exclusion of BS from post-fire stands and that JP could replace BS as the dominant species in these stands. We also aimed to identify the main environmental drivers of regeneration success for each species, in order to understand the underlying mechanisms of the transition from BS to JP and to determine where in the landscape a loss of resilience was more likely to occur.

Study area
The study area is the BS-feathermoss bioclimatic domain of Quebec. It occurs from 488 N to 528 N and covers 412,400 km 2 (Fig. 1). In the western sector of the domain, in the James Bay region, mean annual precipitation is 700 mm, mean annual temperature is À0.658C with an average growing season of 165 days. In the North Shore region in the eastern sector of the domain, annual precipitation is higher (1000-2000 mm), due to maritime influence of the Atlantic Ocean and the Gulf of St Lawrence. Mean annual temperature is À1.58C, and the growing season is approximately 150 days (Saucier et al. 2009). The southwestern sector of the domain is a plain slightly sloped towards James Bay. The central area is a landscape of moderate-altitude hills (400-600 m) with a few higher massifs (1000 m). In the eastern sector of the domain, the landscape is uneven and mountainous; the altitude ranges between 550 and 900 m. The western edge of the study area is within the Little Clay Belt, characterized by clayey and loamy superficial deposits. The remainder of the study area is located on the Precambrian Shield, where superficial deposits are mostly glacial tills. All of the study sites are located outside the permafrost area. Most of the soils are ferro-humic podzols with a thick organic layer (5 to 40 cm; Robitaille andSaucier 1998, Soil Classification Working Group 1998). BS is the dominant species of this domain. It occurs mainly in mono-specific stands, where tree cover generally exceeds 40%. JP is the second species in importance; it is more frequent in the western and central portions of the domain, where it forms mixed stands with BS on well drained sites with coarse textured soils (Saucier et al. 2009). The fire cycle varies from 100 years north of central Quebec, to more than 1000 years in the southeastern part of the domain (Chabot et al. 2009).

Sampling design
A systematic nested design was used (Fig. 2). The MRN fire database was used to identify sites burned by 2005 to 2007 wildfires that were less than one km away from a passable road. We first visited these fire sites in 2009 and selected 14 fire sites with large patches of forest that had not been affected by salvage logging. We made preliminary assessments of burn depths, damage to the canopy and conifer regeneration. Preliminary data exploration indicated that both burn severity and proportion of JP in pre-fire canopy varied across topography and drainage gradients. Therefore, we established transects on toposequences, in order to capture a large variability of slope positions and drainage with contrasting burn severity and stand composition. Within each burn, two to four transects were set up; the number of transects depended on the availability of mature unsalvaged stands. Each transect consisted of several 20 3 20 m plots. Plot topographic positions were classified as crest, upslope, midslope, depression, downslope or toe (Fig. 2). A total of 14 fire sites, 46 slopes and 133 plots were inventoried throughout the BS-feathermoss bioclimatic domain (Table 1). Data were collected in 2010 between June 15th and August 15th. Each plot was systematically cross-ruled in order to establish nine 1.5 m-diameter circular subplots at the four corners and inside the plot. Then, an inner square of 100 m 2 was established at the center of the plot (Fig. 2). Measurements made in the inner square or in the nine subplots were averaged at the plot scale before further analyses.

Measurements
All the variable abbreviations and definitions are presented in Table 2.  Table 1.  Pre-and post-fire conifer densities.-In every inner 100 m 2 quadrat, total tree density of the pre-fire stand (treesÁm À2 ) and species-specific density (JP trees or BS treesÁm À2 ) were determined. Trees larger than 2 cm in diameter and taller than 2 m in height were included in this census, even fallen trees if they were rooted in the quadrat before they fell. However, trees for which charring and decay state indicated prefire death were not considered. BS and JP seedlings present in the inner 100 m 2 quadrat were also counted, in order to compare pre-and post-fire species densities. To evaluate possible changes in species dominance after fire, we used compositional difference (CD) as our response variable (modified from Ilisson and Chen 2009). Because post-fire stands were younger than 25 years, CD was calculated as species density of the post-disturbance stand (i.e., seedling density) subtracted from that of the pre-disturbance stand (i.e., burned adult-tree density) divided by total density of pre-fire stand ( Table 2). Variations of CD were related to environmental parameters in the light of three a priori hypotheses: ''Regeneration success is mainly driven by: post-fire seed rain (hypothesis #1) post-fire seedbed characteristics (hypothesis #2) local climate and post-fire weather (hypothesis #3).'' Parameters related to seed rain.-By ''seed rain'' we refer to the quantity of viable seeds dispersed on post-fire sites. The size of the aerial seedbank is correlated with size class distribution in mature stands of BS and JP Johnson 1994, Johnstone et al. 2009). Stand age was estimated to ensure that trees were mature at the time of fire event. On each transect, nine dominant burned trees were chosen and cored at breast height (1.3 m) to obtain 5.4 mm diameter cores for stand-age assessment. Cores were mounted in wooden blocks and sanded with 220 and 320-grit papers prior to counting rings, using a microscope with a 403 objective lens. We measured diameter at breast height (1.3 m, Dbh in cm) of all trees inventoried in the inner 100 m 2 square, and species-specific diameter was used to represent variations in seed production.
Post-fire seed availability of serotinous conifers can also be affected by fire behavior, because crown fires that prevail in these ecosystems can be sufficiently intense to kill seeds by heating or to cause cone ignition and destruction, affecting both seed quantity and seed viability (Kasischke et al. 2008, Johnstone et al. 2009). An aerial severity index was determined on the basis of three criteria: scorch height, presence and abundance of secondary twigs and bark detachment. Ratings were attributed according to the magnitude of damage for each of the three criteria and were then added to obtain the overall aerial severity index for the plot. Distance to the closest unburned patch was also included in models to account for seed dispersal from remnant living trees in the neighborhood of our plots Johnson 2000, Johnstone et al. 2009). Geographical coordinates of the closest patch of living trees were recorded with a GPS and distance between this point and the plot was measured with GIS software (ArcGIS 10, ESRI 2011, Redlands, California, USA).
Burn severity and post-fire seedbeds.-Soil borings were carried out with an auger to determine site drainage class and to take samples of the deepest mineral horizon (C horizon) of each plot. Theses samples were transported to the lab in a cooler and soil texture was assessed by decantation (McKeague 1978).
Ground cover type was noted in each of the nine 1.5 m-diameter subplots. Cover type classes included: unburned organic matter, charred organic matter, exposed mineral soil, woody debris and rocks. The proportion of different seedbeds was then calculated for each plot. Conifer seedlings were also identified and counted in the nine subplots in order to assess species-specific establishment on different seedbeds. Residual organic layer thickness was measured in each of the nine 1.5 m-diameter circular subplots. Measurements were averaged in order to obtain mean residual organic matter thickness for the plot (ROM in cm). Pre-fire organic layer depth and consumed organic matter thickness were reconstructed. Adventitious roots grow above the initial root collar of BS trees into the organic soil/moss layer that surrounds the tree (LeBarron 1945). The distribution of these roots is clearly visible on burned trees for many years after fire. Therefore, height measurements of small (,2 mm diameter) adventitious roots on BS boles above the remainv www.esajournals.org ing soil organic material can be used to reconstruct pre-fire organic matter depth (Boby et al. 2010). We measured the distance separating the highest adventitious root and surface of mineral soil in centimeter (H i ) in 5 subplots out of 9 (not all of the nine subplots included BS boles). We used the equation of Turetsky et al. (2011) to estimate pre-fire organic matter depths (linear regression model: F 1, 317 ¼ 3378; R 2 ¼ 0.93; p ¼ 0.0001; slope ¼ 1.03 6 0.02; intercept ¼ 4.56 6 0.35). Consumed organic matter thickness at the plot scale (COM, in cm) was then calculated as the difference between the estimated pre-fire depth and ROM.
Local climatic conditions and post-fire meteorological data.-We used the BioSIM software (Régniè re 1996) with daily precipitation data from the provincial climate data base (Ministè re du Développement Durable, de l'Environnement, de la Faune et des Parcs du Québec, MDDEFP, 2012) to generate monthly total precipitation for years 2005 to 2009 in our sites. The BioSIM software estimates meteorological data for each plot by interpolation and extrapolation from weather records of the height nearest meteorological stations to each plot location. We first compared rainfall in 2005, 2006 and 2007 to the normal over the 1981-2010 period. Then, total amount of precipitation (TsP), for summer months (June, July and August) following fire events, and for summer months of the two following years, were included as explanatory variables in regeneration models. On each transect, slope inclination and aspect were determined with a compass and a clinometer. Slope position was recorded as a 10-level categorical variable: flat sites with inclination ,2% were classified as ''crest or shoulder'' or ''toe or depressions'' and sloped sites with inclination .2% were sorted according to slope aspect (''N'', ''NE'', ''E'', etc.). The greatest number of plots was observed in the West-facing slope category; therefore ''W'' was used in analyses as the reference level for comparisons among slope positions. We used inclination, aspect and latitude data to calculate the ''heat load index'' of a site, based on equation 2 in McCune and Keon (2002).

Data analysis
Comparison of BS and JP regeneration.-We first compared the regeneration of the two species using a linear mixed model: CD was the response variable, species was the only fixed effect, and fire sites, toposequences and plots were random effects. The nlme package was used (3.1-101 version, Pinheiro et al. 2012) in R statistical software (2.15.0 version, R Development Core Team 2012, Vienna, Austria).
Identification of environmental factors explaining individual species regeneration.-The regeneration of each species was then analysed separately. Model selection was run to identify the combination of environmental variables that explained the greatest amount of variation of CD. For each species, 12 linear mixed models with the same Included with a linear term and a quadratic term in candidate models for BS and JP regeneration. à Included with a linear term and a quadratic term in candidate models for BS regeneration.
v www.esajournals.org random effects (fires and toposequences) were defined (Appendix: Table A1). The fixed effects were different combinations of variables related to the three hypotheses described above (influence of seed rain, seedbed characteristics and local climate). The sets of candidate models comprised 6 simple models related to a single hypothesis and 6 comprehensive models which were combinations of all the simple models. An interaction term between the variable related to seedbed quality and the variable related to climate was added to every comprehensive model since these factors influence the amount of water available to seeds and seedlings (Duchesne and Sirois 1995).
The numerical explanatory variables were centered prior to analysis. Two variables were included with a quadratic term. First, the quadratic effect of aerial severity was added in both species candidate models, as medium severity fires are expected to result in greater seed rain, because of possible seed damage during high-severity fires and granivory after low-severity fires (Ahlgren 1966, Greene et al. 2004. ROM was also included with a quadratic term in BS models to account for the fact that both very thin (mineral soil with thin humus) and very thick (unburned Sphagnum) ROM are expected to result in high BS establishment rates (Greene et al. 2004, Veilleux-Nolin andPayette 2012). For JP, only the linear term of ROM was included in candidate models. We performed model selection using the package AICcmodavg in R (version 1.24, Mazerolle 2012) and Akaike's information criterion (AICc) and Akaike weights (w) were used to identify the most parsimonious model (Burnham andAnderson 2004, Mazerolle 2006) (Appendix: Tables A2 and A3). A weighted average of the estimates, incorporating model uncertainty, was calculated. We checked model fit of the most parsimonious comprehensive model and calculated the predicted values, standard errors and confidence intervals presented in figures using Monte-Carlo simulations (Gelman and Hill 2007).
Variations of burn severity across the landscape.-We analysed how the proportion of JP in pre-fire stand, slope position and heat load influenced burn severity using linear mixed models: the dependant variables were pre-fire organic layer, percentage of organic layer consumed by fire, and ROM. The effect of pre-fire organic layer thickness on consumed organic layer was also investigated. Again, fire sites and toposequences were included as random effects.

Comparison of regeneration of BS and JP
Pre-fire stand characteristics.-Prior to fire events, 10% of the studied burned plots were dominated by JP. Pre-fire density of BS averaged 0.28 6 0.16 stemÁm À2 (mean 6 SD), whereas prefire density of JP averaged 0.04 6 0.07 stemÁm À2 . Proportion of BS in pre-fire canopy averaged 84 6 22%, whereas JP proportion averaged 14 6 22%.
Species had a significant effect on variations of CD, with JP having a significantly higher regeneration than BS. Mean CD for BS was negative (À0.40 6 1.76), with high inter-plot variability. CD values were negative in 84% of the plots. When CD values were averaged at fire site scale, 13 burns out of 14 displayed negative CD values (Fig. 1). Mean CD value was positive for JP (þ0.69 6 1.421) with smaller variability than for BS. JP CD was positive (!0) in 90% of the plots. CD values were also averaged at fire scale and were strictly positive in 13 fire sites and null in one fire site (Fig. 1).

Principal environmental factors influencing conifer regeneration
Model selection.-When BS regeneration was analysed, the best model given the set of 12 candidates was a comprehensive one combining hypothesis #1, #2 and #3, with an Akaike weight of 0.84 (Appendix : Table A2). This model included Dbh of burned BS trees, ROM, TsP, and the interaction between the last two parameters. For JP regeneration, the only model that had a Di value ,2 was the simple model with slope position as the unique fixed effect (Appendix: Table A3). This model, corresponding to v www.esajournals.org hypothesis #3, had an Akaike weight of 0.68. In both analyses, the most parsimonious models had high rankings, but rather than basing conclusions solely on interpretation of these best models, we used multi-model averaging, as it is recommended when w i 0.90 (Burnham and Anderson 2004). Model-averaged estimates (MAE) and 95% confidence intervals are presented in Appendix: Table A4. Parameters for which confidence interval excludes 0 are considered to have non-null effect on variations of species CD.
Parameters related to seed rain.-Mean age estimations for pre-fire stands ranged from 60 years to 181 years (Table 1) indicating that canopy trees had reached sexual maturity at the time of fire event. Model averaging revealed that both BS and JP regeneration increased with increasing Dbh of pre-fire trees. The other parameters related to hypothesis #1, aerial severity and distance to the closest unburned trees, did not influence spruce nor pine regeneration since confidence intervals around modelaveraged estimates of these parameters included zero (Appendix : Table A4).
Burn severity.-In the studied plots, fire consumed on average 50% of the organic matter (between 8 and 82%, 614%) and the resulting ROM averaged 17 6 8 cm. Both pre-and post-fire organic layers were thicker on toe or depression positions and decreased with increasing proportion of JP in pre-fire canopy ( Table 3). Percentage of organic layer consumed was significantly higher in pre-fire plots with shallow organic layers, high proportions of pre-fire JP and on SWfacing slopes (Table 3).
Post-fire seedbeds.-Charred organic matter was the most frequent seedbed, covering on average 86% of the plot area (16% of the plot area was charred organic matter covered by Polytrichum and 70% of the area was charred organic matter with litter). Unburned organic matter was almost exclusively composed of patches of unburned Sphagnum. Sphagnum, mineral soil and woody debris and rocks, accounted for respectively 8%, 3% and 3% of post-fire seedbeds (Fig. 3). JP regeneration decreased with increasing ROM (Appendix : Table A4). Establishment rates of JP on the different seedbeds were consistent with relative seedbed frequencies, as 59% of the seedlings were found on the most abundant seedbed: bare charred organic matter. However, establishment rates on mineral soil and Sphagnum were high (respectively 9 and 16%) compared to their relative importance in the studied plots (Fig. 3). Sphagnum, despite the fact that it only accounted for a small proportion of total seedbeds, represented the preferential substrate for BS regeneration: 47% of all the recorded seedlings for this species had established on Sphagnum (Fig. 3). In contrast, charred organic matter without Polytrichum supported less than a third of the total BS seedlings, despite the fact that it was the dominant seedbed (Fig. 3).
Local The interaction between ROM and TsP influenced BS regeneration levels (Appendix: Table  A4). Rainfall of the three post-fire summers had no effect on spruce regeneration in sites where ROM was less thick than 25 cm (Fig. 4). In sites where ROM was very thick (.25 cm), BS regeneration increased with increasing amount of rain received the first three summers after fire events. Slope position was the only variable included in the best model for JP regeneration. In the studied pre-fire sites, the highest JP density was observed on N-facing slopes. After fire, JP density increased in all site categories, and was still maximal on N-facing slopes (Table 4). Apart from this peak of regeneration on northern expositions, regeneration on sloped sites showed a gradual response to Aspect, with higher CD values on S and SW-facing slopes (Table 4). When model-averaged estimates were computed to compare each slope position with W-facing v www.esajournals.org sites, only N-facing sites and toes and depressions had CI excluding zero. Following fire, BS density decreased in all site categories except for S-facing slopes, toes and depressions, where seedling densities were respectively higher and about the same as adult densities (Table 4). Slope position was not included in the best model for spruce regeneration, but CD values were significantly higher on toes and depressions.

Potential loss of resilience of BS stands
Fires during the years 2005-2007 were associated with an important change in stand structure, since a decrease in density of the dominant tree species was observed in 84% of the plots. Chrosciewicz (1976) suggested that a uniformly-spaced stocking of 7.4 stemsÁm À2 is required for development of commercial BS stands; this threshold was not reached in any of our plots. Seedling densities were much below those noted by Johnstone et al. (2004) (0.5-0.6 stemÁm À2 with up to 8 stemsÁm À2 ) and by St-Pierre et al. (1992) (0.83, 0.79 and 0.35 stemÁm À2 ). The decrease in stand density was associated with a change in species dominance, with 90% of pre-fire plots dominated by BS and more than 60% of plots dominated by JP after fire. Pine seedling density was similar to that measured by St-Pierre et al. (1992) (1.5; 1.6 and 0.22 stemÁm À2 ) and lower but in the range of values predicted by the equation of Greene and Johnson (1999) relating regeneration density to pre-fire basal area (0 to 3.1 stemsÁm À2 ; mean of 0.56 6 0.74). The high JP seedling density, considering its low proportion in pre-fire canopy, was not sufficient to offset the density decrease associated with low spruce regeneration. Chrosciewicz (1970) suggested that 0.2 JP seedlingÁm À2 was required for adequate commercial development of JP stands; a threshold met in only 1/3 of the post-fire plots. Hence, the ecosystem regenerating after fire in our study sites was dissimilar to the original closed-crown BS-dominated stand, and was characterized by low stem density and JP dominance. This establishment pattern resembles the transition observed in the taiga by Lavoie and Sirois (1998). We focused on the early stages of re-vegetation, since previous studies showed that most of BS and JP seedlings establish within three post-fire growing seasons (Black and Bliss 1980, Lavoie and Sirois 1998). Following fire in mixed spruce-pine stands, there may be an initial successional phase of JP dominance, followed by mortality of pine and return to spruce dominance (Carleton 1982, Lesieur et al. 2002, Lecomte 2005. BS return is generally associated with divergent growth rates of BS and JP following concomitant establishment of the two species (Van Cleve and Viereck 1981, Carleton 1982, however in the present study we did not observe this establishment of spruce. Seed could be dispersed from proximate nonserotinous BS (Carleton 1982) but several feedbacks are likely to preclude further spruce establishment. Low stand density and JP dominance usually favour the formation of a lichen mat (Payette et al. 2000) and a shrub layer of Kalmia angustifolia, both of which can prevent seedling recruitment (Fisher 1979, Mallik 1987, Jasinski and Payette 2005. Thus, our observations of the initial regeneration patterns suggest that the new ecosystem is not a transient successional stage and that the 2005-2007 fire events caused a loss of resilience of the studied spruce stands. However, longer term studies of species dynamics are still necessary to confirm that the new open JP stand is a true alternative state, i.e., an ecosystem that persists under the Notes: * P , 0.05, ** P , 0.01, *** P , 0.001. Estimates are presented when the effect is significant. Sites on W-facing slopes are the reference levels for comparisons between different slope positions. v www.esajournals.org same environmental conditions as the pre-fire BS-dominated stand (Jasinski and Payette 2005).

Seed rains of the two species
For obligate seeders such as JP and BS, seed availability can be a strong limit to post-fire establishment (Greene and Johnson 1999, Jasinski and Payette 2005, Johnstone et al. 2009). Previous studies demonstrated that fire return interval was a good predictor of post-fire relative abundance of the two conifers because of divergent dynamics of seed production between species (Le Goff and Sirois 2004). Regeneration inventories were carried out 3 to 5 years after fire, when the great majority of seeds produced by on-site burned trees had been dispersed. Indeed, post-fire seed release from BS and JP has been shown to be brief, with a pulse in the first and second growing seasons after fire (Charron andGreene 2002, St-Pierre et al. 1992).
Regeneration of both species increased with Dbh of pre-fire mature trees, a surrogate for seed production (Greene and Johnson 1999). Dbh was included in the best model for BS but not for JP. In the range of age of our stands when they burned (60-181 years), seed production of BS could be more variable than that of JP, because sexual maturity occurs later for BS (Cayford and McRae 1983). Neither aerial severity, nor distance to the closest unburned tree patch affected regeneration of the two species. These variables accounted respectively for seed destruction and seed inputs from living seed sources. Studies of post-fire white spruce regeneration have shown that continuous long distance dispersal from external unburned plots could constitute nonnegligible background inputs for stand regeneration (Johnstone et al. 2010b). BS and JP both have contagious dispersal modes (Filion andMorin 1996, Greene andJohnson 2000). If living trees are at distances greater than 200-300 meters from a burn, as was the case for most of our plots, seed inputs from external sources may have been too low to significantly impact postfire regeneration. Observation of this phenomenon will require longer post-fire monitoring (10 years) than in the present study. Finally, we hypothesize that the effect of seed destruction and seed inputs from external green trees on regeneration could be masked by other strong environmental constraints on seedling establishment, such as unfavourable seedbeds and dry weather.

Fire regime and weather in 2005, 2006 and 2007
Lack of good seedbeds.-Thick charred organic matter dominated all post-fire plots. Black colour v www.esajournals.org and low density of this substrate cause high temperatures and low water retention (Miyanishi and Johnson 2002). This negatively affected regeneration of both species, with a less important effect on JP than on BS. Indeed, 73% of pine seedlings occurred on charred duff. In contrast, spruce seedling establishment was low on charred substrate thicker than a few cm, which is consistent with previous observations (Chrosciewicz 1976, Charron andGreene 2002, Veil-Fig. 4. Differential effect of TsP (total amount of precipitation three summers after fire event) on black spruce regeneration in sites of varying ROM (residual organic matter thickness). Symbols represent real data, solid lines represent model predictions, and dashed lines represent 95% confidence intervals. Table 4. Pre-and post-fire characteristics of plots (mean with SD in parentheses) according to topographic position. Mean values presented here do not include four plots that were outliers for spruce regeneration analysis (one W-facing slope plot was discarded); or for pine regeneration analysis (two NW-facing slopes and one SW-facing slope were discarded).
v www.esajournals.org leux-Nolin and Payette 2012). However, Johnstone and  and Shenoy et al. (2011) found that spruce regeneration increased with ROM in Alaska.
JP seedlings have a faster and longer radicle elongation, improving critical access to water resources located deeper in the underlying mineral soil (Thomas and Wein 1985). The inability of spruce to establish on the dominant seedbed likely caused its exclusion in most of the sites. BS post-fire regeneration is generally high on mineral soil, a substrate almost absent in our 14 burns (Charron andGreene 2002, Greene et al. 2004). Sphagnum has also been identified as a good substrate for germination and survival due to its high water holding capacity (Greene et al. 2004, Veilleux-Nolin andPayette 2012). High densities of BS seedlings were indeed observed on patches of unburned Sphagnum, but this will probably cause self-thinning mortality. Moreover, Sphagnum is a poor substrate for future development due to its low nutritional quality (Oechel andVan Cleve 1986, Munson andTimmer 1989).
Low burn severity.-The study sites were characterized by diverse topographic, edaphic and vegetation conditions but all shared the common dominant characteristic of thick ROM (mean of 16.8 6 8.2 cm). ROM recorded in the present study was in the same range of what was measured by Veilleux-Nolin and Payette (2012) in Quebec (between 10.1 and 23.2 cm), but exceeds mean ROM measured in Alaska by Johnstone and Kasischke (2005) and Boby et al. (2010) (respectively 6.6 cm and 8.2 6 1.0 cm). We hypothesize that the observed low burn severity was related to the important thickness of the prefire organic layers. Reconstructed pre-fire organic layers of the studied plots were deeper than those calculated for sites in Alaska by Turetsky et al. (2011) (32.4 6 8.5 cm vs. 19.5 6 7.13 cm). This may reflect the longer fire-return intervals in Quebec (Chabot et al. 2009, Barrett et al. 2011). Depths of burn (14.6 6 4.8 cm) were in the same range of what was observed in Alaska by Kasischke and Johnstone (2005) (mean of 23.1 cm) and by Turetsky et al. (2011) for early and late season burns (respective mean COM of 11.9 and 21.7 cm). However, on average, 50% of total organic matter was burned, resulting in thick ROM and a quasi-absence of exposed mineral soil. An important factor limiting burn severity is ground fuel moisture, which increases with depth of pre-fire of organic layers that have an insulating effect (Kane et al. 2007). Percent organic layer consumed was indeed significantly affected by pre-fire organic matter depth (Table  4).
Ground fuel moisture and organic matter accumulation vary across the landscape with lower moisture levels on S-facing slopes and flat uplands that receive more radiation than Nfacing slopes and flat lowlands (Kane et al. 2007). The percentage of organic matter removed was higher on SW-facing slopes, where ground fuels were likely drier at the time of fire. Moreover, ROM was higher in cold and moist sites on toes of slopes and depressions. Variation in insolation and ROM among different topographic positions resulted in more important JP regeneration rates on S and SW-facing slopes than on NE slopes and toes or depressions, because shallower ROM on the former resulted in easier access to mineral soil water. Slope position was not included in the best model for BS, as BS regeneration rates were uniformly low across the landscape, except in toes and depressions where unburned Sphagnum was abundant (Appendix : Table A4).
Pre-burn organic layers were thinner and more severely burned in sites where pre-fire proportions of JP were high, which is consistent with results from Duchesne and Sirois (1995). JP and BS can be found on a variety of site conditions in boreal forest landscapes, but JP stands usually dominate in the driest uplands characterized by low fertility, coarse-textured soils (Cayford and McRae 1983), while moderately-to poorlydrained cold sites are usually occupied by pure BS stands (Viereck 1983). Biomass production and moss accumulation are less important in JP stands than in BS stands (Gower et al. 1997, Rapalee et al. 1998. Such species-associated differences in site quality and ecosystem functioning can result in a difference in burn severity and post-burn organic layer thickness, as observed in the present study. The 10 plots sampled on N-facing slopes were characterized by a high percentage of JP in the pre-fire stands, and consequently displayed moderate ROM thickness in spite of reduced sun exposure. Regeneration of JP peaked on these northerly slopes where accessible soil water resources and protection from solar exposure resulted in optimal v www.esajournals.org protection from drought stress. Finally, ground fuel moisture evolves during the fire season with maximum humidity levels shortly after snowmelt that decrease as fuel progressively dry during spring and summer (Kane et al. 2007, Viereck 1983. All sites were burned by large early season fires which were abundant in years 2005 to 2007 (88% of area were burned before July 31st, during these 3 years, vs. 50% in the 1973-2009 period). This early season timing may have been an additional cause for low severity in our sites. However, in Quebec, Veilleux-Nolin and Payette (2012) found that variations in ROM were unrelated to fire season, with low severity recorded in the four summer burns as well as in the nine spring burns studied. In Alaska, Turetsky et al. (2011) showed that the relationship between fire season and burn severity did not hold during large fire years, as deep burning occurred through the entire season.
Low precipitation in the first growing seasons after fire In , 2006, several summer months were marked by low precipitation in our sites (BioSIM, Environment Canada 2010, MDDEFP 2012. TsP only influenced BS regeneration in plots where ROM was very thick, i.e., in plots with high Sphagnum cover, suggesting a hierarchy in environmental drivers of BS regeneration. Seedbed quality was the strongest constraint; on plots where only charred organic matter was available, regeneration of BS was too low to observe any effect of rainfall, which was only detectable once the seedbed constraint was relaxed on Sphagnum microsites.

Specificity of large fire years
Results of the present study highlight the vulnerability of a widespread fire-adapted species to a combination of low burn severity and dry weather in sites burned during large fire years. Large fire years are timely opportunities to study the initial regeneration phase in diverse environmental conditions as the greater extent of area burned provides numerous potential study sites. Finding the same range of conditions would be more challenging after small fire years, nevertheless, it would be interesting to compare regeneration patterns of the studied stands with observations from sites burned during interme-diate or small fire years. Several previous studies provide empirical data on BS establishment after wildfires in Quebec. For most of these, sampling sites were located within 1 to 4 fire sites that burned during large fire years (St-Pierre et al. 1992, Greene et al. 2004, Jayen et al. 2006). Veilleux-Nolin and Payette (2012) inventoried 13 fire-sites, 10 of which were burned during large fire years (2005,1991) and the three remaining were burned during a small fire year (1999). In two of the 1999 sites, seedling density exceeded pre-fire BS density, while in the third, establishment rates were very low. We hypothesize that abiotic stress and environmental conditions unsuitable to spruce regeneration are not exclusive to large fire years and could occur punctually during ''normal'' fire seasons. However, large fire years such as [2005][2006][2007] in Quebec are characterized by extended drought periods, which increase the likelihood of unfavourable post-fire conditions and the probability of regeneration failure.

Future vegetation changes in the spruce-moss bioclimatic domain
Our findings support hypotheses suggesting that extreme events and climate variability could affect ecosystems to a greater degree than mean trends in climate (Gutschick andBassiriRad 2003, Jentsch et al. 2007), and that gradual species responses to directional climate warming are likely to be overshadowed by abrupt transitions triggered by weather-disturbance interactions , Johnstone et al. 2010b). If the response of the closed-crown spruce forest of Quebec to climate warning was gradual, the spruce-moss forest should theoretically expand progressively northward; Girard et al. (2008) demonstrated that the opposite is happening. In the present study, we inventoried 14 burns across a 600 km-west to east transect, covering a wide range of site conditions and variation in regional climate. We assumed that observations in these sites represented a reliable overview of the vegetation response to fires that occurred throughout Quebec during these three years, burning 1.2 M ha of forest across the province. Large proportions of forest burned during large fire years increase the probability to observe broad-scale loss of resilience compared to years when ecosystem change can be distributed across v www.esajournals.org the landscape and limited to the patch scale. Therefore, the frequency of major fire years could have a decisive influence on the rate of vegetation response to climate change in the closedcrown spruce forest.

SUPPLEMENTAL MATERIAL APPENDIX
Notes: Estimates are weighted according to model ranking presented in Tables A2 and A3. Estimates are in boldface when 95% CI excludes zero, meaning there is strong evidence of an effect of the parameter on species regeneration. v www.esajournals.org