Episodes of opposing survival and reproductive selection cause strong fluctuating selection on seasonal migration versus residence

(a) Study system and data collection

The focal shag population breeds on the Isle of May (IoM) National Nature Reserve, Scotland (56°11′ N, 2°33′ W). These shags are typically socially monogamous, rearing a single brood per year with biparental care (electronic supplementary material, S1). Since 1997, more than 17 000 chicks (approx. 80% of all those hatched) and more than 900 additional adult recruits have been ringed with uniquely coded colour rings (field-readable to 150 m with a telescope), generating a breeding population of individually marked adults. During ten breeding seasons (April–July 2009–2018, ‘summers’), virtually all nest sites on IoM were monitored through frequent, systematic checks (electronic supplementary material, S1). Colour-ringed nest owners were systematically identified and hence classified as breeders, and sexed through vocalizations and/or genotyping. For active nests (mean 533 year⁻¹, range 388–821), and hence their associated owners, breeding success was recorded as the number of chicks fledged (0–4), or recorded as unknown in cases with uncertain success (approx. 1% of nests; electronic supplementary material, S1). In addition, regular (approx. 3 week⁻¹) resighting surveys at roost sites on IoM were undertaken to identify ringed adults that apparently did not attempt breeding or could have failed early (hereafter early-failed/non-breeders). These individuals were assigned breeding success of 0 fledglings, affecting 2–18% of all resighted adult females and 1–30% of males (means 5%, electronic supplementary material, S1). Due to the intensive ringing, comprehensive nest monitoring and high overall breeding season resighting probability (mean 0.95 during 2010–2018, range 0.90–0.98 [24]), annual breeding success was assigned for a very high proportion of the total adult population. Ringing and nest monitoring were licensed by British Trust for Ornithology and NatureScot.

Because shags have partially wettable plumage and hence must return to shore every day to dry, ringed individuals can be resighted at coastal locations throughout the non-breeding season (winter). Hence, throughout each winter (September–February) during 2009–2018, major roost sites on IoM and across the known winter range of migrant IoM shags (eastern and northern Scotland) were surveyed approximately every two weeks and resightings of colour-ringed individuals recorded (electronic supplementary material S1) [24,25]. Since breeding dispersal from IoM is very rare [24,26], these winter surveys allowed individuals to be directly classified as residents when resighted on IoM, and as migrants when resighted elsewhere. These resightings also effectively inform on true survival, with virtually no confounding permanent emigration.

(b) Model design

To estimate survival and reproductive selection on migration versus residence, we devised a discrete-time multi-event capture–recapture model that considers hidden transitions between individual states and imperfect observation of these states. States are defined as locations (i.e. residence on IoM, versus migratory areas elsewhere; electronic supplementary material S1) and breeding outcomes (i.e. breeding status and number of fledglings). The state transition process thus represents seasonal movement and survival, and subsequent breeding success. The observation process represents spatially and temporally varying resighting effort, and uncertainty in breeding success assessment. The model thereby allows robust probabilistic inference on the (partially observed) full-annual-cycle sequence of individual phenotypes (migrant or resident), and hence migration-dependent survival and breeding success (i.e. survival and reproductive selection; figure 1a).

• Download figure
• Open in new tab
• Download PowerPoint

To maximize the use of available year-round resighting data to make probabilistic inference on individuals' winter locations, and hence phenotypes and resulting selection, we divided each annual cycle (y; one breeding season to the next) into five capture-resighting occasions (o), comprising the breeding season and four subsequent winter occasions (figure 1a; electronic supplementary material S1) [24]. In each breeding season, new adults enter the dataset and all alive individuals are assumed to be in the residency area (on IoM). Through the four subsequent winter occasions (figure 1a), alive individuals can be in the residency area or a migratory area (corresponding to resident and migrant states). At each occasion, alive individuals can be seen where they are located or not seen, according to the occasion- and location-dependent resighting probability (p). To model spatial heterogeneity in the observation process, we defined multiple migrant states, including a ‘ghost area’ encompassing sites with no resighting effort (i.e. an unobservable state with p = 0, electronic supplementary material S1) [24]. Between occasions, survival probability (ϕ) is sex, time- (i.e. occasion by year) and migration-dependent (i.e. all migrants versus residents). Individuals can move between residency and migratory areas, according to probabilities of departing from residency (ε), moving to a particular migratory area conditional on departure (δ), returning from a migratory area to the residency area (ω) and switching between migratory areas conditional on not returning (σ) [24]. The model was parameterized with interacting sex-, location- and time-dependent in movement and resighting probabilities. However, σ was set constant across locations and time, because switching between migratory areas between winter occasions was rarely observed [24]. Since the data did not suggest any major sex bias in migrants' destinations, δ and σ were modelled as sex-independent. Since an individual cannot be a migrant during the breeding season, parameters were constrained such that individuals can only move from or remain in the residency area between occasion 1 (breeding season) and 2 (first winter occasion) and can only move to or remain in the residency area between occasion 5 (late-winter occasion) and 1 (figure 1a). Full details of non-breeding season model structure and parameterization are in [24].

