Predator–prey dynamics of bald eagles and glaucous‐winged gulls at Protection Island, Washington, USA

Abstract Bald eagle (Haliaeetus leucocephalus) populations in North America rebounded in the latter part of the twentieth century, the result of tightened protection and outlawing of pesticides such as DDT. An unintended consequence of recovery may be a negative impact on seabirds. During the 1980s, few bald eagles disturbed a large glaucous‐winged gull (Larus glaucescens) colony on Protection Island, Washington, USA, in the Salish Sea. Breeding gull numbers in this colony rose nearly 50% during the 1980s and early 1990s. Beginning in the 1990s, a dramatic increase in bald eagle activity ensued within the colony, after which began a significant decline in gull numbers. To examine whether trends in the gull colony could be explained by eagle activity, we fit a Lotka–Volterra‐type predator–prey model to gull nest count data and Washington State eagle territory data collected in most years between 1980 and 2016. Both species were assumed to grow logistically in the absence of the other. The model fits the data with generalized R 2 = 0.82, supporting the hypothesis that gull dynamics were due largely to eagle population dynamics. Point estimates of the model parameters indicated approach to stable coexistence. Within the 95% confidence intervals for the parameters, however, 11.0% of bootstrapped parameter vectors predicted gull colony extinction. Our results suggest that the effects of bald eagle activity on the dynamics of a large gull colony were explained by a predator–prey relationship that included the possibility of coexistence but also the possibility of gull colony extinction. This study serves as a cautionary exploration of the future, not only for gulls on Protection Island, but for other seabirds in the Salish Sea. Managers should monitor numbers of nests in seabird colonies as well as eagle activity within colonies to document trends that may lead to colony extinction.


