Ecography

ECOGRAPHY
Research
Landscape host abundance and configuration regulate periodic
outbreak behavior in spruce budworm Choristoneura fumiferana
Louis-Etienne Robert, Brian R. Sturtevant, Barry J. Cooke, Patrick M. A. James, Marie-Josée Fortin,
Philip A. Townsend, Peter T. Wolter and Daniel Kneeshaw
L.-E. Robert (http://orcid.org/0000-0002-4406-4079) (louis-etienne.robert@umontreal.ca) and P. M. A. James, Dépt de sciences biologiques, Univ. de
Montréal, Montréal, QC, Canada. – B. R. Sturtevant, Inst. for Applied Ecosystem Studies, Northern Research Station, USDA Forest Service, WI, USA.
– B. J. Cooke, Ontario Ministry of Natural Resources and Forestry, Sault Ste. Marie, ON, Canada. – M.-J. Fortin, Dept of Ecology and Evolutionary
Biology, Univ. of Toronto, Toronto, ON, Canada. – P. A. Townsend, Russell Labs, Madison, WI, USA. – P. T. Wolter, Iowa State Univ., Dept of Natural
Resource Ecology and Management, Ames, IA, USA. – D. Kneeshaw, Centre d’étude de la forêt (CEF), Univ. du Québec à Montréal, Montréal, QC, Canada.
Ecography Landscape-level forest management has long been hypothesized to affect forest insect
41: 1556–1571, 2018 outbreak dynamics, but empirical evidence remains elusive. We hypothesized that
doi: 10.1111/ecog.03553 the combination of increased hardwood relative to host tree species, prevalence of
younger forests, and fragmentation of those forests due to forest harvesting legacies
Subject Editor: Mikko Mönkkönen would reduce outbreak intensity, increase outbreak frequency, and decrease spatial syn-
Editor-in-Chief: Hanna Tuomisto chrony in spruce budworm Choristoneura fumiferana outbreaks. We investigated these
Accepted 13 November 2017 hypotheses using tree ring samples collected across 51 sites pooled into 16 subareas
distributed across a large ecoregion spanning the international border between Ontario
(Canada), and Minnesota (USA). This ecoregion contains contrasting land manage-
ment zones with clear differences in forest landscape structure (i.e. forest composition
and spatial configuration) while minimizing the confounding influence of climate.
Cluster analyses of the 76-yr time-series generally grouped by subareas found within
the same land management zone. Spatial nonparametric covariance analysis indicated
that the highest and lowest degree of spatial synchrony of spruce budworm outbreaks
were found within unmanaged wilderness and lands managed at fine spatial scales in
Minnesota, respectively. Using multivariate analysis, we also found that forest compo-
sition, configuration, and climate together accounted for a total of 40% of the variance
in outbreak chronologies, with a high level of shared variance between composition
and configuration (13%) and between composition and climate (9%). At the scale of
our study, climate on its own did not explain any of the spatial variation in outbreaks.
Outbreaks were of higher frequency, lower intensity, and less spatially synchronized in
more fragmented, younger forests with a lower proportion of host species, with oppos-
ing outbreak characteristics observed in regions characterised by older forests with
more concentrated host species. Our study is the first quantitative evaluation of the
long-standing ‘silvicultural hypothesis’ of spruce budworm management specifically
conducted at a spatio-temporal scale for which it was intended.
Keywords: spruce budworm, harvest disturbance, landscape ecology, forest management

legacies, dendrochronology, outbreak synchrony
––––––––––––––––––––––––––––––––––––––––
© 2017 The Authors. Ecography © 2017 Nordic Society Oikos
www.ecography.org
1556
16000587, 2018, 9, Downloaded from https://nsojournals.onlinelibrary.wiley.com/doi/10.1111/ecog.03553 by Universitaet Du Quebec A Montrea, Wiley Online Library on [22/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Introduction control (Jones 1977, Hassell et al. 1999). What was missing
from these early models was robust empirical evidence for the
Outbreaks of boreal forest insects are a major forest distur- hypothesized effect of the forest on the budworm.
bance that have significant ecological and economic conse- In parallel, the ‘silvicultural hypothesis’ of spruce bud-
quences. Insect outbreaks affect biogeochemical processes and worm management posited that the impact of the spruce
can influence nutrient cycling, hydrology, and atmospheric budworm in the 20th century was accentuated by ‘man’s
carbon cycling, all of which can affect climate at a global influence in the forest’ – i.e. regional-scale commercial har-
scale (Kurz et al. 2008, Dymond et al. 2010). Insect out- vesting of spruce (Picea spp.), and its replacement by bal-
breaks also affect forest landscape structure (i.e. stand species sam fir – that concentrated preferred host abundance and
composition and spatial configuration) with consequences increased outbreak extent and intensity (Blais 1965, 1983).
for forest succession (Baskerville 1975), timber production Baskerville (1975) argued that the reciprocal feedback
(Chang et al. 2012), and the risk of future disturbances such between the budworm and the forest was strong, although
as fire (James et al. 2017). Miller and Rusnock (1993), who accepted the impact of
The notion that reciprocal feedback between forest and budworm on the forest and that broad-scaled concentra-
insects at least partially drives periodic population eruption tions of mature balsam fir were more vulnerable to budworm
was central to the earliest formalized concepts of forest-insect damage, emphasized the lack of supporting evidence for any
ecosystem dynamics (Holling 1973), and continues to be influence of such concentrations of balsam fir on budworm
a central theme in ecosystems ecology (Raffa et al. 2008). outbreak dynamics. More recent empirical studies have
An important role of the forest on outbreaks has been dem- shown that: 1) budworm defoliation consistently tends to be
onstrated in multiple defoliator systems, including the jack lower when hardwoods are abundant (Bergeron et al. 1995,
pine budworm Choristoneura pinus (Volney and McCullough Su et al. 1996, Campbell et al. 2008, Gray 2013), and 2)
1994), forest tent caterpillar Malacosoma disstria (Roland this effect is due, in part, to an increased abundance of key
2005), gypsy moth Lymantria dispar (Haynes et al. 2009), parasitoid species associated with this hardwood component
and winter moth Operophtera brumata (Wesołowski and (Cappuccino et al. 1998, Quayle et al. 2003). This has led
Rowiński 2006). Although the top-down effects of natu- to a broadening of the silvicultural hypothesis to include the
ral enemies typically underlie cyclic outbreak behaviour effect of forest composition on natural enemy interactions.
(Berryman 2002, Turchin 2003), there is also evidence Still missing is convincing evidence for the penultimate effect
that some of this outbreak behaviour may be attributable of forest landscape structure (i.e. the age, amount, and spatial
to bottom-up host-plant effects (Price et al. 1984, Nealis arrangement of host vs nonhost tree species) on the dynamics
and Lomic 1994, White 2012). For example, impact of the of spruce budworm outbreaks.
European spruce bark beetle Ips typographus is determined by The relative importance of top-down versus bottom-
available deadwood and its spatial distribution (Lausch et al. up processes in budworm outbreak dynamics has been a
2013). Another example is the larch budmoth Zeiraphera subject of intense empirical debate for nearly a century
griseana in the European Alps, where gradients in density of (Sturtevant et al. 2015, Pureswaran et al. 2016). The idea of
host tree species induced epicenter formation and travelling forest-driven feedback to budworm outbreaks was dismissed
waves within outbreaks (Johnson et al. 2004). Further, the by Royama (1984), who argued that the core population oscil-
ability of natural enemy communities to respond to popula- lation underlying periodic outbreaks is produced by delayed
tion irruptions appears to be linked to forest landscape con- density-dependent mortality due to natural enemies, and that
ditions (Roland and Taylor 1997, Cappuccino et al. 1998), spatial outbreak synchrony is produced by regionally-corre-
suggesting there are important interactions among trophic lated factors such as weather (i.e. the Moran effect (Moran
levels affecting outbreak dynamics. 1953, Royama 2005)). However, more recent research sug-
The spruce budworm Choristoneura fumiferana is the gests that complex cycling of outbreaks may emerge from a
most damaging forest insect defoliator in North America, combination of both top-down and bottom-up reciprocal
with balsam fir Abies balsamea and white spruce Picea glauca feedback processes (Cooke et al. 2007, Régnière and Nealis
being especially vulnerable (MacLean 1980, Hennigar et al. 2007, Sturtevant et al. 2015).
2008). This native insect is a strong flier capable of long Considering the full weight of the evidence, both
distance dispersal (Greenbank et al. 1980) and best known empirical and theoretical, one would expect that reducing
for the broad-scale spatial synchrony, comparatively low host forest cover might affect the frequency and intensity
frequency, and long duration of its outbreaks (Blais 1983, of budworm outbreak cycles, through its effects on local
Royama 1984, Williams and Liebhold 2000, Jardon et al. budworm dispersal success (Nealis 2016) and on natural
2003). enemies (Royama 1984, Eveleigh et al. 2007). One might
The spruce budworm inspired much of the foundational also expect a reduction in host to reduce spatial synchrony,
thinking about forest-insect feedbacks as the cause of cross- through its effect on long-range dispersal success (Régnière
scale eruptive dynamics (Ludwig et al. 1978). Indeed, the early and Lysyk 1995, Peltonen et al. 2002) and cycle amplitude
‘multiple equilibrium’ models of spruce budworm outbreaks (Barbour 1990, Cooke et al. 2007). This would represent a
hypothesized that abundance of host foliage was a principle coherent multivariate response of outbreak characteristics
bottom-up driver enabling budworm escape from top-down that are functionally related to one another, and modified
1557
through changes in forest landscape structure. Robert et al. the effects of forest landscape structure on the characteristics
(2012) used tree-ring data from Ontario and Minnesota of a single outbreak (Bergeron et al. 1995, Su et al. 1996,
to show that spruce budworm outbreaks during the 20th Candau and Fleming 2005, Campbell et al. 2008, Gray
century tended to affect fewer trees in managed landscapes 2013), our goal is to examine the long-term (i.e. over mul-
than in unharvested wilderness areas. In particular, they tiple outbreaks) effects of forest landscape structure on out-
found that outbreak intensity was reduced in areas where break dynamics at scales both relevant to that of dispersing
the spatial pattern of forest harvesting was diffuse rather spruce budworm (i.e. 10’s to 100’s of km; (Greenbank et al.
than spatially aggregated. These results were qualitatively 1980, Anderson and Sturtevant 2011) and consistent with
consistent with a forest landscape structure feedback to out- the silvicultural hypothesis more generally (Miller and
break dynamics. Outbreaks within older, less fragmented Rusnock 1993).
forests (wilderness) exhibited infrequent, spatially synchro-
nous outbreaks. In contrast, younger, managed forests con-
taining a greater proportion of hardwood species exhibited Methods
outbreaks that were more frequent, but less synchronized
– all in accordance with the predictions of the silvicultural Study area
hypothesis.
In this paper, we extend the study of Robert et al. (2012) The Border Lakes Landscape (BLL) contains a large
to quantitatively investigate the effect of forest landscape (~20 000 km2) ecoregion that straddles the border
structure on spatio-temporal variation in budworm outbreak between Minnesota (USA) and Ontario (Canada) (Fig. 1;
dynamics. We hypothesized that attributes of forest landscape (Sturtevant et al. 2014)) at the transition between the Great
structure related to landscape management practices, includ- Lakes-St Lawrence mixed-wood and boreal forest regions.
ing forest composition (host abundance, hardwood content, The full extent of the BLL includes a 50 km buffer to address
forest age) and configuration (disturbance rates, host con- potential edge-effects related to insect movement (Anderson
nectivity), influence budworm outbreak behavior (i.e. cycle and Sturtevant 2011). Forest composition is mixed ‘near
intensity, cycle frequency, and spatial synchrony) accord- boreal’ (Heinselman 1973) with a high proportion of boreal
ing to the expectations of a multi-variate response outlined tree species (e.g. jack pine Pinus banksiana, black spruce Picea
above. We further wished to distinguish the effects of forest mariana, white spruce Picea glauca, balsam fir A. balsamea,
landscape structure on outbreaks from the effects of climate birch Betula papyrifera, trembling aspen Populus tremuloides),
(Gray 2008, 2013). Although other studies have examined as well as several species near the northern limit of their range,
0 25 50 100 Km
OWN
OEN
OWS
WWN WNN
OES
WEN
MWN
WWS
WES MEN
MWS MCN WSS
MES Lake
Superior
MCS
Legend
Wilderness
Minnesota
Ontario
Border Lakes
Ecoregion
Figure 1. The Border lakes landscape (BLL) study area located at the border between Ontario (Canada) and Minnesota (United States).
Points represent sampling sites for dendrochronological reconstructions of outbreaks that are grouped by the subarea unit of analysis.
Subarea labels: 1st letter = management zone (Ontario, Wilderness, Minnesota); 2nd letter = area (East, West, ‘Central’, South, North); 3rd
letter = subarea (North, South).
1558
such as white pine Pinus strobus, red pine Pinus resinosa, and recent large wind and fire disturbances. Forest composition
red maple Acer rubrum. (i.e. tree species basal area) was also mapped via remote sens-
The entire BLL shares a common early logging history. ing by Wolter et al. (2008) ca 2002 (Supplementary material
Region-wide logging started at the end of the 19th century Appendix 1 Fig. A1B) and also ca 1985 (Wolter unpubl.)
and focused on selective harvest of the ‘big pines’ (P. strobus, using archived Landsat TM imagery. Wavelet analysis applied
P. resinosa). Logging patterns diverged sharply between to the 2002 dataset (James et al. 2011) indicated that man-
Minnesota and Ontario approximately 70 yr ago, coinciding agement zones differed somewhat in terms of the basal
with the advent of mechanized harvesting on both sides of area of spruce budworm hosts (spruce and fir), especially
the border. Cut-block sizes in managed Canadian forests are between Minnesota and Wilderness (Supplementary material
an order of magnitude larger in size than managed American Appendix 1 Fig. A1B). The managed parts of the landscape
forests, though the rates of harvest have been historically (Ontario and Minnesota) contained similar forest types, but
similar (James et al. 2011, Sturtevant et al. 2014). Between had more early-successional, deciduous forest relative to the
these managed regions lies an approximately 10 000 km2 Wilderness (James et al. 2011).
wilderness area that includes Quetico Provincial Park in
Ontario and the Boundary Waters Canoe Area Wilderness Sampling design
in Minnesota, where no timber harvest has occurred since
the early 1970s (Heinselman 1996), referred to henceforth as We reconstructed spruce budworm outbreak histories using
‘Wilderness’. Forest fires have decreased over the last century tree-ring analyses of white spruce tree cores sampled from
(Ahlgren and Ahlgren 2001), enhancing succession toward 51 sites distributed throughout the BLL (Fig. 1). White
spruce and fir. Periodic defoliation by spruce budworm often spruce was selected because it is as susceptible as balsam fir
kills balsam fir in this region (Frelich and Reich 1995). to defoliation by spruce budworm, but has a higher prob-
ability of surviving attack (Hennigar et al. 2008). We tar-
Forest landscape structure of the study area geted mesic sites that contained five to fifteen old, large
canopy white spruce ( 30 cm DBH) on flat or mid-slope
Divergence in the spatial patterning of disturbances between topographic positions. Two cores were taken from each
management zones in the BLL has been documented spruce tree at a height of 1 m. Between 5 and 10 trees were
by a time-series of Landsat-derived land cover maps from sampled per site (Table 1). Groups of three neighbouring
(1975–2000, (Wolter et al. 2012) (Supplementary mate- sites were then aggregated into ‘subareas’ to provide suffi-
rial Appendix 1 Fig. A1A). Prior analyses of these data by cient sample size (n  15 trees) for outbreak reconstruction
Sturtevant et al. (2014) showed that: 1) disturbance patch- (Fig. 1). In total we created 16 subareas that represent a strati-
size distributions in the BLL were the most consistent fied sampling of the BLL based on forest management region
through time in Minnesota managed zones, and least con- ( Wilderness = ‘W’, Minnesota = ‘M’), lon-
sistent spatially and temporally within the Wilderness; 2) gitude (East = ‘E’, Central = ‘C’, West = ‘W’), and latitude
Minnesota cut-blocks were both small in size and more (North =’N’, South = ‘S’). The goal of this stratification was
diffuse than the aggregated, large-sized cut-block patterns to account for climatic gradients across the study area. The
in Ontario; 3) natural, stand-replacing disturbances in the minimum distance between subarea centroids was 15 km and
Wilderness were comparably rare, with the exception of a few the maximum distance was 250 km. A minimum of fifteen
Table 1. Summary statistics of the dendrochronological reconstruction for all subareas produced by the program COFECHA.
Sensitivity
Zone Area Subarea Label No. of trees Time span Mean SD Inter-tree correlation
Ontario Eastern North OEN 18 1849–2005 0.25 1.13 0.496
Eastern South OES 22 1883–2005 0.226 1.024 0.472
Western North OWN 19 1895–2005 0.28 1.157 0.604
Western South OWS 17 1926–2005 0.272 1.159 0.571
Wilderness Eastern North WEN 16 1924–2006 0.279 1.029 0.661
Eastern South WES 13 1902–2006 0.25 1.105 0.506
North North WNN 20 1906–2005 0.268 1.356 0.643
South South WSS 24 1897–2005 0.277 1.026 0.604
Western North WWN 21 1836–2006 0.259 0.92 0.426
Western South WWS 15 1880–2006 0.253 0.826 0.594
Minnesota Center North MCN 21 1915–2005 0.208 1.351 0.636
Center South MCS 23 1928–2005 0.207 1.413 0.583
Eastern North MEN 22 1836–2005 0.238 1.403 0.368
Eastern South MES 22 1872–2005 0.249 1.075 0.609
Western North MWN 22 1903–2005 0.277 1.363 0.63
Western South MWS 22 1920–2005 0.263 1.508 0.61
1559
non-host coniferous trees (i.e. jack pine and red pine;  30
Table 2. Summary of potential and selected predictor variables to describe variation in SBW outbreak attributes at each sample subarea (Fig. 1; see text for details). Dates refer to
the time periods for which this predictor was available. Predictors indicated as being calculated at multiple spatial scales were calculated at using multiple extracted buffers of
Landsat MSS and Landsat TM

cm DBH) per subarea were also cored to serve as a control for
local climatic variation. All cores were stored in plastic straws
McKenney et al. 2007

