Ecosystem connectivity and trophic subsidies of sandy beaches

. Ecological connectivity can influence the distributions of diversity and productivity among ecosystems, but relationships among multiple marine ecosystems remain relatively uncharacterized. Sandy beaches are recipient ecosystems that support coastal food webs through deposits of drift macrophytes (wrack), and serve as test cases for exploring within- seascape connectivity. We present results from the first comprehensive survey of geographic and temporal patterns of wrack cover and composition on beaches along the North Central Coast of California and test the role of local donor ecosystems and physical factors in predicting wrack distribution. We surveyed wrack at 17 beaches in August 2010, and monthly at a subset of 10 beaches for 13 months. We estimated explanatory variables of (1) local donor ecosystem cover (kelp forests, rocky intertidal, and bays and estuaries), (2) biomass transport, and (3) beach morphology. Regression analyses were used to evaluate relationships among the cover of six key wrack categories and the explanatory variables above, for two time periods. We found persistent geographic variation in wrack composition and detected significant relationships between wrack cover and cover of local donor ecosystems for five of the six wrack categories ( Nereocystis , Zostera , Postelsia , mixed red algae, and mixed brown algae). Transport mechanisms (wind exposure, swell exposure) or attributes of the recipient ecosystem (beach width, beach slope) explained additional spatial variation for three of the six wrack categories ( Zostera, Phyllospadix , and mixed red algae). Our results support the concept of considering ecological connectivity (particularly the role of donor ecosystems upon which recipient ecosystems rely) in the design and management of protected areas.


