DOI: 10.5061/dryad.3ffbg79p8
title: A unified meta-ecosystem dynamics model: Integrating herbivore-plant subwebs with the intermittent upwelling hypothesis

description: 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.
description: <strong>Study Systems</strong>
 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.
 <strong>Environment</strong>
 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). <em>Upwelling Regimes and Shelf
 Width</em>. 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 &lt; 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) <em>Water
 Temperature. </em>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). <em>Nutrients and
 Phytoplankton.</em><strong> </strong>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.  <em>Prey
 Recruitment.</em><strong> </strong>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. <strong>Similar Biota and
 Functional Roles among Systems</strong> 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). <strong>Community Patterns among
 Upwelling Regimes</strong>
 <em>Abundances.</em><strong> </strong>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).
 <em>Limpet Size.</em> 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. <strong>Comparative Experiments
 Testing Top-Down Control Among Upwelling Regimes </strong>
 <em>Experimental Timing.<strong>
 </strong></em>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). <em>Experimental Design.
 </em>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). <strong>Data
 Analysis</strong> 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 &gt;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).
description: Excel or similar spreadsheet.