Accepted Ar tic le Beak morphology predicts apparent survival of crossbills: due to selective survival or selective dispersal? David Gómez Blanco1, Simone Santoro2, Antoni Borrás3, Josep Cabrera3, Juan Carlos Senar3 and Pim Edelaar2 1 Dept of Biology, Lund Univ., Lund, Sweden 2 Dept of Molecular Biology and Biochemical Engineering, Univ. Pablo de Olavide, Seville, Spain 3 Behavioural and Evolutionary Ecology Research Unit, Museu de Ciències Naturals de Barcelona, Barcelona, Spain Corresponding author: David Gómez Blanco, Dept of Biology, Lund Univ., Lund, Sweden. E-mail: david.gomez_blanco@biol.lu.se Decision date: 08-Oct-2019 This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process, which may lead to differences between this version and the Version of Record. Please cite this article as doi: [10.1111/jav.02107]. ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Abstract Dozens of morphologically differentiated populations, subspecies and species of crossbills (genus Loxia) exist. It has been suggested that this divergence is due to variation in the conifer cones that each population specialises upon, requiring a specific beak size to efficiently separate the cone scales. If so, apparent survival should depend on beak size. To test this hypothesis, we undertook multievent capture-recapture modelling for 6,844 individuals monitored during 27 years in a Pyrenean common crossbill L. curvirostra population in a forest of mountain pine Pinus uncinata. Apparent survival was indeed related to beak width, resulting in stabilizing selection around an optimum that was close to the observed mean beak width, indicating that local crossbill beak morphology is adapted to the conifer they feed upon. Both natural selection (selective mortality) and selective emigration of maladapted individuals may explain our findings. As is often the case in capture-recapture analyses but rarely recognised, we could not formally decompose apparent survival into selective mortality versus selective permanent emigration. Nonetheless, there are several indications that selective permanent emigration should not be fully excluded. First, natural selection by itself would have to be unusually strong compared to other empirical estimates to create the observed pattern of apparent survival. Second, the observed mean beak width was a bit lower than the estimated optimum beak width. This can be explained by immigration of crossbills with smaller beaks originating from southern populations, which may subsequently have left the study area permanently in response to low food intake. This is in line with a detected transient effect in the data, yet apparently little influx from crossbills from northern Europe. When permanent emigration is phenotypically selective this will have ecological and evolutionary consequences, so this possibility deserves more attention in general. Keywords: Loxia curvirostra, Pyrenees, apparent survival, local adaptation, matching habitat choice, selection of the environment ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Introduction Organisms are constantly facing ecological challenges due to variation in their biotic and abiotic environments, and the interaction between phenotype and environment modifies their survival and/or reproductive success (Darwin 1859). It is therefore adaptive to have a better match between phenotype and environment, in order to improve ecological performance. Several processes may achieve this improved match, operating both at the individual and population level (Edelaar and Bolnick 2019). One process is the classical local adaptation of populations (e.g. adaptive tracking and adaptive population genetic divergence), based on differences in fitness. If individual variation in phenotypic match with the environment is caused by heritable traits, natural selection causes adaptation to the environment across generations at the population level (Manly 1985, Endler 1986). However, in variable environments the observed phenotype distribution often lags behind shifts in the optimum. Therefore, a second well-known, widespread and important process is matching the phenotype to the environment by adaptive phenotypic plasticity (Black 1993, Przybylo et al. 2000, Ghalambor et al. 2007). A third potential process to deal with environmental variation is selection of the environment, by which an individual selects an environment that matches well with its phenotype, e.g. matching habitat choice (Edelaar et al. 2008). And finally, a fourth process that can improve the phenotype-environment match is adjustment of the environment (Edelaar and Bolnick 2019; not considered further in this paper because it is unlikely to operate a priori in our system). These processes can operate simultaneously to improve fitness, but their distinction is important for our understanding of the ecological processes underlying the patterns we observe in nature. By default, the route to increase the fit of the phenotype to the environment is ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le population genetic adaptation and differentiation through natural selection on nonplastic phenotypes. However, plasticity and selection of the environment are favoured as environmental heterogeneity increases and their costs are low relative to individual benefits (Scheiner 2016, Edelaar et al. 2017). The common crossbill Loxia curvirostra is a sparrow-sized songbird specialized in feeding on seeds enclosed in conifer cones. They are characterized by their thick beaks with curved and crossed bill tips that are used to extract seeds from between the scales of closed cones (Benkman and Lindholm 1991). As is typically the case in bird taxonomy, crossbill taxa have historically been defined by distinct geographic (allopatric) distributions, biometrics and plumage colouring (Cramp and Perrins 1994). In more recent times, morphologically, vocally and genetically distinct sympatric populations (if not incipient species) of crossbills are recognised (Groth 1993, Benkman 2016, Parchman et al. 2016). This seems to be linked to the availability of and adaptation to different conifers. Benkman (1993) showed that there is not a single beak phenotype that allows feeding optimally on all cone types. Given that the beak is a nonplastic trait (Summers et al. 2007) that develops before immatures start feeding on cones, adaptive plasticity is not an option to increase a crossbill’s food intake on a given cone type. Benkman (2003) subsequently showed for a resident crossbill population that return rate (a combination of survival and recapture) was correlated with bill morphology, with highest values for those birds that had the highest intake rate on the local cones. Therefore, by assuming return rates reflect apparent survival, divergent selection should act on bill traits among crossbill taxa feeding on different conifers. The use of alternative resources then would be the ultimate cause of adaptive radiation in this species complex. ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le In Spain crossbills are irregularly distributed in areas with coniferous forests (Borrás and Senar 2003). Currently, two subspecies are recognized: Loxia curvirostra curvirostra (the same subspecies as occurring across most of Eurasia), and L. c. balearica, from Mallorca (Balearic Islands) and according to some authors (see Alonso et al. 2006) also south-eastern continental Spain. L. c. curvirostra from northern Spain feeds almost exclusively on the seeds of mountain pine Pinus uncinata, scots pine P. sylvestris and black pine P. nigra (Clouet 2000, Borrás and Senar 2003). Loxia c. balearica feeds on Aleppo pine P. halepensis, more abundant in SE Spain and the only pine on the Balearic Islands (Alonso et al. 2006). The divergence in the bill morphology of crossbill populations in Spain and a reduction in gene flow between them appears to be due to the differential use of resources, i.e. isolation by ecology (Edelaar et al. 2012). Aleppo pine produces larger cones than mountain and scots pine (Castroviejo et al. 1986), but has relatively long and weak scales. The need to separate the shorter and stronger scales of mountain, scots and black pine cones therefore seems to cause the beak of N Spanish crossbills to evolve to be relatively shorter but wider than those of the birds in S Spain and Mallorca (Alonso et al. 2006, Borrás et al. 2008, Mezquida and Benkman 2010, Edelaar et al. 2012). However, such an effect of beak morphology on adaptation has not been shown, only assumed to be general based on the results of Benkman (2003). The first objective of our study is therefore to investigate if beak morphology affects apparent survival, and if so, to determine which size of beak is locally adaptive. For this we used an extension of the Cormack-Jolly-Seber capture-recapture model (Lebreton et al. 1992) to analyse data collected over 27 years from a Pyrenean crossbill population. The use of a statistical ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le capture-recapture model allows us to decompose observed return rates into the underlying apparent survival and recapture rates, therefore providing more accurate and unconfounded estimates of apparent survival. It was not beforehand obvious that the crossbill population that we studied would be locally adapted, since mountain pine has a very restricted distribution, only occurring in the subalpine zone of the Pyrenees and western Alps. Small ranges are typically associated with smaller population sizes, making populations more vulnerable to environmental and demographic stochasticity as well as maladaptive gene flow. However, mountain pine shows relatively little temporal and small-scale spatial variation in cone production (Senar et al. 1993, Clouet 2000, Borrás and Senar 2003), so this likely facilitates local adaptation and morphological specialisation since resource stability within the range of a population enhances the evolutionary life span of specialist populations (Björklund et al. 2013, Parchman et al. 2018). Irrespective of any relationship between beak morphology, ecological performance, and apparent survival, the specialisation of and divergence among sympatric crossbill taxa is surprising, because population contact and geographical overlap are normally assumed to prevent or erase divergence, due to homogenizing gene flow (Hendry et al. 2001). Such potential for homogenizing gene flow is particularly large in crossbills, due to large spatio-temporal variation in food availability. Many of the other conifers that crossbills feed on vary greatly from year to year in cone production, with booms and busts occurring across large areas (Clouet 2000, Borrás and Senar 2003). In years of low local cone production, crossbills therefore move to other areas in search of better cone crops of the same type of pine (or a suitable alternative). Flocks of crossbills ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le moving to areas with other pine species after local crop failures have often been observed (Senar and Borrás 1985, Clouet 2000, Borrás and Senar 2003, Borrás et al. 2011), sometimes in the form of massive invasions involving millions of birds. Indeed, the crossbills of northern Eurasia may well be the most dispersive passerine on Earth, with average natal and breeding dispersal distances of over 2,100 km (Newton 2006, Marquiss et al. 2008). These movements of individuals should increase the mixing of populations, effective gene flow, and the homogenization of genetic groups. However, alternatively, these movements could bring or keep together phenotypically similar individuals if, after evaluating different areas, they settle in the same area that best fits their ecological needs given their phenotype, i.e. via habitat choice. Matching habitat choice and its associated selective movement could actually promote, instead of preventing, differentiation of populations (Marquiss and Rae 2002, Siepielski and Benkman 2005, Edelaar et al. 2008, Holt and Barfield 2008, Armsworth and Roughgarden 2008, Edelaar and Bolnick 2012, Bolnick and Otto 2013, Berner and Thibert-Plante 2015), because matching habitat choice depends on a direct comparison of local performance (i.e. the phenotype-environment match) across different options. Benkman (2017) recently provided empirical support for matching habitat choice driving the evolution of bill size in the same location for which a link between bill size and return rate was reported (Benkman 2003). Hence, it remains unclear to what extent adaptation and ecological divergence in crossbills is only due to natural selection (e.g. selective survival) in local populations, and/or due to selective dispersal among populations (selection of the environment). That is because even though these two processes are opposites (in the first the phenotype is matched to the environment, in the second the environment is matched to the ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le phenotype), they can result in the same pattern of local phenotype-environment match. This uncertainty in driving processes is also the case for any study on ecological population divergence in which selective dispersal and settlement cannot be excluded, or in fact for any study on local adaptation. Nonetheless, this latter possibility is often ignored. Therefore, as a second objective, we evaluate the possible roles of natural selection (mortality) versus selection of the environment (permanent dispersal) in explaining patterns of apparent survival in our Pyrenean crossbill population. Material and methods Study population and data collection Capture of crossbills took place in the Port del Comte mountain range in the Catalan Pre-Pyrenees, approximately 100 km NW of Barcelona. The only pine occurring here is Mountain pine Pinus uncinata. Mist nets were placed where birds came down to drink or feed on minerals. Captured individuals were marked with individually numbered aluminium rings and released at the capture site. Sex and age were determined according to Svensson (1992). We considered three age groups. Juveniles are those birds that have not yet moulted into adult plumage, so their sex is difficult to determine visually (although see Edelaar and Terpstra 2004, Del Val et al. 2014). Yearlings are birds in adult plumage but with some remainder of juvenile plumage (incl. some greater coverts and flight feathers, which are typically retained until the next full moult). Adults are individuals that have gone through a complete moult. Biometric traits were measured by a single observer, incl. beak width (base of the lower mandible, ± 0.1 mm) and beak height (height of closed mandibles at the nostrils, ± 0.1 mm). All measurements showed high repeatability: beak width, r = 0.93, p < 0.001, beak height r ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le = 0.98, p < 0.001 (N = 10 in both cases, Borrás et al. 2008). Repeatability was measured by an analysis of variance ANOVA on two measures per individual (Lessells and Boag 1987, Bailey and Byrnes 1990). The capture-recapture dataset The birds were ringed and recaptured between 1988 and 2014. Given that recaptures were relatively scarce, we combined all the observations from May to October of the same year into a single capture occasion. This yields higher recapture rates, which give more precise and less biased estimates of state transition parameters like survival rate (Hargrove and Borland 1994, O’Brien et al. 2005). Nonetheless, the resulting dataset consisted of a total of 6,844 marked crossbills and only 272 (4%) recaptures. Furthermore, due to the lack of biometry for some individuals, two reduced datasets were constructed for birds with measures of beak width or beak height (respectively, 949 with 82 (9%) recaptures and 1,929 with 123 (6%) recaptures of crossbills marked and measured). For the few occasions when a crossbill was measured multiple times, its first measurement was used in the analyses. Multievent capture-recapture modelling Multievent models were constructed and adjusted to the data using the software ESURGE 2.1.2 (Choquet and Nogue 2011). The multievent framework allows one to use all the data by distinguishing between the events, i.e. the individual states as they have been observed in the field (often incomplete or uncertain), and the underlying biological states of the individuals, which must be inferred (Pradel 2005). We defined three underlying biological states: dead=†, female alive=F, male alive=M, and four possible events: 0=non-detected, 1=female-detected, 2=male-detected, 3=sex unknown-detected. ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Age groups have been coded within the “headed format” (see Choquet and Nogue 2011) data file as: (1) juvenile, (2) yearling and (3) adult (details on the probabilistic framework are given in Supplementary material Appendix 1). The parameters (initial state, state transition and event) of multievent models are calculated simultaneously from the full encounter histories by maximum likelihood (Choquet et al. 2009a). Goodness of fit Since no goodness of fit (GOF) tests are available for multievent models, we tested the fit of the Cormack-Jolly-Seber model (CJS) using U-CARE ver. 2.3.4 (Choquet et al. 2009b). The CJS is a model that allows for unspecified temporal variation in the probabilities of apparent survival and of recapture (Lebreton et al. 1992). The CJS model only distinguishes two biological states (alive, dead) and two possible events (capture, no-capture). The GOF analysis was performed on our three groups that were defined according to the age at first capture. U-CARE allows testing whether the model fits the data and investigating the potential causes of lack of fit (if any). One reason for lack of fit is known as the “transient effect”, when the apparent survival probability in the first interval after marking is substantially lower than in the subsequent intervals (test component 3.SR). A transient effect can be caused by the presence of individuals which use the study area sporadically and/or by individuals with systematically lower survival than others (Duriez et al. 2009). Another cause of lack of fit is known as trapdependence, when individuals captured (or recaptured) in the previous occasion have a different probability of being recaptured than the other individuals (test component 2.CT). When the global GOF test (which merges all the components together) is significant, the extra-binomial sources of variation have to be accounted for in the ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le structure of the capture-recapture models and/or an over-dispersion coefficient (c-hat) must be used in the analyses to make the model selection more conservative. For our data, we did not find evidence of trap-dependence (Z= -0.63; P = 0.53). However, we did find a transient effect (χ2 = 86.44; df = 144; Z= 3.27, P <0.001). When the analyses were performed separately for the three age groups, we still found statistically significant transient effects for juveniles (χ2 = 28.85; df = 22; P <0.01) and yearlings (χ2 = 8.26; df = 14; P = 0.04), but less so for adults (χ2 = 8.53; df = 18; P = 0.13). In order to control for the detected transient effects, we included an age effect (three age classes) on survival in the global model (i.e. the most complex model, see below). Model selection We first modelled the entire dataset, independent of whether morphological measures were taken. In our global model (the initial full model, before simplification) we included temporal variation in all the parameters (initial state, survival, capture probability and sex-assignment), including in interaction with sex and age for some parameters when the data allowed it (i.e. models not over-parameterised). We followed a step-down model selection approach (Lebreton et al. 1992) based on the quasi Akaike Information Criterion corrected for small sample size (QAICc) (Burnham and Anderson 2004) that consisted of sequentially determining the best structure of the parameters in this order: Sex-Assignment, Capture probability, Initial State, and Survival (see Santoro et al. 2016 for an analogous approach). To this aim, we ran for the first parameter (sexassignment) all the nested models for the entire dataset, retained the structure for this parameter given by the model with the lowest QAICc, and repeated the same procedure ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le with the next parameter (Capture probability), etcetera. As more data is involved, using the entire data set gives us the best demographic estimates and the best assessment of which factors affect each parameter. Next, we modelled the reduced datasets with the individual covariates “beak width” and “beak height”, in order to test if beak morphology influenced apparent survival. For this, the structure of the simplified final model based on all data was used as the global model, and the step-down model selection procedure was repeated. Then, after this dataset-specific model selection, we used the simplified final model with the lowest QAICc as the null model to test for additive linear and/or quadratic effects of beak width and beak height on apparent survival. We considered three possible models to test the effect of the individual covariate (beak width or height): only a linear effect (linear directional selection), only a quadratic effect (stabilising or disruptive selection around the population mean), and both effects (stabilising or disruptive selection around a value different from the population mean). Models including polynomial effects of the third order, i.e. cubic effects, did not improve model fit, so were not further considered. Results Model selection with the complete dataset The probability of sex assignment varied over time and was, on average, higher for males (0.83, 95%CI: 0.74-0.89) than for females (0.52, 95%CI: 0.48-0.56). The probability of recapture was independent of time and sex, and low (0.05, 95% CI: 0.040.07). There was a slight but non-significant majority of males among first-time captured individuals (0.53, 95% CI: 0.49-0.56). We also did not find temporal variation ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le in the initial state, i.e. the probability that a first-time captured individual was a male (0.47, 95%CI: 0.44-0.51). For apparent survival, the best-supported model included the additive effects of age, sex and time (see Table 1 and Fig. 1). On average (model 11 in Table 1), juveniles had the lowest survival (females, 0.30, 95% CI: 0.23-0.39; males, 0.33, 95% CI: 0.25-0.42), followed by adults (females, 0.45, 95% CI: 0.39-0.50; males, 0.48, 95% CI: 0.43-0.53) and yearlings (females, 0.55, 95% CI: 0.41-0.67; males, 0.58, 95% CI: 0.44-0.70). The effects of beak size on apparent survival The model selection procedure with the two reduced datasets (individuals with beak measures available) gave similar results as for the whole dataset, except for (due to reduced sample size) a lack of a sex effect on apparent survival and on sex assignment (Table A1). Subsequent testing of the additive effect of beak size on apparent survival led to different results for beak width versus beak height (Table 2). For beak width, the best model included the linear and quadratic effect of the covariate (see Fig. 2) (coefficient on logit scale for the linear term = 0.22; 95% IC -0.07 – 0.50; SE = 0.14, coefficient for the quadratic term on logit scale = -0.35; 95% IC -0.70 – 0.01; SE = 0.18) closely followed by the model only including the quadratic effect (ΔQAIC = 0.19; coefficient = -0.33; 95% IC -0.66 – 0.00; SE = 0.17). The optimal beak width with respect to apparent survival for crossbills in the study area was estimated at 11.43 mm (Fig. 2), with the probability of apparent survival decreasing away from this optimum. ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le For beak height, the best model did not include the covariate, but it had only a moderate model QAIC weight (0.45). Overall, models including a quadratic effect of beak height (indicating stabilising selection) obtained a cumulative model QAIC weight of 0.37, and for models including the linear effect this was 0.26. Discussion Does beak morphology affect survival? Our first objective was to investigate if crossbill beak size affects apparent survival, and if so, to determine which kind of beak is locally adaptive. Our findings suggest that beak morphology (beak width) of crossbills indeed is related to their apparent survival in the study population. The estimated values of apparent survival rates should be interpreted with caution as they could suffer some bias due to low recapture rate (Hargrove and Borland 1994, O’Brien et al. 2005). However, low recapture rate does not affect the effects of covariates (here: beak size) other than in reducing statistical power: it does not create linear or quadratic effects were there are none. This supports an earlier finding by Benkman et al. (2003) who showed, for a resident crossbill population in North America feeding on a different pine species, that the return rate of marked individuals depended on beak morphology (especially beak height). We estimated that the population optimum beak width in our population is 11.43 mm, with apparent survival decreasing away from the optimum. Crossbills utilising mountain pine, as they do in our study area, have the largest beaks among Spanish crossbills (Edelaar et al. 2012). The fact that the observed average beak (Supplementary material Appendix 3 Fig. A1) size is close (slightly smaller, see below) to the estimated locally optimal beak size, supports the idea that this crossbill population benefits from having ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le such a large beak, and is in that respect locally adapted. While we do not have data that links beak morphology to feeding performance (Benkman 1993), it is likely that the need to have such large bills is related to the very thick cone scales of mountain pine, requiring more force to separate them (Mezquida and Benkman 2010). Finally, we found that the probability of apparent survival varies according to individuals’ age. Apparent survival of yearlings and adults greatly overlap (Fig. 1), but juvenile birds have lower survival. This suggests that, as in other species, juveniles experience higher mortality because of their inexperience (Kershner et al. 2004, Suedkamp Wells et al. 2007). Alternatively, juveniles might also be more prone to disperse than older individuals as has been described for several taxa (Matthysen 2012). Causes of stabilizing selection By using capture-recapture models, we separated recapture (the observation process; potentially confounded by individual traits, including morphology) from apparent survival (the underlying biology). Our second (and more challenging) objective was to further separate the possible roles of natural selection (mortality) versus selection of the environment (permanent emigration) in explaining the variation in the observed apparent survival as related to beak morphology. A pattern of apparent survival as found by us (Fig. 3) is generally interpreted as stabilising natural selection acting through selective mortality of locally maladapted individuals. The alternative, that locally maladapted individuals permanently leave the study area and therefore are accounted as “locally dead”, is far less often considered (if not generally ignored). Given that our estimates proceed from a uni-site modelling approach, it is not possible to directly discriminate between true mortality versus permanent emigration (as is true for most capture-recapture studies): both processes can yield the exact same pattern, since local ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le survivors need to survive and stay present in the study area. Indeed, previous studies on apparent survival in crossbills (Senar et al. 1993, Benkman 2003, Alonso and Arizaga 2012, Santisteban et al. 2012, Benkman et al. 2014, Benkman 2017) were all done in single study areas, and virtually none of them (with the exception of Benkman 2017) suggested permanent emigration to explain individual or annual variation in apparent survival rates. Despite the difficulty distinguishing between natural selection (here: selective mortality) and selection of the environment (here: selective permanent emigration), we discuss a few important indications that suggest to us that natural selection is not the only or more likely explanation. First, when comparing our survival estimate for quadratic selection (β = -0.36) with those compiled by Kingsolver and Diamond (2011) (Fig. 3), our estimate would be a rather unusually strong measure of stabilising natural selection. Said otherwise, an unusual number of selective deaths would have to occur because of a phenotypic trait. Therefore, it appears probable that selective dispersal of locally maladapted individuals out of the study area has also contributed in order to produce this high value. The same kind of observation and logic led Benkman (2017) to conclude the same, in that case, to explain a pattern of apparent directional selection for larger bills. Our parallel observation, based on a more correct capture-recapture methodology, reinforces his conclusions. Second, the estimated optimal beak width (11.43 mm) is somewhat higher than the observed average beak width (11.02 mm) in the population (Table 2 and Supplementary material Appendix 3 Fig. A1). We simulated (Supplementary material Appendix 2) how often this observed average beak width produced a lower survival rate than the predicted optimal beak width, based on the confidence intervals obtained from our capture-mark-recapture analyses. This occurs ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le in 87.6% of simulated cases. Hence, there is fairly strong support that the studied population has a mean beak width which is too small for what provides the highest apparent survival. To the extent that our population indeed has a smaller than optimal mean beak width, this is in line with the immigration of individuals with smaller beaks originating from other populations. Given that Spanish crossbills using other pines (Scots pine, Black pine, Aleppo pine) which are growing relatively nearby have smaller beaks on average (Senar et al. 1993, Edelaar et al. 2012, Borrás et al. 2008), it is likely that at least a proportion of the birds we captured were birds from other areas and thereby this reduced mean beak size. In view of the high mobility of crossbills, it is also likely that these birds subsequently would return to their original habitat or move on in search of better places. Such behaviour is in line with the observation of a transient effect in our recapture data (see Methods). We also explored whether such influxes into our study population could be due to immigration from northern Eurasia, when the availability of resources is limited (Davis 1964, Herremans 1988, Cramp and Perrins 1994, Summers et al. 1996, Edelaar and Terpstra 2004). However, in our study area, years with high numbers of birds caught did not coincide at all (data not shown) with years in which invasions appeared to be occurring in central Europe (based on Marquiss et al. 2012, and the websites Trektellen.org, Waarneming.nl, Euro Bird Portal and Vogelwarte.ch). This suggests that any influxes would be coming from more nearby, probably Spain itself. In view of these indications, and in view of what is known about the biology of crossbills (high mobility, ease for crossbills to estimate individual local adaptedness via food intake rate, ease of recognising different pine species from the air), we feel that it is hard to exclude that morphologically maladapted crossbills were caught and marked ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le by us, but subsequently left the study area permanently because of low food intake rates. How large this effect of selection of the environment on apparent survival is relative to natural selection cannot be established in studies such as ours that have only local recapture data: it might be very small, but in the absence of actually knowing that non-recaptured individuals are dead it might also be the main contributor to the observed pattern of apparent survival. We feel that our study system is not unique in this respect: selective immigration or emigration associated with phenotype and ecological performance could occur in many organisms, as long as they have at least some control over their dispersal events. Biological knowledge and indirect evidence can help distinguish somewhat between selective dispersal and selective mortality, but ultimately different study designs (e.g. monitoring across ecologically heterogeneous areas) and data types (e.g. true deaths or radio-/satellite-tracking) are required for this. Data availability statement Data available from the Dryad Digital Repository: (XXXXX et al. 2019). Acknowledgements – We thank the people who helped to collect the data, especially to Xevi Colomè, Toni Cabrera and Jose Molina. Funding – This study was financially supported by the Ministry of Economy and Competitivity (CGL-2016-79568-C3-3-P to JCS and PECGL2013-49460-EXP and CGL2016-79483-P to PE), with support from the European Regional Development Fund FEDER. During the realization of this study, SS has been financed by a Juan de la Cierva formation grant (FJCI-2015-24579). ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le References Alonso, D., Arizaga, J., Miranda, R. and Hernández, M. Á. 2006. Morphological diversification of common crossbill Loxia curvirostra populations within Iberia and the Balearics. – Ardea 107: 94–99. Alonso, D. and Arizaga, J. 2012. The impact of vagrants on apparent survival estimation in a population of common crossbills (Loxia curvirostra). – J. Ornithol. 154: 209–217. Armsworth, P. R. and Roughgarden, J. E. 2008. The structure of clines with fitnessdependent dispersal. – Amer. Nat. 172: 648–657. Bailey, R. C. and Byrnes, J. 1990. A new, old method for assessing measurement error in both univariate and multivariate morphometric studies. – Syst. Biol. 39: 124– 130. Benkman, C. W. 1993. Adaptation to single resources and the evolution of crossbill (Loxia) diversity. – Ecol. Monogr. 63: 305–325. Benkman, C. W. 2003. Divergent selection drives the adaptive radiation of crossbills. – Evolution 57: 1176–1181. Benkman, C. W. 2016. The natural history of the south hills crossbill in relation to its impending extinction. – Amer. Nat. 188: 589–601. Benkman, C. W. 2017. Matching habitat choice in nomadic crossbills appears most pronounced when food is most limiting. – Evolution 71: 778–785. Benkman, C. W. and Lindholm, A. K. 1991. The advantages and evolution of a morphological novelty. Nature 349: 519–520. ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Benkman, C. W., Parchman, T. L., Favis, A. and Siepielski, A. M. 2003. Reciprocal selection causes a coevolutionary arms race between crossbills and lodgepole pine. – Amer. Nat. 162: 182–194. Benkman, C. W., Colquitt, J. S., Gould, W. R., Fetz, T., Keenan, P. C. and Santisteban, L. 2005. Can selection by an ectoparasite drive a population of red crossbills from its adaptive peak? – Evolution 59: 2025–2032. Berner, D. and Thibert-Plante, X. 2015. How mechanisms of habitat preference evolve and promote divergence with gene flow. – J. Evol. Biol. 28: 1641–1655. Björklund, M., Alonso, D. and Edelaar, P. 2013. The genetic structure of crossbills suggests rapid diversification with little niche conservatism. – Biol. J. Linn. Soc. 109: 908–922. Black, A. R. 1993. Predator induced phenotypic plasticity in Daphnia pulex: Life history and morphological responses to Notonecta and Chaoborus. – Limnol. Oceanogr. 38: 986–996. Bolnick, D. I. and Otto, S. P. 2013. The magnitude of local adaptation under genotypedependent dispersal. – Ecol. Evol. 3: 4722–4735. Borrás, A. and Senar, J. C. 2003. Piquituerto común Loxia curvirostra. – In: Martí, R. and del Moral, J. C. (ed.), Atlas de las aves reproductoras de España. Dirección General de Conservación de la Naturaleza-Sociedad Española de Ornitología, pp. 588–589. Borrás, A., Cabrera, J. and Senar, J. C. 2008. Local divergence between Mediterranean crossbills occurring in two different species of pine. – Ardeola 55: 169–177. Borrás, A., Senar, J. C., Cabrera, J. and Cabrera, A. 2011. Trencapinyes Loxia curvirostra. – In: Herrando, S., Brotons, L., Estrada, J., Guallar, S. and Anton, ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le M. (ed.), Atlas dels ocells de Catalunya a l'hivern 2006-2009. Institut Català d'Ornitologia and Lynx Edicions. Burnham, K. and Anderson, D. 2004. Model selection and multimodel inference. – Springer New York. Castroviejo, S., Laínz, M., López-Gonzalez, G., Montserrat, P., Muñoz-Garmendia, F., Pauci, J. and Villar, L. 1986. Flora Ibérica. Plantas vasculares de la Península Ibérica e Islas Baleares. – Real Jardín Botánico, C.S.I.C. Choquet, R., Rouan, L. and Pradel, R. 2009a. Program E-SURGE: a software application for fitting multievent models. – In: Thomson, D. L., Cooch, E. G. and Conroy, M. J. (ed.), Modeling demographic processes in marked populations, pp. 845–865. Choquet, R., Lebreton, J. D., Gimenez, O., Reboulet, A. M. and Pradel, R. 2009b. UCARE: Utilities for performing goodness of fit tests and manipulating CAptureREcapture data. – Ecography 32: 1071–1074. Choquet, R. and Nogue, E. 2011. E-SURGE 1.8 User’s Manual. – CEFE, Montpellier. Clouet, M. 2000. The breeding biology of the common crossbill Loxia curvirostra in the Central Pyrenees. – Bird Study 47: 186–194. Cramp, S. and Perrins, C. M. 1994. Handbook of the birds of Europe, Middle East and North Africa. Volume VIII. – Oxford Univ. Press. Darwin, C. 1859. On the origin of species by means of natural selection, or the preservation of favoured races in the struggle for life. – John Murray, London. Davis, P. 1964. Crossbills in Britain and Ireland in 1963. – Br. Birds 57: 477–501. Del Val, E., Negro, J. J., Garrido-Fernández, J., Jarén, M., Borràs, A., Cabrera, J. and Senar, J. C. 2014. Seasonal variation of red carotenoid pigments in plasma of wild crossbill males Loxia curvirostra. – J. Ornithol. 155: 211–218. ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Duriez, O., Sæther, S. A., Ens, B. J., Choquet, R., Pradel, R., Lambeck, R. H. D. and Klaassen, M. 2009. Estimating survival and movements using both live and dead recoveries: A case study of oystercatchers confronted with habitat change. – J. Appl. Ecol. 46: 144–153. Edelaar, P. and Terpstra, K. 2004. Is the nominate subspecies or the common crossbill Loxia c. curvirostra polytypic? I. Morphological differences among years at a single site. – Ardea 92: 93–102. Edelaar, P., Siepielski, A. M. and Clobert, J. 2008. Matching habitat choice causes directed gene flow: a neglected dimension in evolution and ecology. – Evolution 62: 2462–2472. Edelaar, P., Alonso, D., Lagerveld, S., Senar, J. C. and Björklund, M. 2012. Population differentiation and restricted gene flow in Spanish crossbills: not isolation-bydistance but isolation-by-ecology. – J. Evol. Biol. 25: 417–430. Edelaar, P. and Bolnick, D. I. 2012. Non-random gene flow: an underappreciated force in evolution and ecology. Trends Ecol. Evol. 27: 659–665. Edelaar, P., Jovani, R. and Gomez-Mestre, I. 2017. Should I change or should I go? Phenotypic plasticity and matching habitat choice in the adaptation to environmental heterogeneity. – Am. Nat. 190: 506–520. Edelaar, P. and Bolnick, D. I. 2019. Appreciating the multiple processes increasing individual or population fitness. – Trends Ecol. Evol. 34: 435–446. Endler, J. A. 1986. Natural selection in the wild. – Princeton Univ. Press. Genovart, M., Pradel, R. and Oro, D. 2012. Exploiting uncertain ecological fieldwork data with multi-event capture-recapture modelling: An example with bird sex assignment. J. Anim. Ecol. 81: 970–977. ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Ghalambor, C. K., McKay, J. K., Carroll, S. P. and Reznick, D. N. 2007. Adaptive versus non-adaptive phenotypic plasticity and the potential for contemporary adaptation in new environments. – Funct. Ecol. 21: 394–407. Groth, J. G. 1993. Evolutionary differentiation in morphology, vocalizations, and allozymes among nomadic sibling species in the North American red crossbill (Loxia curvirostra) complex. –Univ. of California Press. Gómez Blanco, D., Santoro, S., Borrás, A., Cabrera, J., Senar, J. C. and Edelaar, P. 2019. Data from: Beak morphology predicts apparent survival of crossbills: due to selective survival or selective dispersal? – Dryad Digital Repository, . Hargrove, J. W. and Borland C. H. 1994. Pooled population parameter estimates from mark-recapture data. – Biometrics 50: 1129–1141. Hendry, A. P., Day, T. and Taylor E. B. 2001. Population mixing and the adaptive divergence of quantitative traits in discrete populations: a theoretical framework for empirical tests. – Evolution 55: 459–466. Herremans, M. 1988. Measurements and moult of irruptive common crossbills (Loxia curvirostra) in Central Belgium. – Gerfaut 78: 243–260. Holt, R. D. and Barfield, M. 2008. Habitat selection and niche conservatism. – Isr. J. Ecol. Evol. 54: 295–309. Kershner, E. L., Walk, J. W. and Warner, R. E. 2004. Postfledging movements and survival of juvenile eastern meadowlarks (Sturnella magna) in Illinois. – Auk 121: 1146–1154. Kingsolver, J. G. and Diamond, S. E. 2011. Phenotypic selection in natural populations: What limits directional selection? – Am. Nat. 177: 346–357. ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Lebreton, J., Burnham, K. P., Clobert, J. and Anderson D. R. 1992. Modeling survival and testing biological hypotheses using marked animals: A unified approach with case studies. – Ecol. Monogr. 62: 67–118. Lessells, C. M. and Boag, P. T. 1987. Unrepeatable repeatabilities: a common mistake. – The Auk, 104: 116–121. Manly, B. F. J. 1985. The statistics of natural selection on animal populations. – Springer Netherlands, Dordrecht. Marquiss, M. and Rae, R. 2002. Ecological differentiation in relation to bill size amongst sympatric, genetically undifferentiated crossbills Loxia spp. – Ibis 144: 494–508. Marquiss, M., Hobson, K. A. and Newton, I. 2008. Stable isotope evidence for different regional source areas of common crossbill Loxia curvirostra irruptions into Britain. – J. Avian Biol. 39: 30–34. Marquiss, M., Newton, I., Hobson, K. A. and Kolbeinsson, Y. 2012. Origins of irruptive migrations by common crossbills Loxia curvirostra into northwestern Europe revealed by stable isotope analysis. – Ibis 154: 400–409. Matthysen, E. 2012. Multicausality of dispersal: a review. – In: Clobert, J. M., Baguette, M., Benton, T. G. and Bullock, J. M. (ed.), Dispersal ecology and evolution, pp. 3–18. Oxford Univ. Press. Mezquida, E. T. and Benkman, C. W. 2010. Habitat area and structure affect the impact of seed predators and the potential for coevolutionary arms races. – Ecology 91: 802–814. Newton, I. 2006. Movement patterns of common crossbills Loxia curvirostra in Europe. – Ibis 148:782–788. ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le O’Brien, S., Robert, B. and Tiandry, H. 2005. Consequences of violating the recapture duration assumption of mark-recapture models: a test using simulated and empirical data from an endangered tortoise population. – J. Appl. Ecol. 42: 1096–1104. Parchman, T. L., Buerkle, C. A., Soria-Carrasco, V. and Benkman, C. W. 2016. Genome divergence and diversification within a geographic mosaic of coevolution. – Mol. Ecol. 25: 5705–5718. Parchman, T. L., Edelaar, P., Uckele, K., Mezquida, E. T., Alonso, D., Jahner, J. P., Summers, R. W. and Benkman C. W. 2018. Resource stability and geographic isolation are associated with genome divergence in western Palearctic crossbills. – J. Evol. Biol. 31: 1715–1731. Pradel, R. 2005. Multievent: an extension of multistate capture-recapture models to uncertain states. – Biometrics 61: 442–447. Pradel, R. 2009. The stakes of capture–recapture models with state uncertainty. – In: Thomson, D. L., Cooch, E. G. and Conroy, M. J. (ed.), Modeling Demographic Processes in Marked Populations, pp. 781–795. Pradel, R., Maurin-Bernier, L., Gimenez, O., Genovart, M., Choquet, R. and Oro, D. 2008. Estimation of sex-specific survival with uncertainty in sex assessment. – Can. J. Stat. 36: 29–42. Przybylo, R., Sheldon, B. C. and Merilä, J. 2000. Climatic effects on breeding and morphology: Evidence for phenotypic plasticity. – J. Anim. Ecol. 69: 395–403. Santisteban, L., Benkman, C. W., Fetz, T. and Smith, J. W. 2012. Survival and population size of a resident bird species are declining as temperature increases. – J. Anim. Ecol. 81: 352–363. ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Santoro, S., Green, A. J. and Figuerola, J. 2016. Immigration enhances fast growth of a newly established source population. Ecology 97: 1048–1057. Scheiner, S. M. 2016. Habitat choice and temporal variation alter the balance between adaptation by genetic differentiation, a jack-of-all-trades strategy, and phenotypic plasticity. – Am. Nat. 187: 633–646. Senar, J. C. and Borrás, A. 1985. Crossbill irruptions. – BTO News 141: 6. Senar, J. C., Borrás, A., Cabrera, T. and Cabrera, J. 1993. Testing for the relationship between coniferous crop stability and common crossbill residence. – J. Field Ornithol. 64: 464–469. Siepielski, A. M. and Benkman, C. W. 2005. A role for habitat area in the geographic mosaic of coevolution between red crossbills and lodgepole pine. – J. Evol. Biol. 18: 1042–1049. Suedkamp Wells, K. M. S., Ryan, M. R., Millspaugh, J. J., Thompson, F. R. and Hubbard, M. W. 2007. Survival of postfledging grassland birds in Missouri. – Condor 109: 781–794. Summers, R. W., Jardine, D. C., Marquiss, M. and Proctor, R. 1996. The biometrics of invading common crossbills Loxia curvirostra in Britain during 1990–1991. – Ringing Migr. 17: 1–10. Summers, R. W., Dawson, R. J. G. and Phillips, R. E. 2007. Assortative mating and patterns of inheritance indicate that the three crossbill taxa in Scotland are species. – J. Avian Biol. 38: 153– 162. Svensson, L. 1992. Identification guide to european passerines. Fourth edition. – BTO. Supplementary material (Appendix JAV-02107 at ). Appendix 1–3. ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Figure Legends Figure 1 – Ninety-five % CIs for annual apparent survival rate over the study period and according to age class (juveniles=red squares; yearlings=blue circles; adults=black triangles). For graphical clarity, the (minor) differences between males and females are not represented (estimates from model 12 in Table 1). ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Figure 2 – Ninety-five % CIs for model-averaged annual apparent survival rate according to beak width (mm) and age class (a, juveniles; b, yearlings; c, adults). ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Figure 3 – Probability densities for estimated quadratic selection gradients from the present study versus the literature. The light blue curve represents the probability density according to our survival estimate for quadratic selection (β = -0.36, 95% credible intervals: -0.70 – -0.01). The grey curve is adapted from the review by Kingsolver and Diamond (2011). ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Table Legends Table 1. Model selection with the complete dataset. Abbreviations: np, number of parameters; QAICc, Quasi-Akaike Information Criterion corrected for over-dispersion; ΔAICc, AICc differences between models; wAICc, Akaike weight (support of the current model with respect to the candidate set of models); t, time effect; a, age effect, c, constant effect. The best model selected for each model selection step is in bold. No. Initial Sex Assignment Capture State Survival np Deviance QAICc ΔAICc wAICc Modeling sex assignment probability 1 sex . t sex . t t a . sex + t 164 17185.35 17521.13 9.82 0.00 2 sex + t sex . t t a . sex + t 138 17229.81 17511.31 0.00 0.00 3 t sex . t t a . sex + t 137 17247.60 17527.02 15.71 0.00 4 sex sex . t t a . sex + t 112 17845.46 18073.07 561.76 0.00 5 c sex . t t a . sex + t 111 18004.79 18230.34 719.03 0.00 Modeling capture probability 2 sex + t sex . t t a . sex + t 138 17229.81 17511.31 39.04 0.00 6 sex + t sex + t t a . sex + t 113 17261.36 17491.04 18.78 0.00 7 sex + t t t a . sex + t 112 17262.10 17489.72 17.45 0.00 8 sex + t sex t a . sex + t 88 17295.70 17473.93 1.66 0.05 9 sex + t c t a . sex + t 87 17296.09 17472.27 0.00 0.12 Modeling Initial state probability 9 sex + t c t a . sex + t 87 17296.09 17472.27 0.15 0.12 10 sex + t c c a . sex + t 61 17349.05 17472.12 0.00 0.12 Modeling survival probability 10 sex + t c c a . sex + t 61 17349.05 17472.12 2.37 0.12 11 sex + t c c a + sex + t 59 17350.75 17469.75 0.00 0.41 12 sex + t c c a+t 58 17353.45 17470.42 0.66 0.29 13 sex + t c c a . sex 36 17440.00 17512.37 42.62 0.00 14 sex + t c c a + sex 34 17444.00 17512.34 42.59 0.00 15 sex + t c c a 33 17445.16 17511.48 41.72 0.00 ‘This article is protected by copyright. All rights reserved.’ Accepted Ar tic le Table 2. Model selection with the effect of (a) beak width and (b) beak height on apparent survival. Abbreviations: np, number of parameters; QAICc, Quasi-Akaike Information Criterion corrected for over-dispersion; ΔAICc, AICc differences between models; wAICc, Akaike weight (support of the current model with respect to the candidate set of models); t, time effect; a, age effect; c, constant effect; ß, linear effect of the covariate; ß2, quadratic effect of the covariate. The best model selected for each model selection step is in bold. (a) Beak width No. Sex Assignment Capture Initial State Survival np Deviance QAICc ΔAIC wAIC 1 c c a+t 42 2556.43 2644.09 2.76 0.11 c (a + t) + ß 43 2554.55 2644.38 3.04 0.09 44 2549.32 2641.33 0.00 0.42 43 2551.69 2641.53 0.19 0.38 t Modeling beak width effect 5 6 t c t 7 t (b) c c (a + t) + ß + ß 2 2 c c (a + t) + ß Beak height No. Sex Assignment Capture Initial State Survival np Deviance QAICc ΔAIC wAIC 1 c c a+t 50 5107.70 5210.25 0.00 0.45 c (a + t) + ß 51 5107.70 5212.35 2.10 0.16 52 5106.50 5213.26 3.01 0.10 51 5106.63 5211.28 1.04 0.27 t Modeling beak height effect 2 3 4 t t t c c c c c (a + t) + ß + ß (a + t) + ß 2 2 ‘This article is protected by copyright. All rights reserved.’