Spatially varying plankton synchrony patterns at seasonal and interannual scales in a well-connected shelf sea

Abstract

Spatial population synchrony, defined as spatial covariation in population density fluctuations, exists across different temporal and spatial scales. Determining the degree of spatial synchrony is useful for inferring environmental drivers of population variability in the wake of climate change. In this study, we applied novel statistical methods to detect spatial synchrony patterns of Calanus finmarchicus on the Northeast U.S. Shelf at multiple spatiotemporal scales using unevenly distributed data. Our results reveal that C. finmarchicus subpopulations connected by advection are not necessarily in synchrony, indicating that the degree of synchrony is likely influenced by heterogeneity of local habitats. In addition, regionally synchronous environmental conditions (e.g., sea surface temperature) may not play as significant a role in influencing subregional population dynamics as was previously hypothesized. Overlooking the spatial heterogeneity of synchronous patterns at different time scales could lead to erroneous inferences of potential environmental drivers responsible for C. finmarchicus variability.

---

Scientific Significance Statement

Spatial population synchrony is the property of correlating population fluctuations through time. Detecting patterns of spatial synchrony is often challenging in dynamic ocean systems, but it is useful for inferring environmental drivers of population variability in the wake of climate change. This study presents a novel approach to assess the spatial synchrony of Calanus finmarchicus on the Northeast U.S. Shelf, revealing previously overlooked synchrony patterns. Our findings suggest that spatially proximate subpopulations are not necessarily homogeneous, and the spatial scale of synchrony varies with the temporal scale of interest. We further find that subpopulations connected via advection and experiencing synchronous temperature conditions are not in synchrony, indicating that heterogeneity of local habitats likely plays a more important role in driving population variability.

---

Climate change is rapidly altering marine ecosystems, resulting in spatial and temporal variability of vital organisms (Grieve et al. 2017). One ecosystem that is particularly vulnerable to the impacts of climate change is the Northeast U.S. Shelf (NES), which is experiencing faster-than-average warming, altered circulation patterns, reduced salinity, and shifts in the spatiotemporal distribution of resident organisms (McHenry et al. 2019). Of prominent interest on the NES is the lipid-rich copepod Calanus finmarchicus, as this species is an important climate indicator and plays a key role in linking lower and higher trophic levels (Runge et al. 2015; Ji et al. 2021). However, identifying and modeling the environmental drivers responsible for C. finmarchicus population variability remains a significant challenge (Holt et al. 2014; Freer et al. 2022).

One promising approach to better understand the processes influencing copepod population dynamics is by analyzing coincident abundance fluctuations of spatially distinct populations, known as spatial synchrony (Liebhold et al. 2004). For decades, synchrony has been widely studied in terrestrial (Liebhold et al. 1996; Post and Forchhammer 2004; Lindström et al. 2012) and freshwater systems (Myers et al. 1997; Lodi et al. 2014), and in more recent years, detecting synchrony in marine planktonic systems has been useful for inferring environmental drivers of population variability at various spatiotemporal scales (Defriez et al. 2016; Sheppard et al. 2019). Multiple mechanisms can influence patterns of synchrony, including dispersal, trophic interactions, and climatic forcing. Understanding the scale of synchrony allows researchers to hypothesize the environmental factors that impact population dynamics. For instance, seasonal synchrony can be attributed to smaller-scale, homogenous forces, while interannual synchrony can indicate that larger-scale forces may dictate long-term population variability (Liebhold et al. 2004). In addition, changes in synchrony patterns can lead to ecological problems such as abnormal phenological shifts, trophic mismatch, decreased community stability, and even extinction (Post and Forchhammer 2004; Abbott 2011; Doiron et al. 2015). Therefore, by studying the spatial patterns and scale of synchrony, researchers can surmise the relative importance of various environmental factors that influence population dynamics, thus helping project future population trajectories and assess potential ecological consequences.

Various statistical techniques have been developed to detect patterns of spatial synchrony. Previous studies have utilized spline correlogram and wavelet analysis to investigate the spatial synchrony of plankton species and infer what factors may be contributing to observed population fluctuations (Defriez et al. 2016; Sheppard et al. 2019). These techniques require input data to be uniformly distributed across time and space, which is not easily accomplished when sampling oceanic systems due to accessibility restraints and resource limitations. Although it is possible to interpolate unevenly distributed data onto a uniform grid, this makes an inherent assumption of spatial and temporal coherence.

