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

Quantifying temporal variation in sex-specific selection on key ecologically relevant traits, and quantifying how such variation arises through synergistic or opposing components of survival and reproductive selection, is central to understanding eco-evolutionary dynamics, but rarely achieved. Seasonal migration versus residence is one key trait that directly shapes spatio-seasonal population dynamics in spatially and temporally varying environments, but temporal dynamics of sex-specific selection have not been fully quantified. We fitted multi-event capture–recapture models to year-round ring resightings and breeding success data from partially migratory European shags (Phalacrocorax aristotelis) to quantify temporal variation in annual sex-specific selection on seasonal migration versus residence arising through adult survival, reproduction and the combination of both (i.e. annual fitness). We demonstrate episodes of strong and strongly fluctuating selection through annual fitness that were broadly synchronized across females and males. These overall fluctuations arose because strong reproductive selection against migration in several years contrasted with strong survival selection against residence in years with extreme climatic events. These results indicate how substantial phenotypic and genetic variation in migration versus residence could be maintained, and highlight that biologically important fluctuations in selection may not be detected unless both survival selection and reproductive selection are appropriately quantified and combined.


Introduction
Quantifying temporal variation in the strength and direction of sex-specific selection on ecologically relevant phenotypic traits is central to understanding eco-evolutionary dynamics [1][2][3][4]. This is because the forms and magnitudes of variation in selection will shape the maintenance of genetic and phenotypic variation, and shape the rate and direction of adaptive evolutionary change [2,3,5]. Temporal variation in selection will thus fundamentally affect population responses to varying and changing environmental conditions. In particular, both fluctuating selection and sexually antagonistic selection, respectively defined as episodes of selection acting in opposite directions within short ecologically relevant periods or between the sexes, can help maintain genetic variation and alter time frames for adaptation [3,5]. Yet, temporal dynamics of sex-specific selection on key traits in wild populations have still rarely been quantified [2].
Empirical evidence of temporally fluctuating selection is particularly scant, once sampling variance is accounted for [4,6]. Moreover, we commonly lack insights into how fluctuations are caused by environmental variation, even though such impacts are central to eco-evolutionary processes and outcomes [2,7,8].
In general, selection on any trait can operate through differential survival and/or differential reproduction in relation to phenotype, yielding survival selection and/or reproductive selection [9,10]. These selection components could act in the same or opposite direction, generating either strong or weak net selection within years [11]. Further, the relative strength and direction of selection through each fitness component could vary among years, potentially generating net fluctuating selection. Moreover, depending on sex-specific responses to underlying environmental variation, fluctuations in net selection could be synchronized or opposing across females and males. Comprehensive studies aiming to quantify temporal variation in the overall magnitude and direction of selection should therefore quantify sex-specific temporal variation through both fitness components, and through their combined effects.
This ambition necessitates explicit estimation of temporal sequences of sex-specific selection within as well as among years. In many systems, reproduction (and resulting reproductive selection) occurs within discrete seasons, while survival selection could occur at any time, and might not coincide with reproductive selection. Such sequential selection episodes could have complex compound effects. For example, strong survival selection acting during a non-breeding season will leave fewer individuals of particular phenotypes available to breed subsequently, reducing the degree to which opposing reproductive selection could reverse the direction of net annual selection [12]. Further, carry-over effects of non-breeding-season phenotypes on subsequent reproduction could cause additional components of time-lagged indirect selection [13]. However, most studies of temporal variation in selection focus on either survival or reproduction [2,6], and/or do not estimate combined effects of different selection components across seasons. Fitness measures that combine survival and reproduction to quantify individual contributions to annual population growth have been developed [10,14,15], but are still rarely applied to estimate selection [16,17].
One key phenotypic trait that could directly link ecological and evolutionary dynamics is seasonal migration (hereafter 'migration'), defined as reversible individual movements between locations across seasons. Migration allows individuals to exploit seasonally varying resources and avoid unfavourable conditions, and directly determines individuals' seasonal locations and resulting spatio-seasonal population distributions [18,19]. Moreover, phenotypic expression of migration versus residence commonly varies among individuals within populations, creating opportunity for selection. Specifically, forms of 'partial migration', where some individuals remain resident at their breeding location during the non-breeding season while other individuals are seasonal migrants, occur in diverse amphibian, reptile, fish, bird and mammal populations [18,[20][21][22]. Sympatric-breeding individuals with migrant and resident phenotypes are then spatially segregated in the non-breeding season. Episodes of strong seasonal selection could then arise due to spatial variation in non-breeding season environmental conditions that causes differences in survival and/or subsequent reproduction, which may be modulated by sex-specific environmental tolerances and/or constraints on the reproductive success [18,23]. Quantifying among-year variation in sex-specific selection on migration versus residence is therefore central to understanding how spatio-temporal environmental variation could drive micro-evolution of migration and hence drive micro-evolution of spatio-seasonal population dynamics and distributions. Yet, to date, such variation in selection has not been fully quantified.
Progress requires quantifying non-breeding season phenotype (resident or migrant), survival and subsequent reproduction of numerous females and males across multiple years within a sympatric-breeding partially migratory population. This can be achieved through large-scale year-round resightings of marked individuals designed to determine individuals' non-breeding season locations, coupled with subsequent reproductive monitoring. However, since not all individuals' locations and reproduction will typically be observed at all times, inference of selection requires advanced full-annual-cycle capture-recapture models that account for the resighting process and resulting partial observation of individuals' phenotypes and uncertainty in survival and breeding outcomes. Recent analyses in European shags (Phalacrocorax aristotelis, hereafter 'shags') demonstrated strong survival selection against residence in both sexes within two of nine non-breeding seasons containing extreme late-winter storms (i.e. extreme climatic events, 'ECEs'), with weak selection or neutrality otherwise [24]. However, the degree to which such temporal variation in survival selection could be overriden by sex-specific reproductive selection manifested through carry-over effects acting in subsequent breeding seasons, potentially generating overall fluctuating selection, has not been quantified.
Accordingly, we fitted multi-event capture-recapture models to year-round resightings and breeding success data from adult shags to jointly estimate annual sex-specific reproductive selection alongside survival selection. We then combined these estimates to quantify among-year variation in overall selection on migration versus residence through annual adult contribution to population growth, explicitly tested for fluctuating selection, and examined whether variation and fluctuations were broadly synchronized across females and males. We thereby quantify how components of seasonal selection, including selective episodes associated with ECEs, can drive strong fluctuating and/or sex-specific selection on a key phenotypic trait that shapes spatio-seasonal population dynamics.

