Summary

1. Although recent decades have seen much development of statistical methods to estimate demographical parameters such as reproduction, and survival and migration probabilities, the focus is usually the estimation of parameters for individual species. This is despite the fact that several species may live in close proximity, sometimes competing for the same resources. There is therefore a great need for new methods that enable a better integration of demographical data, e.g. the study of synchrony between sympatric species, which are subject to common environmental stochasticity and potentially similar biotic interactions.

2. We propose a mark–recapture statistical model that uses random effect terms for studying synchrony in a demographical parameter at a multi-species level, adapting a framework initially developed to study multi-site synchrony to this novel situation. The model allows us to divide between-year variance in a demographical parameter into a ‘synchronous’ component, common to all species considered, and species-specific ‘asynchronous’ components, as well as to estimate the proportion of each component accounted for by environmental covariates.

3. We demonstrate the method with data from three colonially breeding auk species that share resources during the breeding season at the Isle of May, Scotland. Mark-resight information has been collected since 1984 for Atlantic puffins Fratercula arctica, common guillemots Uria aalge and razorbills Alca torda marked as breeding adults. We explore the relationship between synchrony in the species’ survival and two environmental covariates.

4. Most of the between-year variation was synchronous to the three species, and the same environmental covariates acted simultaneously as synchronising and desynchronising agents of adult survival, possibly through different indirect causation paths.

5.Synthesis and applications. The model proposed allows the investigation of multi-species synchrony and asynchrony in adult survival, as well as the role of environmental covariates in generating them. It provides insight into whether sympatric species respond similarly or differently to changes in their environment, and helps to disentangle the sources of these differences. The estimated indices of synchrony/asynchrony can facilitate the generation of further hypotheses about similarities/differences in these species’ ecology, such as the potential overlap of wintering areas. The method is readily applicable to other species, ecosystems and demographical parameters.

Introduction

The monitoring of demographical parameters is generating a wealth of valuable information for ecology and conservation and recent decades have seen a corresponding proliferation of statistical models for analysing these types of data. However, the potential to integrate different types of data has not been fully exploited, with the majority of these models targeted at analysing single demographical parameters for individual species (Lebreton et al. 1992; Williams, Nichols & Conroy 2002), although some approaches such as integrated population modelling (Besbeas, Freeman & Morgan 2005) jointly estimate several parameters in a single-species analysis. Data from several species have recently been combined in models to study population trends (Sauer & Link 2002) or species richness (multi-species occupancy models, Russell et al. 2009) but a move away from single-species single-location to more encompassing approaches is still largely overdue.

Species exist within the context of communities and ecosystems, and when populations of different species are sympatric they are exposed to biotic interactions and the same abiotic environment (Begon, Townsend & Harper 2006). Some species may react in a similar way to their common environment, showing synchrony in population trends or in the temporal variation of some demographical parameters such as survival. The underlying cause of synchrony between species is usually not clear, with hypotheses suggesting shared stochastic effects, such as weather (Hawkins & Holyoak 1998) and the response to common predators (Raimondo et al. 2004). The study of the species in a community in isolation may lead to only a partial understanding of their ecology or even to incorrect conclusions.

Synchrony between sympatric populations of different species has received less attention compared with synchrony between allopatric populations of a single species (Raimondo et al. 2004). The relatively few multi-species examples to date typically address synchrony in abundance through the study of time series of population size (Swanson & Johnson 1999; Raimondo et al. 2004) and are often dedicated to understanding mechanistic predator–prey interactions (New 2009). In general, investigating the mechanisms underlying population change is a difficult task when studying time series of abundance alone, and the incorporation of demographical parameters such as survival, reproductive success and dispersal probabilities is often key in understanding such mechanisms (Loison et al. 2002).

