Determining the relative influence of biotic and abiotic processes in structuring communities at local to large spatial scales is best understood using a biogeographic comparative-experimental approach. Using this approach, previous work suggests that intertidal community dynamics (top-down and bottom-up effects) vary unimodally along an upwelling-based productivity gradient, termed the Intermittent Upwelling Hypothesis (IUH). Evidence consistent with the IUH comes from the sessile invertebrate/predator (SIP) subweb in certain rocky intertidal communities, but whether this pattern extends to macrophyte/herbivore (MH) subwebs is unknown. Here we ask: Are MH subwebs also structured as predicted by the IUH? What is the relative importance of herbivory and predation in structuring these communities? Under what conditions do ecological subsidies like nutrients or propagule production drive community dynamics? And are omnivorous interactions important? We hypothesize that MH subwebs are driven by a new construct, the Grazing-Weakening Hypothesis (GWH), which states that MH interactions weaken monotonically with increasing nutrients, with strong (weak) herbivory and low (high) macrophyte productivity at low (high) nutrients. We explored local-to-large spatial scale dynamics of both subwebs using a biogeographic comparative-experimental factorial field experiment testing joint and separate effects of herbivores and predators between two continents. Experiments at ten sites ranging across from persistent upwelling to persistent downwelling regimes ran for 26-29 months in Oregon and California, and New Zealand South Island. For the MH subweb, results were consistent with the GWH: herbivory declined and macrophytes increased with increasing nutrients. As expected, results for the SIP subweb were consistent with the IUH: predator effect size was unimodally related to upwelling. Overall, herbivory explained more variation in community structure than did predation, especially in New Zealand. Omnivory was weak, sessile invertebrates outcompeted macrophytes, and ocean-driven subsidies provided the basic template driving ecosystem dynamics. We propose a unified Meta-Ecosystem Dynamics Model (MEcoDynaMo) combining MH and SIP results: with increased upwelling, sessile invertebrates and underlying dynamics vary unimodally (as in the IUH), while herbivory decreases and macrophytes generally increase. While this model was based on research in temperate ecosystems varying in upwelling regime, its wider applicability remains to be tested.

Study Systems

Research was done on wave-exposed rocky benches on the OR and northern CA coast located in the California Current Large Marine Ecosystem (LME), and NZ, located in the New Zealand Shelf Large Marine Ecosystem (Fig. 2). Experiments were done in OR and NZ, and supplementary data on nutrients and macrophyte cover were obtained in CA by the PISCO consortium. Tidal ranges were 3.0 to 3.5 m in OR and CA and 2.1 (east coast) to 3.4 m (west coast) in NZ (Appendix S1: Table S1). OR/CA experiences mixed semi-diurnal tides (i.e., successive low or high tides are of unequal height, with one “good” low tide per day during spring tide periods), while NZ tides are semi-diurnal (daily low and high tides are of similar height). Except for Nine-Mile Beach in NZ, sites used in this study have lengthy histories of continuous ecological research, beginning in 1983 in OR, the 1990s in CA, and 1994 in NZ.

Environment

Experiments in NZ and OR were conducted between October 2004 to October 2007 and May 2005 to August 2007, respectively. CA and NZ data on nutrients and chlorophyll-a were obtained in 1999-2002 and have been collected in OR continuously since 1999. The research included environmental and ecological measurements and identically designed and conducted experiments in both LMEs, with initiation offset by 6 months so that each was begun in the respective spring in each LME (Appendix S2: Fig. S1).

Upwelling Regimes and Shelf Width. Oceanic conditions in the two LMEs have been detailed elsewhere (e.g., Menge et al., 1997, 1999, 2003, 2011, 2015; Menge & Menge, 2013; Schiel, 2004, 2011; Stevens et al. 2021). Briefly, the OR/CA coast is characterized by intermittent upwelling, with a reversing, conveyor belt-like cross-shelf flow pattern (Appendix S1: Table S2). Upwelling events typically deliver temporally varying pulses of cold, nutrient-rich water nearshore and then via subsequent cross-shore flows to offshore shelf waters, while the alternating downwelling-driven events deliver warmed, nutrient-poor surface waters to the nearshore. Periods between upwelling cessation and downwelling onset are termed “relaxations” and are also conditions that can deliver materials to the coast as the uneven sea level generated by upwelling processes evens up with shoreward water movement. In this study, upwelling was stronger at southern OR and CA sites, and weaker at the four central OR sites.