Methods (a) Study system and data collection
The focal shag population breeds on the Isle of May (IoM) National Nature Reserve, Scotland (56°11 0 N, 2°33 0 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 royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 288: 20210404 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 −1 , 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 −1 ) 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).
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 royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 288: 20210404 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 (B =) according to migrationdependent 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 P 4 n¼0 x 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 (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 p 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 royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 288: 20210404 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 earlyfailed/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(BS p ) ¼ P 4 n¼1 z p x p,n 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.  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].

Results
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

(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 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 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 amongyear 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).

Discussion
Eco-evolutionary dynamics will partly depend on the magnitude and between-sex synchrony of temporal variation in the strength and direction of selection on key ecologically relevant traits that shape spatio-temporal population dynamics. We demonstrate episodes of strong, and strongly fluctuating, selection on one such key trait, seasonal migration versus residence, through annual fitness in adult European shags (figure 4). Fluctuations were broadly synchronized across males and females and arose because strong reproductive selection against migration in several years contrasted with episodes of strong survival selection against residence in 2 other years, coupled with weakened reproductive selection in the subsequent summers (figures 2 and 3). While the underlying components of both survival and reproductive selection varied substantially among years, including episodes of approximate neutrality, neither fluctuated strongly (strictly defined as a change in direction). Other empirical studies testing for fluctuating selection commonly consider only one component or the other, or do not quantitatively combine them into a single annual fitness measure (e.g. [2,4,6]). Our analyses demonstrate the value of doing so, since the presence of overall strong fluctuating selection would not otherwise have been evident. Predictions regarding micro-evolutionary outcomes based on single-component estimates of selection might consequently be misleading.
In our system, the two notable episodes of strong survival selection against residence were associated with ECEs that occurred during late winter in 2012-13 and 2017-18 [24]. Here, prolonged periods of strong onshore wind, rain or cold, which reduce foraging efficiency and incur high thermoregulatory costs in shags [30], primarily impacted the residency area, probably causing increased mortality in residents [24]. Location, and hence migrant versus resident phenotype, would then directly affect survival probability, constituting direct selection. However, because these phenotypes are only expressed during the non-breeding season, any reproductive selection against migration must constitute time-lagged indirect selection. This could reflect 'carry-over' effects, for example, if an individual's migrant versus resident phenotype affects its condition and/or capability to acquire or retain a breeding site, which then affects its subsequent breeding success [31][32][33]. Reproductive selection against migrants was also weakened, or perhaps even reversed in females, in the summers following the ECE winters (2013 and 2018), particularly in older individuals (electronic supplementary material S6). Such patterns could arise if residents that survived through selective events were in poorer condition (e.g. due to poorer foraging conditions [34]). This effect could be more pronounced in females, because females' foraging efficiency is more negatively impacted by strong winds than that of males [30]. Temporal variation in selection may therefore predominantly reflect underlying environmental conditions. However, formally demonstrating such links, and the intermediate physiological and/or behavioural mechanisms, is generally challenging and rarely achieved [2,7,8,35]. Our results set up valuable opportunities to attempt such analyses once longer timeseries, including multi-dimensional environmental data at appropriate spatio-temporal scales, can be assembled [24]. Nevertheless, our current results indicate that environmental variation including ECEs can generate major reversals of the direction of net selection in annual fitness, encompassing both direct and time-lagged indirect components. Depending on the frequency of ECEs, and on underlying environmental versus additive genetic (co)variances and forms of phenotypic plasticity, such fluctuations may contribute to maintaining  royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 288: 20210404 additive genetic and phenotypic variation [3], yet alter the optimal genotype and phenotype and constrain rates of adaptive micro-evolution [36][37][38]. Classic examples of strongly fluctuating selection on other traits in other systems have also been linked to extreme environmental variation. For example, drought-induced changes in food supply changed the direction of selection on beak morphology in Darwin's finches [39]. Seasonal migration versus residence is a good candidate trait for such effects because sympatric-breeding migrant versus resident phenotypes are, by definition, completely spatially segregated in the non-breeding season and hence subject to different environmental conditions. However, although studies across diverse taxa have now quantified differences in components of survival and/or reproduction between such residents and migrants [18,22], none have explicitly quantified selection through composite measures of annual fitness. Further, some previous studies pooled data across years to generate sufficient sample sizes, precluding estimation of temporal variation selection, such as elk (Cervus elaphus) [40], moose (Alces alces) [41] and European blackbird (Turdus merula) [42]. Others showed little or no temporal variation and no strong evidence of fluctuating selection, albeit across a few years and individuals (e.g. red-spotted newt (Notophthalmus viridescens) [23], skylark (Alauda orvensis) [43] and pronghorn (Antilocapra americana) [44]). Our findings (which encompass the whole adult population) broadly concur with earlier analyses in the same study system, which used a subset of resighted individuals with known winter locations during 2009-2012 to show that resident shags had consistently higher breeding success than migrants (figure 3) [31]. Our longer timeseries now shows temporal variation in selection in subsequent years. In general, very short time frames are a common limitation across studies of variation and fluctuations in selection (median 3 years [2]) and may miss biologically important fluctuations caused by infrequent environmental perturbations.

(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 multievent 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.
Ethics. Bird ringing was licensed by the British Trust for Ornithology ( permit numbers A400 and A4607). Access to field sites and annual research licences from NatureScot (MON/RP181 for 2018).
Funding. This study was funded by Joint Nature Conservation Committee and Natural Environment Research Council (grant nos NE/M005186/1, NE/R000859/1 and NE/R016429/1).