Latitudinal variation in the constitutive and inducible defenses of a canopy-forming rocky intertidal seaweed

A long-standing theory in biogeography is that stronger biotic interactions at lower latitudes select for better-defended phenotypes. However, greater environmental variability at temperate latitudes may also shape defensive strategies by increasing temporal or local spatial variation in consumer pressure and thus selecting for greater phenotypic plasticity, or inducible defenses. Distinguishing between inducible and constitutive defense strategies is therefore necessary to test for latitudinal defense gradients, but also for understanding how species interactions and community dynamics vary across ecological and evolutionary scales. 

We investigated latitudinal variation in antiherbivore defenses of a cosmopolitan rocky intertidal seaweed (Fucus vesiculosus) against a similarly common and abundant grazer (periwinkle snails, Littorina littorea). We used a multi-phase common garden experiment with seaweeds from three different regions along the US Atlantic coast and snails from the southern- and northern-most study sites. We manipulated and measured snail grazing on seaweeds in a series of no-choice and two-choice assays to identify regional differences in seaweed functional defenses (i.e., reduced grazing). 

Across all assays, grazing rates declined with seaweed latitude. Prior grazing reduced the palatability of southern, but not northern, seaweeds. Changes in seaweed nutritional content (C:N) and phlorotannins (a putative chemical defense) correlated with induced but not constitutive functional defenses.  

Our results indicate that the constitutive anti-herbivore traits of F. vesiculosus increase with latitude and negatively covary with defense plasticity. This result suggests that the selective pressure of herbivory is stronger for northern seaweed populations, while southern populations may face a different set of tradeoffs leading to defense plasticity, such as increasing environmental stress.  

The strength of trophic interactions plays an important role in community dynamics and food web stability. Our findings add to a growing literature highlighting the importance of ecological context in shaping trophic interactions and suggests that estimates and comparisons of interaction strength need to consider spatiotemporal variation in prey defenses. Defensive traits that vary with latitude or along environmental stress gradients may be particularly important for predicting the effects of climate change on trophic interactions and their consequences for community dynamics and ecosystem function. 

Methods
Organism Collections- We collected F. vesiculosus seaweeds from nine sites spanning >750km of NW Atlantic coastline and three regions in and around the Gulf of Maine (Fig. 1). Seaweed individuals (one 10-20cm thallus and its holdfast) were collected from the mid-intertidal zone, removing each from the substrate at its holdfast with a metal scraper. Individual holdfasts were separated by >0.25m to avoid sampling genetically identical thalli (N= 80-100 individuals site-1). We avoided seaweeds that were reproductive or visibly damaged. Upon collection, we froze subset of seaweeds (n= 10 site-1) at -20°C for later analysis of ‘in situ’ tissue chemistry. L. littorea snails were collected from one southern and one northern field site (N= 400 site-1; shell length= 15-25 mm). All snails and seaweed were maintained under ambient light and seawater conditions in shallow flow-through tanks (water temperature range 19-21°C). Snails were kept in separate tanks and fed ad libitum with non-experimental sympatric Fucus until 24h before starting experiments.  

Relaxation Phase – Seaweeds were cleared of all epiphytes and epifauna and maintained in the absence of herbivores for 10d beginning 21August 2019. This “relaxation phase” allows for defenses that may have been induced in the field to relax or dissipate (Underwood et al., 2002), which can take as little as 4d (Rohde & Wahl 2008). After 10d, samples of “relaxed” seaweeds (n= 6-10 site-1) were frozen as above, while the remainder were used in the induction phase (n= 56 site-1).  

Induction Phase - Pairs of relaxed, sympatric seaweed individuals were maintained in independent flow-through “induction units” in the presence or absence of herbivorous snails for 12d. Each induction unit consisted of a 1.2L mesh-sided container submerged within a 2.5L bucket receiving its own supply of running seawater. A mesh panel divided the inner container to separate the two individual seaweeds. Induction units were assigned randomly to one of three treatments: Southern Inducer Snails, Northern Inducer Snails, or No-Snail Control (n= 7, 7, or 14 site-1, respectively; N= 252 induction units). Units assigned to inducer snail treatments received 3 snails from the respective source population in one side of the container, allowing the other seaweed to serve as an autogenic growth control (Long et al. 2013). The second seaweed in the no-snail controls ensured that total seaweed biomass was similar across all induction units and treatments (see Fig. S1).  