We propose a statistical model for studying the variation of a demographical parameter at a multi-species level, through the use of random effects. Between-year variance in the demographical parameter is divided into a ‘synchronous’ component, which represents the common response of all species considered, and ‘asynchronous’ components, specific to each species, and we estimate the contribution of environmental covariates to each of these components. The model is based on that presented by Grosbois et al. (2009) for studying the variation of adult survival for a single species at a multi-population scale, although it is conceptually different in its interpretation, and further we relax the variance structure in the model to accommodate differences among species. In this article, we demonstrate an application of the method to explore synchrony in adult survival using 25 years of individual mark-resight data for three seabird species, the Atlantic puffin Fratercula arctica (L.), the common guillemot Uria aalge (Pontoppidan) and the razorbill Alca torda L., collected at the breeding colony on the Isle of May, southeast Scotland. These three auks (Alcidae) have broadly similar life histories and ecology (Gaston & Jones 1998). Birds from breeding populations on the Isle of May show largely overlapping distributions throughout the year (Wernham et al. 2002) and are thus likely to be exposed to similar environmental stochasticity. Consequently, we would expect some degree of synchrony in their response in terms of the temporal variation of demographical parameters. Adult survival probabilities for Isle of May puffins, guillemots and razorbills have previously been analysed separately (Harris et al. 1997; Harris, Wanless & Rothery 2000; Crespin et al. 2006), but to date no attempt has been made to integrate survival data for these species, and in particular, to look for synchronising and desynchronising agents.

Materials and methods

Mark-resight data

Mark-resight information was collected for 543 Atlantic puffins (hereafter puffin), 831 common guillemots (hereafter guillemot) and 153 razorbills at the Isle of May (56°11′N, 2°34′W), southeast Scotland. As with many seabirds, annual adult survival probabilities of puffins, guillemots and razorbills are normally high (Harris et al. 1997; Sandvik et al. 2005). Birds visit land only for breeding and while puffins nest in burrows, guillemots and razorbills lay eggs directly on narrow cliff ledges. Outside the breeding season, auks from the Isle of May disperse over broad areas of the North Sea (Wernham et al. 2002), and during the breeding season they eat similar prey, mainly small, lipid-rich, shoaling fish such as the lesser sandeel Ammodytes marinus and sprat Sprattus sprattus (Daunt et al. 2008). Between 1984 and 2007, breeding birds in front of permanent hides were marked with unique colour rings and resightings of these birds took place each year up to 2008. Once they have bred, individuals of all species rarely breed more than a few metres from where they were marked (MPH, pers. obs.), so resighting effort was mainly focussed on these areas although regular searches were also made in all nearby areas.

Multi-species synchrony model

Using the standard open population capture–mark–recapture/resight models for estimating apparent survival and recapture/resight probabilities (reviewed in Lebreton et al. 1992), likelihood functions can be constructed individually for each of the species involved in the model. Following standard notation, we denote resight probability in year Y as p(Y) and annual apparent adult survival probability from year Y to Y + 1 as Φ(Y). Both resight and survival probabilities can then be modelled to depend on explanatory variables. To allow the study of synchrony in survival probabilities, we followed the framework presented by Grosbois et al. (2009) and introduced random year effects in the relationship of survival with covariates as follows:

(eqn 1)

ΦS(Y) is the apparent adult survival from year Y to Y + 1 for species S. Survival is related to covariates and random effects through the logit link function, although alternative link functions are possible. The relationship with covariates is handled through fS(.), a species-specific function of ns covariates cSi. The function could be, for example, a linear regression or a nonparametric relationship. For the year random effects, δ(Y) is a random term which is common to all species considered and only depends on the year, and εS(Y) is a random term that depends on the year and species. The δ and εS terms are assumed to be independent and have distributions and respectively, with no correlation between terms. We extend the approach of Grosbois et al. (2009) so that the ‘year × species’ random terms can have different variances for the different species (i.e. are species-specific) and so the between-year variance in survival unexplained by the covariates can be differently partitioned for different species. The δ term corresponds to the amount of between-year variation (unexplained by the covariates, if present) that is synchronous to all species considered, while the εS terms characterise the species-specific (asynchronous) components. Note also that Grosbois et al. (2009) use a single common covariate that takes different values for each colony, while in this study each common covariate has the same value for all the species considered (as the geographical area is the same), but each species might have a different combination of covariates.

Assuming independence between the data for the different species, the overall likelihood function for all species together can be written as the product of the individual likelihoods. This is similar to the way in which likelihood components corresponding to different demographical parameters are combined for a single species in an integrated population modelling framework (Besbeas, Freeman & Morgan 2005). In the proposed model, the species-specific likelihood components share at least the common random term.

Once the model parameters have been estimated, a species-specific intra-class correlation coefficient can be calculated based on the variances of the random terms as an index of synchrony in adult survival:

(eqn 2)

This quantity represents the synchrony of species S with the rest of the species: the amount of between-year variance for species S (either total or unexplained by the covariates, if present) that is accounted for by the common random term δ(Y). When is large compared with , then ICCS is large and the between-year variation for that species is then mostly synchronous with the other species.