Data source
and were later mounted and sanded using progressively finer
grit (80, 150, and 220).
Landsat TM
Landsat TM
Landsat TM
Landsat TM
Landsat TM
Outbreak history reconstruction
BIOSIM
Budworm outbreaks were reconstructed from spruce tree ring
data using the methods described in Robert et al. (2012).
Site-level chronologies were built using tree-ring widths
Multi-scale
measured using a Velmex uni-slide measuring table with an
Yes
Yes
Yes
Yes
Yes
Yes
Yes
Yes
Yes
Yes
Yes
No
No
No
No
accuracy of 0.001 mm connected to a computer (Velmex
Incorporated, Bloomfield, New York, USA). Widths were
cross-dated using the program COFECHA (Holmes et al.
Average of 1980–2000, 5 yr increment

1986) which was also used to locate missing or false rings and
to test for measurement errors. Chronologies were then aggre-
Candidate landscape factors affecting budworm dynamics

gated into their respective subareas as described above. Index
chronologies were calculated using the program ARSTAN
Time period
(Cook 1985, Holmes 1992) with a cubic smoothing spline
Average of 1971–2000
to detrend the series and remove age-related growth trends
(Cook 1985). Spline parameters were set to a 50% frequency
1985 and 2002

1985 and 2002
1985 and 2002
1985 and 2002

