1. Introduction
Highlights:
・ The dynamics of forest tent caterpillar outbreaks in Alberta shifted at the turn of the 20th century from simple, cyclic, synchronous behavior to complex, multi-frequential, asynchronous behavior. This shift coincided with the emergence of a latitudinal gradient in outbreak stability and periodicity reflecting a similar gradient in the periodicity of 20th century winter temperatures.
・ Instrumental weather records in northwestern Canada begin around the start of the 20th century, which is when we observed a shift in forest tent caterpillar outbreak dynamics. There is no evidence of forest structure in northern Alberta having changed suddenly at this time.
・ Asynchronous eruptive behavior in western Canada is not congruent with synchronized cycling reported in eastern Canada. In Alberta, incipient outbreaks of forest tent caterpillar are frequently terminated by Arctic cold air anomalies that tend to be regional in scope.
・ A decrease in polar vorticity at the start of the 20th century has led to increased meridional jet stream flow and more frequent arctic weather anomalies, possibly accounting for an observed loss of synchronous decadal periodicity in outbreak occurrence.
・ Climate change-driven changes in insect disturbance probabilities have profound consequences for forest disturbance risk management under rising uncertainty.
The forest tent caterpillar [Malacosoma disstria (Hbn.)] is a periodic defoliator of broad-leaved trees throughout North America [1] . In Ontario, eastern Canada, outbreaks cycle regularly with a periodicity of 13.00 ± 0.95 yr [2] . Curiously, outbreaks in the prairie provinces of western Canada have exhibited poorly synchronized, aperiodic fluctuations [3] [4] . Recurring outbreaks are thought to be the result of a deterministic host-parasitoid interaction that is subject to stochastic meteorological perturbations [5] [6] [7] [8] . But it is not clear to what extent each of these forces is responsible for the induction of population cyclicity and synchrony and whether they might also be responsible for local departures from perfectly synchronized, cyclic behaviour.
Trembling aspen (Populus tremuloides Michx.) is the principal host-tree-species consumed by forest tent caterpillar larvae [6] and tent caterpillar defoliation is the major factor accounting for large-scale temporal variation in aspen ring-widths, outweighing, by far, any effect of moisture limitations caused by drought [9] . Outbreaks of forest tent caterpillar can thus be inferred from careful measurements of aspen ring-widths [10] .
Using tree-ring records, we reconstructed the long-term history of forest tent caterpillar outbreaks in Alberta, western Canada in order to determine if there was a time in the past when forest tent caterpillar outbreaks were once as cyclic and as synchronous as they are now in Ontario [11] . We sought to identify any temporal shifts and/or spatial gradients in the dynamics of forest tent caterpillar outbreaks that could clearly be associated with similar shifts or gradients in climate. Knowing that winter temperature is a key factor affecting egg survival and influencing population dynamics in northern climates [12] , we were specifically interested in any shifts in dynamics that might have resulted from climate warming as North America emerged from the “Little Ice Age” of the late 19th century. We show that the dynamics of caterpillar outbreaks shifted abruptly at the start of the 20th century, from periodic and synchronous to aperiodic and asynchronous, reflecting similar patterns in modern-day winter temperature variability.
2. Methods
2.1. Sample Collection and Preparation
Two to three trembling aspen (Populus tremuloides Michx.) trees were felled in each of 246 plots. Plots were arranged in a 3 × 3 grid of plot clusters. This grid covered 93,500 km2 of forested land in nine municipalities across central Alberta (Figure 1). A section was taken at the base of each tree. Sections were dried at 40˚C for several days, sanded, and polished. Annual growth rings were resolved using a 10 - 60× magnification stereo-microscope with a variety of light sources. Each ring was dated and measured to the nearest 0.01 mm using an ocular micrometer. Sectional chronologies were cross-dated by examining the cross-correlation among ring width profiles, particularly with regard to high-frequency fluctuations in ring width.
2.2. Outbreak Reconstruction
Tree rings were scored as “normal” or “white”, light-coloured rings being indicative of past defoliation [13] . Annual ring widths were measured along a single radial transect across the stem section. Ring-width chronologies from each stem were detrended to remove non-stationarity due to tree-age effects and local stand effects. Two successive cubic splines with 50% frequency responses of 50 and 25 years were used to remove the low-frequency variation in growth [14] . Detrended ring-width series were then averaged across stands and across municipalities to produce local and municipal-scale mean chronologies. Detailed methods of sample preparation, cross-dating, and data processing are described and discussed elsewhere [10] .
2.3. Patterns of Regionalization
K-means cluster analysis of local aspen ring-width chronologies was used to partition the study area into regions of similar dynamics. With 9 grid points we chose k = 9 clusters to test the hypothesis that outbreak dynamics were more similar within municipalities than between municipalities. To detect any spatial differences in outbreak periodicity spectral analysis was applied to cluster mean chronologies. All analyses were performed using S-plus (Insightful Inc.).
2.4. Nonlinear Stability Analysis
A two-stage method of stochastic, nonlinear, autoregressive modeling and stability analysis [15] was used to determine whether outbreak oscillations over the period 1910-1994 were chaotic. In the first stage of analysis, all nine municipal aspen ring width chronologies spanning the period 1910-1994 were modeled by nonlinear neural net regression. A single hidden-layer model with a variable size
![]()
Figure 1. Location of the study area, in Alberta, Canada. Squares indicate locations of tree-ring plots, arranged in a 3 × 3 grid.
of one to four hidden units (k = 1 to 4) was parameterized using a third-order time-lag (embedding dimension d = 3). With an embedding dimension of three, the use of four units implies 1 + k (d + 2) = 21 model parameters, yielding an over-parameterized model that is highly likely to fit the data, but at an increased risk of falsely detecting chaos. For all chronologies, the single-unit models provided the lowest generalized cross-validation. In the second stage of analysis, the FUNFITS package (https://www.stat.ncsu.edu/information/library/mimeo.archive/ISMS__2289.pdf) was used to analyze the stability of the fitted nonlinear neural net models. Global Lyapunov exponents were estimated for time-series of fitted values derived from single-unit neural net regression models of municipal aspen ring width chronologies. Estimation was based on first-order Taylor expansion. Briefly, these Lyapunov exponents measure the exponential rate of divergence away from an attractor given a perturbation of arbitrarily small size. Positive Lyapunov exponents indicate extreme sensitivity to perturbation and rapid divergence from a trajectory, which is the hallmark of chaos.
2.5. Stationarity of Dynamics
The three longest mean ring-width chronologies were split into two series using 1910 as the pivot year (sensitivity analysis shows that our results are insensitive to choice of pivot year). Spectral analysis was used to examine periodicity in outbreaks before and after 1910. Cross-correlations were computed for pairs of chronologies, before and after 1910. These were compared using a paired t-test to determine if outbreak synchrony had changed over time. Cross-spectral analysis of pairs of chronologies was used to determine if changes in correlations were a result of a loss of coherence in a particular frequency range. Nonlinear stability analysis was performed on the early portions of the three long chronologies, and results here were compared to results from the latter part of the chronologies to determine if stability properties had changed over time.
3. Results
3.1. Outbreak Patterns
Six of the municipal aspen ring-width chronologies originated between 1880 and 1920 and three originated before 1850. The three long chronologies revealed no fewer than 23 episodes of growth suppression that spanned three or more years and were therefore likely due to tent caterpillar defoliation (Online Appendix Supplementary Figure A1). Many of these episodes were evident in all nine municipal chronologies, suggesting periodic outbreaks tended to be synchronized to some extent. Outbreak synchrony was, however, imperfect. Some chronologies were missing some episodes of growth reduction, and some episodes varied in timing among municipalities.
Cluster analysis of all 246 stand-level ring-width chronologies over the period 1936-1994 indicated a latitudinal gradient in dynamics (Figure 2(a)). In the northern region, fluctuations were weakly periodic, with outbreaks occurring at super-decadal intervals (Figure 2(f)). In the southern region, fluctuations were more periodic, with outbreaks occurring at roughly decadal intervals (Figure 2(h)). In the central region, the spectrum of variability was bi-modal, with both super- and sub-decadal periodicity, and a distinct lack of decadal periodicity (Figure 2(g)). Imperfect synchrony of province-wide outbreaks was therefore a result of regional differences in outbreak periodicity. So not only were local outbreaks not phase-synchronized, they also were not homogeneous in frequency.
3.2. Nonlinear Stability
We wanted to know whether the high-frequency pattern of alternating low- and high-amplitude fluctuations exhibited in chronologies from municipalities of the southern region (Figure 2(e)) might be a result of nonlinear, chaotic dynamics [16] [17] . We were especially interested in the possibility that large-scale outbreak asynchrony might be a result of municipal-scale chaos [18] rather than regional stochasticity.
![]()
Figure 2. Results from k-means cluster analysis of stand-level aspen ring width chronologies over the period 1936-1994 (k = 9). Top: (a) cluster numbers assigned to stands, mapped in geographic space; (b) close-up view of Cooking L. municipal grid. Clustering with k = 3 suggests the nine clusters can be aggregated hierarchically into three larger regions constituting a clear latitudinal gradient in FTC outbreak dynamics, shown as three coloured regions in (a). Bottom: Mean chronologies for each of nine clusters derived from k-means cluster analysis plotted in the temporal and frequency domains, and grouped according to latitudinal regions of similarity in (c), (d), (e). Corresponding regional spectra shown in (f), (g), (h). Gray bands on spectra indicate 5- to 10-year range in periodicity. Super-decadal outbreaks in the north give way to sub-decadal outbreaks in the south. Outbreaks in the central region display both low- and high-frequency components.
There was a clear geographic pattern (Figure 3(a)) in the predictability and stability of modeled outbreaks (Figure 3(b), dark circles): northeastern municipalities exhibited fluctuations that were highly stochastic, but highly stable, whereas southwestern municipalities exhibited fluctuations that were less stochastic and
![]()
Figure 3. Geographic patterns in stochasticity and stability of FTC outbreak dynamics at nine municipalities in Alberta: (a) locations of municipalities; (b) R2 values from single-unit neural net regression models plotted against global Lyapunov exponents derived from stability analysis (black circles: 1910-1994 data from nine municipalities; white circles: 1850-1910 data from three exceptionally long chronologies). Dotted line divides municipalities into two groups based on geography (a), and dynamics (b). Over the period 1910-1994 highly stochastic (R2 < 0.3) and stable dynamics (GLE < 0) were prevalent in the northeast, while less stochastic and less stable dynamics were prevalent in the southwest. Outbreak dynamics from 1850-1910 at Peerless L. (P˚) and Conklin (C˚) were no more stable or predictable than after 1910, whereas dynamics at Whitecourt from 1850-1910 (W˚) were more stable and more predictable than after 1910 (see Figure 4).
less stable. None of the chronologies were found to be chaotic, suggesting regional asynchrony in outbreaks was not a result of municipal-scale chaos. The simplicity and clarity of this spatial gradient in stability and stochasticity suggests that the pattern is robust and should therefore have a simple explanation.
3.3. Shift in Dynamics
A visual inspection of the three longest municipal ring-width chronologies (Whitecourt, Peerless L., Conklin) suggested that the dynamics of tent caterpillar outbreaks may have shifted around the turn of the century. Spectral analysis confirmed that, from 1850 to 1910, outbreaks occurred at 6- to 7-year intervals (Figure 4(a)), and from 1910 to 1994, they occurred at both 5-year and 20-year intervals (Figure 4(b)). Cross correlations among pairs of chronologies indicated that outbreak synchrony was higher prior to 1910 (r = 0.59 ± 0.02) than after (r = 0.41 ± 0.02). A paired t-test indicated this difference was significant (one-tailed t = 15.6; d.f. = 2; P = 0.002). Cross-spectral analysis showed that the
![]()
Figure 4. Time-series analysis of long chronologies from three municipalities over two time periods, 1850-1910 (left panels) and 1910-1994 (right panels). Spectra for municipal chronologies shown for (a) early and (b) late period. Cross-spectral coherence among pairs of municipal chronologies shown for early (c, e, g) and late (d, f, h) period, along with correlation (Pearson’s r) between pairs of partial ring-width chronologies. Shaded bars highlight variability in the 5- to 10-year range of periodicity. Spectra smoothed by Daniell smoother with spans = {3, 3}. Dotted lines on coherency estimates indicate 95% confidence intervals.
reduced correlation between pairs of chronologies after 1910 was due to a loss of coherence in the sub-decadal frequency range (Figure 4(c), Figure 4(e), Figure 4(g) vs. Figure 4(d) Figure 4(f) Figure 4(h). Thus a shift from periodic, synchronous behavior to multifrequential, asynchronous behavior appears to have occurred around 1910.
Was this shift in periodicity and synchrony associated with a reduction in stability? As before, nonlinear neural net regression was used to model the long chronologies from Whitecourt, Peerless L., and Conklin, both before and after 1910 (Figure 5). Stability analysis of the early-period models indicated that the global Lyapunov exponents at Peerless L. and Conklin were roughly equivalent before and after 1910 (Figure 3, light vs. dark circles). At Whitecourt, global Lyapunov exponents increased substantially after 1910, indicating a sudden decline in stability, but not to the point of chaos.
![]()
Figure 5. Results from nonlinear neural net regression analysis of early and recent portions of three long chronologies (thick lines: original data, as in supplementary Figure A1; dotted lines: fitted values).
3.4. Regional Environmental Stochasticity
We examined spatial and temporal variability in mean January temperatures in Alberta in a search for patterns that might explain the emergence of regionalized outbreak dynamics during the last century. We discovered that winter climate patterns are regionalized in a remarkably similar manner. In northern Alberta, January temperature fluctuations during the 20th century were dominated by super-decadal variability, whereas in southern Alberta there was an additional sub-decadal component of variability in the 5- to 10-year range of periodicity that was not shared with northern Alberta (Figure 6). Furthermore, super-decadal peaks in mean January temperature in northern Alberta (Figure 6, Fort Vermilion) corresponded with the occurrence of large-scale outbreaks in 1924, 1942, 1962, and 1982 throughout northern Alberta.
![]()
Figure 6. January monthly mean temperature 1913-1985 (left panels) along a latitudinal gradient in Alberta. Corresponding spectra (right panels) smoothed by Daniell smoother with spans = {3, 3}. Arrows indicate periodicity in the 0.1 to 0.2 y-1 frequency range (10 to 5y). Horizontal and vertical lines on periodograms indicate bandwidth and 95% confidence interval, assuming white null spectrum.
4. Discussion
4.1. Robustness of Reconstructions
A major shift in the patterns of trembling aspen ring widths in Alberta appears to have occurred around 1910. The small sample sizes and wide error bars on the long chronologies prior to 1910 are a concern, especially for the Peerless L. data. However two facts are worth noting. First, it is despite the smaller sample sizes that between-site cross-correlations were found to be higher prior to 1910. Largersample sizes prior to 1910 would likely increase the similarity of estimated mean ring widths. Second, though the pre-1910 mean ring width estimates for Peerless L. were subject to exceptionally high sampling error, the between-site cross-correlations involving Peerless L. were no smaller than the cross-correlation that did not involve Peerless L. (Figure 4). Larger sample sizes at Peerless L. would likely decrease error while having little influence on the periodic pattern in ring widths. In summary, it is highly unlikely that sub-decadal coherence among early chronologies is due to random chance alone. It is more likely a result of synchronized, sub-decadal outbreaks of forest tent caterpillar.
The error bars on cross-spectral coherency estimates (Figure 4) were sufficiently wide that coherency differences before and after 1910 were not statistically significant. This, however, only serves to illustrate the limited descriptive power of spectral decomposition when the data are only weakly periodic, as in the 1910-1994 ring widths. The change in outbreak periodicity and synchrony after 1910 was so profound that the effect could not be revealed using a simple cross-spectral test for a change in outbreak frequency.
Additional proof that defoliator outbreaks once occurred with high frequency can be found in the white ring counts (Supplementary Information Figure A1, bars), which are a reliable, but conservative indicator of insect defoliation [13] . Of the ~10 cycles in aspen growth that occurred between 1840 and 1910, the last three (1888-1891, 1896-1898, 1900-1904) were accompanied by high white ring counts (in one, two, and two of the three municipalities respectively). This further supports our contention that similarities among ring-width chronologies prior to 1910 are not spurious.
If climate were the primary factor limiting tree growth, this shift in ring-width patterns would be puzzling. Climate warming over the 20th century has caused trembling aspen in Alberta to become increasingly sensitive to moisture limitation [19] [20] , yet aspen ring widths have become decreasingly correlated. The only logical explanation for the loss of spatial synchrony in ring widths is that some limiting factor other than moisture has changed, from synchronous and periodic to asynchronous and aperiodic. We argue that this shift in ring-width patterns occurring at the start of the 20th century must be due to a shift in the dynamics of forest tent caterpillar outbreaks.
4.2. Insect Population Dynamics and Climate Change
The cause of this shift in outbreak periodicity and outbreak synchrony is unknown. It cannot be a result of anthropogenic changes in forest structure (e.g. [21] ) because it is only within the last 40 years that the forests in and around the three long-chronology municipalities have been disturbed by large-scale forest activity and oil and gas development. It could be a result of climate change because cold winter temperatures can kill tent caterpillar eggs [22] , and can assist in the collapse of localized outbreaks [12] .
The Albertan climate has, on average, been warming since record-keeping began in the 1880’s [19] , so it is not immediately obvious how this might lead to greater episodes of winter mortality of caterpillars during emerging outbreaks. We suggest that the answer to this paradox lies in understanding how changes in polar vorticity under global climate warming may result in more altered jet stream flows, including increasingly meridional flows of wider amplitude, and thus a higher incidence of Arctic weather anomalies during the boreal winter [23] - [30] . An increase in winter weather variability under climate warming might have destabilized caterpillar population dynamics by increasing the frequency of heavy overwintering mortality, thus reducing the degree of periodicity associated with the predator-prey oscillation. This would explain the temporal trend toward lower stability at Whitecourt (the southernmost municipality of the three long chronologies) and the spatial trend toward lower stability in southwestern Alberta (Figure 3).
Additionally, the weak low-frequency variability observed in January temperatures in Alberta could act as a stochastic period-forcing mechanism [31] . This might explain the bidecadal outbreak cycle in northern Alberta and the higher frequency, sub-decadal and decadal outbreaks in southern Alberta (Figure 2). Notably, the super-decadal peaks in mean January temperature in northern Alberta (Figure 6, Fort Vermilion) corresponded closely with the occurrence of large-scale outbreaks in 1924, 1942, 1962, and 1982 throughout northern Alberta, lending this hypothesis some credence. The possibility of stochastic periodic forcing is made even more likely by the fact that the impact of stochastic winter temperature perturbations is not entirely density-independent; caterpillar populations are more likely to be impacted by winter cold when eggs are laid in the shrub layer and in frost pockets, places they tend to avoid at low density [12] .
Of course, if outbreak periodicity were forced by an exogenous agent, this would violate the assumptions of the nonlinear neural net regression models, and would make the results of stability analysis difficult to interpret. However, it wouldn’t make the geographic pattern in municipal model parameters any less meaningful. The geographic gradient in apparent stochasticity and apparent stability, if difficult to interpret, is nonetheless robust, and therefore requires explanation. We conjecture that super-decadal temperature cycles in northern Alberta have disrupted the endogenous sub-decadal outbreak cycle that once occurred there, and that the absence of such low-frequency climatic variability in southern Alberta has allowed the sub-decadal outbreak cycle there to persist. We are currently exploring this hypothesis through experimentation, process-modeling and simulation.
Similar reconstructions of multicentury, regional-scale outbreaks of conifer-feeding Lepidoptera have suggested a shift towards more frequent and more severe outbreaks induced by anthropogenic changes in forest structure [32] [33] . Our results are more in agreement with small-scale reconstructions of Pandora moth outbreaks [34] , which indicate a reduction in outbreak frequency, severity and synchrony occurring in the 20th century. What is unique about our study is that we have implicated an ecological mechanism that has strong empirical precedence. In the prairie provinces of west central Canada outbreaks of forest tent caterpillar are frequently terminated prematurely by winter weather anomalies that result in high overwintering mortality of eggs [12] .
Why do tent caterpillar outbreak dynamics today seem to differ between eastern and western Canada? The answer does not seem to lie in the intrinsic dynamics of the system, as studies from Alberta yield population growth parameters [8] that are very similar to studies from Ontario [11] . This is a system whose intrinsic dynamics, across Canada, are only weakly nonlinear. Our results suggest the difference in extrinsic dynamics may be attributable to a differential change in climate-driven winter weather variability arising at the start of the 20th century. Our hunch is that winter weather in the continental climate of the boreal plains region of Alberta is now far more variable than in Ontario, and thus tends to act as a strong forcing agent on tent caterpillar cycling. In Ontario, winter temperatures are neither as harsh nor as variable, consequently the system is free to cycle and to synchronize.
4.3. Influence of Forest Structure and Relevance to Forest Management
Although it has been demonstrated that spatial variation in forest tent caterpillar outbreak duration across Canada is determined largely by anthropogenic changes in forest landscape structure affecting the configuration of host trembling aspen on the landscape [8] [21] , there is no reason to suppose the species composition of the northwestern boreal forest changed dramatically or suddenly through the 19th and early 20th centuries. For example, a commercial forest industry did not become established in this part of northern Alberta until the late 20th century. Although there were major transitions from grasslands to aspen parkland in southcentral Alberta (resulting from grass fire suppression) and from spruce to aspen in central Alberta (resulting from selective harvesting of conifers), there is no evidence of such drastic changes in northern Alberta, in the areas where our long chronologies originate (Whitecourt, Peerless Lake, Conklin). For this reason, we surmise that the reported early 20th century shift in tent caterpillar dynamics is not attributable to any change in forest structure that might have increased or reduced the proportion of host trembling aspen in the landscape. This leaves us with the important conclusion that climate change―including not just climate warming, but also changes in climate variability―has the potential to significantly alter forest disturbance probabilities through time.
Managers of boreal mixed woods will want to consider carefully how best to manage a forest resource that is subject to probabilistic insect disturbances whose probabilities are changing in uncertain ways, including large and uncertain spatial variations in apparent responses to climate change. In other major forest-pest systems, the risks of established, indigenous disturbance agents are rising at the northern edge of their range, but declining in the southern part of their range. Key examples include the mountain pine beetle (Dendroctonus ponderosae Hopkin) in western North America [35] and the spruce budworm (Choristoneura fumiferana Clem.) in eastern North America [36] [37] . Predicting how quickly insect outbreak ranges will shift in response to ongoing climate change is a major research challenge. Forest managers therefore need to be prepared to manage resources in the face of considerable, even irreducible, uncertainty.
5. Conclusions
The dynamics of forest tent caterpillar populations in Alberta appears to have shifted from cyclic to asynchronous around the turn of the 20th century. The cause of this shift in dynamics is not known with certainty, but the pattern is reminiscent of the sort of change one might expect from an increase in the variability of winter temperatures. The changes observed in the early 20th century in Alberta are consistent with the contrasting observations from the late 20th century between Alberta and Ontario: more variable winter weather leads to asynchronous and aperiodic outbreak dynamics.
A fundamental issue in population ecology is the role of nonlinear feedback in explaining complex dynamic behaviour [38] . Our results suggest periodic forcing also should not be ignored. That the interaction between stochastic and deterministic nonlinear processes can be a rich source of dynamic behaviour is widely recognized by ecologists [39] [40] [41] [42] . However, our study shows such interactions may also be important for understanding the behaviour of linear and weakly nonlinear systems.
This is the first example of a spatially structured population exhibiting a clear and sudden switch from simple to complex dynamic behaviour in response to global climate change. This illustrates why it is so challenging to predict how ecosystems will respond to climate warming. A slow change in climate brings changes in weather and weather variability, and because it is on this short time scale that many ecosystem processes tend to function, climate change may result in rapid and non-intuitive ecosystem responses. Predicting the impacts of climate change on ecosystem function is going to be challenging. For small, well-studied parts of an ecosystem, it may be possible to make predictions, if, as in this study, we support 1) long-term, large-scale, process-oriented studies of well-studied populations in natural environments, 2) continued environmental monitoring, and 3) the development of more sophisticated analytical methods. For whole ecosystems, which include many lesser-studied components, we are much less optimistic. Given the apparent lack of predictability in these disturbance-driven systems, we recommend modern approaches to the management of probabilistic risks.
Acknowledgements
We thank A.R. Ives and an anonymous reviewer for commenting on an earlier draft. Field studies were assisted by J. Brzustowski, L. Louden, J. Litke, K. Sabourin and E. Hunt. Weather data were obtained through the Global Historical Climate Network database maintained by the National Oceanic and Atmospheric Administration.
Statement of Authorship
BC designed the study, performed the research, and wrote the manuscript. JR contributed revisions to an early draft.