| INTRODUC TI ON
After years of decline, bald eagle (Haliaeetus leucocephalus) populations throughout North America rebounded in the latter part of the twentieth century following tightened protection, reduction in the use of lead shot by hunters, and the outlawing of pesticides such as DDT (Hipfner et al., 2012;Watson, Stinson, McAllister, & Owens, 2002). This recovery has provided one of the great success stories of the conservation movement (Millar & Lynch, 2006).
Nowhere has recovery been more pronounced than in the Pacific Northwest of North America where inland waterways such as the Salish Sea, Columbia River, and scores of smaller lakes and streams provide ideal perching, hunting, and nesting opportunities for these raptors (Elliott, Elliott, Wilson, Jones, & Stenerson, 2011;Stinson, Watson, & McAllister, 2001;Watson, 2002;Watson et al., 2002).
An unintended consequence of bald eagle recovery has been the negative impact on seabirds, which are already stressed by overfishing, gill netting, and habitat destruction (Atkins & Heneman, 1987;Blight, Drever, & Arcese, 2015). Although populations of some seabirds may be declining to historic levels (Elliott et al., 2011), local populations of seabirds such as common murres (Uria aalge) may be threatened (Parrish, Marvier, & Paine, 2001). Numbers of salmon and other fish traditionally eaten by wintering bald eagles have plummeted in recent years, affecting eagle survival and possibly resulting in their shift to other food sources (Elliott et al., 2011).
Seabirds always have formed part of the diet of eagles (Stalmaster, 1987), but increasing numbers of eagles and concurrent prey fish shortages have resulted in increased eagle foraging on waterfowl, a cause for concern among ornithologists (Elliott et al., 2011;Hipfner et al., 2012;Moul & Gebauer, 2002;Parrish et al., 2001;Sullivan, Hazlitt, & Lemon, 2002;Vennesland & Butler, 2004;White, Heath, & Gisborne, 2006). Potential impacts of bald eagle populations on marine food-web structure appear to be due to resident eagles, rather than overwintering eagles, and the rates at which they consume seabirds as prey (Harvey, Good, & Pearson, 2012).
Bald eagles can impact seabirds both directly and indirectly (Hipfner et al., 2012;Parrish et al., 2001). The most obvious direct effect is the killing and eating of adults, juveniles, and eggs (DeGange & Nelson, 1982;Hayward, Galusha, & Henson, 2010;Hayward, Gillett, Amlaner, & Stout, 1977). A second direct effect is the extra expenditure of energy needed for nesting or feeding in the presence of eagles Parrish et al., 2001).
An indirect effect results when disturbances displace breeding adults from their nests and expose unprotected eggs and young to other predators (Hayward et al., 2010). A second type of indirect effect involves changes in distribution patterns in response to the presence of eagles. For example, diving waterbirds in the Strait of Georgia moved away from inshore waters, and dabbling ducks formed larger aggregations inshore and were more vigilant, in response to increased eagle presence (Middleton, Butler, & Davidson, 2018).
From 1900 to the early 1980 s, breeding populations of glaucouswinged gulls (Larus glaucescens) markedly increased in the Georgia Basin of the Salish Sea, British Columbia. By 2010, populations had declined to about 50% of peak levels (Blight et al., 2015;Sullivan et al., 2002). A study that incorporated more southern areas of the Salish Sea also reported overall declines from 1975 to 2007 (Bower, 2009 has functioned as a breeding center for marine birds since at least the 1940s (Power, 1976). A large glaucous-winged gull colony had become established by the early 1960s (Richardson, 1961). Today, Rhinoceros auklets (Cerorhinca monocerata), glaucous-winged gulls, pigeon guillemots (Cepphus columba), and harbor seals (Phoca vitulina) breed there in large numbers. Adult auklets, gulls, and guillemots, as well as the eggs and chicks of gulls and the afterbirths and pups of seals, all serve as food for nesting and visiting eagles (Cowles, Galusha, & Hayward, 2012;Hayward, 2009;Hayward et al., 2010).
During the 1980s, few eagle disturbances of the gull colony on Protection Island's Violet Point were noted, and from 1980 to 1993 gull nest numbers increased by 37% (3,796-5,189; Table 1).
Beginning in the 1990s, however, a dramatic rise in bald eagle activity over and within the colony was observed (Galusha & Hayward, 2002;Hayward et al., 2010), with a significant decline in numbers of breeding gulls at the site (Cowles et al., 2012). Bald eagles constitute the only significant source of interspecific predation on the gulls in this colony .

| Classic Lotka-Volterra predator-prey model
In their classic paper on Canada lynx (Lynx canadensis) and snowshoe hare (Lepus americanus) populations, Elton and Nicholson (1942) used 100-year records from the Hudson Bay Company on the numbers of pelts purchased from trappers. The classic predator-prey cycles of theoretical ecology (Figure 1a), often illustrated with lynx-hare data, are produced by the Lotka-Volterra predator-prey ordinary differential equation model (Henson, 2012;Lotka, 1925;Volterra, 1926) Here the "prime" denotes the derivative with respect to time, G and E refer to numbers or densities of prey (gulls, in this context) and predators (eagles), a > 0 is the per capita growth rate of the prey population in the absence of the predators, b > 0 is the per capita decline rate of the predator population in the absence of prey, α > 0 is the predation rate (the probability per unit time that a given prey individual will be taken by a given predator), and β > 0 is the conversion rate of prey into predators.
The classic predator-prey model (1) has two major deficiencies.
First, the prey population grows exponentially, without bound, in the absence of predators; and second, the predator population declines exponentially to extinction in the absence of the prey. Neither of these scenarios is feasible in most ecological communities because population growth is always eventually bounded by self-limitation, and predators usually can switch prey and hence do not decline to extinction with the removal of a single prey species. This is true for bald eagles, which are considered opportunistic foragers (Buehler, 2000).

| Gull-eagle predator-prey model
To examine the relationship between the Violet Point gull colony and eagle activity in terms of a predator-prey interaction, we modified the Lotka-Volterra predator-prey model (1) to include a multiple prey base for eagles and self-limitation terms for both gulls and eagles. In particular, we used the Lotka-Voterra-type ordinary differential equation model (Ricklefs, 1990) where G and E are the numbers of gull and eagle pairs, respectively, rather than individual animals as in the original Lotka-Volterra model.
Here r > 0 and s > 0 are the inherent per capita growth rates for gulls and eagles, respectively, at small population sizes, and r/K > 0 and s/C > 0 are their rates of self-limitation. The parameter α > 0 denotes the predation rate of eagles on gulls, and β is the conversion rate of gulls into eagle births. In the absence of the other species (when α = β = 0), each species grows logistically with carrying capacities K > 0 for gulls and C > 0 for eagles.
Model (2)  (1) (2) • If r > αC, then both Ḡ and Ē are positive in Equation (3), and so the coexistence state (3) is biologically feasible. In this case, the equilibria (0, 0), (K, 0), and (0, C) are unstable and the coexistence equilibrium (Ḡ,Ē) is stable. This equilibrium is either a stable spiral or a stable node. That is, gulls and eagles either approach the coexistence equilibrium through damped predator-prey oscillations ( Figure 1b), or else they approach equilibrium in a nonoscillatory fashion. In the latter case, early transient dynamics may resemble predator-prey oscillations, but the oscillations do not persist ( Figure 1c).
• If r < αC, the coexistence equilibrium (3) is not biologically feasible (because the equilibrium number of gulls Ḡ is negative). The equilibria (0, 0) and (K, 0) are still unstable, but the equilibrium (0, C) in which gulls are absent and eagles are at carrying capacity C is now stable. That is, gulls approach extinction, whereas eagles approach their carrying capacity C. Early transient dynamics may resemble predator-prey oscillations before gulls eventually go extinct ( Figure 1d).
The biological interpretation of these alternatives is the following.
The number r is the inherent net reproductive rate of gulls, and the number αC is the rate at which gulls are taken by C eagle pairs. If the inherent net reproductive rate of gulls is larger than the rate at which gulls can be taken by C eagle pairs, then gulls and eagles both survive and approach a positive coexistence equilibrium. If, however, the inherent net reproductive rate of gulls is smaller than the rate at which gulls can be taken by C eagle pairs, then gulls go extinct and eagles approach their carrying capacity C.

| Gull nest count data for Violet Point, Protection Island
We used glaucous-winged gull nest count data collected between 1980 and 2016 at a large breeding colony on Violet Point, Protection Island National Wildlife Refuge, Washington (48°07′40″N, 122°55′3″W), which lies at the eastern end of the Strait of Juan de Fuca in the Salish Sea ( Figure 2). The Violet Point colony is populated by glaucous-winged gulls and glaucous-winged gull × western gull (L. occidentalis) hybrids (Bell, 1996(Bell, , 1997. Most of these hybrids resemble glaucous-winged gulls more than western gulls Moncrieff, Megna, Hayward, & Henson, 2013); hence, we refer to these birds collectively as glaucous-winged gulls.
Gull counts from 1980 through 2014 were carried out in squad fashion as described in Galusha, Vorvick, Opp, and Vorvick (1987).

| Occupied eagle territory data for Washington State
We obtained data for the number of breeding territories occupied

| Counts of nesting and non-nesting eagles on Protection Island
We obtained annual maximum numbers of subadult and adult bald eagles observed simultaneously on Protection Island from 1993  (Table 1)

| Washington eagle territories as a proxy for eagles on Protection Island: correlation analysis and Poisson regression
Numbers of occupied eagle territories in Washington State were larger and relatively less noisy than numbers of eagles observed on Protection Island. Hence, for model fitting we wished to use a scaled version of the statewide eagle data as a proxy for the Protection Island eagle data. To determine whether we could use the numbers of occupied eagle territories in Washington State as a proxy for eagle activity on Protection Island, we performed a correlation analysis on the number of eagles observed on Protection Island and the number of occupied eagle territories in Washington State for the eight years these data overlapped (1993-1998, 2001, 2005;

| Model parameterization
Ecologists must consider several factors when fitting theoretical models to time series data. Populations are subject to both demographic and environmental stochasticity resulting in "process noise" (Dennis, Ponciano, Lele, Taper, & Staples, 2006).
Inaccuracies in the estimates of population densities also result in measurement error (Carpenter, Cottingham, & Stow, 1994).
Although methods have been developed to deal with these factors (e.g., Valpine & Hastings, 2002), the data demands for these methods are sometimes prohibitive. In particular, the data used here raise challenges for model parameter estimation. There are gaps in the gull and eagle time series, and in some years, data are not available for both species. The sample sizes differ for gulls and eagles, and the numbers for the two species are different in magnitude. These difficulties preclude the application of techniques based on autoregressive time series that underlie many of the methods and software commonly in use (Bolker et al., 2013).
To estimate the parameters in model (2), we first scaled the state variables G and E, as well as the data, by dividing by the observed standard deviations σ g and σ e , respectively (that is, Ĝ = G∕ g and Ê = E∕ e ), to scale the data to comparable magnitudes. It follows that the derivatives are Ĝ� = G � ∕ g and Ê� = E � ∕ e , and so one can rewrite model (2) in terms of the scaled variables: To fit model (4) to the scaled data, we used the ode45 differential equation solver in Matlab ® to produce predicted model trajectories from 1980 to 2016. We treated the (scaled) initial conditions ĝ 0 and ê 0 as parameters to be estimated. Given a vector of parameter estimates = ĝ 0 ,ê 0 ,r,K, ,s,C, , we computed residuals on the log scale to account for environmental stochasticity, which (4) F I G U R E 4 Histograms of the equilibria for both species, based on 2,000 bootstrapped parameter estimates F I G U R E 3 Observed data and model predictions. Observed data (symbols) and predictions of models (2) and (11)  is approximately additive on the log scale (Cushing, Costantino, Dennis, Desharnais, & Henson, 2003;Dennis, Munholland, & Scott, 1991): where ĝ t ( ) and ê t ( ) are predicted values obtained by numerically integrating model (4) from year 1980 to 2016, using parameters θ, the parameters ĝ 0 and ê 0 as initial conditions, and the observed standard deviations σ g and σ e for gulls and eagles, respectively. We obtained best fit parameters ̃ by minimizing the sum of the root mean squares (RMS) of the residuals as a function of θ, where n g and n e are the number of residuals for gulls and eagles, respectively, using the fminsearch downhill search algorithm in Matlab ® .
The fitting method described above is based on a number of considerations. First, the two time series are not paired; in many years, estimates are available for only one of the two species.
Therefore, we cannot view the data for every year as a traditional bivariate observation which would allow a traditional sum of squared errors, and we also cannot use one-step predictions in computing residuals. Second, each species has different numbers of observations, so the sum of squared residuals must be scaled by the number of observations. Otherwise, the parameter estimates would bias the species with the most observations. Third, we cannot estimate parameters using separate RMS values because the equations are coupled, so the fit of one species affects the fit of the other. Fourth, the overall magnitudes of the two species differ; hence, we scaled the data by the standard deviations so that the two terms in the RMS equation would be commensurate and the parameter estimates would not be biased in favor of fitting the species with overall higher numbers.
We performed diagnostic analyses of the gull and eagle residuals to check for independence and normality. We plotted the residuals as a function of time and examined normal quantile-quantile plots for departures from normality. We computed first-and second-order autocorrelations of the gull and eagle residuals that were separated by one and two years, respectively, and tested these correlations for significance. We also computed the Shapiro-Wilk test statistic for normality (Shapiro & Wilk, 1965).

| Goodness-of-Fit
We used a generalized R 2 to check the goodness of fit of the scaled model (4) to the scaled data: Here RMS M is the fitted root mean square using model (4) as the predictor and using Equation (5) to compute the residuals, whereas RMS T is the sum of the root mean squares using the central tendency (mean) of the data as the predictor and using the following equation to compute the residuals in Equation (6) for RMS T : We also computed the adjusted goodness-of-fit R A 2 by where p = 8 is the number of estimated model parameters. In general, R A 2 is smaller than R 2 because it takes into account the number of estimated parameters and penalizes the goodness of fit as p increases. A version of Equation (9) is used in multiple regression models, but in that case one uses −1 instead of −2 (Zar, 2009). Here we have two means (for gulls and eagles) instead of one mean, so we reduce the degrees of freedom by one more unit. Equations (7) and (9) represent a "generalized" coefficient of determination (Anderson-Sprecher, 1994), not the traditional value used in linear regression.

| Confidence intervals for parameters
Once a deterministic model has been fitted to population time series data, bootstrapping methods can be used to obtain confidence intervals for the estimated parameters (Dennis, Desharnais, Cushing, Henson, & Costantino, 2001;Falck, Bjornstad, & Stenseth, 1995).
We randomly sampled, with replacement, from the model residuals 1 , 2 , ⋯ , n g and 1 , 2 , ⋯ , n e to create sets of surrogate residuals (9) R A 2 = 1 − n g + n e − 2 n g + n e − p − 2 1 − R 2 , F I G U R E 5 Scatter plot of the 2,000 estimates of r versus α with the coexistence cutoff condition r = αC for the point estimate of C appearing as a dashed line. Estimates below that line lead to gull colony extinction. The dotted lines are r = αC using the lower and upper 95% CI bounds for parameter C. The area between the dotted lines could be considered an "uncertainty parameter region" for coexistence * 1 , * 2 , ⋯ , * n g and * 1 , * 2 , ⋯ , * n e . The time order of the residuals was ignored when sampling (Dennis et al., 2001;Falck et al., 1995).
The surrogate residuals were used to create surrogate data: For each surrogate data set, we estimated point parameters, using the method explained in section 3.6. This process was repeated n S = 2,000 times using an independent random sampling of the original residuals for each iteration. If the fminsearch algorithm did not converge to a solution within 1,000 functional evaluations, these steps were repeated for a new set of surrogate data. This occurred at a rate of 21.6%. The lack of convergence in some of the bootstrap realizations was due to the fact that the time series are relatively short and, consequently, the overall number of residuals is small, frequently leading to sets of resampled residuals that negatively impact the rate of convergence for the minimization algorithm.
We independently repeated the analyses several times with only trivial variations in the results to verify that 2,000 repetitions were adequate and that the nonconvergent bootstrap realizations were not a problem.
This procedure yielded a set of bootstrapped parameter estimates * 1 , * 2 , ⋯ , * n s that should reflect the variation one would see in the best fit parameters assuming the model (4) is valid, and the observed residuals from the model are random effects with no autocorrelation or cross-correlation.
The 95% confidence intervals for the point parameter estimates were obtained by ranking the parameter estimates for the surrogate data sets and computing the 2.5th and 97.5th percentiles (Dennis et al., 2001). For the fitted point estimates of the parameters ( Washington are predicted to equilibrate at 823 territories, which is equal to the predicted carrying capacity C ( Table 2). The goodness-of-fits indicated that the model explains at least 77% of the variability in the data (R 2 = 0.819 and R A 2 = 0.772). At the predicted equilibrium of 823 eagle territories, Equation (11) Figure A2.

| R E S U LT S
A scatterplot matrix of the bootstrapped parameter estimates, excluding β, which was always close to zero, is shown in Figure A3.
The diagonal plots are histograms showing the distribution of the 2,000 estimates for each parameter. None of these histograms suggest unusual properties for the distributions such as high skewness or multiple modes. The off-diagonal scatter plots show the pairwise relationships between the parameters for the 2,000 estimated parameter vectors. These plots can reveal strong or unusual dependencies between parameter estimates, as is the case for parameters r and α. This suggests that large estimates of gull population growth rates coincide with large estimates of gull predation rates, and vice versa.
In 13.4% of cases, the vector of parameter estimates predicted gull colony extinction. Predicted equilibria were derived for each of the 2,000 bootstrapped parameter vectors, and plotted for both species (Figure 4). The scatter plot of the 2,000 estimates shows r versus α with the coexistence cutoff condition r = αC for the point estimate of C appearing as a dashed line ( Figure 5). Estimates below that line lead to gull colony extinction. The dotted lines are r = αC using the lower and upper 95% CI bounds for parameter C. The area between the dotted lines could be considered an "uncertainty parameter region" for coexistence.
Although our deterministic model cannot predict a time to extinction, we can compute the amount of time it takes gull numbers to fall below a threshold in those cases (267 of 2,000) for which the bootstrapped parameter vectors predict extinction. If we define the threshold as 10% of the estimated carrying capacity for gulls, the estimated mean for the year in which extinction occurs is 2039 with a 95% confidence interval of (2022, 2059). For a threshold of 5% of K, the mean is 2073 with a 95% confidence interval of (2040, 2122).

| D I S C U S S I O N
We have demonstrated a strong dynamic relationship between the bald eagle population in Washington State and numbers of glaucous-winged gull nests on Protection Island's Violet Point colony. (10) (11) ln PI Eagles = 0.5074 + 0.003701 WA Eagle Territories This relationship exhibits a Lotka-Volterra-type dynamic that, at the point estimates of the parameters, predicts long-term coexistence and equilibrium for the two species. The model does not predict predator-prey oscillations as depicted originally by Lotka (1925Lotka ( , 1932 and Volterra (1926). It is notable, however, that the model predicts gull colony extinction for some parameters within the 95% confidence intervals about the point estimates.
The carrying capacity estimate for Washington eagle occupied territories (823)  Eagles have nested on Protection Island since at least the 1920s (Cowles & Hayward, 2008), and one or two eaglets were raised from a nest located on the island during many years since the early 1980 s (Hayward et al., 2010). The impact of bald eagles on seabirds has become increasingly apparent as populations have recovered and some marine fish populations on which they feed have declined (Anderson, Bower, Nysewander, Evenson, & Lovvorn, 2009;Anderson, Lovvorn, Esler, Boyd, & Stick, 2009;Stick & Lindquist, 2009;Therriault, Hay, & Schweigert, 2009). Along the west coast of North America, bald eagles have been implicated as being responsible for declines in local populations of common murres (Uria aalge; Parrish et al., 2001;Hipfner, Morrison, & Darvill, 2011), double-crested cormorants (Phalacrocorax auritus, Chatwin, Mather, & Giesbrecht, 2002;Harris, Wilson, & Elliott, 2005), pelagic cormorants (P. pelagicus, Chatwin et al., 2002;Harris et al., 2005;Carter, Hebert, & Clarkson, 2009), great blue herons (Ardea herodias; Vennesland & Butler, 2004), western grebes (Aechmophorus occidentalis; Bower, 2009), and glaucouswinged gulls (Hayward et al., 2010;Sullivan et al., 2002). The declines have been dramatic in some places. For example, a 131-km stretch along the coast of Oregon formerly supported more than 380,000 breeding pairs of common murres, but successful reproduction by these birds today is virtually nonexistent. Entire colonies have been abandoned. Murres that remain in colonies harassed by bald eagles typically give up on the breeding process before completing the nesting season (Hipfner et al., 2012). Similar effects on seabirds have been noted in Northern Europe where populations of white-tailed eagles (H. albicilla) have rebounded from declines (Hipfner et al., 2012). Hayward, unpublished data). Rhinoceros auklets (Cerorhinca monocerata) remain abundant breeders on Protection Island (Pearson, Hodum, Good, Schrimpf, & Knapp, 2013), although in 2001 they comprised the most common remains beneath an active bald eagle nest on the island (Hayward et al., 2010). Our informal observations suggest that auklets are preyed upon mostly during predawn hours when auklets leave their nest burrows to forage.
Although the model predicts coexistence at the point parameter estimates, the 95% confidence intervals include the possibility of extinction. Of the 2,000 bootstrapped parameter vectors, extinction was predicted in 13.4% of the cases. When we limit our prediction to the 1,788 cases in which the elements of the bootstrapped parameter triplets (r,α,C) are inside their 95% confidence intervals, then in 11.0% of the cases gull colony extinction was predicted. Thus, for some parameter values within the 95% confidence intervals, our model predicts that the Violet Point gull colony on Protection Island will disappear. If we define the threshold for extinction as 10% of the estimated carrying capacity for gulls, the estimated mean for the year in which extinction occurs is 2039 with a 95% confidence interval of (2022,2059). For a threshold of 5%, the mean is 2073 with a 95% confidence interval of (2040, 2122). In fact, extinction did occur on Colville Island, located 33 km north of Protection Island.
Bald eagles were known to disturb and prey on these gulls during the 1970 s , although it is unknown whether this was the cause of colony abandonment.
Extinction of the Protection Island gull colony would impact the local ecosystem in a variety of ways. The effects on vegetation would be pronounced. Gulls physically alter vegetation in their colonies through trampling, digging of nest scrapes, collection of nest material, and disturbance during boundary disputes; they chemically alter the soil through defecation and regurgitation of nondigestible components of food (Ellis, Fariña, & Witman, 2006;Lindborg, Ledbetter, Walat, & Moffett, 2012;Sobey & Kenworthy, 1979); and the decomposition of adult and juvenile gull carcasses on breeding colonies contributes nutrients to the soil (Emslie & Messenger, 1991;Lord & Burger, 1984). Thus, extinction of the gull colony would result in significant changes in the vegetation and in organisms that depend on that vegetation (Sobey & Kenworthy, 1979), and in the loss of a nutrient subsidy to the waters surrounding Protection Island (Hutchinson, 1950;Leentvaar, 1967;McColl & Burger, 1976).
Extinction also would eliminate a significant local food source for bald eagles, although gulls are not their only island food (Hayward et al., 2010), and bald eagles may gradually relocate in response to dwindling historic sources of prey (McClelland et al., 1994).
Factors other than the direct and indirect effects of eagles also may impact gull populations. Gulls are scavengers and gull populations increased dramatically during most of the twentieth century Duhem, Roche, Vidal, & Tatoni, 2008;Kadlec & Drury, 1968;Sullivan et al., 2002). Closure of landfills in this and other regions worldwide has been associated with sharply reduced gull populations (Payo-Payo et al., 2015). Indeed, the 1992 closure of the Coupeville Landfill (Anonymous, 2001), a popular feeding site for gulls located 19 km northeast of the Violet Point gull colony (Schmidt, 1986), was followed by a 10-year decline in gull nest counts on Violet Point (Figure 3).
Declines in forage fish populations (Blight et al., 2015;McKechnie et al., 2014;Therriault et al., 2009) may have played a role in gull declines. Gulls nest only along the edges of dune grass (Leymus mollis), so increases in cover by this plant could impact the size of the colony.
Dune grass cover increased from 2.5 ha (14% of Violet Point) in 1980 to 6.6 ha (39% of Violet Point) by 2009 (Cowles et al., 2012). Considerable area suitable for nesting, however, remains unoccupied suggesting that dune grass is not a limiting factor. Increasing sea surface temperatures (SSTs) that occur with El Niño events and climate change have been implicated as a factor that increases gull egg cannibalism and decreases gull colony reproductive output . We do not know whether the effects of increasing SSTs interact in some way with eagle effects on the gull population. Although these various confounding factors, which are not included explicitly in our model, undoubtedly contributed to the decline in gull numbers, it is important to note that the gull dynamics nevertheless are well predicted by the model. This  Table A1) are nearly identical to those in Table 2.
Another issue with the parameter estimation and bootstrapping methods is a lack of independence of the model residuals. This occurs because we are directly fitting the model to time series data using OLS. An alternate approach to parameter estimation is to take advantage of transitions between consecutive model states using conditional least squares (CLS). For this method, one uses parameter values and the observed population numbers at a given time to predict the population values in the following time interval (see, e.g., Dennis, Desharnais, Cushing, & Costantino, 1995). One repeats this procedure for all time steps. The best set of parameter estimates are ones that minimize the sum of the squared deviations between the observed numbers and the one-step predictions. This approach takes advantage of the Markov assumption implicit in the ordinary differential equation model: Future states depend only on the present, not the past. Since the data used for parameter estimation are conditional one-step transitions, the "observations" are, by assumption, independent. If the population data are not spaced evenly in time, however, then one cannot assume that the random one-step deviations are identically distributed, since the variance will, in general, depend on the length of the time step. Also, if one does not have observations for all the state space variables at the same times, then one-step predictions are not possible. Given that both situations exist for the gull and eagle data, we were unable to use the CLS method.
More complex methods, such as Bayesian state space modeling and Gibbs sampling, might mitigate these data deficiencies (Clark & Bjørnstad, 2004). Nevertheless, our OLS and bootstrapping procedures provided parameter estimates with a good visual fit to the data and model residuals that showed no evidence of autocorrelation or deviations from normality. The techniques we used are not specific to the predator-prey system we analyzed and could prove useful in other situations where the ecological data provide similar challenges.
Our parameter estimation and model predictions for gulls on Protection Island are based on statewide eagle data. This is because these data are more frequent and consistent over time than the eagle observations for Protection Island. However, the Poisson regression can be used to predict eagle numbers on the island based on the statewide data. We repeated our parameter estimation procedure using the predicted island eagle numbers from the nonlinear Equation (4). The estimated parameter values appear in the Table   A1. The overall fit was slightly poorer with generalized coefficients of determination of R 2 = 0.789 and R A 2 = 0.735. A graphical comparison of the model fits for the regression-predicted eagle numbers versus statewide data ( Figure A4) suggests that the former underestimates observed gull numbers for the time period 1984-1997.
A graphical analysis of the residuals supports this observation and suggests some consistent departures from normality ( Figure A5).
Moreover, the analysis for the regression-predicted eagle numbers does not take into account the error associated with the regression predictions, which, based on Figure A1, could be substantial. For these reasons, we feel that the parameter estimation and model analyses based on the statewide eagle data are more reliable.

| CON CLUS IONS
We have shown that the dynamics of a glaucous-winged gull colony on Protection Island from 1980-2016 can be explained by the number of occupied bald eagle territories in Washington with generalized R 2 = 0.82. This supports the hypothesis that the rise and decline in gull numbers observed on Protection Island are due largely to the decline and recovery of the bald eagle population. We also have shown that, with 95% confidence, the long-term dynamic predictions include coexistence but also the possibility that the gull colony will disappear, as occurred on Colville Island.
This study serves as a reminder that the necessary and successful management of one species can have direct and dramatic effects on other species; and it illustrates the uncertainty of those effects. It serves as a cautionary exploration of the future, not only for gulls on Protection Island, but for other seabirds in the Salish Sea. In particular, managers should monitor the numbers of nests in seabird colonies as well as the eagle activity within the colonies to document trends that may lead to colony extinction.

ACK N OWLED G M ENTS
We thank Jennifer Brown-Scott, Lorenz Sollmann, and Sue Thomas,

CO N FLI C T O F I NTE R E S T
None declared. All authors contributed to the writing of the manuscript. All authors gave final approval for publication.

DATA ACCE SS I B I LIT Y
The data are available in Table 1  Thus, the coexistence equilibrium is positive if and only if the inherent growth rate r of the gull population is greater than the rate at which gulls can be taken by C eagle pairs.
F I G U R E A 1 Poisson regression relationship between number of bald eagles on Protection Island and occupied bald eagle territories in Washington State, from Equation (11) in the main text

S TA B I LIT Y O F TH E EQ U I LI B R I A
The stability of each equilibrium pair is determined by the process of linearization. A comprehensive overview of linearization and stability analysis is found in Henson (2012). Linearization involves computing the Jacobian matrix at the equilibrium pair and finding its eigenvalues. If both eigenvalues of Jacobian matrix are negative real numbers, then the equilibrium is called a stable node and nearby solutions approach it in a nonoscillatory fashion. If the eigenvalues are a complex conjugate pair with negative real part, then the equilibrium is called a stable spiral, and nearby solutions approach it in an oscillatory fashion. If at least one of the eigenvalues is positive or if the eigenvalues are a complex conjugate pair with positive real part, then the equilibrium is unstable.
At equilibrium (0, C), the Jacobian has one eigenvalue that is always negative, λ = −s, and one eigenvalue that can be either positive or negative, λ = r − αC. If condition (A1) holds, the second eigenvalue is positive, so the equilibrium (0, C) is unstable. Note that the coexistence equilibrium solution (Ḡ,Ē) is positive only if (0, C) is unstable.
At the coexistence equilibrium (Ḡ,Ē), the eigenvalues are If (A1) holds, then both eigenvalues are negative (if real) or have negative real part (if complex). Thus, if the coexistence equilibrium (Ḡ,Ē) is positive, then it is stable. Whether it is a stable node or stable spiral depends on the parameter values.

S U M M A RY
If r > αC, then all four equilibrium pairs are biologically feasible. The equilibria (0, 0), (K, 0), and (0, C) are unstable, and the coexistence equilibrium (Ḡ,Ē) is positive and stable. The system will approach coexistence either with or without damped oscillations, depending on whether the eigenvalues are complex or real.
If r < αC, the coexistence equilibrium is not biologically feasible. The equilibria (0, 0) and (K, 0) are unstable and the equilibrium (0, C) is stable. In this case, the gull population will go extinct and the eagle population will approach its carrying capacity.

APPENDIX B Graphical Analyses of Model Residuals
We conducted a graphical analysis of the model residuals defined in Equation (5). These residuals are the log-scale deviations from the predicted gull and eagle numbers, obtained using the model (4) with the point estimates of the parameters and the observed values after dividing by their respective standard deviations. The first two panels of Figure A2 show the residuals plotted as a function of time. There is no evidence TA B L E A 1 Point parameter estimates and coefficients of determination for the three methods of parameter estimation (PI = Protection Island).
F I G U R E A 3 Scatterplot matrix of the 2,000 bootstrapped parameter estimates, excluding β, which was always close to zero of systematic departures from zero, although the residuals for gulls are smaller in the earlier years. The last two panels of Figure A2 show normal quantile-quantile plots of the gull and eagle residuals. The dashed line is the normal distribution expectation, and the solid line is the reference line connecting the first and third quartiles. There is no evidence in these plots of large systematic deviations from normality.

APPENDIX C Parameter Estimation for Decoupled Equations
For the full model (2) in the main text, the parameter β is the conversion rate of gulls into eagle births. Our estimate of β = 1.876 × 10 −23 is effectively zero. This might be expected since we are using statewide data for eagles, while gull population numbers are local to Protection Island. If we assume a priori that β =0, then we can decouple the equations for eagles from the equation for gulls and estimate the parameters for the eagle population separately. The predicted densities for the eagle population can then be used to find best fitting parameter estimates for gulls. Two advantages of this approach are (a) the population densities do not need to be scaled by the standard deviations to arrive at commensurate numbers for eagles and gulls, and (b) there is one less parameter to be estimated.

F I G U R E A 4
Time series plots of observed (circles) and model-predicted numbers for (a) glaucous-winged gulls and (b) bald eagles based on the regressionpredicted eagle numbers for Protection Island. The dashed curve in panel A is the prediction for gulls based on the statewide data F I G U R E A 5 Time series plots of model residuals for (a) glaucous-winged gulls and (b) bald eagles based on the regression-predicted eagle numbers for Protection Island. Normal quantilequantile plots of model residuals for (c) glaucous-winged gulls and (d) bald eagles. Compare to Figure A2