1985 and 2002
response cut-off of 60 yr, which is standard for spruce bud-
worm outbreak reconstruction (Boulanger and Arseneault
2004, Bouchard et al. 2006, Robert et al. 2012). Between 13
and 24 trees (mean = 19.8 trees; SD = 3.2 trees) were used to
1975
2000
develop chronologies for each subarea (Table 1).
At the level of an individual tree, outbreaks were defined
CoheHst[SpatialScale][Year]
CoheFir[SpatialScale][Year]
when a tree-ring growth reduction was observed for more
Decid[SpatialScale][Year]
POld[SpatialScale][Year]
AWMOpn[SpatialScale]
PFor[SpatialScale][Year]
than 5 yr and at least one year showed a reduction greater
AwmPtch[SpatialScale]
Hst[SpatialScale][Year]
Fir[SpatialScale][Year]
Abbreviation
than 1.28 standard deviations from the mean chronology ring
POpn[SpatialScale]
PDist[SpatialScale]
width (Boulanger and Arseneault 2004, Bouchard et al. 2006).
Outbreak histories, based on the percentage of trees classi-
fied as affected by an outbreak at a given year and scale, were
WintTmp
SumTmp
ClimInd
reconstructed using these parameters and the program OUT-
SpTmp
BREAK (Holmes and Swetnam 1996). OUTBREAK has
increasing radii (0.5, 1, 5, 10, 20, 30 km) centered on each subarea.
been previously used in spruce budworm outbreak reconstruc-

tion and can detect outbreaks while avoiding confounding
Unit
Mean
Mean
Mean
Mean
Mean
Mean
Mean
Mean
Mean
Mean
Index
consequences of other defoliating insects (Boulanger and Arse-
%
%
neault 2004, Bouchard et al. 2006). No other defoliators are

known in the region that are capable of causing growth reduc-
Area-weighted patch size of openings
tions over multiple years such that they would be identified

Configuration Area-weighted patch size of forest
COHESION of spruce-fir patches
as an outbreak. Outbreak detection with the program OUT-

Proportion of forest disturbance
Spruce budworm growth index
BREAK was more efficient when applied solely to host spe-

Proportion of forest opening
Min. summer temperature
cies (Bouchard et al. 2006). Non-host chronologies were only

COHESION of fir patches
Spruce and fir basal area
Min. winter temperature

Min. spring temperature
Variable
Proportion of old forest
used to validate detected outbreaks by visual comparison of

Deciduous basal area
Composition Balsam fir basal area
the growth pattern (not shown). The percentage affected trees

Proportion of forest
by outbreaks each year at either the site or subarea scales were

disturbance
used as the response variable for cluster analysis, redundancy

analysis, and analysis of spatial synchrony (described below).
Predictor variables
Forest landscape structure
Forest landscape structure was spatially quantified in two
Category
Climate
ways, each using separate sets of multi-temporal Landsat