In each breeding season (occasion o = 1), adults that survived the previous winter transition to one of six possible breeding states conditional on whether they were resident or migrant in the preceding late-winter occasion (o = 5; figure 1b). Specifically, individuals first become breeders (B), or conversely transition to the early-failed/non-breeder state ($\mathrm{B}/$) according to migration-dependent breeding probability ζ. Breeders then produce n fledglings (0 ≤ n ≤ 4) and thus transition to the corresponding states (B_n), following the migration-dependent set of nest outcome probabilities χ_n (with ${\sum }_{n=0}^4{\chi }_n=1$; figure 1b). The model was parameterized with interacting migration- and year-dependence in ζ and χ_n. Because there are as yet no capture–recapture methods allowing the breeding outcome to be modelled as a joint state for two paired individuals, the breeding outcome was modelled for one sex at a time (hereafter ‘focal sex’) and hence treated as independent from the other sex. Corresponding observation events for focal sex individuals are resighted as early-failed/non-breeder ($\mathrm{B}/$), resighted as a breeder with the success of n fledglings (B_n) or unknown success (B_?), or not resighted (Ø) (figure 1b). We assume that, since nests are exhaustively monitored (electronic supplementary material S1), all breeders are resighted (i.e. p_B = 1 for each breeder state). Breeder states are consequently either recorded with certainty or unknown (i.e. a breeder that produced n fledglings can only have observation event B_n or B_?), following state-dependent assignment probability α_Bn (figure 1b). However, early-failed/non-breeders can be resighted or not, with resighting probability ${\mathit{p}}_{\mathrm{B}/}$. All surviving non-focal sex individuals transition into a single live state unlinked to reproduction and are resighted with probability p [24]. However, winter observations of non-focal sex individuals inform movement parameters of both sexes and therefore improve precision and accuracy of estimates for the focal sex.

(c) Data

We used 42 322 year-round resightings to compile individual capture-resighting histories (i.e. sequences of observation events) of 2147 known-sex adult shags that bred on IoM at least once during 2009–2017 (19 011 resightings of 1108 females; 23 311 resightings of 1039 males). Each individual was first assigned to the residency area in the summer of its first observed breeding attempt during 2009–2017, then assigned as observed in an area or unobserved in each subsequent occasion, with a specific breeding observation in summer (electronic supplementary material S1).

Breeding season events comprised 2569 and 3004 direct observations of breeding success (and 31 and 25 unknown success) for females and males, respectively. The success of each individual's first observed breeding attempt during 2009–2017 is excluded from current analyses. This is because previous winter location (and hence migrant versus resident state) cannot be inferred for individuals that entered the dataset in summer 2009 (before winter resightings started; 340 females, 382 males) or were ringed as breeding adults during 2010–2017 (187 females, 127 males). Other individuals originally ringed as chicks entered the dataset at recruitment, typically aged three years (563 females, 508 males). However, individual pre-recruitment histories cannot be included without further assumptions, model developments and data regarding natal dispersal and recruitment processes. Since first breeding attempts were necessarily excluded, we did not aim to test general hypotheses regarding age-specific breeding success. Accordingly, we retained individuals of known and unknown ages (ringed as chicks and adults, respectively) in the data. However, to confirm that estimated differences in breeding success between residents and migrants, and hence apparent reproductive selection, were not simply due to correlated effects of age (e.g. if younger individuals were independently likely to migrate and to have low breeding success), we fitted further multi-event models that included basic age structure in transitions to breeding state (electronic supplementary material S6).

(d) Model analyses

The model was built and analysed using Stan, a Bayesian probabilistic programming language using Hamiltonian Monte Carlo, with package rstan [27] in R v. 3.6.3 (code in electronic supplementary material S2) [28]. Objective (uninformative) uniform priors were used for all parameters (electronic supplementary material S2). Posterior predictive checks showed no major discrepancies between the data and posterior predictions, implying good model fit (electronic supplementary material S3). Complete details on posterior samples, including for elementary model parameters, are in electronic supplementary material S5 and [29].

