A global geography of synchrony for marine phytoplankton

Aim 
Spatial synchrony in plankton is imperfectly understood yet may have far-reaching implications, for example for carbon export to the deep ocean. Several techniques have been used to describe patterns of spatial synchrony, from correlation coefficients to spectral methods. Some studies have used temporally extensive data sets to identify causes of synchrony. This study instead uses the exceptional spatial extent provided by remotely sensed data to describe, for the first time as far as we know, geographical patterns of synchrony in marine phytoplankton. We use these patterns to illuminate drivers of synchrony and of its geography. 
 
Location 
The oceans. 
 
Time period 
2003–2015. 
 
Major taxon 
Chlorophyll a-containing phytoplankton. 
 
Methods 
Synchrony in chlorophyll a concentrations is mapped globally. Spatial statistics and model selection are used to illuminate main statistical determinants of synchrony and of geographical patterns in synchrony. 
 
Results 
The first main result is that there is a pronounced and previously unmapped geography of synchrony for phytoplankton. For instance, synchrony was highest in the open ocean, specifically in gyres, and lowest in coastal regions. Spatial modelling provided the second main result that synchrony in sea surface temperature (SST) was a major statistical determinant of chlorophyll synchrony in both the Pacific and Atlantic Oceans, indicating a strong Moran effect, although possibly an indirect and/or complex one. In the Pacific Ocean, this effect depended on the time-scales on which synchrony was assessed, providing our third result, which is that synchrony of phytoplankton and its geography can be time-scale specific. Synchrony of surface solar irradiance was not associated with synchrony of chlorophyll. 
 
Main conclusions 
To our knowledge, this study is the first to map geography of synchrony in marine plankton. We showed that this geography is pronounced. Geographical patterns illuminated determinants of synchrony. The geography of synchrony is a major phenomenon that has been little explored.

