Forager (predator) abundance may mediate feeding rates in wading birds. Yet, when modeled, feeding rates are typically derived from the purely prey-dependent Holling Type II (HoII) functional response model. Estimates of feeding rates are necessary to evaluate wading bird foraging strategies and their role in food webs; thus, models that incorporate predator dependence warrant consideration. Here, data collected in a mangrove swamp in Puerto Rico in 1994 were reanalyzed, reporting feeding rates for mixed-species flocks after comparing fits of the HoII model, as used in the original work, to the Beddington-DeAngelis (BD) and Crowley-Martin (CM) predator-dependent models. Model CM received most support (AICc wi = 0.44), but models BD and HoII were plausible alternatives (ΔAICc ≤ 2). Results suggested that feeding rates were constrained by predator abundance. Reductions in rates were attributed to interference, which was consistent with the independently observed increase in aggression as flock size increased (P < 0.05). Substantial discrepancies between the CM and HoII models were possible depending on flock sizes used to model feeding rates. However, inferences derived from the HoII model, as used in the original work, were sound. While Holling's Type II and other purely prey-dependent models have fostered advances in wading bird foraging ecology, evaluating models that incorporate predator dependence could lead to a more adequate description of data and processes of interest. The mechanistic bases used to derive models used here lead to biologically interpretable results and advance understanding of wading bird foraging ecology.
Estimates of feeding rates are essential to understand predator-prey relationships and trade-offs among foraging strategies in wading birds (Charnov 1976; Kushlan 1978; Draulans 1987; Erwin 1985, 1989). Also, estimates of feeding rates are necessary to quantify the impact of wading birds on fish communities in nursery potential studies. These studies assess potential factors limiting colonization of juvenile life stages to nursery habitats and their demographic implications (Gilmore et al. 1983; Miranda 1995). Estuaries are productive systems (Odum et al. 1982; Hamilton and Snedaker 1984) that may serve as nursery habitat for many species of fish (Austin 1971; Gilmore et al. 1983). Nursery habitats provide juvenile fish shelter from adult predatory fish and high quality growth habitat (Zale and Merrifield 1989; Aliaume et al. 1997; Aliaume et al. 2000; Zerbi et al. 2001). Access to such habitats include the shallow reaches of mangrove swamps. These swamps may expose the fish community to other predators and wading birds are among the most prominent (Lopez et al. 1988; Kushlan et al. 1985; Kent 1986a; Ramo and Busto 1993).
The most prevalent model used to estimate feeding rates in the general ecological and waterbird literature is Holling's (1959) Type II functional response (e.g. Draulans 1987; Piersma et al. 1995; Skalski and Gilliam 2001; Durant et al. 2003; Stillman et al. 2003; Gillings et al. 2007). Despite its seminal role in advancing our understanding of predator-prey interactions (Brown 1991), Holling's model is a purely prey-dependent model; in the model, feeding rates are unaffected by predator densities and competition among predators occurs only via depletion of prey. However, feeding rates might be affected by predator interference or changes in prey behavior in a predator-dependent context. These factors may mediate intake rates in wading birds and are central to hypotheses predicting patch use and revisitation rates and foraging flock size dynamics (e.g. Kushlan 1978; Kent 1986b; Erwin 1989; Norris and Johnstone 1998). Skalski and Gilliam (2001) reviewed and compared several functional response models that permit assessing the influence of such factors on intake rates. We believe that these alternative models warrant consideration in future wading bird foraging ecology studies.
We illustrate our point with data collected at a mangrove swamp in southwestern Puerto Rico in January–March 1994 (Miranda 1995; Miranda and Collazo 1997a, b). A stated objective of the study was to assess the nursery potential of the Boquerón Wildlife Refuge for snook (Centropomus spp.) and tarpon (Megalops atlanticus), and thus, special attention was given to estimating per capita feeding rates. Estimates were used to determine mortality rates and the number of juvenile fish removed by mixed-species flocks of Great Egrets (Ardea alba), Tricolored Herons (Egretta tricolor) and Snowy Egrets (E. thula). On average two Little Blue Herons (E. caerulea) occurred per flock, but these were not considered in the assessment because they fed almost exclusively on fiddler crabs (Uca spp.) and because the sampling gear used to estimate prey density was not suitable to sample crabs (Miranda 1995, Miranda and Collazo 1997a, b). We re-assessed the support in the data for the Holling's Type II model (HoII), as used in the original study (Miranda 1995), together with two alternative models that include predator interference. One was derived independently by Beddington (1975) and DeAngelis et al. (1975), henceforth the BD model; the other was developed by Crowley and Martin (1989), henceforth the CM model. In essence we assessed three competing hypotheses, one predicting changes in feeding rates as a function of prey density and two others on the basis of both prey density and predator abundance. A major distinction between the BD and CM models is that the former assumes that time devoted to prey handling and interference are exclusive activities; the CM model does not (Skalski and Gilliam 2001). The inclusion of these models in the candidate set was justified because aggressive behavior (interference) among members of a flock increased with flock size (Miranda 1995). We used an information-theoretic model selection process to compare the relative support among models in the set (Burnham and Anderson 2002). We also report prey density and number of interactions among members of foraging flocks to describe the foraging and social context of the study. Finally, we discuss the appropriateness of using feeding rates as derived by Miranda (1995) to assess the impact of waders on the fish community at the Boquerón Wildlife Refuge and the value of considering multiple models when evaluating predictions about wading bird foraging strategies.
METHODS
Study Site
The Boquerón Wildlife Refuge is located in southwestern Puerto Rico, east of Boquerón Bay, Cabo Rojo (18°00′N, 67°08′W). The refuge is approximately 3 km2 and comprised mostly of shallow water lagoons bordered by Red (Rhizophora mangle), Black (Avicenia germinans) and White mangrove (Laguncularia racemosa). Data were collected in the southwestern portion of the refuge (Fig. 1). The area consists of 83 ha of mangrove forest demarked on the east by dikes that enclose most of the estuary. Thirty percent of the 83 ha consisted of open pools of an average depth of 11 cm, with low canopy and nematophore cover (i.e., cover ≤20%; Miranda 1995; Miranda and Collazo 1997a).
Data Collection and Analyses
Data were collected between January and March 1994. This period was typified by seasonal flooding and movement of juvenile fish into the shallow reaches of mangrove swamps (Miranda and Collazo 1997a; Aliaume et al. 1997, 2000). We monitored 11 mixed-species flocks, ranging in size from 5–70 individuals. A total of 103 individuals of Great and Snowy Egrets and Tricolored Herons were observed from blinds, approximately 10 m from birds, between sunrise and 10.00 h, using 8 × 40 binoculars. Observation bouts did not exceed 75 min. Before starting every sampling period, we recorded the number of wading birds. A foraging flock was defined as two or more individuals foraging with interindividual distances ≤5 m (Wiggins 1991; Master et al. 1993). We started observations by arbitrarily selecting a focal individual from the “center” of the flock. Focal birds of all species were chosen randomly and observed for 5 min (Draulans 1987), recording the number of successful attempts and number of aggressive encounters. Aggressions were defined as overt antagonistic behavior (e.g. pecking, chase, threat posture). Data were collected on as many individuals as possible. We cannot rule out the possibility that an individual may have been observed more than once during an observation period, but can assert that none were observed consecutively. Wilcoxon tests were used to compare average number of aggressive encounters among species. Linear regression was used to determine if the average number of aggressions per mixed-species flock increased with flock size. Data for the latter analysis were log-transformed to meet homogeneity of variance assumptions (Levene's test, P > 0.05).
Figure 1.
Boquerón Wildlife Refuge, Puerto Rico showing study site on the southwestern quadrant of the refuge. Foraging observations were made in the shallow reaches of the mangrove swamps between the east and west dikes, January—March 1994.
Prey density (m-2) estimates were obtained immediately following each observation period. Prey consisted of fish and shrimp, with fish being ≥5 times more numerous than shrimp in samples (Miranda and Collazo 1997a). Prey was sampled using a modified 1-m2 throw-trap. Prey density estimates were consistent among sites, but likely biased low due to predator aversion by the prey and removal (Miranda and Collazo 1997a). The same sampling technique was used to describe the foraging environment throughout the study area and use patterns of waders by comparing used and non-used sites (N = 38; Miranda and Collazo 1997a). A “used” site was defined as a location where two or more wading birds actively fed. A non-used site was randomly located at approximately 50 m from the used site, staying within similar habitat features. Paired estimates were compared using a Wilcoxon test because data did not meet homogeneity of variance assumption.
We fitted data to three functional response models using PROC NLMIXED (SAS 2002). The HoII model (Holling 1959) is defined as f (N, P) = aN/(1+bN); where a (units: 1/time) and b (units: 1/prey) are expressions of capture rate and handling time (handing time = b/a) and N is prey density. The BD model is defined as f (N, P) = aN/(1+bN+c (P-1)); where c (units: 1/predator) is an expression of the magnitude of predator-predator interference and P is predator abundance (i.e., number of waders per flock). In this equation, we used P-1 because a focal predator is assumed not to be interfering with itself (Skalski and Gilliam 2001). The CM model is defined as f(N, P) = aN/((1+bN) (1+c(P-1)); where c and P have the same meaning as in the BD model. The BD model and the CM model differ in how they assume interference to occur. In the BD model interference is depicted as occurring via encounters between predators, with encounters resulting in “time wasted.” Interference enters into the equation in a way analogous to handling time: just as time spent handling prey is subtracted from searching in the HoII model, the BD model introduces a term that depicts time lost from searching due to interference. In the BD model, a predator cannot be both handling prey and experiencing interference simultaneously. Hence, holding P constant and letting N become large, interference is modeled as becoming negligible (Skalski and Gilliam 2001). The CM model depicts interference differently. Rearranging the CM model as [aN/((1+bN)]/ (1+c (P-1)), the CM model can be viewed as the HoII model divided by an interference term. Hence, the interference is depicted as a penalty even at arbitrarily high prey density. Thus, a difference between the BD model and the CM model is that the asymptotic feeding rate as N increases does not depend on P in the BD model, but the asymptote does depend on P in the CM model. The BD model and the CM model both reduce to the HoII model for P = 1.
We report estimates of feeding rates (i.e., successful attempts) per capita for mixed-species flocks, as did Miranda (1995). In the analyses, flocks were the primary sampling unit, not individuals within a flock. As such we decoupled and accounted for two sources of variation; that due to a random effect (flocks, N = 11) and the residual variance due to individuals within a flock. The latter is influenced by the fact that behavior within flocks is correlated (e.g. Petit 1987). We recognize that modeling species-specific feeding rates is valuable to discern the relative predation pressure exerted by species as well as to gain insights on mixed-species flock foraging ecology, but the number of flocks per species did not permit analyses. Miranda (1995) and Miranda and Collazo (1997a) showed that physical parameters (e.g. salinity, temperature, depth) between paired prey sampled sites throughout the study area did not differ significantly. Thus we felt no need to adjust prey density by an index of vulnerability (Kersten et al. 1991; Gawlik 2002). Support in the data for models in the candidate set was assessed using the Akaike's Information Criterion (Burnham and Anderson 2002). Models with ΔAICc ≤ 2 were considered models with highest support.
We summarized model results in three ways. First we plotted the observed mean and predicted successful feeding attempts/5 min derived from the Hollings (1959) functional response model and from the alternative models. Second, we depicted the interplay between feeding rates and three flock sizes. We arbitrarily used flock sizes of 13, 26 and 52 individuals, which represented the average flock size (26; Miranda 1995) and half and twice the average flock size for the lower and upper values. Prey density used to generate feeding rates were within the range recorded in the field (range = 35–375 m-2). Third, we evaluated the implications of different model results on the assessment conducted by Miranda (1995). This assessment consisted of two parts. First, we compared estimates of feeding rates using the HoII and the alternative model with highest support (AIC) using PROC NLMIXED (SAS 2002). Input values were the field-derived average estimates of mixed-species flock size (26) and prey density at used sites (131 m-2). Estimates were compared using a Z-test. Second, we plotted the difference between estimates of fish removed/capita obtained from the alternative model with highest support (CM or BD) and HoII model across a range of flock sizes. The rationale was to explore how inferences about wader impact on fisheries resources could be affected by varying flock sizes if only purely prey-dependent models are considered to estimate feeding rates. For this assessment we ran a simulation using 103 pairs of randomly selected of values of prey density and flock sizes within the range of observed values in the field. The sample size for this simulation was equal to the total number of observations made in the field (Miranda 1995). Predicted feeding rates were generated using PROC NLMIXED (SAS 2002). We report means (±SE). We used a P value of 0.05 to determine test significance.
RESULTS
Average prey density in used sites (131.4 ± 17.0 m-2) was significantly higher than in non-used sites (17.0 ± 5.1 m-2) (χ2 = 42.97, P < 0.001). Average prey density at sites where foraging observations were made was 180.81 ± 38.79 m-2 (range 35–375). The average size of flocks at those sites was 25.63 ± 6.20 individuals (range 5–70 individuals; N = 11 flocks). On average Great Egrets were the most abundant species in flocks (45% ± 0.05 of individuals), followed by Snowy Egrets (29% ± 0.05), Tricolored Herons (16% ± 0.05) and Little Blue Herons (8% ± 0.02).
The average number of aggressions per individual (1.90 ± 0.19) increased with flock size (F1,10 = 3.95, P < 0.05). Great Egrets (2.18 ± 0.22) and Tricolored Herons (2.04 ± 0.33) exhibited significantly higher mean aggressive encounters than Snowy Egrets (1.20 ± 0.44; 
= 9.89, P = 0.01). There were almost no interactions among Little Blue Herons (0.20 ± 0.12). More than half (56%) of all individuals displayed aggressive behavior towards other birds. Aggression was nine times more frequent (91.7%) towards conspecifics that towards other species.
Variation in feeding rates in mixed-species flocks was best explained by the CM model (AIC c wi = 0.44; Table 1). Plausible alternatives were the BD model (AIC c wi = 0.36), followed by the purely prey-dependent HoII model (AIC c wi = 0.20). Model parameters (coefficients) are listed in Table 2. Greater model support for the CM and BD models suggested that feeding rates were influenced by predator abundance. The negative influence of flock size on feeding rates was especially evident in two instances in the field data (Fig. 2). Both were when prey density was between 309–310 and 355–375, but flocks doubled in size within those prey densities. We also depicted the constraining influence of larger flocks on predicted feeding rates (Fig. 3). For a given prey density, say 155 m-2, the predicted intake would be about 30% lower for a flock of size 52 as compared to a flock of 26 birds. Feeding rates predicted by the HoII model at the mean prey density of 131 m-2 (used sites) was 5.49 ± 3.34, similar to the CM-derived estimate of 5.15 ± 2.96 at the same mean prey density and mean flock size of 26 individuals (Z = 0.17; P > 0.05). Our simulations showed that the same HoII model underestimated the number of prey removed/capita as flocks decreased in size or progressively tended to overestimate predicted removal rates as foraging flocks increased in size (Fig. 4).
Table 1.
Model selection results of three functional response models used to estimate per capita feeding rates of waders at the Boquerón Wildlife Refuge, Puerto Rico, Model support evaluated using ΔAICc, AICc weights (wi) and model likelihood (ML). K is the number of model parameters, Functional response models were: CM = Crowley-Martin (1989); BD = Beddington-DeAngelis (1975); HoII = Holling Type II (1959)
DISCUSSION
Wading birds in the mangrove swamps of the Boquerón Wildlife Refuge conformed to previously documented patterns of habitat use and foraging behavior. Waders foraged in sites with highest prey density as documented recently by Gawlik (2002) in an experimental setting and Stolen (2006) in salt marshes in central Florida. In a close parallel to our study, Kushlan (1976) reported that Wood Storks (Mycteria americana) in Florida foraged in sites where fish concentrations averaged 141 fish.m-2 as compared to 10 fish.m-2 in unused sites.
Re-analyses of Miranda's (1995) data yielded higher support for the CM model (AIC c w i = 0.44), followed closely by the BD model (AIC c w i = 0.36). Support for these models suggested that predator abundance influenced feeding rates, and we interpreted these results as the effects of interference. This interpretation was consistent with the increase in aggressive behavior with increasing flock size reported here and in previous studies (Takita and Dean 1984; Kent 1986b). We speculate that Great Egrets and Tricolored Herons contributed more to the support of the CM and BD models than Snowy Egrets as they exhibited greater aggressive behavior. Holling's (1959) model was also a plausible alternative (ΔAIC c ≤ 2). Thus, model averaging is an alternative and appropriate way to express feeding rates (Burnham and Anderson 2002). The equivocal evidence for predator interference might be rooted on unique interactions between field conditions and flocks at the Boquerón Wildlife Refuge. Admittedly, the number of flocks under observation was low, and thus may have prevented better discrimination among models. For this reason, we cannot overemphasize the importance of careful study design and sampling effort in future studies to adequately assess the evidence in the data for competing models. It is also important to bear in mind that true replication is measured by the number of flocks sampled, not by the number of focal individuals within flocks observed in the field. This is necessary to avoid inferences based on pseudoreplicated data (Hurlbert 1984).
Table 2.
Parameters (SE) of three functional response models used to estimate feeding rates of wading birds at the Boquerón Wildlife Refuge, Puerto Rico. Parameters a (units: 1/time), b (units: 1/prey) and c (units: 1/predator) are expressions of capture rate, handing time (b/a) and interference, respectively. Parameters s2u and s2e account for group (flock) and individuals within flock variance, respectively. Functional response models were: CM = Crowley-Martin (1989); BD = Beddington-DeAngelis (1975); HoII = Holling Type II (1959). df = 10; * = significant at 0.05.
Figure 2.
Predicted and mean observed successful feeding attempts/5 min derived from the Crowley-Martin (1989; saCM), Beddington-DeAngelis (1975; saBD) and Hollings (1959; saHoII) functional response models. Input parameters were the number of mixed-species foraging flocks and prey density (fish + shrimp) recorded at the Boquerón Wildlife Refuge, Puerto Rico, 1994.
Miranda (1995) assessed the potential impact of wading birds on the fish community at the Boqueron Wildlife Refuge by estimating the number of fish removed/hr by an average flock size (26) after adjusting for their dietary composition (Miranda and Collazo 1997b). Feeding rates were estimated at the average prey density of 131 m-2. At this prey density, the HoII (5.49 ± 3.34) predicted similar feeding rates than the CM (5.15 ± 2.96) model at the same prey density and predator abundance of 26 individuals. We conclude that because Miranda (1995) did not extrapolate feeding rates to other flock sizes or estuaries, his estimates of fish removal and inferences about wading bird impact on the fish community were sound. Our simulation showed that had evidence for a predator interference effect been stronger; failing to account for it could have led to substantial under- or overestimates of prey removed/capita depending on the flock size used to estimate feeding rates. Average feeding rates, however, will likely be of limited value to model the impact of wading birds on fish communities as studies increase in complexity (e.g. multiple estuaries). This is because systems under study will vary temporally, spatially or as a result of management actions. Circumstances at the Boquerón Wildlife Refuge represent a plausible example. A major factor undermining the refuge's nursery potential was its levees and dikes, built to improve waterfowl hunting opportunities in the 1970s (Miranda 1995). Such structures restrict fish movement to and from the estuary, increasing prey vulnerability (Morales and Pacheco 1986; Crowder et al. 1997). Increasing access to the estuary is a reasonable management option and could fundamentally change the dynamics of predators-prey interactions. These circumstances favor modeling approaches that can accommodate such dynamic states (e.g. stochastic-dynamic) and stress the value of evaluating multiple functional response models to ensure an adequate description of processes of interest.
Figure 3.
Predicted successful feeding attempts/5 min as a function of flocks of 13, 26 and 52 individuals, Prey density ranged from 35 to 350/m2. Predictions were based on the Crowley-Martin (CM, 1989) functional response model, Model coefficients were derived from fitting the model to field data collected at the Boquerón Wildlife Refuge, Puerto Rico, 1994.
Figure 4.
Difference between estimates of predicted fish removed/capita for mixed-species flocks obtained from the Crowley-Martin (1989; CM) and Holling Type II (1959; HoII) functional response models. Predicted values were obtained from fitting models to randomly selected values of flock size and prey density from within the range of observed values in the field. Mixed-species flocks ranged in size from 5 to 70 individuals and prey densities from 35 to 375/m2.
Foraging flocks have been of interest to ecologists because they may lead to higher foraging efficiency (e.g. Power 1985; Cezilly and Hafner 1990; Wiggins 1991; Master et al. 1993; Stolen 2006). The possibility has been the impetus to quantify tradeoffs between flock size and feeding rates (e.g. Kushlan 1978; Stolen 2006). Because interference might influence this process (e.g. Russell 1978; Kent 1986b; Petit 1987), predator-dependent models provide a suitable framework to test predictions about flock size dynamics. In this vein, Erwin (1989) advanced predictions about how temporal changes in encounter rates affect patch revisitation schedules and feeding rates for four different prey types (e.g. schooling fish). Pertinent to our discussion here is that evidence to support or reject predictions is based on estimates of feeding rates, and these can be influenced by interference (Erwin 1989). Moreover, parameters such as attack and depletion rates could be used as complementary evidence because they serve to gauge patch quality. We argue that predator-dependent models warrant consideration because, if supported by the data, feeding and attack rates are directly derived from models. Appropriate model selection would also lead to more accurate estimation of depletion rates (e.g. removal/capita). Finally, while our questions focused on the possible influence of predator interference on feeding rates, questions about facilitation are also of interest as noted above and could be explored with the BD and CM models. For example, an extension of these models, such that c is treated as a function of predator abundance (e.g. ([c0+ c1(P-1)](P-1)) rather than as a constant as it was in this work, would allow exploration of interference/facilitation across flock sizes.
Purely prey-dependent models, general linear models, and the use of indirect measures (e.g. giving-up-densities or times) have been helpful in advancing our understanding of wading bird foraging ecology (e.g. Erwin 1985; Draulans 1987; Kersten et al. 1991; Gawlik 2002). We suggest that the predator-dependent models used here are a useful benchmark to achieve further progress. This is because modeled parameters are germane to predictions about wading bird foraging ecology and to applied questions such as assessments of nursery potential. Furthermore, the mechanistic bases behind the derivation of models lead to biologically interpretable results (Skalski and Gilliam 2001).
ACKNOWLEDGMENTS
We thank M. Sanderson, J. Cruz, M. Laboy, A. Martinez, W. Miranda, and J. Morales for assistance in the field. The work was supported by a US Fish and Wildlife Service Federal Aid Grant (F-33) through the Puerto Rico Department of Natural and Environmental Resources. We thank C. Arellano for statistical advice, and R. M. Erwin, T. R. Simons and an anonymous referee for reviewing the manuscript.