We derived posterior distributions for compound quantities of biological interest that are not elementary parameters, thus synthesizing key effects while propagating associated uncertainty. Annual survival probability (Φ_s) for possible phenotypic sequences (s) of seasonal residence versus migration through the annual cycle (figure 1a) was calculated as the product of survival probabilities of focal migratory phenotypes π (R or M) across the five successive occasions within each year (ϕ_π,o). For current purposes, we focused on two stereotypical and biologically relevant sequences: ‘full-winter migration’ (leaving the residency area by September and returning next breeding season, s = R-M-M-M-M, hence Φ_RMMMM,y = ϕ_R,1ϕ_M,2ϕ_M,3ϕ_M,4ϕ_M,5), and ‘full-winter residence’ (s = R-R-R-R-R, Φ_RRRRR,y = ϕ_R,1ϕ_R,2ϕ_R,3ϕ_R,4ϕ_R,5). Model estimates show that these sequences are the two most frequently realized: across 2009–2018, the posterior mean of the probability of full-winter residence ranged from approx. 20–50% (grand mean 30%) and full-winter migration from approx. 10% to 30% (grand mean 20% [24]). Each alternative path was unlikely: grand mean range 0–7% [24]. The two stereotypical sequences also capture key biological variation because major differences in Φ_s among sequences are driven primarily by late-winter survival; all sequences ending as migrant (or resident) in late-winter occasion have similar survival probabilities [24]. Based on model estimates across years, approximately 35–75% of individuals were residents in the late-winter occasion (grand mean 60%) while 25–65% of individuals were migrants (grand mean 40%). There was no strong sex bias in the proportion of migrants, or hence of sexual dimorphism in late-winter location (female-male difference ranged −0.18–0.07, grand mean −0.04; [24]).