Each seaweed’s tissue mass was measured at the beginning and end of the induction phase (mean ± SE= 1.068 ± 0.014g, N= 504). Individuals were blotted dry with paper towels to remove all visible moisture prior to weighing. Inducer snails (n= 3 unit-1) were also blotted dry with paper towels and weighed, as a group, for total biomass per unit (Bunit). At the end of the induction phase and after measuring final mass, two apical tips were clipped from each seaweed for use in choice assays (Sotka et al., 2002, Taylor et al., 2002), with additional tips sampled and frozen for later analysis of “post-induction” tissue chemistry.  

The total amount of seaweed grazed by snails within a given unit (Gunit) was estimated as Si(Af/Ai)-Sf, where Si and Sf are the initial and final mass, respectively, of the seaweed with snails, while Ai and Af are those of the paired autogenic control (Long et al., 2013). Because northern snails were slightly larger than southern snails (mean ± SE individual mass, 3.82 ± 0.05 vs. 3.16 ± 0.04 g), we also calculated mass-adjusted, per capita daily grazing rates (Gpc) normalized for an individual 3.5g snail using the formula 3.5*(Gunit/Bunit)/Dunit, where Bunit is the total biomass of the three snails in the unit, and Dunit is the induction phase duration (11.6 – 12.0 d). We analysed Gunit and Gpc with separate linear mixed models (“LMMs”) that included seaweed region, inducer snail population, and their interaction as fixed effects. Seaweed collection site (within region) and the site x inducer interaction were included in the model as random intercepts to account for the nested experimental design. Seaweed growth rates in the no-snail control treatment were calculated as (Cf-Ci)/Ci where Cf and Ci are the final and initial mass, respectively. We calculated then analysed the average growth rate in each induction unit using a similar LMM that included region (fixed effect) and seaweed collection site (random effect). Seaweeds from 5 units were excluded from analyses and choice assays: two units experienced water flow issues and mortality and three units were missing seaweeds. 

Choice Phase – We conducted two types of choices assays to assess differences in relative palatability (1) between regions (“inter-region” assays) and (2) between induction treatments within each region (“intra-region” assays). Choice assays were conducted in the same flow-through setup as the induction phase but with single-compartment mesh containers (“choice units”). Each choice unit contained a pair of apical tips from seaweeds exiting the induction phase and either three northern or three southern “chooser” snails (Fig. S1). Seaweed collection sites within each region were pooled and seaweeds were chosen at random for each choice replicate. (1) Inter-region assays offered snails a choice between one lower- and one higher-latitude tip from no-snail control seaweeds. Three paired region choices (lower/higher: South/Central, South/North, Central/North) were crossed with the two chooser snail populations (N, S) in a 3 x 2 design with n= 12 choice units for each choice x chooser combination (N= 72 inter-region choice units). (2) Intra-region assays offered snails choice between two seaweed tips from the same region but different induction treatments: one tip was from a seaweed in the no-snail control treatment and the other was from a seaweed in either the N- or S-inducer snail treatment. Two paired control/inducer choices (Control/N Inducer, Control/S Inducer) for each of the three seaweed regions were offered to N and S chooser snails in a 2 x 3 x 2 factorial design. Seaweeds from two induction units (S seaweed in the N inducer treatment) could not be included in choice assays because high grazing rates during induction removed most or all apical tissues. In addition to mortality losses described above (see “Induction Phase”), this led to an unbalanced design with 7-12 replicates for each inducer x region x chooser treatment combination (N= 138 intra-region choice units).  

In both choice assays, the different choice combinations (e.g., South/Central, Control/N-induced, etc.) were replicated by randomly selecting one seaweed (i.e., induction unit) from each of the two corresponding regions or treatments. Source sites within each region were pooled and selected randomly. One apical tip from each seaweed/induction unit was then placed in each of two choice units, which were then randomly assigned to receive either northern or southern chooser snails (see Fig. S1 for an example). In this way, we ensured that any among-unit differences (i.e., variation among individual seaweeds) within each induction treatment were distributed uniformly between the two chooser snail populations.  