Previous studies have suggested that regional warming or North Atlantic Oscillation may play a key role in synchronizing allopatric subpopulations of C. finmarchicus in the North Atlantic (Perry et al. 2004). However, before we can ascertain whether these or other drivers are influencing C. finmarchicus variability, it is necessary to assess the underlying patterns and scales of spatial synchrony in dynamic shelf regions like the NES. Here, we introduce a new approach to detect spatial synchrony in unevenly distributed spatiotemporal survey data. Our specific objectives are to (1) investigate the small-scale seasonal and interannual spatial synchrony patterns of C. finmarchicus, (2) compare the seasonal and interannual synchrony of subpopulations in two or more spatially distinct locations that are connected via advection, and (3) hypothesize potential drivers that may be contributing to these observed patterns for future research.

Methods Data

Plankton and sea surface temperature (SST) data come from the Marine Monitoring Assessment and Prediction (MARMAP) plankton survey (1977–1987) and the Ecosystem Monitoring (EcoMon) program (1988 to present) collected by the U.S. National Oceanic and Atmospheric Administration (NOAA) Northeast Fisheries Science Center (NMFS/NEFSC). Plankton samples were obtained approximately every 2 months at fixed or randomly selected stations using a 61-cm bongo net with a 333 μm mesh size towed at a maximum depth of 200 m. Refer to Sherman (1980), Meise and O'Reilly (1996), and Kane (2007) for more survey protocols. The NES region was separated into 46 distinct strata based on bathymetry (Fig. 1a), and we assigned each sample of the MARMAP/EcoMon data set to a stratum based on sampling location. Although the MARMAP/EcoMon data set is comprehensive, it is irregularly distributed across space and time (Fig. 1b).

Fig. 1. The spatial area of the Northeast U.S. Shelf separated into 46 distinct strata (a). Data for this region are comprehensive, but samples are unevenly distributed across space and time (b). Stratum coordinates were provided by Harvey Walsh from the NOAA NEFSC.

Statistical analysis

The statistical analysis had two goals: (1) to assess seasonal and interannual spatial synchrony of C. finmarchicus within each stratum and (2) to assess seasonal and interannual synchrony between two or more strata. In both cases, the analysis was based on fitting generalized additive models (GAMs), which were developed by Hastie and Tibshirani (1986) as an alternative to multiple linear regression (Hastie and Tibshirani 1986). The advantage of GAMs is that they do not require the specification of the functional forms of the dependence of the response variable on the regressors, and they have proven to be effective in analyzing complex ecological systems (Guisan et al. 2002; Drexler and Ainsworth 2013; Pedersen et al. 2019). Before fitting to the GAM, we translated the approximate bi-monthly sampling dates of the MARMAP/EcoMon data set to day-of-year. Throughout our statistical method, we aggregated all data sampled from within each stratum and utilized a standard identity link function to simulate abundance. All GAM analyses were performed using the “mgcv” R package, and we will retain the notation of that package (Wood 2017).

The general, unrestricted GAM considered here is:[Image Omitted. See PDF]where $R_{\mathit{ij}}$ is density of C. finmarchicus in sample $i$ ($i=1,2,\dots ,n_j$) from stratum $j$ ($j=1,2,\dots ,46$); ${\beta }_{\mathit{ij}}$ represents the intercept; ${\mathrm{u}}_{\mathit{ij}}$ and ${\mathrm{v}}_{\mathit{ij}}$ are the latitude and longitude of this sample; ${\mathrm{DoY}}_{\mathit{ij}}$ and ${\text{Year}}_{\mathit{ij}}$ are the day-of-year and year, respectively, of this sample; and ${\epsilon }_{\mathit{ij}}$ is a normal error with mean 0 and variance ${\sigma }^2$. Here, $s_{1j}$, $s_{2j}$, and $s_{3j}$ are smooth functions representing the direct effects of location, day-of-year, and year on $Y_{\mathit{ij}}$; and ${\mathit{ti}}_{1j}$, ${\mathit{ti}}_{2j}$, and ${\mathit{ti}}_{3j}$ are tensor interaction smooths representing interactions between the regressors. For example, ${\mathit{ti}}_{1j}$ represents a smooth interannual change in the seasonality of log abundance. Following standard practice, we defined density as $\left[\left(\text{abundance}/{\mathrm{m}}^3\right)+1\right]$ and applied a log transformation both to ensure positivity and to stabilize variance, which, as is typical with such data, increases with the mean.

