A global geography of synchrony for terrestrial vegetation

geography of synchrony for marine phytoplankton, and used that geography to infer statistical environmental determinants of synchrony. Here we determine if terrestrial vegetation (measured by the Enhanced Vegetation Index, EVI) also shows a geography of synchrony, and we infer determinants of EVI synchrony. As vegetation is the basis of the terrestrial food web, changes in spatio-temporal vegetation dynamics may have major consequences. Methods: Synchrony in terrestrial vegetation is mapped globally. Spatial statistics and model selection are used to identify main statistical determinants of synchrony and of geographic patterns in synchrony. Results: The first main result is that there is a pronounced and previously unrecognized geography of synchrony for terrestrial vegetation. Some areas such as the Sahara and Southern Africa exhibited nearly perfect synchrony, whereas other areas such as the Pacific coast of South America showed very little synchrony. Spatial modeling provided the second main result that synchrony in temperature and precipitation were major determinants of synchrony in EVI, supporting the presences of dual global Moran effects. These effects depended on the timescales on which synchrony was assessed, providing our third main result: synchrony of EVI and its geography are timescale specific. Main conclusions: To our knowledge, this study is the first to document the geography of synchrony in terrestrial vegetation. We showed geographic variation in synchrony is pronounced. We used geographic patterns to identify determinants of synchrony. This study is one of very few studies to demonstrate two separate synchronous environmental variables driving synchrony simultaneously. The geography of synchrony is apparently a major phenomenon that has been little explored.