All choice assays began on 13 September and lasted 3-6d. Units were inspected daily to visually assess grazing on apical tips. Assays were terminated on a unit-by-unit basis when either tip appeared ~50% consumed, or after a maximum of 6d (>70% of units crossed the 50% threshold by day 6). The flexible assay duration minimized potentially confounding effects of preferred resource limitation (i.e., grazers consuming all of a preferred tip) or seaweed trait changes (i.e., defense relaxation or induction) during the choice assay.  

Seaweed tips were individually weighed immediately before and after each assay to calculate tissue mass grazed by snails (Gtip= initial-final mass) and standardized per capita daily grazing rates as Gpc= 3.5g*(Gtip/Bunit)/Dunit, where Dunit is the choice unit-specific assay duration (3-6d). Grazing rates on each tip were analysed with LMMs that included choice unit, seaweed source site, and induction unit as random intercepts. Grazing rates were square root transformed prior to analysis to satisfy assumptions of normality and homoscedasticity. Analysis of inter-region choice data included seaweed tip relative latitude (low or high), chooser snail population (N or S), seaweed region combination (S-N, S-C, C-N), and their interactions as fixed effects. Intra-region analyses included seaweed tip induction treatment (induced or control), inducer snail population (i.e., induced tip exposed to N or S grazers during the induction phase), chooser snail population (N or S snail grazing during the choice assay), and all interactions as fixed effects. We then conducted joint contrast tests on estimated marginal means (R package emmeans; Lenth 2023) to test for main and interactive effects of induction treatment and inducer snail population within each region x chooser combination.  

Tissue Chemistry - Seaweed apical tips from different experiment phases were freeze dried, homogenized, and ground into fine powders. We extracted phlorotannins from pre-weighed subsamples (~15mg) following Kurr & Davies (2019) and spectrophotometrically analysed liquid extracts following Pavia & Toth (2000) using Folin-Ciocalteu reagent with phloroglucinol (1, 3, 5-trihydroxybenzene, Sigma-Aldrich P3502) as a standard. An additional subsample of Fucus powder (~1.5mg) was weighed on a semi-microbalance (± 0.001mg) and passed through a Fisons NA 1500 Series 2 Elemental Analyzer (Costech Analytical, Inc.) to estimate carbon and nitrogen mass based on curves produced using acetanilide standards (C8H9NO, Costech). C:N ratios and phlorotannin concentrations (proportion dry mass) of samples collected before (in situ, relaxed) or after the induction phase were analysed with LMMs and GLMMs, respectively, that included region, relaxation state (or induction treatment), and their interaction as fixed effects plus site and site x relaxation state (or site x induction treatment) as random intercepts. GLMMs of phlorotannin concentrations (proportions) were fit with a beta distribution and logit-link. 

Field Surveys – We conducted quadrat surveys of herbivore and fucoid assemblages in the late spring/early summer 2020 at 7 of the 9 collection sites (Fig.1a). Quadrats (0.5 x 0.5 m, n= 10 site-1) were haphazardly placed in the mid-intertidal zone at each site, avoiding tide pools and maintaining >1m between quadrats. For each quadrat, we estimated canopy % cover for each algal species using a 25-point-intercept grid and recorded identity and density (no. 0.25m-2) of fucoid holdfasts and herbivores. We measured the maximum frond lengths of fucoids located within all or part of the quadrat (up to 30 individuals quadrat-1). Densities of the two most common fucoids (F. vesiculosus and Ascophyllum nodosum) were analysed with negative binomial GLMMs (log-link function) that included the region, species, and their interaction as fixed effects, with site and the site x quadrat interaction as random intercepts. Frond lengths and herbivore densities (Littorina littorea and L. obtusata) were analysed with similar GLMMs. Percent cover data were logit-transformed for analysis with a similar LMM.  

Statistical Analyses – All statistical analyses were performed in R (R Core Team 2024). Models were fit via REML with type 3 contrasts to account for unequal sample sizes. Fixed effects were evaluated with marginal Wald c2 (GLMMs) F tests on Kenward-Roger df (LMMs). Diagnostic tests and residual plots were used to verify assumptions of homoscedasticity and normality for LMMs, with data transformations and/or weighted variance structures applied as necessary following Zuur et al. (2009). We evaluated significant interactions (p < 0.05) using linear contrasts on estimated marginal means and adjusted contrast p-values for the false discovery rate. Details of each model and analysis are given with corresponding results in the supporting information (Tables S1-6).