Spatial synchrony cascades across ecosystem boundaries and up food webs via resource subsidies

The Santa Barbara Channel is a highly productive coastal region with nearshore rocky reef habitat dominated the foundation species giant kelp (Macrocystis pyrifera). This region has mixed coastline types but is largely rocky intertidal and sandy beach. Because of the strong ecological connectivity between kelp forests and sandy beaches (17, 18, 23), sandy beach monitoring is conducted at five study sites within the Santa Barbara Channel as part of the Santa Barbara Coastal Long Term Ecological Research (SBC LTER) program (Fig. 3B). These sites have been surveyed monthly since 2009, with the exception of April to June 2020, and include Arroyo Quemado (34${}^{°}$28’13.3”N, 120${}^{°}$07’09.6”W), Isla Vista (34${}^{°}$24’33.4”N, 119${}^{°}$52’27.0”W), East Campus (34${}^{°}$24’38.6”N, 119${}^{°}$50’31.3”W), Arroyo Burro (34${}^{°}$24’11.0”N, 119${}^{°}$44’38.3”W), and Santa Claus Lane (34${}^{°}$24’31.4”N, 119${}^{°}$33’06.7”W).

Surveys of giant kelp wrack abundance (cover) and dry beach width at the five study sites are conducted monthly at low tide ($\le$ 0.76 m) along three permanent cross-shore transects. Wrack cover is measured as the total length of giant kelp intersecting the measurement transect for a 1 m wide band from the back beach boundary (i.e., cliff base) to the upper swash limit (17). Wrack is classified as fresh or old and is identified by part of the kelp plant (blade, stipe, or holdfast). We focus our analyses on total blade and stipe wrack because blades and stipes comprise the majority of wrack and are consumed at far higher rates than holdfasts. Wrack cover along replicate transects was averaged by site for each monthly survey. Dry beach width is measured along each transect as the distance from the back beach boundary to the 24-h high tide strandline. The dry beach width measurement is representative of the available beach area for wrack deposition and retention. We used monthly giant kelp wrack data and associated dry beach width data for the period 2009 to 2019 (60).

Data on the biomass of live giant kelp are derived from Landsat satellite imagery using a spectral unmixing algorithm and an empirical relationship relating kelp canopy cover to biomass (34, 35). We obtained monthly data at the spatial resolution of 100 m coastline segments and spatially aggregated the data at two spatial scales: “local,” i.e., averaged within a 1.5-km radius of the beach site, and “regional,” i.e., averaged over all locations along the mainland coast from Point Conception to Carpenteria and San Miguel, Santa Rosa, and Santa Cruz islands (approximately 280 km coastline).