INTRODUCTION
Understanding the dynamics of terrestrial vegetation biomass and production is both an interesting and important topic: vegetation is constantly changing over a range of temporal and spatial scales (Ichii et al., 2002;Zhu et al., 2016); and it is a key component of the global carbon cycle (Beer et al., 2010;Wieder et al., 2015) that is tightly coupled with climate because it directly affects land-atmosphere heat and moisture fluxes (Meir et al., 2006). Of course vegetation is also the base of the terrestrial food web and hence the dynamics of vegetation biomass or production is implicated in essentially every area of pure and applied ecology. Therefore, changes in details of vegetation dynamics, including spatio-temporal aspects of those dynamics, may have far reaching consequences.
One incompletely understood aspect of spatio-temporal vegetation dynamics is their spatial synchrony. Spatial synchrony is the phenomenon whereby geographically separate population time series (or, in this context, vegetation biomass or production time series) fluctuate partly in 4 unison. Spatial synchrony has been observed even in populations separated by hundreds or thousands of kilometers (Post & Forchhammer, 2004;Liebhold et al., 2004), across a very wide variety of taxa including protists, insects, fish, birds, mammals, and many others (Hanski & Woiwod, 1993;Myers et al., 1997;Liebhold et al., 2004). One of the primary mechanisms that has been cited to account for synchrony is the presence of spatially synchronized environmental factors which drive population dynamics, thereby inducing synchrony in the populations. This is known as the Moran effect (Moran, 1953). The Moran effect is one of the main causes of synchrony (Lande et al., 1999;Liebhold et al., 2004), but historically it was difficult to convincingly show that Moran effects operate in specific scenarios and to identify the environmental drivers (Liebhold et al., 2004;Abbott, 2007). This was partly because the historically most common statistical descriptors of synchrony can show similar patterns for Moran effects and for other causes of synchrony (Ranta et al., 1999;Abbott, 2007).
We previously mapped global geographic variation in patterns of synchrony in ocean phytoplankton (Defriez & Reuman, In review). We thereby provided evidence that a Moran effect, operating through synchronized sea surface temperatures or through synchrony of an environmental variable highly correlated with sea surface temperature such as nutrient availability, and possibly acting through complex oceanographic mechanisms, is a major driver of phytoplankton synchrony globally. The main goal of the present study is to apply the same statistical techniques to map the geography of synchrony in terrestrial vegetation, and then to infer determinants of synchrony in vegetation. Temperature and precipitation are two important climatic variables affecting productivity (Nemani et al., 2015;Clinton et al., 2014);and Koenig (2002) found synchrony in both of these factors over large spatial scales (up to 5000km).
Vegetation dynamics may be similarly synchronized. Shetakova et al. (2016) investigated tree growth in two contrasting forest biomes and found that large scale synchrony responded to 5 climate warming. But the detailed geography of synchrony in terrestrial vegetation and the extent to which that geography reflects geographic patterns of synchrony in temperature and precipitation variables are unknown.
Multiple methods have been used to describe patterns of spatial synchrony, ranging from standard methods based on correlation coefficients (Hanski & Woiwod, 1993;Bjørnstad & Falck, 2001) to spectral and wavelet methods (Grenfell et al., 2001;Vasseur & Gaedke, 2007;Keitt, 2008;Sheppard et al., 2015) and matrix regressions and others (Haynes et al., 2013). Many of the prior studies that have statistically illuminated the processes driving synchrony (e.g., Sheppard et al., 2015;Shetakova et al., 2016) have used temporally extensive data sets. Here, and following Defriez & Reuman (In review), we adopt a different approach, taking advantage of the unparalleled spatial coverage (but limited temporal extent) provided by remotely sensed data. We use data describing geographic patterns and variation in vegetation, globally, as measured through the enhanced vegetation index (EVI), process them so as to quantify and map the phenomenon of synchrony, and then make comparisons to geographic patterns in potential causal factors of synchrony. Drivers of synchrony should have statistically similar spatial patterns to the geographic patterns of synchrony revealed by the EVI data. We use spatial linear models to compare geographic patterns in EVI synchrony with putative drivers.
A key feature of several of the studies cited in the previous paragraph is their attention to the timescale dependence of population dynamics in synchrony, through the use of spectral methods.
We also use spectral methods. Spectral methods allow the decomposition of synchrony into the frequencies, or timescales (timescale here indicating the reciprocal of frequency) at which it occurs, thereby showing which frequencies contribute most to synchrony. Synchrony on one timescale can be independent of synchrony on another timescale, and this independence can obscure analysis of synchrony by correlation based methods (Figure 1; Keitt, 2008;Sheppard 6 et al., 2015;Defriez et al., 2016). In addition, synchrony on longer timescales may be more important than short-timescale synchrony because it is more likely to affect longer-lived consumers (Sheppard et al., 2015;Defriez et al., 2016). We believe timescale-specific approaches to synchrony are under-applied, and that this has limited our understanding of the causes and consequences of synchrony (Sheppard et al., 2015;Defriez et al., 2016).
Past researchers have typically considered Moran effects resulting from just one environmental driver at a time (Batchelder et al., 2012;Sheppard et al., 2015). However, it is possible, in principle, to have two or more distinct simultaneous Moran drivers of synchrony.
These multiple Moran effects may, a priori, reinforce or counteract each other. Here we investigate the possibility that Moran effects from both land surface temperature and precipitation environmental drivers are simultaneously important for the synchrony of vegetation. We also include other possible covariates of synchrony: average EVI density (areas with more vegetation may, a priori, exhibit systematically more or less EVI synchrony), average temperature, average precipitation, average altitude, extent of variation in altitude, latitude, and average wind speed.
The main questions asked here are: Q1) What regions of the terrestrial realm exhibit high degrees of regional synchrony in vegetation density (as measured by EVI), and what areas exhibit low synchrony? Q2) When regional synchrony is decomposed into long-and short-timescale components, do maps differ in their main features, i.e., is the geography of synchrony timescale-specific? Q3) What are the main statistical determinants of synchrony in vegetation density, as inferred from its geography? Do determinants of long-and short-timescale synchrony differ? Q4) Do patterns of synchrony and statistical determinants of synchrony differ between major land masses? And finally Q5) Is there evidence for dual Moran effects contributing to synchrony in vegetation density? We hypothesize that: H1) Moran effects will be detectable via our approach, and will comprise some of the major determinants of vegetation synchrony and its geography. Our past research (Sheppard et al., 2015;Defriez et al., 2016) indicates timescale-specific structure is a common feature of synchrony, so we also hypothesize that H2) Moran effects will be timescale-specific, and therefore the geography of synchrony will be timescale specific. H2 would be supported if short-timescale components of environmental synchrony are spatially associated with short-but not long-timescale components of vegetation synchrony, or if long-timescale components of environmental synchrony are spatially associated with long-but not short-timescale components of vegetation synchrony. This paper is the first explicit exploration we are aware of into the detailed geography of synchrony in terrestrial vegetation, and is also the first time the geography of synchrony has been used to infer determinants of synchrony in the terrestrial vegetation.