Coastal currents interact with shelf width to modify retention of subsidies along the coast (Kirincich et al., 2005). The shelf is relatively wide in the Cape Perpetua (CP) or “central” region (sites SH and TK) and relatively narrow at all other sites (Menge et al., 2015), generating weak and retentive currents in the CP region and stronger, less retentive offshore currents elsewhere (see bathymetry in Fig. 2).

In NZ, East coast sites are dominated by persistent downwelling while West coast sites experience intermittent upwelling (Menge et al., 2003; Appendix S1: Table S2). Shelf width also differs between coasts, narrow at the East coast sites and wide at the West coast sites (Menge et al., 2003). It thus seems likely that the West coast sites have retentive oceanic regimes for the same reason those at CP; the wide shelf likely weakens upwelling currents. On the East coast, downwelling flows occur shoreward, so would also transport particles toward the coast, at least in surface waters.

We quantified upwelling regimes using two measures, a high resolution (to 0.25o [28 km] latitude) dataset used in analyses of OR and NZ data in which we used the Ekman cross-shelf transport component (e.g., see Close et al. 2020), and in analyses including CA data, the globally available Bakun upwelling index (Bakun, 1975). Units of both are m3 water * sec-1 * 100m of coastline. These metrics were highly correlated (p < 0.0001, adj. R2 = 0.92) indicating using these alternatives should give comparable results. Although an updated upwelling index (CUTI, Jacox et al., 2016) is available for California Current sites, it has not yet been extended to NZ. We calculated the monthly average upwelling index at each site for the period of study in each LME. Upwelling intermittency was calculated as in Menge & Menge (2013; p. 290 and Appendix B: Figs. B2, B3) and shelf width was determined as detailed in Menge et al. (2015) for both OR and NZ. Data for 100m isobaths were obtained using the ETOPO1 database (see Appendix S1: Table S2)

Water Temperature. Seawater temperature was quantified using Onset HOBO TidBit ® (Onset Computer Corp., Bourne, MA, USA) replicate temperature loggers (n = 2 or 3) placed in situ in the low intertidal zone at each site. Loggers recorded at 15-minute (OR) to hourly (NZ) intervals, which were averaged to daily, then to monthly values for analysis. We used a detiding program (Menge et al., 2003) to separate air from water temperatures (hereafter termed intertidal sea surface temperature or ISST).

Nutrients and Phytoplankton. Nutrients (NO3 + NO2) and phytoplankton (proxied by Chlorophyll-a, hereafter Chl-a) were quantified using replicate (n = 3-5) bottle samples taken during low tides (see details in Menge et al. 1997, 1999, 2003, 2015). We used data averaged across upwelling months in OR (April-September in OR, October-March in NZ). Nutrient data were only available from 1999-2002 for NZ so we limited our comparisons of nutrients to upwelling months in those years in both NZ and OR. 

Prey Recruitment. We considered prey recruitment to be an external (environmental) input in our analyses (see below). We quantified barnacle and mussel recruitment using standard collectors to test for associations between recruitment and experimental communities. Barnacle recruitment was determined using 10 x 10 cm PVC plates covered with Saf-T-Walk ® (3M), a rubbery rugose tape used as an anti-slip surface on boats. Such collectors have been widely used (e.g., Farrell et al., 1991; Broitman et al., 2008; Dudas et al., 2009; Menge, 1992; Menge et al., 1999, 2003, 2010, 2011, 2015; Menge & Menge, 2013, 2019; Navarrete et al., 2005, 2008; Pfaff et al., 2011, 2015), and provide an index of barnacle recruitment. Mussel recruitment was similarly determined using plastic mesh balls (Tuffies), also in wide use (references in previous sentence). Collectors (n = 5 replicates each) were placed in the lower mid intertidal zone near the experiments at all sites in both LMEs and replaced monthly in OR and monthly to quarterly in NZ. Samples were processed in the laboratory by counting recruits under a microscope. Barnacle recruits could usually be identified to species, but Oregon mussel recruits were difficult to identify to species so likely included both Mytilus californianus and M. trossulus. New Zealand mussel recruits of P. canaliculus and A. ater maioriensis were distinguishable but M. galloprovincialis and Xenostrobus pulex are very similar. In analyses, we combined all acorn barnacle and mussel recruits as “barnacle recruits” and “mussel recruits”, respectively. We used data from the recruitment seasons (~April-November in OR, October-March in NZ) for 2003 through 2006.

Similar Biota and Functional Roles among Systems

