can be easily observed at the scale of an ecosystem (i.e. local = ecosystem). Regional spatial variation describes the differences across regions in the regional-average ecosystem state, whereas local spatial variation describes variation of average ecosystem state across all ecosystems (Wagner et al. 2007). Temporal variation can be decomposed into two components: coherent temporal variation (also called spatial synchrony), which describes the variation in ecosystem state across all ecosystems that is synchronous through time (Walter et al. 2017;Shestakova et al. 2018); and local temporal variation, which describes temporal variation where ecosystems vary independently through time (Kincaid et al. 2004;Wagner et al. 2009;Vidal et al. 2017). Past studies have quantified some of these variance components and their drivers. However, no study has quantified all components of variation in a single analysis at macroscales, which is needed to provide a deeper understanding of the relationship between spatial and temporal variation and their drivers at broad scales.

The past extensive research on ecological variation of species distributions and abundances provides a rich starting point to study ecosystem material and energy at macroscales.

Table 1 describes key propositions (i.e. confirmed generalisations, sensu Scheiner & Willig 2011) of spatial and temporal variation and its controls that are well-accepted in ecology, labelled as 'assumed'. For example, spatial variation increases as the spatial extent of the study increases, and is controlled by different processes at different spatial extents. Temporal variation increases following perturbation and prior to regime shifts. Coherent variation can result from synchronous temporal responses to broad-scaled drivers, such as regional weather patterns; and local temporal variation can result from ecosystem-specific characteristics mediating how individual ecosystems respond to drivers over time.

Here, we test four propositions of variation in ecosystem properties (labelled as Pr1-Pr4, with supporting citations in Table 1) that are specific to the temporal and spatial scales of this studyannual estimates of ecosystem properties over two decades (temporal) and estimates of ecosystem properties in tens to hundreds of ecosystems nested within 63 regions at the macroscale (spatial). (Pr1) Spatial variation includes both local and regional components, which is based on past studies that have found variation at both local and regional scales. (Pr2)

Table 1 Propositions defined as confirmed generalisations (sensu Scheiner and Willig 2011) for ecological spatial and temporal variation at macroscales and the annual temporal scale. We include propositions that are most relevant to this study, those that have extensive evidence from the literature (labelled as 'assumed'), and those that have less support from the literature and that we evaluate in this study. Scale refers to both grain and extent. Macroscales are spatial extents at the range of sub-continent to continent (i.e. 1000's of km's) Proposition Evaluated in study?

Evidence to date (selected sources) Spatial variation at macroscales a. Observed patterns of spatial variation in ecosystem properties are dependent on the scale of observation Assumed Sawyer (1989) and Horne & Schneider (1995) b. As spatial extent increases from regional to macroscales, the total spatial variation of ecosystem properties increases Assumed Meentenmeyer & Box (1987), Wiens (1989) and Horne & Schneider (1995) c. Spatial variation in ecosystem properties includes both local and regional components Yes (Pr1) This study; Wickham et al. (2005), Park et al. (2010), Fergus et al. (2011), Cheruvelil et al. (2013), Read et al. (2015) and Lapierre et al. (2018) Temporal variation at macroscales d. Observed patterns of temporal variation in ecosystem properties are dependent on the scale of observation Assumed Horne & Schneider (1995) e. Coherent temporal variation in ecosystem properties among ecosystems should decrease as spatial extent increases, and so should be lower at macroscales than finer spatial extents Assumed Stoddard et al. (1998), Liebhold et al. (2004), Huttunen et al. (2014), O'Reilly et al. (2015), Lottig et al. (2017) and Oliver et al. (2017) f. Annual temporal variation is dominated by local (sitespecific) variation relative to coherent temporal variation. (i.e. it takes a large and sustained perturbation to move a system in a directional way) Yes (Pr2) This study; Odum et al. (1979), Underwood (1991), Kincaid et al. (2004), Lottig et al. (2017), Oliver et al. (2017) and Shestakova et al. (2018) Relative amounts of spatial and temporal variation at macroscales g. Spatial variation (local plus regional) exceeds annual temporal variation Yes (Pr3) This study; Larsen et al. (2001), Lottig & Carpenter (2012), Cheruvelil et al. (2013), Oliver et al. (2017) and Vidal et al. (2017) Controls of spatial and temporal variation at macroscales h. Similarity in the spatial structure of ecological properties and their drivers may indicate common spatial scaling and the possibility of linkages Assumed Grieg- Smith (1979), Schneider & Piatt (1986), Legendre & Fortin (1989), Wiens (1989), Horne & Schneider (1995), Broitman & Kinlan (2006) and Lapierre et al. (2018) i. Coherent temporal variation (i.e. synchrony) can be a result of broad-scale exogenous drivers, such as climate, and so can be thought of a measure of 'regional temporal variation' Assumed Liebhold et al. (2004), Rusak et al. (2008), Vogt et al. (2011), Walter et al. (2017) and Shestakova et al. (2018) j. Local temporal variation is a result of site-specific characteristics that influence an individual site's response to drivers through time