To summarize the distribution of breeding success across all focal sex individuals (including defined breeders and early-failed/non-breeders) dependent on residence or migration in late-winter occasion, we derived the expected number of fledglings per individual (hereafter ‘expected breeding success', E(BS_π) where π denotes resident or migrant phenotype): effectively $E(B{S}_{\pi })={\sum }_{n=1}^4{\zeta }_{\pi }{\chi }_{\pi \mathit{,}\mathit{n}}\mathrm{n}$ (electronic supplementary material S5). Underlying full probability distributions of breeding outcomes are summarized in electronic supplementary material S5.

We computed annual fitness as the expected contribution to population growth (E(C_s)) of adult residents and migrants encompassing survival probability from year y to y + 1 followed by breeding success in year y + 1, for a given annual phenotypic sequence s ending with phenotype π, such that E(C_s) = Φ_s(1 + ½E(BS_π)). This measure represents the expected direct contribution of an adult alive following a given breeding season to the population immediately following the next breeding season (i.e. a post-breeding census), conditional on being resident or migrant. It comprises the expected contribution of an individual itself and half its expected number of offspring (given that all offspring have two parents), both conditional on annual survival. It is broadly analogous to the pre-breeding census formulation used elsewhere [14].

To quantify sex-specific selection on residence versus migration, we computed full posterior distributions of the differences (Δ) in Φ_s, E(BS_π) and E(C_s) between residents and migrants within each sex and year. We assessed evidence for each difference through the posterior probability that it was positive (Pr(Δ > 0)). Pr(Δ > 0) values close to 1 or 0 indicate substantial evidence for positive or negative differences, respectively (and hence for selection in one direction or the other), while values close to 0.5 indicate no clear evidence for selection in either direction. To explicitly test for variation in components of selection between sexes and years, we computed the difference in the resident-migrant difference (Δ_Δ, and corresponding Pr(Δ_Δ > 0)) between females and males within each year, and between every pair of years within each sex (electronic supplementary material, S4). Because the ordering of years within pairwise comparisons and resulting directionality is arbitrary when summarizing across several comparisons, we report values as the distance from the [0,1] boundaries rather than as the absolute values (Pr(Δ_Δ)’; electronic supplementary material, S4). Accordingly, Pr(Δ_Δ)’ values close to 0 indicate substantial evidence for a difference in selection between two focal years (in a given direction), while values close to 0.5 indicate no clear evidence in either direction. These measures of differences in selection combine differences in magnitude and direction. To explicitly test for fluctuating selection (i.e. differences in direction), we calculated the posterior probability of a sign change in selection between every pair of years within each sex (electronic supplementary material S4) [4]. Pr(Δ_Δ±) values close to 1 indicate strong evidence that selection acted in opposite directions between two focal years, whereas values close to 0 indicate no evidence that selection differed in direction.

(a) Annual survival selection

Estimated annual adult survival probabilities (Φ_s) were notably lower in 2012–13, 2013–14 and 2017–18 than in the other six study years (figure 2). There were clear episodes of strong survival selection against full-winter residence in two of these years (2012–13 and 2017–18; figure 2). Here, Φ_s was approximately 0.2 higher for full-winter migrants than residents in both sexes (figure 2; electronic supplementary material S5). These episodes of mortality and survival selection were previously noted to coincide with late-winter ECEs comprising severe storms [24].

• Download figure
• Open in new tab
• Download PowerPoint

By contrast, there was no evidence of substantial differences in Φ_s between migrants and residents, or hence of survival selection, in the other 7 years (figure 2; grand mean of Δ was 0.00 in females and 0.02 in males; electronic supplementary material, S5). Survival selection in 2012–13 and 2017–18 exceeded that in the other years in both sexes (figure 2), as confirmed by pairwise comparisons (Pr(Δ_Δ)’ was less than 0.03 and less than 0.11 for 2012–13 and less than 0.01 and less than 0.12 for 2017–18, in females and males, respectively; electronic supplementary material S5). Further, survival selection in 2012–13 did not differ from that in 2017–18 (Pr(Δ_Δ)’ was 0.45 in females and 0.20 in males; electronic supplementary material S5). There were no marked differences in selection among the other 7 years, except for some evidence that selection against resident males differed between 2015–16 and other years apart from 2009–10 and 2010–11 (figure 2, Pr(Δ_Δ)’ ranged 0.05–0.09; electronic supplementary material S5).

There was no evidence of strong or consistent sex-specific survival selection in any year (figure 2; Pr(Δ_Δ > 0) ranged 0.11–0.81 across years; electronic supplementary material S5). Hence, overall, there was strong among-year variation in survival selection, driven by two extreme years with low survival and strong selection against residents in both sexes. However, there was little definitive evidence of fluctuating survival selection, defined as sign changes, in either sex (Pr(Δ_Δ±) ranged 0.04–0.80 in females and 0.04–0.89 in males across all pairs of years; electronic supplementary material S5).

(b) Reproductive selection

Expected breeding success (E(BS_π)) varied substantially among years, and there was evidence of strong reproductive selection against late-winter migrants in multiple years (figure 3; recall that for E(BS_π), phenotype π reflects location in late-winter occasion o = 5; figure 1). Specifically, there was strong evidence of selection against migrants in both sexes in 2011, 2012 and 2014, and in females in 2016. For these years, posterior means of Δ ranged 0.26–0.64 and 0.28–0.45 fledglings in females and males, respectively (figure 3; electronic supplementary material S5).

• Download figure
• Open in new tab
• Download PowerPoint

Evidence of reproductive selection was weaker in other years, but estimated effects were typically in the same direction (i.e. against late-winter migrants, figure 3). However, there was weak evidence for reversed selection (i.e. against late-winter residents) in females in 2013 and 2018 (figure 3; electronic supplementary material S5), which are the 2 years that also showed strong survival selection against residence (figure 2). Consequently, pairwise comparisons confirmed that reproductive selection varied among years in females. In particular, there was strong evidence that selection was different in 2013 and 2018 compared to other years (figure 3; Pr(Δ_Δ)’ ranged 0.00–0.14 and 0.01–0.26 respectively; electronic supplementary material S5), but not to each other (Pr(Δ_Δ)’ = 0.37). Among the other 7 years, there was some evidence that selection was different in 2011 and 2012 compared to later years (Pr(Δ_Δ)’ ranged 0.05–0.26), but not otherwise (Pr(Δ_Δ)’ ranged 0.22–0.45). However, there was only weak evidence of fluctuating reproductive selection in females; Pr(Δ_Δ±) ranged 0.72–0.82 when comparing 2013 with other years excluding 2017 and 2018, but ranged 0.01–0.64 otherwise (figure 3; electronic supplementary material S5).

Evidence of among-year variation in reproductive selection was weaker in males. Here, selection against migrants differed between 2012 and later years (figure 3; apart from 2014, Pr(Δ_Δ)’ ranged 0.07–0.10; electronic supplementary material S5), but did not differ between any other years (Pr(Δ_Δ)’ ranged 0.14–0.48; electronic supplementary material S5). There was consequently no strong evidence of fluctuating reproductive selection in males (Pr(Δ_Δ±) ranged 0.04–0.51). Finally, even though there was stronger evidence of among-year variation in reproductive selection in females than males, there was no strong evidence of sex-specific reproductive selection on residence versus migration in any year (figure 3; Pr(Δ_Δ > 0) ranged 0.13–0.91; electronic supplementary material S5).

(c) Selection through annual fitness

Analyses of the composite measure of annual fitness (E(C_s)) showed clear evidence of episodes of strong fluctuating selection on full-year residence versus migration (figure 4). Specifically, there was strong evidence of sign changes in selection in 2012–13 and 2017–18 compared to several other years: 2009–10, 2010–11 and 2011–12 in females, and 2011–12 and 2013–14 in males (Pr(Δ_Δ±) for these pairwise years ranged 0.92–0.99 in females and 0.94–1.00 in males; electronic supplementary material S5). However, 2012–13 and 2017–18 did not differ from each other (Pr(Δ_Δ±) was 0.01 in both sexes; electronic supplementary material S5). There was no strong evidence of sex-specific overall selection on residence versus migration within any year (Pr(Δ_Δ > 0) ranged 0.22–0.89; figure 3; electronic supplementary material S5). Variation and fluctuations in selection were consequently broadly synchronized across both sexes.

• Download figure
• Open in new tab
• Download PowerPoint

The fluctuations were caused by the underlying episodes of strongly opposing survival and reproductive selection that occurred in different years. Specifically, E(C_s) was substantially higher in migrants than residents in both sexes in the 2 years with strong survival selection (2012–13 and 2017–18; figure 4; electronic supplementary material S5). Conversely, E(C_s) was substantially higher in residents than migrants in years with strong reproductive selection, most notably in 2009–10 to 2011–12 and 2015–16 in females and in 2011–12, 2013–14 and 2014–15 in males (Δ ranged 0.15–0.34 and 0.15–0.32 in females and males, respectively; figure 4; electronic supplementary material S5).

(a) Estimation and implications

In general, estimates of selection can be biased by missing or error-prone phenotypic and/or fitness data [45,46], and evidence of varying and fluctuating selection should be evaluated given sampling variance [4]. These challenges are ubiquitous but commonly ignored and come to the fore when focal phenotypes are not always directly observed with certainty [47]. Uncertainty is inevitable for resident versus migrant phenotypes inferred from resighting data. Our multi-event analyses accounted for key uncertainties by modelling missing phenotypic and survival data due to resighting failure (including due to migration outside surveyed locations), and missing reproduction data due to non-breeding and unobserved nest outcomes. The Bayesian implementation allowed straightforward computation of posterior distributions of between-sex and among-year differences and sign changes that evidence varying, fluctuating and cross-sex synchrony in selection. This approach differs somewhat from standard regression approaches to estimating selection, and recent extensions, designed for readily observable continuous traits. For example, we did not estimate variance in selection through temporal random effects [4,8]. Such estimation is not straightforward in our case, since annual survival and expected breeding success and annual fitness are all derived parameters. Further, variances may be poorly estimated across relatively few years, and assumptions regarding Gaussian distributions may be violated (electronic supplementary material S5). Our analyses illustrate how variation and fluctuations in selection can be robustly quantified in such (common) circumstances.

The variation and fluctuations in selection on the defined migrant versus resident phenotypes estimated in adult shags could appreciably affect phenotypic dynamics and underlying additive genetic variation. Since both phenotypes are frequently expressed, there is considerable phenotypic variance on which selection can act. Indeed, such variance could be partly maintained by fluctuating selection, and the apparent lack of sexual dimorphism in migration is consistent with the cross-sex synchrony in selection. Nevertheless, further steps that consider how the effects of fluctuating selection are propagated across years and life-history stages are required to fully consider the eco-evolutionary consequences. In a moderately long-lived species such as shags, selection through adult survival is likely to impact strategy-specific population growth rates at multi-year scales more than selection through reproduction. However, such effects will also depend on the degree of phenotypic plasticity following selection episodes, and on forms of selection acting before recruitment, which remain to be quantified in our system. Moreover, while our focus on individual full-winter migration versus residence captures considerable annual phenotypic and fitness variation, subtle forms of selection might act on the full diversity of possible phenotypic sequences of seasonal residence versus migration [24,48] and on the joint phenotypes of breeding pairs [31]. Future analyses with further data and methodological developments, including analyses of paired capture-resighting histories, will allow us to quantify these components and thereby reveal the full form and consequences of varying and fluctuating selection on liability for migration.