Data
The Enhanced Vegetation Index (EVI) and land temperature emissivity data sets from the  (Huete et al., 2002;Sims et al., 2015;Sjöström et al., 2009). Compared to the Normalized Difference Vegetation Index (NDVI) it minimizes the confounding effects of soil background, atmosphere and canopy density (Huete et al., 1997(Huete et al., , 2002Xiao et al., 2003). The pixels designated 'lowest quality' or below in the quality assurance flags and any estimated pixels (pixel reliability 4) were removed and the data were converted to 1 ∘ resolution by taking averages. Data from the Aqua and Terra satellites were averaged.
Daily precipitation values for the same time period were obtained from the Global Precipitation Climatology Project (GPCP) (Huffman et al., 2001) as a 1 ∘ daily dataset compiled from satellite observations and rain gauge measurements. The daily precipitation data were averaged to monthly time series.
Annual time series for EVI, land surface temperature (LST) and precipitation were obtained from monthly time series for each 1 ∘ grid cell as long as there were no more than 5 missing months in an individual time series and no more than 1 missing month from any individual year.
Otherwise the grid cell was omitted. Before annualising, missing months were replaced by the average value for that month for that time series. Time series were linearly detrended and their mean was subtracted before further analysis.
Altitude data were downloaded from the British Oceanographic Data Centre  Kalnay et al. (1996)). Data were provided as velocity components, u representing the east-west component and v the north-south component of wind. Annual data were downloaded at 2.5 ∘ resolution over the time period 2001 to 2014 and were re-gridded to 1 ∘ resolution. Some 1 ∘ cells fell entirely within a 2.5 ∘ cell and were given the value of the larger cell. For 1 ∘ cells which crossed more than one 2.5 ∘ cell an average was taken. For each 1 ∘ cell an average of all annual values was calculated for the period of study.
Wind speed was calculated from the u and v components with the formula u 2 +v 2 . All datasets were downloaded September 2015.

Correlation-based synchrony
For EVI, temperature and precipitation (separately), for each 1 ∘ grid cell, Spearman's correlation was calculated between the annual time series of that focal cell and time series of each of the other cells within a 500 km radius. These values were averaged to produce a synchrony value for the focal cell. Global maps of the strength of synchrony (out to 500 km) were thereby produced for EVI, temperature and precipitation. These spatial variables will be referred to as EVI synchrony, temperature/temp synchrony, and precipitation/precip synchrony in subsequent spatial modeling. Because of the averaging over values for cells within 500 km of the focal cell, our synchrony maps are maps of regional synchrony. Justification for the specific choice of 500 km is in Appendix S1 in Supporting Information. We computed synchrony maps for distance bands 500-1000 km and 1000-1500 km to test sensitivity of patterns to the choice of distance band.
Spearman's correlation was used as not all time series were normally distributed and no single transformation was able to normalize all time series simultaneously.