Annual temporal variation is dominated by local variation, which is based on studies that have suggested that large and sustained perturbations are required to move the average state of a system in a directional way, thus it is more likely that ecosystems will vary independently through time. (Pr3) Spatial variation exceeds annual temporal variation, which is based on studies that find variation in ecological properties across ecosystems or sites is larger than temporal variation of ecosystems or sites. (Pr4) Each component of spatial and temporal variation is controlled by different drivers, which is based on evidence from prior studies conducted at individual spatial or temporal scales or studies that consider some components of variation, but not others.

We evaluated evidence for these four propositions at macroscales using observations across both space and time in thousands of ecosystems (i.e. lakes) in a 1 800 000 km 2 study extent in the Northeastern and Midwestern US. This study extent includes a wide range of lakes and broad gradients in climate, geomorphology, hydrology and land use. We quantified four components of variation in five ecosystem properties of lakes and examined the drivers of that variation. We found support for all four propositions; and, our study informs future work to scale up local-scaled observations to regional and continental scales, to extrapolate from well-studied to unstudied ecosystems, and to design effective monitoring programs that considers both spatial and temporal variation.

We fit Bayesian hierarchical models to thousands of observations to: quantify four components of spatial and temporal variation as proportions of the total variation in ecosystem properties, quantify the effects of ecological drivers on each component of variation and calculate the proportion of variation explained by the drivers. We fit separate models for each ecosystem property. The four components of variation that we quantified are two spatial variance components (regional and local) and two temporal components (coherent and local). We also estimated a residual component that is variation unaccounted for by the other sources, including measurement errors. We selected candidate drivers a priori based on understanding about the controls on lake ecosystem properties for each type of variance component. Our model explicitly accommodates the hierarchical structure of our data, where observations are nested within lakes, lakes are nested within regions, and observations are also nested within years. The random effects in our models (i.e. the model error structure) help to accommodate this data structure and the dependencies (e.g. the lack of statistical independence) that may result from this structure.

The study area is a lake-rich zone of the US including 17 upper Midwest and Northeastern US states. It includes a total of 51 101 lakes and reservoirs ≥ 4 ha in a north temperate climate zone (Soranno et al. 2017). There is a wide range of ecological, geological, climatic, hydrologic and land use characteristics and a wide range of lake types (Fergus et al. 2017); lake watershed land-use ranges from all forested to almost all agricultural. There is a large regional variation in lake and geographic features (Wagner et al. 2011;Cheruvelil et al. 2013;Lapierre et al. 2018). In this study, we delineate the study area into 63 regions using 4-digit hydrologic units (U.S. Geological Survey 2013; hereafter, regions), which have been found to capture the above observed regional spatial variation similarly well compared to other regionalisation frameworks in a subset of our study lakes (Cheruvelil et al. 2013).