sensor data (Table 2): 1) forest configuration was quantified
1560
using discrete land use land cover time series maps at five-year COHESION were calculated using FRAGSTATS (McGa-
intervals from 1975–2000 (i.e. 6 points in time) (Wolter et al. rigal et al. 2012).
2012); and 2) forest composition was quantified using use
maps of spatially continuous estimates of basal areas (m2 ha–1) Climate
for balsam fir, combined budworm host species (i.e. fir and We assessed the potential effect of climate on outbreaks
spruce), and deciduous tree species. Composition was quan- within our study area using mean minimum seasonal (win-
tified at only two points in time (1985 and 2002) due to ter, summer, spring) historical (1901–2000) temperature
the availability of forest plot data necessary to develop the for each subarea (McKenney et al. 2007) (Table 2). Mini-
underlying tree species models (Wolter et al. unpubl. and mum seasonal temperatures have been previously shown to
Wolter et al. (2008)). influence spruce budworm outbreaks (Swetnam and Lynch
Because we do not know a priori the scales at which bud- 1993, Candau and Fleming 2005, 2011). We also examined
worm outbreaks respond to forest conditions, we quantified the influence of a spruce budworm population growth index
our forest landscape structure predictor variables at multiple (Table 2, Supplementary material Appendix 1 Fig. A1C) on
scales using a range of sample radii (i.e. 0.5, 1, 5, 10, 20, 30 outbreaks averaged between 1970 and 2000. This generation
km), where each predictor variable was measured within a survival probability index is based on a model of budworm
circle defined by a given sample radius for each site, and site- physiological development and in response to climatic varia-
level variables were then averaged across sites within a subarea tion (Régnière et al. 2012).
(approximately 3 sites per subareas).
Both forest composition and configuration were charac- Statistical analysis
terized using multiple metrics, where the specifics of each
metric were contingent on the corresponding data available Time-series clustering
at different points in time. Forest composition was repre- Using complete linkage cluster analysis (Keogh and Kasetty
sented using the basal area of combinations of budworm host 2003) applied to time series of raw percentage affected trees,
and non-host species, and the proportion of either forest area we identified subareas that exhibited similar historical out-
(circa 1975) or ‘old forest’ area in 2000 (see below). We cal- break dynamics. Prior to clustering, time series were trimmed
culated the average basal area of balsam fir, spruce combined so that each year examined with a subarea time series con-
with fir, and the combination of all deciduous tree species for tained information on at least five individual trees, resulting
each sample radius for both 1985 and 2002. We measured in series spanning the common interval 1928–2005. Cluster-
the proportion of forest in 1975 and then the proportion ing was based on the Euclidean distance between time-series
of pixels that were classified as mature forest in 1975 and (Keogh and Kasetty 2003) using the ‘hclust’ function in R.
remained so in 2000 as our estimate of old forest coverage in The statistical significance of the identified clustering was
2000. Given the constraints of the available datasets, these assessed using distance-based redundancy analysis (dbRDA;
latter variables were the best available surrogates for forest (Legendre and Legendre 2012)) as implemented using the
age at each time period, with clear limitations limiting the ‘adonis’ function in the vegan (Oksanen et al. 2017) package
degree of inference linking forest age to budworm outbreak in R. Here, the matrix of Euclidean distances was used as
dynamics. the response variable and cluster identification was used as a
Forest configuration was quantified using landscape pat- predictor.
tern metrics that capture the spatial connectivity of bud- To determine the degree of periodicity and the primary fre-
worm host patches as well as the amount and size of forest quency of oscillation in outbreak chronologies we conducted
openings (i.e. forest fragmentation). Connectivity of fir and spectral analyses on the cluster-level time-series (function
spruce-fir patches was calculated for each sample radius ‘spec.pgram’ in R; Daniell kernel; smoothing parameter = 3).
based on basal area for each of the two map dates (1985 and
2002) using the COHESION landscape metric (Schumaker Spatial synchrony
1996). Forest openings (non-forest) were represented using We characterized the spatial synchrony of site-level outbreak
‘disturbed’ and ‘open’ land classifications. ‘Disturbed’ land time-series (n = 51) using spatial nonparametric covariance
included pixels that were classified as mature forest dur- functions (SNCF), as implemented in the ncf package in R
ing one 5-yr interval, but were then classified as open, (Bjørnstad 2009). SNCF measures the covariance amongst
grass, brush, or regenerating in the next year interval pairs of sites close together versus those further apart, and
(Sturtevant et al. 2014), representing the disturbance rate. typically produces smooth decay patterns for pairs of sites
‘Open’ land combined all non-forest land cover types (i.e. further apart. SNCFs were computed for all pairs of sites
ignoring water and open wetlands) (Sturtevant et al. 2014) within each management zone (M, W, O), and for all
with disturbed cover types to estimate cumulative fragmen- sites within the study area (i.e. a global SNCF). We used
tation of forest. We calculated the area-weighted mean patch site level data instead of subareas despite relatively fewer
size (AWMPS) and proportion of area for both disturbed replicates per site (5–15 trees) to maximize the number of
and open categories across all sample radii averaged over pairwise comparisons within management zones. We pre-
all 5-yr intervals between 1975 and 2000. AWMPS and dicted that, for any given distance class, there would be
1561
greater synchrony (i.e. covariance) within the less-disturbed, Data deposition
more host-rich Wilderness (W) than in the managed forest
landscapes (M, O). Data available from the Dryad Digital Repository:  http://
dx.doi.org/10.5061/dryad.cj9sp  (Robert et al. 2017).
Multivariate analysis
Constrained ordination (RDA) was used to model how for-
est landscape structure and climate affect variation in subarea Results
outbreak chronologies. Our multivariate response matrix
consisted of the annual % of trees affected by budworm out- Outbreak temporal patterns
break, 1928–2005, for each subarea (rows; n = 16; Fig. 1).
The overall chronology contained four complete, major out-
Potential predictor variables (total n = 101; Table 2), were of
break cycles between 1850 and 2005, yielding a region-wide
three types: 1) forest configuration (n = 48); 2) forest compo-
estimate of cycle periodicity of 39 yr (Fig. 2). The first (peak
sition (n = 48); and 3) climate (n = 5). A greater number of
~ 1872, cycle I) and last (peak ~ 1991, cycle IV) of these four
predictors are included for forest configuration and composi-
cycles were intense, and synchronous, affecting more than
tion due the multiple spatial scales at which these predictors
65% of host trees at the cycle peak. While the first outbreak
were calculated, and in some cases multiple dates (Table 2).
(cycle I) was not used in the analysis because it was defined
Prior to analysis, all predictors were standardized to zero mean
by very few sample trees, an outbreak at this time is consis-
and unit variance. We first identified a parsimonious subset of
tent with previous studies (Blais 1954). The second cycle (II)
predictors within each predictor type using forward selection
came in discrete pulses with variable timing (peaks ~ 1910,
and adjusted R2 as the selection criterion (Dray et al. 2007,
1916, and 1931, spanning a long interval, but never rising
Oksanen et al. 2017). Forward selection was implemented
to affect more than 40% of trees in any one year). The third
using the ‘forward.sel’ function in the ‘packfor’ package in R.
cycle (III; peak ~ 1961), was more synchronous in its timing,
Once predictor subsets were identified for each type, we
but affected just over 40% of trees at the cycle peak.
combined them into a final ordination model. We further
reduced the retained predictors by removing those that were
Clustering of outbreak time-series
identified as collinear based on the variance inflation factor
(VIF). Collinearity among predictors can result in biased The 16 subarea chronologies (see Supplementary material
models and can confound coefficient direction and interpre- Appendix 2 for raw subarea-level chronologies) grouped into
tation. Predictors with a VIF greater than 10 were removed four significant clusters on the basis of raw outbreak chronol-
from the model (Zuur et al. 2009). VIFs were calculated using ogies of percent trees affected (Fig. 3A) (dbRDA; F = 2.44;
the ‘vif.cca’ function in the vegan package in R (Borcard et al. p  0.05). The clusters corresponded largely to management
2011). The relative importance of each predictor type (i.e. zone, although they were also oriented along a southwest to
forest composition, configuration, and climate) was assessed northeast axis (Fig. 3B). Cluster 1 was comprised of a sin-
using variance partitioning (Peres-Neto et al. 2006). Variance gle subarea in the central part of the study area, in north-
partitioning was undertaken using the ‘varpart’ function in ern Minnesota, where just two, long-lasting outbreaks were
the vegan package in R (Oksanen et al. 2017). observed (Fig. 4C). Cluster 2, restricted to the southwestern
Our analysis is the first to relate forest landscape struc- part of the study area, and entirely in Minnesota, exhibited
ture to budworm outbreak dynamics over multiple out- irregular, high-frequency, low-amplitude behaviour that was
breaks. This objective presents two important challenges. poorly synchronized within-cluster (Fig. 4D). Clusters 3 and
First, individual outbreaks may respond more strongly
to forest landscape structure at the time of each outbreak % affected trees of
than the aggregate response of multiple outbreaks to either Border Lakes Master Chronology
100 350
average structure or structure at any specific point in time. Outbreak
I IV 300
Sample II
Second, data on forest landscape structure becomes more
% affected trees
80
Sample depth
Depth 250
limited as one goes back in time; the most comprehensive III
60 200
remote sensing data is associated with the most recent bud- 150
40
worm cycle. To better understand how outbreak timing 100
and the availability of spatial data affect our ability to infer 20
50
outbreak drivers, we repeated the RDA procedure outlined 0 0
above for three shorter windows of roughly equal length, 1850 1900 1950 2000
each spanning one outbreak cycle: 1928–1957 (cycle II, Year
28 yr), 1958–1983 (cycle III, 26 yr), 1984–2005 (cycle IV, Figure 2. Master chronology for the full BLL study area. The grey
22 yr). According to our logic above, we predicted that the section was excluded from formal analyses due to low sample size at
variance explained by the predictors applied to the latest the level of the sample unit (i.e. subarea). Sample depth (red line)
cycle (i.e. IV) would be greater than either the aggregate, refers to the number of trees upon which the disturbance chronol-
multi-cycle analysis or disaggregated analyses applied to ogy (percent affected trees; black line) is based. Cycle peaks associ-
either of the earlier cycles. ated with four major cycles (I–IV) are indicated with arrows.
1562
(A) Cluster dendrogram curves), or within the three management zones (Fig. 5A, B,
C, black curves). As expected, spatial covariance within any
MCN
given distance class was higher (i.e. more synchronous) in
300
the less-disturbed, more host-rich Wilderness (Fig. 5B) than
Height
MCS in the other two areas (Fig. 5A, C). This was especially true
MWS
MWN
OES
100 for distances greater than 40 km, where covariance actually
WNN
WES
OEN
MEN
OWN
WEN
OWS
increased with distance in the Wilderness (Fig. 5B).
WWS
WWN
WSS
0
MES
1 2 3 4 Forest landscape structure and climate influence on
(B) outbreak dynamics
0 25 50 100 Km
The final RDA model applied to the full chronology (1928–
OWN
2005) contained seven variables distributed among the three
predictor types (Table 3). Forest configuration was included
OEN
OWS WNN as the proportion of forest disturbance at scales of 30 and
WWN 5 km. Forest composition was included as the basal area of
WEN
OES balsam fir at scales of 5 and 1 km in the year 2005, and the
MWN
MEN proportion of old forest in 2005 at 10 and 20 km. Climate
WWS WES
MWS MCN
was included as the budworm growth index (Fig. 6).
WSS
MES Lake In total, this model captured 40% of the variance in out-
Superior break pattern among the 16 subareas over 76 yr. 20% of the
MCS
Legend total variance was attributable to forest configuration and
Wilderness
Minnesota forest composition alone, representing roughly half of the
Ontario
Border Lakes
variance explained. The other half of the explained variance
Ecoregion was shared among predictors, with half of that (13%) shared
between forest composition or configuration, and most of
Figure 3. (A) Results of the cluster analysis (n = 4) of the percentage
the other half (9%) shared between forest composition and
of trees affected by defoliation in the 16 subareas. (B) The four
clusters overlaid in geographical space. See Table 1 for subarea climate. Climate alone accounted for none of the variation
acronyms. (Fig. 6).
In analysing each outbreak cycle independently (cycles II–
4 each experienced three periodic outbreaks, with relatively IV; Fig. 2), we found that, as expected, our models captured
high synchrony (narrow confidence interval) in cluster 3 and a greater proportion of variance as we moved forward in time.
lower synchrony (wider confidence interval) in cluster 4 (Fig. During cycle II (1928–1957; Fig. 7A) 14% of the variance in
4A, B). Clusters 3 and 4 were slightly out-of-phase, with outbreak pattern could be explained by forest configuration;
cluster 3 cycles consistently occurring three years later than no other predictors were retained. During cycle III (1958–
cluster 4 cycles, as evidenced in a cross-correlation tabula- 1983; Fig. 7B) 46% of the variance in outbreak pattern was
tion conducted a posteriori (rlag 3 = 0.891). This minor shift explained by forest composition and configuration. During
in phasing simply indicates some potential directionality in cycle IV (1984–2005; Fig. 7C) 66% of the variance in out-
the expansion of outbreaks (e.g. travelling waves) that cannot break pattern was explained by a combination of all three
be further evaluated without analysis of spatial dynamics at classes of predictors, with 40% of the variance explained
scales larger than the study area. by forest landscape structure (inclusive) and the remainder
Spectral analyses applied to the mean chronologies for shared between climate and forest landscape structure vari-
the clusters showed that clusters 1, 3, and 4 exhibited a sin- ables, with none of the variance attributed to climate alone.
gle spectral peak at 40, 33, and 33 yr, respectively (Fig. 4).
Cluster 2, in contrast, exhibited non-stationary behavior,
including two major cycles at the start and end of the series Discussion
(in the 1930s and 1990s), and four minor cycles occurring at
regular intervals in between the major cycles. This resulted in Our results indicate that historical outbreak dynamics were
weak periodicity in the 12–15 yr frequency range. Over the affected by differences in forest landscape structure that
interval 1940–1990 the % of trees affected in cluster 2 never resulted from different land management legacies in the BLL.
rose higher than 40%. We found that the greater the basal area of host tree species
and the proportion of older forest on the landscape, the lower
Outbreak spatial synchrony the frequency of outbreaks and the higher their amplitude,
resulting in outbreaks that are more intense, more extensive,
Spatial covariance in outbreak occurrence generally decreased of longer duration, and of greater synchrony.
with increasing geographic distance between sites, whether We did not expect an important effect of climate on
measured at the scale of the whole study area (Fig. 5, red outbreak behavior because we predicted that climate was
1563
Cluster 4 (OWN,WES,OES,OEN,MEN)
100
80
60
40
20
Cluster 3 (OWS,WEN,WNN,MES,WSS,WWN,WWS)
100
80
60
40
20
% affected trees
0
Cluster 1 (MCN)
100
80
60
40
20
0
Cluster 2 (MCS,MWN,MWS)
100
80
60
40
20
0
1940 1960 1980 2000
Year
Figure 4. Combined chronologies for each of the four clusters shown in Fig. 3, listed from the southwest (bottom) to the northeast (top).
Black = average percentage of trees affected for the cluster. Dashed gray = 95% confidence intervals for the cluster. Red = average percentage
of trees affected for the entire study area (Fig. 1). Dotted horizontal line indicates the 25% threshold defining outbreaks in prior study
(Robert et al. 2012) – but this line is used for reference only and not applied in any analysis. Difference between red and black curves shows
the deviation of a cluster from the broader BLL temporal pattern. See Table 1 for subarea acronyms, and Supplementary material Appendix
2 for raw subarea-level chronologies. Spectral peak from spectral analysis is the following for each cluster: cluster 4 and 3: 33 yr, cluster 2:
12–15 yr, cluster 1: 44 yr.
not sufficiently variable within our study area. However, higher population growth rates, whereas progressively interior
we learned that climate does vary somewhat with distance sites and sites south of our study area may be too warm for
from Lake Superior such that subareas that are closer to the optimal reproduction (Régnière et al. 2012) (Supplementary
lake have cooler springs and summers (Fig. 7C). Indeed, the material Appendix 1 Fig. A1C). Although some of the varia-
spruce budworm population growth index indicates that tion in outbreak histories was due to climate, this variance
the cooler climate closest to Lake Superior could produce was consistently shared with, and therefore indistinguishable
1564
(A) 1.0 Cycle frequency, amplitude, and synchrony
Ontario
0.8 We predicted that outbreak cycle frequencies, amplitudes,
0.6 and synchrony would be related to one another in a sys-
0.4
tematic way, and that this multi-variate response would be
related to forest landscape structure, and this is precisely what
0.2
we observed. Outbreak behavior within clusters 3 and 4 in
0 Wilderness and in Ontario was consistent with most other
(B) 1.0 tree-ring based studies of spruce budworm outbreaks: out-
Wilderness
0.8 breaks were periodic, and evenly spaced, with 33 yr between
Correlation
0.6
cycle peaks (Fig. 4A, B). They were well-synchronized at the
scale of both clusters (Fig. 5A, B), which were extensive, span-
0.4 ning 75% of the subareas (Fig. 3B: 7 subareas in cluster 3, 5
0.2 subareas in cluster 4). Such strong periodicity and synchrony
0 are also consistent with Royama’s (1984) theory of harmonic
(C) 1.0 oscillations in the spruce budworm system.
Minnesota The contrast in observed outbreak behavior between clus-
0.8
ters 1 and 2 was striking even though both clusters are in
0.6 Minnesota. Here, a greater diversity of cycling behavior was
0.4 observed relative to the other zones. We found evidence for
0.2
high-amplitude, low-frequency, 40-yr cycles in tiny cluster
1 (Fig. 4C: 1 subarea), and low-amplitude, high-frequency
0
cycling in cluster 2 (Fig. 4D: 3 subareas). It is possible that
0 20 40 60 80 our statistical test of climatic effects on outbreak patterning
Distance (km) was weak, and that the anomalously low amplitude and high
Figure 5. Spatial non parametric covariance functions (SNCF) for frequency of irregular cycling within Minnesota’s cluster 2 is
the different management zones in the BLL. Each management consistent with the environment being too warm to promote
zone (black) is represented individually on separate graph (A, B, C) regular, high-amplitude cycling. However, this would not
with their 95% confidence interval (dashed lines) and overlaid with explain the opposite pattern being observed in Minnesota’s
the SNCF calculated for the entire BLL (red). adjacent cluster 1.
What might explain these divergent patterns within a
single management zone is spatial variation in forest frag-
from, that due to forest landscape structure. Thus, at the scale mentation of host forests in Minnesota, including the
of our uniquely structured study area, it appears that forest decreased abundance and contagion of host cover. This in
landscape structure is more important than climate in shap- turn would be wholly consistent with the growing number of
ing of spruce budworm outbreak dynamics. results that have emerged in the last two decades supporting
Table 3. Summary of variables selected for inclusion in the final RDA models at each outbreak time window. Multiple time periods and/or
spatial scales selected for a given variable indicated by multiple values separated by the ‘/’ character.
Temporal window Variable family Variable selected Time period Spatial scale
Full Composition Balsam fir basal area 2002/2002 1 km/5 km
Proportion of old forest 2000/2000 10 km/20 km
Configuration Proportion of forest disturbance Time-series average 5 km/30 km
Climate Spruce budworm growth index Time-series average Subarea average
1928–1957 Configuration Area-weighted patch size of forest disturbance Time-series average Subarea average
1958–1983 Composition Spruce and fir basal area 2002 30 km
Proportion of forest 1975 0.5 km/5km
Proportion of old Forest 2000 5 km
Configuration Proportion of forest disturbance Time-series average 5 km
1984–2005 Composition Fir basal area 2002 1 km
Proportion of forest 1975 0.5 km
Proportion of old forest 2000 30 km
Configuration Proportion of forest disturbance Time-series average 1 km/30 km
COHESION of fir patches Time-series average 1 km
Climate Min. spring temperature Time-series average Subarea average
Min. summer temperature Time-series average Subarea average
1565
Configuration Climate
1928-2005 MCN
.05 .02
15 (Full Chronology)
PDist30km .13
OWN .09
10
WES PDist5km .15
OES
RDA2
5
OEN Composition
ClimInd WNN OWS Residual R2 = 0.60
0 MWS
Fir1km2002 MEN MCS
WWNMES
WEN
−5 Fir5km2002
WWS
WSS
POld20km2000
−10 POld10km2000 MWN
−30 −20 −10 0 10 20 30

RDA1
Figure 6. Redundancy analysis (RDA) and variance partitioning relating dissimilarity in the annual percent trees affected to independent
variables chosen by forward selection for each of our three hypotheses: 1) forest configuration; 2) forest composition; and 3) climate. See
Table 1 for subarea acronyms in black font. Inset bubble diagrams indicate the independent and shared variance explained by each class of
variables (see Table 2 for predictor variable descriptions, where the color label corresponds to the class of variables in the bubble diagram).
The RDA is applied to the full chronology is contrasted with separate RDAs applied to each of the three major outbreak periods (Fig. 7).
the ‘silvicultural hypothesis’, which suggests that outbreak Pearson correlations between balsam fir basal area and decid-
dynamics may be stabilized (i.e. lower amplitude and higher uous basal area were low (r = 0.08 and r = 0.34 at 1 and 5 km
frequency) through the removal of susceptible host forest radius, respectively). It is possible that balsam fir is the limit-
from the landscape. ing factor in this region, relative to the geographic locations
The contrast in outbreak synchrony among management where negative relationships between hardwood content and
zones indicates that forest management legacies influence budworm damage have been documented (i.e. southeastern
spatial synchrony (Fig. 5). This result is interesting given the Canada; (Cappuccino et al. 1998, Campbell et al. 2008)).
known long-range dispersal capacity of the spruce budworm Further research will be required to identify which of the
(Greenbank et al. 1980, Anderson and Sturtevant 2011). mechanisms are responsible for the interesting asynchronous
Our study purposefully selected contrasts in forest manage- high-frequency oscillations observed in cluster 2.
ment legacies at spatial scales large enough to overcome the
large scale spatial correlation of outbreak dynamics created Contrast with other studies
by the cycle synchronization process. West and central areas
of Minnesota that contained the least amount of budworm The patterns we report compare well to other reported pat-
host exhibited the lowest degree of synchrony of the three terns in the spruce budworm dendroecology literature, but
management zones (i.e. cluster 2; Fig. 4). In contrast, the are perhaps more directly interpretable because we specifi-
Wilderness, where host was most abundant and forests were cally chose a spatially structured landscape with clear con-
least disturbed, exhibited the highest degree of synchrony trasts forest landscape structure associated with divergent
(Fig. 4–6). Ontario, which had similar rates of disturbance land management legacies (Sturtevant et al. 2014). On the
as Minnesota but where disturbances were more aggregated one hand we observed a relatively stable cycle that was fairly
in space, was intermediate in its degree of synchrony (Fig. 5). well synchronized throughout the bulk of the study area,
Notably, outbreak dynamics of subareas from Minnesota that which is consistent with (Boulanger et al. 2012), in their
contained more host (i.e. eastern Minnesota; Fig. 6) tended small-scale, multi-century study in southern Quebec. On the
to follow the dynamics of the broader study area (i.e. clusters other hand, we also observed persistent asynchrony at finer
3 and 4; Fig. 3) – demonstrating the importance of treat- temporal and spatial scales, as did Jardon et al. (2003), in
ing forest landscape structure as a continuous variable that their large-scale study from across Quebec. We found that
does not necessarily follow political boundaries (James et al. spatial synchrony varied the most among management zones
2011). with different forest landscape structures. Our study sup-
In our study, the concentration and connectivity of ports Jardon et al. (2003)’s assertion that a study area needs
preferred budworm host, rather than deciduous content, to be extensive to reliably estimate outbreak frequency. Our
emerged as the important compositional factors in our RDA observation of a wide range of outbreak intensities (Fig. 2, 4)
analyses (Fig. 6). This result suggests that it is the effect of is further consistent with a small-scale, multi-century study
host forest on dispersal success, and not the effect of non- in Maine (Fraver et al. 2007) that described a relatively low
host forest on natural enemy communities, through which outbreak frequency (~67 yr), but with alternating major and
cycling behavior is modified. While one might speculate that minor cycles of defoliation intensity, where only the major
these factors could be negatively correlated with each other, cycles were counted as ‘outbreaks’.
1566
(A) MCN Configuration Climate
15 1928-1957
.14
10
5 OWS
WESOWN
MWS WWN
0 MES WWS
WNN
AWMptch1km WSSWEN Composition
MEN OEN Residual R2 = 0.86
OES
MWN
−5
MCS
−30 −20 −10 0 10 20
(B) MESPFor0.5km1975 Configuration Climate