Although synchrony is well examined in plankton and other organisms, aspects of synchrony are still poorly understood. For instance, spatial variation in the strength of synchrony, here called 'geography' of synchrony, has seldom been studied directly. The main goal of the present study was to undertake a new exploration and statistical description of geographical variation in patterns of phytoplankton synchrony in the ocean. A benefit of this exploration is that details of the geography of synchrony can be used to illuminate determinants of phytoplankton synchrony and to relate synchrony to extensive existing literature on drivers of plankton dynamics in the ocean (Williams & Follows, 2011).
Past common statistical descriptors of synchrony typically cannot reveal the causes of synchrony. The most commonly used past methods plot correlations between pairs of population time series in different locations versus distances between the locations (Bjørnstad & Falck, 2001;Koenig, 2006). Such plots indicate the strength of synchrony and how it declines with increasing distance between populations. Where the pattern of this decline for the populations closely resembles the analogous pattern for a putative environmental cause of the population synchrony, it has been argued that a Moran effect is supported (Myers, Yule, & Jones, 2015). This is a limited approach, however, as declines in correlation with distance resulting from all three main mechanisms of synchrony (listed above) can be similar (Ranta, Kaitala, & Lindstrom, 1999;Kendall, Bjørnstad, Bascompte, Keitt, & Fagan, 2000;Abbott, 2007), and even when a Moran effect is the only synchronizing mechanism the decline in environmental synchrony with distance need not clearly match the decline of population synchrony (Abbott, 2007). The basic problem is that the decline of correlation with distance is an insufficiently detailed description of synchrony.
Recent work has developed new methods for obtaining more detailed descriptions of synchrony, which have sometimes revealed mechanisms of synchrony. For example, Haynes, Bjørnstad, Allstadt, and Liebhold (2013) used matrix regression to calculate the relative importance of several factors potentially synchronizing outbreaks of gypsy moths. The spline correlogram method of Bjørnstad and Falck (2001) was recently extended to enable synchrony to be calculated as a function of distance while correcting for a second co-variate (for example, altitude or climatic variables; Gouveia, Bjørnstad, & Tkadlec, 2016). Using this method, Gouveia et al. (2016) identified lower local synchrony in common vole populations in more highly forested regions.
Spectral methods have been used to reveal synchronous dynamics at various temporal scales within plankton populations (Vasseur & Gaedke, 2007;Keitt, 2008) and to map virus spread (Grenfell, Bjørnstad, & Kappey, 2001;Viboud et al., 2006). Grenfell et al. (2001) identified travelling waves of measles outbreaks from focal start points using correlograms and wavelet analysis. Using wavelets, Sheppard et al. (2015) discovered that spatial synchrony of aphid first flights in the U.K. is caused by synchrony of winter temperatures, and that changes in winter temperature synchrony attributable to changes in the North Atlantic Oscillation (NAO) caused major changes in aphid synchrony around 1993. All the above techniques allow for a more detailed description of synchrony than is possible using plots of correlation decline with distance between sites. Among other benefits, detailed descriptions can facilitate inference of causes of synchrony.
A key feature of several of the studies cited in the previous para-  Keitt, 2008;Sheppard et al., 2015;Defriez et al., 2016). We believe that time-scale-specific approaches to synchrony are under-applied, limiting our understanding (Sheppard et al., 2015;Defriez et al., 2016).
Our approach here is similar to the above approaches (Cazelles & Stone, 2003;Keitt, 2008;Haynes et al., 2013;Sheppard et al., 2015) in providing detailed descriptions of synchrony, but our approach is novel in that the description is geographical. We provide information on the geography of synchrony in phytoplankton in the world's oceans and make comparisons to the geography of potential causal factors. If the Moran effect is the dominant cause of synchrony, we may expect geographical patterns of phytoplankton synchrony and synchrony in the operating environmental variables to be spatially similar. Dominant factors affecting phytoplankton population growth and which may cause synchrony through Moran effects include light, nutrients and temperature. We use geographical variables representing sea surface temperature (SST), SST synchrony, irradiance, irradiance synchrony, currents and wind; and we consider covariates including average chlorophyll abundance and depth. Our descriptions of synchrony are time-scale specific. Our methods take advantage of spatially extensive but temporally limited remote sensing data. Our approach is descriptive and correlative, but much is known about what drives plankton dynamics, and one of our goals is that our statistical results will help to extend existing understanding about the drivers of plankton dynamics to apply to improved understanding of the drivers of synchrony.
Our two main candidate environmental determinants of phytoplankton synchrony (SST and irradiance, selected in part because of spatially extensive data availability) are strongly associated with each other and with nutrient availability (Marañ on, Cermeño, Latasa, & Tadonl ek e, 2012; Agawin, Duarte, & Agusti, 2000), and these three factors can affect plankton dynamics in direct and indirect ways and via different mechanisms in different parts of the ocean (Williams & Follows, 2011); thus it is unlikely that our correlative approach will produce precise and authoritative statements about Moran-effect causes of synchrony. For instance, any apparent SST effects may be attributable to real effects of nutrients and light, coupled with associations between SST and nutrients, and between SST and mixed-layer depth, which affects nutrient and light availability. Nevertheless, much is known about how drivers of plankton dynamics vary in importance spatially, and it should be possible to link patterns that emerge from our analysis of synchrony to that literature. We explore some initial links in the Discussion.   (d), misleadingly suggesting lack of important synchronous phenomena. Note that normalized cospectra (e,f; see Methods) reveal that positive synchrony at one frequency is masked by negative synchrony at the other. In practice, exact cancelation is unlikely, but asynchrony at some frequencies may nevertheless strongly conceal important synchrony at other frequencies DEFRIEZ AND REUMAN | 3 ocean basin? Much research has demonstrated the importance of Moran effects, so we hypothesize that (H1) Moran effects will be detectable via our approach and will comprise some of the main determinants of phytoplankton synchrony and its geography. Our past research (Sheppard et al., 2015;Defriez et al., 2016) indicates that time-scale-specific structure is a common feature of synchrony, so we also hypothesize that (H2) Moran effects will be time-scale specific, and therefore, the geography of synchrony will be time-scale specific.
H2 would be supported if short-time-scale components of environmental synchrony are spatially associated with short-but not longtime-scale components of phytoplankton synchrony, or if long-timescale components of environmental synchrony are spatially associated with long-but not short-time-scale components of phytoplankton synchrony. This study is the first explicit exploration we are aware of into the detailed geography of synchrony in marine phytoplankton.