Single-stratum analysis

We assessed seasonal asynchrony within stratum $j$ by testing the null hypothesis ${\mathrm{H}}_0:{\mathit{ti}}_{2j}=0$ that the seasonal cycle does not depend on location within the stratum. Similarly, we assessed interannual asynchrony by testing the null hypothesis $H_0:{\mathit{ti}}_{3j}=0$, that interannual variability does not depend on location within the stratum. We tested these null hypotheses against the alternative hypothesis $H_1$ that seasonal/interannual variability does depend on location within a stratum, as in Eq. 1.

We tested $H_0$ using the quasi-likelihood ratio ($\mathrm{qLR}$) statistic:[Image Omitted. See PDF]where ${\mathrm{RSS}}_0$ and ${\mathrm{RSS}}_1$ are the residual sums of squares for the null and alternative models, respectively. This is not a true likelihood ratio statistic because the model is not fit by maximum likelihood.

The following bootstrap procedure was used to assess the significance of the $\mathrm{qLR}$ statistic. First, the $n$ residuals from the fitted null model (${\widehat{Y}}_{\mathit{ij}}^0$) were sampled with replacement. These bootstrapped residuals are denoted by ${\epsilon }_{\mathit{ij}}^{*}$, where $i=1,2,\dots ,n$. With these residuals, we formed a new set of values for $Y_{\mathit{ij}}^{*}$.[Image Omitted. See PDF]

Second, we refit the unrestricted and null models to this bootstrap sample and recorded the value of the $\mathrm{qLR}$ statistic. Finally, we repeated the bootstrap procedure 100 times and approximated the observed significance level (or p-value) by the proportion of times the value of the statistic exceeded the observed value. Supporting Information S1 provides further details and the results of a power study of this single-stratum analysis (SSA) method.

Multiple-strata analysis

Seasonal asynchrony between two strata $j$ and $k$ for which synchrony could not be rejected by the SSA was assessed by testing the null hypothesis $H_0:s_{2j}=s_{2k}$; ${\mathit{ti}}_{1j}={\mathit{ti}}_{1k}$ of a common seasonal cycle in every year against the general alternative $H_1$. Similarly, interannual asynchrony was assessed by testing the null hypothesis $H_0:s_{3j}=s_{3k}$; ${\mathit{ti}}_{1j}={\mathit{ti}}_{1k}$ of a common interannual trend in every day of the year against the general alternative $H_1$. As before, these tests were based on $\mathrm{qLR}$, with significance assessed by the residual bootstrap. Refer to Supporting Information S2 for further details and a power study of this multiple-strata analysis (MSA) method. Since C. finmarchicus is primarily a cold-water species, we focus on using the MSA to analyze the synchronicity of subpopulations in the Georges Bank (GBk) and Gulf of Maine (GoM) regions (Grieve et al. 2017). In addition, in order to investigate a potential environmental driver of the observed C. finmarchicus spatial synchrony patterns, we utilized the MSA to analyze interannual SST synchronicity by changing the GAM response variable from $Y_{\mathit{ij}}$ to ${\mathrm{SST}}_{\mathit{ij}}$ (units of °C) in Eq. 1.

Novelty

Our approach differs from other methods developed for sparse and unevenly distributed data (Thorson et al. 2018; Siple et al. 2020), as we do not apply synchrony metrics on interpolated time-series fit from autoregressive models. Instead, the SSA tests for the significance of the interaction between temporal variability at locations in close spatial proximity, and the MSA assesses how seasonal/interannual variability differs across distinct locations. For both methods, our bootstrap approach to test for significance gives a simple binary result of synchrony or asynchrony (we use the terms “synchronized/synchronous” for the case in which the null hypothesis could not be rejected). In addition, including the ${\mathit{ti}}_{1j}$ term in the MSA null hypothesis allows for the detection of asynchrony between strata that appear to have synchronous average seasonal or interannual curves, but have asynchronous seasonal curves across different years. Conversely, the power of the MSA is lower when testing two strata with asynchronous average seasonal or interannual curves, but they each have a seasonal cycle that systematically advances across years (see Supporting Information S2). This allows for analysis of how seasonality changes across years and spatial locations, which is not commonly considered in synchrony studies using unevenly distributed data (Thorson et al. 2018; Thorson 2019; Siple et al. 2020).

