#### Data sets

Four different data sets were analyzed. The first set includes fur-counts of the eleven animal populations within a given year. The three other sets were obtained by augmenting the fur records by one climatic record in each case. The three records we use are the North Atlantic Oscillation (NAO) index, as defined in [28], the cold-season Niño-3 sea-surface temperatures (SSTs), and the mean surface-air temperature of the Northern Hemisphere.

The time series of eight animal populations (bear, beaver, fox, lynx, marten, otter, wolf and wolverine) are 98-year long, extending from 1752 to 1849; those of fisher, mink and muskrat start in 1767 and are thus only 83-year long. All the fur-counts are provided in Appendices 1 through 11. The main results reported in this paper were obtained using the eight species with 98-year long records. When adding the three species with shorter records to the previous eight, the results are very similar (not shown). All the climatic data were available throughout the 1752–1849 time interval.

The NAO index represents a suitably normalized difference in sea-level pressure between the Açores High and the Icelandic Low; it determines the strength of the westerly winds over the North Atlantic Ocean and is correlated with several climatic patterns over the adjacent land masses [28]. We calculated the NAO index by using 3-month averages of the monthly sea-level pressure data provided by [29] for Iceland and Gibraltar, from December 1658 on; the monthly data are available in ref. [30] and our 3-month box-car averages, as used in the present study, appear in Appendix 12. This data set is well correlated with other paleoclimatic proxy indicators of NAO [31], as well as with modern instrumental records, when available.

The cold-season Niño-3 SSTs are a reliable index of the phase of the El Niño-Southern Oscillation (ENSO); ENSO dominates interannual climate variability in the Tropical Pacific [32] and has important effects on other parts of the world, including North America [16,17]. The Niño-3 SSTs and the Northern Hemisphere temperature (NHT) index are taken from [33]; the values used in this study are provided in Appendices 13 and 14, respectively.

#### Principal component analysis

Our data analysis includes two steps: principal component (PC) analysis and spectral analysis. PC analysis has been extensively used in the biomedical sciences [34,35], as well as in climate dynamics [36]. Its main purpose here is twofold: (i) data reduction, *i.e.*, separating the signal from the noise in each data set; and (ii) the identification of significant correlations among the time series. We refer to the previously cited publications [34-36] for technical details.

#### Advanced spectral methods

Once the population-count or climatic signal was isolated, we used a battery of spectral analysis tools to determine the main periods it contains. Singular-spectrum analysis (SSA) is based on eigenvalue-eigenvector decomposition of a time series' lag-covariance matrix [37,38]. Given a series of length *N*, and a maximum lag *M*, the eigenvectors are data-adaptive basis functions for the representation of the series and are called empirical orthogonal functions (EOFs), by analogy with the PC analysis of meteorological fields. The eigenvalues are the associated variances λ_{k}, ordered from largest to smallest. When two eigenvalues are nearly equal, and the corresponding pair of (odd and even) EOFs are in phase quadrature, they may capture, subject to statistical significance tests, an anharmonic (*i.e., *not sinusoidal) oscillation of possibly nonlinear origin [8,12].

SSA can decompose a short, noisy time series into a (variable) trend, periodic oscillations, other statistically significant components that are aperiodic, and noise. The projection of the time series onto an EOF yields the corresponding principal component (PC) of length *N*-*M*+1. Reconstructed components (RCs) are series of length *N *that are obtained by the least-square fitting of their lagged copies, at lag 0, 1, ..., *M*-1 to the projection of the original series and its copies onto a given EOF or a set of EOFs [5,12]. In the life sciences, SSA has already been applied to the analysis of continuous zooplankton records and their environment [39], the possible connections between ENSO and cholera [40], biophysical [41] and neurophysiological problems [42], among others.

The multi-taper method (MTM) is designed to reduce the variance of spectral estimates by using a small set of tapers rather than the unique data taper or spectral window used by Blackman-Tukey methods [43]. MTM has the additional advantage of being nonparametric, in that it does not prescribe an *a priori *(*e.g., *autoregressive) model for the process generating the time series under analysis. A set of independent estimates of the power spectrum is computed, by pre-multiplying the data by orthogonal tapers which are constructed to minimize the spectral leakage due to the finite length of the data set. The optimal tapers or "eigentapers" belong to a family of functions known as discrete prolate spheroidal sequences (DPSS) and defined as the eigenvectors of a suitable Rayleigh-Ritz minimization problem [44]. Averaging over this set of spectra yields a better and stabler estimate – *i.e., *one with lower variance – than do single-taper methods.