To evaluate the effect of the environmental covariates in generating synchrony and asynchrony between species survival, two models are compared, following the method in Loison et al. (2002) and Grosbois et al. (2009). Both models include random effects δ(Y) and εS(Y) but one of them does not have covariates (only a separate intercept for each species). We can define the residual variance of δ, , for the model with covariates, and the total variance, when there are no covariates. The same can be done for the variance of the εS terms: , . Based on these values, a set of coefficients can be calculated:

(eqn 3)

(eqn 4)

Cδ and the Cs coefficients measure the contribution of the environmental covariates to the interspecific synchronous δ and asynchronous εs components of the between-year variances, respectively.

Environmental covariates for survival

Environmental covariates are known to influence demographical parameters in many species and have been shown in some cases to be responsible for interspecific synchrony (Hawkins & Holyoak 1998). For North Atlantic seabirds, studies of survival often include covariates related to two oceanographic factors, the North Atlantic Oscillation (NAO) and the temperature at the sea surface (SST).

Heterogeneity in resight probability

Before fitting the data in combination, we assessed the goodness-of-fit (GOF) of the general fully time-dependent Cormack–Jolly–Seber (CJS) model {Φ(t)p(t)} with program RELEASE (Burnham et al. 1987), for each species individually. In this model, both survival and resight probabilities are allowed to vary from year to year. The GOF was very similar for all species studied. The general CJS model fits the data poorly, due mostly to the 2.C component (guillemot: χ2 =173·49, d.f. = 22; puffin: χ2 =129·9, d.f. = 22; razorbill: χ2 =55·35, d.f. = 17; all P-values <0·001), which indicates heterogeneity in resight probability (trap dependence), an effect that has been reported already for puffins at the Isle of May (Harris et al. 2005). Component 3.SR fitted well for all three species (P-values >0·9), therefore showing no evidence of individual heterogeneity in survival probability, as noted in previous analysis of these species from the Isle of May (Harris, Wanless & Rothery 2000; Harris et al. 2005; Grosbois et al. 2009).

The trap dependence in resight probability detected for the three species was taken into account in the synchrony models by adding a 1-year trap dependence structure as follows:

(eqn 5)

For each species S, the resight probability for individual i in year Y, pS(i,Y), depends through a logit link on a year-specific resight probability rS(Y) and an additive term aS that is only included if the individual was resighted in the previous occasion. This is achieved by using the indicator function TS(i,Y) that can be seen as an individual covariate for each capture occasion. Thus, TS(i,Y) = 1 if bird i was caught in year Y − 1, and 0 otherwise. The species-specific terms aS represent the amount of 1-year trap dependence for each species studied.

Analysis of the auk data

We applied the method outlined above to investigate the amount of synchronisation in adult survival for the three auk species at the Isle of May. For simplicity, we used the same covariates for all three species, but this is in general not a restriction and species-specific covariates could be considered. The vector cov = {c1, c2, c3, c4} = {wNAO_0, wNAO_1, SST_0, SST_1} in the models hereafter refers to the four covariates together. All covariate time series (1984–2008) were standardised prior to inclusion in the models by subtracting the mean of the series and dividing by its standard deviation. We verified that the covariates did not have high correlation. For adult survival, we considered a logit link function and a linear regression, with the aforementioned set of four standardised covariates cj and corresponding species-specific regression coefficients βjS:

(eqn 6)

All models considered in the following sections had fully time-dependent resight probability with 1-year trap dependence modelled as explained in Eqn (5) and these are denoted ‘p(t + a)’.

As the Bayesian approach is more flexible for handling random effects than the classical maximum likelihood framework (Barry et al. 2003), we conducted our study within a Bayesian framework with markov chain monte carlo (MCMC) sampling. All models were programmed in WinBUGS (Spiegelhalter et al. 2003). The code used to fit the models can be found in Appendix S1. After a burn-in of 100 000 samples, the MCMC chains were run for 150 000 iterations (with a thinning of 3). Convergence was assessed with the Gelman–Rubin statistic calculated as modified by Brooks & Gelman (1998), after starting three chains with dispersed initial values for all variables. The statistic suggests that convergence had been achieved after 100 000 samples. Uninformative priors were used for all variables: regression coefficients βiS ∼ U(−5,5); standard deviation of the δ and ε random terms σx ∼ U(0,3); year-specific component of resight probabilities rS(Y) ∼ N(0,104); trap dependence coefficients aS ∼ U(−5,5). We conducted a prior sensitivity study for the random effect variances by specifying conventionally used vague inverse-gamma priors as an alternative to uniform priors.