Results SSA results

Through our SSA, we find that there are different patterns of seasonal and interannual spatial asynchrony on the NES (Fig. 2). Most of the seasonal synchrony occurs within the GoM, with some seasonal synchrony in the coastal and offshore strata of the Mid-Atlantic Bight (MAB) and Southern New England (SNE) regions (Fig. 2a). Many mid-shelf strata corresponding to depths of approximately 20 to 60 m are seasonally asynchronous. Interannually, many strata on the NES are synchronized (Fig. 2b), with approximately 85% of the 46 strata exhibiting interannual spatial synchrony. Almost half of the strata that are not synchronized appear on GBk (i.e., Strata 29, 30, and 32).

Fig. 2. Seasonal (a) and interannual (b) spatial synchrony patterns from the SSA. Shaded strata indicates spatial synchrony of data collected within that stratum.

MSA results

Comparing seasonal variability across Wilkinson Basin (WB) and Jordan Basin (JB), the average seasonal curves ($s\left(\mathrm{DoY}\right)$ partial effect curve), along with the sum of the average seasonal curve and the $\mathit{ti}\left(\mathrm{DoY},\text{Year}\right)$ interaction term, are displayed in Fig. 3b. The overall average trends of these seasonal curves coincide for most of the year, so through our MSA, we failed to reject the null hypothesis, indicating that these strata are synchronous. Conversely, there is a statistically significant difference in the average seasonal effect curves for C. finmarchicus subpopulations in WB and Northwest Georges Bank (NWGBk), so for these strata, we were able to reject the null hypothesis to assume these strata are asynchronous (Fig. 3c).

Fig. 3. Strata for the GoM and GBk (a), seasonal MSA results for WB and JB (b), and seasonal MSA results comparing WB and NWGBk (c). The thick, dark lines represent the average seasonal GAM-effect (sDoY), while the thin, light lines represent the sum of the average seasonal GAM-effect and the interaction between DoY and year (sDoY+tiDoYYear) to demonstrate seasonal variability across different years.

Using all available data, the three basins of the GoM (WB, JB, and Georges Basin [GBn]) are interannually asynchronous. Figure 4b shows that the abundance peaks in JB and GBn trail the WB peaks on the order of a few years, particularly in the first 20 years of the data. Furthermore, the average interannual GAM effect curves ($s\left(\text{Year}\right)$) using data from the spring months (April, May, and June) appear to be synchronous, but due to variability of the $\mathit{ti}\left(\mathrm{DoY},\text{Year}\right)$ interaction term in JB, these strata are still asynchronous. In addition, testing just WB and GBn, these strata are interannually asynchronous using all data (Fig. 4c), as the interannual abundance fluctuations exhibit compensatory dynamics for several decades (e.g., 1990–2000 and 2010–2019). Yet, the fluctuations become synchronous when only using spring data (Fig. 4f), so the synchronicity of WB and GBn is season dependent. Moreover, comparing the interannual cycles of abundance using all data (Fig. 4d) and just spring data (Fig. 4g) for WB and NWGBk gives the result of asynchrony in both cases, which demonstrates that areas close in spatial proximity and connected through advection are not necessarily spatially synchronized in the dynamic NES. However, SST fluctuations in the basins of the GoM (Fig. 5b) as well as WB and NWGBk (Fig. 5c) are synchronous.

Fig. 4. Strata for the GoM and GBk (a) along with interannual MSA results comparing combinations of the basins of the GoM and NWGBk using all data (b–d) and using just spring data from April, May, and June (e–g). The thick, dark lines represent the average interannual GAM-effect (sYear), while the thin, light lines represent the sum of the average interannual GAM-effect and the interaction between DoY and year (sYear+tiDoYYear) to demonstrate interannual variability across different seasons.

Fig. 5. Strata for the GoM and GBk (a) and interannual MSA results for SST in the basins of the GoM (b) and for WB and NWGBk (c). The thick, dark lines represent the average interannual GAM-effect (sYear), while the thin, light lines represent the sum of the average interannual GAM-effect and the interaction between DoY and year (sYear+tiDoYYear).

