The genetic regulation of avian migration timing: combining candidate genes and quantitative genetic approaches in a long-distance migrant

Plant and animal populations can adapt to prolonged environmental changes if they have sufficient genetic variation in important phenological traits. The genetic regulation of annual cycles can be studied either via candidate genes or through the decomposition of phenotypic variance by quantitative genetics. Here, we combined both approaches to study the timing of migration in a long-distance migrant, the collared flycatcher (Ficedula albicollis). We found that none of the four studied candidate genes (CLOCK, NPAS2, ADCYAP1 and CREB1) had any consistent effect on the timing of six annual cycle stages of geolocator-tracked individuals. This negative result was confirmed by direct observations of males arriving in spring to the breeding site over four consecutive years. Although male spring arrival date was significantly repeatable (R = 0.24 ± 0.08 SE), most was attributable to permanent environmental effects, while the additive genetic variance and heritability were very low (h2 = 0.03 ± 0.17 SE). This low value constrains species evolutionary adaptation, and our study adds to warnings that such populations may be threatened, e.g. by ongoing climate change.


Introduction
Rising spring temperatures have advanced most avian annual cycle events, but the magnitude of this recent change differs between traits and species (Rubolini et al. 2007). For instance, the timing of breeding has advanced more than spring migration (Both and Visser 2001;Weidinger and Král 2007;Pearce-Higgins and Green 2014;Tomotani et al. 2018), and the latter has advanced more in short-distance migrants compared to long-distance migrants (Usui et al. 2017;Lehikoinen et al. 2019;Koleček et al. 2020). Such observations suggest that phenological shifts are often driven by plastic responses of birds to local temperatures and vegetation phenology-environmental cues that are more available to sedentary species and short-distance migrants (Helm et al. 2013;Clark et al. 2014). In contrast, temperatures at the non-breeding sites of long-distance migrants are not very indicative of phenology at their breeding sites (but see Saino and Ambrosini 2008). Therefore, these species largely rely on their genetically based internal circannual pacemaker that tells them when to start the migration (Gwinner 1996;Helm et al. 2013;Åkesson et al. 2017).
Long-distance migrants can also modify the speed of migration according to temperatures or rainfalls along their route (Ahola et al. 2004;Hüppop and Winkel 2006;Both 2010;Haest et al. 2020), but their arrival to a breeding site is largely determined by the onset of migration (Stanley et al. 2012;Ouwehand and Both 2017;Briedis et al. 2019;Schmaljohann 2019). Thus, unlike rapid plastic responses by short-distance migrants, the advanced phenology of long-distance migrants may be mainly driven by evolutionary changes to their internal clocks (Pulido 2007a; Communicated by Suvi Ruuskanen. Schmaljohann 2019). Both plastic changes and evolutionary responses help birds to cope with changing environmental conditions. Plastic responses are faster and thus may serve as a "first aid" (Charmantier et al. 2008;Vedder et al. 2013), but because plasticity has its own costs and limits (DeWitt et al. 1998;Auld et al. 2010), ultimately an evolutionary response is needed for the long-term sustainability of a population (Gienapp et al. 2013;Gienapp and Brommer 2014).
One prerequisite for an evolutionary response to a changing environment, e.g. the current rapid climate change, is a sufficient amount of genetic variation upon which selection may act. Laboratory experiments have demonstrated rapid evolutionary responses in the timing of migration to artificial selection (Pulido et al. 2001). However, such findings need to be verified in the wild, where the environmental variation may mask genotypes from the action of selection. In other words, we need to know the amount of additive genetic variation scaled by natural amounts of phenotypic variation, i.e. heritability of the trait measured in the field. The heritability of breeding times has been often studied, and found to vary both between and within wild populations (reviewed by Liedvogel et al. 2012). In contrast, only a handful of studies have considered the heritability of migration timing. These have usually found substantial heritability, suggesting a high evolutionary potential for this phenological trait in wild populations (Potti 1998;Møller 2001;Teplitsky et al. 2011;Arnaud et al. 2013;Tarka et al. 2015). This was recently confirmed in a population of the pied flycatcher (Ficedula hypoleuca), which considerably advanced its internal annual rhythm in as few as 21 years (Helm et al. 2019).
Although heritability is a useful measure for estimating the possible rate of adaptation to environmental change, it tells nothing about proximate mechanisms underlying this genetic adaptation. To reveal how the genetic variation is translated to differences between phenotypes that are visible to selection, one needs to identify the involved genes and the function of their protein products. One possibility is to adopt a candidate gene approach, which looks for orthologs of genes with known function in genetic models, like Drosophila or mice, in the studied non-model organism (Fitzpatrick et al. 2005). However, these genetic models are usually short-lived, which prevents studying circannual rhythms on them. Fortunately, an emerging view suggests that circannual rhythms are closely linked to circadian rhythms, the latter having a very well described molecular basis (Dunlap 1999;Bell-Pedersen et al. 2005). Steinmeyer et al. (2009) found several circadian genes (CLOCK, NPAS2, ADCYAP1, CREB1) with allelic variation caused by a short repetitive sequence (microsatellite) in the blue tit (Cyanistes caeruleus). CLOCK and its paralog NPAS2 are part of the core circadian oscillator (Reick et al. 2001;DeBruyne et al. 2007). CREB1 is best known for its role in memory formation and learning (Silva et al. 1998), while ADCYAP1 is a regulator of metabolism including digestion, respiration and immune function (Sherwood et al. 2000). However, both these genes are also involved in the regulation of circadian rhythms (Hannibal et al. 1997;Kako and Ishida 1998;Lonze and Ginty 2002;Nagy and Csernus 2007). For brevity, we will label these four genes as "circadian", but it should be borne in the mind that some of them have many additional functions.
Along with the regulation of daily rhythms, these four circadian genes are also suspected to have roles in the regulation of annual cycles (Visser et al. 2010;Liedvogel et al. 2011;Helm et al. 2013;Merlin and Liedvogel 2019). Correlations between individual circadian genotypes and the timings of breeding (Liedvogel et al. 2009;Caprioli et al. 2012;Bourret and Garant 2015), moult (Saino et al. 2013), and migration (Mueller et al. 2011;Peterson et al. 2013;Bazzi et al. 2015;Ralston et al. 2019) were found in several avian populations. In contrast, other studies did not find such within-population associations between individual genotypes and the timing of either migration (Contina et al. 2018;Parody-Merino et al. 2019) or breeding (Liedvogel and Sheldon 2010;Dor et al. 2012;Chakarov et al. 2013). While studies of migration are not usually very powerful due to the difficulty of tracking birds along the route, the large sample sizes of the studies of breeding phenology make their negative findings robust. Thus, evidence for the effects of circadian genes on avian phenology is mixed, and this topic warrants further study (Ralston et al. 2019).
To add to the knowledge on whether circadian genes also regulate circannual cycles, we tested the association between the circadian genotype at four candidate genes (CLOCK, NPAS2, ADCYAP1, CREB1) and the phenology of migration in the collared flycatcher (Ficedula albicollis), a nocturnal long-distance migrant. First, we determined the timing of all stages in the annual cycle of the flycatchers by tracking with light-level geolocators. Second, we directly recorded the spring arrivals of individual birds to their breeding site over four consecutive years. Our set of candidate genes likely represents only a small fraction of all genes involved in the regulation of avian phenology. Consequently, negative results of these tests would not conclusively prove a minor role of genes in the regulation of annual cycles. Therefore, we also employed quantitative genetic methods and decomposed phenotypic variance in the timing of spring arrival to its causal components (see also Liedvogel et al. 2012). A high additive genetic component combined with a low effect of the investigated candidate genes would indicate that other genes retain high evolutionary potential in the studied population. A finding of a low heritability of migration phenology, on the other hand, would add to warnings that ongoing climate change may threaten populations of long-distance migrants (Radchuk et al. 2019;Koleček et al. 2020).

