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 confined to western Quebec at first 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–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 modified 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 first 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 first 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 first application was changed after field trials showed that a first application at peak fourth-instar followed by a second application at peak fifth-instar larvae (late timing) protected almost twice as much foliage as the standard early timing (first 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, specifically 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 fir-white birch (Betula papyrifera Marsh.) and balsam fir-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 fir, and white spruce as accompanying species .
The balsam fir-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 fir 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 fir and white birch are dominant, whereas xeric sites of lower elevation are characterized by mountain maple (Acer spicatum Lamb.) and balsam fir (Grandtner 1972; Robitaille and Saucier 1998). Hydric sites in turn are dominated by balsam fir 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 fir, red spruce and black spruce (Robitaille and Saucier 1998).
The balsam fir-white birch domain is located in the southern part of the boreal forest and the dominating species are balsam fir 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 fir-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 fir, white birch and white spruce whereas balsam fir, 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 fir and black spruce, balsam fir 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) number of years of treatment (insecticide applications) (0–13 years), (2) stand composition (pure stands of balsam fir, balsam fir- spruce ssp. or balsam fir – hardwood mixed stands), (3) drainage quality (good [subhydric with seepage (classes 31 and 41], moderate [mesic (classes 20 and 30)] or bad [ xeric, hydric (classes 10, 40, 50, 51, 60, 61)] (for more details on drainage classes see MFFP 2019)), and (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 fir 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 classified according to their status (dead or alive). Furthermore, 3–6 balsam fir 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 fir, 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 significant. 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 fit 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 fixed effects assuming binomial distribution. Differences between the levels of significant fixed effects were verified 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 fir 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 overfitting (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 shuffled) 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 fir mortality with stands characteristics (age, altitude, drainage, total density, host, balsam fir, 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 significance of the random forest model was tested with 1000 permutations of mortality in the “rfUtilities” package (Evans and Murphy 2016). The significance 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 cross-validation (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 fitted (Meyer et al. 1999) separately to evaluate the efficacy 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 fir and spruce) and balsam fir 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):
M = c ‒ a/(1 + exp(‒k(N-b)),
Where M is either host or balsam fir 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) (inflection point) and k affects the steepness of the curve at its inflection point. The logistic model is symmetric around the midpoint b (Meyer et al. 1999). The c parameter represents total potential host or balsam fir 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 fit the models.