Predator-prey feedback in a gyrfalcon-ptarmigan system?

Specialist predators with oscillating dynamics are often strongly affected by the population dynamics of their prey, yet they do not always participate in a predator-prey cycle. Only those that exert strong population regulation of their prey do so. Inferring the strength and direction of the predator-prey coupling from time series therefore requires contrasting models with top-down versus bottom-up predator-prey dynamics. We examine such population-level coupling using multivariate autoregressive models. The models translate several hypotheses for the joint dynamics of population densities of the Icelandic gyrfalcon Falco rusticolus, and its prey, the rock ptarmigan Lagopus muta. The dynamics of both species are likely not only linked to each other but also to stochastic weather variables acting as confounding factors on the joint dynamics. The classical MAR(1) model, used most often in ecology, predicts that the times series exhibit predator-prey feedback (i.e., Granger causality): the predator helps to explain prey dynamics and the prey helps to explain predator dynamics. Weather, in the form of spring temperature, influences gyrfalcon population growth but not ptarmigan population growth, despite individual-level evidence that ptarmigan chicks can be strongly affected by weather. MAR(2) models, allowing for species to cycle independently from each other, further suggests alternative scenarios where a cyclic prey influence its predator but not the other way around; such bottom-up models produce a better fit but less realistic cross-correlation patterns. Simulations of MAR(1) and MAR(2) models further demonstrate that the top-down MAR(1) models are most likely to be misidentified as bottom-up dynamics than vice-versa. We therefore conclude that predator-prey feedback in the gyrfalcon-ptarmigan system is very likely, though bottom-up dynamics cannot be excluded with certainty. We finally discuss what sort of information is needed to advance the characterization of joint predator-prey dynamics in birds and other vertebrates.


Introduction
. Rather, real predator-prey systems are constantly buffeted by 28 outside forces, be those climatic or biotic variables unaccounted for (i.e, other players in the 29 interaction web). The study of Vik et al. (2008) reports at best around 55% of prey variance 30 in log-densities explained by both prey and predator densities; this therefore leaves ample 31 room for other factors to influence hare dynamics (Barraquand et al., 2017; see also Vucetich 32 et al. 2011 on ungulate-wolf systems). In birds, contrasted feedback structures (bottom-up 33 or top-down) were found in goshawk (Accipiter gentilis) -grouse dynamics, depending on 34 the grouse species considered (Tornberg et al., 2013), with marked effects of weather forces. 35 To gain a better appraisal of the strength of top-down regulation in the field, compared 36 to other drivers of herbivore dynamics (see Sinclair, 2003, for a discussion in mammals), 37 the list of predator-prey systems to which stochastic models of interacting populations are 38 fitted to time series needs to increase. Having a number of reference predator-prey systems, 39 whose dynamical structure (e.g., top-down vs. bottom-up) have been vetted by time series 40 analysis, will also help evaluating how future food web models should be structured to obtain

Statistical models
Multivariate AutoRegressive (MAR) models have been used to assess the strength of predator-81 prey coupling (Ives et al., 2003;Vik et al., 2008). Let us denote the ln-transformed predator 82 density p t = ln(P t ) and n t = ln(N t ) the ln-tranformed density of the prey; the log transfor-83 mation is useful to transform log-normal into Gaussian noise. These ln-densities are then 84 centered and stacked into a vector x t = (x 1t , x 2t ) = (n t , p t ) . The dynamics of the MAR(1) 85 model, with one timelag, are then written as a forced recurrence equation (eq. 1), 86 x t+1 = Bx t + Cu t + e t , e t ∼ N 2 (0, Σ) where B is an interaction matrix that characterizes the effects of net interactions on popu-87 lation growth of the 2 species, C describes the effect of environmental covariates u t on the 88 population growth rates of predator and prey, and e t is a Gaussian bivariate noise term. 89 We considered both a model without interactions where B =  in https://github.com/fbarraquand/GyrfalconPtarmigan_MAR) in suggested an optimal 110 lag order of 2 (BIC, HQ) or 3 (AIC, FPE). Because 2 lags are enough to model indepen-111 dently cycling populations of period up to 10 years and more (Royama, 1992), and MAR (2) 112 models are already parameter-rich, we considered a maximum of 2 lags in MAR models. The

113
MAR(2) model can be written as The independent cycling model has diagonal matrices B (1) and B (2) . The full model To model also an 116 asymmetric and nonreciprocal effect from the cyclic prey to its predator, we used the following

121
Models without environmental covariates