| Data
Chlorophyll a, SST and photosynthetically active radiation (irradiance) data were used from the Sea-viewing Wide Field-of-View Sensor (SeaWIFS) and MODerate resolution Imaging Spectroradiometer (MODIS) Aqua projects available at http://oceancolor.gsfc.nasa.gov/. Monthly data were downloaded for January 2003 to December 2015 from level 3 standard mapped images. Data were downloaded on a 9 km resolution (1/68) and remapped to 18 resolution by averaging over all pixels within a 18 cell.
Annual time series for chlorophyll, SST and irradiance were obtained from monthly time series for each 18 grid cell as long as there were no more than five missing months in an individual time series and no more than one missing month from any individual year. Otherwise, the grid cell was omitted. Before annualizing, 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. Owing to cloud cover and light limitations during winter months, data at high latitudes are limited, and this analysis is predominantly done on data between 508 N and 508 S. Also because of missing data, the analysis did not include a few other areas, such as the northern Indian Ocean and an area off the western coast of Africa.
Bathymetry data were obtained from the British Oceanographic Data Centre (http://www.bodc.ac.uk/data/online_delivery/gebco/) at one arc minute resolution, and were regridded to 18 resolution. Averages were taken of all cells below sea level that fell within a given 18 cell. gov/psd/). Data were provided as velocity components, with 'u' referring to the east-west component and 'v' referring to the north-south component. Annual data were downloaded at 2.58 resolution over the time period 2003-15 and were re-gridded to 18. Where a 18 cell fell entirely within a 2.58 cell, the value was used. Where a 18 cell crossed more than one 2.58 cell, an average was taken. For each 18 cell, an average value was calculated for the period of study. All data were downloaded on 5 April 2016. Only time series that were complete after the above processing steps were used in our analyses. Time series were linearly detrended and their mean was subtracted before computing synchrony.

| Correlation-based synchrony
For chlorophyll, SST and irradiance (separately), for each 18 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 thus produced for SST, chlorophyll and irradiance.
Spearman correlation was used as not all time series were normally distributed and no single transformation was able to normalize them simultaneously. Synchrony was primarily investigated over 0-500 km, but analogous maps of synchrony were produced for distances 500-1000 and 1000-1500 km for comparison. For descriptive purposes and because such plots are standard, the shape of the decline of synchrony with distance was plotted using nonparametric correlation functions (Bjørnstad & Falck, 2001), separately for chlorophyll a, SST and irradiance in each ocean basin, using the Sncf function in the ncf package in R. All analyses were done using R v3.1.3. 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 before calculating the normalized cospectrum.

| High-and low-frequency synchrony
For chlorophyll, SST and irradiance (separately), for each 18 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 500 km. The normalized cospectra were then integrated over 'high' (0.25-0.5 cycles/year) frequencies (time-scales of 2-4 years) and 'low' (0-0.25 cycles/year) frequencies (time-scales exceeding 4 years), and average values for each frequency band within 500 km were then computed. Via these definitions, high-and low-frequency synchrony for a cell summed to the all-frequencies synchrony computed in the previous section. Our use of the normalized cospectrum is described with formulas in Appendix S1, the dividing frequency 0.25 cycles/year is justified in Appendix S2, and a method is described in Appendix S3 for producing 95% confidence thresholds for values on maps of synchrony.

| Statistical modelling
Our overall approach was one of model selection based on an information criterion. We fitted multiple spatial linear models of chlorophyll synchrony with data and used information criteria to determine the relative support for each model. Different combinations of predictors were used in different models. The predictors that appeared in bestsupported models were considered the most important determinants of chlorophyll synchrony and geographical patterns therein, among the predictors for which data were available. This is a standard approach (Burnham & Anderson, 2002).

Models were run separately for the Pacific, Atlantic and Indian
Oceans. For each ocean, spatial linear models were run with the following: (a) chlorophyll synchrony as the response variable and with SST synchrony and irradiance synchrony potentially included as explanatory variables; (b) high-frequency chlorophyll synchrony as the response variable and with high-and low-frequency SST synchrony and highand low-frequency irradiance synchrony potentially included as explanatory variables instead of total (all frequencies) SST and irradiance synchrony; and (c) low-frequency chlorophyll synchrony as the response variable and again with high-and low-frequency SST and irradiance synchrony potentially included as explanatory variables instead of total SST and irradiance synchrony. Chlorophyll synchrony, and high-and low-frequency chlorophyll synchrony were logit transformed (values were first linearly mapped from the interval 21 to 1 onto the interval 0 to 1) before models were run so that they were not limited by 21 and 1. In total, globally, 22123 18 spatial grid cells were used. Four types of model were considered. First, standard linear models y5X 1 b1E 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 y5X 1 b1WX 2 d1E 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 3 N matrix of weights encoding the geographical neighbourhoods of each location.
The term WX 2 d represents neighbourhood or spatially 'lagged' effects of the explanatory variables in X 2 on the response variable, y, and the estimated parameters d represent the strengths of these effects for each explanatory variable in X 2 . Third, the spatial error model y5X 1 b1 u was considered; and fourth, the spatial Durbin error model y5X 1 b1 WX 2 d1u (Elhorst, 2010)  The variance inflation factor (VIF; Kutner, Wasserman, & Neter, 2004) was used to test for collinearity, and non-negligible collinearity was detected between SST and irradiance (but not between SST synchrony and irradiance synchrony). Following the guidelines of Dormann, Elith, and Bacher (2013), principal component analysis of these two variables was performed, and the variables were replaced by the first and second principal components, labelled 'PC1' and 'PC2'. After this replacement, the VIF indicated no non-negligible collinearity among predictor variables (Supporting Information Tables S1-S3), again following Dormann et al. (2013) in assessing collinearity. An increase in PC1 corresponded to increases in both SST and irradiance, whereas an increase in PC2 corresponded to an increase in irradiance and a decrease in SST.
To limit the number of models fitted, the total number of variables allowed in any one model was restricted to five or fewer. Without this simplification, the total number of models for the all-frequencies analysis alone would have been 118,096 for each ocean basin, a computationally prohibitive number also considered excessive for model selection (Burnham & Anderson, 2002). We sought main determinants of synchrony, so we considered only simple models. To test whether results were sensitive to this choice, we repeated the model selection exercise for the Pacific Ocean, for all frequencies, allowing models with DEFRIEZ AND REUMAN | 5 six or fewer variables. For OLS models with lagged variables and 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 effects). Neighbourhood effects without local effects 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. The BIC was used instead of the Akaike information criterion (AIC) because it is known that AIC favours more complex models (Burnham & Anderson, 2002) and we sought main determinants of synchrony. The importance of a given variable as a predictor of chlorophyll 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, 2002). Residual plots were produced for top models to check model assumptions.
BIC approaches can identify which of several models is best supported (Burnham & Anderson, 2002) but do not 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 chlorophyll synchrony explained by the model. Second, we compared, via BIC and ANOVA, our best (lowest BIC) model from each model comparison exercise with the null model y5X 1 b1u, for X 1 a column of ones 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).