1958-1983 PFor5km1975 WENHst30km2002
5 WNN
POld5km2000 WSS
WWS WWNWES OES
MEN OEN
0 .15
RDA2
OWS .31
MCN
OWN
MCS PDist5km
−5 Composition
MWS Residual R2 = 0.59
−10
MWN
−20 −10 0 10 20
(C) 10 OWN
OES Configuration Climate
1984-2005
PDist30km
WES .10
5
OEN MCN PDist1km .16
WNN .34 .03
PFor0.5km1975 MCS
MWS
0 OWSSumTmp .06
Fir1km2002 WEN SpTmp
MEN Composition
−5 MES MWN Residual R2 = 0.34
CoheFir1km1985 WWN
WSS
POld30km2000
WWS
−20 −10 0 10 20
RDA1
Figure 7. Redundancy analysis (RDA) and variance partitioning relating dissimilarity in the annual percent trees affected to independent
variables chosen by forward selection for each of our three hypotheses: 1) forest configuration; 2) forest composition; and 3) climate. See
Table 1 for subarea acronyms in black font. Inset bubble diagrams indicate the independent and shared variance explained by each class of
variables (see Table 2 for predictor variable descriptions, where the color label corresponds to the class of variables in the bubble diagram).
The RDA is applied to each of the three major outbreak periods (A–C).
Previous investigators have reported on large-scale pat- synchronization can break down due to strong transitions in
terns in outbreak synchrony as a function of distance, and forest landscape structure.
have inferred that the Moran effect is a ubiquitous phenome- Looking beyond the spruce budworm, there is evidence
non in both spruce budworm (Williams and Liebhold 2000) for reciprocal insect–host plant feedback operating in other
and other forest insect systems (Peltonen et al. 2002). More periodic defoliator systems, namely: forest tent caterpillar
recent study of synchrony in the spruce budworm demon- in boreal Canada (Roland 1993, Roland and Taylor 1997);
strates that different causal factors – including bottom-up gypsy moth in the northeastern United States (Haynes et al.
factors such as host tree cone production – can contribute 2009); jack pine budworm in central North America (Nealis
to the Moran effect at different times in the outbreak cycle, and Lomic 1994); larch budmoth in the Swiss alps (John-
and concluded that synchrony is unlikely to be caused by any son et al. 2004); and winter moth in Fennoscandia (Ims et al.
one factor (Bouchard et al. 2017). While their study focused 2004). Collectively, this suggests that the research commu-
on the effect of high frequency fluctuation of factors related nity should consider re-evaluating the idea that population
to food availability and budworm survival on synchrony, our outbreak cycles are driven strictly by reciprocal feedback
study focused on the effect of slow changes in forest com- occurring in either the upper or lower trophic levels of a tri-
position and structure on cycling behavior. Our results sug- trophic (host plant–herbivore–natural enemy) interaction.
gest that in spite of the Moran effect, the process of outbreak Wherever both trophic levels seem to be implicated in cycling
1567
behaviour, there may be some opportunity for manipulating relevant forest landscape structure variables through time.
long-term outbreak dynamics by manipulating landscape Alternatively, a stronger effect of forest landscape structure
structure, which is in accordance with the silvicultural on outbreaks moving forward through time (Fig. 7) might
hypothesis. also be the consequence of increasingly divergent landscape
structure as the legacies of forest management activities accu-
Study limitations mulated over time (Sturtevant et al. 2014). That is, budworm
outbreak behavior may simply be tracking forest landscape
While tree-ring data is a good proxy to quantify the temporal structure as it diverges among management zones through
dynamics of budworm outbreaks within a given location, it is time. Without access to high quality data of population
not a perfect proxy. Indeed, the use of the percentage of trees abundance and landscape structure at earlier time periods, we
with radial growth sufficient to indicate defoliation damage cannot yet distinguish between these two plausible explana-
may not capture some of the fine-scaled variation in populations for the results from the outbreak-specific analyses, and,
tion dynamics that may be occurring at the tree-level. None- indeed, they are not mutually exclusive.
theless, a comparison with aerial defoliation surveys suggests
that tree ring data do capture larger scale variation in tree Conclusions
damage associated with budworm outbreaks (Bouchard et al.
2006). Clearly, intensive population studies provide the We found that spatial variation in forest landscape structure,
best local understanding of the underlying processes that including species composition and spatial configuration
may not be detectable using impact proxies such as tree ring influences spruce budworm outbreak cycle frequency, inten-
and aerial survey data (Royama 1984). There are, however, sity, and synchrony. In areas with reduced fir abundance and
no practical methods available to quantify such fine-scaled connectivity, spruce budworm outbreaks had low intensity
population dynamics at the spatiotemporal scale necessary to and high frequency, and were not well synchronized. The
evaluate landscape-scale feedbacks of the forest on budworm opposite was observed in areas where host was more abun-
outbreaks. dant, as in the relatively undisturbed Wilderness. This result
Our study is the first to explicitly examine effects of forest is consistent with the ‘silvicultural hypothesis’ that increasing
landscape structure on budworm outbreak dynamics across severity of spruce budworm outbreaks over the last century
multiple outbreak cycles. A challenge facing any approach resulted from past management activities that enhanced fir,
such as ours is the absence of detailed historic data on forest in this case by demonstrating the reverse situation is true –
conditions during earlier cycles. The availability of Landsat- i.e. a reduction in landscape-level host mitigates the cycling
derived spatial data extends back only to the early 1970s. behavior of budworm outbreaks. This paradigm of strong for-
Further, the earlier Landsat Multispectral Scanner sensors est-insect reciprocal feedback is generally more accepted for
did not have the spectral resolution necessary to distinguish bark beetle systems in western North America (Whitehead
tree species composition (Moore and Bauer 1990) which was and Russo 2005) and western Europe (Temperli et al. 2013),
enabled later by launch of the Landsat Thematic Mapper sen- but until now has been considered controversial for defolia-
sors in 1982 and 1984. However, we know that small-sized tors such as the spruce budworm (Miller and Rusnock 1993,
and spatially-diffuse clearcut harvest disturbances have been Muzika and Liebhold 2000).
applied consistently in Minnesota since the 1930s (White and Relationships between forest landscape structure and out-
Host 2003, Sturtevant et al. 2014), and that the Wilderness break patterns appear to be intertwined with weak responses
had a lower disturbance rate during the last century than to a climatic gradient, which highlights the challenge of iden-
surrounding managed forests (Heinselman 1996). Forestry tifying the unique contributions of forest landscape structure
operations in northern Ontario have been traditionally on spruce budworm outbreaks. Our results also did not
road-limited and applied to lands leased from the provincial directly support the notion that deciduous content affected
government (Rempel et al. 1997) – so there is no reason to outbreak behavior, a surprising result given the mixed for-
suspect the spatial pattern of mechanized forestry operations est conditions of the BLL. Nevertheless, ongoing land man-
prior to 1975 were markedly different from those observed in agement and consequent land-use changes at large scales has
the 1970s and 1980s. Nonetheless, balsam fir is clearly sensi- the potential to have significant unintended effects on spa-
tive to budworm disturbance in the region (Frelich and Reich tial ecological processes. Continued sustainable stewardship
1995), and it is likely that this host species is more dynamic of forest ecosystems requires that more attention be given to
than the otherwise dominant legacy of forest management spatial legacies in forest landscape structure as these legacies
activities because spruce budworm is simultaneously killing affect insect outbreak dynamics over broad spatial scales.
the fir (James et al. 2011). Our results support the synthetic view that while budworm
We found that the analysis specific to cycle IV (i.e. the outbreaks are generally cyclic, these cycles may be muted in
most recent; Fig. 7) provided a more robust estimate of host-deficient forests. To the extent that the budworm itself,
the forest response than both a) the aggregated analysis of through stand-replacing disturbances, is responsible for sig-
all cycles, and b) analysis of earlier cycles (cycle II and III) nificant changes in forest species composition and forest
(Fig. 6, 7). This greater robustness is likely due to the increased spatial configuration, our study implies that there is a recip-
quantity and quality of remote sensing data used to generate rocal feedback relationship between the budworm and forest
1568
that occurs at the landscape scale, as hypothesized decades Blais, J. 1965. Parasite studies in two residual spruce budworm
ago (Baskerville 1975). This provides considerable qualitative (Choristoneura fumiferana (Clem.)) outbreaks in Quebec.
validation to the earliest models of spruce budworm dynam- – Can. Entomol. 97: 129–136.
ics, which hypothesized that forest renewal was a key process Blais, J. 1983. Trends in the frequency, extent, and severity of spruce
affecting spruce budworm cycling (Jones 1977). Even if, as budworm outbreaks in eastern Canada. – Can. J. For. Res. 13:
539–547.
hypothesized by Royama (1984), reciprocal feedback with Borcard, D. et al. 2011. Canonical ordination. – In:
natural enemies is a primary source of outbreak cycling behav- Borcard, D. et al. (eds), Numerical ecology with R. Springer,
iour, it appears that the budworm-forest reciprocal feedback pp. 153–225.
cannot be ruled out as a significant component of the full Bouchard, M. et al. 2006. Forest dynamics after successive spruce
complexity in outbreak patterning. Additional study in dif- budworm outbreaks in mixedwood forests. – Ecology 87:
ferent ecological contexts is needed to better quantify the 2319–2329.
influence of forest landscape structure on budworm outbreak Bouchard, M. et al. 2017. Bottom-up factors contribute to large-
cycling as well as how forest landscape structure may indi- scale synchrony in spruce budworm populations. – Can. J. For.
rectly affect top-down control by natural enemies. Addition- Res.  https://doi.org/10.1139/cjfr-2017-0051 .
ally, we think there would be enormous value in examining Boulanger, Y. and Arseneault, D. 2004. Spruce budworm outbreaks
in eastern Quebec over the last 450 years. – Can. J. For. Res.
the strength of the insect–host relationship in other defoliator
34: 1035–1043.
systems to determine if reciprocal forest–insect feedback is a Boulanger, Y. et al. 2012. Dendrochronological reconstruction of
universal feature of forest defoliator outbreak systems. spruce budworm (Choristoneura fumiferana) outbreaks in
southern Quebec for the last 400 years. – Can. J. For. Res. 42:
Acknowledgements – We thank Nate Aspelin and François Larouche 1264–1276.
for their help in the field and in the laboratory. We also thank Campbell, E. M. et al. 2008. The severity of budworm-caused
Patrick Tobin, Yan Boulanger, and two anonymous reviewers for growth reductions in balsam fir/spruce stands varies with the
helpful comments on the manuscript. hardwood content of surrounding forest landscapes. – For. Sci.
Funding – This project was funded through a grant from the USDA 54: 195–205.
Cooperative State Research, Education and Extension Service Candau, J.-N. and Fleming, R. A. 2005. Landscape-scale spatial
(CSREES) Managed Ecosystems program (2005-35101-16342) distribution of spruce budworm defoliation in relation to
to BRS, PAT, DK, BJC, and MJF, a grant from the Sustainable bioclimatic conditions. – Can. J. For. Res. 35: 2218–2232.
Forest Management Network of Canada to DK, BRS, and BJC, Candau, J.-N. and Fleming, R. A. 2011. Forecasting the response
and through a grant from the US Endowment for Forest and of spruce budworm defoliation to climate change in Ontario.
Communities ( www.usendowment.org ; project # E2014- – Can. J. For. Res. 41: 1948–1960.
0260) to BRS, BJC, and PMAJ, all of which were used to support Cappuccino, N. et al. 1998. Spruce budworm impact, abundance
LER, and the National Fire Plan supporting BRS. and parasitism rate in a patchy landscape. – Oecologia 114:
236–242.
Chang, W.-Y. et al. 2012. Economic impacts of forest pests: a case
References study of spruce budworm outbreaks and control in New
Brunswick, Canada. – Can. J. For. Res. 42: 490–505.
Ahlgren, C. E. and Ahlgren, I. 2001. Lob trees in the wilderness: Cook, E. R. 1985. A times series analysis approach to tree ring
the human and natural history of the Boundary Waters. – Univ. standardization. – Univ. Arizona Press.
of Minnesota Press. Cooke, B. J. et al. 2007. Insect defoliators as periodic disturbances
Anderson, D. P. and Sturtevant, B. R. 2011. Pattern analysis of in northern forest ecosystems. – In: Johnson, E. A. and
eastern spruce budworm Choristoneura fumiferana dispersal. Miyanishi, K. (eds), Plant disturbance ecology: the process and
– Ecography 34: 488–497. the response. Elsevier Academic Press, pp. 487–525.
Barbour, D. A. 1990. Synchronous fluctuations in spatially Dray, S. et al. 2007. Packfor: forward selection with permutation.
separated populations of cyclic forest insects. – In: Watt, A. – R package ver. 0.0-7.
D. et al. (eds), Population dynamics of forest insects. Intercept, Dymond, C. C. et al. 2010. Future spruce budworm outbreak may
pp. 339–346. create a carbon source in eastern Canadian forests. – Ecosystems
Baskerville, G. L. 1975. Spruce budworm: super silviculturist. 13: 917–931.
– For. Chron. 51: 138–140. Eveleigh, E. S. et al. 2007. Fluctuations in density of an outbreak
Bergeron, Y. et al. 1995. Balsam fir mortality following the last species drive diversity cascades in food webs. – Proc. Natl Acad.
spruce budworm outbreak in northwestern Quebec. – Can. J. Sci. USA 104: 16976–16981.
For. Res. 25: 1375–1384. Fraver, S. et al. 2007. Dendrochronological reconstruction of
Berryman, A. 2002. Population cycles: the case for trophic interac- spruce budworm outbreaks in northern Maine, USA. – Can. J.
tions. – Oxford Univ. Press. For. Res. 37: 523–529.
Bjørnstad, O. 2009. ncf: spatial nonparametric covariance func- Frelich, L. E. and Reich, P. B. 1995. Spatial patterns and succession
tions. – R package ver. 1.1-3,  https://CRAN.R-project.org/ in a Minnesota southern-boreal forest. – Ecol. Monogr. 65:
package=ncf . 325–346.
Blais, J. R. 1954. The recurrence of spruce budworm infestations Gray, D. R. 2008. The relationship between climate and outbreak
in the past century in the Lac Seul area of northwestern Ontario. characteristics of the spruce budworm in eastern Canada.
– Ecology 35: 62–71. – Clim. Change 87: 361–383.
1569
Gray, D. R. 2013. The influence of forest composition and climate MacLean, D. A. 1980. Vulnerability of fir-spruce stands during
on outbreak characteristics of the spruce budworm in eastern uncontrolled spruce budworm outbreaks: a review and discus-
Canada. – Can. J. For. Res. 43: 1181–1195. sion. – For. Chron. 56: 213–221.
Greenbank, D. O. et al. 1980. Spruce budworm (Lepidoptera: McGarigal, K. et al. 2012. FRAGSTATS v4: spatial pattern analysis
Tortricidae) moth flight and dispersal: new understanding from program for categorical and continuous maps. – Univ. of
canopy observations, radar, and aircraft. – Mem. Entomol. Soc. Massachusetts, Amherst.
Can. 112: 1–49. McKenney, D. et al. 2007. Customized spatial climate models for
Hassell, D. C. et al. 1999. A mathematical analysis of Jones’s site Canada. – Frontline Tech. Note i108.
model for spruce budworm infestations. – J. Math. Biol. 38: Miller, A. and Rusnock, P. 1993. The rise and fall of the silvicultural
377–421. hypothesis in spruce budworm (Choristoneura fumiferana) man-
Haynes, K. J. et al. 2009. Spatial analysis of harmonic oscillation agement in eastern Canada. – For. Ecol. Manage. 61: 171–189.
of gypsy moth outbreak intensity. – Oecologia 159: 249–256. Moore, M. M. and Bauer, M. E. 1990. Classification of forest
Heinselman, M. L. 1973. Fire in the virgin forests of the Boundary vegetation in north-central Minnesota using Landsat
Waters Canoe Area, Minnesota. – Quat. Res. 3: 329–382. Multispectral Scanner and Thematic Mapper data. – For. Sci.
Heinselman, M. L. 1996. The boundary waters wilderness 36: 330–342.
ecosystem. – Univ. of Minnesota Press. Moran, P. A. P. 1953. The statistical analysis of the Canadian Lynx
Hennigar, C. R. et al. 2008. Differences in spruce budworm cycle. – Aust. J. Zool. 1: 291–298.
defoliation among balsam fir and white, red, and black spruce. Muzika, R. M. and Liebhold, A. M. 2000. A critique of silvicultural
– For. Sci. 54: 158–166. approaches to managing defoliating insects in North America.
Holling, C. S. 1973. Resilience and stability of ecological systems. – Agric. For. Entomol. 2: 97–105.
– Annu. Rev. Ecol. Syst. 4: 1–23. Nealis, V. G. 2016. Comparative ecology of conifer-feeding spruce
Holmes, R. L. 1992. Dendrochronology program library. budworms (Lepidoptera: Tortricidae). – Can. Entomol. 148:
– Laboratory of Tree-Ring Research, Univ. of Arizona, Tucson, S33–S57.
AZ. Nealis, V. G. and Lomic, P. V. 1994. Host-plant influence on the
Holmes, R. and Swetnam, T. 1996. Program OUTBREAK user population ecology of the jack pine budworm, Choristoneura
manual. – Laboratory of Tree-Ring Research, Univ. of Arizona, pinus (Lepidoptera: Tortricidae). – Ecol. Entomol. 19:
Tucson, AZ. 367–373.
Holmes, R. L. et al. 1986. Quality control of crossdating and Oksanen, J. et al. 2017. Package “vegan”. – R packge ver. 2.4-3,
measuring: a user manual for program COFECHA. – Labora-  https://CRAN.R-project.org/package=vegan .
tory of Tree-Ring Research, Univ. of Arizona, Tucson, AZ. Peltonen, M. et al. 2002. Spatial synchrony in forest insect
Ims, R. A. et al. 2004. Do sub-Arctic winter moth populations in outbreaks: roles of regional stochasticity and dispersal. – Ecology
coastal birch forest exhibit spatially synchronous dynamics? – J. 83: 3120–3129.
Anim. Ecol. 73: 1129–1136. Peres-Neto, P. R. et al. 2006. Variation partitioning of species data
James, P. M. et al. 2011. Two-dimensional wavelet analysis of spruce matrices: estimation and comparison of fractions. – Ecology 87:
budworm host basal area in the Border Lakes landscape. - Ecol. 2614–2625.
Appl. 21: 2197–2209. Price, P. W. et al. 1984. A new ecology: novel approaches to inter-
James, P. et al. 2017. Lagged cumulative spruce budworm defoliation active systems. – Wiley.
affects the risk of fire ignition in Ontario, Canada. – Ecol. Appl. Pureswaran, D. S. et al. 2016. Paradigms in eastern spruce budworm
27: 532–544. (Lepidoptera: Tortricidae) population ecology: a century of
Jardon, Y. et al. 2003. Périodicité et synchronisme des épidémies debate. – Environ. Entomol. 45: 1333–1342.
de la tordeuse des bourgeons de l’épinette au Québec. – Can. Quayle, D. et al. 2003. Forest composition, host-population
J. For. Res. 33: 1947–1961. density, and parasitism of spruce budworm Choristoneura
Johnson, D. M. et al. 2004. Landscape geometry and travelling fumiferana eggs by Trichogramma minutum. – Entomol. Exp.
waves in the larch budmoth: landscape geometry and travelling Appl. 107: 215–227.
waves in the LBM. – Ecol. Lett. 7: 967–974. Raffa, K. F. et al. 2008. Cross-scale drivers of natural disturbances
Jones, D. D. 1977. Catastrophe theory applied to ecological prone to anthropogenic amplification: the dynamics of bark
systems. – Simulation 29: 1–15. beetle eruptions. – Bioscience 58: 501–517.
Keogh, E. and Kasetty, S. 2003. On the need for time series data Régnière, J. and Lysyk, T. 1995. Population dynamics of the spruce
mining benchmarks: a survey and empirical demonstration. budworm, Choristoneura fumiferana. – In: Armstrong, J. A. and
– Data Min. Knowl. Discov. 7: 349–371. Ives, W. G. H. (eds), Forest insect pests in Canada. Canadian
Kurz, W. A. et al. 2008. Mountain pine beetle and forest carbon Forest Service, Ottawa, ON, pp. 95–105.
feedback to climate change. – Nature 452: 987–990. Régnière, J. and Nealis, V. G. 2007. Ecological mechanisms of
Lausch, A. et al. 2013. Spatio-temporal infestation patterns of Ips population change during outbreaks of the spruce budworm.
typographus (L.) in the Bavarian Forest National Park, Germany. – Ecol. Entomol. 32: 461–477.
– Ecol. Indic. 31: 73–81. Régnière, J. et al. 2012. Predicting insect distributions under
Legendre, P. and Legendre, L. F. 2012. Numerical ecology. climate change from physiological responses: spruce budworm
– Elsevier. as an example. – Biol. Invasions 14: 1571–1586.
Ludwig, D. et al. 1978. Qualitative analysis of insect outbreak Rempel, R. S. et al. 1997. Timber-management and natural-
systems: the spruce budworm and forest. – J. Anim. Ecol. 47: disturbance effects on moose habitat: landscape evaluation. – J.
315–332. Wildl. Manage. 61: 517–524.
1570
Robert, L.-E. et al. 2012. Effects of forest management legacies on Temperli, C. et al. 2013. Cross-scale interactions among bark
spruce budworm (Choristoneura fumiferana ) outbreaks. – Can. beetles, climate change, and wind disturbances: a landscape
J. For. Res. 42: 463–475. modeling approach. – Ecol. Monogr. 83: 383–402.
Robert, L.-E. et al. 2017. Data from: Landscape host abundance Turchin, P. 2003. Complex population dynamics: a theoretical/
and configuration regulate periodic outbreak behavior in spruce empirical synthesis. – Princeton Univ. Press.
budworm Choristoneura fumiferana. – Dryad Digital Repository, Volney, W. J. A. and McCullough, D. G. 1994. Jack pine budworm
 http://dx.doi.org/10.5061/dryad.cj9sp . population behaviour in northwestern Wisconsin. – Can. J. For.