Starting from the full model {Φ(cov + δ + ε)p(t + a)}, we constructed all of the eight combinations of up to three of the arguments of survival (covariates, ‘year’ random term δ, ‘year × species’ random terms εS), or none at all. In the cases when covariates were removed, a species main effect was kept through a species-specific intercept. For brevity, we did not attempt a formal model-reduction exercise in terms of reducing the number of individual covariate terms required. We fitted each of the resulting eight models (Table 1) to the auk mark-resight data. The models were ranked in terms of their Deviance Information Criterion (DIC), a Bayesian analogue of AIC (Spiegelhalter et al. 2002) that balances model fit and complexity. It is calculated as DIC = D(θ) + 2pD, where D(θ), the deviance when using the mean of the posterior distribution of the parameters, is penalised by twice the effective number of estimated parameters pD. DIC is available directly in WinBUGS, with the best model being the one with the lowest DIC value. Although its use is controversial in the context of hierarchical models (Spiegelhalter et al. 2002; Barry et al. 2003; Millar 2009), note that the model ranking does not affect the analysis of synchrony.

Table 1. DIC values for the different models compared

Model

DIC

ΔDIC

‘cov’ refers to the set of four covariates (wNAO_0, wNAO_1, SST_0 and SST_1). ‘S’ refers to species main effect (intercept only). ΔDIC is the Deviance Information Criterion (DIC) increment compared to the model with lowest DIC.

Φ(cov + δ + ε)p(t + a)

1104·2

0

Φ(S + δ + ε)p(t + a)

1105·1

0·9

Φ(cov + δ)p(t + a)

1108·3

4·1

Φ(cov + ε)p(t + a)

1110·6

6·4

Φ(S + ε)p(t + a)

1111·7

7·5

Φ(S + δ)p(t + a)

1139·8

35·6

Φ(cov)p(t + a)

1153·6

49·4

Φ(S)p(t + a)

1202·2

98·0

Simulation study

We used simulation to study the performance of the proposed method in fitting a set of data derived from known parameters. We selected the full model structure {Φ(cov + δ + ε)p(t + a)} from the previous section and chose parameter values based on the best model obtained in the Isle of May auk study, to stay within ecological realism. Mark-resight data were generated 20 times (the processing time required for the MCMC sampling is prohibitive for a much larger simulation study), matched to values estimated for the three auk species, with the same number of animals as in the real data set. For each species, each of these data sets differed only in the value of the survival probabilities, as the random effect terms (both common and species-specific) that were added to the linear relationship were generated independently with same variance for each simulated data set. The rest of the parameters were kept unchanged. The model {Φ(cov + δ + ε)p(t + a)} was fitted to the 20 data sets using WinBUGS (50 000 MCMC iterations after a burn-in of 100 000). We used the medians of the posterior distributions for each parameter to calculate the bias, and then averaged over the 20 data sets.

Results

Data analysis

According to the DIC values (Table 1), the best of the eight possible models fitted was the full model {Φ(cov + δ + ε)p(t + a)}, where survival depended on the set of covariates but had also common (δ) and species-specific (εS) random terms. As expected, the models with covariates outperformed the corresponding models with only species main effect (intercept). In both cases, with and without covariates, the inclusion of any kind of random effects gave a substantial improvement in terms of DIC, and having both common and species-specific random terms was better than having either in isolation.

Prior sensitivity was tested for the best model using alternative priors. In particular, the use of inverse-gamma priors for the random effect variances appeared to be slightly more informative than specifying uniforms for their standard deviation, and the posterior distributions were sensitive to the choice of the gamma distribution parameters, as has been noted in previous studies (Royle 2008). This was particularly the case for razorbills, the species with least data. These results support the selection of uniform priors for these parameters.

Concentrating on the full model {Φ(cov + δ + ε)p(t + a)}, estimated survival probabilities (Fig. 1) differed substantially for the three species, although most values remained relatively high, as is typical for long-lived seabirds. Note that the size of the 95% credible intervals reflected the amount of data available for each species, being wider for razorbill (153 birds) and very narrow for guillemots (831 birds). Survival was relatively stable over the years for guillemots, showed wider variation for razorbills, with pronounced peaks in a few particular years, whereas estimates for puffins were intermediate.