Discussion Small-scale spatial heterogeneity on the NES

Even though the NES is a well-connected shelf sea, C. finmarchicus subpopulations close in spatial proximity are not necessarily homogeneous due to the dynamic nature of small-scale environmental conditions throughout the shelf. Our SSA results clearly corroborate that subpopulations in many mid-shelf strata in MAB, SNE, and GBk cannot be represented by a single seasonal or interannual curve. The spatial heterogeneity of population variability is difficult to detect with statistical confidence when data points are sparse and irregular. Hence, asynchronous spatial patterns could be obscured by assuming coherence at certain spatiotemporal scales.

Detecting synchrony patterns at multiple temporal scales allows us to understand and infer potential drivers of population variability. For instance, the prevailing pattern of seasonal asynchrony within mid-shelf strata (i.e., Strata 5, 8, 11, 15, 16, 19 in Fig. 2a) suggests that subpopulations on the coastal sides of the MAB and SNE could be influenced by different environmental drivers than subpopulations on the shelf-break side. It is thus unrealistic to link seasonal subpopulation variability at different locations with a single, shelf-wide, environmental driver as some previous studies have suggested (Conversi et al. 2001; Licandro et al. 2001; Greene et al. 2003). Instead, the seasonal cycle could be influenced by multiple drivers operating at variable spatiotemporal scales, including spatially heterogeneous seasonality of hydrography and productivity cycles, strong cross-shelf gradients of bottom-up and top-down forcing, and a southward decline of advective inflow from upstream source populations (Manning 1991; Pershing et al. 2010).

At the interannual scale, subpopulations are synchronized within most strata throughout the NES (Fig. 2B). This indicates that similar environmental drivers may be responsible for interannual C. finmarchicus variability at a spatial scale equal to or larger than the size of an individual stratum. However, subpopulations in a small portion of the NES (e.g., Strata 25, 29, 30, 32, and 33 near the GBk region, Fig. 2b) are not synchronized. These strata are located in areas between the shallow bank and the deep GoM, where both drivers and population responses may vary at a smaller spatial scale than the rest of the NES region. A small shift of the population supply pathway from the GoM to the GBk coupled with different in situ population source/sink dynamics, could lead to asynchronous population variability within these strata. Using stratum area size as a benchmark, the spatial scale of synchrony at the interannual temporal scale appears to be larger than that at the seasonal scale, suggesting that corresponding environmental drivers are also interannually synchronized at or beyond the individual stratum scale.

Time-scale dependency and advection

At a larger spatial scale, we use the MSA to assess the scale-dependent variability of subpopulations that are connected via advection. A specific hypothesis to test is whether population variability is more likely influenced by environmental factors when the relevant timescales of the two match. This idea has been examined through the linear tracking window hypothesis proposed by Hsieh and Ohman (2006) and a related concept has been applied to synchrony through the study of ecological “detuning” (Hsieh and Ohman 2006; Haynes et al. 2019). In the GoM and GBk, the relevant timescales of advection and population growth are both on the order of several months, so if advection was the primary driver influencing population dynamics among the basins of the GoM and from WB to NWGBk, we would expect to see lagged seasonal asynchrony with the upstream source subpopulation peaking a few months before the downstream subpopulation and overall synchronous interannual fluctuations (Lynch et al. 1998; Luo et al. 2021). Based on the physical circulation of the GoM, modeling studies (i.e., Lynch et al. 1998; Miller et al. 1998; Li et al. 2006) have shown that individuals from JB can get advected south into GBn and southwest into WB. However, our results at all time scales of analysis reflect that advection is likely not the primary driver of population abundance. We find that subpopulations in WB and JB are seasonally synchronous, and based on a power analysis of our method, the seasonal MSA is sensitive enough to detect asynchrony at a time lag of 20 d (with a power greater than 90%), which is shorter than the seasonal time lag associated with advection (Supporting Information S2). Therefore, the MSA is robust to detect lagged seasonal asynchrony at the relevant timescales of advection and population growth. In addition, since all three basins are interannually asynchronous, this, along with our result of seasonal synchrony, suggests that other environmental factors may be playing a more significant role in influencing subpopulation abundances in this region. Furthermore, adjusting the time period of analysis by using only spring data for all three basins, the average interannual fluctuations appear more synchronous. However, the spring populations in the basins are still asynchronous due to high seasonal variability in JB (Fig. 4e). These fluctuations are not similarly reflected in WB and GBn, because the spring data for these basins give the result of synchrony (Fig. 4f). This finding further substantiates that subpopulations in these basins likely depend more on local rates of growth than on advective transport from JB.