Roland, J. 1993. Large-scale forest fragmentation increases the Res. 24: 502–510.
duration of tent caterpillar outbreak. – Oecologia 93: 25–30. Wesołowski, T. and Rowiński, P. 2006. Tree defoliation by winter
Roland, J. 2005. Are the “seeds” of spatial variation in cyclic moth Operophtera brumata L. during an outbreak affected
dynamics apparent in spatially-replicated short time-series? An by structure of forest landscape. – For. Ecol. Manage. 221:
example from the forest tent caterpillar. – Ann. Zool. Fenn. 42: 299–305.
397–407. White, M. A. and Host, G. E. 2003. Changes in disturbance
Roland, J. and Taylor, P. D. 1997. Insect parasitoid species respond to frequency, age and patch structure from pre-Euro-American to
forest structure at different spatial scales. – Nature 386: 710–713. settlement to present in north-central and northeastern
Royama, T. 1984. Population dynamics of the spruce budworm Minnesota. – Minnesota Forest Resources Council.
Choristoneura fumiferana. – Ecol. Monogr. 54: 429–462. White, T. C. 2012. The inadequate environment: nitrogen and
Royama, T. 2005. Moran effect on nonlinear population processes. the abundance of animals. – Springer Science and Business
– Ecol. Monogr. 75: 277–293. Media.
Schumaker, N. H. 1996. Using landscape indices to predict habitat Whitehead, R. J. and Russo, G. L. 2005. Beetle-proofed lodgepole
connectivity. – Ecology 77: 1210–1225. pine stands in interior British Columbia have less damage from
Sturtevant, B. R. et al. 2014. Forest recovery patterns in response mountain pine beetle. – Pacific Forestry Centre.
to divergent disturbance regimes in the Border Lakes region of Williams, D. W. and Liebhold, A. M. 2000. Spatial synchrony of
Minnesota (USA) and Ontario (Canada). – For. Ecol. Manage. spruce budworm outbreaks in eastern North America. – Ecology
313: 199–211. 81: 2753–2766.
Sturtevant, B. R. et al. 2015. Modeling insect disturbance across Wolter, P. et al. 2008. Remote sensing of the distribution and
forested landscapes: insights from the spruce budworm. – In: abundance of host species for spruce budworm in northern
Perera, A. H. et al. (eds), Simulation modeling of forest Minnesota and Ontario. – Remote Sens. Environ. 112:
landscape disturbances. Springer, pp. 93–134. 3971–3982.
Su, Q. et al. 1996. The influence of hardwood content on balsam Wolter, P. T. et al. 2012. Greater Border Lakes Region land cover
fir defoliation by spruce budworm. – Can. J. For. Res. 26: classification and change detection. – Newtown Square, PA,
1620–1628. U.S. Dept of Agriculture, Forest Service, Northern Research
Swetnam, T. W. and Lynch, A. M. 1993. Multicentury, regional- Station,  http://dx.doi.org/10.2737/RDS-2012-0007 .
scale patterns of western spruce budworm outbreaks. – Ecol. Zuur, A. et al. 2009. Mixed effects models and extensions in ecology
Monogr. 63: 399–424. with R. – Springer Science and Business Media.
Supplementary material (Appendix ECOG-03553 at  www.

ecography.org/appendix/ecog-03553 ). Appendix 1–2.
1571