IntroductIon
Ecosystems have often been defined and managed as self-contained systems; however, there is growing recognition that most ecosystems are porous to some degree and that some rely on connectivity with other ecosystems to fuel their food webs (Loreau et al. 2003, Kool et al. 2013, Menge et al. 2015. The term "spatial subsidies," coined by Polis et al. (1997), frames this process of allochthonous inputs that enrich seemingly discrete systems, ranging from transport of nutrients or biota within terrestrial and aquatic ecosystems to subsidies across their interface. Notable examples include marine inputs supporting terrestrial food webs on desert islands (Polis and Hurd 1996), salmon as vectors of nutrient transport from highly productive coastal ocean ecosystems to less productive lentic and riparian ecosystems (e.g., Polis et al. 2004, Moore et al. 2011, contribution of riparian zones to stream food webs Dietrich 2002, Sabo andHagen 2012), and the subsidies of marine macrophyte production to fuel secondary production in deep-ocean habitats (Harrold et al. 1998, Vetter and Dayton 1998, britton-Simmons et al. 2012 and sandy beaches , Polis and Hurd 1996. importantly, this process is largely donor-controlled (i.e., determined by variation in productivity and delivery of nutrients, detritus, or prey from the donor ecosystem) and thus reliant on spatially and temporally variable biotic or abiotic transport vectors Polis 2004, Spiller et al. 2010). Additionally, subsidies are often pulsed, creating a complex interaction of subsidy and community responses across time (Sears et al. 2004). Patterns of spatial subsidies may be particularly notable in systems with low in situ productivity that are adjacent to highly productive ecosystems (Persson et al. 1996, Polis et al. 1997, have high ratios of subsidy resource to trophically equivalent ambient resource (Marczak et al. 2007), or have high perimeter:area ratios, creating large and potentially permeable boundaries (Polis and Hurd 1996, Whitman et al. 2004, Marczak et al. 2007.
Sandy beaches are classic examples of ecosystems dependent on trophic subsidies and are underappreciated centers of invertebrate production and biodiversity that provide trophic support for fishes and shorebirds (Schlacher et al. 2008, Defeo et al. 2009), making them ideal systems for studying this dynamic process. beaches have little to no autochthonous primary production and therefore rely on donor inputs, such as ocean phytoplankton and drift macrophytes (e.g., McLachlan and brown 2006). A variety of marine and coastal macrophytes (as well as some carrion and wood) can represent major subsidies to beach ecosystems in the form of beach-cast wrack deposits. Some categories of macrophytes are highly nutritious, such as kelps, which can provide significant trophic subsidies for secondary production of invertebrates (Polis and Hurd 1996, ince et al. 2007, Crawley et al. 2009, baring 2014. Others, such as seagrasses, are less palatable to invertebrates, but can provide some nourishment and important physical habitat structure for fauna on exposed beaches (Jȩdrzejczak 2003, Mews et al. 2006, Heck et al. 2008. in turn, wrack-associated invertebrates are important food sources for a range of coastal species, such as shorebirds (including the western snowy plover, a Federally listed threatened species), seabirds, marine mammals, and fishes (bradley and bradley 1993, Anderson and Polis 1998, Stapp et al. 1999. Additionally, a portion of this wrack subsidy-fueled production can be transferred to terrestrial ecosystems via terrestrial invertebrates (Polis and Hurd 1995, Whitman et al. 2004, Paetzold et al. 2008, Mellbrand et al. 2011) and vertebrates (e.g., lizards, barrett et al. 2005, rodents, Stapp and Polis 2003, birds, Nielsen et al. 2013.
in addition to directly fueling secondary production on beaches, wrack deposition plays a number of key roles in ecological processes on beaches. Wrack deposited in storm conditions provides nutrients, traps windblown sand, and may contain coastal strand and dune plant seeds, setting the stage for plant recruitment and dune formation (Dugan and Hubbard 2010). Wrack accumulations may also act as "metabolic hotspots" of nutrient processing, driving a key ecosystem function of beaches (Coupland et al. 2007, Spiller et al. 2010, and thus, wrack accumulations have great indirect effects as well. Given the importance of wrack in fueling diverse sandy beach ecosystems and processes, it is important to understand the spatial and temporal variation in cover and composition of beach wrack, the biotic and abiotic determinants of that variation, and more generally, the patterns of connectivity between offshore donor ecosystems and the recipient beach ecosystems. The enormous variation in abundance and composition of wrack on beaches can be influenced by both physical and biological drivers at various spatial and temporal scales. There are a large number of potential causal factors, which can be organized by considering their influence on the steps that lead to presence of wrack on beaches: (1) availability of donor biomass, v www.esajournals.org LiEbOWiTz ET AL.
(2) detachment and transport of donor material, and (3) deposition and retention of donor material in the recipient system. Availability of wrack biomass to the beach is proximally controlled by the habitat extent, productivity, and phenology of a macrophyte source. Detachment of macrophytes can be controlled by characteristics of the wrack species (e.g., strength of attachment, profile in flow, phenology, and seasonal senescence) or storm-driven forces leading to high orbital velocities and detachment, and attributes of the substratum to which they are attached (Polis et al. 1991, blanchette 1997, Koehl 1999, Gaylord et al. 2008. Transport and dispersal of the detached macrophytes is related to the buoyancy of the macrophyte and driven by a wide variety of factors, such as winds, waves, and currents (Kirkman andKendrick 1997, Whitman et al. 2004) which vary seasonally, and tides (Orr et al. 2005). Deposition and retention on the beach may be controlled by physical characteristics of the beach itself, such as beach slope, width, length, and substrate type (Orr et al. 2005, Revell et al. 2011, Gómez et al. 2013, or the characteristics of the macrophytes, such as buoyancy, size, form, life stage, and palatability (Koop and Field 1980, Ochieng and Erftemeijer 1999, McLachlan and brown 2006, Duong and Fairweather 2011, Oldham et al. 2014. The coast of California has large stretches of ecologically and economically important sandy beaches, with a varied mosaic of productive nearshore habitats, making it an ideal region to study ecosystem connectivity and subsidies. Additionally, the state recently established a network of 119 marine protected areas (MPAs), designed to function as a spatial network to support conservation and of marine species, especially invertebrates and fishes that use different ecosystems at different life history stages. The network design, intended to support organisms that have open populations at a local scale and exhibit metapopulation dynamics, did not fully account for ecosystem connectivity and metaecosystem dynamics, which involve detrital, energetic, and material flows (Loreau et al. 2003, Massol et al. 2011. Sandy beach food webs may be especially sensitive to changes in management of adjacent donor ecosystems that influence the amount and composition of macrophyte wrack supplied to beach ecosystems. Adaptive, ecosystem-based management of MPAs should be informed by a broader understanding of meta-ecosystem dynamics, the spatial and temporal patterns of connectivity among donor and recipient ecosystems, the processes that drive the patterns of connectivity, and the scale at which these processes occur (e.g., Menge et al. 2015).
Here, we examine ecological questions of regional ecosystem connectivity, which have implications for adaptive management of the network of MPAs along this diverse and highly utilized coast. Specifically, we surveyed spatial and temporal patterns of beach wrack distribution and modeled those patterns in relation to local donor ecosystems and physical factors that can influence detachment, transport, deposition, and retention of macrophytes on sandy beaches. Although patterns of wrack subsidies have been explored in recent years within other regions, this is the first study to couple several wrack categories with multiple donor ecosystems. We ask two main questions: (1) What are the patterns of beach wrack cover and composition across space and time? (2) What are the predictive correlates of the spatial patterns of cover of six different wrack categories from three prominent donor ecosystems?

Study region
We conducted the first large-scale study of sandy beaches along the North Central Coast (NCC) region of California, in conjunction with the establishment of a network of 25 MPAs and marine managed areas in the region in May 2010. The NCC study region extends from Pigeon Point (37°10.55′ N, 122°23.41′ W) to Point Arena (38°57.35′ N, 123°44.50′ W) and is located in one of the four highly productive upwelling regions of the world. The NCC region is characterized by strong prevailing northwesterly winds that peak during spring and summer, a relaxation season in late summer where strong winds are less prevalent, and a winter storm season (Largier et al. 1993, Hickey 1998, García-Reyes and Largier 2012. The study region has large expanses of sandy beaches (51% of the 592 km of shoreline), and a mosaic of productive ecosystems including kelp forests, shallow rocky reefs, rocky intertidal, v www.esajournals.org LiEbOWiTz ET AL. and large estuarine embayments with numerous small rivers.

Approach
We compared spatial and temporal variation in beach wrack cover among NCC beaches and then used stepwise multiple linear regressions to relate these results to estimates of local cover and proximity of donor ecosystems, macrophyte detachment and transport variables, and metrics of beach morphology (variables summarized in Table 1). We conducted separate analyses for geographic patterns of wrack cover for six key categories of wrack, with each wrack category modeled at two temporal categories: (1) August, the month with peak average wrack deposition, and (2) June-December, the period of high average wrack deposition. We capitalized on data sets generated by multiple projects within the North Central Coast Marine Protected Area baseline Program (www.Oceanspaces.org), combined with contextual data from the National Data buoy Center (www.ndbc.noaa.gov) and a wave model (Simulating Waves Nearshore (SWAN); Storlazzi et al. 2005).

Geographic and temporal distribution of wrack composition and cover
To characterize geographic variation in wrack composition and cover among beaches across the NCC study region, we surveyed 17 beaches once in August 2010 (Table 1; Appendix S1). We surveyed a subset of 10 of those beaches monthly from May 2010 to July 2011 (excluding July 2010) to generate a more robust estimate of geographic variation and determine how well the instantaneous (August) estimate of spatial variation reflected more persistent differences among Notes: Additional descriptors include survey frequency (focal beaches were surveyed regularly, from May 2010 to July 2011; all others were surveyed once, in August 2010), beach type, beach orientation, latitude, and longitude. Swell exposure, beach width, and beach slope (WTO) are listed for all beaches from their August 2010 samples. See Table 2 for descriptions of variables.
v www.esajournals.org LiEbOWiTz ET AL. beaches (Table 1), as well as to characterize seasonal variation in wrack composition and cover within and among beaches. Percent cover (as a metric of abundance) was estimated for each wrack category, using a line-intercept procedure on each of three transects that extended perpendicular to the shoreline, from the edge of terrestrial vegetation or the bluff to the lowest intertidal level exposed by receding waves at each location (swash). The transects were assigned to locations within the first 100 m of shoreline from the access point using a random number table and a distance-measuring wheel. One edge of the track of a distance-measuring wheel was used to define a reference line for enumerating wrack cover and beach width. The extent and presence of each category of macrophyte was recorded along the reference line using size categories (1 mm to 8 m) yielding total wrack cover on the transect line and then expressed as m 2 of wrack per meter of shoreline by wrack category for each transect (m 2 /m of coastline, hereafter referred to as m 2 /m). These surveys characterized the spatial and temporal distributions of cover (m 2 /m) of a suite of wrack categories, of which six were dominant (Table 1): Nereocystis luetkeana, Zostera marina, Postelsia pal maeformis, Phyllospadix spp., mixed red algal species (a mixture of low intertidal and subtidal species including Mazzaella spp., Chondracanthus exasperatus, Cryptopleura/Hymenena, Poly neura, Microcladia spp., Plocamium cartilagineum, Ptilota/ Neoptilota, Mastocarpus spp., Halosaccion glandi forme, Erythrophyllum delesserioides, Prionitis spp., Neorhodomela larix and various articulated coralline algae), and mixed brown algal species (primarily Stephanocystis osmundacea, Desmarestia spp., and fragments of various kelps that could not be unambiguously identified).

Geographic variation in wrack composition and cover: donor ecosystem correlates
We identified three main donor ecosystems for our analyses: (1) kelp forests, (2) rocky intertidal zones, and (3) bays and estuaries. To determine the relative contribution of donor ecosystems in explaining observed distributions of the cover and composition of wrack among beaches, we developed a set of metrics that quantified the cover or proximity of each donor ecosystem to recipient beaches. We used multispectral aerial imagery to estimate the spatial extent of each donor ecosystem, which was acquired, processed, and interpreted by Ocean imaging (Oi) through the MPA baseline Program (see Svejkovsky 2013 for detailed methods). These data sets were categorized into intertidal substrates (20 spectral classifications), estuarine substrates (11 spectral classifications), and offshore kelp bed cover extent using supervised maximum likelihood and unsupervised iso cluster classification techniques in Esri inc. intertidal imagery was available for the entire NCC region, but kelp forest canopy cover was only available from Point Arena to San Francisco, creating a data gap south of San Francisco. The most recent kelp cover assessments prior to this, conducted by the California Department of Fish and Wildlife in 2009, were imported to explore replacing missing values, but given the high variation in yearly kelp cover, we decided to maintain those southern points as missing data for kelp cover assessments.
The Oi aerial imagery classifications were used to calculate metrics for kelp forest, rocky intertidal, and bay and estuarine ecosystems, as well as for two habitats within the rocky intertidal (Table 2). Kelp forest ecosystem cover was identified by a unique kelp spectral signal, and we calculated a metric of m 2 areal cover within 1 km radius of the beach wrack sampling site. Rocky intertidal ecosystem cover was identified as a "mixed red-brown" signal, and we developed a metric of m 2 cover within a 1 km radius of the wrack sampling site for this ecosystem as well. Two habitats within the rocky intertidal ecosystem were identified: (1) Phyllospadix habitat was identified as a unique spectral signal, with a metric of m 2 cover within a 1 km radius, while (2) a binary metric (presence or absence) of probable Postelsia habitat within the 1 km radius was developed by identifying rocky outcrops within the mixed red-brown category, and groundtruthing potential Postelsia habitat with the wrack survey project leads (one of whom [KJN] has done extensive fieldwork on Postelsia in the region). bay and estuarine ecosystems were characterized by two metrics: The first was a continuous metric developed using GiS distance tools to measure the linear shoreline distance from each recipient beach to the mouth of the nearest estuary or bay (proximity in m), and the second was a binary classification (presence or absence) of a v www.esajournals.org LiEbOWiTz ET AL. bay or estuary opening at the coastline within a 1 km radius of the wrack sampling site.
Characterizations were unique for each donor ecosystem, but we used a common spatial scale at which each donor ecosystem cover (areal extent, proximity, or presence/absence) was estimated. bull kelp (Nereocystis luetkeana) canopy cover and wrack were used for the selection of the common spatial scale, as it was the most thoroughly surveyed donor ecosystem, with the best estimates of areal canopy cover across the study region and period, to allow precise analyses. The total Nereocystis canopy cover (m 2 ) was calculated using the Oi shapefiles of Nereocystis canopy cover, and GiS software (ESRi ArcMap 10.2) buffer zone statistics. The scale with the strongest relationship between offshore Nereocystis cover and beach wrack cover of Nereocystis was determined by first estimating the areal extent of total canopy cover of Nereocystis forests within a set of buffer zones around each point where the wrack was sampled (length of radii of buffer zones = 0.5, 1, 3, 5, and 7 km). Then, we evaluated the strength of the relationships between these estimates of the total cover of offshore Nereocystis canopy for each of the five buffer zones and the cover of Nereocystis wrack on adjacent recipient beaches (individual analyses conducted for radii of 0.5, 1, 3, 5, and 7 km). These analyses used Pearson's correlation to test the relationship of the five radii, on two wrack data sets: (1) 17 beaches surveyed in August 2010, and (2) for the highwrack period averages (June-December 2010, described in Geographic variation in wrack compo sition and cover: regression analyses) of the 10 focal beaches surveyed monthly from 2010 to 2011. in August 2010, there were significant relationships among offshore Nereocystis canopy cover and beach wrack for radii of 0.5 km (r = 0.61, P = 0.02) and 1 km (r = 0.56, P = 0.04), but not for radius of 3, 5, or 7 km. in high-wrack periods, the pattern was similar and strongest for 0.5 km (r = 0.97, P < 0.0001) and 1 km radii (r = 0.91, P = 0.002), but it was also significant, although slightly weaker for radii of 3 km (r = 0.89, P < 0.01), 5 km (r = 0.88, P < 0.01), and 7 km (r = 0.86, P < 0.01). As the 1 km radius had the strongest statistical relationship outside of the 0.5 km radius, and previous studies of rafting in giant kelp (Macrocystis) have shown that drift kelp can move several km in the days before washing ashore or sinking (Harrold et al. 1998, Hobday 2000, the 1 km radius was chosen for subsequent analyses, and the remaining donor ecosystem cover and proximity metrics were calculated at this spatial scale for consistency (Table 2A). See Appendix S1 for maps of the beaches, the 1-km buffer zones, and the aerial imagery used for the analyses.

Geographic variation in wrack composition and cover: abiotic correlates
To explore the spatially explicit abiotic correlates with wrack cover and composition, we developed variables to represent marine macrophyte detachment and transport (swell and wind exposure), as well as deposition and retention (beach length, width, and slope; Table 2A). Swell exposure was estimated using 30 yr of seasonally averaged orbital velocities from the SWAN model (Storlazzi et al. 2005). We used the summeraveraged orbital velocities to extract average values within a 1 km radius of each sampled beach. The orientation of the beach was measured as compass degrees of the shore-normal line for each beach site (0/360° = North) and was a presumed proxy for exposure to prevailing northwest winds. The following transport data types were collated and assessed for utility in these models, but they were not available at the appropriate spatial scales, time frames, or locations for spatially explicit analysis: significant wave height (NDbC, NOAA), ocean surface currents (HF radar, iOOS), rip currents (time series satellite imagery, Google Earth), and local wind and tide data. We defined beach type as long (>1 km expanse of continuous sand between rocky outcrops) vs. pocket (<1 km stretch of continuous sand). beach width was measured from the lower edge of terrestrial vegetation (or the bluff, if no vegetation was present) to the lowest intertidal level exposed to breaking waves (the swash zone). beach slopes were measured at the hightide strand line (HTS; a metric of the daily hightide water level) and water table outcrop (WTO; the upper bound of saturated sand where the subaerial water table reaches the beach surface).

Geographic variation in wrack composition and cover: regression analyses
Spatial analyses were conducted on two sets of data to evaluate possible correlates of wrack cover and composition at two time periods. The v www.esajournals.org LiEbOWiTz ET AL. first set of models examined the wrack cover and composition data from August 2010, which included the largest number of beaches (n = 17), had the highest average wrack cover (2.7 m 2 /m) compared to the averages of the other monthly samples (0.06-2.6 m 2 /m), and included representation from five of the six major wrack categories assessed here (data set referred to as " All-beaches_ August"). The second set of models examined the potential spatial correlates of wrack distribution for the ten focal beaches with monthly samples, during the period of high wrack cover (June-December 2010; referred to as "Focal-High-Wrack"). This second set of models yielded more data-rich cover assessments, but it was limited to a smaller sample size of 10 focal beaches. Therefore, the two methods provided insights on the strength of different variables contributing to geographic variation in wrack composition and abundance, and are both presented for comparison. For all statistical analyses, data were tested for normality, and those indicating transformations (all response variables and offshore Nereocystis cover) were natural-log-transformed to meet assumptions of normality.
For both of the above modeling data sets ( All-beaches-August and Focal-High-Wrack), we developed general linear regression mixed models using AiC c forward entry model selection (SAS software [version 9.4]; SAS institute inc., Cary, North Carolina, USA). For both modeling data sets, individual regression models were developed for each of the six wrack categories as the dependent variable: (1) Nereocystis luetkeana, (2) Zostera marina, (3) Postelsia palmae formis, (4) Phyllospadix spp., (5) mixed red algae, and (6) mixed brown algae. Each model selection process considered the same suite of independent abiotic variables (swell exposure, wind exposure, beach width, HTS and WTO beach slope, and beach length [pocket vs. long]), and all potential donor ecosystems' variables for that wrack category. For example, Nereocystis wrack could only originate from kelp forest cover, and therefore, kelp forest was the sole donor ecosystem entered into the model. However, mixed red algae could originate from either rocky intertidal or the understory of kelp forests, and therefore, both donor ecosystem variables were included in the model selection process. Relationships for each significant predictor variable were plotted using non-log-transformed values for visualization.

Temporal variation in wrack composition and cover: contextual data
To consider the abiotic context of the monthly patterns of wrack composition and abundance, we collated available data sets with high temporal resolution. This precluded most data with high enough spatial resolution to model these monthly distribution patterns across beaches. Additionally, due to autocorrelation and only 1 yr of temporal data, it was not possible to develop statistical models for these monthly patterns. However, to visualize the physical context of wrack transport, we developed monthly averages for wind and swell exposure. Wind and swell variables were generated from the National Data buoy Center (NDbC) stations of Point Reyes (Station PRYC1, May-December 2010) and Point Arena (Station ANVC1, January-July 2011) merged, as none of the local buoys had complete records for the study period. We used these data to calculate 7-d moving averages for maximum wind speed, the average southern wind speed, and the average northern wind speed (Fig. 2c). The 7-d moving averages for orbital velocity ( Fig. 2d) were calculated from the same merged buoy data, with the following equation: where H is wave height, T is wave period, d is water depth (20 m), s is height above substrate (1 m), k is wave number (=2π/L), and L is wavelength (Denny 1988). This calculation of orbital velocity was selected to represent wave energy (swell exposure) as it integrates the influence of both wave height and wave period, and therefore indicates the stresses faced by subtidal algae (kelp and subcanopy algae) in the donor ecosystem.

Geographic and temporal patterns of wrack cover and composition
The distribution of wrack cover on the study beaches was highly variable, both geographically ( Fig. 1) and temporally (Fig. 2). Geographic variation in combined cover of the six wrack categories ranged three orders of magnitude (0.01-11.3 m 2 /m) across All-beaches-August (Fig. 1b), and there was no apparent systematic latitudinal gradient in wrack cover. beaches with high wrack cover were found throughout the region, with two beaches in the north (Cb and ST), one site in the middle of the region (DOR), and another in the south (Mb) (Fig. 1b); beaches with the lowest observed cover were often adjacent to those with the highest, such as Ab vs. Cb, and STG vs. DOR. These patterns of local variation in wrack cover were similarly reflected in the averaged Focal beach data (Fig. 1c).
There was marked variation in wrack composition. Mixed brown algae (especially Lessoniopsis littoralis) and Nereocystis were prevalent constituents of wrack at northern sites, while Zostera and Phyllospadix were prominent throughout the central section of the study region, and red algae and Nereocystis were prominent in the south (Fig. 1b, c). Across the All-beach-August data, the largest contributors to wrack cover were mixed brown algae (0.84 m 2 /m), followed by Zostera (Cb), the average cover of mixed brown algae was the fourth highest regionwide (0.32 m 2 /m). Average wrack composition in the Focal beaches was again similarly reflective of the All-beach-August data (Fig. 1b, c). Temporal patterns of wrack cover and composition were highly variable as well (Fig. 2a, b), with regional wrack cover (averaged across the Focal beaches) varying more than an order of magnitude among months, ranging from a January low of 0.07 m 2 /m to an August high of 2.7 m 2 /m. Composition also shifted temporally among wrack categories, with mixed red algae

Correlates with geographic variation in wrack cover and composition
The linear regression models generated to explain geographic variation in cover of the six wrack categories found that the local abundance of the donor ecosystem, attributes of the recipient beach, and metrics of transport processes all contributed to explaining geographic patterns of wrack cover, but the relative contribution of these categories varied both among the different wrack categories, and within each wrack category across the different data sets (Table 3). Generally, the models with strongest explanatory power were associated with the Focal-High-Wrack data.
Kelp forest ecosystem cover was the only significant explanatory variable for the models of geographic variation for Nereocystis wrack cover, as well as for the models of mixed brown algae wrack cover. For Nereocystis wrack, this was the case for both models (Fig. 3a, Table 3). Explanatory power was greatest for the Focal-High-Wrack model (adjusted R 2 = 0.80, P = 0.002). For mixed brown algae wrack, the strongest model was generated for the Focal-High-Wrack data (Table 3; adjusted R 2 = 0.94, P < 0.001), followed by the All-beaches-August data (Fig. 3b, Table 3; adjusted R 2 = 0.36, P = 0.015).
Rocky intertidal ecosystem metrics provided explanatory power for models of geographic variation of the cover of intertidal wrack categories (mixed red algae, Postelsia, and Phyllo spadix), while physical variables played a large role in these models as well (Table 3). For mixed red algae, the best model for the Allbeach-August data was predicted by the rocky intertidal (mixed red-brown) cover within 1 km (adjusted R 2 = 0.29, P = 0.015; Fig. 3c). However, for the Focal-High-Wrack data, it was predicted by a negative relationship with beach width (adjusted R 2 = 0.50, P = 0.013). For Postelsia, the only significant model was for the Focal-High-Wrack data, which was predicted by the presence of Postelsia habitat within a 1 km radius (adjusted R 2 = 0.57, P = 0.007; Fig. 3d). For Phyllospadix wrack cover, the only significant model was for the All-beach-August data (Table 3; adjusted R 2 = 0.43, P = 0.007), which included a negative relationship with maximum Table 3. Linear regression models predicting spatial cover of each of six wrack categories (column 1), assessed using two data sets: (1) All-beach-August data set, the full complement of 17 beaches with values from August 2010, and (2)  Notes: The independent variables included the associated metrics of spatial extent or proximity of the donor ecosystems (see Table 1), and the suite of physical metrics for physical drivers (swell exposure, wind exposure) and beach morphometrics (beach width, slope [HTS or WTO], and categorization as a long vs. pocket beach). swell exposure (Fig. 3e) and a positive relationship with WTO beach slope.
Estuarine and bay ecosystem metrics were the strongest explanatory variables for the models of geographic variation for Zostera wrack cover in the All-beaches-August data, while the Focal-High-Wrack model was predicted only by abiotic drivers (Table 3). The All-beach-August model (Fig. 3f) found the cover of Zostera wrack positively related to the presence of a bay or estuary within 1 km (adjusted R 2 = 0.41, P = 0.004). The Focal-High-Wrack model was predicted by a negative relationship with swell exposure (adjusted R 2 = 0.57, P = 0.007).

dIscussIon
We conducted a large-scale study to explore local and regional connectivity among the major donor ecosystems of macrophyte wrack (kelp forests, and bays and estuaries, and rocky intertidal ecosystems), and recipient sandy beach ecosystems. The abundance of marine macrophyte wrack is known to provide important energetic and nutrient subsidies to beach ecosystems ), but the scale and predictability of connectivity among specific donor and recipient ecosystems have been not been quantified. Our results suggest that strong regional patterns found in wrack subsidies to beaches are a result of the combination of local processes and larger-scale geographic patterns in the distribution of major subtidal and intertidal donor ecosystems. The strong temporal signal of composition and abundance of wrack on NCC beaches was influenced by both physical factors, and the phenology of individual macrophytes.

Spatial scales of variation in wrack cover and composition
Our surveys of beach wrack cover and composition revealed several general patterns. First, average wrack cover could vary over three orders of magnitude among adjacent beaches. This large range in wrack cover is consistent with studies that cited enormous variability in local deposition of wrack biomass, ranging from 360 to 2900 t·km −1 ·yr −1 in Western Australia (Hansen 1984, cited by Kirkman andKendrick 1997), 550 to 2660 t·km −1 ·yr −1 in Patagonia (Piriz et al. 2003), and 0.41 to 46.4 kg/m standing stock in southern California . Second, there was little latitudinal gradient in wrack cover evident at the scale of the entire study region. Thus, the data suggest that processes influencing the quantity of wrack cover on a specific beach involve local factors, which act at small spatial scales. These processes have been documented in other studies, which found local factors such as beach morphometrics (Orr et al. 2005, Duong and Fairweather 2011, Revell et al. 2011, Gómez et al. 2013 or the characteristics of the macrophytes such as buoyancy (Hobday 2000, Oldham et al. 2014) influenced wrack deposition on beaches.
While the abundance of wrack demonstrated a local scale of variation, the major constituents of beach wrack exhibited a larger regional scale of variation, broadly reflecting the regional patterns of distribution of the donor ecosystems. Geographically, the northern area of the NCC study region is characterized by relatively continuous rocky shores and subtidal reefs that support high production of the mixed brown algae and Nereocystis, which contributed disproportionately to wrack composition in this section of the coast. The central section of the study region is characterized by extensive embayments and estuarine ecosystems (e.g., San Francisco bay, bodega bay, bodega Harbor, Tomales bay, Drakes bay, Drakes Estero, Estero Americano) that support extensive beds of Zostera, the species which comprised a large proportion of the wrack on nearby beaches. Further south, beaches more removed from these sources of Zostera appear to receive a mix of wrack categories from rocky intertidal ecosystems (e.g., Postelsia, Phyllospadix), and subtidal kelp forests along that coast. Taken together, the composition of wrack across the study region reflects the distribution and relative cover of donor ecosystems that vary latitudinally across the study region (see Fig. 1 and Appendix S1 for maps to visualize these patterns). in addition to donor ecosystem distribution, physical metrics of macrophyte detachment, transport, and deposition/retention played a role in explaining patterns of wrack cover in beaches, discussed further below.

Spatial correlates with wrack composition and cover
Regression models were able to explain 31-94% of the variance in geographic distribution of wrack among beaches of the NCC, with local abundance, presence, or proximity of the relevant donor ecosystems by far the most prominent correlates. The strongest associations stemmed from the Nereo cystis forests as a donor ecosystem, which was correlated with wrack cover of both Nereocystis and mixed brown algae. Additionally, the local cover of rocky intertidal ecosystems was the strongest predictor for mixed red algae wrack cover, Postelsia habitat predicted Postelsia wrack, and the proximity to estuarine ecosystems predicted the cover of Zostera wrack. Variation in Phyllospadix wrack was not explained by the local availability of its donor ecosystem, suggesting it may have been more susceptible to physical drivers, or because our estimate of donor ecosystem abundance only captured intertidal, and not subtidal, Phyllospadix cover. Therefore, the largest contribution to discerning patterns in distribution was the local proximity or presence of the donor ecosystem. While this seems intuitive, the local cover or proximity of donor ecosystems could have played a lesser role if transport processes or beach morphology factors were dominant, or the scales of analyses were inappropriate. Thus, our analyses provide strong evidence for the importance of productive local donor populations in connecting sandy beaches to essential trophic subsidies.
The physical metrics describing macrophyte detachment, transport, and deposition/ v www.esajournals.org LiEbOWiTz ET AL. retention added explanatory power to some of the models of spatial variation in wrack cover and increa sed their relative contributions for the Focal-High-Wrack data sets. Zostera exhibited a negative relationship with swell exposure, which indicated a habitat association with lower energy systems. Phyllospadix displayed a negative relationship with swell exposure as well, although it shows a weak positive relationship with the WTO slope within the model, suggesting that it was cast high up on the beach at peak high tide by waves and deposited there. interestingly, we did not detect a significant relationship between wrack distribution and pocket vs. long beaches within these multivariate models, although others have found that sheltered pocket beaches had the highest amount of wrack (barreiro et al. 2011, Duong andFairweather 2011). The complex physical interactions impacting wrack transport and retention are only now being modeled in detail (Oldham et al. 2014).

Temporal patterns and correlates of beach wrack cover and composition
The cover of wrack on sandy beaches across the NCC exhibited enormous temporal variability as well, in keeping with highly temporally variable wrack distribution patterns elsewhere. Even at a daily scale in which loading is relatively limited, Dugan et al. (2011) found that deposition varied over an order of magnitude, from 0.1 to 5.6 kg wet wt·m −1 ·d −1 among sampling sites and dates on Southern California beaches. While most studies found high variability, some identified clear seasonal patterns, often with later summer or fall senescence (Stenton-Dozey and , Piriz et al. 2003, Revell et al. 2011, Gómez et al. 2013, while others saw more stochastic variability and nonseasonal interactions with the energy level of the system (barreiro et al. 2011, Goncalves andMarques 2011).
We found strong evidence for seasonal wrack deposition. For some macrophyte species with an annual life history and phenology, like Nereocystis, this would be expected. Nereocystis grows from early spring to fall and generally senesces in the winter, often dislodged in the first large winter storm (Springer et al. 2010), due to detachment and transport by ocean swell and winds (Harrold et al. 1998). Perennial species examined here (e.g., Phyllospadix) might be expected to exhibit a less pronounced seasonal signal. The intensive nature of beach ecosystem monitoring precluded more frequent surveys of wrack, which can be deposited or removed quickly by storms, tides, and consumption ). More frequent surveys may have documented higher frequency variability. However, the consistent patterns of negligible wrack from January to May, and variable high covers across the other months on all beaches, suggest that seasonality dominates across this region. Additionally, studies in southern California (Revell et al. 2011, Dugan andHubbard 2016) and ongoing studies in the region just north of Point Arena and up through Oregon found similar patterns of seasonality in wrack cover (Reimer 2014; K. J. Nielsen, J. E. Dugan and D. M. Hubbard, personal observation).
Our consideration of associations with the temporal variation of wrack cover at a regional scale could only include factors with sufficient temporal resolution; therefore, despite the importance of local donor ecosystem cover in the spatial assessments, our temporal considerations could only examine abiotic covariates, while accounting for phenology. For example, Nereocystis biomass that accumulates during the growing period (May-September, a period of high nutrient availability) becomes available for deposition as wrack in fall, due to susceptibility to breakage by the high wave energy that arrives with strong early-winter storms. This delivers large quantities of wrack, which align with a seasonal peak in wrack cover on the beaches (Fig. 2). With much of the biomass removed at the onset of storms, far less biomass is available for removal and delivery to beaches later in winter, consistent with the observed reduction in wrack cover later in the winter storm season. Thus, an understanding of phenology coupled with a visual comparison of seasonal patterns of wrack cover and these two transport processes suggests that seasonal patterns of standing biomass accumulation, removal from donor habitat, and delivery to beach habitat collectively help explain temporal patterns of wrack cover on beaches.

Trophic implications and management of pulsed and highly connected seascapes
Ecological subsidies can be vital to the productivity and diversity of recipient ecosystems by stabilizing populations and food webs v www.esajournals.org LiEbOWiTz ET AL. (Anderson et al. 2008), particularly through the delivery of temporally pulsed subsidies from multiple sources (Huxel et al. 2002, Anderson and and with different nutrient profiles (zimmer et al. 2004). beaches exemplify this conceptual framework; however, the consumption and turnover of wrack subsidies on beaches can vary greatly depending on the fauna, environmental conditions, and wrack category. The proportion of wrack subsidy consumed may also vary widely with the fluctuations in supply; for example, the large biomass of wrack delivered during high-deposition periods can swamp consumers, while during low-wrack deposition periods, the relative consumption of wrack can be rapid, with almost every palatable piece of wrack consumed each night. During times of scarcity, wrack supply may affect intertidal consumer populations more strongly than when supply is high. Consumption rates in southern California ranged from 87% of deposition per day for the giant kelp, Macrocystis, to a negligible 1% of deposition per day for Phyllospadix , and a study in british Columbia found similar rates, ranging from 90% consumption of deposition per day for Nereocystis to a low of 10% consumption of deposition per day for Phyllospadix (Mews et al. 2006). This variation in consumption rates among wrack categories, combined with the high temporal variability in deposition, highlights the extreme variability in energy and nutrient subsidies to this donordependent ecosystem. biodiversity and functional species redundancy in adjacent donor ecosystems may be necessary to buffer the impacts of highly pulsed subsidies from multiple sources on function and higher trophic levels in recipient ecosystems like beaches , and the drivers of temporal dynamics in these systems will benefit from further study. The abundance of macrophyte wrack is a strong correlate of the diversity and abundance of macroinvertebrates and in turn the abundance and diversity of shorebirds on beaches , Nielsen et al. 2013, Dugan and Hubbard 2016, presumably because the energy and habitat material of wrack fuel and support the base of the food web.
This study highlights the importance of multiple proximal donor ecosystems when considering ecosystem connectivity and resource subsidies.
There is a growing emphasis on the importance of looking beyond traditional definitions of ecosystems to consider ecological subsidies and connectivity, viewing multiple ecosystems as a "meta-ecosystem" to be considered as a larger whole (Loreau et al. 2003). Our results underscore the importance of incorporating landscape ecology methods developed in terrestrial realms into marine ecology and conservation (Carr et al. 2003, Jelinski 2015, Young and Carr 2015. Such awareness may help optimize management consideration of all factors that influence the condition of ecosystems, as exemplified by those that rely upon subsidies (McLachlan et al. 2013).
Adaptive, ecosystem-based management of MPAs and recent natural resource policies conceptually acknowledge the importance of managing entire ecosystems and protecting networks of connectivity among ecosystems, instead of focusing on particular populations or species at risk (Crook et al. 2015). However, the data, scientific tools and models, and the socioeconomic incentives for examining populations are much better established, and thus have remained a primary focus in planning and assessment of management actions, despite the emphasis on ecosystems in the language of the policies and laws that motivate the work. The tools and models for understanding and predicting the effects of management actions and connectivity on particular populations exist because the mechanisms driving the patterns of population dynamics are much better understood than those connecting ecosystems through the flow of energy, materials, nutrients, and detritus. The emerging theory and study of meta-ecosystems and their dynamics illustrates how important connectivity and subsidies between ecosystems can be in driving the structure and dynamic functioning of a single ecosystem. in this study, we provide a foundation for understanding the drivers, scale and degree of connectivity among donor ecosystems (kelp forest, rocky intertidal, and estuarine), to a recipient ecosystem, sandy beaches.
in particular, effective conservation and management of nearshore marine ecosystems (e.g., kelp forests) is an important component of ecosystem-based management for sandy beaches, as they provide a large portion of the energetic base of beach food webs , and the migrating and resident shorebirds that feed on them. Our findings also imply that sandy beach ecosystems will be sensitive to the vulnerability of kelp forests and intertidal macrophytes to climate change and commercial take (Steneck et al. 2002, Springer et al. 2010, Thompson et al. 2010, Sundblad et al. 2011, burnaford et al. 2014, and other factors that influence the productivity, extent, or proximity of their important donor ecosystems. Resource management policies and the design of MPAs should consider patterns of ecosystem connectivity to ensure protection of the donor ecosystems upon which recipient ecosystems rely.