| RE SULTS
Plotting all-frequency chlorophyll synchrony (500 km radius; Figure 2a) provided answers to the first question of the Introduction (Q1). There is a pronounced geography of synchrony. Geographical variation in strength of synchrony is strikingly strong, spanning almost the entire range of possible values from 0 to 1. Highest chlorophyll synchrony values tended to occur in the centres of oceans, and preferentially at lower latitudes. In coastal areas, synchrony was typically much reduced, for instance, off the coasts of China and Japan and around the islands of Oceania. The largest strongly synchronized area was in the middle of the Pacific Ocean. Patterns of synchrony were broadly similar for 500-1000 and 1000-1500 km (Supporting Information Figure S1). Figures   S2-S4) showed declines typical of such plots and did not reveal detail as did maps of the geography of synchrony.

Plots of synchrony versus distance (Supporting Information
Plotting observed high-and low-frequency synchrony (Figure 2c,d, respectively) answered Q2; the geography of synchrony was frequency specific. For instance, synchrony along the equator in the Atlantic and Pacific Oceans occurred principally at high frequencies. There was also a patch of synchrony off the coast of Chile, mostly at high frequencies.
At many places in the central Pacific Ocean away from the equator, synchrony occurred predominantly at low frequencies.
Model-averaged predictions of chlorophyll synchrony (all frequencies) from the 95% confidence set of models from individual ocean basin analyses (Figure 2b) were reasonable first approximations.  porting Information Tables S4-S6) and were always far superior to the null model y5X 1 b1u, for X 1 a column of ones, according to BIC and ANOVA (Appendix S8). Models explained considerably less variation in low-or high-frequency chlorophyll synchrony than in all-frequencies synchrony (Supporting Information Tables S4-S6), as expected given the greater variability intrinsic to estimates of frequency-specific quantities. For low-frequency synchrony in the Atlantic Ocean, the worst    Information Table S5).
DEFRIEZ AND REUMAN | 7 case, models explained none of the geographical variation in chlorophyll synchrony (Supporting Information Table S5).
Statistics showing the relative importance of variables in determining geographical patterns of synchrony (Tables 1-3) answered Q3. Log chlorophyll is an important statistical determinant of chlorophyll synchrony across all ocean basins and frequencies and has a negative association; regions with more chlorophyll tended to have reduced synchrony. SST synchrony is an important predictor of chlorophyll syn-chrony (all frequencies) in the Atlantic and Pacific Oceans, with positive coefficients in both cases, confirming hypothesis H1; more SST synchrony is associated with more chlorophyll synchrony, indicating apparent Moran effects. However, we emphasize again that this association need not be directly causal, in that SST might serve as a surrogate for nutrient effects or might be an indicator of complex processes (Discussion). SST synchrony is also positively associated with chlorophyll synchrony in the Indian Ocean, although summed BIC weights indicated that it is less important as a predictor there than in the Atlantic and Pacific Oceans (Table 3). For the Pacific Ocean, all-frequencies analysis, allowing six-variable models showed essentially the same results as  Table S5).