Study populations and fieldwork
We conducted this study in two close (35 km) nestbox collared flycatcher populations in Czechia. Study plots were located in a mixed forest dominated by 40% beech (Fagus sylvatica) and 40% oak (Quercus petraea) at the site Dlouhá Loučka site (49° 50′ N, 17° 13′ E, 340-500 m asl) and in a > 90% oak forest (Quercus petraea) with some pines (Pinus silvestris) at the site Velký Kosíř (49°32′N, 17° 03′ E, 300-400 m asl). The collared flycatcher is a long-distance nocturnal migrant with non-breeding residency areas about 7000 km away from the European breeding sites (Briedis et al. 2016). Non-breeding areas of Dlouhá Loučka and Velký Kosíř breeding populations are overlapping in southern Africa including Democratic Republic of Congo, Angola, Zambia and Botswana (Briedis et al. 2018). Any potential small-scale differences between non-breeding areas of these two breeding populations cannot be distinguished due the inherit limits of geolocation accuracy (Lisovski and Hahn 2012). The species easily adopts nest-boxes for breeding, and each of our populations hosts about 100 pairs of collared flycatcher. Other species breeding in our nest-boxes are the great tit (Parus major), blue tit and nuthatch (Sitta europaea). Nest-boxes are attached to trees about 160 cm above ground, and have inner dimensions of 22.5-25.5 × 11 × 12 cm (height × width × depth) and a nest entrance diameter of 32 mm.
We recorded the spring arrival of male collared flycatchers to the Velký Kosíř site from 2013 to 2016, taking advantage of their habit of exploring potential cavities for breeding immediately after their arrival to a breeding site. We equipped all nest-boxes with traps, and caught newly arriving males at intervals from 1 to 7 (mean = 2.55) days during the early part of the breeding season. Shorter intervals between catching sessions were used at the peak of male arrival, while longer ones were used at times when only a few males arrived. We used the mean time between the first capture of a male and the previous catching session as an estimate of the arrival date. These estimates were highly correlated with the true arrival dates known from geolocators (r = 0.95, n = 16, Supplementary Fig. 1). Each male was measured and blood sampled by tibial venipuncture, and its age was determined from plumage and categorized as young (in the second year of life) or old (in the third year of life or older). Note that females cannot be efficiently caught by this method since it is the male who enters the cavity first as part of its display to females. For further details about male catching see the Supplementary material.