Figure 1.

Estimated apparent adult survival from model {Φ(cov + δ + ε)p(t + a)} for (a) puffin, (b) guillemot and (c) razorbill at the Isle of May. The point estimates are the median of the MCMC samples for each variable, obtained with WinBUGS. Vertical bars show 95% credible intervals. Survival probabilities from the fully time-dependent model {Φ(t)p(t + a)}, estimated with WinBUGS for each species separately, are shown as a dotted line.

The trap dependence coefficients (aS in Eqn 5) were all positive for the three species (Table 2) and therefore the probability of seeing a bird was higher if it was seen the previous year. Using the estimates of aS and rS, we calculated the estimated resight probabilities for the three species, for the case when a bird was seen the year before, and for when it was not (Appendix S2).

Table 2 shows the regression coefficients for the full model, {Φ(cov + δ + ε)p(t + a)}. Most of the point estimates were below zero, denoting a negative relationship between adult survival and the covariate represented. Note that some of the 95% credible intervals spanned both sides of 0. In the particular case of 1-year time-lagged SST for razorbill, the corresponding beta was very close to zero, indicating a lack of strong influence of that covariate on razorbill survival. The fact that some of the regression coefficients corresponding to the time-lagged versions of wNAO and SST were far from zero indicated that they also had an indirect effect on adult survival, acting possibly through the food chain (Harris et al. 2005; Sandvik et al. 2005).

Interspecific synchrony (ICCS) and the fraction of variation accounted for by the covariates for each species (Cδ and Cs terms) were calculated from the estimates of the full model {Φ(cov + δ + ε)p(t + a)} and the ‘species main effect’ model {Φ(S+δ+ε)p(t+a)} (Table 3).

Table 3. Estimated residual and total variance of the common (δ) and species-specific (εS) random effect terms and inter-class correlation (ICCS) coefficients

Interspecific synchronous variance component

Species-specific asynchronous variance component

Inter-class correlation ICCS

The fraction of between-year variance in survival accounted for by the climatic variables (Cδ and CS) was calculated based on the estimated variances. ‘Species 1’ refers to puffins, ‘species 2’ to guillemots and ‘species 3’ to razorbills. 95% credible intervals are shown in parentheses.

Model Φ(S + δ + ε)p(t + a) (total variances)

= 0·386 (0·066, 0·885)

= 0·191 (0·017, 0·628)

ICC1 = 0·667 (0·173, 0·965)

= 0·137 (0·008, 0·487)

ICC2 = 0·735 (0·245, 0·982)

= 0·202 (0·005, 0·849)

ICC3 = 0·665 (0·117, 0·987)

Model Φ(cov + δ + ε)p(t + a) (residual variances)

= 0·288 (0·091, 0·711)

= 0·036 (0·000, 0·346)

ICC1 = 0·894 (0·304, 0·999)

= 0·079 (0·001, 0·377)

ICC2 = 0·787 (0·350, 0·996)

= 0·082 (0·001, 0·660)

ICC3 = 0·785 (0·205, 0·998)

Fraction of variation accounted for by the climatic covariates

Cδ = 0·256

C1 = 0·810

C2 = 0·425

C3 =0·595

For the model with covariates, the residual variances of the species-specific random terms were all substantially lower than that of the common random term which is also noticeable when looking at the estimates of the random terms for each year of the study (Fig. 2). ICCS values were consequently high, which suggested that most of the variation unexplained by the environmental covariates was synchronous to the three species.

Figure 2.

Value of the random effect terms (on the logistic scale) estimated for each year by the best model {Φ(cov + δ + ε)p(t + a)}. Both common random terms δ(Y) and the species-specific random terms εS(Y) for each species are shown.

In the ‘species main effect’ model, all and variances increased compared with the model with covariates, to accommodate the extra variation created by the lack of covariates. The species-specific variances increased more, in proportion, and therefore the ICCS values decreased to below 75%.

