Many management policies consider the conservation of one group of species independent of their interactions with other species of conservation concern. However, populations are frequently linked in natural systems (e.g., predators and prey; Stenseth et al., 1997) and understanding how changes in one taxon affect another remains a key challenge in conservation and ecology. We may expect interactions between taxa to be strong in aquatic systems where plants and seaweeds act as foundation species (sensu Dayton, 1973) for a diversity of organisms that utilize macrophytes for habitat and/or as a primary food source (Bakker et al., 2016). Within these systems, strong linkages between macrophytes and their consumers can arise when consumers depend on the macrophytes, but also exert strong top-down control on macrophyte populations through extensive tissue damage and biomass removal (reviewed in Poore et al., 2012; Wood et al., 2012, 2017; Bakker et al., 2016). In cases where both the primary producer and consumer are a focus of conservation, recognition of these reciprocal linkages can inform the effective management of both taxa.

In temperate systems, scientists and managers use coastal waterfowl to motivate the protection of seagrass beds, given an understanding that seagrass provides high-value resources which supports migratory and breeding activities (Sedinger et al., 2011; Schamber et al., 2012). Climate change, over-hunting, and habitat destruction all threaten waterfowl populations (e.g., Ward et al., 2005). Declines in Zostera spp. (hereafter Zostera) can reduce waterfowl carrying capacities or necessitate a behavioral shift to other sites or resources. In the 1930s, an epidemic of seagrass wasting disease (Labyrinthula zosterae) removed the majority of Zostera marina from North Atlantic coastlines in America and Europe, followed by a collapse in brant (Branta bernicla sp.) populations (Addy and Aylward, 1944). Although brant were thought to be obligate consumers of Z. marina (reviewed by Muehlstein, 1989), the geese have since recovered with a shift in diet to agricultural lands (Moore et al., 2004; Inger et al., 2006a,b). The Pacific brant (Branta bernicla nigricans) largely remains an obligate Z. marina consumer during migration and over winter, but will occasionally consume other marine primary producers such as Ulva (Hereu and Jorgensen, personal observation) and feed in Arctic salt marshes while nesting in the summer (Ward et al., 2005). Local losses of Zostera may be compensated through behavioral flexibility, which enables brant to use other stopover or wintering locales (Sedinger et al., 2006) rather than necessitating dietary flexibility or resulting in population-level declines. Nevertheless, the local consequences of Zostera losses appear dramatic, as predicted by models relating brant carrying capacity to Zostera abundance (Stillman et al., 2015) and evidenced by the decline in brant geese populations in Morro Bay after Z. marina essentially disappeared in the mid-2000s (Harencar et al., personal communication).

While conservation of avian herbivores includes consideration of their seagrass food, the linkage between seagrasses and their herbivores is typically not considered in the management of temperate seagrass populations (Maxwell et al., 2016). Several waterfowl species consume Zostera (Figure 1), but their top-down effects on seagrass bed dynamics is often assumed to be minor due to the birds’ temporally limited residence and rapid productivity in Zostera (Valentine and Heck, 1999; Ganter, 2000; Valentine and Duffy, 2007). Studies of top-down effects by grazing birds on seagrasses typically aim to estimate the impact of herbivory as reductions of standing stock biomass or production (Thayer et al., 1984; Bakker et al., 2016) and results show that annual Zostera productivity can be two or three times greater than peak standing stock biomass (McRoy, 1966, 1970; Penhale, 1977; Thom and Hallum, 1990; Kaldy, 2006). Furthermore, grazing pressure is neither uniform temporally nor spatially, and total impacts of waterfowl on annual production can be minimal (as low as 3%; Thayer et al., 1984). Thus, the perspective first proposed by Kikuchi and Peres (1977) that most vegetative material in temperate seagrass beds is consumed through detrital pathways remains widely held.