High-and low-frequency synchrony
Following ( To obtain the normalized cospectrum of two time series, one starts with their cospectrum (a standard method, see Brillinger, 2001), and normalizes by dividing by the geometric mean of the variances of the time series. Because the integral of the cospectrum of two time series is their covariance (Brillinger, 2001), this normalization ensures that the integral of the normalized cospectrum of the time series is their Pearson correlation coefficient. To make the integral equal the Spearman correlation used in the previous section, time series values were replaced by ranks prior to calculating the normalized cospectrum.
For EVI, temperature and precipitation (separately), for each 1 ∘ grid cell, normalized cospectra were calculated between the rank time series of that focal cell and the rank time series of each other cell within 500km. The normalized cospectra were then integrated over 'high' (0.25-0.5 cycles/year) frequencies (timescales of 2 years to 4 years) and 'low' (0-0.25 cycles/year) frequencies (timescales exceeding 4 years) and average values for each frequency band within 500km were then computed. The resulting six spatial variables will be referred to as high-(respectively, low-) frequency EVI synchrony, high-(respectively, low-) frequency temperature/temp synchrony, and high-(respectively, low-) frequency precipitation/precip synchrony. Because the integral of the normalized cospectrum over the whole range of frequencies (0-0.5 cycles/year) equals the correlation, high-and low-frequency synchrony values for a cell for a given variable summed to the all-frequency values of the previous section. The terms total or all-frequencies EVI/temp/precip synchrony will sometimes be used to refer to the variables of the previous section, to specifically contrast them with frequency-specific quantities.
Our use of the normalized cospectrum is described with formulas in Appendix S2, the dividing frequency 0.25 cycles/year is justified in Appendix S3, and a method is described in Appendix S4 for producing 95% confidence thresholds for values on maps of synchrony.

Statistical modeling
Spatial linear models were run with: 1) EVI synchrony as the response variable and temperature and/or precipitation synchrony potentially included as explanatory variables (among other potential explanatory variables, see below); and 2) high-and 3) low-frequency EVI synchrony as the response variable and high-and low-frequency temperature and precipitation synchrony potentially included (among others, see below) as explanatory variables instead of total temperature and precipitation synchrony. Total, high-, and low-frequency EVI synchrony were linearly mapped from the interval -1 to 1 onto the interval 0 to 1 and then logit transformed before models were run so that they were not limited by -1 and 1. Four types of model were considered. First, standard linear models y=X 1 β+ε were considered, where y is the response variable (a column matrix with one entry for each grid cell for which data were available) and X 1 is the standard design matrix, with columns containing explanatory variables. The specific columns depended on the explanatory variables included in a given model (see below). Second, models with spatially "lagged" variables y=X 1 β+WX 2 δ+ε were considered (LeSage, 2014). If N is the number of grid cells for which data were available (this is the length of the column y) then W is an N×N matrix of weights encoding the geographic neighbourhoods of each location. The term WX 2 δ represents neighbourhood, or spatially "lagged" effects of the explanatory variables in X 2 on the response variable, y, and the estimated parameters δ represent the strengths of these effects for each explanatory variable in X 2 . Third, the spatial error model y=X 1 β+u was considered; and fourth, the spatial Durbin error model y=X 1 β+WX 2 δ+u (Elhorst, 2010) was considered, where for both models the equation u=λWu+ε implicitly defines u, and u represents spatially autocorrelated residuals. See Appendix S5 for further details on the models.
Variables which could enter in a model are listed in abbreviation in table 1 and are explained here. Temporally averaged EVI (abbreviation 'EVI') in the focal cell, as well as temporally averaged land surface temperature ('temp') and precipitation ('precip') variables were used, as was average altitude within the focal cell (log transformed as described in the Data section, abbreviation 'log(altitude)'), standard deviation of log(altitude) values over grid cells within the 500 km radius disk centered at the focal cell ('sd log(altitude)') and average wind speed within a focal cell ('wind speed'). The strengths of synchrony of temperature and precipitation ('temp synchrony', 'precip synchrony'), as computed above, were also used to test for Moran effects.