The biota in each LME consisted of similar groups including some of the same genera, but few species occurred in both LMEs (Appendix S1: Table S3). Each LME had similar functional groups of macrophytes, sessile invertebrates (mytilid mussels, acorn and stalked barnacles), grazing herbivores (limpets and chitons), and predators (whelks, and sea stars). Further, each functional group consisted of morphologically similar species; both LMEs had larger and smaller species of mussels and barnacles, limpets and whelks of similar size and appearance, and similar apex sea star predators known to be dietary generalists and to consume similar prey (for details see Paine, 1971; Menge et al., 1994, 1999, 2003; Navarrete & Menge, 1996; Schiel, 2004, 2011; Novak, 2010, 2013). A few taxa occurred in, or were abundant in, one system but not the other. For example (1) surfgrass (e.g., Phyllospadix spp.) does not occur in NZ (but did not colonize our OR plots); (2) gooseneck barnacles Pollicipes polymerus are abundant in OR, but comparable stalked barnacles Calantica villosa and C. spinosa were sparse at NZ sites and mostly restricted to cryptic habitats; (3) the large chiton Katharina tunicata is abundant in OR but similar-sized chitons (e.g., Eudoxochiton nobilis) were sparse at our NZ sites; and (4) herbivorous fish and crabs can have effects at some sites in NZ (e.g., Rilov & Schiel, 2006; Taylor & Schiel, 2010). We have never observed crabs at our NZ sites, nor have we seen evidence of their activity or that of fish (e.g., broken or cracked mussel shells or bitten off barnacles), so feel it unlikely they had effects in our experiments. Crabs are sparse at our OR/CA sites, and possible fish predators seem restricted to calmer waters, and again, we’ve seen no signs of predation by these taxa at our sites. Despite these differences, the similarities remain striking and researchers familiar with one system quickly become comfortable working in the other.

In this study, because few species occurred in both systems, we focused on interactions at the functional-group or “trait-based” level. The taxa listed in Appendix S1: Table S3 were sorted into ten functional groups including crustose algae, filamentous algae, turf-forming algae, blade-forming algae, large brown algae, chthamaloid barnacles, balanoid barnacles, gooseneck barnacles, mussels and “other” (consisting mostly of sea anemones).

Community Patterns among Upwelling Regimes

Abundances. To provide the broader community context for the experiments, at all sites we quantified community structure, defined as abundance of sessile and mobile species, using three methods. In NZ, we conducted horizontal (i.e., parallel to the water’s edge) transect-quadrat surveys for sessile organisms and small mobile organisms (e.g., Menge, 1976; Lubchenco & Menge, 1978). In OR and CA, we conducted vertical transect-quadrat surveys (i.e., perpendicular to water’s edge, spanning the intertidal zones). In the horizontal transect-quadrat method, 0.25 m2 quadrats divided into 0.04 m2 subquadrats were placed at three-meter intervals along 30 m transect tapes placed in the center of the low-mid zone (i.e., the ~ 1 m-wide zone just below mussel beds) parallel to the water’s edge. Abundance, quantified as percent cover of all sessile taxa, was estimated from photographs using the subquadrats (each covering 4% of the 100% of the quadrat area) as a guide to facilitate abundance scaling to the 1% level. Mobile organisms were then counted in each 0.25 m2 quadrat to provide density estimates. Mobile organism abundance estimates in vertical transect-quadrat photographic surveys (OR, CA) were done similarly; although transects ran from the low-mid zone to the high zone, we used only data from the low zone (i.e., the ~ 1 meter below the mussel bed). In OR and NZ, we also used larger belt transects in areas of sea star habitat (typically just below mussel beds or in channels) to quantify sea star abundance (e.g., Menge et al., 2011, 2016). Belt transects consisted of replicated 2 x 10 or 2 x 5 m plots placed parallel to the water’s edge in which all sea stars were counted, weighed and measured. Here we present averages of percent cover data across 2015-2020 (6 surveys, OR and CA) and 1995-2016 (11 surveys, NZ), small mobile species (including limpets and chitons) density data from 2006-2016 (11 surveys, OR) and 2008-2016 (3 surveys, NZ), and sea star density data from 2005-07 (3 surveys, OR), and 2000-2019 (4 surveys, NZ).

Limpet Size. Because herbivore size could be related to herbivore impact as well as herbivore density, we quantified limpet size using two methods. In NZ, we measured species-specific limpet length directly in the field at each site using rulers and calipers. All accessible limpets (some were in crevices) in the vicinity of each transect were measured with the goal of measuring at least 200 individuals of each species. In OR, we measured limpets in photographs of experimental treatments (see below) using rulers placed at the edge of each plot as a scale. We used photographs taken in the second year of the experiment (i.e., after limpets recolonizing plots had grown to sizes like those of limpets in the vicinity) to ensure that sizes were representative.