Geolocator data
We used light-level geolocators to study the timing of six migratory stages in collared flycatchers (departure from the breeding site, crossing the Sahara in autumn, arrival to the non-breeding site, departure from the non-breeding site, Sahara crossing in spring, and arrival to the breeding site). We deployed 69 geolocators (model GDL2.0 with 7 mm light stalk, Swiss Ornithological Institute) on 33 adult males and 36 adult females at Dlouhá Loučka when they took care of their broods during the late nestling stage in 2013. At Velký Kosíř, we equipped adults with geolocators upon their spring arrival (139 males) or during the nestling stage (18 males and 8 females) in 2014. Each geolocator weighed about 0.6 g (ca 4.6% of the adult body mass). Traditionally, geolocators weighing up to 5% of the adult mass have been considered to cause little trouble to birds wearing them (Portugal and White 2018), and this has recently been confirmed by an extensive meta-analysis (Brlík et al. 2020).
We retrieved 29 geolocators at Dlouhá Loučka and 30 at Velký Kosíř in the following years, but our final sample size for the timing of annual cycle events varied from 26 (arrival to the breeding site) to 41 (Sahara crossing in autumn and arrival to the non-breeding site) due to the failure rate of the devices. The lower retrieval of loggers at Velký Kosíř site (30/165 = 18%) compared to Dlouhá Loučka site (29/69 = 42%) was likely caused by different methods of geolocator deployment. At Dlouhá Loučka, we deployed loggers on birds breeding in nest-boxes in forest patches that we knew to be very attractive to flycatchers. Birds breeding in such high-quality forest patches were likely to be in superior condition and have high fidelity to these attractive breeding sites. On the contrary, most loggers at Velký Kosíř site were deployed on birds before they established territories. Many of these birds emigrated from our plots after deployment of the loggers-not necessarily due to this manipulation-and bred in natural cavities or became floaters. They also had a smaller likelihood of surviving to the next breeding season simply because of the longer time interval. These methodological as well as actual site differences should not adversely affect our results because we statistically controlled for site differences in all analyses (see below).
Geolocators recorded ambient light intensity on an arbitrary scale ranging from 0 to 63 units, corresponding to 0 and approximately 3500 lx, respectively. We used a threshold approach (Lisovski et al. 2020) to estimate individual migration timing from the recorded light data. First, we identified sunrise and sunset times within the light-level data using 'GeoLocator' software (Swiss Ornithological Institute) and setting the light intensity threshold to 1 unit on the arbitrary scale. Second, we determined the start and end of migratory periods within each dataset using the 'changeLight' function (parameters: q = 0.85 and minimum duration of stationary/ stopover periods = 3 days) of the R package 'GeoLight' v 1.03 (Lisovski and Hahn 2012). Timing of the Sahara crossing was determined by a manual inspection of raw daily light recordings to identify days when prolonged periods of uninterrupted maximal light intensities (63 arbitrary units) were recorded. Such patterns of uninterrupted maximal light recordings are conspicuous in the raw light-level data and are characteristic of non-stop diurnal flights when birds cross large ecological barriers, like the Sahara Desert (Adamík et al. 2016, including data from the collared flycatcher). We used the first day of Sahara crossing (the total duration of the crossing is 2-3 days in most cases) in all calculations.

