The biogeography of irruptive insect herbivores is determined by host availability and climate conditions. As such, outbreak distributions are sensitive to climatic change, especially across large latitudinal gradients. Here, we investigate the outbreak distributions of two understudied defoliators, hemlock sawfly (Hymenoptera; Neodiprion tsugae) and western blackheaded budworm (Lepidoptera; Acleris gloverana), that have both recently impacted the greatest land area recorded across the Pacific coastal temperate rainforest since the establishment of aerial survey programs. We compiled polygon-based estimates of insect damage collected by aerial observers, forest inventory, and downscaled climatic data to develop gridded estimates of bioclimatic conditions across the extent of the Pacific coastal temperate rainforest, including the continental United States, British Columbia, and Alaska. We leveraged these data to develop ensemble machine learning models with the goal of predicting the outbreak distribution of each insect. In this manuscript we: (1) describe the historical patterns of defoliator outbreaks, (2) identify and describe climatic conditions associated with outbreaks in both species, and (3) assess whether historic outbreaks have tracked geographic shifts in climate conditions across the region. We demonstrate that outbreaks of hemlock sawfly and western blackheaded budworm have been observed across the Pacific coastal temperature rainforests of North America in each decade since the establishment of the Canadian and United States aerial survey programs. The distribution of outbreaks by both insects were best explained by host availability, a limited range of spring, summer, and winter temperatures, and minimum precipitation. Finally, we demonstrate that outbreaks have tracked the poleward shift in suitable climate over the last century. This study establishes a baseline understanding of the climatic constraints and biogeographic patterns of historic sawfly and budworm outbreaks across the Pacific coastal temperate rainforest and emphasizes the overarching importance of climate in driving the irruptive dynamics of these defoliator species.

Data

Aerial surveys

We compiled aerial detection survey data describing the extent of tree damage from insects and diseases. These data were collected by aerial observers in fixed-wing aircrafts by the United States Forest Service (i.e., Alaska and the western Continental United States (CONUS); United States Detection Survey Data) and the British Columbia (BC) Ministry of Forests (British Columbia Overview Survey Data). We calculated the presence/absence of outbreak populations of hemlock sawfly and western blackheaded budworm by intersecting annual polygon and point detection survey data with a 0.25 km² grid. Grid cells that contained any visible defoliation were listed as affected. Aerial detection data were available for different time periods in each region of our analysis (Alaska: 1989-2022, BC: 1960-2022, CONUS: 1997-2022).

Forest Structure

