Ecacy of Insecticide Aerial Spraying Programs to Reduce Tree Mortality During a Spruce Budworm Outbreak (1967 ‒ 1992) in the Quebec Province.

Spruce budworm (Choristoneura fumiferana) is the most important pest of conifer forests in Eastern North America. The main approach to protect Quebec forests against this defoliator is aerial spraying of insecticides. Despite the crucial role of aerial spraying of insecticides in the global forest protection strategy, little is known about its real impact on tree mortality reduction or the stand characteristics that may affect its long-term ecacy. We evaluate the ecacy of the protection program implemented between 1967 and 1992 in reducing tree mortality during the previous outbreak in Quebec and how its ecacy is affected by stand characteristics such as drainage quality and stand age. We established 422 plots in Eastern Quebec based on the following parameters: insecticide application (0-13 years), stand composition, drainage quality, and stand age at the beginning of the previous outbreak. Our results show that unprotected plots exhibited in average between 18 and 32.6% higher mortality than protected plots. Insecticide ecacy was affected by drainage quality. For example, protected plots established on sites with bad drainage quality exhibited 49% less balsam r mortality than unprotected plots. Furthermore, immature and mature protected stands exhibited a similar reduction in balsam r mortality (32.2 and 32.1% less mortality respectively) compare to unprotected stands. Our results may be useful for decision-makers not only to decide where and when apply insecticides during a spruce budworm outbreak but also, they may help to determine the pertinence of continuous protection during longer than usual outbreaks depending on protection goal. The specic goals of this study were to document the ecacy of the on to classify according and to determine and


