Biogeographic features mediate marine subsidies to island food webs

Abstract

Although marine subsidies often enrich terrestrial ecosystems, their influence is known to be context‐dependent. Additionally, the multitrophic impact of marine subsidies has not been traced through food webs across physically diverse islands. Here, we test predictions about how island characteristics can affect marine enrichment of food web constituents and how nutrients flow through island food webs. To evaluate enrichment and trace marine nutrients across food webs, we used stable isotopes of soil, flora, and fauna (n = 4752 samples) collected from 97 islands in British Columbia, Canada. Island area was the strongest predictor of enrichment across taxa; we found that samples were more ¹⁵N‐rich on smaller islands. Enrichment declined with distance from shore but less so on small islands, implying a higher per‐unit‐area subsidy effect. These area and distance‐to‐shore effects were taxon‐specific, and nearly twice as strong in basal food web groups. We also found that increases in δ¹⁵N correlated with increases in %N in basal trophic groups, as well as in songbirds, implying biologically relevant uptake of a potentially limiting nutrient. Path analysis demonstrated that subsidies in soil flow through plants and detritivores, and into upper‐level consumers. Our results reveal an interplay between island biogeography and marine subsidies in shaping island food webs through bottom‐up processes.

INTRODUCTION

MacArthur and Wilson's (1967) The Theory of Island Biogeography provides a framework for predicting equilibrium species richness according to island isolation and size, but marine inputs can impose deviation from these basic predictions (Barrett et al., 2003; Obrist et al., 2020; Polis et al., 1997). Marine subsidies connect terrestrial and aquatic ecosystems through the movement of materials (i.e., nutrients, detritus, and organisms), and alter the dynamics of recipient communities (Anderson et al., 2008; Polis et al., 1997). They can occur either through passive, donor-controlled processes, such as wind and waves depositing macroalgae on shores, or through active, animal-mediated means, including egestion and excretion by consumers, but also carcass deposition (Subalusky & Post, 2019). Since such cross-boundary transfers often involve the movement of limiting resources, they have the potential to facilitate increases in recipient habitat productivity. In fact, a meta-analysis across several biomes showed that subsidized ecological communities have 36%–57% larger populations (Montagano et al., 2019). This effect is amplified on smaller islands due to higher perimeter-to-area ratios, which facilitate greater per-unit-area effects of marine subsidies than on larger islands (Polis & Hurd, 1996). Accordingly, the Subsidized Theory of Island Biogeography (SIB) suggests more variable consumer responses on subsidized small islands due to trade-offs between increased resource availability and associated competition (Anderson & Wait, 2001). As a result, subsidized small islands are ideal model ecosystems to test the potential effects of marine subsidies on multilevel ecological communities.

Physical, spatial, and temporal characteristics of both subsidies and recipient habitats determine the magnitude and type of effect on recipient communities. Subsidies fluctuate in amount, quality, timing, and duration, eliciting variable effects on ecosystems (Piovia-Scott et al., 2019; Subalusky & Post, 2019). In one study, larger amounts of spawning salmon significantly increased riparian per-shrub berry production in Canadian streams (Siemens et al., 2020), which potentially affects avian and mammalian consumers (Traveset & Willson, 1998). Another study found that changing the seasonal timing of subsidy addition to plants resulted in a shift from bottom-up effects for herbivores to top-down suppression by predators (Gratton & Denno, 2003). Furthermore, physical characteristics of recipient habitats affect rates of delivery and retention of subsidies (Wiens et al., 1985). For instance, island beaches with less rocky substrate (i.e., more sand, gravel, and pebble substrate) have higher rates of wrack retention (Wickham et al., 2020). The degree of isolation between habitats also matters—lizard populations on small, isolated islands responded to artificial seaweed inputs more strongly than more connected populations in plots on large islands (Wright et al., 2020). A recent experiment evaluated the effects of both the characteristics of the subsidy and the recipient habitat on very small islands (4–157 m²). This study found varied predator effects on invertebrate prey depending on wrack subsidy amount, pulse frequency, and island size (Piovia-Scott et al., 2017, 2019), suggesting a connection between island characteristics and marine subsidies in island food webs. Overall, both the characteristics of the subsidy and those of the recipient habitat determine how a subsidy will influence the biotic components of an ecosystem, but these two components have rarely been studied in unison. Additionally, although much of the seminal work linking trophic ecology with marine subsidies occurred on islands (Anderson & Polis, 1998; Polis et al., 1997; Polis & Hurd, 1996; Stapp et al., 1999), these studies did not explicitly investigate the role of recipient island characteristics.