Data on daily maximum significant wave heights corresponding to each beach site were obtained from the CDIP MOP v1.1 model (http://cdip.ucsd.edu/MOP_v1.1/) and were aggregated to a monthly time step, matching the beach wrack data, by averaging. In brief, the model combines hourly empirical wave height and direction measurements with swell propagation and hindcast models (61). Detailed methods for computation of the wave intensity data are given in ref. 37. We analyzed first-differenced time series of wave intensity to better correspond to the mechanisms linking wave intensity to kelp wrack. Following a summer peak in kelp biomass, wave intensity tends to increase through the fall and winter. Initially, relatively small increases in wave intensity can dislodge kelp leading to an increase in wrack; however, once waves have begun to dislodge kelp increasingly strong waves are required to further dislodge kelp and increase wrack.

Data on the abundances of shorebirds at the same sites are collected concurrently with wrack sampling (60). Observers identify (generally, to species) and count the number of birds on 1 km along-shore transects. We considered shorebirds to include American Avocet, Baird’s Sandpiper, Black-bellied Plover, Black-necked Stilt, Black Turnstone, unidentified Dowitchers, Dunlin, Greater Yellowlegs, Killdeer, Least Sandpiper, Lesser Yellowlegs, Long-billed Curlew, Long-billed Dowitcher, Marbled Godwit, Pectoral Sandpiper, Red-necked Phalarope, Red Knot, Ruddy Turnstone, Sanderling, Semipalmated Plover, Short-billed Dowitcher, Snowy Plover, Spotted Sandpiper, Surfbird, Wandering Tattler, Western Sandpiper, Whimbrel, Willet, and Wilson’s Plover. We totaled the abundances of all shorebird species to obtain one aggregate shorebird abundance for each beach and survey event.

Data on monthly average daily air temperatures were obtained at the 5 beach sampling locations from the PRISM dataset (62). Air temperature effectively accounts for seasonal variation in shorebird abundances associated with their seasonal migrations; additionally, local variations in bird migration timing are associated with air temperature (38, 39).

We quantified spatial synchrony in kelp wrack using wavelet mean fields, which assess spatial synchrony over all sites as a function of time and timescale (27). Statistical significance of spatial synchrony was determined using wavelet phasor mean fields, a closely related technique for which a null model is available, enabling significance testing (27). Statistical significance was assessed relative to 1,000 sets of random phasors, representing a null hypothesis of no synchrony except that arising by chance. Detailed descriptions of these mean field methods, and other analysis methods, and pedagogical figures illustrating their use, are provided in SI Appendix.

To examine the synchronizing roles of kelp biomass, waves, and beach width, we used spatial wavelet coherence (27) and wavelet linear models (33). Spatial wavelet coherence tests timescale-specific bivariate relationships between two spatiotemporal variables, giving, for each timescale of interest, a magnitude corresponding to the strength of relationship between the two variables and a phase indicating the temporal offset between peaks in the oscillations of the two variables occurring on the timescale (27). When one variable is an external environmental driver and the other is a biological variable, a statistically significant spatial wavelet coherence relationship indicates that the environmental driver, or possibly some third variable related to both the environmental driver and biological response, is a mechanism of synchrony in the biological response (27). We tested for bivariate spatial wavelet coherences between kelp wrack and local kelp biomass, channel-wide kelp biomass, waves, and beach width (4 drivers $\times$ 3 timescale bands = 12 tests). To infer how synchronizing mechanisms change across timescales, we tested for coherence on 2 to 8, 8 to 16, and 16 to 60 mo timescale bands. Significance testing was performed using the “fast” method introduced by Sheppard et al. (50), with 10,000 surrogate (i.e., appropriately randomized) datasets used to test for significance. Because multiple factors drive kelp wrack subsidies, we followed the methods of Sheppard et al. (33) in constructing multivariate wavelet linear models using the combination of all candidate environmental drivers shown to be coherent ($P<0.1$) with kelp wrack, individually for each timescale band. If both local and channel-wide kelp biomass were coherent with kelp wrack, the one with the lower $P$-value was selected for the multivariate model. We then applied the wavelet Moran theorem and synchrony attribution theorem of ref. 33 to estimate the fraction of synchrony in kelp wrack explained by these models at a given timescale band.

To test whether synchrony in kelp wrack cascades across trophic levels to produce spatial synchrony in the abundance of shorebirds, we tested for spatial wavelet coherence between shorebird abundance and wrack cover (3 timescale bands = 3 tests). We subsequently applied the wavelet Moran and synchrony attribution theorems to wavelet linear models that included air temperature and wrack cover as predictors. Air temperature was included as a predictor in the wavelet linear model to account for the seasonal cycle in bird abundances arising from their migrations, which improved model diagnostics testing assumptions underpinning quantification of synchrony explained and contributions of driver variables to synchrony. All computations were done in R, using the wsyn package, available on the Comprehensive R Archive Network.

Given the potential for multiple testing to produce false positives, we report (above) the numbers of significance tests performed, in order to aid interpretation of strength of evidence for our conclusions. For instance, of the 12 tests in Table 1, a type I error rate of $0.05$ would tend to yield $0$ or $1$ false positives if the null hypothesis of no relationship between variables were true; instead, we observed $7$ significant results using a threshold of $0.05$, a substantially greater number. Note that our $P$-values are additionally conservative because they assume that the power spectrum (timescale structure) of response and driver variables arise independently, though in a true driver–response relationship, the driver would, under reasonable linearity assumptions, make the power spectrum of the response variable more similar to its own (51).

J.A.W., T.W.B., L.W.S., K.C.C., D.C.R., and M.C.N.C. designed research; J.A.W., K.A.E., J.E.D., D.M.H., T.W.B., and V.A.K. performed research; J.A.W., L.W.S., and D.C.R. contributed new reagents/analytic tools; J.A.W. analyzed data; and J.A.W. and K.A.E. wrote the paper.

The authors declare no competing interest.