122
The predator-prey time series and the MAR(1) model one-step ahead predictions are pre-123 sented in Fig. 1 There is a time-lag l P for the ptarmigan (0 or 1 year) and l G = 5 (always) for the gyrfalcon: weather is expected to have such delayed effects on the gyrfalcon counts because of age structure. April weather is considered for gyrfalcon as it is the critical period for reproduction, and it is always included in models from row 3 and below. Models from rows 3 to 7 considered May temperature for ptarmigan, log(precipitation), or both. The models of rows 8 and 9 considered instead July and June temperatures as environmental variables for ptarmigan.
Based on the comparison of AICc and BIC between the full 2 × 2 interaction matrix and the diagonal matrix model (null model), the full model was favored (Table 2). The 131 model with environmental covariates did not lead to substantially better fit or very different 132 biotic interaction parameters (Table 2) than the full model, and the weather effects were not 133 consistent (Table 3), save for those of delayed April temperature on predator growth (see 134 below).

135
Additional Granger causality testing using the MAR(1) model revealed a two-way recip-136 rocal feedback, though the Wald test was weakly statistically significant (at the 0.1 level) 137 due to the low sample size (i.e., not by ecological standards but compared to other fields The addition of environmental covariates did not improve significantly model fit (Table 2).

143
The coefficients were mostly non-significant, as illustrated by the model including both 144 temperature and log(precipitation) (0 is included within CIs for environmental C matrix 145 coefficients, Table 3). The model with both precipitation and temperature was deemed 146 over-parameterized by the information criteria. It is likely that an effect of 5-year delayed 147 temperature on predator growth is present as this effect was found positive, large and nearly 148 statistically significant at 95% (Table 3). Although this weather effect on the predator does 149 not seem to improve significantly the predictive ability of the model.  Table 3: Species 1 is ptarmigan and species 2 gyrfalcon. May variables only affect species 1 while April variables, delayed by 5 years (we model the effect of variables at t − 4 on growth between t and t + 1), affect only species 2's population growth.
possible that such a negative effect is present, but it does not appear clearly with the current 156 dataset.

157
We also fitted models where winter weather affects ptarmigan growth (Appendix S2), to 158 test the idea that the survival of first-year chick might be lower in harsher winters, but again 159 we did not find consistent effects of weather on ptarmigan growth.

161
The MAR(2) models showed uniformly better fit than the MAR(1) models ( populations (i.e., diagonal interaction matrices B (1) and B (2) ) and the bottom-up predator-166 prey model (see Methods), assuming an independently cycling prey and a predator whose 167 dynamics is forced by its prey, were the better-ranking models (  (2) null + temperature uses April temperature with a 5 year delay, which affects the predator only -this model adds temperature to the list of potential drivers for predator dynamics, as it was found marginally significant in previous MAR(1) analyses. (Fig. 2). The examination of time series plots is however difficult because the simulated time 171 series are relatively short and noisy. 172 We therefore simulated 100 datasets using the fitted models (Fig. 3) and examined their  affect the growth of the adult segment of the population. However, prey population growth 187 was not affected by any of the weather covariates considered, neither in spring nor in winter.

188
This is a surprising find, which we discuss below.

189
If we had stopped the analyses to the MAR(1) model, which is customary in ecology (e.g.  the optimal lag order was 2; 3 according to AIC). MAR(2) models were therefore found to 196 realize a better trade-off between parsimony and fit than MAR(1) models (Table 4,   We found that for n = 35, a simulated MAR(1) F was recovered respectively 52%, 56%, 231 64% for AIC, AICc, and BIC, respectively. By contrast, a simulated MAR(2) BU model 232 was recovered 98%, 98% and 97%. With the current length of our dataset, the important 233 message from these simulations is that based on AIC or BIC, we are much more likely to 234 mistake a fully interacting predator-prey system for a bottom-up system than the reverse.

235
These percentages were all very close to 90% and 100% in the case of n = 100 (Appendix 236 S3).

237
From our simulation experiments, we can derive three lessons. First, from an ecological 238 viewpoint, given that the full interaction MAR(1) model both predicts the cross-correlation 239 pattern better and is the most likely to be misidentified as MAR(2) BU, we should not give 240 too much weight to the better (lower) AIC and BIC scores of the MAR(2) BU model. It is 241 fairly likely that top-down prey regulation and therefore reciprocal predator-prey feedback 242 is at work here. Second, from a more statistical viewpoint, it is informative to notice that 243 whether MAR(1) or MAR(2) models are better identified is parameter-specific: sometimes 244 a MAR(1) model will be more likely to be correctly classified (as in the simulation study 245 of Lütkepohl, 2005), sometimes a MAR(2) will (our case study). Therefore, whether an 246 ecological scenario is better identified than another one from time series is likely to be 247 context-specific, and not simply dependent on the lag order. The corollary being that new

253
Given that the data presented here is collected once a year for the most part, and that 254 it is not feasible to census the population much more frequently with current means (other 255 technologies would be necessary, such as camera traps or DNA-based evidence), it is un-256 likely that we will get the time series near 100 years within acceptable time frames for 257 management of both populations (i.e., conservation of the gyrfalcon and sustainable hunting 258 management of ptarmigan). Therefore, differentiating unequivocally between the bottom-up 259 and predator-prey feedback scenarios will likely require other type of models and data. We 260 still view the MAR(p) approach as useful, however, as a means to delineate likely scenarios 261 to investigate further, and check for important abiotic drivers that need to be considered.

262
Why the ptarmigan population growth is not affected -at the NE Iceland scale -by weather 263 variables is also a puzzling question for further study.

264
Mechanistic modeling might help further understand the effect of drivers on ptarmigan 265 dynamics during certain phases of the cycle. For instance, are ptarmigan declines mainly 266 driven by its predator (i.e., the gyrfalcon is responsible in large part for the declines) or mainly by other causes such as parasites? Rough estimates of predation demonstrate why 268 this question is intrinsically difficult. Around 100 adult predator pairs can be found near peak 269 abundance on the NE Iceland ptarmigan management zone, and to these correspond about 270 100 000 ptarmigan individuals at best (Sturludóttir, 2015). One might think, given these 271 numbers, that the predators are unlikely to make their prey decline. Ptarmigan, however, 272 have a slow, long-period cycle (Fig. 1). Therefore, they decrease at worst by ≈ 20 000 in a Iceland, using averaged state variables and covariables, and yet that weather locally affects 294 the survival of ptarmigan chicks, as additional data seem to suggest: Nielsen et al. (2004) 295 found that mean windspeed and mean precipitation in June-July explained a considerable 296 part of the variance in chick production.

298
Our results have implications for other studies on birds and more generally vertebrates 299 with relatively slow life histories (compared to e.g., plankton sampled many times a year).

300
Using long time series by ecological standards (34 years), we found some evidence of re-301 ciprocal predator-prey feedback in this cyclic predator-prey system, without being able to 302 exclude nonetheless more bottom-up predator-prey dynamics. MAR(p) models with p = 1, 2 303 described well this system as a forced oscillator, although the unexplained noise was gener-304 ally stronger than weather effects, which may point to other important biotic factors driving 305 the dynamics, such as parasites. Simulations of the fitted models revealed than unequivocal 306 inference of bottom-up versus reciprocal predator-prey coupling (i.e., including top-down 307 predator influence on prey) would require about a century of time series data. We therefore 308 think that additional demographic data (e.g., through capture-recapture, genetics,...) should 309 always be considered in conjunction to counts taken once or twice a year, if one of the goals 310 of a monitoring study is to infer interactions between the populations of different species.   We also considered minimum temperature but this did not alter the following results.

407
The two above mentioned winter weather variables were inserted in place of spring 408 weather variables for ptarmigan into a MAR(1) model. The estimated parameters are re-409 produced in Table S1 and the Information Criteria, with previous models for comparison, in 410   Table S2. None of the models are able to significantly improve the fit, although it is possible 411 that a weakly statistically significant effect of winter temperature exists.  Table S1: Coefficients for biotic and abiotic effects on population growth. Species 1 is ptarmigan and species 2 gyrfalcon. Winter variables only affect species 1 while April variables, delayed by 5 years (we model the effect of variables at t − 4 on growth between t and t + 1), affect only species 2's population growth. Effects of winter variables are depicted in italics.  Table S2: Comparison of model selection criteria for MAR(1) models. MAR(1) 'null' indicates a diagonal B matrix while MAR(1) 'full' indicates a full 2x2 interaction matrix. Models including temperature (third row and below) effects on growth rates take the form x t+1 = a + Bx t + Cu t + e t , e t ∼ N 2 (0, Σ). Here the environmental vector u t = (T t−l P , R t−l P , T t−l G , R t−l G ) T , with T the temperature and R log-rainfall. There is a timelag l P for the ptarmigan (0 or 1 year, depending on the month) and l G = 5 for the gyrfalcon. IC scores for the two winter models are depicted on the last two rows . prec. = precipitation (rain and snow).  Table S1: Frequency of correct identification of MAR(1) full and MAR(2) bottom-up models for two time series lengths.