The dependence of synchrony on timescale and geography in freshwater plankton

Spatial synchrony is defined by related fluctuations through time in population abundances measured at different locations. The degree of relatedness typically declines with increasing distance between sampling locations. Standard approaches for assessing synchrony assume isotropy in space and uniformity across timescales of analysis, but it is now known that spatial variability and timescale structure in population dynamics are common features. We tested for spatial and timescale structure in the patterns of synchrony of freshwater plankton in Kentucky Lake, U.S.A. We also evaluated whether different mechanisms may drive synchrony and its spatial structure on different timescales. Using wavelet techniques and matrix regression, we analyzed phytoplankton biomass and abundances of seven zooplankton taxa at 16 locations sampled from 1990 to 2015. We found that zooplankton abundances and phytoplankton biomass exhibited synchrony at multiple timescales. Timescale structure in the potential mechanisms of synchrony was revealed primarily through networks of relationships among zooplankton taxa, which differed by timescale. We found substantial interspecific variability in geographic structures of synchrony and their causes: all mechanisms we considered strongly explained geographic structure in synchrony for at least one species, while Euclidean distance between sampling locations was generally less well supported than more mechanistic explanations. Geographic structure in synchrony and its underlying mechanisms also depended on timescale for a minority of the taxa tested. Overall, our results show substantial and complex but interpretable variation in structures of synchrony across three variables: space, timescale, and taxon. It seems likely these aspects of synchrony are important general features of freshwater systems.

correlated environmental variability (i.e., Moran effects; Moran 1953), dispersal between populations (Bjørnstad et al. 1999;Liebhold et al. 2004), or mobile or synchronized predators that synchronize prey populations (Ims and Andreassen 2000;Sundell et al. 2004). Understanding and differentiating among these causal mechanisms of spatial synchrony has been a longstanding goal in ecology (Liebhold et al. 2004). Synchrony and our understanding of it has ramifications for conservation and other applications of population biology. For instance, populations that are synchronously low in abundance are at greater risk of extinction (Heino et al. 1997;Earn et al. 2000), and pest populations that are synchronously abundant constitute an outbreak or epidemic (Williams and Liebhold 2000;Okland et al. 2005).
Standard methods for assessing spatial synchrony typically include calculating the mean of correlation coefficients between all pairs of locations (Buonaccorsi et al. 2001) and plotting correlation coefficients against geographic distances between locations. The typical expectation is that correlation decreases with increasing distance because closer populations are more likely to experience similar biotic or abiotic conditions and to exchange dispersers (Bjørnstad et al. 1999;Koenig 2002). Researchers have sometimes attempted to ascertain drivers of synchrony by visually comparing the patterns of distance decay in correlation for population data and putative drivers. This approach has provided useful information, but it has limitations. Foremost, multiple drivers can produce similar patterns of distance decay (Kendall et al. 2000;Abbott 2007). The use of distance-decay relationships to assess synchrony may obscure the potential for uncovering more complex spatial patterns that are likely to occur in natural systems (Gurevitch et al. 2016;Walter et al. 2017). Furthermore, approaches based on standard correlation focus on in-phase relationships (i.e., relationships with no time lag) and conflate potentially distinct patterns of dynamics occurring on different timescales (Defriez et al. 2016;Sheppard et al. 2016;Defriez and Reuman 2017a;Defriez and Reuman 2017b). Only recently have studies begun to investigate more complex spatial and timescale structure in synchrony.
Several recent studies have identified mechanisms that can create detailed geographic patterns of synchrony (e.g., Powney et al. 2012;Haynes et al. 2013;Gouveia et al. 2016), here referred to as the "geography" of synchrony (Walter et al. 2017). For example, synchrony in primary production varied spatially in both terrestrial and oceanic systems, which was ultimately linked to Moran drivers (e.g., temperature and rainfall) with synchrony that had similar geographic patterns (Defriez and Reuman 2017a;Defriez and Reuman 2017b). Spatial structure in the synchrony of Moran drivers has been the primary cause investigated thus far for the geography of synchrony, although several other processes also could produce geographic structure in synchrony, including spatial variation in density-dependent processes or dispersal (Walter et al. 2017). There is an important conceptual distinction to be made, here, between mechanisms of synchrony itself (e.g., Moran effects, dispersal), which are well studied, and mechanisms of geography of synchrony, which are often, but not necessarily, related to mechanisms of synchrony themselves. For instance, a Moran driver can be the mechanism of synchrony in a population variable, whereas it is the spatial variation in the synchrony of that driver that is the mechanism of the geography of population synchrony. Other mechanisms of the geography of synchrony operate by modifying the effectiveness, in a spatially structured way, of a separate synchronizing driver. For instance, spatial variation in density dependence can modify the effectiveness of Moran drivers, producing a geography of population synchrony even if there is no geographic pattern in the synchrony of the Moran driver. These ideas are explored in detail by Walter et al. (2017), but, in general, investigations of the geography of synchrony are still in their infancy. There is much to be learned, for instance, on how geography of synchrony may manifest for different taxa, ecosystems, spatial scales, and other axes of ecological variation.
In addition to detailed spatial structure in synchrony, populations can also exhibit synchrony with important timescale structure (Fig. 1a-d), i.e., synchrony can manifest itself differently on some timescales than others (Vasseur and Gaedke 2007;Defriez et al. 2016;Sheppard et al. 2016;Walter et al. 2017). Here, the timescale of a periodic oscillation refers to its period, the reciprocal of its frequency. Ecological time series usually are a combination of fluctuations occurring on a variety of timescales (e.g., annual and decadal), and fluctuations on some timescales can be synchronized across space differently than fluctuations on other timescales .
Timescale structure of synchrony has been used to facilitate inference of the drivers of synchrony. For example, nearly 80% of aphid synchrony at timescales greater than 4 yr was explained by synchrony of winter temperatures occurring at the same timescales ), but much less shorter timescale synchrony was explainable. We believe that approaches leveraging timescale-specific information in synchrony have been underapplied in ecology, limiting our understanding of this aspect of synchrony. The possibility that a particular mechanism of spatial synchrony may be operating in a given system primarily at specific timescales is also underexamined, and that oversight probably has limited our ability to determine mechanisms of synchrony. Furthermore, timescale structure in synchrony and its drivers may vary spatially, such that different locations are synchronized at different timescales (Fig. 1e,f ), potentially because different drivers of the geography of synchrony occur at specific timescales.
Synchrony in freshwater plankton is well studied via distance-decay and correlation approaches (e.g., Rusak et al. 1999;Rusak et al. 2008;Seebens et al. 2013;Lodi et al. 2014) which conflate components of synchrony on different timescales (Fig. 1). Only a small minority of studies have tested whether and why longer timescale, interannual fluctuations, specifically, may be spatially synchronized (Keitt and Fischer 2006;Vasseur and Gaedke 2007;Winder and Cloern 2010;Carey et al. 2016). Most prior studies of freshwater plankton have instead focused on seasonal fluctuations in abundance that relate to life history variation among species (Sommer et al. 1986). Thus, causal mechanisms that may induce spatial synchrony over multiannual time periods are less explored. A rich literature on plankton population and community dynamics suggests many environmental factors that could drive Moran effects on interannual timescales (Dodson et al. 2005;Yan et al. 2008;Shurin et al. 2010), some of which have recently been shown to affect synchrony on longer timescales in marine plankton (e.g., Defriez et al. 2016). Shorter timescale fluctuations (Fig. 1c,d) could include annual phenological patterns in zooplankton, whereas longer timescale fluctuations may occur due to distinct Moran effects (e.g., synchrony in pH) operating at longer timescales. This study provides new information by examining the longer, interannual timescales rather than the better known intraannual and phenological patterns.
Prior studies of synchrony in freshwater plankton also have generally examined synchrony between different water bodies (Rusak et al. 1999;Rusak et al. 2008;Vogt et al. 2011;Pandit et al. 2016), rather than within a single water body. Using separate water bodies limits comparative investigations of drivers of synchrony (i.e., dispersal, predation, and Moran effects), as dispersal among lakes is unlikely to be sufficient to synchronize their plankton dynamics. Investigations of synchrony of plankton within a single water body permit more direct comparisons of mechanisms driving synchrony, making it possible to infer which mechanism of synchrony dominates because of the geography of synchrony that is observed. For instance, Moran effects vs. current-induced dispersal may be expected to yield distinct geographies of synchrony. Furthermore, geographies of plankton synchrony may vary by timescale, possibly revealing how mechanisms of synchrony also vary by timescale. For example, abundance of a zooplankter at some locations within a lake may have correlated fluctuations at longer timescales driven by Moran effects, whereas populations at a different set of locations may be synchronized at shorter timescales by dispersal (Fig. 1e,f ). Few analyses of synchrony using data from single bodies of water have been performed, and those studies that have been done (Lansac-Tôha et al. 2008;Seebens et al. 2013;Lodi et al. 2014) included limited examination of synchrony mechanisms with little to no focus on decomposing synchrony by timescale and/or geography.
We assess the geographic and timescale structure of synchrony in freshwater plankton using long-term data collected . Combining y1 with y3 and y1 with y4 gives two time series (c), here imagined to be obtained from measurements at different locations, which are synchronized on long timescales, but antisynchronized on short timescales. The reverse is also possible (d). The timescale-specific structure of synchrony cannot be detected with correlation coefficients, which are 0 for both c and d, because contributions from different timescales are conflated by correlation methods and cancel (see Defriez et al. 2016, for further discussion). Dependence of geography of synchrony on timescale can also be modeled with our simple example by plotting time series arbitrarily in Cartesian space, creating a simple demonstration of geographies of shortand long-timescale synchrony which differ (e, f ). This figure is an augmented version of a figure from Defriez and Reuman (2017a). from multiple locations in Kentucky Lake, a large reservoir in western Kentucky, U.S.A. (Fig. 2). Reservoirs provide a good environment in which to investigate geographic patterns of synchrony for plankton because the diversity of habitats that reservoirs contain, including both lotic and lentic areas (Thornton et al. 1990), creates the potential for geographic structure in abiotic conditions (i.e., Moran effects). Moreover, directional movement of water in more lotic areas creates potential for highly directed dispersal, which would tend to homogenize connected population phenomena. Water flow also can impede dispersal, making locations on opposite sides of a strong current less mutually accessible to dispersers. Each of these examples has the potential to generate geographic structure in synchrony. We previously documented strong overall spatial synchrony of multiple biotic and abiotic variables that were measured in Kentucky Lake (range of cross correlation values = 0.49-0.82). We also showed that synchrony (as measured by correlation) did not decline systematically with geographic distance between locations (Anderson et al. 2018; Supporting Information File S1), in essence because sites on opposite sides of the lake can be close together or far apart, but are always limnetically separated by the channel, and sites on the same side can be close together or far apart but are always relatively connected by water flow patterns. Here, we leveraged relatively new approaches to assess timescale and geographic structure in plankton synchrony. We examined two specific questions within each of these overarching areas of inquiry: (Q1a) Does strength of spatial synchrony depend on timescales of analysis? (Q1b) What is the relative importance of different potential mechanisms of synchrony on different timescales? (Q2a) Do geographies of synchrony depend on timescale? And (Q2b) does the relative importance of mechanisms that are likely driving observed geographies of synchrony depend on timescale? Recently developed analytical approaches applied here help unravel complex structure and mechanisms Walter et al. 2017), resulting in the first or one of the first detailed and systematic examinations of geographic and timescale-specific aspects of synchrony and its potential mechanisms in freshwater plankton.