Despite this prevailing view of minor top-down impacts of waterfowl on temperate seagrasses, reductions in plant standing crop by waterfowl could be significant where waterfowl biomass per unit area is high, where alternative foraging areas are scarce, or when waterfowl mobility between sites is limited (Thayer et al., 1984; Wood et al., 2012, 2017). Several case studies highlight that birds have contributed to large-scale losses of seagrass beds. Many seagrass-consuming birds dig up and eat the rhizomes or shoot meristems, which can damage the bed and inhibit recovery (Thayer et al., 1984). In San Francisco Bay, California, migrating Canada geese (Branta canadensis) remove all the eelgrass shoots in a shallow bed upon their arrival in the fall. Experiments showed that unless the geese are excluded, the bed must recover by seedling recruitment each spring (Kiriakopolos, 2013). Reliance on seeds is risky; Canada geese grazing led to local extinction of eelgrass in a New England bed where seedling recruitment was minimal (Rivers and Short, 2007). In New Zealand, intense grazing by black swans (Cygnus atratus) can remove 19–20% of the annual Zostera production, resulting in a 43–69% decrease in the standing biomass during the following growing season (Dos Santos et al., 2012, 2013). These studies highlight the potential for waterfowl to exert strong top-down pressure on seagrass ecosystems, but the generality of such phenomena is unknown and there remains little systematic understanding of the conditions under which bird herbivory alters the structure, function, or dynamics of temperate eelgrass beds.

Here we address the assumption that linkages between vertebrate herbivores and temperate seagrasses are weak by reviewing the literature to explicitly quantify the strength of waterfowl interactions with seagrass of the genus Zostera at temperate latitudes (30–60° north and south of the equator). Specifically, we examine the prevalence, scope, impact, and consequences of waterfowl herbivory on Zostera by addressing the following research questions:

We conclude our analysis with a broader view of system-wide consequences of waterfowl herbivory on seagrass ecosystems, which extend beyond reductions in seagrass abundance. Finally, we consider the implications of these reciprocal linkages for the conservation of both waterfowl and seagrass species as well as identify important avenues for future research in this area.