Frequency specificity of Moran effects was not evident in the Indian
Ocean (Table 3). Moran effects of irradiance synchrony were not apparent in any ocean basin for any frequency band (but see Discussion).
In addition to frequency specificity of apparent SST Moran effects, there were also other differences by frequency band in what factors were important determinants of synchrony, confirming H2. For instance, in the Pacific Ocean PC1 was an important predictor of highbut not low-frequency chlorophyll synchrony, as was PC2. Absolute u current was an important predictor of low-but not high-frequency synchrony.
Our results also reveal a degree of heterogeneity across ocean basins in the statistical determinants of synchrony, although with certain determinants operating universally or commonly, answering Q4.  In the Atlantic and Pacific Oceans, we tried ancillary models that included absolute latitude as a predictor; these were no better by BIC than the models in Supporting Information Tables S4 and S5. In the Indian Ocean, models with absolute latitude were substantially better than those of Supporting Information

| DISCUS SION
Here, we used the excellent spatial coverage of satellite chlorophyll data to map the geography of phytoplankton synchrony and to find main statistical determinants of geographical variation in synchrony.
Perhaps our most striking result is the simple demonstration of the existence of a geography of synchrony and the strength of this phenomenon; spatial variation in synchrony encompasses essentially the entire biologically reasonable range of values of the synchrony statistic (0 to 1). Regional weather cannot be manipulated, but there is substantial variation in weather and its synchrony. Hence, large-scale analyses allowed us to examine how chlorophyll synchrony changes as physical conditions vary with location, thereby using inference to make up, in part, for the impossibility of experimental manipulation (Abbott, 2007).
We found that across the Pacific and Atlantic Oceans, synchrony in SST was an important statistical determinant of chlorophyll synchrony.
In the Pacific Ocean, this effect was frequency specific. Other drivers of synchrony also depended on frequency band.  & Dutkiewicz, 2002). Thus in the subtropical North Atlantic Ocean, greater spring mixing leads to enhanced chlorophyll, but in the subpolar North Atlantic Ocean, it can lead to reduced chlorophyll (Dutkiewicz et al., 2001). But Follows and Dutkiewicz (2002) provide evidence that interannual fluctuations in the strength of mixing eclipse regional variation in subtropical but not subpolar areas, suggesting enhanced phytoplankton synchrony in subtropical relative to subpolar regions. This is exactly what we observed (Figure 2a). Signorini, Franz, and McClain (2015) explore similar ideas. SST is known to be a reasonable indictor of mixed layer depth in the gyres (Signorini et al., 2015), and SST is very spatially synchronous, possibly further explaining both the accentuated phytoplankton synchrony we observed in the North Atlantic gyre and the apparent SST Moran effects detected by our statistics for the Atlantic Ocean. This possible mechanism for some of the synchrony patterns we observed implicates light availability, whereas we found that synchrony maps for surface irradiance were unassociated with synchrony maps for phytoplankton. This is not a contradiction, because h c depends on water clarity as well as surface irradiance.