Ecography - 2017 - Robert - Landscape Host Abundance and Configuration Regulate Periodic Outbreak Behavior in Spruce

Uploaded by

Copyright:

Available Formats

Ecography - 2017 - Robert - Landscape Host Abundance and Configuration Regulate Periodic Outbreak Behavior in Spruce

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Ecography - 2017 - Robert - Landscape Host Abundance and Configuration Regulate Periodic Outbreak Behavior in Spruce

Uploaded by

Copyright:

Available Formats

Keywords: spruce budworm, harvest disturbance, landscape ecology, forest management

Landsat MSS and Landsat TM

Landsat MSS and Landsat TM

McKenney et al. 2007

Average of 1980–2000, 5 yr increment

Average of 1980–2000, 5 yr increment

Candidate landscape factors affecting budworm dynamics

1985 and 2002

1985 and 2002

than 1.28 standard deviations from the mean chronology ring

been previously used in spruce budworm outbreak reconstruc-

neault 2004, Bouchard et al. 2006). No other defoliators are

tions over multiple years such that they would be identified

COHESION of spruce-fir patches

as an outbreak. Outbreak detection with the program OUT-

Spruce budworm growth index

BREAK was more efficient when applied solely to host spe-

Min. summer temperature

cies (Bouchard et al. 2006). Non-host chronologies were only

Min. winter temperature

Proportion of old forest

used to validate detected outbreaks by visual comparison of

the growth pattern (not shown). The percentage affected trees

by outbreaks each year at either the site or subarea scales were

used as the response variable for cluster analysis, redundancy

ways, each using separate sets of multi-temporal Landsat

−30 −20 −10 0 10 20 30

(B) MESPFor0.5km1975 Configuration Climate

Supplementary material (Appendix ECOG-03553 at  www.

You might also like