The fraction of the synchronous variance accounted for by the set of covariates (Cδ) was around 26%, that is, about one-fourth of the variation that is synchronous to the three auk species was explained by components of the climate related to wNAO and SST. Climate is acting to some extent as a synchronising agent in the survival of puffins, guillemots and razorbills but there is still about 75% of synchronous variation that is not explained by these covariates. The environmental covariates were also responsible for a large part of the asynchronous variation, as shown by the values of the CS coefficients. For puffins and razorbills, the values were very high (≈81% and 60%, respectively), implying that most of the between-year variation asynchronous to the other auk species was related to these climatic covariates. For guillemots on the other hand, less than half of the asynchronous variation in adult survival was explained by these covariates. Thus, it appears that the same climatic factors can act simultaneously as synchronising and desynchronising agents for adult survival of these species at the Isle of May. There is some indication that both wNAO and SST can act indirectly on survival (Harris et al. 2005). It is therefore possible that the oceanographic effects reflected in wNAO and SST can act through different indirect causation paths, some of them affecting the three species in synchrony, some others affecting them differently or only affecting some of the species.

Simulations

We obtained the average over the 20 simulated data sets of the median value of each parameter (Appendix S3). Bias was calculated as the average over the 20 simulations of the absolute value of the difference between the point estimate (median) and the true value. It was generally small for the regression and trap dependence coefficients. The largest values appeared with the estimation of the variance of the random effects. In relation to the species-specific random terms, it is worth-noting that as expected the largest bias was associated with the species with least data (razorbill, 153 marked individuals), whereas the smallest corresponds to guillemots (with 831 birds). These differences disappeared when the simulations were repeated with 831 individuals for each of the species. Bias in survival estimates was in almost all cases below 3% and was again in general largest for razorbills (smallest data set) and smallest for guillemots (largest data set).

Discussion

This article presents a model, fitted using Bayesian methodology, for studying synchrony in adult survival between several species, and the contribution of environmental covariates as synchronising and desynchronising agents, adapting the framework used for a multi-population study by Grosbois et al. (2009) to the multi-species situation. This method does not directly shed light into the typically complex mechanisms that underlie the observed synchronisation or desynchronisation between different species, but it can be used to provide insight into community dynamics and to point out further avenues of investigation in terms of environmental covariates.

Auk survival at the Isle of May

The survival estimates obtained in our study with the best model {Φ(cov + δ + ε)p(t + a)} are consistent with previous analyses of the three species individually (Harris, Wanless & Rothery 2000; Harris et al. 2005). However, estimates of a species’ survival from a more integrated study have the potential for borrowing strength from the rest of the ensemble, with the consequent gain in precision. In this study, some of the regression coefficients seem to point to the existence of indirect environmental effects, possibly through the food web, as noted in Sandvik et al. (2005): regression coefficients were negative for SST with no delay and others with 1-year lag were not zero. Some of the estimated regression coefficients were low and had 95% credible intervals that included zero, possibly pointing to a lack of a strong influence of the corresponding environmental covariates on that particular species’ survival. We did not attempt a systematic covariate selection process prior to the modelling as the primary aim at this stage was to develop the statistical model for studying multi-species synchrony and demonstrate the potential of this framework.

There was a significant proportion of variance not explained by our set of covariates, which indicates that there is scope for further investigation. This may include the existing environmental covariates with longer time-lags (Harris et al. 2005) or averaged over different periods of the year or broader areas in which auks overwinter (Sandvik et al. 2005). Biotic covariates, like prey stock estimates (Harris et al. 1997), could also be considered, as well as nonlinear or nonparametric relationships with the covariates (Gimenez et al. 2006). These covariates will be the object of further exploration of this data set, with a focus on the ecology of these auk species. Our study lays the methodological groundwork for this.

Future developments of the framework

A number of interesting generalisations can be considered for the framework presented by Grosbois et al. (2009) for the multi-colony case and extended in this study for multi-species synchrony. First, the framework of using species-specific and common random effect terms could be adapted to other demographical parameters, as already suggested by Grosbois et al. (2009) for the multi-population situation. The natural next step would be to consider synchrony in several demographical parameters by analysing them together and potentially incorporating time series of abundance, in an integrated population modelling framework (Besbeas, Freeman & Morgan 2005). The joint likelihood of such analysis would extend over demographical parameters and different species, an analysis that to our knowledge has not been done to date. Apart from the inherent benefits of the integrated modelling, the partition of variation into synchronous and asynchronous components could be carried out for the different demographical parameters, in the same fashion as was done here for adult survival. The synchrony/asynchrony of the response to the environmental covariates could be simultaneously assessed across species for different life-history traits. Conversely, synchrony in different demographical rates could be studied for a single species, investigating for example if juvenile and adult survivals are synchronous and if climate contributes to this effect.