Similarly, even though WB and NWGBk are physically connected, these strata are asynchronous at both the seasonal and the interannual time scales. Surface currents in the GoM flow from WB to GBk, yet it is clear that seasonal subpopulations in NWGBk peak in abundance about a month before subpopulations in WB. This suggests that advective population supply from WB to GBk does not play as important a role as was previously hypothesized (Meise-Munns et al. 1990). Therefore, despite the close spatial proximity and strong connectivity of these strata, the asynchronous fluctuations in abundance between these strata could be attributed to either different seasonal drivers or differing responses to the same seasonal drivers.

SST and the Moran effect

Moran's (1953) theorem states that synchrony/asynchrony in allopatric populations is driven by synchrony/asynchrony in environmental variables (Moran 1953). SST has been suggested as a key driver for C. finmarchicus population variability (Hare and Kane 2012; Grieve et al. 2017). Our SST synchronicity test for the basins of the GoM and WB/NWGBk gives the result of synchrony (Fig. 5), while conversely, C. finmarchicus interannual abundances are not synchronized (Fig. 4b,d). These results imply that interannual variability in the basins likely does not depend entirely on SST, and instead, the degree of C. finmarchicus synchrony in the basins and on GBk is more likely influenced by internal growth through season-specific heterogeneity of local habitats. This finding is consistent with theoretical studies examining the interplay between synchrony, advection/dispersal, and environmental correlation. It is widely accepted that dispersal increases spatial synchrony by decreasing local population variability at long timescales (Kendall et al. 2000; Abbott 2011; Luo et al. 2021). Yet, in order to fully understand the impact of dispersal on population dynamics, it is important to consider dispersal in the context of other physical and biological factors (Yang et al. 2022). For instance, Kendall et al. (2000) and Yang et al. (2022) found that the synchronizing effect of dispersal tends to be weaker when environmental fluctuations are positively correlated. In the GoM and GBk, both dispersion through advection and environmental correlation are acting on subpopulations of C. finmarchicus. However, our finding of interannual spatial asynchrony suggests that neither of these factors are substantially driving subpopulation variability at the interannual scale. Therefore, a combination of other environmental forces are more likely contributing to asynchronous interannual fluctuations, such as differing bottom-up or top-down processes or changes in shelf-slope exchange due to a shift in large-scale circulation patterns (Greene et al. 2003; Ji et al. 2021).

Conclusions

With climate change projected to increase subregional environmental variability, identifying spatial asynchrony patterns is an effective method for analyzing how populations may be experiencing different environmental conditions or exhibiting varying responses to the same environmental conditions (Seidov et al. 2021). Our novel statistical methods allow us to detect synchrony patterns using uneven spatiotemporal observation data. Our findings suggest that spatially proximate subpopulations are not necessarily homogeneous and that subpopulations connected via advection and experiencing synchronous temperature conditions are not in synchrony. This highlights the importance of considering small-scale influences on spatiotemporal variability instead of generalizing over large domains, which may lead to erroneous conclusions.

Acknowledgments

The authors acknowledge funding support from the Department of Defense (DoD) through the National Defense Science & Engineering Graduate (NDSEG) Fellowship Program, the Northeast U.S. Shelf Long-Term Ecological Research (NES-LTER) Project (NSF OCE-1655686), and NOAA Projects (NA21OAR4170379 and NA22OAR4310557). The MARMAP/EcoMon data set was provided by the U.S. National Marine Fisheries Services through D. Richardson and H. Walsh. The authors thank Dr. Casey Youngflesh and Dr. Justin Suca for their useful input in the early stages of the model development, and Dr. Jeffrey Runge and Dr. Greg Britten for insightful discussions about the manuscript content. The authors also acknowledge the editor, associate editor, and two anonymous reviewers for their constructive comments to improve this work.