Comparative Experiments Testing Top-Down Control Among Upwelling Regimes

Experimental Timing. In OR, the experiments testing predator x herbivore effects started in May 2005 at six sites, two each in three regions (from north to south, Capes Foulweather [hereafter northern region], Perpetua [central region], and Blanco [southern region]. In NZ, experiments began in October 2004 at four sites, two each on East and West coasts. Depending on the site, and how quickly an endpoint was reached (i.e., when little change occurred during several consecutive monitoring visits), experiments ran for 21-26 months in OR and 22-29 months in NZ (Appendix S2: Fig. S1). Although sampling was often monthly, for easier analysis and comparison, we averaged samples by season (termed a “sample period”) in each LME with OR sample periods lagged 6 calendar months after NZ sample periods. In sample periods 2-7, seasons were summer, fall, winter, spring, summer and fall. Because recovery was slow on the NZ East coast, the last two sample periods (8 and 9) were six and 12 months long (Appendix S2: Fig. S1).

Experimental Design. All experiments were established in the lower-mid zone, transitional between typically mussel-dominated mid zones and algal-dominated low zones. Plots ranged in tidal height from about 0.8 to 1.8 m (Appendix S1: Table S1). Reflecting the different tidal range between NZ east and west coasts, experiment heights were slightly lower at east than at west coast sites (Appendix S1: Table S1). The basic unit was a 25 x 25 cm square experimental plot, which was cleared of all biota using scrapers and wire brushes. Oven cleaner (NaOH) was then applied to remove most remaining algal and animal tissue. Treatments were (1) controls that allowed herbivore and predator access (abbreviated = +H+P), (2) predator exclusion (herbivores present, predators absent = +H-P), (3) herbivore exclusion (herbivores absent, predators present = -H+P), and (4) consumer exclusion (herbivores absent, predators absent = -H-P). Three control treatments, all +H+P, were included. To control for stainless-steel fences and antifouling paint (see below), we used partial (two-sided) fence controls (FC) and partial (strips with gaps) paint control barriers (PC). Plots with stainless steel lag screws at each corner served as marked plot controls (MP). MPs, FCs, and PCs allowed entry by limpets, chitons, whelks, sea stars, and crabs. The efficacy of the paint exclusion and fence methods to exclude or allow entry by limpets, the main herbivore group, was high (Appendix S2: Figs. S2, S3). Limpets were abundant in control (+H+P) and predator effect treatments (+H-P), and low in herbivore effect and consumer effect treatments (–H+P and –H-P). The design was blocked, with each block including all treatments, and with five replicate blocks per site. Blocks were spaced out on the shore, with 3-10 meters between each block.

Predator effects (+H-P) were tested by fastening four-sided stainless steel mesh fences or “open cages.” Fences hindered entry by certain benthic predators (whelks, sea stars; but not crabs or fish) but allowed grazers (limpets, chitons) to enter (by crawling underneath cage edges) or recruit to the cages. Those few whelks entering cages (species shown in Appendix S1: Table S3) were removed at each monitoring visit. The only sea stars found in OR cages were one Pisaster ochraceus at SH (once), and occasional Leptasterias sp. at FC. No sea stars were found in NZ cages during monitoring visits.

Herbivore effects (-H+P) were tested by attaching a 3-5 cm wide strip of Z-Spar® (Pettit Paint A-788 Splash Zone marine epoxy) around each plot, smoothing the edges to ensure a continuous smooth surface was available to grazers, and applying a layer of copper-based antifouling paint (Coastal Copper 250 Ablative Antifouling Bottom Paint, 1st Marine Products) (e.g., Fig. 3). Prior research has shown that this effectively excludes limpets but does not impede whelks, sea stars, or crabs (Cubit, 1984; Menge, 2000; Freidenburg et al., 2007; Guerry, 2008; Guerry et al., 2009; Guerry & Menge, 2017). Herbivorous snails (e.g., species shown in Appendix S1: Table S3) are also not deterred by the paint, but at all OR and NZ sites were either sparse (Tegula, Lunella, Melagraphia, Diloma) and/or very small and unlikely to have much effect (littorines, Risselopsis). Consumer effects (-H-P) were tested by combining paint strips and stainless-steel fences (e.g., Fig. 3). Fences were attached to the rock and then surrounded by the Z-spar/paint barriers.

After initiation, abundances in plots were determined visually in the field or photographically during each monitoring visit by the senior author. Because prior experiments showed that changes typically were initially rapid and slowed over time, the frequency of visits varied, ranging from monthly during the first months of the experiment and extending to three (OR) to twelve (NZ) months. During monitoring visits, we removed invading grazers or predators from exclusion plots, repainted barriers as needed, and repaired damaged cages. In total, we analyzed 4,711 samples across all combinations of LME (n = 2) x Region (n = 5) x Site (n = 2 per region) x Treatments (n = 6) x Replicates (n = 5) x Sample dates (7 to 9 in NZ, 17 to 19 in OR, depending on region).

Data Analysis

To organize our results and categorize datasets, we have assembled each according to region, inclusive dates, the sites used, and in which figure(s) and/or table(s) analyses each dataset was used into Table 1. The table heading lists all sites used including core sites in each region and sites from which supplementary data were collected.

Data were analyzed using JMP v16.1.0 (SAS Institute, Inc., 2021), PRIMER 7, and PERMANOVA+ for PRIMER (Anderson & Gorley, 2008). We used natural log-transformed data (ln (x+1)) in JMP analyses and square root transformations in PRIMER analyses. We chose not to use arc-sine transformations because percent covers could range >100% due to algal layering or overgrowth. Factors in the PERMANOVA model included the nested spatial effects crossed with the time and treatment effects, including LME (large scale effects of OR and NZ as a random factor), region nested within LME (mesoscale effects of cape or coast as a random factor), and site nested within region (local-scale effects of site, often as a random factor), predator treatment, herbivore treatment, time point and their interactions. Because large LME effects masked the experimental effects of interest, we ran subsequent PERMANOVAs for each LME (as above without LME). Estimates of variance components converted into percent variance explained (see Table 2 heading) were used to examine the relative importance of each factor and interaction. We visualized temporal changes in experimental communities using nMDS (non-metric Multidimensional Scaling) to plot centroid-based vectors of averaged consumer treatments among sites through time show community development patterns.

Inclusion of block and all three control treatment in the overall analysis did not allow the model to run due to insufficient degrees of freedom. We therefore conducted preliminary analyses testing for differences among the three control treatments and found no differences among them either as main effects or in interactions (Appendix S1: Table S4). Similarly, we tested block (i.e., replicate) effects and found it was not significant as a main effect or in most interactions except (the Site(Region(LME)) x Block interaction; p = 0.0001, 11.6% of the variance; Appendix S1: Table S5). This decision likely made our relevant analyses conservative, with some factors identified as not significant when in fact they might have been.

We used two methods to estimate effect sizes of top-down control, assess upwelling regime influences on these effect sizes, and determine if results vary with procedure. In Method One, we conducted nMDS on functional group abundances by site and treatment (averaging MP, PC, and FC to obtain a single value for the +H+P treatment) and saved x and y coordinates for centroids of each by treatment and replicate (control [+H+P], herbivore effect [-H+P], predator effect [+H-P] and consumer effect [-H-P]). We used these coordinates to calculate vector distances between control and exclusion treatments which served as our estimate of effect size. We visualized the association between all site-by-treatment combinations in 2D space using nMDS and overlaid environmental metrics to assess the relationship between each site-by-treatment combination and the environment. Finally, we regressed herbivore, predator and consumer effect sizes against upwelling to determine if the experimental results were consistent with expectations of the combined IUH/GWH model (Fig. 1).

In Method Two, we summarized temporal variation in community structure in each treatment and site by calculating how individual functional groups varied during the experiment. After square-root transforming percent cover of all functional groups, we standardized each to its maximum value, then took the average across the groups by each replicate, site, treatment and sample period to calculate a community index. We then examined linear, quadratic and cubic regressions (by site, treatment) for the average community index vs. time and found that cubic regressions provided the best fit to the data. Next, we used the second-order coefficient (i.e., the slope) of the cubic regressions to estimate the effect sizes of herbivore exclusion, predator exclusion and consumer exclusion, and plotted the slopes against the average Bakun upwelling index at each site. We also plotted the third-order coefficient vs. upwelling to determine the rate of deceleration or acceleration of the effects of each consumer group. To visualize which taxa were most closely associated with patterns observed in the cubic regressions, we used nMDS and plots of temporal changes of each functional group.

We used Principle Coordinates Analysis (PCO) to ordinate the environmental regimes (upwelling, intermittency, shelf width at 100m dept, ISST, mussel and barnacle recruitment, limpet density and shell length) to determine which environmental factors were most closely associated with each site. We assessed the association between average experimental results using nMDS with an overlay of environmental factors and distance-based linear models (DistLM).

Excel or similar spreadsheet.