13
Either total or low-and high-frequency versions of these variables were used, as described above.
Variables could enter as local effects or, in models with lagged variables could enter as neighbourhood effects, the latter specified by 'lag' in table 1.
The variance inflation factor (VIF) was used to test for collinearity. VIF indicated only negligible collinearity among predictor variables (table S1 in Supporting Information), following the recommendations of Dormann et al. (2013) for assessing collinearity.
To limit the number of models fitted and because we sought main determinants of synchrony, the total number of variables allowed in any one model was restricted to 5 or fewer. For models with spatially lagged variables and for spatial Durbin error models, explanatory variables were allowed to enter the models either as part of X 1 (local effects) only, or as part of both X 1 and X 2 (both local and neighbourhood/lagged effects). Neighbourhood effects without local effects (i.e., putting a variable in X 2 but not X 1 ) were not considered.
For each model, the Bayesian information criterion (BIC) was calculated and BIC weights were computed to determine the top 95% confidence set of models. BIC was used instead of the Akaike information criterion (AIC) because it is known that AIC tends to favor more complex models (Burnham & Anderson, 2004) and we sought the main determinants of synchrony.
Importance of a given variable as a predictor of EVI synchrony (total, low-or high-frequency) was measured by summing BIC weights across models that included the variable. Signs of model-averaged coefficients were computed (Burnham & Anderson, 2004). Residual plots were produced for top models to check that model assumptions were reasonable.
Spatial statistical models were fitted and model selection was performed using the global data.
But if determinants of synchrony differed markedly by continent, the global analysis would represent an average of different processes. To diagnose whether important differences occurred 14 (Q4 from the Introduction), analyses were also run separately for Eurasia, Africa, North America, South America and Oceania. See Appendix S6 for methodological details.
BIC approaches can identify which of several models is best supported by the data (Burnham & Anderson, 2004), but do not, on their own, indicate whether any of the models was objectively good. We provided three pieces of information to assess model fit. First, we computed R 2 , the fraction of spatial variation in EVI synchrony explained by the model. Second, we compared, via BIC and ANOVA, our best (lowest BIC) model from each model comparison exercise to the null model y=X 1 β+u, for X 1 a column of 1s and u spatially autocorrelated residuals. All analyses were done using R v3.1.3. Spatial models were fitted using the R package spdep (Bivand & Piras, 2015). and 1000-1500km ( Figure S1 in Supporting Information), but the geographic variation in strength of synchrony at these distances was not as great, as may be expected because averages are computed across more space.