Study site
Kentucky Lake forms the terminal and largest mainstem reservoir on the Tennessee River in western Kentucky, U.S.A. (length ≈ 300 km, width ≈ 2 km, surface area ≈ 650 km 2 , mean depth ≈ 6 m; Fig. 2). The reservoir is considered mesotrophic and is well mixed by wind and water currents. It has a very short retention time, averaging less than 30 d (Bukaveckas et al. 2002;Yurista et al. 2004), making the reservoir functionally more riverine than lacustrine. Similar to other mainstem impoundments, the original channel is the deepest part of the reservoir (max depth of ≈ 21 m in summer), and much of the inundated surface area is the old flood plain. Water depth varies by approximately 2 m between winter and summer. Discharge from the dam, about 24 km downstream from the study sites, averages approximately 40,000 m 3 s −1 , but varies substantially with power generation and flood control situations.

Data collection
The ongoing Kentucky Lake Monitoring Program (KLMP) regularly collects samples of physiochemical and plankton variables in Kentucky Lake at 16+ sites every 16 d during the spring-fall months and every 32 d during the winter months ). See Bukaveckas et al. (2002)  Program (KLMP) on Kentucky Lake, KY. Colored symbols represent sampling locations for plankton and abiotic variables (pentagons, squares, circles, and triangles) and crappie (stars). Crappie were sampled over several locations within two embayments, but data were combined to produce a single time series for crappie. of sampling protocols. For this analysis, we focused on a 26-yr period  for 16 primary sampling sites (Fig. 2) that occur within an embayment arm on the western shore (N = 4 sites; mean depth = 3.9 m), in embayment mouths on the western shore that primarily drain agricultural land (N = 6 sites; mean depth = 6.5 m), in embayment mouths on the eastern shore that drain the Land Between the Lakes National Recreation Area (primarily forested land, N = 3 sites; mean depth = 7.5 m), and sites along the original river channel (N = 3 sites; mean depth = 16.9 m).
Zooplankton communities were sampled using a 15-liter Schindler-Patalas trap (fitted with a 243 μm sieve) that was lowered to 5 m below the surface or half of the maximum water depth (whichever was shallower) and then retrieved. Samples were collected in triplicate at each site on each visit, and total numbers of individuals of each species were calculated across samples to estimate total abundance per species and sampling occasion (number of individuals per 45 liters). Only two workers conducted the zooplankton identification over the entire survey period (1990-2012 and 2012-2015 were their respective periods of employment). Analyses were focused on the most abundant zooplankton taxa: calanoid and cyclopoid copepods, Bosmina longirostris, Daphnia lumholtzi, Daphnia retrocurva, Diaphanosoma birgei, and Leptodora kindtii. See Williamson and White (2007) for a list of the dominant copepod species observed in Kentucky Lake, though here counts of copepods were by Order. Analyses here represent data from approximately 7800 zooplankton samples. We assume that L. kindtii and cyclopoida are primarily predaceous, calanoida are primarily omnivorous, and the cladocerans other than L. kindtii are primarily herbivorous.
Fourteen environmental parameters were selected for this analysis that are known to be important to plankton dynamics (Dodson et al. 2005;Rusak et al. 2008;Shurin et al. 2010) and that potentially could produce Moran effects. At each sampling site, water temperature, dissolved oxygen, conductivity, and pH were recorded at 1 m intervals throughout the entire water column using a YSI multiparameter sonde (Yellow Springs Instruments). Water samples were collected for laboratory analyses from 1 m below the surface and 1 m above the lake bottom using a 2-liter Kemmerer sampler. One-liter subsamples of water were filtered through 1.2 μm, 2.5 cm Whatman GF/C glass fiber filters and stored on ice prior to chemical analysis. Chlorophyll a was extracted from the water collected 1 m below the surface using pigment extraction and spectrophotometric methods and provided an estimate of phytoplankton biomass (APHA 1989). The following parameters were determined from both near-surface and bottom water samples: dissolved and total phosphorus (DP, TP), soluble reactive phosphorous (SRP), ammonium (NH 4 ), silicon dioxide (SiO 2 ), and dissolved and total nitrogen (DN, TN). See Bukaveckas et al. (2002) and Yurista et al. (2004) for detailed descriptions of the lab techniques used to determine concentrations of each chemical. Secchi depth was recorded at each sampling site using a 20 cm Secchi disk. Discharge rates (m 3 s −1 ) from Kentucky Lake Dam, approximately 25 km down current from the study area, were taken from Tennessee Valley Authority River Reports. Annual (1990Annual ( -2015 catch-per-unit-effort abundance estimates of age-1 crappie (counts included both black [Pomoxis nigromaculatus] and white crappie [P. annularis]) were obtained from the Kentucky Department of Fish and Wildlife Resources. Crappie data were not collected in the same way as the other biotic and abiotic data. Instead, these data were collected in only two embayments (Fig. 2) in relatively standardized locations within each embayment. For data quality reasons and because two locations were not enough for a spatiotemporal comparison of crappie and plankton data, crappie data were combined into a single time series representing a coarse annual estimate of lake-wide average abundance of crappie. Young age classes of crappie feed partially to exclusively on zooplankton (Ellison 1984;Martin 2012).

Data preparation
Measurements listed above were used to construct a database, which had no missing values, for the abundances of seven zooplankton taxa, phytoplankton biomass, zooplanktivore abundance (crappie), and five physical variables (discharge, specific conductance, pH, Secchi depth, and water temperature). We assumed an abundance of 0 for zooplankton taxa when sampling occurred but no individuals of that taxon were captured. A nearly complete database of eight chemical variables (NH 4 , DN and TN, DP and TP, SRP, SiO 2 , and dissolved O 2 ) was also constructed. TN and TP had several missing values in 1993, which we filled with mean values over sampling trips carried out at the same time of year for a given site over years other than 1993 for which data were available. We also discarded values of TN and TP and DN and DP from all sites on two sampling dates, SRP from two sampling dates at two individual sites, and NH 4 from one site in 1 yr, due to anomalously high values; these were not filled with means of data from other years because sufficient sampling occurred during the affected years so that annual means could simply be constructed from the other, unaffected sampling occasions.
We generated annual mean values (using the arithmetic mean) for each variable using data from April to November, which approximated the growing season, when maximum depths occur in Kentucky Lake and when KLMP performed sampling at 16-d intervals. We then had complete annual time series for all variables and locations ( Fig. 3; Supporting Information File S2). Seasonal plankton dynamics are much more commonly studied (Sommer et al. 1986(Sommer et al. , 2012 than dynamics on interannual timescales, although strong and ecologically meaningful patterns of synchrony are present on these longer timescales. Because they are relatively little studied, we focused on the longer timescales. We log(x + 1)transformed annual mean values for plankton variables to improve normality. For ancillary analyses, for each plankton variable, we also computed the means of the log(x + 1)transformed abundances from the separate sampling occasions within a location and year to compute an alternative annual value for that variable, location, and year. This alternative to the arithmetic mean used in the main analyses was similar to the geometric mean but could cope with zeros, and was used to discover whether outliers that may have influenced arithmetic means more than geometric means could have changed results. Each site's annual time series for each biotic or abiotic variable was linearly detrended to remove any longitudinal trends that might obscure synchrony patterns (Buonaccorsi et al. 2001). Detrending involved regressing against year for each location, extracting the residuals, and then dividing by the standard deviation of those residuals.
We tested whether partitioning data collected at different depths into surface samples (within 2 m of the surface) and Log(x + 1) abundance of plankton separated by sampling site on Kentucky Lake. Colors represent the different limnetic groupings (red squares = embayment, orange plus signs = west shore, blue circles = channel, green triangles = east shore). Each point is the annual mean value per site for data collected from April to November. All units are the total number of individuals per 45 liters, except for phytoplankton which is μg/L. deeper water samples (> 2 m depth) affected results. We found similar results for most analyses, using either shallow or deep variables, so we only report results for analyses using data averaged over all depths.

Analysis overview
Distinct statistical methods were used to answer each of the main questions from the Introduction (Q1a,b; Q2a,b; Table 1). For Q1a, we examined timescale structure in the synchrony of each abiotic and biotic variable using wavelet phasor mean field magnitude plots. For Q1b, we used timescale-specific spatial wavelet coherence methods to test for timescale differences in potential drivers of synchrony (environmental Moran drivers or species interactions). For Q2a, we used cross-wavelet analysis to measure geographic variation in synchrony in a timescale-specific way, and then tested whether the spatial structure of synchrony varied by timescale using Mantel tests. And for Q2b, we used matrix regression methods to determine the likely relative importance of different mechanisms of the geographic structure of synchrony, making comparisons between results obtained for different timescales. Some of these methods are recently developed, but mathematical details and prior ecological applications of all methods are published elsewhere (Sheppard et al. 2013;Sheppard et al. 2016;Sheppard et al. 2017;Walter et al. 2017), so we have not repeated the mathematical content. We instead describe each method below from an operational viewpoint, indicating the general idea and goal of each analysis, and the input data that each method requires and how to interpret output. See Table 1 and Supporting Information Files S3-S5 for further details and heuristic examples. All methods were implemented in R (R Core Team 2018).

Q1a: Timescale structure of synchrony
We characterized timescale-specific synchrony between locations by generating wavelet phasor mean fields (Supporting Information File S3) for each of our biotic and abiotic variables (except for crappie and discharge, for which suitable spatial data were not available). For each variable, we computed the continuous Morlet wavelet transforms of the time series for individual sites. The wavelet transform detects and decomposes the variational content of a time series by time and timescale (Addison 2002). The technical specifications for the mother wavelet and for the spacing of wavelet scales were exactly the same as those of Sheppard et al. (2016). For each variable, we combined the wavelet transforms for all sites to generate a wavelet phasor mean field magnitude plot (Supporting Information File S3), which visually depicts the strength of spatial synchrony among sampling locations as a function of time and timescale. While visually similar to plots of wavelet power (e.g., Winder and Cloern 2010;Carey et al. 2016), wavelet phasor mean field magnitude plots display distinct information, namely the degree of phase consistency (a measure of the strength of spatial synchrony) among locations as a function of time and timescale (Supporting Information File S3). Wavelet transforms are less reliable for longer timescales and times closer to the beginning or end of available data (Addison 2002), so the wavelet phasor mean Table 1. Summary framework of our research questions, data inputs, analytical methods, locations of associated results in the paper, and citations that describe methodological details. For all questions, data inputs were log(x + 1)-transformed annual mean abundances of zooplankton or chlorophyll a biomass from our 16 locations, linearly detrended (see Methods section). These data were then transformed into the "Transformed data" inputs each method required. See Methods section for details.

Question
Transformed data Method Result Citation  field, being based on the wavelet transform, has the same property. Our plots are therefore "scalloped" (i.e., trimmed at the edges, more so for longer timescales) to omit unreliable values (see Fig. 4). We tested for nominal significance of synchrony for each time and timescale (without correcting for testing multiple times and timescales) by comparison with a distribution of phasor mean field values consistent with a null hypothesis of no synchrony between sites (Sheppard et al. 2013).
Q1b: Timescale-dependent mechanisms of synchrony We used the spatial coherence technique Sheppard et al. 2017), which measures the degree to which two spatiotemporal variables measured at the same locations and times (e.g., phytoplankton and temperature time series at all Kentucky Lake sampling locations over the period 1990-2015) are related to each other, in a timescalespecific way. In a general sense, spatial coherence is a statistical method that tests whether patterns of spatial synchrony are similar among two spatiotemporal variables. Whereas measures of spatial synchrony detect in-phase relationships between time series of the same variable measured in different locations, spatial coherence detected relationships, which are consistent over space and time, between two different spatiotemporal variables. Technically, it shows whether the variables exhibit consistent phase differences and correlated magnitudes Fig. 4. Wavelet phasor mean field magnitudes for zooplankton abundance and phytoplankton biomass data. Higher values indicate stronger spatial synchrony between sites at specific times (x-axis) and timescales (y-axis). Timescale units are years. Contour lines indicate times and timescales that exhibited nominally significant synchrony at p < 0.05 (solid lines) and p < 0.01 (dashed lines), without correction for testing at multiple times and timescales. Right panels for each taxon show time-averaged synchrony by timescale and summarize the timescale structure in synchrony. of oscillation through time and across space as a function of timescale (Supporting Information File S4). Given two spatiotemporal variables, spatial coherence can be interpreted as being similar to a timescale-specific correlation coefficient between the variables, but with the further advantage of detecting a specific phase relationship between the variables, if one occurs, without having to specify in advance what the expected phase difference is (Supporting Information File S4). For instance, in-phase and antiphase relationships are both detected, as are any time-lagged relationships, with equal effectiveness regardless of the lag. Traditional correlation, on the other hand, only detects in-phase relationships between variables. Spatial coherence is superior to correlation methods for our purposes because time-lagged relationships can produce correlations of 0, falsely suggesting no relationship between variables, and because spatial coherence decomposes potential relationships between variables by timescale.
We tested for significant spatial coherence between each pair of spatiotemporal variables for which we had data, including all zooplankton taxa, phytoplankton biomass, and all environmental variables, separately for "short" and "long" timescales (see below), thereby generating two timescalespecific networks of significant spatial coherences. Network nodes were variables, and connections between nodes were significant spatial coherences between variables. Thus, networks revealed patterns by which variables were related, for short and long timescales. See Supporting Information File S4 for information on significance testing of spatial coherences, which followed the methods of Sheppard et al. (2017).
In addition to spatial coherences between pairs of spatiotemporal variables, we also included in our networks coherence results with crappie and lake discharge time series, although spatially resolved data for these variables were not available. The same crappie or discharge time series was associated with all plankton sampling locations in the lake for these comparisons, constituting an assumption that crappie and discharge variables, as spatially averaged and whole-lake quantities, respectively, were related in the same way to other variables at all locations.
We here and henceforth use the term "short timescales" to refer to the 2-4-yr timescale band; the term "long timescales" refers to > 4-yr timescales. One cycle every 4 yr was chosen as the dividing line between short and long timescales because it was exactly half the Nyquist frequency (1 cycle every 2 yr for annual data), and because it is a boundary between persistent and antipersistent dynamics (i.e., successive values are more similar or dissimilar, respectively) in Fourier components, as measured with lag-1 autocorrelation . Annual or intra-annual timescales were not considered in our analyses; they corresponded to frequencies that were above the Nyquist frequency for our annualized data. Seasonal fluctuations in plankton are already well studied (Sommer et al. 1986; Winder and Cloern 2010), as explained above.
As opposed to correlation-based approaches that assume inphase relationships, spatial coherence allows for any kind of phase relationship, and output of the technique includes the average phase that occurred. Phase information was incorporated into our networks and used for comparing the short-and long-timescale networks. When two variables were significantly spatially coherent across a timescale band (short or long timescales), we calculated the mean phase difference between the variables over the band, in units of π, to determine if oscillations had approximately in-phase (φ = −0.25 to 0.25), quarter-cycle (φ = 0.25 to 0.75; OR −0.25 down to −0.75), or antiphase (φ = 0.75 to 1.00; OR −0.75 down to −1.00) relationships, and we associated this information to the corresponding network link.
Our networks can reveal differences between mechanisms of synchrony in Kentucky Lake operating at short vs. long timescales because the networks reflect the mechanisms of synchrony that operated; if mechanisms were the same across timescales, then our two networks should be similar, but substantial differences between networks suggest differences in mechanisms. First, direct transmissions of synchrony between the variables we have measured (e.g., Moran effects on zooplankton) has the potential to be reflected in our networks. If variable A causally influences variable B and induces synchrony in B, then A and B will typically be spatially coherent . If this effect occurs at both short and long timescales (so that mechanisms of synchrony are similar across timescales for these variables), there will be a corresponding link in both networks; if it occurs only for one timescale band (reflecting different mechanisms across timescales), there will be a link in one network and not the other. Second, joint effects whereby an unmeasured variable induces synchrony in each of two or more of our measured variables will also be reflected in our networks. In that case, the measured variables should typically be spatially coherent for the same timescale band(s) for which the unmeasured variable acts on both of them. So they will be correspondingly linked in one or both of our networks, again demonstrating dissimilar or similar mechanisms of synchrony, respectively. Thus, although networks of significant coherences do not imply that synchrony has been transmitted causally along all network links (joint effects can occur), differences between networks can reveal dissimilarities in the propagation of synchrony among variables and in the potential mechanisms of synchrony on short vs. long timescales. Type I and II errors, which are inevitable in the networks we construct, will result in inevitable differences between long-and shorttimescale networks. Differences between networks that result only from these types of error cannot result in major network differences, which, if they occur, will indicate differences by timescale in the potential mechanisms of synchrony. So, we will compare short-and long-timescale networks to detect major differences in their structure as opposed to small, random differences that may have come about through Types I or II statistical errors.

Q2a: Geography of synchrony at different timescales
The geography of synchrony for a variable sampled at 16 locations can be studied in detail using 16 × 16 matrices, the ijth entries of which measure the strength of synchrony between sampling locations i and j. This approach (Haynes et al. 2013;Pandit et al. 2016) retains all geographic information about synchrony, because all pairwise synchrony values are retained, in contrast to traditional distance-decay approaches for which values are effectively averaged across all pairs of sites separated by similar distances. Our strategy for comparing shortand long-timescale geographies of synchrony was to compute synchrony matrices separately for each timescale band and then to compare how similar or different the two matrices are using Mantel tests. A significant Mantel correspondence between the matrices indicated similarity of geography of synchrony across timescales. Nonsignificant Mantel results indicated that geographies were no more similar across timescales then expected by chance.
We represented the timescale-specific strength of synchrony between time series in locations i and j using the real part of a power-normalized cross-wavelet transform (Grinsted et al. 2004). This quantity, which is a real-valued function of timescale, is a timescale-specific analogue of the correlation of two time series (Supporting Information File S5). High values correspond to a consistently close-to-zero phase difference between oscillations. Values range from −1 to 1. Entries of the short-timescale (respectively, long-timescale) synchrony matrix were averages across short (respectively, long) timescales of the real part of the power-normalized cross-wavelet transform.

Q2b: Mechanism of geography of synchrony at different timescales
We used matrix regression (Lichstein 2007;Haynes et al. 2013) with recently developed model selection tools based on cross-validation and resampling techniques (Walter et al. 2017) to examine potential mechanisms of geographic structure in synchrony at different timescales. Recall that these are distinct from but can be related to mechanisms of synchrony itself (see Introduction and Walter et al. 2017). Briefly, short-and long-timescale synchrony matrices of Q2a (previous paragraph) were used as response variables in matrix regressions, carried out separately for each timescale and for each plankton variable, with predictors including matrices for a variety of possible mechanisms of geography of synchrony. Matrix regression model selection techniques indicate different levels of support for different predictor matrices and therefore for different potential mechanisms, thereby answering our question of which mechanisms are important at each timescale, though only to the extent that causality can be determined using observational datasets. This model-selection approach to inference of mechanisms in ecology is now standard (Burnham and Anderson 2002;Clark and Gelfand 2006). We consider three of the mechanisms (abbreviated M) of geography of synchrony listed by Walter et al. (2017): geography in the synchrony of Moran drivers (M mor ), spatial heterogeneity in density-dependent processes (M dens ), and spatial structure in dispersal (M disp ). These are labeled mechanisms A, B, and D in Walter et al. (2017). We also included Euclidean distance between sampling locations, because locations that are farther apart are expected to be less synchronous for multiple reasons.
Predictor matrices for M mor were short-and long-timescale synchrony matrices for environmental drivers. As for the plankton variables, short-and long-timescale synchrony between two locations in an environmental variable (e.g., temperature, pH, and nutrients) was quantified using the real part of the power-normalized cross-wavelet transform (Supporting Information File S5). Creating separate short-and long-timescale synchrony matrices for each of our 12 environmental variables would have saturated model selection procedures with too many predictors, resulting in overly complex models and infeasible computation times. Instead, and following Haynes et al. (2013), we did a principal components analysis of our abiotic variables, where PC1 and PC2 explained 52% of the variation in abiotic conditions (Supporting Information File S6). We then created short-and long-timescale synchrony matrices for the time series corresponding to PC1 and PC2. Suitable spatially sampled data for lake discharge and crappie abundance were not available and thus were not used in these analyses.
For each taxon, a matrix M dens serving as a rough surrogate for differences in the density dependence affecting the taxon at the different sampling locations was generated. We produced this dissimilarity matrix for each zooplankton abundance variable and for phytoplankton biomass using two metrics which were computed for each variable and sampling location. First, we averaged each log(x + 1)-transformed abundance time series over time for each site and plankton taxonthis may correspond to the carrying capacity of the site for the taxon, and as such may reflect density dependence. Second, we estimated the slope from a regression of log(x + 1) t + 1 abundance against log(x + 1) t abundance for each time series and taxon, a measure which has classically been used in studies of density dependence (Pollard et al. 1987). For each taxon, the 16 pairs of values obtained for the 16 sampling locations were arrayed in a two-dimensional Euclidean space, and M dens was generated for the taxon by computing distances in this space between pairs of points. This joint measure was considered a surrogate for differences in density dependence between sites, a factor known to influence levels of synchrony (Liebhold et al. 2006;Walter et al. 2017). Differences in density dependence might be expected due to spatial heterogeneity in resources or other niche requirements; there are many possible reasons for spatial variation in density dependence, details are discussed by Liebhold et al. (2006).
We represented M disp by generating a matrix of hypothesized dispersal likelihoods between sites, based on our qualitative knowledge of large-scale water flow patterns in the lake (Supporting Information File S7, Fig. G1). The main dispersal assumptions we followed in constructing the matrix were as follows. We assigned a 0 (easy dispersal) to the ijth matrix entry when site j was downstream of i and both locations were in the main channel and when i and j were both within the embayment. We assigned 1 (medium dispersal) for i and j such that plankton would have to move from the side of the lake into the main channel to get to j, a downstream site on the same side of the lake. Upstream and cross-channel movements were typically assigned a 2 (difficult dispersal). The matrix is displayed in Supporting Information File S7, Table G1. The matrix constructed in this way was not symmetric, but our synchrony matrices were symmetric, and synchrony likely depends more on the average strength of dispersal in both directions between two sites than on dispersal in either of the directions taken singly. So, we averaged the matrix constructed above with its transpose to give a matrix which hypothesizes overall between-site dispersal connectivity, with high values corresponding to low connectivity (Supporting Information File S7, Table G2), and the notation M disp refers to the transpose-averaged matrix. Thus, pairs for which dispersal was easy in one direction (0, e.g., upstream to downstream in the main channel) and hard in the other (2) have an intermediate final matrix entry value (1, the mean of 0 and 2), whereas sites on opposite sides of the main channel, which were characterized by difficult dispersal (2) in both directions even for sites directly across the lake from each other, received a final matrix entry value of 2. A previous analysis (Anderson et al. 2018) found that our proxy matrix for dispersal connectivity was still an important predictor of spatial patterns of synchrony, even after controlling for the patterns of spatial synchrony of numerous environmental variables, providing indirect support that the matrix sufficiently characterizes broad movement patterns in KY Lake.
The spatial dispersal patterns reflected in M disp comprise only a hypothesis about dominant movement tendencies, based on knowledge of system hydrology. Therefore, our approach can provide evidence rather than certainty as to whether dispersal is or is not a mechanism of synchrony. However, in systems such as Kentucky Lake which have clear and strong habitat structure, hypotheses about dominant dispersal patterns are reasonable, and detailed measurements of dispersal of the kind that would be needed for an empirical connectivity matrix are rare. So the evidence our approach can provide with regard to dispersal as a mechanism of synchrony and its geography is valuable. Reasonable hypotheses about habitat structure and dispersal pathways such as ours are likely to be available in many systems, e.g., in riverine, marine, and structured terrestrial systems (Bunnell et al. 2010;Powney et al. 2012). We revisit this topic in the Discussion section. We conducted separate model comparisons for each zooplankton species and phytoplankton biomass at each timescale. For each taxon and timescale, we compared models containing all combinations of the five different predictor variables (M mor of PC1, M mor of PC2, M dens , M disp and Euclidean distance; N = 31 models considered for each response matrix). We did not test for synchrony across timescales (e.g., longtimescale M mor of PC1 synchrony was not considered as a predictor of short-timescale zooplankton synchrony); this was because Defriez and Reuman (2017a) and Defriez and Reuman (2017b) found that Moran transmission of synchrony appeared commonly within long-and short-timescale bands but not between them. To perform model selection, we could not use standard approaches based on the Akaike or Bayesian Information Criteria (AIC, BIC) because matrix regression is based on resampling (Lichstein 2007) instead of likelihood, and AIC and BIC ordinarily are for likelihood approaches (Burnham and Anderson 2002). We instead employed a recently developed model selection framework for matrix regressions using leave-n-out cross validation and randomization procedures (Supporting Information File S8; Walter et al. 2017), similar to other recently used approaches (Peterman et al. 2014;Kastens 2015). This method generates model rankings and model weights for each model, and predictor importance weights for each predictor, similar to those produced using standard AIC/BIC approaches. Model weights range between 0 and 1 and quantify the relative support of data for each considered model. Predictor importance weights, which are also between 0 and 1, quantify the importance of each predictor matrix for explaining the geography of synchrony encoded in a response variable synchrony matrix. Additional details of model selection are in Supporting Information File S8.
Model selection approaches reveal the relative merits of models but do not provide information on whether any of the models considered are objectively adequate. For this purpose, we also computed the coefficient of determination (R 2 ) of the top-ranked model for each response matrix, and evaluated that model in a standard hypothesis-testing approach (Lichstein 2007) to determine whether included variables were significant. Because matrix regressions typically result in R 2 values that are less than standard multiple regression models (Mortelliti et al. 2015), we concluded that none of the models we considered were adequate only when the best model according to the model selection procedure had R 2 < 0.1. When predictors of highest importance weight were not significant in the top model, we concluded that model selection with the available variables had not revealed the mechanisms of geography of synchrony. We determined signs of regression coefficients for top models to see whether variables had a positive or negative effect on synchrony.

Q1a: Timescale structure of synchrony
Wavelet phasor mean field magnitude plots showed time and timescale structure in synchrony for phytoplankton biomass and for abundance of all zooplankton species (Fig. 4). For most zooplankton species and phytoplankton biomass, synchrony was present for some times and timescales (red coloration/higher values in Fig. 4) and absent or undetectable for others (blue coloration/lower values in Fig. 4), showing timescale structure that was substantially variable across the taxa examined. For instance, the majority of L. kindtii synchrony occurred at timescales around 4 yr, whereas the majority of phytoplankton synchrony peaked around both 4 and 8 yr. (peaks in Avg Sync, sides of each panel in Fig. 4). In some Fig. 5. Spatial coherence results indicate probable differences by timescale of mechanisms of synchrony in Kentucky Lake. Connections between variables indicate they are significantly coherent for (a) short timescales (< 4 yr) and (b) long timescales (> 4 yr). Arrowheads point in the hypothesized direction of the effect. Arrowheads are hypotheses based on biological reasoning, are not provided by the spatial coherence technique, and are not used in assessing similarities and differences between the two networks. Lines with double-ended arrows indicate potential bidirectional effects. Black lines indicate approximately in-phase relationships, blue lines indicate quarter-cycle relationships, and red lines indicate antiphase relationships (see Methods section). Precise p-values and phases are tabulated in Supporting Information File S10. All p < 0.001 would remain significant after p-value correction for multiple testing. Far more tests were significant than expected from the type I error rate (Supporting Information File S10).
cases, taxa were consistently synchronous across all or almost all timescales (e.g., D. retrocurva). The results were visually very similar using an averaging technique similar to geometric mean (see Methods section) as opposed to the arithmetic mean, indicating outlier data points were unlikely to be the cause of the absences of synchrony for some times and timescales. Nearly all abiotic variables also showed timescalespecific synchrony (Supporting Information File S9). Average strength of synchrony for abiotic variables was greater than that of plankton variables, as is typically the case.

Q1b: Timescale-dependent drivers of synchrony
Because the Kentucky Lake ecosystem is complex and our data are observational, our spatial coherence networks reflect rather than reveal precise causal pathways by which synchrony may cascade through the Kentucky Lake food web (see Methods section). Nevertheless, our networks differ substantially by timescale, evidence that the mechanisms of synchrony are likely to differ by timescale (Fig. 5). For instance, all but one significant spatial coherence between abiotic and biotic variables occurred at long timescales. And only 7 of the 21 significant spatial coherences between plankton variables we observed were shared across timescales (Supporting Information File S10). Though Type I and Type II errors may have produced or prevented some arrows on Fig. 5, error rates would not have been sufficient to account for differences of this magnitude between the long-and short-timescale networks. Thus, we can say with some confidence that there are differences between these networks, and probably therefore also differences between the mechanisms of long-and shorttimescale synchrony. Spatial coherence networks also had some similar features across timescales. Spatial coherences between plankton variables were more often significant than were coherences between abiotic and biotic variables: 13 of 28 possible plankton-plankton connections were significant for short timescales, whereas there was only 1 significant of 112 possible connections between plankton and abiotic variables or crappie; for long timescales, there were 8 of 28 plankton-plankton connections and 10 of 112 possible connections between plankton and abiotic variables and crappie ( Fig. 5; Supporting Information File S10). When separating out variables into approximate trophic levels (phytoplankton, herbivorous zooplankton, omnivorous/predaceous zooplankton, crappie), the highest percentage of potential connections between levels that materialized as significant spatial coherences, for both timescales, was between predaceous zooplankton to other predaceous zooplankton (three out of three potential connections at each timescale). The majority of the strongest (p < 0.001) coherences were between predaceous zooplankton and herbivorous zooplankton for both timescale bands (Fig. 5). At both long and short timescales, among the strongest spatial coherence values were those between Calanoida and D. retrocurva, Cyclopoida and D. retrocurva, and Cyclopoida and B. longirostris; and these relationships were all approximately in-phase ( Fig. 5; Supporting Information File S10). Among plankton variables, 11 of the 13 significant coherences at short timescales and all eight coherences at long timescales were approximately in-phase. Of the 11 significant spatial coherences between abiotic and plankton variables across both timescales, 5 of the phase relationships were in-phase, 5 were quarter-cycle phase relationships, and 1 was antiphase. There were no significant coherences with the top predator (crappie) for either short or long timescales. All significant spatial coherences between taxa were between a taxon in one trophic level and another in the same level or in a level immediately above or below the first one-there was no coherence which "skipped" a trophic level.

Q2a: Geography of synchrony at different timescales
There was a positive association between short-and longtimescale synchrony matrices that was significant for six of the eight taxa when using a 5% significance threshold, and five of the eight taxa using a 1% threshold (mean Mantel's r for taxa with p < 0.05 = 0.44; Table 2). The positive relationships indicate that long-and short-timescale geographies of synchrony were similar for most taxa. The nonsignificant case(s) indicate that long-and short-timescale geographies of synchrony can also be no more related than expected by chance.

Q2b: Mechanisms of geography of synchrony at different timescales
Using model selection (see Methods section), multiple models of synchrony matrices received at least some support for each zooplankton taxon and for phytoplankton biomass, as is often the case in model selection. Models with the highest weights tended frequently to have only a single predictor, especially at long timescales (Table 3; Supporting Information File S11).
There was a fairly high degree of consistency across timescales in the potential mechanisms of the geography of synchrony that were supported by data, according to our model selection approaches. This result is consistent with the result that geographies of synchrony were related across timescales in six or seven of our plankton taxa ( Table 2). The predictor variable with the highest weight was the same across long and short timescales for seven of eight taxa, the exception being L. kindtii, one of the two taxa that did not show significant correlation between long-and short-timescale synchrony matrices at 95% confidence level (Table 2). For the other taxon with nonsignificant correlation in Table 2, B. longirostris, our model selection methods supported M dens as the most important predictor of the geography of synchrony for short timescales, but our models were not adequate and model selection did not reveal the mechanisms of the geography of synchrony (see Methods section) for long timescales. These results may reflect differences across timescales in mechanisms of the geography of synchrony for B. longirostris if the mechanisms that truly operated at long timescales for B. longirostris were mechanisms other than the ones we considered. In sum, our results seem to support the conclusion that geographies of synchrony and their mechanisms were similar across timescales for most of the taxa we considered, but probably differed for one or two of these taxa.
As may be expected given the physiological and ecological variety represented by our eight taxa, mechanisms of the geography of synchrony supported by our analyses were diverse. They encompassed all of the mechanisms we investigated, though perhaps with M disp and M dens better represented than M mor . Considering only taxa for which the top model had R 2 > 0.1 (Mortelliti et al. 2015), spatial heterogeneity in synchrony of an environmental variable (M mor ), spatial heterogeneity in density dependence (M dens ), and spatial heterogeneity in dispersal (M disp ) were all included in top-ranked models at both short and long timescales. Euclidean distance was only included as a predictor in a top-ranked model in two cases, indicating that when predictors that actually represent mechanisms are included in models, Euclidean distance, which is only a correlated quantity with several geographic mechanisms and does not directly represent any single mechanism, can cease to be well supported by data. Each variable also exhibited the expected relationship to synchrony (positive for M mor , and negative for M dens and M disp ; Table 3), supporting the reasonableness of our method of determining mechanisms. M disp and M dens most commonly had the highest predictor importance weight, and most commonly appeared as predictors in top-ranked models. Table 3. Importance weights (see Methods section) for each predictor of the spatial structure of synchrony for each of the seven zooplankton taxa and phytoplankton biomass. Bold values indicate the predictor with the highest importance weight, and italicized values indicate variables included in the top model for that response. The superscript * indicates that more than one top model were equally supported (Supporting Information File S11) for that taxon; in those cases, only information for the model with a higher R 2 is shown. Symbols in parentheses indicate the sign of the coefficient (significantly positive [+], significantly negative [−], or not significant [ns]) of the predictor(s) in the top model. R 2 values are for the top model. Predictors are spatial heterogeneity in dispersal (M disp ), spatial heterogeneity in density-dependent processes (M dens ), geography in the synchrony of PC1 (M mor -PC1) and PC2 (M mor -PC2), and Euclidean distance. The top model for six of the eight taxa explained relatively high amounts of variation at short timescales (R 2 > 0.1) and five of the eight taxa at long timescales. The average across taxa of the R 2 of the top model was higher at short (mean AE SD, R 2 = 0.24 AE 0.11) compared to long timescales (mean AE SD, R 2 = 0.19 AE 0.07). Models of D. retrocurva had the highest R 2 at short timescales, whereas the model R 2 for phytoplankton biomass was highest at long timescales.

Discussion
We found geographic and timescale structure in the patterns of spatial synchrony of phytoplankton biomass and zooplankton abundance in Kentucky Lake. We found that many zooplankton taxa and phytoplankton biomass showed synchrony between locations that varied by timescale. Furthermore, spatial coherence networks among variables (see Methods section) differed by timescale, providing evidence that mechanisms of synchrony differed by timescale. Spatial coherences between taxa were strongly significant, while coherence between plankton and environmental Moran drivers were weakly (long timescales) or not (short timescales) supported. Geographic structure in synchrony was related across timescales for most of the taxa investigated, and potential mechanisms of the geographic structure of synchrony were also found generally to be similar across timescales. One or two taxa seemed to show timescale differences in the geography of synchrony and its mechanisms. We furthermore observed that data-supported drivers of the geography of synchrony were diverse across taxa, including three of the four mechanisms of geography of synchrony described by Walter et al. (2017) that we were able to investigate with our data and techniques, though M dens and M disp appeared disproportionately active as mechanisms in this system. The mechanisms we examined tended to better explain spatial patterns of synchrony than did Euclidean distance (Table 3), a standard covariate used in studies of synchrony. Our understanding of synchrony and its mechanisms can be substantially improved by examining geographic and timescale structure.

Mechanisms of synchrony
In answering Q1b and Q2b, we observed some support for two of the three main hypothesized mechanisms of spatial synchrony: environmental Moran effects and dispersal. Moran effects were revealed through both geographic approaches and spatial coherence. Geography of synchrony in the environment (PC1 and PC2) was included in the top model of the geographic variability in synchrony, for two of the eight taxa considered using matrix regression, at each timescale; and associations between environmental and population synchrony matrices were positive in each case (Table 3). This indicated that greater environmental synchrony between locations was associated with greater population synchrony between those locations and strongly suggests that Moran effects were contributing to the synchrony itself in these cases. We remind the reader here that determining mechanisms of synchrony and mechanisms of its geography are conceptually distinct but related ideas (see Introduction), and the argument here is essentially that the geography of the synchrony of an environmental variable is unlikely to be significantly positively associated with the geography of the synchrony of a population variable (as we observed, Table 3) unless the environmental variable helps drive the synchrony in the population variable through Moran effects. Multiple variables contributed to PC1 and PC2 (Supporting Information File S7), leaving some ambiguity as to which abiotic variable(s) produced the effect. Individual putative Moran effect variables (e.g., temperature, pH, or nutrient levels) were weakly supported as mechanisms of synchrony on their own at either timescale using spatial coherence techniques (Fig. 5). Overall, the two analytical methods combined provide consistent results regarding the strength of Moran effects: matrix modeling showed M mor was never the most important mechanism and was only included in top models along with other variables; and individual putative Moran effects were only weakly supported by spatial coherence tests. Dispersal may have been a more important potential mechanism of synchrony than Moran effects in this system.
Dispersal as a mechanism of synchrony has been demonstrated less frequently in natural systems than other mechanisms, at least partly due to difficulties in quantifying dispersal rates (Nathan 2001). However, in theory and controlled experiments, dispersal has been shown to strongly increase synchrony (e.g., Ims and Andreassen 2005;Fox et al. 2011). Our characterization of dispersal was indirect, but was a reasonable hypothesis for this system. Dispersal in our system is almost certainly dominated by water movement, as opposed to typical between-lake dispersal pathways, mediated by animals or air-borne movements (Havel and Shurin 2004). Our hypothesized dispersal matrix was highly consistent with synchrony patterns on both timescale bands for calanoid and cyclopoid copeods and D. lumholtzi. Such results would be unlikely if our characterization of dispersal did not correspond approximately with actual dispersal. Dispersal probably is an important mechanism of synchrony in this system, and spatial heterogeneity in dispersal is probably a mechanism for geographies of synchrony. The approach of constructing hypothesized dispersal or connectivity matrices and relating them to synchrony has been successfully applied for other taxa (Bunnell et al. 2010;Powney et al. 2011;Powney et al. 2012).
We found timescale structure in synchrony and mechanisms of synchrony, but generally found consistency across timescales in the geography of synchrony and its mechanisms. This was not an inconsistency in our results; it makes sense in light of the distinction between mechanisms of synchrony and of its geography mentioned in the Introduction. A mechanism of the geography of synchrony, such as geographic variation in density dependence (M dens ), can overlay and produce similar geographic variation in synchrony that was caused in different ways. Geographies of synchrony can also differ, even when the underlying causes of synchrony itself are similar. For instance, calanoid copepods and D. retrocurva were strongly significantly coherent with each other at short and long timescales (Fig. 5), suggesting interactions between these taxa that couple their dynamics, and/or common external influences, both of which suggest similar mechanisms of synchrony for the species. However, the geographies of synchrony of these two taxa differed (Table 3). This can occur, even with a common external driver of synchrony, if some mechanism such as spatial variation in density dependence modified spatial patterns of synchrony differently in the two taxa. Table 3 suggests spatial variation in density dependence may have played a larger role in the geography of synchrony of D. retrocurva than it did for calanoid copepods, though future research is still needed to substantiate this hypothesis and may be an interesting avenue to pursue.

Interspecific interactions
We observed very strong, in-phase spatial coherence between several species of zooplankton (Fig. 5), indicating fluctuations between some pairs of species occurred similarly on some timescales. This result is similar to other studies on timescale-specific dynamics of zooplankton (Vasseur et al. 2014;Buttay et al. 2017). Spatial coherences were especially strong between herbivorous and predaceous zooplankton on both timescale bands. While it is possible that interactions among these species helped drive synchrony in some of them, it seems likely we would then have observed antiphase or quarter-phase relationships between some predaceous and herbivorous taxa (e.g., copepods and D. retrocurva). Instead, nearly all taxa fluctuated in phase with one another, making it seem more likely that a common, unmeasured environmental or biotic driver caused spatial coherence by simultaneously affecting all these taxa. These results are consistent with the observations of Vasseur et al. (2014) that synchronous dynamics of co-located zooplankton taxa are much more common than compensatory or otherwise phase-shifted dynamics.
An additional way of testing for possible species-interaction-based drivers of synchrony used crappie abundance. Age-1 crappie are likely to be less impacted by water currents than zooplankton, meaning they could, in principle, act as a mobile synchronizing predator. Crappie were not coherent with any zooplankton taxon. Age-1 crappie either do not substantially affect zooplankton abundance fluctuations or our measure of crappie abundance was a poor representation of predator abundance for zooplankton. Annualizing the data, as we did for reasons explained previously, also may have limited our ability to detect interspecific interactions (predation, competition, or herbivory) using our methods, as dynamical effects of these processes could be occurring principally at subannual timescales.

Spatial scales and elementary visualization of geography of synchrony
Our study shows that geography of synchrony can occur on small spatial scales. Haynes et al. (2013) found geographic variation in the synchrony of gypsy moths over a large section of the northeastern U.S.A., and Gouveia et al. (2016) observed geographic structure in patterns of synchrony in voles over the entire Czech Republic. Similarly, geographic variation in the synchrony of primary production has been described across the entire globe in both oceanic and terrestrial systems (Defriez and Reuman 2017a;Defriez and Reuman 2017b). Complementing these large-scale patterns, we found strong support for spatial variation in synchrony within a small part of a single lake spanning only~50 km 2 .
Geography of synchrony can be directly pictured, but not explained, via the classic and very commonly used nonparametric cross-correlation plot. The extreme residual variation that is essentially ubiquitous feature of such plots (e.g., Supporting Information File S1, Fig. A1) indicates that pairs of locations separated by similar distances show variable strengths of synchrony, indicating the potential for some geographic factor other than distance to explain synchrony. Synchrony was generally higher in Kentucky Lake between locations within limnetic groupings (groupings are inside Ledbetter embayment, in the embayment mouths on the west, in those on the east, and in the channel-colors in Fig. 2) than it was between sites in different limnetic groupings. And generally correlations were independent of geographic distance within a limnetic grouping (Supporting Information File S1, Fig. A1). As our results demonstrate, statistically testing for geographic structure can improve our understanding of synchrony and population dynamics (Walter et al. 2017).

Implications for monitoring
Our results have implications for long-term monitoring programs of freshwater systems. It has been proposed that monitoring a single location that acts as a "sentinel" within a lake or reservoir can provide an understanding of ecological dynamics in the whole body of water (Baines et al. 2000;Lansac-Tôha et al. 2008). Our results suggest otherwise, at least for reservoirs similar in size to our sampling area. Monitoring multiple locations that encompass flow gradients, bathymetric variability and diverse limnetic zones, as the KLMP was designed to do, should typically give a much more complete picture. This is particularly true for reservoirs, in which lakeriver hybrid characteristics (Thornton et al. 1990) increase the potential for habitat structure.
Additional methodological considerations should influence the design of future long-term monitoring projects so spatial synchrony and spatiotemporal dynamics generally can be better studied in plankton. For example, it is unclear how vertical migration patterns of plankton may impact analyses, such as ours, of the synchrony of plankton measured at a fixed depth. Levine et al. (2014) provide evidence that seasonal depth preferences of B. longirostris may have caused single-depth Schindler-Patalas-trap measurements such as ours to miss autumnal peaks of the species, though Levine et al. only examined 1 yr of the data analyzed here. Such difficulties could have impacted our results so that the synchrony patterns we describe are a conflation of patterns of population synchrony and synchrony of depth preference behavior of the species. Decoupling the influences of these two phenomena would require data that is vertically resolved, as well as spanning multiple locations, species, and decades. To the best of our knowledge, such a dataset does not exist for freshwater plankton. Examining the possible geographies of synchrony in three dimensions (one of depth and two horizontally) in aquatic systems would be a novel and potentially fascinating topic for future studies. Such studies may already be possible in the marine realm, where long-term sampling that is both horizontally and vertically resolved may be more common.
The spatial arrangement of our sampling sites within Kentucky Lake, established decades ago, may also aid or hinder the ability of analyses such as ours to detect some mechanisms of synchrony. In our system, the proximity of the Ledbetter sites to each other may not make them "distinct enough" to be considered "truly separate populations" by some ecologists. In our view, the degree to which sampling efforts in different locations reflect "separate populations" is a quantitative spectrum, rather than a binary, and is precisely one of the things being studied when one measures the synchrony between the populations. Nevertheless, documentation of actual dispersal or movement rates of organisms between sampling locations may help directly quantify the relationship between synchrony and the degree of population mixing (Michels et al. 2001), and would further illuminate the importance of dispersal.

Conclusion
Numerous studies have documented taxonomic variability in the strength of freshwater plankton spatial synchrony (Rusak et al. 1999;Lansac-Tôha et al. 2008;Rusak et al. 2008;Vogt et al. 2011;Seebens et al. 2013;Lodi et al. 2014;Pandit et al. 2016), but many previous studies were unable to explore multiple mechanisms, timescale structure, or geography of synchrony in explaining this taxonomic variability. We demonstrated that synchrony has an important timescale structure and a geography that will both complicate and enrich efforts to compare synchrony among species and to understand the causes and consequences of synchrony. We suggest that if future studies of synchrony use techniques such as we have demonstrated and take into account timescale and geographic patterns of synchrony, it will improve understanding. For example, traditional methods using simple correlations to measure synchrony in D. birgei and D. lumholtzi would suggest similar strengths of synchrony and no distance-decay relationships (Supporting Information File S1, Fig. A1); using our methods, we were able to identify distinct timescale-specific drivers of synchrony and its geography for these and other taxa.