We used a curated lake ecosystem database called the Lake Multi-scaled Geospatial and Temporal Database of the Northeast US (LAGOS-NE; Soranno et al. 2017). Lake nutrient, primary productivity, and clarity data were obtained from LAGOS-NE LIMNO v1.087.1.

We analysed five common measures of lake ecosystem properties: two total nutrient measures (total phosphorus, TP and total nitrogen, TN); one dissolved nutrient measure (nitrate, NO 3 ); one surrogate for pelagic primary production in lakes (chlorophyll a concentration, Chl); and one measure of the light environment in lakes, water clarity, as measured by Secchi depth (clarity). The lakes were sampled at least once per year from 1990 to 2011 during the summer stratified season (15 June-15 September). This time period is the season of maximum primary and secondary production in lakes and the time of dominant ecological interactions in lakes. The number of lakes with any measurement for one of these five properties range from 3560 (NO 3 ) to 7601 (clarity) lakes (Fig. 1; Table S1). Lakes differed in the frequency of sampling through time, ranging from being sampled only one year to up to 22 years (Table S1). Our analytical approach is designed to accommodate such an imbalanced design to make full use of all available data (see below). Further, a prior study of this data set has shown that known biases, such as the fact that larger lakes are over-sampled relative to smaller lakes, did not substantially change statistical distributions of the lake response variables (Stanley et al. 2019). Nevertheless, biases are inherent in any data compilation and may influence the outcomes. However, with large sample sizes such as ours, the effect of such biases on outcomes is likely to be less than for smaller data sets.

Lake ecological driver data were obtained from LAGOS-NE GEO v1.05 (Soranno & Cheruvelil 2017), except for monthly climate data, which are available from Collins et al. (2018b). Data are described in Soranno et al. (2015Soranno et al. ( , 2017)). We selected potential ecological drivers based on evidence from past studies of lakes (Table 2). For example, because coherent temporal variation may be more sensitive to seasonal climate signals, we included seasonal climate drivers for coherent temporal variation, whereas, we included monthly climate drivers for local temporal variation. We selected candidate ecological drivers and their appropriate scales of quantification a priori based on ecological knowledge and evidence that a driver controls a component of variation, particularly from published studies (Collins et al. 2017;Lottig et al. 2017;Oliver et al. 2017;Lapierre et al. 2018). Table 2 describes four categories of ecological drivers to explain lake variation that were made up of multiple metrics, often measured at different scales suited to the variance components: climate, atmospheric deposition, hydrology and land use/cover. For climate data, because we were modelling interannual variation, we used monthly and seasonal data.

We used a Bayesian hierarchical model to decompose the total variation in each of the five ecosystem properties into their four components, as well as a residual. The model was as follows:

where y ijkl is the log e -transformed observation i, from lake j, in year k, and region l, for an ecosystem property response variable. The parameter b 0 is the intercept and µ 1j , µ 2k , µ 3jk , µ 4l are the random effect for each variance component and ǫ ijkl is the residual error, that are assumed normally distributed with a mean of zero and variance r 2 x . Diffuse priors were used for all parameters: b 0 ~N (0, 1000) and r x ~U(0, 5). Three parallel Markov chains were run, each starting with a random value. Each chain was run for 10 000 iterations, discarding the first 5000 samples. This resulted in 15 000 samples used to summarise posterior distributions. We assessed convergence using the Brooks-Gelman-Rubin statistic (Brooks & Gelman 1998) and trace plots. We estimated the posterior means and 95% credible intervals for all estimated sigma parameters (Fig. S1), and we examined the model fit by examining diagnostic plots and calculating the RMSE for each model (Fig. S2).

We included covariates in the variance decomposition model as follows:

where y ijkl and b 0 are as defined above and b z , a z , c z and w z are regression slope parameters for covariates (X) quantified at each spatial and temporal scale. Diffuse priors, as described above were used for b 0 and r x . We performed predictor variable selection using the horseshoe shrinkage prior on regression slope parameters (b z , a z , c z and w z ; Carvalho et al. 2010).