We compiled forest structure data from raster-based estimates of hemlock (spp.) and spruce (spp.) basal area over Alaska and CONUS (per acre; 240 m raster; published in 2012 [Individual Tree Species Parameter Maps [US]) and polygon-based forest inventory records published in 2020 (Vegetation Resources Inventories [BC]). Forest inventory data from the United States were aggregated by upscaling 240 m raster-based estimates of basal area (ft²/acre) to 0.25 km² and transforming units to m²/ha. Inventory data from BC were aggregated by summing the area-weighted mean (per hectare) of forest inventory polygons that intersected each 0.25 km² cell. We used total host basal area (hemlock + spruce; m²/ha; log-transformed) for all models.

Climate

We compiled decadal climate data from raster-based estimates of annual and seasonal climate variables (1901-1910…2011-2020; 800 m raster; ClimateNA; Wang et al., 2016). We selected 14 candidate climate variables: Hogg’s climate moisture index; precipitation as snow (cm); mean annual temperature (˚C); mean annual precipitation (cm); number of frost-free days; spring, summer, autumn, and winter mean temperature (˚C); spring, summer, autumn, and winter precipitation (cm); and minimum winter temperature (˚C). Decadal climate data were intersected with the 0.25 km² grid and averaged per grid cell. Climate variables were selected based on hypothesized relationships and tested for inclusion in the final model (Supporting Information). Ultimately, we selected minimum winter temperature, mean spring temperature, mean summer temperature, total summer precipitation, and total winter precipitation for inclusion.

Models

Ensemble Model

We constructed classification ensemble (i.e., ‘model stacking’) models with the ranger implementation of random forest and the XGBoost implementation of ‘extreme’ gradient boosting. We tuned model hyperparameters based on a grid approach, although neither individual modeling approach is particularly influenced by hyperparameter selection. Based on tuning, we used Gini index as the split criterion for the random forest and grew 300 trees. We fit our gradient boosting models as probability boosted trees with the following parameters: 300 rounds (number of trees); a max depth of 9; η of 0.25 (step size shrinkage), and λ of 0.3 (L2 regularization term). Training data were selected for each ensemble model by randomly selecting 63% of observed presences per decade of observed sawfly and budworm defoliation and an equal number of absences to construct a balanced model. Predictions from each individual model (the random forest and gradient boosting) were averaged to compute a final prediction. We assessed model performances with confusion matrix metrics based on an independently sampled dataset. Since we report model performance across both the full dataset and an independent representative sample, we do not report out-of-bag performance metrics. Further, we used 1-AUC (area under the receiver-operator curve) dropout loss (Bücker et al. 2022) to assess variable importance, and accumulated local effects (Apley and Zhu 2020) to visualize relative mean predictions.

Model Workflow

We constructed a series of separate ensemble models to predict the biogeographic distribution of hemlock-spruce forestlands and each insect in this study. First, we developed a simple model to predict the presence/absence of hemlock-spruce forestlands for data visualization purposes. We provide the methods and results of this model in Supporting Information. Second, we developed models to predict the spatial extent of outbreak populations of hemlock sawfly and western blackheaded budworm. Models were fit by selecting 63% of presences per decade (1961-1970 onward) and an equal number of absences. Models were constructed 10 times for each insect species and average predictions and performance metrics are reported. Model selection was performed to select the best predictors and predictors were checked for multicollinearity (Spearman pairwise correlation < 0.7). Model selection details are provided in Supporting Information. Finally, based on our model predictions (decadal predictions from 1901-1910…2011-2022) and the thermal climate envelope based on our model results, we explored how the latitudinal range of areas suitable for western blackheaded budworm outbreaks has shifted over the duration of this study. We decided to focus on budworm outbreaks because we have a better historical record of the cyclical outbreaks in this species.

Modeling Decisions and Caveats

Analyses were restricted to hemlock-spruce defoliation attributed to either hemlock sawfly or western blackheaded budworm. This decision likely results in an underestimate of defoliation because western blackheaded budworm defoliation of spruce is similar to damage caused by other common defoliators, such as western spruce budworm, and thus may be incorrectly attributed; and there are several instances of unspecified ‘unknown’ or ‘defoliator’ reports in forests containing western or mountain hemlock. Further, the reported extent of sawfly and budworm defoliation across southeast Alaska is likely an underestimate due to the limited aerial survey coverage and frequent low cloud ceilings that obscure possible defoliation at higher elevations (see Supporting Information for more details).

Model Confidence

Constructing spatiotemporally explicit models that predict the likelihood of insect outbreaks across large geographic extents is challenging, especially across areas with low spatial and unequal temporal coverage. We are most confident in our model predictions for Southeast Alaska, Vancouver Island, and the Olympic mountains in Washington. We are the least confident in our model predictions across the coastal mountains of British Columbia and at higher elevations for several reasons: 1) the coastal mountains in BC are sparsely populated and aerial detection survey coverage is extremely sparse; and 2) forestlands across these areas are among the wettest included in our study area. As a result, we are not as confident in predicted relationships across areas where precipitation is extremely high.

Packages-All analyses were conducted in R (version 4.2.0). Integral packages included: tidyverse (wrangling; (Wickham et al. 2019), terra (geospatial formatting; (Hijmans 2024), sf (geospatial formatting; (Pebesma 2018), stars (geospatial formatting; (Pebesma and Bivand 2023), and mlr (machine learning framework; (Bischl et al. 2016).