In this study, we use a suite of hierarchical models to evaluate how island characteristics (island area, mean slope, and shore-cast macroalgal biomass [i.e., wrack]) and distance from shore affect the enrichment of island taxa across six trophic levels. We use stable isotopes of carbon and nitrogen to represent enrichment by marine subsidies on 97 islands on the central coast of British Columbia, Canada, where we collected samples of soil, plants, herbivores (weevils), detritivores (isopods), carnivores (ground-dwelling beetles), and upper-level consumers (songbirds). The main conduits of marine inputs to these islands are wrack accumulation, sea spray, and river otter feces, urine, and discarded prey items. Evidence of enrichment by these sources can be tracked through elevated δ¹⁵N, δ¹³C, or both. We predicted island characteristics would have weaker effects on higher trophic levels because these organisms are more mobile (both within and among islands), and feed over larger, more heterogeneous areas. To determine whether subsidy inputs resulted in biologically relevant uptake of nitrogen, we also evaluated the associated increase in %N alongside enrichment of ¹⁵N.

We paired this analysis with a structural equation modeling approach to gain insight into the dominant pathways of subsidies through island food webs—specifically, to determine whether subsidies move through food webs through plant fertilization, direct consumption of subsidies by consumers, or both. Direct consumption is observable as elevated values of δ¹³C and δ¹⁵N in higher trophic levels but not terrestrial plants, while fertilization is observable as enrichment of ¹⁵N without simultaneous ¹³C enrichment in plants and subsequent consumers. These predictions rely on the reality that plants obtain nitrogen from the soil but carbon from the atmosphere. We predicted that subsidy effects would be primarily bottom-up: Subsidies enter the food web through the soil, and signals of marine influence would flow up the food web through plants, detritivore, herbivore, and carnivore invertebrates, and possibly even into upper-level consumers. We tailor this prediction by considering the types of subsidies considered here (i.e., wrack, feces, urine, and food scraps). Stronger upper-level consumer pathways likely exist in salmon-subsidized ecosystems (e.g., bears in Helfield & Naiman, 2006) or in ones with emerging aquatic invertebrates (e.g., birds, bats, and spiders in Recalde et al., 2020), where the subsidy is more attractive to upper-level consumers. Collectively, testing these predictions allowed us to assess how island characteristics can affect both the enrichment of individual components of food webs and the flow of nutrients through island food webs, and, as a result, tie together key theories in island biogeography and food web ecology.

METHODS Study area