Candidate genes
DNA was extracted from blood using DNeasy (r) Blood & Tissue kits (Qiagen). We genotyped flycatchers at four candidate genes (CLOCK, NPAS2, ADCYAP1 and CREB1) for circadian/circannual cycles (Steinmeyer et al. 2009). Primers for the amplification of NPAS2 and ADCYAP1 were the same as in Steinmeyer et al. (2009). CREB1 was amplified using the forward primer from Steinmeyer et al. (2009) and a modified reverse primer, AGA ATA ACG CAG CCC AGA GC, from Bourret and Garant (2015). CLOCK primers were adopted from Caprioli et al. (2012). Forward primers for the amplification of CLOCK, NPAS2, ADCYAP1 and CREB1 were labelled with the dyes 6FAM, VIC, PET and NED, respectively. All loci were amplified in a single multiplex PCR using Type-it ® Microsatellite PCR kits (Qiagen) following the manufacturer's protocols. Annealing temperature was 53 °C. PCR products were mixed with a GeneS-canTM-500 LIZ ® Size Standard (Applied Biosystems) and their size was resolved using fragment analysis in a 3130xl Genetic Analyzer (Applied Biosystems). Genotypes were scored with GeneMarker ® 1.9 (Softgenetics).

Pedigree
At the Velký Kosíř site, most nestlings and adults were ringed each year since 1998. This provides an extensive social pedigree for this population. In some years, paternity analysis was conducted on parts of the population (2001-2002, Krist et al. 2005-2009, Krist and Munclinger 2011 or on the whole population (2013, Edme et al. 2017). There was no case of maternal error, so all maternal links in the pedigree should be accurate. In total, extra-pair paternity was detected at a rate of 23.5%, which is comparable to other populations of this species (Hungary: 20.6%, Rosivall et al. 2009;Sweden: 15%, Sheldon and Ellegren 1999). We corrected the social pedigree using the genetic information, which was available for 18% of the offspring born between 1998 and 2015 (Edme et al. 2019).
Consequently, after this correction, less than 20% of paternal links probably remained erroneous, giving in a total of less than 10% errors in the full pedigree (maternal and paternal links together). This inadequacy should not have excessive impact on our ability to estimate quantitative genetic parameters, since pedigree-based animal models are robust to even higher rates of paternity errors. A simulation study showed that 40% of paternal errors may cause an underestimation of heritability by 20% (Charmantier and Réale 2005). Similarly, Firth et al. (2015) reported a 15% underestimation of heritability with 12.5% extra-pair paternity and non-random mating. Although these results highlight the need for caution when interpreting heritabilities from systems with frequent paternal errors, they also suggest that the underestimation of heritability is not overwhelming.

Data analyses
Allele frequencies were calculated in Cervus 3.0 (Kalinowski et al. 2007). All other statistical analyses were done in R 3.6.2 (R Core Team 2019). We used mean allele length as the main variable describing individual genotype (see also Liedvogel et al. 2009;Mueller et al. 2011;Bourret and Garant 2015). This approach relies on two assumptions. First, that there is not much dominance on the locus. This assumption is supported by observations that phenotypic variation in complex traits such as behaviour is usually more affected by additive compared to dominant or epistatic effects (Hill et al. 2008;Wolak and Keller 2014). This was also the case of CLOCK allele length and the timing of breeding in the blue tit (Liedvogel et al. 2009). Second, this method assumes a linear effect of the genotype on the phenotype. Support for this assumption comes from studies that have found latitudinal clines in circadian genotypes (Johnsen et al. 2007, O'Malley et al. 2010, reviewed in Kyriacou et al. 2008). An alternative approach of using exact genotype as a factor in statistical analysis is free of these assumptions. However, this second approach is less powerful due to the higher numerator degrees of freedom and inability to reveal possible trends with allele length. We used only mean allele length for tests based on small sample sizes of tracked individuals. We used both mean allele length and exact genotypes for analyses based on large sample sizes obtained by direct observations of spring arrival for loci with intermediate genetic diversity (CLOCK and ADCYAP1, Table 1). The NPAS2 locus had low genetic diversity (Table 1), and therefore it was coded only as a categorical factor and its mean allele length was not calculated. In contrast, the CREB1 locus had too high a genetic diversity to enable us to use exact genotype as a factor in analysis (Table 1). Therefore, we also tested the separate effects of the short and the long allele in addition to the main analysis based on the mean allele length at this locus.
The relationships between genotypes at three loci (CLOCK, ADCYAP1 and CREB1) and timing of the six core phases of the annual cycle inferred from geolocators were assessed and visualized in 18 (3 × 6) separate general linear models, with the day of year of the stage (1 = 1st January) as the dependent and mean allele length as the independent variables. We did not test for the effect of the NPAS2 gene because of its low allelic diversity and small sample size of tracked individuals. Three other independent variables were included in each model to reduce residual variation in the day of year (sex: male vs. female, age at deployment: young vs. old and site/year: Kosíř/2014 vs. Loučka/2013). Male age was determined from plumage characteristics while that of females from ringing data, as all females included in this study had been previously marked. Note that the factor site/ year comprises several different sources of variation-site characteristics (e.g. differences between elevations, forest types), methodological differences (deployment of loggers in pre-breeding vs. breeding stage) and years (2013/2014 vs. 2014/2015 season). By its inclusion into the model, we control all these sources of variation in the timing of annual cycles. However, as we are unable to separate these sources of variation from one another, we will not try to interpret the site/year factor even if it is statistically significant.
Direct observations of male spring arrivals to the breeding site were conducted in four consecutive years. To test for the potential effects of candidate genes, we first averaged values of repeatedly sampled individuals. Before averaging, year and age effects were removed from the data by means of using residuals from an ANOVA of arrival date on these two factors (Supplementary Table 1). Average residual date (response variable) was then predicted by the mean length of allele (CLOCK, ADCYAP1 and CREB1) and/or by the exact genotype (CLOCK, ADCYAP1 and NPAS2) locus in a linear model. The contribution of each data point (individual mean arrival date) was weighted by the number of observations (years) used for its calculation.
Phenotypic variance in male spring arrival dates from direct observations was decomposed to its causal components by fitting an animal model in the R package Asreml-R, version 3 (Butler 2009). Male spring arrival date was the dependent variable, age (young or old) the fixed independent variable. The random part of the model included four random effects: year, male identity not linked to the pedigree, male identity linked to the pedigree, and identity of the rearing nest. Male identity linked to the pedigree estimates the additive genetic component, male identity not linked to the pedigree estimates the permanent environment component, and identity of the nest estimates common environment effects (Wilson et al. 2010). All these random effects conveyed useful information, and the resulting model was not over-fitted as indicated by the singularity test (is Singular?) implemented in the lme4 package (Bates et al. 2015). Before entering the animal model, pedigree was pruned to contain only informative individuals (i.e. those that were either phenotyped or linked to at least two phenotyped individuals) in the R package Nadiv (Wolak 2012). The pruned pedigree retained 379 phenotyped males and their 669 relatives. Phenotyped males included, for example, 143 father-son pairs and 37 full-sibling groups raised in the same nest (31 pairs and 6 trios). Repeatabilities and their standard errors were calculated in the Rptr package (Stoffel et al. 2017).

Results
In total, 407 individuals that were either geolocator-tracked or had a record of their spring arrival to the Velký Kosíř breeding site were genotyped at the four candidate genes. CLOCK alleles varied by multiples of three base pairs. We sequenced CLOCK homozygotes with the length of alleles 120 bp, 123 bp, 126 bp and 129 bp. These alleles differed due to CAG repetition on its 5′ end and coded 10-13 glutamins in polyglutamine chain (polyQ) of the coded Clock protein. We therefore label them as Q 10 -Q 13 . ADCYAP1 and CREB1 were more variable compared to CLOCK, while NPAS2 was the most conservative locus (Table 1). There was an excess of homozygotes at the CREB1 locus, which together with the highest rate of amplification failure among loci suggests the presence of null alleles. A minor part of Table 1 Allele frequency, observed (HObs) and expected (Hexp) heterozygosity, the chi-square test of Hardy-Weinberg equilibrium and the estimated proportion of null alleles for four candidate genes (n = 407 individuals) The Hardy-Weinberg test was not performed for the NPAS2 gene because of its low allelic variability K number of alleles, nt individuals not successfully typed at the particular gene, F(null) estimated frequency of null alleles the allelic variation at this locus thus remained undetected and the results have to be interpreted with greater caution. Tracking flycatchers using geolocators revealed that males were ahead of females in most stages of the annual cycle (Table 2), although this difference was not statistically significant. Old individuals were somewhat delayed when compared to young ones (Table 2), but note that young birds in this case were already experienced in migration and were recorded on their second migration to and back from Africa using geolocators. Site/year had the strongest effect on the timing of migration. Birds from the Loučka population tracked during the 2013/2014 season migrated before those from the Kosíř population tracked during the 2014/2015 season (Table 2). In contrast to these consistent effects of covariates, the effects of candidate genes on the timing of annual events were relatively weak (Table 2, Fig. 1). We found only one significant relationship (earlier arrival of individuals with longer CLOCK alleles to their non-breeding sites, Table 2). However, this relationship was also not strong enough to remain significant after accounting for multiple testing (Rice 1989; p value boundary for statistical significance after sequential Bonferroni adjustment is 0.05/18 ≈ 0.003).
Catching of arriving males at the Kosíř breeding site revealed that this stage of the annual cycle was delayed by about a week in young compared to old males (Supplementary Table 1). Note that spring arrival to the breeding site was the final stage of the first migration of young males, which is in contrast to the data from geolocators that covered the second migration of these young birds (see above). The mean timing of spring arrival differed only slightly (up to 5 days, Supplementary Table 1) over the four studied years, although this result was statistically significant. In contrast, male spring arrival was unrelated to allele length or genotype for any of the candidate loci. CLOCK mean allele length: − 0.29 ± 0.24 (estimate ± SE), F 1,368 = 1.40, p = 0.237; CLOCK genotype: F 8,361 = 1.65, p = 0.109; NPAS2 genotype: F 3,366 = 0.83, p = 0.476; ADCYAP1 mean allele length: 0.30 ± 0.28, F 1,367 = 1.13, p = 0.288, ADCYAP1 genotype: F 12,356 = 0.71, p = 0.740; CREB1 mean allele length: − 0.03 ± 0.10, F 1,352 = 0.09, p = 0.761 (Fig. 2). The last relationship was closely similar if the response variable was the length of the short or the long allele instead of the mean length: CREB1 short allele length: − 0.03 ± 0.10, F 1,352 = 0.14, p = 0.715, CREB1 long allele length: − 0.01 ± 0.08, F 1,352 = 0.09, p = 0.864.
In total, we recorded 502 arrivals of 372 males to the Velký Kosíř breeding site by catching them in nest-boxes. 266 males were caught only in 1 year, but we had repeated records for 106 individuals (86 were caught in two years, 16 in three years and 4 in all four years of the study). These males often arrived at the breeding locality at different times in consecutive years (Fig. 3). Despite this variation, we also detected a modest but significant repeatability in the timing of their arrival (r = 0.24 ± 0.08 SE). We further decomposed the variance in male spring arrival to its causal components with the animal model. We found the repeatability to be explained mainly by permanent environmental effects (relative variance component: 0.194 ± 0.712 SE), while the additive genetic component, and consequently heritability, was much lower and insignificant (h 2 = 0.031 ± 0.170 SE, Table 3).

Discussion
We did not detect any strong relationships between genetic variation at the four circadian genes and the timing of any phase of the annual cycle in geolocator-tracked collared flycatchers. These negative findings, however, should be treated with caution as they are based on tens of individuals and thus the tests have small statistical power. Nevertheless, we found the same negative result in the dataset with hundreds of individuals, for which we directly observed their spring arrivals to the breeding site. Taken together, our results suggest that none of the genetic variation at any of the studied circadian genes plays an important role in circannual rhythmicity in this migratory species. However, we cannot exclude the possibility that epigenetic variation of these genes may affect the timing of migration (see Saino et al. 2017). Furthermore, our quantitative genetic analysis revealed only moderate repeatability (R = 0.24 ± 0.08 SE) and very small and insignificant heritability (h 2 = 0.03 ± 0.17 SE) in the timing of spring arrival to the breeding site. The true heritability of spring arrival may be somewhat higher due to paternal errors in our pedigree. However, after adjustment for this inadequacy (up to 20% increase of h 2 , see Charmantier and Réale 2005), the heritability of spring arrival was still very small (h 2 = 0.036) and suggests a low evolutionary potential for this important phenological trait.
Previous studies have sometimes confirmed associations between individual circadian genotype and phenology, for instance in the timing of breeding (Liedvogel et al. 2009;Caprioli et al. 2012), dispersal (Chakarov et al. 2013), migration (Mueller et al. 2011;Bazzi et al. 2015;Saino et al. 2017;Ralston et al. 2019) or moult (Saino et al. 2013;Bazzi et al. 2017). However, other studies have not found any such relationships (Dor et al. 2012;Liedvogel et al. 2012;Peterson et al. 2013;Contina et al. 2018;Romano et al. 2018;Parody-Merino et al. 2019, this study), or found them only in an interaction with environmental variables, such as breeding density (Bourret and Garant 2015). Many of these studies were based on small sample sizes because of the difficulty to obtain phenological data throughout the entire avian annual cycle (e.g. by tracking individuals using geolocators). Moreover, even those studies that Table 2 Results of general linear models relating the timing of six phenological stages as inferred from geolocators (day of the year, 1 = 1st January) to mean allele length at three candidate genes The Mean alelle length Fig. 1 Scatterplots showing the relationships between mean allele length at three candidate loci and the timing of six annual cycle events (day of the year, 1 = 1st January) determined from geolocators. BS breeding site, non-BS non-breeding site, Sahara timing of the Sahara crossing. Lines are driven from simple OLS regressions of the day of the year (response variable) on the length of the allele (predictor). No covariates were controlled in these models. For regression statistics controlled for effects of covariates see Table 2 reported significant associations between circadian genotypes and phenology usually did this only for some of the tested loci and only for a subset of the genetic traits that were considered (mean allele length, length of longer allele, length of shorter allele, heterozygosity, exact genotype and methylation level). On the other hand, only a limited number of loci that are involved in the regulation of circadian rhythms were examined in these field studies. Most research was done on the four microsatellite loci also considered in this study, as their allelic variation is easily determined by molecular methods. Less attention has been paid to circadian genes with a single-nucleotide polymorphism (e.g. AANAT, PERIOD2, CKIε; Steinmeyer et al. 2009), although their alleles also do not distinguish migratory and sedentary species (Ramos et al. 2017). Taken together, evidence for a causal link between the timing of daily and annual cycles is relatively weak (reviewed by Ralston et al. 2019;Parody-Merino et al. 2019). This is also in line with a recent study in the great tit, where genomic selection caused phenological shifts without any correlated change in the endogenous daily cycles (Verhagen et al. 2019). This could indicate that these two cycles are not genetically correlated (Verhagen et al. 2019) and that circadian genes cannot explain much of the phenological variation. However, migratory traits often have substantial genetic variation and heritability (Pulido 2007a). Consequently, there is a need to look for other loci that would explain the variation in key annual cycle events (Contina et al. 2018). This is currently achievable thanks to the development of next-generation sequencing methods (Stapley et al. 2010). Whole genome sequencing and high-density single-nucleotide polymorphism (SNP) chips can reveal sites of divergent selection and thus allelic variation between migratory phenotypes (Liedvogel et al. 2009;Merlin and Liedvogel 2019). The latter method has recently been employed in a study of Vermivora warblers whose migratory phenotypes correlated with SNPs at a new candidate gene, VPS 13A (Toews et al. 2019). The function of this gene is unknown in birds, but the authors hypothesize that it may be involved in the processing of metabolic products that arise when the warbler is migrating (Toews et al. 2019). Similarly, migratory phenotypes of the willow warbler (Phylloscopus trochilus) were most divergent at loci related to the metabolism of fatty acids (Lundberg et al. 2017) and those of Swainson's thrushes (Catharus ustulatus) at the loci TMEM192, with unknown function, and PALLD, with a function in cell signaling (Delmore et al. 2016). It is important to note that different migratory phenotypes came from different populations across a migratory divide. Thus, they represented a suite of intercorrelated traits, such as different breeding and non-breeding areas, in addition to different migratory patterns in space and time. It is therefore not clear whether these candidate genes may also affect the timing of avian migration. Nevertheless, support for this idea includes a recent finding that the laying date of the great tit was most strongly associated with another gene controlling metabolism, thyroglobulin, within a breeding population (Gienapp et al. 2017). Overall, these recent studies suggest that the timing of annual cycles may be more closely linked to metabolism than to circadian cycles. Therefore, it may now be time to look for new candidates in avian phenology.
Technological advances in tracking techniques like GPS and light-level geolocators provided an opportunity to repeatedly follow even small migrants through their whole annual cycle without severely compromising their fitness (Brlík et al. 2020). Repeated tracking has usually revealed a high individual consistency in the timing of migration. This suggests that the start of migration is controlled by an innate program (Alerstam et al. 2006;Lourenço et al. 2011;Vardanis et al. 2011;Stanley et al. 2012;Conklin et al. 2013;López-López et al. 2014;van Wijk et al. 2016). In contrast, birds are rather flexible in the selection of their route (Alerstam et al. 2006;Vardanis et al. 2011;Stanley et al. 2012;López-López et al. 2014;Klvaňa et al. 2018, Lisovski et al. 2021, although the opposite pattern of a consistent route and flexible timing has also been found (Hasselquist   (Both et al. 2016). In such a case, we would expect a decrease of repeatability along the route, with the smallest values at the time of arrival to the breeding ground (Both et al. 2016), as the progress of migration is dependent on environmental conditions en route and these often fluctuate between years (Both 2010;Briedis et al. 2017). A lower repeatability in the timing of migration near the final destination has been found in many studies (Alerstam et al. 2006;Lourenço et al. 2011;Vardanis et al. 2011;Sergio et al. 2014;Fraser et al. 2019, but see van Wijk et al. 2016). However, the absolute difference was often small and so the evidence for such lower repeatability is not particularly strong. Caution is also needed when comparing repeatability between studies (Conklin et al. 2013;Both et al. 2016). For example, lower reported repeatability for the arrival to the breeding site compared to the departure from the non-breeding site might be caused by the fact that the latter values are typically obtained on tracked birds while the former are most often derived by direct observations that are prone to larger measurement errors. On the other hand, studies using direct observations have often used much larger sample sizes than those of tracked birds (reviewed in Both et al. 2016). Our estimate of the repeatability of male spring arrival (r = 0.24 ± 0.08 SE) broadly fits within the values obtained by direct observations at breeding sites (Pulido 2007b;Both et al. 2016) and is actually very close to that of male pied flycatchers (Ficedula hypoleuca, r = 0.27 ± 0.03, Both et al. 2016). As repeatability is usually the upper limit to heritability (Falconer 1993, but see Dohm 2002, both these studies suggest a rather limited evolutionary potential for arrival date in Ficedula flycatchers. We further decomposed the variance in male arrival timing to its causal components and found very low additive genetic variance and heritability (h 2 = 0.03 ± 0.17 SE) and larger but still insignificant common and permanent environment components. Thus, the repeatable timing of the same individual is not much controlled by its genotype. Instead, it is possible that young birds learn how to migrate during their first migration and then repeat this successful tactic in the following years, a mechanism also suggested for common terns (Sterna hirundo; Arnaud et al. 2013) and black kites (Milvus migrans; Sergio et al. 2014). Alternatively, the timing of migration may also be somewhat pre-determined by early condition (e.g. Pulido 2007b), as suggested by the non-zero common environment component. However, these conclusions have to be treated with caution since our heritability estimate has wide standard errors despite having more than a thousand informative individuals in our pedigree. The difficulty of determining heritability in the wild with high precision is mirrored by the scarcity of quantitative genetic studies of avian migration.
We know of only five such studies (Potti 1998;Møller 2001;Teplitsky et al. 2011;Arnaud et al. 2013;Tarka et al. 2015). All of them tested the heritability of avian spring arrival to breeding grounds, as we did, and found heritability ranging from 0.10 to 0.54. Thus, our heritability estimate is the lowest so far. However, many more studies are needed to determine whether different heritability estimates have something to do with the biology of the studied species or if they are merely a product of sampling variation in space and time.
A low evolutionary potential in phenological traits, as we found for the spring arrival of collared flycatchers, would severely limit the rate of adaptive evolution and thus may have negative consequences for populations currently experiencing rapid climate change. On the other hand, a low consistency of individuals in their spring arrival suggests high plasticity in their decisions on when to arrive at the breeding locality, and this variation would also be adaptive if corresponding to the local phenology of vegetation. Unfortunately, we were unable to test for this scenario because our study spanned only four years. In general, phenotypic plasticity indeed helps populations to persist under climate change (Charmantier et al. 2008;Gienapp et al. 2013). However, the length of this persistence depends on the costs of phenotypic plasticity (Chevin et al. 2010) and species life history (Vedder et al. 2013). Adaptive evolution is necessary for the long-term persistence of populations under directional selection caused by global warming (Gienapp et al. 2013;Charmantier and Gienapp 2014;Radchuk et al. 2019). The evolved responses enable animals to appropriately time their life cycle events without the need to maintain costly physiological machinery that is necessary for high phenotypic plasticity (DeWitt et al. 1998;Auld et al. 2010). In addition, the cues upon which birds plastically react may be available to resident species or short-distance migrants, but not to long-distance migrants (Both and Visser 2001). Consequently, long-distance migrants may be unable to plastically adjust their spring arrival to match their breeding to the peak of food supply (Both et al. 2006), which could lead to population declines (Møller et al. 2008;Both et al. 2010;Koleček et al. 2020) and an increased risk of extinction (Radchuk et al. 2019). The application of quantitative genetic methods to a broader spectrum of species differing in life-histories, e.g. migration distance or generation time, would shed useful light on the mechanisms by which birds may adapt their phenology to ongoing climate change.