We used a half-Cauchy prior on both the slope parameterspecific and global (across all slope parameters) variance parameters for the horseshoe prior (slope parameter i |k i , s ~N (0, k 2 i s 2 ), where ki and s are the slope-specific and global variance parameters, respectively, and k i ~half-Cauchy(0, 1) and s ~half-Cauchy(0, 1). Predictor variables measured as proportions were logit-transformed and standardised, while continuous predictors were log e -transformed and standardised.

Three parallel Markov chains were run, each starting with a random value. Each was run for 80 000 iterations, discarding the first 60 000 samples. The remaining 60 000 samples were used to summarise posterior distributions. All models were fitted using the program JAGS (Plummer 2003) using the jags UI function (Kellner 2017) called from within the program R (R Core Team 2018), and are available at Wagner (2019). We assessed convergence using the Brooks-Gelman-Rubin statistic (Brooks & Gelman 1998) and trace plots.

We used several metrics to evaluate model fit and covariate importance. First, we examined plots of residuals versus fitted values, histograms of residuals, quantile-quantile plots and we calculated the RMSE for each model (Fig. S3). Second, we calculated the approximate percentage of the total variation explained (approximate R 2 ) for each spatial and temporal component, which was determined using the unconditional variances from eqn 1 and the conditional variance from eqn 2. For example, calculating the approximate R 2 of the local spatial level l was as follows:

ðÞ (Raudenbush & Bryk 2002). This approximate R 2 calculation provides an intuitive and simple measure of the variance explained; however, it is possible to obtain negative approximate R 2 if the conditional variance is larger than the variance of the data. This result may happen because the between-group variances are a function of variance that occurs at multiple levels of the model (LaHuis et al. 2014). Third, we determined covariate significance by evaluating whether or not the 90% credible interval of the coefficient overlapped with zero. Fourth, we conducted a 10-fold crossvalidation and report the plots of residuals versus fitted values, histograms of residuals, quantile-quantile plots and we calculated the RMSE for each model (Fig. S4).

Although we have included spatial (i.e. local, region) and temporal (i.e. year and local 9 year) random effects in our models that accommodate the hierarchical structure of these datai.e. the lack of statistical impendence that may exist across space and time-we acknowledge that because these models do not explicitly model spatiotemporal dependencies that may exist, posterior uncertainty may be underestimated.

Proposition 1: Spatial variation includes both local and regional components Spatial variation included both local and regional components that were roughly equal, in terms of proportion of the total variation (Fig. 2). On average, 31% of the total variation was regional, and 33% was local, with only small differences across ecosystem properties. For example, the more biological reactive NO 3 had the lowest regional variation (25%) as might be expected; whereas, TP and clarity had the highest (36% and 35%, respectively). Local spatial variation was more similar across ecosystem properties, although TN was the lowest (29%) and TP and clarity were the highest (35%).

Proposition 2: Annual temporal variation is dominated by local variation On average, the proportion of temporal variation was small (Fig. 2) and ranged from 2% (for TP and clarity) to 14% (for NO 3 ). Local temporal variation was made up from 10 to 14% of the total variation, which was 2-5 times higher than coherent temporal variation. These results suggest that at macroscales and for the inter-annual scale across two decades, temporal variation is ecosystem-specific rather than coherent.

Proposition 3: Spatial variation exceeds annual temporal variation For all five ecosystem properties, the proportion of combined spatial variation was four times larger than the proportion of combined temporal variation (Fig. 2). The proportion of total spatial variation ranged from a minimum of 58% of the total variation for NO 3 , to about 70% for TP and clarity. However, the proportion of total temporal variation ranged from a minimum of 11 and 13% for TP and clarity, to a maximum of 20% for NO 3 and TN. Residual variation was on average, 20% of the total variation across all of the ecosystem properties, but was highest for Chl and NO 3 and lowest for TP and clarity. The differences that we document are supported by an analysis of the posterior uncertainty of the differences among each pair (Fig. S5).

Proposition 4: Each component of spatial and temporal variation is controlled by different drivers For all ecosystem properties, one or more ecological drivers accounted for variation in all variance components except coherent temporal; however, the exception was for regional spatial variation for Chl and TN (Fig. 3). In addition, for the clarity model, we could not include drivers to model coherent temporal variation because there was so little coherent variation that adding predictors led to a conditional variance near zero that led to failed model convergence. This lack of drivers for coherent temporal variation was not surprising because it was the lowest variance component for all models. Although local temporal variation was also relatively small, five to seven climatic driver variables accounted for some of that variation, as shown by their small but non-zero effect sizes. For the two spatial components of variation, there were numerous ecological drivers that had non-zero effect sizes for the local spatial variation, and 1-2 drivers that had non-zero effect sizes for the regional spatial variation. It is beyond the scope of this paper to discuss the ecological implications of each driver, however, the sign of the effects and the drivers themselves are consistent with well-established ecological relationships. Uncertainty around the effects for the regional variation was substantially larger than for the local-scale variation, as were the number of non-zero coefficients, reflecting the smaller number of regions in the dataset compared to the number of ecosystems.

The percent variation explained by the spatial drivers was much larger than for the temporal drivers (Table 3). Most of the approximate R 2 values for the drivers of the two temporal components of variation were zero, meaning that the drivers could not explain any of the temporal variation. Exceptions to this result are TN and NO 3 . However, because the credible intervals overlap zero, we interpret these approximate R 2 values as zero. The approximate R 2 values for the drivers of regional spatial variation were larger (81-94%) than those for the drivers of local spatial variation (32-51%). Thus, regional variation was better explained by its spatial drivers despite the larger uncertainty of estimated effects and a fewer number of significant variables compared to the local scale.

For each of the five models, we also calculated an approximate total variation explained by the combined drivers by weighting the approximate R 2 values of each of the variance components (Table 3) by its proportion of total variance from Fig. 2 (and assuming the negative R 2 values are zero). Models for TP and clarity had the highest estimated total R 2 values and we were able to explain nearly half of the variation across thousands of lakes. These variables also have the lowest temporal components of variation. We were only able to explain 33, 36, 39% of the variation in NO 3 , Chl and TN, respectively. For Chl, it is likely that we explained less variation due to it being a biological variable, which has been shown to be subject to more extreme events (Batt et al. 2017), and so is likely more difficult to predict.

Our results help to fill knowledge gaps of the major sources and drivers of spatial and temporal variation in ecosystem properties at macroscales. The four propositions held across all ecosystem properties that we studied, although we observed some differences in dissolved nutrients and Chl compared to total nutrients, which were expected given the labile nature of dissolved nutrients, and the biological nature of Chl. Next, we consider the strength of support for the four propositions based on our results and past studies of ecological variation.

Proposition 1: Spatial variation includes both local and regional components Both local and regional spatial variation are relevant to macroscales. Some past studies emphasise local variation (Kincaid et al. 2004;Fraschetti et al. 2005;Wickham et al. 2005;Read et al. 2015); and others emphasise regional variation (Latham & Ricklefs 1993;Cheruvelil et al. 2013). However, spatial variation in ecosystem properties likely results from a combination of locally and regionally structured geographic drivers that induce both local-and regional-scale spatial structure (Park et al. 2010;dal Bello et al. 2017;Lapierre et al. 2018), an idea supported by our result of equally large local and regional variances. An important next step is to test this proposition in data sets that span increasingly broader spatial extents and that include a wider range of ecological properties from conservative chemical constituents to biotic properties (Kincaid et al. 2004;McGuire et al. 2014).

Our data support the proposition that most annual temporal variation at decadal and macroscales is local (i.e. ecosystemspecific) rather than coherent, as was also found in a study at a smaller spatial extent for fewer lakes (Kincaid et al. 2004). This result is also supported by a trend analysis of these data that showed the majority of lake nutrients showed no trends over two decades, nor did they change synchronously (Oliver et al. 2017). Despite the size of our data set, it is possible that our study lacked power to detect coherent temporal variation. For example, analyses of time series of air temperatures found that datasets that are 10-20 years long have small signal-tonoise ratios and so have trouble detecting trends (Santer et al. 2011). However, the measured amount of annual temporal variation in our data was low, a result confirmed by studies with adequate power (Santer et al. 2011). Thus, even with more data, any coherent variation may still be small and may not be ecologically relevant. Many recent studies conducted within individual regions have found evidence for synchrony in a wide range of biotic and abiotic ecological properties (Ranta et al. 1997;Rusak et al. 2008;Shestakova et al. 2016;Walter et al. 2017;Defriez & Reuman 2017). However, at increasing spatial extents, synchrony likely breaks down, and the importance of non-coherent variation likely increases. Although we detected low to no coherent variation in our study, we expect that it would be higher if we increased the temporal extent of our study to many decades or centuries. Thus, temporal coherent variation needs to be assessed specific to the spatial and temporal footprint of the study and measured at a wide range of scales to help illuminate underlying drivers of temporal variation (Walter et al. 2017;Shestakova et al. 2018).

Proposition 3: Spatial variation exceeds annual temporal variation

The proportion of total spatial variation was 2-5 times greater than the proportion of annual temporal variation. These patterns make sense given the large and heterogeneous study area and the fact that the minimum and maximum values for lake nutrients or productivity across thousands of lakes is generally larger than the minimum and maximum value in an individual lake through time. The variation we have measured at the annual scale is supported by evidence that ecological properties in most ecosystems fluctuate around a relatively stable average state and mostly deviate as a result of rare and large perturbations (Carpenter & Brock 2006;Fraterrigo & Rusak 2008). Due to data limitations, we were not able to quantify within-year variation in individual lakes, which could increase temporal variation. An important future research direction will be to quantify spatial versus temporal variation in macroscale studies that have temporal data that are either finer-resolution or for a longer temporal extent (i.e. many decades or centuries). We show that different ecological drivers were related to different components of spatial and temporal variation, supporting past theoretical and empirical work (Schneider & Piatt 1986;Wiens 1989;Horne & Schneider 1995). Despite the equal amounts of local and regional spatial variation, the ecological drivers explained a much larger proportion of regional versus local variation. This result could be due to our measures of regional drivers being better able to capture Figure 3 The estimated effects of driver variables on each of the four components of variation for five response variables. Filled symbols represent effect sizes whose 90% credible intervals (error bars) do not overlap zero. Symbols with very small effect sizes, but that are filled, do not overlap zero and only appear to do so due to the size of the symbol. Colours as are for Fig. 2.

ecologically relevant processes compared to our measures of local drivers that are either poor proxies of processes that we know matter (such as nutrient loading or runoff), or that were missing entirely because of lack of data at macroscales, such as lake water residence time, groundwater connections, biotic interactions or watershed soils. Further, although we could not explain much of the small amount of annual temporal variation across several decades in our study area, our results suggest that predicting temporal variation at the scale of decades and thousands of ecosystems will be challenging with the commonly available data sets at macroscales. In sum, we were able to explain the most regional spatial variation, followed by local spatial variation, and only some of the small amounts of coherent temporal variation. This type of analysis on the major types of ecological variation and their controls helps to identify gaps in our understanding of ecological variation across broad scales of space and time.

Biological properties at macroscales have been well-studied in biogeography and macroecology and scaling laws have been developed that describe controls on species and organisms at regional to macroscales. Ecosystem properties (i.e. materials and energy) have not been similarly studied at macroscales, but are equally important, particularly for quantifying the contribution of terrestrial and aquatic ecosystems to continental and global carbon cycles (DelSontro et al. 2018;Mitchard 2018;Seekell et al. 2018). Our study describes the importance of capturing adequate spatial variation in ecosystem properties, and that at increasingly broader spatial extents, spatial variation will likely be more important than temporal variation, at least at the temporal extent of several decades. Nevertheless, we were not able to include two important considerations in our study. First, our annual estimates are based on summer values, which is commonly done because it is the period of maximum productivity in most ecosystems, but also because of a lack of data from other seasons (Stanley et al. 2015;Hampton et al. 2017). Second, we were only able to study ecological variation in ecosystem states, as opposed to ecosystem rates or fluxes, which are critical for assessing ecosystem roles in global cycles. Unfortunately, data are also lacking for these measures at macroscales and until they are available, we can assume that spatial and temporal patterns that we have identified for ecosystem state will apply similarly to ecosystem flux.

Our ability to accurately capture spatial and temporal variation in the critical ecosystem properties influencing global cycles will depend on the accumulation of observations across broad gradients in space and time. Such observations will be compiled in large and complex databases, and there is evidence that ecology is increasingly becoming more dataintensive, particularly compared to its historical roots (Hampton et al. 2013;Peters et al. 2014;Elliott et al. 2016;Cheruvelil & Soranno 2018). There is an increasing availability of large publicly accessible data sets that include estimates of ecosystem materials and energy in a range of freshwater, marine and terrestrial systems across a range of spatial and temporal scales (O'Reilly et al. 2015;Henson et al. 2016;Soranno et al. 2017;Anderson-Teixeira et al. 2018;Smith et al. 2018). In our study, we used a modelling approach that combines spatial and temporal observations in a single framework, which could be applied to most, if not all of the above databases to further explore the potentially complex relationships among spatial and temporal variation in ecosystem materials and energy. This framework that standardises the measurement of variation into their component parts can allow for the comparison across vastly different ecosystem types, variables and spatial scales and can lend new insight into mechanisms explaining ecological variation from regional to global scales.

When, where and how ecological properties respond to anthropogenic stressors can often be obscured by ecological variation. Understanding the effect of broad-scale stressors such as climate change on ecological properties at macroscales is difficult because the system is multi-scaled through both space and time, yet we rarely consider multiple components of variation. Our study provides insight into both spatial and temporal variation in ecosystem properties at annual and macroscales. We support much past research that suggests that different drivers control ecological variation at different scales. We also provide evidence that at macroscales, spatial variation will exceed decadal interannual temporal variation. Thus, to effectively extrapolate findings from one site to another, or to scale up ecological properties to continental or global scales, it will be essential to collect observations across many ecological settings to capture the full range of spatial variation in ecosystem states and processes. Doing so for other properties and ecosystems will be especially important to improve our understanding of the role of ecological systems in global cycles, such as carbon and nitrogen, and how they are likely to respond to global change.

Table 3 The percentage of variation explained (approximate R 2 ) for each spatial and temporal component by the driver variables (posterior mean, followed by 95% credible intervals in parentheses). Negative % variance explained is possible with these models and they reflect poor predictive capability and are reported as zero. The estimated total variation is calculated by weighting the approximate R 2 values of each variance components by its proportion of total variance from Fig. 2, and assuming the negative approximate R 2 values are simply zero. For clarity, coherent temporal covariates were not included in the model due to convergence issues, so values of NA are reported Variance component Approximate R 2 TP Chl Clarity TN NO 3 Spatial, regional, % 86 81 90 94 85 (77, 93) (67, 90) (81, 95) (90, 97) (68, 93) Spatial, local, % 51 38 51 48 32 (44, 51) (34, 42) (48, 54) (43, 52) (26, 39) Temporal, coherent, % 0 0 NA 20 0 (À115, 70) Temporal, local, % 0 0 0 0 5 (À7, 16) Estimated total variation, %: 49 36 49 39 33 © 2019 John Wiley & Sons Ltd/CNRS Letter Ecological variation at macroscales 1595