Results
Answering Q2, geographic patterns of synchrony were strongly frequency specific (Figures 2c and 2d), with areas that were strongly synchronized at low frequencies often differing from (though sometimes being near to) areas that were strongly synchronized at high frequencies.
For example, synchrony in the Sahara and in China predominantly occurred at low frequencies.
Synchrony on the Atlantic coast of Brazil was primarily at low frequencies, but synchrony inland from the coast had a larger high-frequency component.
Both globally and for continental regions, spatial linear models of synchrony were fitted and ranked by BIC weight (tables S2-S7), variable importance tables were generated by summing BIC weights (tables 1, S8-S12), signs of model coefficients were tabulated (tables S13-S19) and model-averaged predictions were generated (Figure 2b, S2, S3). Figure 2b shows model-averaged total EVI synchrony as predicted by the top 95% confidence set of global models (by BIC weight). The models generally identified the areas of highest synchrony such as the Sahara, Southern Africa, and Eastern Brazil (compare to Figure 2a). However, predicted synchrony was often lower than observed synchrony when observed synchrony was high, and was higher than observed synchrony when observed synchrony was low: total geographic variation in the strength of synchrony was underestimated by the model. High synchrony in Northern Europe was also not captured by predictions. Although model predictions roughly captured large-scale trends, they were poor at representing fine spatial structure in synchrony. The best model explained 29% of the variation in total EVI synchrony ( ANOVA revealed a highly significant difference between these models (p<0.001).
Models of frequency-specific synchrony explained less of the variation than models for all frequencies, as may be expected given the greater variability intrinsic to estimates of frequency-specific quantities (table S2). But BIC comparisons and ANOVA p-values nevertheless indicated highly meaningful differences between best models and the null model (Appendix S9).
Synchrony in temperature and precipitation are both consistently highly important in statistically explaining all-, high-and low-frequency EVI synchrony ( Results show that both temperature and precipitation effects are frequency specific, confirming H2 from the Introduction. High-frequency precipitation synchrony was a more important predictor of high-frequency EVI synchrony than it was of low-frequency EVI synchrony. Low-frequency temperature synchrony was a more important predictor of low-frequency EVI synchrony than it was of high-frequency EVI synchrony, and likewise for precipitation (table 1).
Given that both synchrony in temperature and synchrony in precipitation were important determinants of total EVI synchrony, it was possible, a priori, that interaction effects were present between the two. As there was only one supported model in the global, all-frequencies analysis (table S2) we compared that model to a model that had interaction effects between temperature and precipitation synchrony variables, but was otherwise the same. The BIC value of the model with interactions was -616 compared to -626 for the best model with no interactions, so there was no evidence for interaction effects.
Also important and also having a positive association with EVI synchrony were mean temperature and log altitude (table 1). These were important predictors for all frequencies and for high and low frequencies, separately. Average EVI itself was only an important predictor of high-frequency EVI synchrony, further demonstrating the frequency specificity of synchrony and its determinants.
Results for the continent-specific analyses (tables 2, S8-S12) were similar to global results, but with some heterogeneity in some determinants of synchrony, answering Q4 from the

Discussion
Our principle result is that there is an important geography of synchrony of terrestrial vegetation, globally. Regional (500 km) synchrony varied enormously, from areas with almost no synchrony to areas with near-perfect synchrony. Dual apparent Moran effects of precipitation and temperature were the major statistical determinants of EVI synchrony and its geographic patterns worldwide, with some variation among continents in the details of these effects and in the importance of other determinants of synchrony. The geography and determinants of synchrony were strongly frequency specific. Inferences leading to our conclusions about likely Moran effects and frequency-specificity of synchrony were effective because they exploited the geography of synchrony. Moran effects from a single environmental driver have been reported (e.g., Batchelder et al., 2012;Sheppard et al., 2015), and studies have combined multiple drivers using PCA (Haynes et al., 2013). However, as far as we are aware, this is the first or one of few studies to provide evidence for two distinct separate synchronous environmental variables driving ecological synchrony in concert. Additional observations on effects of altitude on synchrony and the partial coincidence of high-synchrony areas with arid regions are in Appendix S7.
This work complements Defriez & Reuman (In review), in which we analyzed the main statistical determinants of synchrony in ocean chlorophyll. The key similarities between synchronies of chlorophyll and terrestrial vegetation are that the geography of synchrony was pronounced in both contexts, and Moran effects were statistically supported and frequency specific. There were also differences. Chlorophyll synchrony is highest in areas where chlorophyll itself is low; but EVI often has a positive association with EVI synchrony.
It is well known that temperature and precipitation can be important factors in vegetation dynamics (Ichii et al., 2002;Piao et al., 2006;Clinton et al., 2014). However, our new result that geographies of synchrony in temperature and precipitation are important correlates of the geography of EVI synchrony, while true, does not follow automatically from the earlier knowledge. Defriez & Reuman (In review) found that the geography of synchrony of incident solar irradiance was not statistically related to the geography of synchrony in chlorophyll-a abundance in the world's oceans, even though it is well known that irradiance can be an important factor in chlorophyll-a dynamics. In complex ecosystems, known importance of an environmental factor for ecological dynamics does not necessarily mean the factor produces a Moran effect or a geography of synchrony, for several potential reasons, including: the possibility that the environmental factor itself has a muted geography of synchrony; and the possibility that the factor's influence on vegetation is nonlinear or complex in some important way that varies geographically (see, e.g., Defriez & Reuman, In review, where mechanisms are discussed for the phytoplankton case). The following paragraphs explore some of these complexities further.
Non-monotonic relationships are increasingly important in ecology (Zhang et al., 2015). The carrying capacity of vegetation varies globally due to spatial variation in soils and other factors, and therefore the nature of density dependence in vegetation dynamics also varies spatially. Liebhold et al. (2006) describe a reduction in environmentally caused spatial synchrony resulting from geographic variation in density dependent dynamics. The effect can mediate or produce a geography of population synchrony, and is related to the concepts of the previous paragraph. If spatial variation in density dependence over a region is pronounced, the geography of ecological synchrony need not match the geography of synchrony of an environmental driver.
Although temperature and precipitation are key drivers of primary productivity, which of these is more important can differ among biomes (Nemani et al., 2015). For example, precipitation is often found to be more tightly coupled with vegetation and productivity in arid zones (Zhang, 2005;Fabricante et al., 2009). Although it is difficult to tell from the analyses of this study whether both synchrony in temperature and synchrony in precipitation are driving EVI synchrony everywhere across a given continent or if their relative importance varies, future work might be able to illuminate this question. There is substantial heterogeneity globally in the correlation between EVI and precipitation and temperature intra-annually (Clinton et al., 2014).
Furthermore, in Africa, according to our results, there are two areas of strong synchrony: the Sahara and an area in Southern Africa. Synchrony in the Sahara is predominantly at low frequencies ( Figure 2d) whereas synchrony in Southern Africa is predominantly at high frequencies ( Figure 2c). Both high-frequency precipitation synchrony and high-frequency 21 temperature synchrony were important determinants of high-frequency EVI synchrony in Africa; but only low-frequency temperature synchrony, and not low-frequency precipitation synchrony, was an important determinant of low-frequency EVI synchrony. This suggests that in Southern Africa synchrony in both temperature and precipitation drive synchrony in vegetation density but in the Sahara only synchrony in temperature drives synchrony in vegetation density, despite the tight coupling with precipitation often found in arid regions (Zhang, 2005;Fabricante et al., 2009). The wavelet methods of Sheppard et al. (2015) can illuminate what is causing synchrony with no need to rely on geographic variation in synchrony. Although our time series are probably too short for their wavelet analyses, a Fourier version of the techniques of Sheppard et al. (2015), applied separately to different regions of Africa (and elsewhere) might identify precisely how temperature and precipitation trade off against each other in relative importance as Moran drivers.
We used land surface temperatures (LST) as opposed to surface air temperatures (SAT, measured at 1.5-2m above ground level) because LST is available from satellite measurements and SAT is measured at discrete weather stations. Data products that provide SAT estimates globally based on interpolation between weather stations are available, but interpolation creates artefactual synchrony, so those products could not be used. LST can be affected by vegetation cover and condition (through evapotranspiration), and by soil wetness and therefore precipitation (Jin & Dickinson, 2010;Mildrexler et al., 2011), so causal relationships between our temperature and EVI synchrony variables may be complex and are incompletely illuminated by our correlative statistical approach. It is possible that temperature synchrony causes EVI synchrony through Moran effects, or that EVI synchrony causes temperature synchrony through vegetation effects on LST, or, most probably in our opinion, that the factors jointly affect each other or the causal relationship itself varies geographically. Our models establish statistical determinants of synchrony, a necessary but not sufficient condition for a causal effect of temperature synchrony on EVI synchrony. SAT may be a better variable for improving causal understanding of synchrony in future research, though a modified statistical approach would be needed to account for the fact that direct SAT measurements are only available at weather stations.
Changes in synchrony through time have received recent attention (Sheppard et al., 2015;Defriez et al., 2016;Shetakova et al., 2016;Koenig & Liebhold, 2016), raising the possibility of changes in synchrony superimposed on the geography of synchrony or affecting this geography itself. Shetakova et al. (2016) describe increasing synchrony through time in spatio-temporal forest tree growth data. They concluded that observed increases are not due to increasing synchrony of Moran drivers, but instead to stronger synchronizing influence of the drivers: as climate becomes more extreme and impacts in climatic variation therefore become more widely influential on ecosystems, climate has a more synchronizing influence even when climate need not itself have become more spatially synchronous. In our analysis we did not address the possibility of changes in synchrony with the passage of time -our time series are probably too short to do so. Although we removed trends before calculating synchrony we found mean temperature was an important determinant of synchrony across all frequencies, and has a positive effect. The results of Shetakova et al. (2016) raise the possibility that longer term increases (or decreases) in synchrony may be superimposed on top of the spatial patterns we found. Changes in synchrony may also interact with and modify spatial patterns of synchrony. To examine this possibility one would need data that are extensive both spatially and temporally. Synchrony is ecologically important, as metapopulations displaying increased spatial synchrony have an increased risk of extinction (Heino et al., 1997). It has also been proposed that an increase in spatial correlation may occur in environments before a regime shift (Dakos et al., 2010). As such, our study may act as a baseline for examining future changes in synchrony and its geography, in addition to being one of the first to document the geography of synchrony at all.

Acknowledgments
We Data underpinning this work is available freely online at locations indicated within the methods.
The date data were downloaded for this work is also indicated within the methods.
7 Supporting Information Short Captions S1) Justification for 500 km disks: On the choice to use 500 km disks.
S2) The normalized cospectrum: Some details of the method.  frequencies, from the 95% confidence set of models. c) Observed high-frequency EVI synchrony (timescales of 2-4 years). d) Observed low-frequency synchrony (timescales exceeding 4 years).
The white horizontal lines on the colour legends in a, c, d are approximate 95% significance thresholds compared to a null hypothesis of no synchrony (Appendix S4). Projection is equidistant cylindrical.