We used a systematic literature review in combination with meta-analytic techniques to (1) quantitatively describe the correlation between Zostera abundance and waterfowl abundance, (2) measure the effect size of waterfowl herbivory on Zostera abundance, and (3) compile a list of avian species known to consume Zostera. On November 15, 2016 we searched ISI Web of Science using the terms “Zostera AND (bird^∗ AND (herbiv^∗ OR graz^∗) OR (waterfowl OR goose OR geese OR brant).” To uncover additional relevant literature, we applied a forward and backward search from those citations, consulted earlier reviews (Valentine and Heck, 1999; Olsen, 2015; Wood et al., 2017), and contacted researchers with coastal waterfowl expertise (P. Clausen, Denmark; M. Hori, Japan; A. Olsen, United States). We also included an additional unpublished data set (C. M. Hereu and P. Jorgensen, experiment performed November 2011 to November 2012). For all research objectives, we excluded studies that did not assess natural grazing by waterfowl (e.g., researcher-simulated herbivory by clipping leaves or digging pits).

The analysis of the shared variation in waterfowl and Zostera required at least four time points or locations. Studies that met this criterion ranged widely in the scope of observation; some spatial studies covered bird use of a single tidal flat with heterogeneous Zostera distribution (e.g., Percival and Evans, 1997; Dixon, 2009), whereas others compared tidal flats spanning hundreds or thousands of kilometers (e.g., Clausen and Percival, 1998; Moore et al., 2004). Most temporal data sets had annual time steps, in which birds and Zostera were recorded using similar methods over as much as three decades, but we also included one study with monthly time steps, where other evidence indicated that changes in bird abundance were not due to seasonal migration but rather to behavioral choices among a range of possible habitats (Tinkler et al., 2009). Variability in bird use due to tidal cycle changes in water level is well established (Clausen, 2000) but at too fine of a temporal scale for relevance in this meta-analysis. At large temporal or spatial scales, birds were typically assessed for bird-days (a metric of population-level habitat use at a site) integrated over a season or for peak numbers. When data sets were of smaller temporal or spatial scales, birds were typically tracked for behavior, such as foraging time or fecal dropping rate. Large-scale data sets for Zostera included estimates of area covered (e.g., at the scale of km²), whereas small-scale data sets for Zostera had sample units of percent cover or biomass per area (usually per m²). We calculated effect size using Pearson’s correlation coefficient (r) between waterfowl and Zostera, followed by transformation with Fisher’s Z (Z =0.5*In(1+r1−r)) and a calculated standard deviation derived from sample size (VZ =1n−3). We calculated an effect size for each waterfowl species separately if the study reported multiple species (Percival et al., 1998; Petersen et al., 2008; Balsby et al., 2017).

We applied a linear mixed effects model to the transformed correlation coefficients, with study considered a random effect to account for multiple bird species measured across the same Zostera samples. We assessed the statistical importance of two potential predictor variables: first, the amount of variation in the abundance of Zostera as calculated by the coefficient of variation among samples within a study (standard deviation divided by the mean; CV); and second, whether the study was carried out within a bay or year (defined as a small scale) or across bays or years (large scale). We tested how well the inclusion of the predictor variables explained the heterogeneity in the data using a Cochran’s Q-test, which tests the null hypothesis of homogeneity among samples. Analysis was carried out using R software (R Core Team, 2015) with the function rma.mv in the metafor package (Viechtbauer, 2010).

To examine the top-down effects of waterfowl on Zostera, we only included studies that measured the response of Zostera abundance variables to the presence and absence of waterfowl (i.e., caging experiments, or observational experiments where waterfowl presence and absence varied spatially or temporally). Response variables included metrics that reflect the abundance of Zostera at the plot level (e.g., number of shoots in a 1-m² area, canopy height, percent coverage of vegetation, and aboveground and belowground biomass). Within a study, we extracted data points across multiple response variables and multiple data points. We treated publications that reported results from multiple study sites as independent only if the seagrass beds were non-overlapping (Dixon, 2009). We excluded records if the variance could not be extracted with the information provided in the publication (Madsen, 1988; Nacken and Reise, 2000, and certain data points within Tubbs and Tubbs, 1983). A subset of data points within one study (n = 14 data points in Kiriakopolos, 2013) had a mean and standard deviation value equal to 0 in both the presence and absence of waterfowl, so these points were removed. Effect sizes were calculated as the standardized mean difference (Hedges’d,d =M1−M2SD*pooled,SD*pooled =(n1−1)SD12+(n2−1)SD22n1+n2−2, Hedges and Olkin, 1985; with confidence intervals corrected for heteroscedastic population variances between groups, Bonett, 2009) using the function escalc from the package metafor (Viechtbauer, 2010) in R (R Core Team, 2017). We only retained the maximum effect size from each study for further analysis. We excluded one outlier from Rivers and Short (2007) because including the maximum effect size in regression analyses (see below) yielded standardized residuals greater than two standard deviations away from the mean. Instead, we retained the second largest maximum effect size for this study.

We first used a random-effects model to estimate the overall effect size of waterfowl on Zostera abundance and test for heterogeneity in effect sizes among studies. Analyses used the rma function in R metafor package (Viechtbauer, 2010). To test for potential publication bias among our set of studies, we used a funnel plot to visualize effect size versus standard error and statistically tested for asymmetry using the function regtest in metafor. There was significant asymmetry among our studies (z = -5.99, p-value <0.001; Supplementary Figure S1). Though this asymmetry may result from chance alone given the low number of studies in the dataset, it is plausible that our set of studies does not represent an unbiased sample (Koricheva et al., 2013). The bias may arise from publication bias against non-significant results, system heterogeneity across studies, or that ecologists tend to set-up labor-intensive experiments only when they expect to see a result. With this caveat in mind, we pursued a simple analysis of top-down effects of waterfowl on eelgrass from the data currently reported in the literature.

We explored potential explanatory variables that could account for the observed heterogeneity among the effect sizes using linear mixed-effects models. Explanatory variables included latitude, longitude, waterfowl species, bird density (measured as individuals per hectare), Zostera response variable, Zostera species, area sampled, and the time elapsed between when the experiment began and when the maximum effect was observed. For studies that only reported experiment start dates and observations as a month in a year, we coded the date as the first day of the month and calculated the time elapsed since the experiment began in number of days. Several studies involved either a goose or swan combined with a duck species (e.g., Tubbs and Tubbs, 1983; Madsen, 1988; Gayet et al., 2012) but never a goose and swan together. No study included a duck species exclusively, so for waterfowl species, we further classified whether the study included a goose or a swan species. Similarly, we grouped the Zostera response variable into categories that classified whether the response estimated aboveground biomass, belowground biomass, or total biomass. We extracted data on explanatory variables directly from the study. We tested whether each predictor variable influenced the effect size in separate models by including the predictor term in the rma function described above. No predictors were combined in a single model due to limitations of sample size. We estimated the fit of each individual model to the data using Akaike information criteria corrected for small sample size (AICc) along with an accompanying pseudo R² statistic. We calculated pseudo R² as the explained proportion of heterogeneity in effect sizes from the mixed-effects model with the explanatory variable relative to the random-effects model with no explanatory variable (Raudenbush, 2009).

Bird diets were compiled from the previously described literature review, including backward and forward citation searches from diet data compiled in Olsen (2015). Data were included from studies on gut contents, fecal analyses, and observations of feeding, as well as consulting more general field guides (Poole, 2005). For birds consuming Zostera, their use was assessed in terms of what fraction of their diet contained Zostera at particular observation times. These were categorized as “dominant,” in which Zostera composed 50–100% of diet during at least one season; “frequent,” in which Zostera is 5–50% and unlikely to be a dominant food source; or “incidental” in which <5% or rare observations of consumption occur (Table 3). Where noted in studies, we also recorded the seasonal stage of primary consumption (i.e., migration, overwintering) and the part of the plant consumed (i.e., leaves, rhizomes). See Supplementary Table S2 for details.

We identified 76 papers via Web of Science and an additional seven papers through forward/backward citation searches and expert consultation. This yielded a total candidate study list of 83 papers. From this list, we retained all publications that satisfied the unique criteria we developed for each analysis (see section “Methods”). This resulted in 10 publications used for investigating the relationship between waterfowl and Zostera abundance, 11 publications used for the top-down meta-analysis, and 32 publications used for the diet assessment (see Figure 2 and Supplementary Table S2). Datasets used in the abundance relationship and the top-down analyses included four Zostera species: Z. japonica (one publication), Z. marina (10 publications, including two under the name Z. angustifolia), Z. noltei (nine publications), and Z. muelleri (two publications). Five publications included multiple species of Zostera (Z. noltei and Z. marina). Figure 3 illustrates the distribution of the Zostera genus, the distribution of the primary waterfowl herbivores, and the study locations from the meta-analyses.

Waterfowl and Zostera abundance tended to track each other spatiotemporally (positive correlations, r > 0), but only when variability in Zostera was substantial. Accordingly, Zostera CV was a significant modifier of the waterfowl–Zostera correlation, and the amount of correlation between the two taxa increased with variability in Zostera (Table 1 and Figure 4). The data sets from which correlation coefficients and Zostera CV were calculated are presented in the Supplementary Figure S2. When the coefficient of variation for Zostera was very low, we observed several negative correlations between waterfowl and Zostera abundance in one study (Percival et al., 1998), possibly reflecting depletion by birds. While Zostera variation significantly predicted the correlation, the scale of the study did not (Table 1). However, with the inclusion of these two predictors, Cochran’s Q-test showed no additional heterogeneity in the meta-analysis (QE₁₂ = 14.45, p-value = 0.27).

Effects of waterfowl on Zostera ranged from strongly negative (e.g., an 81% reduction in rhizome biomass from Tinkler et al., 2009) to slightly positive (an 18% biomass increase in rhizome biomass from Dixon, 2009 at one site) with an overall effect size estimated at -3.96 (95% confidence interval of -5.72 to -2.20) (Figure 5). The random effects model without predictors showed a large amount of heterogeneity in effect sizes among the studies (Q₁₃ = 84.78, p-value <0.001), leading us to consider potential predictor variables to explain this heterogeneity (see section “Methods”). Among the models tested, the model of effect size as a function of bird density had the greatest AICc support, but explained only 3.5% of the variance among studies (Table 2). The magnitude of the negative effect of waterfowl on Zostera abundance tended to increase with increasing bird density (Figure 6A), though the slope was not statistically significant (estimate = -0.11, Z = -0.64, p-value = 0.53; note that excluding Madsen, 1988 results in a significant slope with an estimate = -0.95, Z = -3.57, p-value = 0.0004). Other models had very little support, but we report two models with significant predictors that ranked second and third in our model comparison according to the AICc values. Bird arrival and bird type (goose versus swan) both explained some of the heterogeneity in effect sizes (27.2 and 18.5%, respectively), but they largely contained the same information. All geese arrived in autumn with one exception (Kiriakopolos, 2013), and all swans arrived in summer, therefore we cannot distinguish these two models. Overall, studies with goose species had stronger negative effects on Zostera than studies involving swan species (Figure 6B).

The number of avian species documented as consumers of Zostera greatly exceeds the number that have been studied for top-down and bottom-up effects. In total, we identified 39 species and subspecies of waterfowl that included Zostera in their diet at any one location or time point: nine dabbling ducks, 16 diving ducks, six geese (including three subspecies of brant, differentiated by distinct migratory pathways and feeding behaviors), five swans, two rails, and one wader species (Table 3 and Supplementary Table S2). Nearly all identified species are generalist herbivores that consume a variety of plant and non-plant food resources aside from Zostera. Zostera is a dominant component of the diets of six species (including the three subspecies of brant), frequent in the diets of 20 species, and rarely consumed by 13 additional waterfowl (Table 3). Except for the mute swan (Cygnus olor), all the heaviest consumers are represented in our correlation and top-down analyses. The majority of infrequent and rare consumers are dabbling and diving ducks. Comparing by species of seagrass consumed, Z. marina is consumed by the greatest number of bird species (27, many from Cottam, 1939), followed by five species for Z. japonica, three for Z. noltei, one for Z. muelleri, and ten in which only Zostera spp. was noted (Table 3 and Supplementary Table S2). While we had anticipated that different taxa might use different feeding techniques (e.g., grazing of leaves versus grubbing for rhizomes), we found that different sources documented different feeding behaviors for the same species. These changes in feeding behavior might reflect variation in whether leaves or rhizomes have higher nutritional value, or when loss of aboveground biomass necessitates a shift to grubbing for rhizomes (Mathers et al., 1998). Of the 39 avian species consuming Zostera, 13 are noted to interact with the seagrass primarily during a specific season, either migration or overwintering (Supplementary Table S2).

Our review provides a quantitative perspective on the reciprocal linkages between waterfowl species and temperate seagrasses. The published literature provides multiple cases in which waterfowl abundances are closely linked to Zostera abundance, while also showing that herbivorous waterfowl can have strong top-down effects on Zostera. Of the 39 avian species documented to consume Zostera, only a small fraction (12 species) are represented in the studies suitable for inclusion in this meta-analysis. Below we discuss these results and reflect on how the effects of waterfowl herbivory could extend beyond direct consumption. We also discuss the implications of waterfowl–Zostera interactions for conservation and identify areas of future research.

All authors contributed to the development of ideas, collection of literature, and data extraction. CH and PJ provided unpublished data for inclusion in the meta-analyses. JR conducted the correlation meta-analysis, AH conducted the diet assessment, NK and MW conducted the top-down meta-analysis. NK, AH, MW, and JR co-wrote the manuscript with significant input from all authors.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.