Resistance, Resilience or Change: Post-disturbance Dynamics of Boreal Forests After Insect Outbreaks

Understanding and measuring forest resistance and resilience have emerged as key priorities in ecology and management, particularly to maintain forest functioning. The analysis of the factors involved in a forest’s ability to cope with disturbances is key in identifying forest vulnerability to environmental change. In this study, we apply a procedure based on combining pathway analyses of forest composition and structure with quantitative indices of resistance and resilience to disturbances. We applied our approach to boreal forests affected by a major spruce budworm outbreak in the province of Quebec (Canada). We aimed to identify the main patterns of forest dynamics and the environmental factors affecting these responses. To achieve this goal, we developed quantitative metrics of resistance and resilience. We then compared forests with different pre-disturbance conditions and explored the factors influencing their recovery following disturbance. We found that post-outbreak forest dynamics are determined by distinct resistance and resilience patterns according to dominant species and stand composition and structure. Black spruce forests are highly resistant to spruce budworm outbreaks, but this resistance is conditioned by the length of the defoliation period, with long outbreaks having the potential to lead the system to collapse. In contrast, balsam fir forests easily change to a different composition after outbreaks but are highly resilient when mixed with hardwood species. Overall, the severity of the disturbance and the tree species affected are the main drivers contributing to boreal forest resistance and resilience. Our procedure is valuable to understand post-disturbance dynamics of a broad range of communities and to guide management strategies focused on enhancing the resistance and resilience of the system.


Introduction
Natural disturbances are an essential component of forest ecosystems. While they play an important role in forest renewal and diversification, disturbances can generate profound changes in forest composition and structure. With the global environmental change, the occurrence and severity of biotic and abiotic disturbances have increased, affecting forests around the world (Schelhaas and others 2003;Kautz and others 2017). Such alterations in disturbance regimes may change forest responses to disturbances, resulting in altered ecosystems and affecting their function and provision of services (Reyer and others 2015;Seidl and others 2016).
Understanding forest responses to disturbances requires considering the mechanisms that enhance ecosystems to persist through time around an equilibrium state subjected to disturbances (Scheffer and others 2012;Reyer and others 2015;Willis and others 2018). Several studies recognize the important insights gained by considering resistance and resilience ('engineering resilience' or recovery) as related but distinct measurable components of ecosystem responses to disturbances (Lloret and others 2011;Hodgson and others 2015;Nimmo and others 2015). Resistance can be defined as the ability of ecological systems to persist through the disturbance event (Tilman and Downing 1994;Connell and Ghedini 2015), while resilience is the ecosystem's capacity to recover or 'bounce back' after the disturbance has been alleviated (Gunderson 2000). Given the variety of meanings around these broad concepts, it is important to establish procedures to specifically measure these two properties when assessing responses of natural communities to disturbances (Lloret and others 2011;Hodgson and others 2015;Nimmo and others 2015).
The resistance of an ecosystem to a given disturbance can be measured by the change in ecosystem structure and composition caused by the disturbance agent and can be estimated by comparing ecosystem characteristics before and immediately after the disturbance. Measuring resilience requires considering changes over a longer post-disturbance period to evaluate the degree to which an ecosystem returns to pre-disturbance conditions (D?az-Delgado and others 2002; Bagchi and others 2017).
Forest resistance and resilience to disturbances depend on multiple variables acting at different levels of organization. Species responses to disturbances mediate changes in their abundance, affecting the composition and structure of the entire community, and therefore, forest functioning as well (Oliver and others 2015). Different quantitative indices have been successfully applied to identify such variables and to compare forests based on their species responses (Bruelheide and Luginb?hi 2009;Duveneck and Scheller 2016;S?nchez-Pinillos and others 2016). While useful, the simplification inherent in quantitative, univariate indices does not capture the complexity of the entire system (Quinlan and others 2016). Further insights into forest resistance and resilience can be gained by combining quantitative indicators based on cross-scale temporal dynamics with descriptive analyses of post-disturbance successional pathways. This approach allows us to comprehensively compare the responses of different forest types while assessing potential factors affecting resistance and resilience. Moreover, ecosystem collapse can be identified as pathways leading to long-lasting and widespread changes in ecosystem state after exceeding thresholds and tipping points (Reyer and others 2015;Lindenmayer and others 2016;Bagchi and others 2017).
In this study, we developed a method of quantitatively measuring and assessing ecosystem resistance and resilience to disturbances by combining multivariate dissimilarity indices from ecosystem composition and structure with pathway analyses of these descriptors after disturbances. This method can help to answer three questions considered essential in the understanding of ecosystem responses to disturbances (see for example Willis and others 2018): (i) which types of ecosystems are the most resistant and resilient, (ii) what factors contribute to resistance and resilience, and (iii) how prone is an ecosystem to losing resilience? We applied our approach to examine boreal forest dynamics following a severe spruce budworm (Choristoneura fumiferana (Clem.)) outbreak that occurred in the second half of the 20th century in the province of Quebec, Canada (1971-1990re des For?ts de la Faune et des Parcs 2016). For this purpose, data from permanent plots repeatedly measured for Quebec's Forest Inventory (Minist?re des Ressources Naturelles 2013) provided highly relevant information to document the responses of long-lived forests to disturbances over long periods of time. Specifically, our objectives were to (1) quantify and compare resistance and resilience to spruce budworm outbreaks in sites dominated by balsam fir (Abies balsamea, L.) or black spruce (Picea mariana, (Mill.) BSP), two species that often exhibit contrasted biological responses to spruce budworm outbreaks (Hennigar and others 2008); (2) evaluate the influence of stand characteristics (species diversity, tree size structure and composition) on responses of these two host species; and (3) assess the main successional pathways of these boreal forests when affected by outbreaks as well as when considering different pre-outbreak conditions.