Introduction
Forest ecosystems are the habitat for several groups of herbivores which consume foliage biomass without normally producing a negative impact on their hosts (Schowalter et al. 1986). However, some species can reach outbreak population densities and cause severe damage and mortality (MacLean 2016), producing negative consequences not only on forest biodiversity, ecosystem function and services, and loss of carbon sinks (Anderegg et al. 2013(Anderegg et al. , 2015 but also can cause important economic impact for the affected jurisdiction (MacLean 2016). This forces forest managers to develop protection programs to minimize the damage (Wainhouse 2005).
Successful protection programs have to be not only appropriate to the scale, diversity, economic and environmental value of forests but also have to take into account the ecological characteristics of the pest (Waters and Stark 1980; Wainhouse 2005). Consequently, the selection of the appropriate combination of pest management methods requires not only a clear understanding of how these techniques affect tree and stand resistance but also a thorough assessment of system responses to these decisions so that subsequent decisions are increasingly better informed (Wainhouse 2005; Ayres and Lombardero 2018).
Canadian forests suffer periodic outbreaks of pest insects which can produce great damage (MacLean 2016). The spruce budworm is the most important defoliator affecting spruce-r forests in eastern Canada. During the 20th century, three outbreaks of this defoliator occurred starting about 1910, 1940 and 1967, covering 10, 25, and 55 million ha respectively (Blais 1983a). During the rst outbreak (1910), researchers noted that budworm damage was heaviest in r-dominated stands (Tothill 1919;Swaine 1922). Consequently, forest management (e.g. removal of overmature balsam r (Abies balsamea [L.] Mill.) individuals, favor spruce-dominated stands, reduce cutting cycle, etc.) was advocated as a method of controlling budworm outbreaks (Craighead 1924;Westveld 1946). However, the scale and complexity of the silvicultural operations required and the fact that silviculture by itself would not protect the forest from spruce budworm made scientists to look for other options (Miller and Rusnock 1993).
Biological and chemical control were tested to reduce budworm damage but the latter option showed the most promising results (Blais 1973). Large-scale spraying operations of chemical insecticides (DDT) were therefore the control method employed during the second outbreak (1940) in Eastern Canada (Prebble 1975). The main objective of the chemical control program was to keep under control the damage produced by this insect. Aerial application of DDT began in 1952 and achieved its goal by keeping alive most of trees in the treated areas (Blais and Martineau 1960;Webb et al. 1961). As a result, the province of Quebec resorted to chemical control to reduce damage produce by the third budworm outbreak   (Blais et al. 1975;Blais 1977;Dorais et al. 1995). Spraying operations began in 1970s using Fenitrothion as principal insecticide. Until 1973, all areas requiring protection according to spray policy were treated (see methodology section for more details). However, the outbreak area increased so much that Quebec was forced to prioritize regions of high hazard (i.e. stand having high balsam r content) and where the forest industry was a cornerstone of economic activity (Blais 1977). The new protection policy aimed at protecting at least 50% of current year foliage to keep trees alive. This foliage protection target was based on the study of Hardy and Dorais (1976), which used the estimated contribution made by the various ages of foliage to the total photosynthetic capacity of balsam r reported by Clark (1961) to evaluate the capacity of balsam r to recover after spruce budworm defoliation. Their results suggest low balsam r mortality when 38 to 51% of its photosynthetic capacity is maintained by the end of the outbreak.
Spray e cacy is usually assessed as the percentage of spruce budworm population reduction or percentage of current-year foliage protection due to treatment (van Frankenhuyzen et al. 2007). According to Dorais et al. (1995), 80% of treated areas were adequately protected in terms of foliage protection 10 out of 14 years between 1974 and 1987 which demonstrate that the spraying operations conducted in Quebec met the protection goal. Despite the good results in terms of foliage protection and the link between the aforementioned variable and tree mortality (Hardy and Dorais 1976), there exists little information on the long-term e cacy of aerial applications of insecticides in reducing budworm-related tree mortality. For example, an early study reported that sprayed stands with DDT exhibited 11% mortality rate (2% of the basal area) whereas tree mortality reached 44% (44% of basal area) in unprotected stands (Batzer 1973). MacLean et al. (1984) reported that there was little difference in spruce and global (spruce and balsam r combined) mortality between protected and unprotected stands whereas balsam r mortality was lower in protected stands than it was in unprotected stands but this host formed only a minor component of the sampled stands. More recently, Fournier et al. (2010) reported an average reduction in tree mortality of 20.5 m3/ha in protected stands with Btk as compare with unprotected stands after 9 years. However, these studies have either been conducted during an outbreak (MacLean et al. 1984) or restricted to a particular region (Northern Minnesota (Batzer 1973) and Ottawa River Valley (Fournier et al. 2010)) which do not permit to do a complete appraisal of treatment e cacy because tree mortality continues several years after the end of the outbreak (e.g. Blais 1981a) and the fact that tree and stand characteristics may affect host and stand vulnerability (MacLean 1980;Fuentealba et al. 2013) and therefore the long-term e cacy of insecticide treatments.
It is well documented that balsam r is the most vulnerable host to budworm defoliation, followed by white spruce (Picea glauca [Moench] Voss), red spruce (P. rubens Sarg.) and black spruce (P. mariana [Mill.] BSP) (MacLean 1980; Hennigar et al. 2008). At the stand level, variables such as stand age, drainage quality, density, composition and altitude can affect tree and stand vulnerability to budworm attack (e.g. MacLean 1980, 2016). However, the relative importance of the parameters that in uence stand vulnerability is not clear and, therefore, no generalization can be made at the provincial level as to their contribution to increase/decrease tree mortality. It is therefore necessary to have a better understanding of the e cacy of protection programs in terms of tree mortality and how variables that in uence stand vulnerability affect their e cacy. This information will help us to improve our ability to predict the impact of spruce budworm outbreaks by improving the accuracy of prediction models, to increase accuracy of vulnerability classi cation systems such as the one proposed by Blais and Archambault (1982) and better plan silvicultural treatments and spraying operations to reduce the impact of this defoliator on productive forests.
The main goal of this study was to evaluate the e cacy of the protection strategy based on aerial application of chemical and biological insecticides deployed in the Province of Quebec during the previous outbreak of the spruce budworm that occurred between 1967 and 1992. The speci c goals of this study were (1) to document the e cacy of the protection program based on aerial application of chemical and biological insecticides in the province of Quebec in reducing host tree mortality and how its long-term e cacy was affected by stand age and drainage quality (two of the most important stand characteristics that affect stand vulnerability (MacLean 1980)), (2) to classify the different parameters in uencing the stands vulnerability according to their relative importance and (3) to determine the optimum number of years of treatment and the level of protection.

Methodology
Brief description of the control program implemented during the 1967-1992 outbreak in the Province of Quebec Spraying operations began 1970 and treated areas were con ned to western Quebec at rst but operations then covered many regions of the province as the outbreak spread eastward (Blais et al. 1975). During the early 70s, the goal of the spraying operations was reducing the infestation intensity in Western Quebec (1970)(1971)(1972) and eradicating the epicenters in Eastern Quebec (1971) (Dorais 1992). Until 1972, the spray policy consisted in treating stands that had sustained two successive years of severe defoliation. However, the increase in severity observed during the third outbreak together with the fact that the chemical insecticides used (Fenitrothion, aminocarb and mexacarbate) were not as effective as DDT in controlling spruce budworm damage force the province to change the spray policy. Consequently, stands were treated after one year of severe defoliation starting in 1973 (Blais 1977) and the main goal was to slow the spread of the outbreak (Dorais 1992). It was originally intended to treat all stands requiring protection according to the spray policy. This was possible until 1973 when the outbreak area increased so much that it was necessary to prioritize the stands exhibiting high vulnerability to budworm attack (Blais 1977). As a result, the goal of the protection program was modi ed in order to protect at least 50% of current-year foliage to reduce the risk of tree mortality in highly vulnerable areas (Dorais 1992).
The number of insecticide applications, the dosage and the volume used evolved during the outbreak. In 1970 and 1971, all treated areas received one application of insecticide at 1.4 L / ha (0.15 gal. / acre) aimed at the fourth-instar larvae (Blais et al. 1975). From 1972 to 1976 the standard treatment consisted in two application of insecticides at a volume between 0.9 to 1.12 L / ha (dosage: 70-280 g AI / ha). The rst application would take place when 50% of second-instar larvae had emerged whereas the second application targeted the fourth-instar larvae (Blais 1977;Dorais et al. 1995). By 1977, the infestation reached its peak in Eastern Quebec. Consequently, two to three applications at a volume between 0.9 to 1.4 L / ha (dosage: 70-280 g AI / ha) were applied between 1977 and 1980 depending on budworm density with the rst application timed at second-instar emergence. During the 1981-1983 period, the standard treatment consisted in two application of insecticides at 1.4-2.34 L / ha (dosage: 52-210 g AI / ha) (Dorais et al. 1995). However, the timing of the rst application was changed after eld trials showed that a rst application at peak fourth-instar followed by a second application at peak fth-instar larvae (late timing) protected almost twice as much foliage as the standard early timing ( rst application at second-instar emergence followed by a second one at peak fourth-instar larvae) (Blais 1981b). Consequently, the late timing became the norm and early timing was only used when larval densities and tree damage were considered excessive (Dorais et al. 1995). Chemical insecticides were used until 1986 and were then replaced in 1987 by the microbial insecticide Bacillus thuringiensis (Bt). Aerial treatments using Bt were tested since 1971 to evaluate its potential in protecting host foliage (e.g. Smirnoff 1980). These tests allowed to improve the performance and reduce the cost associated to Bt (van Frankenhuyzen 1990; Dorais et al. 1995). Further improvements in the potency of Bt formulations (e.g. Smirnoff 1985) and application technology (van Frankenhuyzen 1990) as well as environmental concerns led the Quebec Ministry of Energy and Resources decided to entirely replace chemical insecticides with this biological insecticide (Dorais et al. 1995).

Study Area
The territory considered in the study was located in the eastern part of the province of Quebec, speci cally in the regions of Gaspésie, Côte-Nord, Québec and Saguenay-Lac-Saint-Jean (Fig. 1). The study area comprised three bioclimatic domains: sugar maple (Acer saccharum Marsh.)-yellow birch (Betula alleghaniensis Britt.), balsam r-white birch (Betula papyrifera Marsh.) and balsam r-yellow birch domains. The sugar maple -yellow birch domain covers the hills bordering the southern Laurentian plateau and the Appalachians (Saucier et al. 2003). Mean annual precipitation is around 1000 mm, with 25% falling as snow. Mean annual temperature ranges from 2.5 to 5°C (Robitaille and Saucier 1998). The dominant species were sugar maple, yellow birch and red maple (Acer rubrum L.), with American beech (Fagus grandifolia Ehrhart), balsam r, and white spruce as accompanying species .
The balsam r-yellow birch domain is a transition zone between the northern temperate zone and the boreal zone. This domain includes the Gaspé peninsula, the Appalachian hills east of Québec region, the Laurentian foothills north of Saint Lawrence river and the lowlands of Saguenay-Lac-Saint-Jean (Saucier et al 2003). This bioclimatic domain exhibit a mean annual temperature of 2.5 ºC and a mean annual precipitation that varies between 900 mm and 1600 mm (Robitaille and Saucier 1998). The vegetation varies according to altitude and topography. On mesic sites, balsam r and yellow birch are the dominant species but white spruce and white cedar (Thuja occidentalis L.) are also present, while sugar maple along with yellow birch are more frequently encountered in valleys. On summits, balsam r and white birch are dominant, whereas xeric sites of lower elevation are characterized by mountain maple (Acer spicatum Lamb.) and balsam r (Grandtner 1972;Robitaille and Saucier 1998). Hydric sites in turn are dominated by balsam r and black spruce with white cedar, black ash (Fraxinus nigra Marsh.) and speckled alder (Alnus incana (Du Roi) J. Clausen) as accompanying species whereas xeric sites are dominated by balsam r, red spruce and black spruce (Robitaille and Saucier 1998).
The balsam r-white birch domain is located in the southern part of the boreal forest and the dominating species are balsam r and white spruce mixed with white birch on mesic sites. This domain covers a strip approximately 150 km wide, between the 48th and 49th parallels of north latitude from the Ontario border to the west part of the Côte-Nord. However, the perimeter of Lac Saint-Jean and of the Gaspé Peninsula belong to the balsam r-yellow birch domain (Saucier et al. 2003). This bioclimatic subdomain exhibit a mean annual temperature of 2.5 ºC and a mean annual precipitation that varies between 1000 mm and 1300 mm (Robitaille and Saucier 1998). Mesic sites are dominated by balsam r, white birch and white spruce whereas balsam r, speckled alder and black spruce in turn are the predominant species on hydric and xeric sites (Grandtner 1972;Robitaille and Saucier 1998). Forests in the Côte-Nord region are dominated by balsam r and black spruce, balsam r being the dominant species in mesic sites whereas black spruce importance increases in xeric and hydric sites and as well as in upland sites. Other species present in the region are white spruce, jack pine (Pinus banksiana Lamb.), white birch, trembling aspen (Populus tremuloides Michx.), tamarack (Larix laricina (Du Roi) K. Koch), yellow birch, and red maple (Robitaille and Saucier 1998).

Study Design
Stand selection was carried out at the end of outbreak using forest maps (inventory of the 2nd decennial) and mapping of the history of aerial application of insecticides in Quebec. The following parameters were considered to randomly select the stands: (1) (4) stand age at the beginning of the previous outbreak. This particular subset of parameters was chosen because it allowed us to include stands that presented different degrees of vulnerability to spruce budworm defoliation and various levels of insecticide treatments..
In each selected stand, one circular plot with an area of 500 m2 (radius 12.62 m) was established. The establishment of circular plots was constrained to areas that presented homogeneity in stand forest composition and hydric regime, where balsam r content comprised at least 25% of total volume, and had not suffered human or natural perturbation other than that produced by the previous spruce budworm outbreak, yielding a total of 422 circular sample plots. In each sample plot, the diameter at breast high (dbh) of all commercially sized trees ≥ 10 cm was measured by class of 2 cm. During the survey, the stems were classi ed according to their status (dead or alive). Furthermore, 3-6 balsam r per plot were selected to measure their age, dbh and height. The exposure, slope class, and water regime were also evaluated.

Statistical analysis
We compared global mortality (host tree volume mortality/total volume), host tree mortality (host tree volume mortality/ total host tree volume), balsam r, white spruce and red-black spruce mortality among treatments using generalized linear models (GLM) that accounted for binomial distributions. Preliminary analysis using likelihood ratio tests (Bates et al. 2015) showed that the candidate random effects (block and blocks nested within region) tested were not signi cant. Plots were classed by level of protection received during the last outbreak (unprotected (no insecticide application) and protected (1-13 years of treatment (insecticide application)), stand age (50 years class (plots between 41 to 60 years) and 70 years class (plots ≥ 61 years) and drainage quality (good, moderate and bad). The glm procedure in R version 4.0.2 (R Core Team 2020) was used to t a model with level of protection (unprotected vs protected), stand age (50 vs 70 years class), drainage quality (good vs moderate vs bad) and their interaction as xed effects assuming binomial distribution.
Differences between the levels of signi cant xed effects were veri ed using pairwise comparisons with Tukey's correction on the marginal means estimated by the model using the emmeans package (Lenth 2020).
To evaluate the impact of stand characteristics on stand vulnerability to spruce budworm attack, host and balsam r mortality were modeled using random forest (RF) method (Breiman 2001). RF is a nonparametric statistical learning technique that builds many decision trees from bootstrap samples of a data set. This technique can incorporate complex nonlinear interactions because is insensitive to outliers or autocorrelation, can handle many noisy variables and is robust to over tting (Breiman 2001). RF also provides measures of variable importance for each candidate predictor (Breiman 2001;Liaw and Wiener 2002). The importance of each predictor variable is measured as the change in prediction accuracy (increase in mean square error), computed by permuting (value randomly shu ed) the variable with out-of-bag data in the RF validation approach (Liaw and Wiener 2002). A larger percent increase in mean square error indicates higher importance of a variable in prediction. RF analysis for regression was conducted using the "randomForest" package (Liaw and Wiener 2002) to examine the relationship between host and balsam r mortality with stands characteristics (age, altitude, drainage, total density, host, balsam r, spruce spp. and hardwood volume). The default of 500 trees (parameter "ntree") was applied as error rate stabilized at 200-250 trees. A random one-third of predictor variables were used to perform data partitioning at each node (parameter "mtry"). The signi cance of the random forest model was tested with 1000 permutations of mortality in the "rfUtilities" package (Evans and Murphy 2016). The signi cance of the relative importance for each of the explanatory variable on each response variable was determined using the "rfPermute" package (Archer 2013). Model cross-validation was completed using leave-one-out crossvalidation (LOOCV) with the "caret" package (Kuhn 2008). In LOOCV, the prediction model is trained on data from all of the observations except one, which is "held out" and used as the test dataset. The process is repeated until all observations have served as the test data set and the accuracy results are averaged.
A 3-parameter logistic model was tted (Meyer et al. 1999) separately to evaluate the e cacy of insecticide applications (the number of years of treatment) and level of protection (number of years of successful protection (50% or more of foliage protection) / number of years of protection required) in reducing host (balsam r and spruce) and balsam r mortality. This model is suitable for quantifying growth phenomenon that exhibit a sigmoid pattern (Fekedulegn et al. 1999). However, given the inverse relationship between the number of years of treatment and level of protection and tree mortality, we used the following variation of the logistic model (i.e. inverse logistic model): Where M is either host or balsam r mortality (expressed as percentage) after the end of the last spruce budworm outbreak, N is either the number of years of treatment or level of protection, a represents the maximum mortality reduced by treatment (asymptote), b determines the number of years of treatment or level of protection at which the curve reaches the midpoint (a/2) (in ection point) and k affects the steepness of the curve at its in ection point. The logistic model is symmetric around the midpoint b (Meyer et al. 1999). The c parameter represents total potential host or balsam r spruce budworm-related mortality in the absence of protection program. We used nls function in R version 4.0.2 (R Core Team 2020) to t the models.

Results
Globally, host and balsam r mortality, varied signi cantly between treatments, age class, drainage quality and treatment × drainage class interaction ( Table 2). Depending on how tree mortality was expressed (global, host or balsam r mortality), unprotected plots exhibited in average between 18 and 32.3% higher mortality than protected plots. Similarly, mature plots (70 years class) exhibited in average 4.5% higher mortality than immature plots (50 years class). Bad drainage, in turns, had the most detrimental impact on tree survival. Indeed, plots having drainage of bad quality exhibited higher mortality than those having drainage of either moderate or good quality whereas no signi cant difference in mortality between plots having moderate and good drainage was observed (Fig. 2). Insecticide performance was also affected by drainage quality (Table 2). Although aerial applications of insecticide substantially reduced mortality in all drainage qualities, its e cacy was higher on plots having drainage of bad quality whereas unprotected plots showed increasing mortality as drainage quality worsened (Fig. 2).  White spruce and black-red spruce content were much lower than balsam r content in the study area (Table 1). White spruce mortality was affected by the interaction between treatment and drainage quality whereas black-red spruce mortality was not affected by any of the treatments (Table 3). Aerial applications of insecticides reduced signi cantly white spruce mortality only in plots having bad drainage quality (20.74 vs 6.71% in unprotected and protected stands respectively). Random forest models, as expected, explained a relative low percentage of the variation (R 2 ) in host and balsam r mortality given the high di culty that exists to relate mortality to stand characteristics (MacLean 1980; Blais 1983b) ( Table 4). Altitude and drainage quality were among the top explanatory variables in both models (Fig. 3). Indeed, altitude had the strongest relationship with host species mortality whereas drainage quality was the most important variable for predicting balsam r mortality (Fig. 3). Other important variable detected in both models were hardwood, balsam r and host content whereas stand density and spruce content importance was not always statistically signi cant (Fig. 3). Logistic growth model tted slightly better host residual content than balsam r residual content as shown by the lower root-mean squared errors (RMSE) ( Table 6). Furthermore, all parameters, except k in the models including level of protection as explanatory variable, were statistically signi cant given that con dence intervals did not include zero (Table 6) (Fekedulegn et al. 1999). However, all models tend to underestimate host and balsam r mortality (Fig. 4). This trend was expected given the high di culty to predict spruce budwormrelated mortality (MacLean 1980, Blais 1983b). Con dence intervals for parameter b (in ection point) shows that insecticide treatments reach their maximum performance in reducing tree mortality after between 5 to 8 years (average of 7 years) of continuous insecticide application ( Table 6). The return on protection investment begins to decline after this period and their contribution to reducing tree mortality becomes marginal after 10 to 12 years (Fig. 4). The models for level of protection suggest that mortality rate reduction reached its peak when level of protection is 23% (con dent interval 19-27%), reaching its equilibrium point at around 40% (Fig. 4).  Table 6 is not available with this version.

Discussion
The previous spruce budworm outbreak  is considered to be the worst outbreak in terms of duration, severity and area affected in the province of Quebec (Blais 1983a). As a result, 13 million ha suffered various degrees of tree mortality, producing estimated losses of 235 million m 3 (Dorais et al. 1995). It was estimated that more than 90% of areas treated with insecticides received adequate protection (treatments were judged to be effective if infestation levels did not exceed 25 larvae per branch) (Dorais et al. 1995) but there exist, to our knowledge, no evaluation of wood losses due to tree mortality in protected stands covering more than one region. Our results show that protected plots exhibited in average 24% balsam r mortality compared with unprotected plots which sustained in average 56.3% balsam r mortality. In other words, unprotected plots lost in average 29.9 m3/ha more balsam r volume due to spruce budworm-related mortality than protected plots. These results are in accordance with previous studies (Batzer 1973;Fournier et al. 2010), con rming that aerial spraying of insecticides can substantially reduce balsam r mortality during spruce budworm outbreaks. The higher reduction in balsam r mortality reported in this study as compared with the reduction reported by Fournier et al. (2010) (29.9 vs 20.5 m3/ha respectively) for the Ottawa River Valley highlights the importance of conducting assessments of the e cacy of protection programs at the end of the outbreak and over large areas to have a better appraisal of treatment e cacy because tree and stand characteristics may affect stand vulnerability to forest pests (MacLean 1980; Fuentealba et al. 2013) and can therefore reduce treatment e cacy by affecting stand resistance in certain regions.
Drainage quality has been found to be signi cantly related to stand vulnerability to spruce budworm attack (e.g. Archambault et al. 1990; Dupont et al. 1991). Our results from control plots are in accordance with the ndings reporting higher balsam r mortality on sites with bad quality drainage and the reduction of stand vulnerability as drainage quality improves (e.g. Archambault et al. 1990; Dupont et al. 1991) (Fig. 2). Interestingly, insecticide treatments reduced tree mortality in all drainage qualities (Fig. 2), refuting the suggestion that sites with bad drainage quality should not be included in the planning of spraying operations given the insigni cant reduction in mortality predicted for these sites after the application of insecticides (Chabot 1990). Protected plots established on sites with bad drainage quality exhibited 49% less balsam r mortality and 36.7% less host mortality than unprotected plots. The high e cacy of insecticide treatments on these sites may be explain by the impact that defoliation has on tree vigour. Following severe defoliation, growth and C levels initially decline (Piene 1980;Vanderklein and Reich 1999;Palacio et al. 2012). Trees are able to compensate to some degree for loss of resources by changing allocation patterns favouring foliage production, increasing the photosynthetic rate of surviving leaves and inducing changes in leaf morphology (Vanderklein and Reich 1999). However, trees growing on sites with bad quality drainage normally have lower energy reserves than those growing on sites with good quality drainage which reduce their ability to produce new foliage, thereby making them more vulnerable to insect damage (Dupont et al. 1991). Consequently, any reduction in defoliation in these trees may have a large impact on their chances to survive continuous defoliation by allowing the affected tree to keep the necessary amount of photosynthetic biomass to maintain its energy reserves and therefore its capacity to grow and produce foliage which is critical for tree survival (Piene 1989). These results con rm that protection programs based on aerial spraying of insecticides are effective in reducing tree mortality in all type of drainage quality during spruce budworm outbreaks.
Stand age is another characteristic that has been identi ed to in uence stand vulnerability to spruce budworm attack. Immature stands normally are less vulnerable than mature stands but mortality in the former stands is much more variable than in the latter (MacLean 1980) which may explain why high mortality rates in immature stands have been reported in some areas (e.g. Blais 1981a). Our results show that balsam r mortality is slightly higher (13.6%) on mature stand than it is on immature stands. It is noteworthy that stand age did not signi cantly affect insecticide treatments e cacy in reducing balsam r mortality which suggest that it is not necessary to assign more intensive protection programs to mature stands. Indeed, both immature and mature protected stands exhibited a similar reduction in balsam r mortality (32.

Importance Of Stand Characteristics On Tree Mortality
Insect outbreak severity (determined by insect population level), extent, and duration (years), as well as the feeding ecology of the insect, host species vulnerability and abiotic factors can affect the severity of the damage produced by forest pests in terms of growth reduction and tree mortality (Anderegg et al. 2015;Blais 1983b;MacLean, 1980MacLean, , 2016. Previous studies on stand vulnerability to spruce budworm that have focused on the impact that stand characteristics have on tree mortality have not shown strong relationships (Blais 1983b;MacLean 1980). According to Raske and Alvo (1986), the weak relationship between stand characteristics and tree mortality is due to the fact that, under outbreak conditions, most area will be severely attacked which result in similar damage level and mortality in all forest types. Prediction of tree mortality on a stand basis using stand characteristics is therefore di cult (MacLean 1980; Blais 1983b) which is in accordance with the low R 2 of our random forest models (Table 5). Consequently, it has been suggested that management and protection programs should instead take into consideration stand vulnerability in order to effectively reduce wood losses through mortality and growth reductions (Blais 1983b). As a result, the development of ecological site classi cation for spruce budworm vulnerable stands was recommended as well as the need to conduct further research to better understand the relationship and interactions between insect population densities and stand and site conditions (Witter 1985).
Despite the low predictive power of our random forest models, they allow us to better understand the relative importance that stand characteristics selected in this study have on mortality. Indeed, altitude had the strongest relationship with host mortality whereas drainage quality was the most important variable for predicting balsam r mortality (Fig. 3). The relative importance of the latter variable is well known (e.g. Archambault et al. 1990;Dupont et al. 1991) and have been included in vulnerability classi cation systems in the past (e.g. Blais and Archambault 1982) whereas the former variable has not been considered to have an important role on tree mortality. In general, stand located at elevations lower than 700 m are expected to sustain higher damage and mortality than those stand located above this elevation (Witter 1985). However, evidence of the opposite response has been observed during the previous outbreak. Blais (1980) reported a balsam r mortality rate two times higher in stand located at elevations higher than 700 m as compared with those stands located at lower elevations which concur with the impact predicted by our random forest models for this variable and the observed balsam r mortality in our plots (30.76 vs 47.61% of balsam r mortality in plots located at elevations lower and higher than 700 m respectively).. This may explain by the fact that balsam r growth at high altitude is adversely affected by growing conditions; the smaller number of buds and reduced shoot development observed in balsam r trees growing at high altitude (Blais 1979) may render trees more vulnerable to spruce budworm defoliation.
In contrast, our random forest models assign a low importance to stand age in mortality prediction which contradict the high relevance that has been attributed to this variable (e.g. MacLean 1980; Witter 1985). This response may stem from the fact that insect outbreaks always have stochastics elements that may produce unusually severe tree mortality regardless of stand characteristics (Blais 1983b; Raske and Alvo 1986; MacLean 2016). In general, immature stands sustain less mortality than mature stands (MacLean 1980; Witter 1985) but mortality rates in immature stands can reach levels as high as those observed in mature stands (e.g. Blais 1981a). This phenomenon may be produced by extremely high larval densities, when the duration the outbreak at the local level is longer than usual (MacLean 2016), by the occurrence of abiotic event, such as droughts, that may stress trees (Anderegg et al. 2013;, by continuous migrations of moths from adjacent areas (e.g. Raske and Alvo 1986) or by the spatial variation in forest landscape structure (i.e. species composition and spatial con guration) which may in uence budworm outbreak frequency, intensity and synchrony (e.g. Robert et al. 2018). Consequently, implementation of effective protection programs requires not only a general knowledge of the impact of stand characteristics on stand vulnerability but also how distinctive features of regional forest can alter its vulnerability. In this regards, our random forest models shed some light on how some distinctive features such as altitude (increasing mortality with altitude) may have a greater impact on tree mortality than previously acknowledge whereas the impact of other stand characteristics such as stand age may be affected by unpredicted events that may increase the severity of outbreaks as it was observed in the previous spruce budworm outbreak (Blais 1983a; Raske and Alvo 1986). Vulnerability classi cation systems should therefore be modi ed accordingly. Furthermore, prediction of stand vulnerability (or expected mortality) should be conducted assuming different levels of population densities (e.g. different outbreak severity levels) and outbreak duration to determine potential damage under realistic, pessimistic and optimistic scenario which is fundamental to evaluate the potential impacts and to develop protection and management plans according to different outbreak severities.

Optimum Level Of Protection
During a spruce budworm outbreaks, host trees sustain high defoliation repeatedly over a period of 5 to 15 years (Candau et al. 1998;Gray et al. 2000) which leads to substantial losses in tree vigour and high mortality unless a protection program is implemented. Our results show that aerial spraying of insecticides show a high performance reducing tree mortality for a period of 5 to 8 years, 7 years being the most likely point where these aerial spraying reach their highest performance (Fig. 4). After this point, the return on protection investment begins to decline, suggesting that further application of treatment would not be warranted unless stands sustain severe defoliation lasting 12 years or more. These results con rm that the control program used in the province of Quebec is suitable to protect trees against spruce budworm attack. However, the duration of the outbreaks may affect its overall e cacy. Furthermore, our models for level of protection suggest that it is not indispensable to achieve the target protection goal (50% or more foliage protection) every year during an outbreak to reduce host tree mortality (Fig. 4). Consequently, insecticide application every 2 years may be a relevant option to protect host species given the adequate level of protection provided by this strategy in the current outbreak and the reduction in the number of applications required (Fuentealba et al. 2019). However, protection during the whole outbreak may be warranted if the main goal is to protect tree growth.

Conclusions
Our results show that the protection program based on aerial application of insecticides implemented in eastern Quebec during the previous outbreak to protect forests against spruce budworm provided a good protection by effectively reducing tree mortality. However, its e cacy was affected by stand characteristics such as drainage quality. Given that insect outbreaks always have stochastics elements that may produce unusually severe tree mortality regardless of stand characteristics (Blais 1983b; Raske and Alvo 1986; MacLean 2016), implementation of effective control programs requires not only a general knowledge of the impact of stand characteristics on stand vulnerability but also how distinctive features of regional forest can alter its vulnerability. In this regard, our results shed some light on how some distinctive features such as altitude (increasing mortality with altitude) may have a greater impact on tree mortality than previously acknowledge. Furthermore, our results may be useful for decision-makers not only to decide where and when apply Btk during a spruce budworm outbreak but also they may help to determine the pertinence of continuous protection during longer than usual outbreaks depending on the protection goals. Finally, the assessment of the protection programs is important because it allows decision-maker to make better informed decisions to protect forests, develop better protection strategies (Wainhouse 2005; Ayres and Lombardero 2018) and to obtain pertinent information to evaluate the pro tability of protection programs based on aerial application of insecticides. Indeed, an ongoing project is testing ve different protection scenarios that were conceived using this information (Fuentealba et al. 2019). Further research will be conducted to evaluate the impact of aerial spray program on tree growth so that the impact on total tree volume protection (mortality + growth) could be estimated. Distribution of plots in the study area. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.

Figure 2
Balsam r mortality, global (host tree volume mortality/ total volume) and host mortality according to drainage quality and level of protection (LSMEANS ± SE). Figure 3