{'fraction_non_alphanumeric': 0.04386864261969118, 'fraction_numerical': 0.01885854537546214, 'mean_word_length': 4.790250044972117, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Seasonality, density dependence and spatial population synchrony

Spatial population synchrony has been the focus of theoretical and empirical studies for decades, in the hopes of understanding mechanisms and interactions driving ecological dynamics. In many systems, it is well-known that seasonality plays a critical role in the density dependence structure of the populations, yet this has hardly received any attention in synchrony studies. Here, we propose a protocol that allows to elucidate deterministic and stochastic sources of spatial synchrony, while accounting for geographic- and season-specific density dependence. We apply our protocol to seasonally-sampled time series of sub-arctic gray-sided voles, known for marked spatial synchrony. Dissociating seasonal density-dependence contributions to the total observed synchrony reveals differential strength and shape of synchrony patterns by season. Mild winter weather reveals to be an important driver of vole spatial synchrony, with lagged effects in the fall. This has direct implications to the future population dynamics of such species when facing climate change.


Introduction
Spatial synchrony -referring to the extent local populations display simultaneous changes across space -is a universal characteristic of geographically distributed populations. The strength and scale of population synchrony, which varies tremendously between species and ecosystems, has been the subject of a large number of theoretical and empirical studies (reviewed by Liebhold et al. (2004); Hansen et al. (2020)). These studies are motivated by their potential to provide unique insights into the mechanisms that drive ecological dynamics across a range of spatial scales (Bjørnstad et al., 1999;Koenig, 1999;Walter et al., 2017). The study of spatial population synchrony is one of the fields within ecology that is, both conceptually and methodologically, most tightly linked to other sciences that also deal with spatio-temporal dynamics (Nareddy et al., 2020;Pérez-García et al., 2021). P.A.P. Moran (1953) developed the first formal theory of spatial population synchrony. Moran's theorem postulates that populations subjected to the same regulatory biotic mechanisms (i.e. log-linear density dependence), and influenced by the same (or perfectly correlated) abiotic environmental variation (e.g. stochastic weather), will display a synchrony that mirrors the synchrony of the environmental variation (Moran, 1953;Hudson and Cattadori, 1999;Hansen et al., 2020). While this theorem has become a cornerstone of the study of population synchrony, Moran himself expressed the need for relaxing some of its restrictive assumptions in order to be more applicable to empirical case studies. Subsequently, many studies have contributed to a "generalization of the Moran effect" (sensu Hansen et al. (2020)) by, for instance, allowing for non-linear density dependence (Blasius et al., 1999;Engen and Saether, 2005), spatially heterogeneous (Royama, 2005;Hugueny, 2006) and temporally autocorrelated environmental variation (Massie et al., 2015), and inclusion of other synchronizing mechanisms (e.g. dispersal Ripa (2000) and trophic interactions Jarillo et al. (2020)). Analytical approaches to elucidate the effect of climatic variation on population synchrony have become particularly timely in the current era of anthropogenic climate change (Sheppard et al., 2015;Koenig and Liebhold, 2016).
Accounting for seasonality was a fundamental aspect highlighted by Moran (1953) when assessing the effect of meteorological conditions on population synchrony. This became clear to him when analyzing population time series of lynx from boreal Canada, which is a region with strikingly different climate in summer and winter. Moran realized that season-specific biotic mechanisms were important, because different demographic parameters are involved in the two seasons (e.g. reproduction only in summer). However, because the lynx population time series were based on only one census per year, Moran was not able to analytically account for seasonspecific population processes (i.e. density dependence). More modern studies of seasonallysampled boreal and arctic rodent populations have shown that marked season-specific density dependence is indeed present and a crucially important determinant of local population dynamics (Hansen et al., 1999;Stenseth et al., 2003;Fauteux et al., 2021). Although seasonality is such a critical aspect of most ecological systems (White and Hastings, 2020), and changing seasonality is one of the most profound consequences of global warming in the northern hemisphere (Xu et al., 2013), we are not aware of any study of population synchrony that has explicitly incorporated seasonality. On the contrary, it has even been argued that one should use yearly averages to remove the influence of seasonality when estimating synchrony (e.g. Dallas et al. (2020)).
The purpose of the present study is to devise a general stepwise analytical protocol, that accounts for season-dependent and geographic context-dependent population processes, to identify which aspects of climatic variation and change are most influential to spatial population synchrony (Fig. 1). We illustrate the applicability and potential of the protocol through a case study of the gray-sided vole (Myodes rufocanus). This boreal-arctic rodent species is renowned for its important role in ecosystem functioning (Boonstra et al., 2016) and multi-annual population cycles (Hansen et al., 1999;Turchin et al., 2000), with suspected impacts of climate change on these cycles (Ims et al., 2008;Cornulier et al., 2013).

Results
Local gray-sided vole abundances were estimated every spring and fall over 21 years based on capture-recapture sampling in northern Norway (Fig. 2 a, b). Nineteen sampling locations (i.e. live-trapping grids) were spaced along a 170 km transect in boreal mountain birch forest and encompassed three predefined geographic regions (R 1 : coast, R 2 : fiord and R 3 ; inland; Fig. 2a) which were expected to influence the density dependent structure of vole population dynamics.
The 21-year population time series encompass five multi-annual cycles, exhibiting profound overall synchrony across the extent of the study area (Fig. 2b). However, despite visible spatial synchrony, and relative temporal stationarity, there is also some variation in timing and amplitude of the cyclic peaks among the localities. This regards especially the spring series, which have lower and more variable abundance estimates than the fall series (Fig. 2b).
Previous studies have demonstrated that local boreal and Arctic vole populations are adversely affected by winter weather phenomena, such as thaw-freeze cycles (Aars and Ims, 2002;Kausrud et al., 2008) and rain-on-snow events (Fauteux et al., 2021). Hence, we derived local time series of the number of days the temperature crossed zero degrees (Celsius), and the total amount of rainfall (mm) during winter (Fig. 2d). The two weather variables exhibit spatial synchrony, with a tendency for milder (more zero crosses) and wetter (more rainfall) climate towards the coastal area.
2.1 Density-dependence structure Following Stenseth et al. (2003), we fitted second-order log-linear autoregressive models to the population time series, according to the density-dependent (DD) models described in Fig. 1 (models II-IV). As we use a Bayesian framework to conduct the data analysis, we selected Bayesian R 2 (Gelman et al., 2019) as a measure of explained variance (i.e. the fit) of the different linear autoregressive models. In general, the models explained more of the abundance variance in the fall than the spring (Fig. 3). The inclusion of geographic region-specific DD parameters (when comparing model II and III; Fig. 3) did not improve the model fit much, suggesting there are small differences in the DD structure between the three geographic regions. However, a large improvement of the model fit was achieved when including including season-specific DD parameters (model IV; Fig. 3), especially concerning the fall abundances. This implies that season-specific biotic interactions are strongly influential components of the overall population dynamics.

Spatial population synchrony
The spatial correlograms, based on the four population metrics (I-IV) outlined in Fig. 1, clearly show that much of the overall spatial population synchrony (Fig. 4, I) is due to a common DD structure across the study area (Fig. 4, II-IV). Moreover, when accounting for season-specific DD (model IV), the synchrony in the residuals drops substantially in comparison to that from models II and III, which only account for annual DD (Fig. 4, IV). The reduction in spatial synchrony due to seasonal DD is particularly sharp for the fall abundances, for which the synchrony between the most distant populations approaches zero. Accounting for the slight differences in density dependence among the three geographic regions provides almost no contribution to the synchrony pattern (i.e., comparison between II and III in Fig. 4).

Weather synchrony vs. population synchrony
The synchrony of both of the weather variables declined steeply as a function of distance between the sampling stations. However, there was more scattering in the cross-correlations in rainfall when compared to the correlations in the zero crosses (Fig. 5a). The synchrony of number of zero crosses was positively and significantly associated to population synchrony corrected for DD structure (model IV) both in fall and spring, while the synchrony in winter rainfall was only related to the population synchrony in the fall (Fig. 5b,c).

Discussion
We have here proposed and exemplified an analytical protocol, that based on time series data, allows for elucidating deterministic and stochastic sources of spatial population synchrony. Potential deterministic sources include density dependence, climatic seasonality and geographic ecological context, while influential stochastic sources are likely weather variables. Spatial covariance in stochastic weather events amounts to the Moran effect provided that the deterministic components of local population dynamics are linear and identical. Nonetheless, under most circumstances, correlated weather events are expected to exert synchronizing effects when the local density-dependent structure is non-linear and spatially heterogeneous (i.e., the generalized Moran effect; cf. Engen and Saether (2005); Royama (2005); Hansen et al. (2020)). Moran (1953) showed that a key step to make "meteorological phenomena show up more clearly" in statistical analyses of population synchrony is to remove the density-dependent structure from the population time series before making further statistical inferences (e.g., by analyzing the residuals of an autoregressive model). Many studies have used Moran's approach to remove serial autocorrelation in order to fulfill the independence requirement for significance tests of synchrony (Buonaccorsi et al., 2001;Liebhold et al., 2004). However, there appears to be a lack of studies that have followed Moran's suggestion to formally analyze whether the scale of synchrony in the population residuals is dependent on synchrony in the weather (but see Grøtan et al. (2005)); i.e., as achieved by step IV in our analytical protocol. Accordingly, Hansen et al. (2020) conclude that there has been an "analytical deficiency" in empirical Moran-effect studies in terms of making formal inferences about how population synchrony is environmentally forced. We show in the present study that by focusing on residuals which by definition depend on an adequate model structure, we draw more accurate inferences regarding the strength and scale of synchrony.
By applying our analytical protocol to bi-annually sampled time series of gray-sided vole populations we demonstrate winter weather contributions to spatial synchrony. We found that both the amount of rainfall and the frequency of mild-spells in winter contribute to spatial synchrony. These two weather variables have previously been found to affect local population dynamics of boreal and arctic vole species by enhancing winter declines (Aars and Ims, 2002;Fauteux et al., 2021). However, the present study is the first to analytically link large-scale spatial synchrony -a phenomenon that appears to be ubiquitous in boreal and arctic cyclic small rodent populations (cf. Stenseth and Ims (1993); Krebs (2013)) -to any form of stochastic environmental forcing; i.e. Moran effects.
An interesting result arising from our analysis is the time-lagged effect of the winter weather on synchrony of fall abundances. Moran (1953) found similar time-lagged weather effects on an annual time-scale for Canada lynx and speculated about which biological mechanisms could be involved. In voles, environmental conditions in the non-breeding seasons may have lasting effects, for instance, by delaying the onset of reproduction and thereby reducing population growth over the summer (Ergon et al., 2001). The combination of direct and lagged effects of winter weather amounts to an enhanced Moran effect. As increased frequencies of rain-on-snow events and thaw-freeze cycles are very likely outcomes of climate warming in boreal and Arctic ecosystems (AMAP, 2017), we predict that the strength and scale of spatial synchrony of rodent populations will change in these ecosystems.
Climatic seasonality is an externally forced oscillator that acts on the dynamics of most natural systems (Fretwell, 1972). Yet both empirical and theoretical studies of ecological dynamics mostly ignore this fact (White and Hastings, 2020). While seasonality has been shown to be a very important component of spatio-temporal disease dynamics (Earn et al., 1998;Grenfell et al., 2001;Moustakas et al., 2018), we are not aware of empirical studies that have explicitly investigated how such seasonal forcing acts on the strength and scale of synchrony in animal population dynamics. Our analytical protocol provides means for filling this knowledge gap. Specifically, the role of seasonality becomes evident by comparing the correlograms of residuals from models with and without seasonal density dependence (i.e. compare correlograms III and IV in Fig. 4). In the case of sub-arctic gray-sided voles, seasonality is evidently an important determinant of the region-scale spatial synchrony. This regards especially the fall abundances, for which both the overall synchrony becomes reduced and the distance effect is enhanced when seasonal density dependence is accounted for. In this case, it appears that the exact nature of such season-specific effects is contingent on the relative magnitude of the spring and fall noise term of the bi-variate autoregressive model (see Appendix C).
The role of seasonality may be a particularly forceful determinant of spatio-temporal population dynamics in species with multivoltine life cycles, like voles. For instance, the length of winter seasons has been found to exert a strong effect on the local vole population dynamics by acting through density dependent structure (Batzli, 1999;Stenseth, 1999;Stenseth et al., 2003;Bierman et al., 2006) and likely also through season-specific noise terms (Vasseur (2007), Appendix C). Hence, it may not be surprising that seasonality also exerts an effect on regional population dynamics (e.g. large-scale spatial synchrony) as here shown for sub-arctic graysided voles. However, as demographic processes are typically season-specific also in univoltine species (Boyce et al., 1999) -including how they are affected by density-dependent and inde-pendent factors -, we believe that our analytical protocol (Fig. 1) will help advance empirical studies of spatial population synchrony for a wide range of species.

Methods
The methods section follows the structure outlined in Fig. 1, with the five steps required to investigate weather effects on the spatial population synchrony of a gray-sided vole population, after accounting for the geographical-and seasonal-DD structure of the population.

Data and Study Area
We use data from a long-term running monitoring program of the rodent community in the region of Porsanger, northern Norway, between 2000-2020. The data collection consisted of a capture-mark-recapture methodology with two trapping days at 19 individual stations, scattered along a linear transect of approximately 170 km of road. Trapping sessions were conducted twice per year, once in late spring after snow melt, and once at the end of the summer, at the end of the vole reproductive season (see Ehrich et al. (2009) for precise trapping specifications). The Porsanger region contains different landscapes and is subject to a strong climatic contrast (in both temperature and precipitation). The different stations can be sorted into m = 3 regions according to their landscape affinities: coastal region (R 1 ), fjord region (R 2 ) and inland region (R 3 ). Stations 1-5 were included in R 1 (n 1 = 5), stations 6-12 were included in R 2 (n 2 = 7) and stations 13-19 were included in R 3 (n 3 = 7). Fig. 2 summarizes spatial features of the study area and data.

Abundance estimation from mark-capture-recapture data
To reduce a potential bias when estimating synchrony (Santin-Janin et al., 2014), we incorporated the sampling error from capture heterogeneity in our estimates of seasonal abundances (Nicolau et al., 2020). Specifically, we fitted a multinomial regression model to the capture history data to estimate the probability of obtaining a given capture history as a function of individual features registered during the live trapping. These features included the weight and sex of the individuals. We also added a random effect for station in the predictor of the regression model. Individual capture probabilities were subsequently estimated by assuming a temporal effect on the capture process (model M th , Otis et al. (1978)). Finally, the individual probabilities were used to estimate seasonal abundances using an empirical Horvitz-Thompson estimator, which is a function of the estimated individual capture probabilities. Denote the resulting estimated log abundances by {X s,t } and {Y s,t }, for spring and fall, respectively, at spatial locations s = 1, . . . , n s and year t = 1, . . . , n t . For the case study, n s = 19 and n t = 21.

Weather variables
To explore the effect of the weather on the spatial synchrony, we should ideally look into the winter snow conditions (i.e, snow depth and ice formation) as they can be considered the most relevant climatic variables affecting rodent population cycles (Hansson and Henttonen, 1985;Hansen et al., 2013;Stien et al., 2012;Fauteux et al., 2021), being an interaction of different weather variables, including temperature and precipitation. As this information was not directly available, we resorted to proxy variables of snow conditions, using the temperature and precipitation estimates from the Norwegian Meteorological Institute between 2000-2020. These measurements correspond to model estimates (not measured at station level) and are prone to large uncertainties, particularly the precipitation (Lussana et al., 2019). For our proxies of snow conditions, we derived two variables: winter zero crosses, as the total number of times the mean daily temperature crossed 0°C during winter (21 Dec -20 Mar); and winter rainfall, as the precipitation sum in days where the mean temperature surpassed 0°C, during winter.

Statistical Framework
We describe the statistical framework to decompose density-dependence contributions into the spatial synchrony of populations, and isolate the weather effects on population synchrony, described in steps 3-5 in Fig. 1.

Density-dependence structure (Step 3)
The general protocol (Fig. 1, step 3) specifies three different models for the DD structure of the estimated time series. Here, we assume that the general function f (.) is linear, describing the log-DD structure in terms of direct and delayed effects up to lag p. Specifically, the three models either include or exclude regional-and seasonal-dependent effects as specified below.
In general, the spatial locations s are assumed to be within a closed geographical region R, which can be partitioned into m mutually exclusive subregions, R = R 1 ∪ . . . ∪ R m . For the gray-sided vole case study, this corresponds to m = 3 regions. The most general model includes both regional-specific and seasonal-specific terms (Fig. 1, model (IV)), and the assumed loglinear dependency structure up to order p can be expressed by X s,t = β r1 Y s,t−1 + β r2 X s,t−1 + . . . + β r,2p−1 Y s,t−p + β r,2p X s,t−p + s,t (1) Y s,t = γ r1 X s,t + γ r2 Y s,t−1 + . . . + γ r,2p−1 X s,t−1 + γ r,2p Y s,t−p + ω s,t where t = p + 1, . . . , n t and s ∈ R r . The terms st and ω st denote individual random environmental noise at each spatial location s for each time point t, while the sets of regionaland seasonal-specific coefficients can be summarized as Θ Rx = {β r1 , . . . , β r,2p } and Θ Ry = {γ r1 , . . . , γ r,2p }.
Simplifications of the given model will yield more simplistic measures of the DD structure. According to the general protocol in Fig. 1, Model (I) corresponds to assuming no DD structure, in which all of the given coefficients are equal to 0. This corresponds to simply using the estimated raw log-abundance series, {X s,t } and {Y s,t }, in further analysis.
Following (Stenseth et al., 2003), we included delayed effects up to order p = 2 for the case study. Model (II) refers to a second-order annual autoregressive processes including coefficients Θ = {β 2 , β 4 , γ 2 , γ 4 } which are neither regional-specific (m = 1; disregarding spatial heterogeneity), nor seasonal-specific (β 1 = β 3 = γ 1 = γ 3 = 0; assuming yearly dynamics). Such AR(2) models are often used in literature (e.g. Turkia et al. (2020); Dallas et al. (2020)). Model (III) is characterized by incorporating regional-specific effects {β r2 , β r4 , γ r2 , γ r4 } m r=1 . This corresponds to AR(2) models which allow for spatial differences in the DD structure which can account for some of the observed synchrony (Hugueny, 2006). Finally, by including the seasonal-specific effects {β r1 , β r3 , γ r1 , γ r3 } m r=1 , we get the bivariate model (IV) which is very similar to a second-order vector autoregressive model (VAR). The difference to a VAR-model, however, is that the time series {X s,t } and {Y s,t } are observed at two different time points in year t, and the fall log-abundances are modeled in terms of the spring observations within the same year. Seasonal-specific DD has been recognized as fundamental to model small rodent population dynamics (Hansen et al., 1999), but to our knowledge seasonal DD contributions to spatial synchrony have not been assessed.

Measuring the scale and shape of spatial population synchrony (Step 4)
To assess the scale and shape of the spatial synchrony, we can consider the spatial correlations of the environmental noise terms in Models (I-IV) as a function of geographical distance. The following analysis is repeated using the four different models for DD structure, specified in section 4.2.1. A major goal is then to understand how the inclusion of regional-and seasonalspecific terms influences the synchrony estimates, i.e., which part of the synchrony is explained by the different DD components.
Define the residual vectors s = ( s,1 , . . . , s,nt ) and ω s = (ω s,1 , . . . , ω s,nt ) for all spatial locations s = 1, . . . , n s and t = p+1, . . . , n t . The contributions to the spatial synchrony are then characterized by the pairwise correlations between vectors within each of the sets { s } ns s=1 and {ω s } ns s=1 . If the associations between these residual series are expected to be linear, the degree of synchrony is typically measured using Pearson's correlation coefficient (Bjørnstad et al., 1999;Liebhold et al., 2004). To model the correlations in terms of geographical distance, let δ i,j denote the Euclidean distance between two stations i and j. In accordance with calculating the spatial correlogram (Bjørnstad et al., 1999;Bjørnstad and Falck, 2001;Liebhold et al., 2004), we discretize the n s (n s − 1)/2 unique distances between stations into distance classes d k , k = 1, . . . , K, where K is the total number of classes. Specifically, a distance class d k is defined by L k < δ i,j < U k , where L k and U k represent the lower and upper bound of the distances within that class, respectively. The corresponding averages of the pairwise correlations {ρ i,j } for distance class d k are then given by where n k is the total number of distances/correlations within distance class d k . The given formulation is analogous to the calculation of Koenig's modified correlogram (Koenig, 1999;Bjørnstad et al., 1999), as the correlations are not centered (zero synchrony is taken as the reference line of the correlogram). For the given case study, we assumed that the distance-class width is U k − L k = 1 for all classes, which corresponds to rounding off the geographical distances to the nearest integer. We used this method to calculate the averaged correlations in (3) as a pre-processing step to reduce random noise in the estimated correlations.
As an alternative to using the non-parametric covariance function (Bjørnstad and Falck, 2001) or other non-parametric estimates of the correlation function (Liebhold et al., 2004), we chose to model the correlations in terms of the distances using the regression model Here, f denotes a smooth underlying function while {ν k } represents zero-mean, independent Gaussian error terms with constant variance. This model is fitted using a Bayesian framework where the function f is assigned a second-order intrinsic Gaussian Markov random field prior (Rue and Held (2005), page 110). The model is scaled according to Sørbye and Rue (2014) and the precision parameter of the model is assigned a penalized complexity prior with parameters U = 0.5 and α = 0.01 (Simpson et al., 2017). Using the methodology of integrated nested Laplace approximation (Rue et al., 2009), both the posterior mean and credible intervals for f are calculated efficiently, without the need of resampling techniques, like Monte Carlo simulation or bootstrapping.

Effects of Weather (Step 5)
Finally, we can use the measures of synchrony accounting for the effects of geographic-and seasonal-dependent DD to investigate potential weather (or other relevant environmental variables) drivers. For this, we can model the set of correlations {ρ k (d k )} K k=1 from model (IV) as a function of the corresponding spatial correlations of different weather covariates, defined by {ρ The availability of such covariates are typically case-specific but should be measured or estimated to represent the same spatial locations and time points used for the logabundance estimates. For the given case study, the relationship between the weather variables (zero crosses and winter rainfall) appeared to be linear, and was thus modeled using simple linear regression models.

Synchrony (residual decomposition and weather)
Distance Distance Distance Population Synchrony  The five main steps of a general methodological protocol to single out the impacts of climatic variation (weather) on spatial population synchrony, by accounting for seasonal and geographical contexts in ecological population processes (density dependence).

Effects of weather
Step 1: Seasonal (spring and fall) sampling of both local populations and a focal weather variable at different locations (crosses). The geographic sampling frame encompasses two regions (R 1 and R 2 ) representing different geographic ecological contexts (e.g. habitats or ecological communities).
Step 2: Season-and region-specific time series of local population abundance estimates resulting from the sampling process, together with time series of the focal weather variable. The estimation of abundance ideally involves separating the observation process and the population process, accounting for detectability.
Step 3: Four alternative models to further analyze spatial population synchrony. (I) corresponds to seasonal abundance estimates (X t and Y t ). (II-IV) correspond to the sets of X t − f (.) and Y t − f (.) from the respective general models for density-dependence, modeling state of the population at time t as a function of previous p states. Model (II) includes only one set Θ of density-dependence parameters with annual time lags (i.e. ignoring seasonal and regional components). Model (III) includes region-specific parameters Θ R , again with annual time lags (i.e. ignoring seasonal components). Model (IV) is a bi-variate model (Stenseth et al., 2003) that includes both geographic-and season-specific parameters Θ R X and Θ R Y .
Step 4: Season-specific synchrony patterns (i.e. scale and shape) of the population (derived from Step 3) and weather metrics (derived from Step 2) as function of distance. The dots are the pairwise cross-correlations of the population metrics and the weather variables, while the lines are estimated correlograms with associated uncertainty intervals (e.g., Bjørnstad et al. (1999)).
Step 5: Estimated effects of weather synchrony on population synchrony. Seasonspecific (Fall and Spring) population synchrony with corrected for seasonal-density dependence and geographic context effects (i.e. residuals from model (IV)) are regressed against the spatial synchrony in the focal weather variable. Illustrations created with Biorender.com.  Associated 95% credible intervals correspond to the bars, with the coefficients associated with each region represented with the respective color (blue for R 1 , green for R 2 and yellow for R 3 ).

C Simulations
To assess when the observed changes in seasonal correlations can occur, we made a short simulation study. We simulated 100.000 pairs of time series, each pair corresponding to 2 populations with a true process described by model IV (see Methods), with regression parameters approximating those of the real data set (section B; parameter values were {β 1 = 0.5, β 2 = 0.1, β 3 = −0.4, β 4 = 0.1, γ 1 = 0.8, γ 2 = 0.1, γ 3 = −0.1, γ 4 = −0.1}. We then compared (a) the true correlation in the noise terms for each pair (which corresponds to a parameter in the simulation; we used r=0.5 for the spring correlation and r=0.3 for the fall correlation), against: (b) the correlations when using year-to-year raw abundances for either spring and fall (model I); and (c) correlations in the residuals of yearly AR(2) for each season separately (model II). Figure 10 displays the noise correlation estimates for each of the approaches, using different variances of spring and fall noise. We can infer that the decrease we have observed for fall densities (i.e., summer season) happens when the winter noise is as large or larger than the summer noise, which is expected be the case in study system (the estimated variances were 0.86 for spring and 0.86 for fall). If the summer noise is larger than the winter noise, the patterns are reversed compared to the case when winter noise is larger than summer noise . Note that correlations based on annual models can be higher than the true ones. This confirms that the observed changes come from considering the seasonal model (IV) instead of the annual models (I-III), rather than artifacts related to the data. This shows that fluctuating (i.e., seasonal) environments can enhance synchrony, but the exact pattern obtained using annual time series will depend on the noise structure and which season is monitored (Vasseur, 2007).  Figure 10: Simulating noise term correlations for either Spring and Fall. X-axis corresponds to the noise terms of either Spring (left) and Fall processes (right). This corresponds to the and ω terms in equations 1-2. The mean true correlations in the noise terms (a) are the black dots; the mean correlations in the raw abundances are the red dots; are the mean correlations in the AR(2) model residuals (c) are the blue dots.