Mann & Lees [4] and Ghil *et al. *[5] show that the spectral resolution of the MTM method equals Δ*f *= min {*f*_{N}/4, *2pf*_{R}}, where *f*_{N }is the Nyquist frequency, *f*_{R }is the Rayleigh frequency, and *p *is a bandwidth parameter; in the present case *f*_{N }= 0.5 cycles.year^{-1}, *f*_{R }= 0.02 cycles.year^{-1}, and we used *p *= 2. This means that Δ*f *= 0.04 cycles.year^{-1}, and so all the periods listed in Table 2 are well separated by the MTM method. The Monte Carlo SSA method only yields a spectral resolution of 1/*M *= 0.11 cycles.year^{-1}, which is marginally useful for separating the shortest periodicities in the Table, but the fact that the peaks coincide in both methods lends further credence to their independent existence.

The SSA-MTM Toolkit was developed by the Theoretical Climate Dynamics group at the University of California, Los Angeles (UCLA) [45] and further improved in collaboration with researchers in Europe and North America [5,46]. It supports four different spectral methodologies: classical Fourier analysis, SSA, MTM, and the maximum entropy method (MEM). The toolkit provides a battery of statistical significance tests for each method, as well as visualisation tools that help compare results between the methods. More complete descriptions of all four methods, as well as comparisons of their features and practical performance can be found in [5]. For the data sets under examination, we found SSA and MTM to be the most useful. The entire SSA-MTM Toolkit, along with the User Guide, is available as freeware [46]; further references on advanced spectral methods and their various applications are also listed there.

#### Heuristic population model

We propose here to link the striking and fairly sudden increase in variance of all the fur-count records to a slight variation of environmental conditions. To do so, consider the very simple model:

of the interaction of a prey population *x *with its predator population *y*. The parameters *r *and *s *represent the birth rate of prey and predator respectively, while α stands for the carrying capacity. The predation function is of Holling type II, *h *is the rate of mortality of the predator, and *k *is the rate of predation exerted on the prey.

All parameters are strictly positive, except *C*, which stands for the effects, negative or positive, of external pressures on the prey population; *C *is zero if *x *= 0 but is otherwise assumed to be independent of *x*. Our model (1) is close to harvesting models used in fisheries management [47-49]. If the hunting pressures imposed on the prey population are too large, then *C *may be small or negative, even if other abiotic parameters are favorable.

Non-dimensionalizing the model (1) gives:

where and τ = *rt*.

The equilibrium associated with system (2) was studied using the CONTENT software [50]. Let us suppose that hunting pressures increased, as it is likely that they did during the period 1752–1848, and so the parameter *e *decreased. To examine the effects of such a decrease, we performed a bifurcation analysis with respect to the parameter *e*. An oscillatory instability occurs at *e *= *e*_{H}: if *e *>*e*_{h}, a single equilibrium, with finite and nonzero *u *and *v*, is stable; for *e *≤ *e*_{h}, this equilibrium becomes unstable and gives rise to a stable limit cycle. The larger the difference *e *- *e*_{h}, the larger the amplitude of the oscillations (see Figure 4). The amplitude is roughly proportional to the square root of the difference *e*_{h }- *e*, as expected in the vicinity of a Hopf bifurcation.

This very simple model, with only one prey and one predator species, shows that, as the external conditions deteriorate even slightly, the amplitude of both the predator's and the prey's population cycles increases. The large change in the oscillations' amplitude for all the population records studied here could therefore be linked to a deterioration in Canada's environmental conditions in the early 1800s. This deterioration might also be linked to increased hunting pressures.

Other results show that harvesting pressure may have a destabilizing effect. Basson & Fogarty [10] have demonstrated that harvesting of adults in an age-structured model may have a destabilizing effect in the sense that complex dynamics may then emerge. They also explained the potential role of predation links in the destabilization; this role is in accordance, at least intuitively, with the presence of multiple prey and predator populations in our records. Gamarra & Sole [11] showed that lynx trapping may be the reason for a switch between a stable regime, with low-amplitude fluctuations of the population, to a less stable regime, with larger fluctuations and more complex dynamics. Our model, while simpler than theirs, still seems to capture the key aspect of this effect.