Few studies address spatial and temporal synchrony simultaneously (but see Swanson & Johnson 1999). For survival, such a situation could be tackled with a multi-species multi-population framework, combining the model proposed by Grosbois et al. (2009) with that of ours:

(eqn 7)

Survival ΦSP(Y) for species S in site P would be related to a species-and-site-specific function fSP(.) of a set of nSP environmental covariates cSPi(Y) and random effects. These would include an overall common term δ(Y), terms specific to species λS(Y) and sites γP(Y), and finally species-and-site-specific terms εSP(Y). With more parameters to be estimated compared with the multi-species or multi-colony cases, we can expect the requirements in terms of amount of data needed to be able to estimate them to increase.

The alternative parameterisation proposed as a generalisation of the multi-population model (Grosbois et al. 2009; eqn 1) can also be adopted in the multi-species framework we present, allowing the incorporation of covariates into the species-specific partition of variance between synchronous and asynchronous components. When mechanistic hypotheses about interspecific relations exist, it could be worth considering applying the framework presented here to models that take into account these interactions explicitly (see New 2009 for an example with predator–prey interaction).

Finally, the application of random effects to study multi-species synchrony could be explored for other types of data beyond mark–recapture. One example is occupancy models (MacKenzie et al. 2006) where detection/non-detection data of an unmarked species are used to estimate the percentage of sampled sites where the species is present, taking into account imperfect detection. In a similar fashion to that in Eqn (1), data from several species sampled at the same sites could be modelled together, adding common and species-specific random effects terms to account for the between-site variation not accounted for by a set of covariates:

(eqn 8)

In this case, ΨS(i) (the probability of site i being occupied by species S) depends on a set of covariates cSj(i) and two random terms. The variance of the common terms δ(i) represents the variation of occupancy across sites that is synchronous to all species considered, while the variances of the species-specific terms εS(i) correspond to the asynchronous components. The derivation of indices of synchrony and the contribution of the covariates in synchronising and desynchronising occupancy across sites is then straightforward. The number of sites is usually large compared with the number of years available in typical mark–recapture studies, facilitating the characterisation of the random effect variances. We note that multispecies occupancy models have already been proposed to study communities (e.g. MacKenzie, Bailey & Nichols 2004; Russell et al. 2009) although not specifically targeted to investigate synchrony in occupancy.

Conclusion

Improved understanding of how the environment synchronises and desynchronises demographical parameters can be of great value in generating ecological hypotheses, especially when coupled with biological knowledge of these species. Links between demography and environmental conditions are complex, with variables acting simultaneously as synchronising and desynchronising agents. For example, in the case of the auks considered here, it is likely that to understand the processes involved, more information will have to be incorporated. The results of synchrony could be related to similarities in wintering grounds, as new research clarifies the picture of where these birds spend the winter months (Harris et al. 2010). Models like the one presented by Grosbois et al. (2009) for multi-populations and its adaptation for multi-species introduced in this paper represent new steps towards more integrative approaches to study demographical parameters. Methods to study multi-species relations are urgently needed given the changing environmental conditions and may play an important role in increasing our understanding of how climate change may affect communities’ composition, as sympatric species react in similar or different ways to changes in their environment.

Acknowledgements

J.J.L.-M. was supported by a grant provided by the Centre for Ecology and Hydrology and the EPSRC National Centre for Statistical Ecology. We thank the many people who have helped with ringing and looking for colour rings on the Isle of May, particularly Mark Newell. Part of the fieldwork was carried out with funding from the Joint Nature Conservation Committee’s integrated Seabird Monitoring Programme. Scottish Natural Heritage allowed us to carry out these studies on the Isle of May National Nature Reserve. We also thank Olivier Gimenez and Marc Kéry for helpful comments and suggestions that improved the article.

Publication History

Supporting Information

Appendix S1. WinBUGS code.

Appendix S2. Resight probability estimates.

Appendix S3. Simulation study results.

As a service to our authors and readers, this journal provides support ing information supplied by the authors. Such materials may be re-organized for online delivery, but are not copy-edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Please note: Wiley-Blackwell is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.

8SJ Burthe, S Wanless, MA Newell, A Butler, F Daunt, Assessing the vulnerability of the marine bird community in the western North Sea to climate change and other anthropogenic impacts, Marine Ecology Progress Series, 2014, 507, 277CrossRef