Study area
Our study area covers the 550,000 km2 continuous boreal forest of the province of Quebec (Canada) found between the 47 th and 52 nd latitudes (Fig. 1). This region comprises two bioclimatic domains as defined by Saucier et al. (2011): the southern area is in the balsam fir-white birch domain, mostly dominated by balsam fir and white spruce (Picea glauca, (Moench) Voss) mixed with white birch (Betula papyrifera, Marsh.); the northern portion is dominated by black spruce growing in monospecific forests or mixed with balsam fir, with abundant hypnaceous mosses and ericaceous shrubs in the undergrowth, constituting the spruce-moss domain.
The spruce budworm (hereafter SBW) is one of the most damaging outbreaking insects in the boreal and sub-boreal forests of eastern North America. SBW populations reach outbreak levels synchronously over large spatial scales, with an average periodicity of 30-40 years (Boulanger and Arseneault 2004).
As the insect mostly eats current year foliage, it takes many years of defoliation to kill a tree. During outbreak periods, continuous defoliation over 5-10 years may result in major mortality events in areas dominated by balsam fir or spruce (MacLean 1980). One of the most severe SBW outbreaks in the last century occurred between 1971 and 1990 (Bouchard and others 2006). That outbreak affected more than 50 million hectares of boreal forest in eastern North America, causing moderate to high mortality in many regions, depending on the abundance of host species and the severity and duration of defoliation (Minist?re des For?ts de la Faune et des Parcs 2016).
Both fir and spruce forests are highly susceptible to SBW attack, but they generally show contrasting responses to outbreaks. Black spruce forests are usually less defoliated than balsam fir or white spruce forests (Hennigar and others 2008) due primarily to a lack of phenological synchrony between budburst and larval emergence (Blais 1957;R?gni?re and others 2012). In contrast, balsam fir has high mortality rates due to high phenological synchrony and fewer years of accumulated foliage. However, because this species is highly competitive and presents abundant seedling banks, outbreaks have been considered to act as a self-regulating mechanism in monospecific fir forests, while more complex pathways have been suggested in mixed forests (Baskerville 1975;Morin 1994;D'Aoust and others 2004).

Individual pathways Representative pathways
Spatial data Temporal data  Figure 1: Location of forest plots in the study area sampled before and twice after outbreaks. The three states were plotted and classified in an ordination space where individual pathways were drawn by joining pre-disturbance, disturbed and final states in chronological order. Representative pathways were the result of grouping plots with the same successive forest types along the three states. Successional pathways and values of quantitative indices informed on post-disturbance forest responses.

Data source and plot selection
We used a combination of temporal and spatial data of forest structure and composition to assess the resistance and resilience of boreal forests affected by SBW outbreaks. Several repeated measures provide the necessary data to assess the magnitude of change, establishing a baseline against which changes can be assessed over time. Moreover, a large spatial scale is needed to capture the effects of environmental variables (Nimmo and others 2015) (Fig. 1).
Quebec's Forest Inventory (QFI) (Minist?re des Ressources Naturelles 2013) has been carried out approximately every 10 years since 1970 in permanent plots densely distributed in the meridional portion of the province (more than 12,000 plots below the 52 nd parallel). Each circular plot covers a total surface of 625 m 2 and consists of three nested subplots in which trees are inventoried depending on their diameter at breast height (dbh). The smallest subplot (40 m 2 ) includes measurements for all trees with a dbh larger than 1 cm. Trees with a dbh larger than 9 cm (including those with dbh > 31) are inventoried in the medium subplot (400 m 2 ). The largest subplot (625 m 2 ) considers measurements for trees with a dbh larger than 31 cm.
Our final data set included 425 forest inventory plots measured before and, at least, twice after the SBW outbreak that occurred in the 1970s and 1980s. From all QFI plots in the study area, we initially selected plots for which the relative abundance of balsam fir and spruce, in terms of basal area (BA) of adult trees (dbh > 9 cm), was more than 25% before the outbreak. We ensured a minimum outbreak impact by selecting plots for which the BA of adult host trees (balsam fir and spruce) decreased by more than 25 percentage points between the pre-outbreak and subsequent inventories. Additionally, we discarded plots with more than 50% of stumps in tree counts to avoid plots with final harvests. Despite using the abundance of adult trees in the criteria of plot selection, all trees were used in the analyses regardless of size.
To identify the outbreak period in each plot and establish a baseline state, we combined the QFI data with the SBW annual defoliation database (SBADD) between 1968 and 2003 (Minist?re des For?ts de la Faune et des Parcs 2016). The outbreak period affecting each plot was defined as the one including both the maximum loss in host species between consecutive inventories and the severe defoliation period recorded in the SBADD. We then characterized three forest states from the QFI: (i) the pre-disturbance state was obtained from the pre-outbreak or early outbreak inventory; (ii) the disturbed state was obtained from the inventory conducted immediately after the end of the outbreak period; and (iii) the final state was obtained from the last inventory available after the disturbed state ( Fig. 1; or Supplemental material for further details).

Characterization of forest plots and construction of pathways
We used a principal component analysis (PCA) to characterize forest structure and composition in the study plots through the 'vegan' R package (Oksanen and others 2017). All time points (pre-disturbance, disturbed, and final states) were jointly considered in the PCA to subsequently classify them by forest type according to a common ordination space. For this, we generated an abundance matrix compiling, for each plot and measurement (i.e. pre-disturbance, disturbed and final states), the species abundance in terms of BA. For the most common species (Abies balsamea, Picea mariana, Betula papyrifera, and Picea glauca), we separately considered three size classes using their 0.33 and 0.66 dbh percentiles as size-class breaks. We included a "No tree" category to account for plots that contained only dead trees in any measurement period. Following recommendations by Legendre and Gallagher (2001), we subjected the matrix to a Hellinger transformation to control for the many zeros in the data due to the large area covered.
We grouped all plots and measurements into forest types through a K-means partition based on the first three PCA axes, using the simple structure index to determine the optimal number of groups. We defined individual pathways for each plot by connecting in chronological order their pre-disturbance, disturbed and final states. Since each state was associated with a forest type resulting from the K-means partition, each individual pathway was characterized by the succession of K-means groups. Then, all plots exhibiting the same combination of K-means groups were combined into a representative pathway characterized by the centroids of the PCA coordinates in each sample period (i.e. pre-disturbance, disturbed and final states). To illustrate the main post-disturbance forest dynamics, successional vectors were constructed in the ordination space by connecting in chronological order the centroids that represent forest structural groups for each of the three defined states (Fig. 1).

Resistance and resilience indices
We estimated resistance and resilience based on the dissimilarities in forest structure and composition through time, rather than focusing on proportional changes in only one ecological variable (e.g. Lloret and others 2011). Resistance (Rt) was calculated as the similarity in forest structure and composition between the pre-disturbance state and the state immediately after disturbance (i.e. the disturbed state).
Therefore, the resistance index measures how stable forest composition and structure are during the disturbance. Because similarity can be calculated by computing the one-complement of common dissimilarity measurements, we defined the resistance index as follows: where D pre−dist is the dissimilarity in forest structure and composition between the pre-disturbance and disturbed states of a given plot. Thus, the resistance index ranges between zero (maximum dissimilarity) and one (identical states).
We calculated resilience (Rs) as a measure of system elasticity by considering the dissimilarity in forest structure and composition between the pre-disturbance and the final states in relation to the disturbance impact (i.e. D pre−dist ): where D pre−f inal is the dissimilarity between the pre-disturbance and final states of a given plot. The disturbance impact (D pre−dist ) was included in the formula to avoid underestimating resilience in low resistant systems in which the initial change is high and overestimating in highly resistant communities (Lloret and others 2011). Our resilience index attains a maximum value (Rs = 1) when the final state equals the pre-disturbance state (D pre−f inal = 0), while negative values represent systems that move away from the initial state (Fig. 1). Values below -1 may represent relatively small initial changes that eventually generate large final dissimilarities.
To estimate dissimilarities in forest composition and structure over time, we applied the framework provided by De C?ceres and others (2013), which incorporates the size structure of individuals in conventional community resemblance measurements. Instead of comparing a given pair of communities using their species abundances only, this procedure uses the size distribution of species in the two communities. Specifically, the abundance of each species (here quantified using BA) was stratified by diametric size classes, and then cumulative values are calculated such that the value in each class corresponds to the cumulative abundance of trees whose diameter is equal to or larger than the class value. This is called the cumulative abundance profile (CAP; De C?ceres and others, 2013

Influence of environmental variables on response to outbreaks
Because forest resistance and resilience are mainly determined by changes in host species, we built exploratory models to assess the influence of environmental factors by considering structural changes in either black spruce or balsam fir for all plots where they were present. We compiled several variables potentially affecting resistance and resilience as suggested by the literature (Table 1, Supplemental material for related literature). We included six variables related to forest structure and composition: tree species diversity (Shannon diversity index), tree size heterogeneity (Gini coefficient of diameters, Lexer?d and Eid 2006), and the relative abundance of balsam fir, black and white spruce, and white birch in terms of BA.
The set of abiotic variables included elevation, site quality index (calculated from the dominant height of host species at 50 years using equations in Pothier and Savard, 1998), presence of organic material in the soil (soil layer 5 to 40 cm thick), and three climatic variables (mean annual temperature, annual precipitation, and length of the growing season) that were obtained using the BioSIM software (R?gni?re and Bolstad 1994) following the procedure described in Terrier et al. (2013).
Since start and end dates of the outbreak vary by site and may not match the assigned pre-disturbance and disturbed states, we considered defoliation duration (i.e. the number of consecutive years recorded in the SBADD) as an additional explanatory variable. Finally, because severity of the outbreak (i.e. relative loss of host species) and time for recovery (years between the disturbed and final states) varied along our set of plots, we included them as covariate to control for undesired effects and to address potential causality due to other factors (see Supplemental material for definition of variables).
We used a beta regression model with a log-log link function to assess the factors affecting the resistance of each host species to SBW outbreaks. To avoid zero and unitary values, the response variable (Rt) was transformed as (Rt · (n − 1) + 0.5)/n, where n is the sample size (Cribari-Neto and Zeileis 2010). In the case of resilience (Rs), we used linear regression models to assess the effects of the explanatory variables on host species. Our resilience index is a proportional measurement and equivalent values can be obtained for different degrees of disturbance impact. To avoid getting the same values in contrasted cases (i.e. small versus large changes following disturbance), we only selected 50% of plots with the lowest resistance values to run the resilience models.
Selection of explanatory variables for optimal models was conducted via a supervised forward process based on differences in corrected Akaike's Information Criterion (∆AIC c ) (Sugiura 1978;Hurvich and Tsai 1989) while avoiding collinearity between variables (see Supplemental material for further details).
All analyses in the study were performed in R 3.3.2 (R Core Team 2017).

Characterization of forest plots before and after the outbreak
Forest plots were clearly discriminated by their composition (balsam fir, black spruce, and hardwoods) and structure (tree size and density) in the PCA. The first three axes of the PCA explained 51.3% of the variation in the abundance matrix. The first axis (27.3% of the variance) separated the two bioclimatic domains of the study area: plots dominated by balsam fir or white birch were located on the right side of the ordination space, whereas those dominated by black spruce were on the left side. The second axis (13.8% of the variance) was associated with deciduous abundance and conifer size structure. The lower quadrants included stands dominated by deciduous species and plots with dispersed conifers of small size (Fig. 2).
The K-means partition classified the three states of the 425 plots into seven forest types according to their composition and structure. To interpret post-disturbance dynamics and characterize each forest type, we explored compositional and structural variables and tested significant differences among forest  (Table 2). Three of these forest types were dominated by balsam fir in monospecific forests (Fir) or mixed with white birch (Fir-Birch) or black spruce (Spruce-Fir).
Monospecific black spruce forests were separated into two different groups: Spruce (high BA) represents old forests with big trees, and Spruce (low BA) represents forests of younger age and lower densities on sites with a low-quality index. The two final forest types corresponded to plots in disturbed and final states dominated by white birch (Birch) or plots without living trees (Treeless).
Because of the criteria we used to select plots, 78.4% of the pre-disturbance states were classified as balsam fir and black spruce forests, including the mixture Spruce-Fir. In the disturbed state, the percentage of plots classified as Fir or Spruce-Fir groups dropped (-25.4 and -11.3 percentage points, respectively), whereas the percentage of plots in groups dominated by white birch (Birch) or without trees (Treeless), increased (+24.3 and +10.8 percentage points, respectively). A similar number of plots dominated by black spruce (Spruce (high BA) and Spruce (low BA)) and mixed Fir-Birch forests were found in pre-disturbance and disturbed states (Table 2).

Quantitative analysis of forest resistance and resilience
At the whole community level (i.e. considering compositional and structural changes in host and non-host species), mixed forests (Fir-Birch, and Spruce-Fir) and Spruce (high BA) forests presented the highest values on the resistance index, whereas pure balsam fir (Fir) and Spruce (low BA) forests had the lowest resistance values (Fig. 3 a). Regarding resilience, we did not find significant differences among forest types (Fig. 3 a).
At the host species level (i.e. when resistance and resilience were calculated from structural changes of either black spruce or balsam fir), black spruce had significantly higher resistance (0.43 ± 0.37; mean ± sd; values considering all forest types where the target species was present), but lower resilience (0.03 ± 0.26) than balsam fir (Rt: 0.33 ± 0.29; Rs: 0.18 ± 0.25). In black spruce, the highest resistance values were found in dense mature forests and mixtures with balsam fir, whereas in low-density forests, black spruce was less resistant (Fig. 3 b). Regarding resilience, we did not find significant differences among black spruce forest types, and 70% of plots had resilience values lower or equal to zero. In contrast to black spruce, resistance values for balsam fir did not significantly vary among forest types. However, balsam fir was more resilient when mixed with white birch and less resilient when black spruce co-dominated the stand (Fig. 3 c).
Models analyzing the influence of environmental factors on black spruce indicate that the magnitude of the outbreak (severity and duration of defoliation) and the length of the growing season were key factors, significantly decreasing and increasing, respectively, resistance (Table 3). In the resilience model, the length of the recovery period was the only significant factor, with a positive effect, but showing a poor goodness-of-fit (R 2 = 0.06) because of the generalized low resilience values for black spruce (Table 3).
Regarding the balsam fir models, lower severity, higher elevation, greater tree size heterogeneity and a lower proportion of black spruce conferred increased resistance (Table 4). Fir resilience was positively related to the time for recovery, the proportion of white birch and plot elevation (Table 4).

Forest dynamics after outbreaks
Post-outbreak forest pathways varied depending on the pre-disturbance forest composition and structure.
Almost half of the dense black spruce forests (Spruce (high BA)) remained unchanged, whereas the rest evolved into either Spruce (low BA) (28.4%) or Treeless communities (10.8%) (Fig. 4 a). Similarly, the vast majority of the Spruce (low BA) plots followed one of two pathways: they did not change (56.8%), or they changed into a Treeless area (32.4%; Fig. 4 b).
In  Figure 4: Successional vectors representing the direction and magnitude of the changes in composition and structure for each pre-disturbance forest type. Representative pathways were drawn from centroids of the pre-disturbance, disturbed and final states of plots with the same successive forest types along the three states. Dots and triangles indicate centroids of pre-disturbance and final states of plots following each pathway. Line colors indicate forest type in the disturbed state. Line width is proportional to the number of plots following each pathway. Numbers represent the percentage of plots ending at the final state represented by the number color (percentage of plots was shown only for the most frequent forest types in the final state for each category). after the outbreak, with most remaining in those classes until the final state (Fig. 4 c). Similarly, 60.8% of the Fir-Birch stands transitioned to pure white birch forests immediately after the outbreak (Fig. 4 d). However, these mixed forests showed a greater capacity to recover than monospecific fir forests, with almost 50% returning to their original compositional group by the final state.
We observed a great variety of post-disturbance pathways for Spruce-Fir forests (Fig. 4 e). The most common pathways implied changes into either class of monospecific black spruce forests (Spruce

Discussion
Combining quantitative indices of resistance and resilience with pathway analyses of post-disturbance community dynamics provides insight into ecosystem responses to disturbances, particularly for coniferous boreal forests affected by insect outbreaks. Our study reveals two contrasting responses of the examined forest types to outbreaks: persistence or collapse, each driven by different mechanisms. While in black spruce forests resistance was strongly determined by the magnitude of disturbance (outbreak severity and duration) and the length of the growing season, in balsam fir forests, outbreak severity, elevation, and stand composition and structure were the main factors affecting resistance and resilience. This lead to different temporal pathways of forest structure and composition that make black spruce forests more prone to collapse into a treeless state when outbreak severity and duration increase.

Combining resistance and resilience indices with successional pathways
Understanding and measuring ecosystem responses to disturbances is critical for adapting management strategies and ensuring the persistence of ecosystem services following disturbance-mediated changes (Puettmann 2011;Ayres and Lombardero 2017). The procedure we applied in this study attempts to provide a deeper understanding of post-disturbance forest dynamics by complementing quantitative indices with descriptive analyses of post-disturbance pathways.
We developed indices that facilitate the comparison of forest responses for different pre-disturbance conditions and to explore the factors conferring resistance and resilience. Our resistance and resilience indices are based on changes in composition and structure of communities, rather than changes in one particular variable (e.g. Lloret and others 2011;Nimmo and others 2015). This makes our approach applicable to a wide range of ecosystems, communities, or populations in which disturbances generate an impact on their composition or structure. In addition, accounting for composition and structure in both quantitative indicators and descriptive analyses facilitates their complementarity. While numeric metrics allowed us to quantify and compare the impact of disturbances in terms of composition and structure for different communities, the descriptive analyses informed on the changes in the same variables along post-disturbance dynamics. A similar analysis to ours was used by Bruelheide and Luginb?hi (2009) to assess floristic changes after severe windthrows. The main difference in our approach lies in the consideration of the community structure as well as composition. By including structural changes, our metrics were able to capture differences between states with the same composition but different size distribution (De C?ceres and others 2013).
The consideration of structure in the analyses has important consequences in the definition of the disturbed and final states. Because our indices are calculated by comparing pre-disturbance, disturbed and final states, the identification of inventories to define these points is a crucial step. In the particular case of biotic disturbances, it is often difficult to define the start and end dates of an outbreak. Even when these dates are known, the potential lack of inventories conducted at the date an event such as an outbreak occurs may lead to errors in the estimation of the real impact of a disturbance. In our study, we defined the disturbed state after the whole impact of the disturbance, and we left forests the maximum period available to recover before the final state. This means that early recovery could have started before the disturbed state was inventoried and that the final state could be characterized by a more advanced developmental or successional stage than the pre-disturbance state. Consequently, resistance and resilience could be over-and underestimated, respectively. In contrast with other similar indices based on species abundance, the use of a resistance index accounting for community structure allowed us to capture differences in individuals' size, minimizing the effects of mismatches between the disturbance period and forest inventories. In addition, we minimized resistance underestimations by fitting as much as possible forest inventory surveys to estimations of insect outbreak period. This is especially useful for systems associated with long temporal development such as forests. However, systems showing rapid post-disturbance dynamics would require a good fit of the inventories to the duration of the disturbance. We ensured that the recovery period (i.e. the period between the disturbed and final states) was shorter than the pre-disturbance forest age. Although the potential resilience of the system could be underestimated due to an incomplete recovery, we were able to detect recovery trends, for instance when small trees of a lost species were present. Nonetheless, we strongly recommend a continuous assessment of resilience over time when data are available (Bagchi and others 2017).
The use of data gathered from a wide gradient of geographic locations was essential to capture multiple pathways, some of which could eventually lead to ecosystem collapse (Taylor and Chen 2011;Lindenmayer and others 2016), and allowed us to define general patterns of post-disturbance dynamics already found in studies at smaller spatial scales (e.g. Bouchard and others 2006).

Black spruce forests: resistance vs collapse
Black spruce forests showed high levels of resistance to SBW, in accordance with the phenological mismatch between budburst and larval emergence (Blais 1957). Importantly, successional pathways towards treeless areas may, however, occur as a consequence of high outbreak severity and long defoliation periods. In this sense, black spruce forests established on poor soils and low-quality sites are particularly vulnerable to collapse (Payette and Delwaide 2003). This might explain why pathways to treeless spaces were more frequent in spruce forests with low basal area, although the effects of forest productivity on spruce responses to SBW outbreaks should be further explored. Furthermore, resilience values close to zero for most plots reinforces the interpretation that long and severe outbreaks lead the system to collapse (Fig. 5). These results agree with previous studies observing alternative stable states in the boreal biome, suggesting that boreal forests can easily shift into lichen woodlands or treeless spaces in response to severe or repeated disturbances and changing conditions (Payette and Delwaide 2003;Scheffer and others 2012).
In this sense, future environmental changes, including range expansion of pests, may accelerate drastic changes towards the collapse of black spruce forests (Pureswaran and others 2015). Positive effects such as increased tree productivity may only be transitory (D'Orangeville and others 2018). Instead, some expected effects of climate change include the northward expansion of SBW, greater overwinter larval survival, phenological synchronization of larval emergence and bud break, and more stressful conditions for host species as a consequence of drought (Weed 2013;De Grandpr? and others 2018). Additionally, the occurrence of cascading severe perturbations (e.g. tree harvesting, insect outbreaks, and fire) may potentially lead to changes in the composition and structure of spruce forests, increasing their risk of collapse (Payette and Delwaide 2003; Pureswaran and others 2015).

Balsam fir forests: the importance of pre-disturbance conditions for resilience
Co-occurring tree species plays an essential role in post-outbreak dynamics of balsam fir forests (Su and others 1996; Bognounou and others 2017). The differential effect of co-occurring tree species on biotic damages suffered by a given host has been widely described in the literature (Jactel and Brockerhoff 2007; Bognounou and others 2017). The abundance of black spruce was negatively related to the resistance of balsam fir, suggesting an associational susceptibility of fir in the presence of spruce. This effect could be due to its later bud flush ensuring that SBW larvae have a food resource for longer in the season, which could maintain populations at higher levels for multiple years. Thus, fir trees re-flushing in subsequent years would be more vulnerable in the presence of spruce (Blais 1957;Despland 2017). In fir-spruce mixtures, the SBW seems to act as a biological filter to reduce fir abundance in favor of spruce (Kneeshaw and Bergeron 1998;Kneeshaw and others 2015).
In contrast to spruce, white birch abundance was positively related to the resilience of balsam fir. It has also been reported at the community and landscape levels (Bergeron and others 1995;Su and others 1996;MacKinnon and MacLean 2003;Bouchard and Auger 2014). Canopy gaps created by outbreaks enable pioneer species such as white birch to invade the space and convert stands into mixed or pure forests of shade-intolerant species (Kneeshaw and Bergeron 1998;Bergeron 2000;Taylor and Chen 2011).
Then, in agreement with previous studies (Morin 1994;Kneeshaw and others 2015), our results point to a facilitative effect of hardwoods on balsam fir regeneration and a progressive transition towards conifer forests (Fig. 5).  Figure 5: Illustrative cup-and-ball diagrams in which resistant states (black spruce) are represented by narrow deep valleys and the resilient states (Fir-Birch) by broad basins of attraction with gradual slopes permitting transitions between states. Bottom diagrams illustrate the way that variables significantly contribute to resistance and resilience of each forest type (gray circles). Arrows around gray circles represent positive factors for resistance. Arrows on solid lines leaving a forest state represent factors contributing to transitions towards a different state. On dashed lines, arrows symbolize factors involved in the recovery of forests pointed by the arrow.

Spruce
Tree size heterogeneity and elevation seem to favor balsam fir forest resistance. Size diversity suggests the presence of trees in different developmental stages that could have less damage than homogeneous mature forests in agreement with previous studies that revealed high mortality rates in mature balsam fir forests (MacLean 1980). In any case, our analyses were markedly exploratory. Further analyses based on pre-defined hypotheses should be done to more accurately test the role of environmental variables on the observed post-disturbance pathways.
Our results support the hypothesis proposed by Baskerville (1975), which states that SBW outbreaks act as a regulating mechanism for boreal forests. As such, when outbreaks occur, most mature balsam fir die, leaving shade-intolerant species and a dense regeneration of spruce and balsam fir the opportunity to grow. Eventually, the competitive balsam fir could dominate mature communities highly susceptible to subsequent outbreaks.

Conclusions and future research
Procedures combining quantitative indices and descriptive analysis of post-disturbance forest dynamics provide powerful tools for both ecosystem management and ecological assessment. Our study highlights the role of some key species in providing stability to the whole forest community, and the potential vulnerability of historically weakly affected forests due to climate-mediated changes in insect outbreak regimes. This approach may be helpful to better understand ecosystem dynamics associated with largescale disturbances, determined by different resistance and resilience patterns depending on the dominant species and composition and structure of the system.
Further research should focus on assessing post-disturbance dynamics in relation to the pre-disturbance stable state. This will require the study of undisturbed system dynamics at large spatial and temporal scales to identify stable states and to compare their dynamics with those of disturbed systems. In addition, potential applications of our procedure considering ecosystem resilience at different time scales may provide important insights into the understanding of post-disturbance ecosystem dynamics. A future challenge to research will be, however, to develop quantitatively integrated procedures accounting for the identity of returning pathways to the pre-disturbance state.