In the summers of 2015, 2016, and 2017, we collected 4752 samples from soil and five taxonomic groups of plants and animals on islands (Figure 1), representing links in food chains across 97 islands in the Central Coast region of British Columbia, Canada (51°26′–52°3′ N and 127°41′–128°28′ W). See Appendix S1: Table S1 for species names of samples collected in both local Indigenous languages (Haíɫzaqvḷa and 'Uik̓ala). The ocean around these islands hosts widespread kelp forests and is one of the most productive marine environments in the world (Steneck et al., 2002). Surveyed islands (0.0001–3 km²) are part of the very wet, hyper-maritime subzone of the Coastal Western Hemlock biogeoclimatic zone that receives over 3 m of precipitation annually (Pojar et al., 1987). There are nearly 1500 islands in this region; the subsample we surveyed was representative of the biogeographical and geomorphological variation of the region and was selected using a two-step clustering method in SPSS statistical software (V23, IBM; results in Appendix S1: Table S1). For this cluster analysis, we considered five descriptors per island: distance from mainland, area, wave exposure according to the British Columbia ShoreZone dataset (Howes et al., 1994), normalized (size-independent) perimeter-to-area ratio, and percentage of area occupied by surrounding landmasses within 500 m of each island. Aside from exposure, these metrics were derived using 2-m resolution WorldView-2 satellite imagery.

FIGURE 1. Stable isotopes of δ13C and δ15N ± SE in soil, plants, and animals on islands on the central coast of British Columbia, Canada. Plots are separated into tertiles according to island size: small (0.00012–0.016 km2), medium (0.016–0.22 km2), and large (0.22–2.88 km2).

Field collections and measurements

To maximize the number of samples, we collected samples from a representative species or group of species from each trophic level that is common and abundant rather than species we explicitly know to be consumed by the next trophic level. Many of these samples were originally collected as discrete components of several separate projects under the umbrella of the “100 Islands Project” at the Hakai Institute (see Davidson et al., 2021; Nijland et al., 2017; Obrist et al., 2020; Wickham et al., 2019, 2020). Additionally, due to the nature of small islands, we were unable to collect samples of every trophic level on every island. This variation in the availability of organisms across study islands is part of the reason we chose not to fit stable isotope mixing models—instead, our modeling approach is inherently able to account for unequal sampling across islands. See Appendix S1: Table S3 for a full list of samples from each island. Additionally, there are different protocols for sample storage, cleaning, and preparation of different tissue types from different organisms. Although varying protocols could lead to biases in stable isotope composition, this is not an issue in our study because, unlike in a mixing model, we compare patterns in responses to environmental covariates and correlations with other taxa rather than raw differences in δ¹³C and δ¹⁵N. All samples within taxa were stored and prepared the same way (below).

Soil and plant samples

We collected soil (n = 973) and foliage samples of salal (Gaultheria shallon, n = 1485) and false lily-of-the-valley (Maianthemum dilatatum, n = 1166) on each study island. These two plant species were selected due to their ubiquity across islands and differences in nutrient regimes—salal is a species that is very tolerant of nitrogen-poor soils, whereas false lily-of-the-valley is a nitrophile (Hocking & Reynolds, 2011; Klinka et al., 1989). On each island, we collected samples immediately adjacent to the shoreline at the furthest point in each cardinal direction and 40 m inland from this point. On islands smaller than 80 m across, we sampled in the approximate middle of the island. We also used soil and vegetation samples collected at each of the 301 random point count locations on the same set of islands, described in Obrist et al. (2020). The distance to the shore of these samples ranged from shoreline to approximately 400 m inland. We ground dried leaf matter (stored in paper envelopes, then dried at 60°C for a minimum of 48 h) using a Wig-L-Bug (Dentsply/Rinn Corp., York, PA, USA) to prepare for stable isotope analysis.

Invertebrate samples

We collected samples of common species of weevils (Coleoptera: Curculionidae, n = 137), isopods (Isopoda: Porcellionidae, Ligiidae, n = 142), and ground beetles (Coleoptera: Carabidae, n = 271) to represent invertebrates from the herbivore, decomposer, and carnivore feeding guilds, respectively. We placed networks of pitfall traps at four nearshore and one interior location on each island. Nearshore locations corresponded to the soil and plant sampling sites. Our traps were filled with 1.5 cm of ~30% propylene glycol solution diluted with seawater, to which surfactant was added. Propylene glycol is nontoxic and appropriately preserves specimens for isotope work (Nakamura et al., 2020). We removed invertebrates from the traps after 3–4 days, rinsed specimens with fresh water, and stored them in a freezer in 95% ethanol until they could be processed. In the laboratory, we identified, cleaned, dried, and ground specimens into a uniform sample using a Wig-L-Bug to prepare for stable isotope analysis. Specimens were selected for each guild based on availability, in the following descending order of preference: Steremnius carinatus was preferred, but if not available, then S. tuberosus for weevils; Porcellio scaber, then Ligidium gracile for isopods; and Pterostichus algidus, P. amethystinus, P. crenicollis, P. lama, and then Scaphinotus angusticollis for ground beetles. Although individual ground beetles were large enough, to ensure adequate material for analysis, we had to combine up to five individuals from each trapping location for isopods and weevils.

Bird samples

To represent upper-level consumers, we used conspecific playback to catch live birds of five common species in mist nets: chestnut-backed chickadees (Poecile rufescens), song sparrows (Melospiza melodia), fox sparrows (Passerella iliaca), Pacific wrens (Troglodytes pacificus), and orange-crowned warblers (Leiothlypis celata). We took both a fecal (n = 241) and feather sample (second right rectrix, n = 337) from each bird. As all five species are primarily insectivorous here on their breeding grounds (Davidson et al., 2015), we considered all of them in combination as representative upper-level consumers in this ecosystem. To prepare for stable isotope analysis, we first soaked feathers for 24 h in a 2:1 solution of chloroform and methanol, dried them at 60°C for 48 h, and used scissors to cut them into pieces to pack into capsules. Feces were dried and packed directly into capsules.

Stable isotope analysis

The main pathways by which marine subsidies enter terrestrial food webs—direct consumption and plant fertilization—can be evaluated using stable isotopes of carbon and nitrogen (δ¹³C and δ¹⁵N). Elevated δ¹⁵N and δ¹³C values in samples imply higher ratios of the heavy isotope, ¹⁵N to ¹⁴N and ¹³C to ¹²C. These samples can be said to be enriched in ¹⁵N and ¹³C. These isotopes can be used to estimate trophic positions of organisms and marine inputs in organisms' diets. Nitrogen enrichment can be tracked through food webs; ¹⁵N fractionates by ~3.4‰ per trophic level (Post, 2002) because metabolic processes favor the lighter isotope of nitrogen, ¹⁴N, resulting in ¹⁵N-enriched animal tissues. Nitrogen enrichment can also indicate marine input, as marine primary producers typically have higher δ¹⁵N values than terrestrial primary producers—for example, in our system, the mean δ¹⁵N value of seaweed wrack was 7.3 ± 1.0 (mean ± SE; Appendix S1: Figure S1), while salal and false lily-of-the-valley had means of −0.3 ± 5.6 and −2.1 ± 6.4, respectively. Conversely, carbon isotopes do not fractionate much (~0.4‰) when passed up the food chain. Accordingly, they better reflect the δ¹³C values of the basal primary producer (Post, 2002). Due to different carbon sources in terrestrial and marine ecosystems, marine photosynthesizers also have higher values of δ¹³C than terrestrial plants (Fry, 2006). In our system, the mean δ¹³C of wrack was −15.4 ± 3.0, salal was −31.4 ± 1.7, and false lily-of-the-valley was −30.4 ± 2.6. Consequently, higher δ¹³C values in a consumer's tissue imply a ¹³C-enriched diet, such as that from a marine origin.

We processed all samples and quantified δ¹⁵N, δ¹³C, and %N at the Pacific Forestry Centre in Saanich, British Columbia, using a Flash 2000 elemental analyzer coupled to a ConFlo IV interface and DELTA Advantage isotope ratio mass spectrometer (Thermo Fisher Scientific, Waltham, MA).

Island characteristics and wrack biomass

We derived island area using 2-m resolution WorldView-2 satellite imagery. The area was measured as land within the outer vegetation boundary excluding the tidal and wash (shore) zones. Barren or nonvegetated areas inside the vegetation boundary were included. The mean island slope metric, however, represents the entire island, including the shore zone. The slope was calculated based on LiDAR where available; where LiDAR was not available, we flew an unmanned aerial vehicle and calculated a ground model with 50-cm resolution. We measured distance to shore as the shortest distance from each sample collection site to the edge of the vegetated land (vegetation boundary) on each island.

We weighed wrack at four locations (one in each cardinal direction) on each island. We placed two 20-m transects parallel to the water—one transect at the most recent high tide line and one at the most recent storm line—and randomly placed three 1-m² quadrats along each (Wickham et al., 2020). We identified each species in the quadrat to calibrate wet weights to dry weights for total wrack biomass per quadrat (Wickham et al., 2019). We calculated the mean wrack biomass of the 24 quadrats measured (4 sides × 2 transects × 3 quadrats) on each island.

Data analysis Enrichment analyses: δ¹⁵N and δ¹³C

To determine how island characteristics can affect the enrichment of island organisms, we first fit separate global linear mixed-effects models (LMMs) to the ratios of heavy nitrogen and carbon isotopes, ¹⁵N and ¹³C (measured as δ¹⁵N and δ¹³C), of samples using the glmmTMB package in R v. 4.0.3 (Brooks et al., 2017; R Core Team, 2021). We fit a separate global model for each taxon, each including island area, distance to shore from which a sample was collected, wrack biomass, the mean slope of the island including the shore zone, and an interaction term between area and distance to shore. Each model also included island as a random effect to account for unequal sampling and because some predictors are on the scale of the sample (e.g., distance to shore), while others are island-specific (e.g., area). We log₁₀-transformed island area, log_e-transformed distance to shore, and square-root-transformed wrack biomass to best linearize their relationships with isotopes prior to scaling and centering all independent variables. Since we had no a priori predictions about the relative importance of our independent variables, we model-averaged to calculate coefficient estimates across all possible subsets of our four predictors and the interaction term using the MuMIn package in R (Barton, 2020). Our main analyses used fecal samples from songbirds to represent a “snapshot” of what was recently consumed by birds (Appendix S1: Tables S3 [δ¹⁵N] and S4 [δ¹³C]), but we repeated all analyses with feather samples (Appendix S1: Tables S5 [δ¹⁵N] and S6 [δ¹³C]). We checked for multicollinearity between predictors in each model using variance inflation factors and checked all model diagnostics using the DHARMa package in R (Hartig, 2020).

Subsidy analysis: %N

To determine whether the enriched samples received a biologically relevant nutrient subsidy (i.e., that samples with more marine-derived nitrogen also contained more nitrogen overall), we fit a LMM to evaluate the effects of δ¹⁵N on %N while accounting for island area, distance to shore, wrack biomass, mean island slope, and the random effect of island. We again model-averaged across all possible subsets of our predictors to obtain coefficient estimates. Here, we used bird feather samples in our main analyses rather than fecal samples. Birds excrete and egest waste through a single cloaca, so fecal samples' %N reflects the high variability of excreted nitrogenous wastes in birds, making it less useful for detecting patterns (Vanderklift & Ponsard, 2003). The results from the averaged model were comparable to those from the global model. Therefore, to determine whether the inclusion of δ¹⁵N improved the fit of the model, we also used the Akaike information criterion corrected for small sample sizes to compete for the global model for each taxon with a null model (one without δ¹⁵N and solely island characteristics) (Appendix S1: Table S7).

Path analysis

To determine the mechanism by which subsidies are incorporated into island food webs, we used the lavaan package in R (Rosseel, 2012) to fit a structural equation model (SEM) of pathways through food webs for each stable isotope tracer (δ¹⁵N and δ¹³C). This analysis used island-level means of stable isotopes for each taxon. This form of path analysis allowed us to evaluate the strength of pathways among components of island food webs, leading through island ecosystems from the soil to upper-level consumers in a holistic, ecosystem-level framework. The SEMs allow for testing of strengths of hypothesized causal relationships among variables, while holding all other variables in the model constant (Pearl, 2012). As such, a causal relationship is implied to exist if the coefficient estimate's 95% CI does not overlap with zero. The SEMs also allow testing of the strengths of both direct and indirect relationships, where indirect paths are quantified by multiplying direct path coefficients. Since we were unable to collect samples of all organisms on all islands, and the fitting of SEMs requires complete datasets, we used the mice package (van Buuren & Groothuis-Oudshoorn, 2011) in conjunction with lavaan to impute missing data using predictive mean matching. Due to their potential to have island-level impacts on stable isotope signatures, we included island area and wrack biomass as predictors to impute the missing data. We chose to include 100 imputations as a conservative approach to the rule of thumb that the number of imputations should be at least equal to the percentage of incomplete cases (White et al., 2011). In our case, we were missing at least one sample from 83% of islands (16/94 complete cases), meaning we needed to impute at least 83 iterations. However, even with this approach, some of the estimates of fraction of missing information, particularly those for upper-level consumers, are higher than 0.5, meaning they should be interpreted with a degree of uncertainty. We used the fecal samples taken from birds (Appendix S1: Tables S8 and S9 for δ¹⁵N, δ¹³C results) but also repeated the analysis using feather samples (Appendix S1: Tables S10 and S11).

We used an iterative decision-making process to build our SEM to trace δ¹⁵N through the food web. We started with a biologically accurate global model of pathways from marine subsidies to upper-level consumers. At each iteration, if the model did not meet the conventional cut-offs (below), we used the modification indices to reselect the subset of pathways that needed to be included to improve fit but maintain biological accuracy. Specifically, a model was required to meet at least three of the following four criteria: a χ² value of at most twice the number of df, a comparative fit index of >0.95, alongside a standardized root mean square residual, and root mean square error of approximation of <0.08 (Hu & Bentler, 1999). Because the model initially had difficulty differentiating between the two plant species, we averaged across all plant samples on an island to represent terrestrial primary producers more generally. Detailed documentation of our entire SEM and decision-making process is available online (www.github.com/debobrist/food-web-ibt). Once the biologically accurate path selection met the standard requirements for model fit, we used the same structure to fit the SEM for δ¹³C, except that the regression models between soil and plants were excluded because we did not expect soil δ¹³C to predict plant δ¹³C.

RESULTS Enrichment analysis: δ¹⁵N

Island characteristics and marine subsidies tended to affect ¹⁵N enrichment across taxa in similar ways, with variation in the strength of the effect depending on the trophic position of the recipient (Figures 2 and 3). Island area was the strongest predictor of δ¹⁵N across soil, island flora, and fauna (Figures 2a and 3a; Relative Variable Importance [RVI] = 1.00 for all taxa [Appendix S1: Table S6]). As predicted, all groups were more enriched in ¹⁵N on smaller islands, an effect that was approximately twice as strong in the two plant species and the herbivorous weevils as in the soil, in the detritivores, and in the upper-level consumers (carnivorous ground-dwelling beetles and songbirds). The second strongest effect was the distance from shore at which a sample was collected; δ¹⁵N decreased with increasing distance from shore (Figures 2a and 3b; RVI = 1.00 for all groups except songbirds). The effect of distance to shore, however, depended on island size in all taxa except for songbirds; enrichment in ¹⁵N decreased more quickly with distance to shore on larger islands than on smaller ones, particularly for salal (Figure 2a). Although all groups were enriched closer to the land–sea interface, this effect was at least twice as strong in the soil and in the two plant species (particularly in false lily-of-the-valley) than in the invertebrate taxa. Mean island slope and wrack biomass had negligible effects on ¹⁵N enrichment of most taxa, aside from the two plant species, in which δ¹⁵N was slightly higher on islands with higher wrack biomass (Figure 2a; Appendix S1: Figure S2a), and detritivorous isopods and carnivorous beetles, where ¹⁵N was less enriched on islands with steeper topographical slopes (Figure 2a; Appendix S1: Figure S2b).

FIGURE 2. Standardized coefficient estimates from taxon-specific model-averaged linear mixed-effects models for (a) δ15N and (b) δ13C. Bars represent 95% CI around estimates.

FIGURE 3. Effects of island area (in square meters) and distance to shore (in meters) on δ15N (a, b) and δ13C (c, d) values of island soil, plants, invertebrates, and songbirds. Shaded areas represent 95% CI around taxon-specific model-predicted means.

Enrichment analysis: δ¹³C

We detected considerable variability in ¹³C enrichment across taxa in response to island characteristics and marine subsidies. Interestingly, δ¹³C in the soil was higher on larger islands, while salal had slightly higher δ¹³C on smaller islands (Figures 2b and 3c; Appendix S1: Table S7). Both plant species and isopod detritivores collected closer to the shoreline had higher δ¹³C (Figures 2b and 3d). The effect of distance to shore was approximately twice as strong in false lily-of-the-valley as in salal. For soil, δ¹³C increased with distance to shore on smaller islands but decreased with distance to shore on larger islands. Contrary to δ¹⁵N, for δ¹³C, the effect of distance to shore was stronger in plants on smaller islands. Both plant species, carnivorous beetles, and songbirds also had lower δ¹³C on islands with steeper slopes (Figure 2b).

Subsidy analysis: %N

For most taxa, %N was not related to island characteristics, but more basal taxa experienced elevated %N in tandem with higher δ¹⁵N. Specifically, %N was higher in soil and plant samples that had higher δ¹⁵N (Figure 4a,c). This effect was approximately three times stronger in false lily-of-the-valley than in salal and in the soil. Percent nitrogen in songbird feathers was also slightly higher in samples that had higher δ¹⁵N (Figure 4a,b). Salal %N was slightly higher on larger islands and further away from the shoreline, whereas %N in the soil was slightly lower on larger islands and higher closer to shore (Figure 4a,c).

FIGURE 4. (a) Standardized coefficient estimates from taxon-specific model-averaged linear mixed-effects models for %N. Bars (a) and shaded areas (b, c) represent 95% CI around estimates. Modeled relationships between δ15N and %N in (b) songbirds, invertebrate beetles, isopods, and weevils, and in false lily-of-the-valley, salal, and soil (c).

Path analysis

Our SEM revealed evidence of the plant fertilization pathway as the dominant mechanism of δ¹⁵N enrichment through island food webs (Figure 5a), while invertebrate detritivore δ¹³C appeared to be driven by plant δ¹³C (Figure 5b). We found evidence that subsidies in the soil in the form of δ¹⁵N strongly predicted enrichment of ¹⁵N in plants (0.93 ± 0.10 [path coefficient estimate ± 95% CI]). This nutrient subsidy was subsequently reflected in isopod detritivores, whose δ¹⁵N signatures were predicted by ¹⁵N enrichment in plants (0.38 ± 0.27) and soil (0.50 ± 0.27). Furthermore, the level of δ¹⁵N in carnivorous beetles was predicted by that of their potential prey source, isopod detritivores (0.95 ± 0.33). We also found strong correlations between the δ¹⁵N of the soil and that of herbivorous weevils (0.83 ± 0.49), but that δ¹⁵N of weevils was not predicted by δ¹⁵N of plants (0.09 ± 0.35). When tracing the stable isotope of carbon, we found evidence that plant δ¹³C predicted δ¹³C in isopod detritivores (0.46 ± 0.25), which subsequently predicted δ¹³C in bird fecal samples (0.62 ± 0.35). We also found a negative correlation between soil δ¹³C and weevil δ¹³C (−0.39 ± 0.33).

FIGURE 5. Structural equation models showing pathways of (a) δ15N and (b) δ13C enrichment through soil, and island plants and animals. Highlighted (darker) paths with associated path coefficients are statistically significant at the 95% level. PhyloPic: Weevil: JCGiron (CC BY 3.0); isopod: B. Lang (Public Domain Mark 1.0); beetle: T. M. Keesey, T. Assmann, J. Buse, C. Drees, A. L. L. Friedman, T. Levanony, A. Matern, A. Timm, and D. W. Wrase (CC BY 3.0); and bird: Troglodytes troglodytes, A. Caravaggi (CC BY-NC-SA 3.0).

DISCUSSION

Our field study explicitly integrated biogeographical characteristics of islands with food web ecology at a large scale of variation in island size, island slope, and marine subsidy. As predicted, we found that ¹⁵N enrichment of soil, primary producers, invertebrate herbivores, detritivores, carnivores, and songbirds was stronger on smaller islands and closer to shore and that these patterns were stronger in lower trophic levels. We found that distance to shore was less important on smaller islands, suggesting higher per-unit-area effects of subsidies and thereby lending support to SIB. Additionally, we found that increased %N corresponded to increased δ¹⁵N in lower trophic levels, implying biologically relevant uptake of a potentially limiting nutrient. Finally, using path analysis, we found evidence that marine subsidies predominantly enrich temperate island food webs via fertilization (rather than direct subsidy consumption).

We showed that island characteristics mediate enrichment by marine subsidies across all six trophic levels of island food webs in our analyses. Island area had the largest effect; ¹⁵N enrichment was higher on smaller islands across all taxa. The strength of this effect depended on the trophic position of the organism—stronger relationships occurred between ¹⁵N enrichment and island area at lower trophic levels. A likely explanation is that lower trophic-level organisms (e.g., plants) uptake nutrients over smaller areas, and are therefore more affected by characteristics of their local environment, while higher trophic-level taxa (e.g., songbirds) forage over larger areas, likely spanning multiple locations both within and among islands. Indeed, a similar study found riparian soil and herbivorous weevils to have higher levels of δ¹⁵N closer to salmon-bearing streams, but δ¹⁵N of more mobile, carnivorous beetles was unaffected (Rammell et al., 2021). We also found that ¹⁵N enrichment drops off more rapidly with increasing distance to shore on larger islands. This effect is likely explained by larger perimeter-to-area ratios on smaller islands, making them more permeable to marine-derived nutrients further inland. To make sense of this relationship, consider an island that is 10 m in diameter. An organism 5 m from shore is 5 m from shore in all directions, whereas on an island that is 100 m in diameter, that organism is 95 m from shore in the opposing direction. This is a key dimension of SIB, which grounds itself in the foundational prediction that smaller islands experience a higher per-unit-area influence of marine inputs (Anderson & Wait, 2001).

By evaluating the relationship between marine-derived nitrogen and %N across constituents of a coastal temperate island food web, we were able to confirm that marine inputs provide a biologically relevant nutrient subsidy. Specifically, we found that the soil, foliage from two plant species, and songbird feathers (but not feces) experience higher %N in tandem with higher δ¹⁵N. This finding implies potential for higher subsidy-fueled productivity. Nitrogen is often the most important growth-limiting soil nutrient, and nitrogen limitation is common in these rainforest ecosystems (Klinka et al., 1989). Elemental composition (i.e., ratios of different nutrients) can vary more in plants than in animals due to plant cells' large central vacuoles, which can store excess nutrients not immediately used for metabolic processes (Sterner & Elser, 2002). Storage of a nutrient beyond what is required for growth is termed “luxury consumption” and often indicates limitation by a different nutrient. Indeed, luxury consumption of nitrogen often indicates phosphorus limitation (Greenwood, 1976). Since egestion, excretion, and scent marking by river otters are all processes driving increases in available nitrogen (Ben-David et al., 1998), it is possible that islands subsidized by river otter fertilization experience a shift in nutrient limitation. This hypothesis is supported by our observation that even salal, which is ericaceous and generally not nitrogen-limited (Bennett et al., 2003), experienced an increase in %N alongside δ¹⁵N. Ericaceous plants have specialized mycorrhizae, which facilitate the acquisition of nutrients from rotting organic matter in nutrient-poor conditions to overcome nitrogen limitation. Similar increases in %N alongside δ¹⁵N have been seen in ericaceous false azalea (Menziesia ferruginea) on salmon (Oncorhynchus sp.) streams in this region of British Columbia (Hocking & Reynolds, 2011). Furthermore, the effect of δ¹⁵N on foliar %N was three times stronger in nitrophilic false lily-of-the-valley, likely because nitrophiles' mycorrhizal associations allow them to more readily capitalize on available nitrogen in the soil's litter layer (Read, 1991).

In our study, aside from songbirds, higher trophic levels showed enrichment of ¹⁵N without higher levels of %N, likely because the ratios between carbon, nitrogen, and phosphorus are far less flexible in animal tissues than in plant tissues (Sterner & Elser, 2002). However, although subsidy effects are not apparent at the tissue level, subsidies elicit variable consumer responses in individuals that can scale up to community-level consequences (Yang et al., 2010). For instance, herbivorous invertebrates tend to have higher rates of larval survival and development in association with higher %N content in foliage (Mattson, 1980). Such enhanced nutrition can lead to higher prey densities for upper-level consumers, even to the point where subsidized herbivores overgraze plant populations (Leroux & Loreau, 2008; Piovia-Scott et al., 2019). Indeed, higher N content in foliage has been associated with multiple pathways of food web change; subsidized plots hosted higher abundance of herbivorous arthropods and predatory lizards in synchrony with higher levels of foliage growth and leaf damage on islands in the Bahamas (Spiller et al., 2010). In addition, food web stability is highly dependent on the trophic position of the subsidy recipient (Leroux & Loreau, 2008; Polis et al., 1997). Low amounts of inputs tend to support food web stability, whereas large, variable pulses of inputs can have a destabilizing effect through stronger trophic cascades (Huxel & McCann, 1998; Leroux & Loreau, 2008; McCary et al., 2021; Recalde et al., 2020). Understanding the pathways by which subsidies enter food webs and knowing who benefits from them can aid our understanding of individual-level to community-level consumer responses.

Finally, we determined that the dominant mechanism of marine nutrient propagation through island food webs was through the plant fertilization pathway. Specifically, our path analysis provided evidence that increased soil δ¹⁵N leads to higher plant δ¹⁵N and that this pattern extends up the food web into upper-level consumers. This hypothesis was also supported by our enrichment analyses, where, on smaller islands and closer to shore, we found enrichment of consumer ¹⁵N without simultaneous enrichment of ¹³C. Enrichment of consumer ¹⁵N but not ¹³C implies that the subsidy enriches soil, litter, and vegetation without direct consumption (Hocking & Reimchen, 2002). Similar conclusions were reached by a study on mouse diet variation on islands in our study area, where indirect effects of marine subsidies on mouse diets were observed through an increase δ¹⁵N of fecal pellets without a simultaneous increase in δ¹³C (Davidson et al., 2021). Due to the lack of a consistent pattern, we infer that variability in δ¹³C across taxa is likely driven by variability in plant δ¹³C. Plants obtain carbon from the atmosphere, but δ¹³C varies according to soil moisture, salinity, and nitrogen availability (Dawson et al., 2002). Additionally, although we found evidence of bottom-up fertilization effects in this ecosystem, some uncertainty remains about the relative contributions from different subsidy sources. Although δ¹⁵N in plants was higher on islands receiving more wrack, wrack biomass appeared to have no effect on δ¹⁵N in the soil nor in any of the consumers. As such, we infer that, moving up the food web, the signal of animal-mediated nutrient deposits from river otters likely overwhelms the effects of wrack fertilization.

CONCLUSIONS

Considering how recipient habitat characteristics affect the enrichment of multitrophic island food webs is an important step to improve our understanding of how subsidies affect dynamic island ecosystems more generally. Our finding that organisms are enriched on smaller islands with higher per-unit-area effects of subsidies lends strong support to rarely tested theoretical predictions of SIB (Anderson & Wait, 2001). Additionally, we determined that plant fertilization was the dominant pathway by which marine subsidies enter temperate island food webs—an effect that we detected propagating up through soil, plants, and invertebrates, and into upper-level consumers. Understanding the pathways by which subsidies enter food webs and knowing who benefits from them are fundamental to our understanding of individual-, population-, and community-level consumer responses.

AUTHOR CONTRIBUTIONS

Debora S. Obrist, John D. Reynolds, Brian M. Starzomski, and Chris T. Darimont conceived of the study; Debora S. Obrist, Christopher M. Ernst, Owen T. Fitzpatrick, Sara B. Wickham, Jeremiah C. Kennedy, Wiebe Nijland, and Luba Y. Reshitnyk collected the data; Debora S. Obrist, Patrick J. Hanly, and Norah E. M. Brown conducted statistical analyses; Debora S. Obrist wrote the first draft; and all authors contributed substantially to revisions.

ACKNOWLEDGMENTS

We are grateful to the Haíɫzaqv and Wuikinuxv First Nations for their support in conducting this research. We also thank all the 100 island field and laboratory technicians for their hard work on this project, the Hakai Institute for field and laboratory support, and Carl Humchitt for his guidance and knowledge. In addition, we would like to thank Jane Pendray for guidance with field season planning, Allison Dennert for sharing plant expertise, and several other members of the Reynolds Lab and the Earth to Ocean Research Group for statistical advice. This work was supported by the Tula Foundation at the Hakai Institute, MITACS, and NSERC Discovery Grants to John D. Reynolds, Brian M. Starzomski, and Chris T. Darimont.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

DATA AVAILABILITY STATEMENT

Data (Obrist et al., 2021) are available from the Hakai Institute's public data repository: https://doi.org/10.21966/8j4j-7039.