| Mechanisms of synchrony
In contrast to the Atlantic Ocean, our results showed no important association between SST synchrony and phytoplankton synchrony in the Indian Ocean, and this appears consistent with known differences between the Atlantic and Indian oceans in drivers of chlorophyll dynamics. Plankton blooms in the Indian Ocean are associated with increased nutrient supply and mixed-layer thickening related not to SST but to monsoon winds (section 7.2.2 of Williams & Follows, 2011).
Winds are attributable to differential heating of land and sea, and resulting blooms exhibit a latitudinal gradient ( fig. 7.12 of Williams & Follows, 2011), which may explain why our ancillary models revealed a negative association between absolute latitude and phytoplankton synchrony.
Local population dynamics can interact with dispersal to generate spatial synchrony over distances an order of magnitude greater than that of dispersal (Gouhier, Guichard, & Menge, 2010). One way to distinguish between the synchronizing effects of dispersal and environmental factors is to study systems in which one mechanism is not possible (for example, Myers et al., 2015), although these systems are uncommon. As dispersal of phytoplankton at the scales of the present study is reliant on physical processes, if dispersal were an important causal factor in chlorophyll synchrony we might expect the wind and/ or current variables to have been important determinants of synchrony in our analysis. But this was generally not the case. This logic does not completely preclude dispersal as a factor, but it does suggest it has secondary importance to Moran effects manifesting via complex mechanisms, as above.

| Other patterns in synchrony
Another statistical determinant of synchrony that was important across all three oceans was log chlorophyll concentration, which was negatively associated with chlorophyll synchrony. Similar results have also been found at sub-annual scales (Uz & Yoder, 2004). This was probably because coastal waters tend to be nutrient rich, promoting phytoplankton growth, and unstable, reducing synchrony, for the same reasons, DEFRIEZ AND REUMAN | 9 namely upwelling, waves and other oceanic activities (Holt, Harle, & Proctor, 2009). Open oceans have the highest levels of water-column stability and synchrony, and typically also the lowest levels of chlorophyll. Differences between coastal and oceanic synchrony have also been observed in zooplankton (Batchelder et al., 2012).
During the period of this study there were three El Niño events (2004-05, 2006-07 and 2009-10). These events would fall into the high-frequency band of our analyses and might explain some of the differences seen between high and low frequencies. El Niño may be part of the cause of the high-frequency synchrony seen along the equator in the Pacific and Atlantic Oceans (Figure 2c). El Niño is associated with increased SST and therefore increased stratification of surface waters over large areas. This should cause changes in chlorophyll over similarly large areas, thereby producing high-frequency chlorophyll synchrony.

| Future work
It seems likely that in most natural systems, synchrony of populations cannot be ascribed to one particular driver, but rather is attributable to a combination of factors (Kendall et al., 2000). For instance, the mecha- here but not explained, or might seek higher-resolution results, or might include factors for which data were not available for our study, such as direct measurements of nutrients or consumers of phytoplankton. Ultimately, we believe that our work helps to open the door on the intriguing approach of directly mapping the geography of synchrony but has taken only the first step through that door. We hope that our study inspires further analyses of the geography of synchrony in phytoplankton. Additional possible avenues for future work are in Appendix S6.