Introduction

Bird migration is one of the most fascinating and well-studied behaviours among birds, including work on the physiological and morphological adaptations required for successful migration and ecological correlates of this behaviour. Considerable interest has also focused on understanding how variation in the migratory phenotype is generated, with several studies demonstrating high variability for various migratory traits using selective breeding studies (e.g. Berthold et al. 1992), displacement experiments (e.g. Perdeck 1958; Chernetsov et al. 2008) and quantitative genetics analyses (e.g. Pulido and Berthold 2010). Nevertheless, the underlying genetic architecture shaping this behaviour remains poorly understood (Liedvogel et al. 2011; Delmore and Liedvogel 2016). One popular approach to enhance our understanding of the molecular basis of migratory traits has been the use of candidate genes for behavioural traits suggested to contribute to variation in migratory phenotype. These candidates are often selected by their molecularly characterised specific function in other (often model) organisms. A candidate gene approach attempts to identify an association between genetic variation in that gene and the focal phenotype, here our focal behaviour is migration. The underlying rationale of this approach is to focus on genetic variation in specific candidate regions of the genome that have been suggested to directly impact the function of the candidate gene and ultimately the target phenotype in other species as well. In the context of migration, the traits receiving most attention are circadian behaviour and personality traits (e.g. exploratory or anxiety-related behaviour); focal candidate genes include CLOCK, ADCYAP1, CREB1 and NPAS2 for circadian rhythm, and DRD4 and SERT for personality traits (see Müller et al. 2011 as one of the pioneer studies in the context of migration).

CLOCK and ADCYAP1 are among the most studied candidate genes in the field of bird migration. The overall pattern for CLOCK variability indicates a latitudinal cline in repeat lengths at the variable region, possibly reflecting local adaptation to seasonal variation at different latitudes. The circadian CLOCK gene is highly conserved among birds throughout most of its sequence with the exception of one C-terminal region that contains a variable poly-glutamine(Q) (poly-Q) motif with variability in the number of glutamine repeats both among and within species (e.g. Johnsen et al. 2007; Liedvogel et al. 2009; Bazzi et al. 2016; but also see Liedvogel and Sheldon 2010; Dor et al. 2011). This region has been suggested to influence the transcription activating potential of the protein, potentially altering rhythms in both physiology and behaviour. The aforementioned pattern of a latitudinal cline in lengths polymorphism has been recovered in several species and in a migration context has been suggested to reflect adaptive features related to migration. Similarly for ADCYAP1, a neuropeptide-coding gene, one 3′ UTR microsatellite length variation has shown associations with migratory behaviour in blackcaps (Müller et al. 2011; Mettler et al. 2015) with longer alleles associated with higher migratory activity. However, the pattern so far lacks consistency across other avian species (e.g. Peterson et al. 2013) in order to confirm its suggested role as regulatory unit of migratory behaviour.

More recently, Delmore et al. (2015a) used next-generation sequencing data to estimate genomic differentiation between two subspecies groups of Catharus ustulatus that exhibit differences in their migratory behaviour (both timing and orientation). They characterised the entire genome to probe for enrichment of candidate genes for migration in areas of elevated differentiation across the genome. The resulting list of candidates included 25 genes that had been identified using a literature search (Ruegg et al. 2014). As predicted, genes from the list of candidates were enriched in areas of elevated differentiation, suggesting selection around these candidates could not only contribute to variation documented in their migratory behaviour, but also help maintain differences between the groups. Delmore et al. (2016) expanded on this work using hybrids and a genome-wide association study to identify one region in particular that is associated with variation in migratory orientation and harbours the CLOCK gene.

Despite several clearly species-specific differences in behavioural traits that make up the migratory phenotype, all migratory birds have to meet a general consensus schedule of key adaptations in order to cope with the challenge of migration. Specifically, migratory birds need a set of navigation and orientation mechanisms to migrate successfully, and must complete several key physiological traits such as hyperphagia and feather moulting in time prior to their migratory journey. Migratory individuals also share common morphological and anatomical adaptations, such as elongated and more pointed wings, lighter bone structure, maybe also brain volume (e.g. Lockwood et al. 1998; Fuchs et al. 2014). Given this common set of features among migratory birds and demonstrations of both their heritability and potentially shared genetic basis (e.g. the aforementioned involvement of CLOCK and ADCYAP1 across species), it has been postulated that there may be a shared set of genes for migration among birds (Berthold 1999).

Here, we test the hypothesis that sequence variation at one or a common set of suggested candidate genes for migration has been exploited to adapt to the requirements of migratory behaviour in different avian lineages using a cross-species comparative approach. Specifically, we focused on each of the 25 candidate genes for migration proposed by Ruegg et al. (2014) and compared patterns of genetic variability separately for each candidate gene between migratory and non-migratory species of birds. In addition to analyses per candidate gene, we also analysed concatenated sequence data from all candidate genes. This work benefited from the wealth of genomic data that have been accumulated in the last half-decade, with the assembly of several draft reference genomes for birds and large-scale initiative of the Avian phylogenomics project (http://avian.genomics.cn/en/index.html; more recently expanding to the B10K project). Our full dataset included 17 obligate migratory species, 32 sedentary or non-migratory species and 21 species with an additional intermediate movement phenotype (e.g. dispersive, partial migrant) (Table S1). We compare patterns of evolutionary divergence of each candidate gene using three different approaches: (1) comparing observed topologies for candidate genes to trees built using phylogenetic relationships with and without distinguishing migratory species from non-migratory species; (2) performing a gene-wide and branch-specific dN/dS analysis to identify if selective pressures on these candidate genes play a role related to migration; and (3) focussing on structural features that previous studies have shown to correlate with migratory behaviour (e.g. microsatellites at ADCYAP1 and CREB1 3′ UTR regions and poly-Q regions in NPAS2 and CLOCK), running linear models between these variants and two predictors of migratory behaviour. If genetic variation at candidate genes included in our study code for differences in migration, we predicted that gene trees based on candidate genes that play a role in shaping migratory behaviour would group migratory species together. Further, selective pressures on those candidate genes with a clear role in shaping the focal phenotype should be picked up in lineages with migratory species and linear models would show strong associations between genetic variation at structural features and predictors of migratory phenotype.

Materials and methods

Migratory phenotype characterisation

We classified each of the 70 species included in our study, according to their migratory phenotype. Our classification was based on a careful literature review of bird guides (Svensson et al. 1999), as well as BirdLife (http://birdlife.org) and Handbook of birds of the World (HBW) (http://www.hbw.com). We defined the following categories: clearly non-migratory (resident, sedentary) (0/R; n = 32), obligate migrant (2/M; n = 17). However, sometimes it is not easy to clearly define a species as either clearly non-migrant or obligate migrant, this is especially true when a migratory trait segregates within a population such as in partial migrants where only some individuals of the population migrate, consequently we added a third category (1; n = 21). This category includes partial migratory species (i.e. not all individuals of the population migrate), and species that exhibit other kind of migration-independent movement behaviour (e.g. dispersal, homing, foraging flights). For partial migrant species we used additional information of the individual used for generating the reference, such as date and geographical origin of sample collection, in order to clearly define migratory phenotype and grouped that individual/species accordingly whenever possible.

Genome sequences, extraction and alignment

We downloaded genome sequences and annotations for most of the species used in our study from the NCBI database (Supplementary Table S1). We further included genome sequences of five additional migratory species here: Siberian stonechat Saxicola maurus (Van Doren et al. 2017), Swainson’s thrush Catharus ustulatus (Delmore et al. 2015a), European blackcap Sylvia atricapilla (Delmore et al., in preparation), Willow warbler Phylloscopus trochilus (Lundberg et al. 2017, accepted), and Greenish warbler Phylloscopus trochiloides (Irwin et al. 2016).

Once we had sequences data and annotation for all of the species, we used the Bedtools (Quinlan and Hall 2010) getfasta module to extract genomic sequences for each of the 25 candidate genes for every species. Sequences for unpublished genomes or genomes without annotations (for details see Table S1) were generated using Blastn and chicken cDNAs from the Ensembl database as a reference. All genomic sequences were aligned with MAFFT (Katoh and Standley 2013) and manually edited in AliView. Coding (CDS) sequences were also obtained from a multiple alignment of the genomic sequences and Ensemble cDNA sequences (including untranslated regions, UTRs) for the flycatcher, chicken, and Zebra Finch. Only sequences covering 50% or more of the chicken genes were considered for further analysis.

Statistical analyses

Topological comparisons

Evolutionary trees were constructed for each candidate gene using whole genomic sequences and cDNA as reference, using a Neighbour Joining approach in MEGA v5.2 (Tamura et al. 2011). The reliability of the trees was evaluated performing a bootstrap analysis of 1000 replicates with the Kimura 2 Parameters model. To visualise the trees we used Figtree (http://tree.bio.ed.ac.uk/).

We compared the pattern of evolutionary divergence of these gene trees with three different hypothetical scenarios: divergence driven by phylogeny (‘phylogenetic topology’); divergence constrained by migration, i.e. different migratory phenotypes clustering in separate braches, while keeping the evolutionary relationship of the phylogenetic topology within each branch, (‘migratory phenotype topology’); and random divergence (‘random topology’). These comparisons were carried out for (a) the full dataset including all three migratory phenotypes; (b) a restricted dataset only contrasting exclusively obligatory migratory and completely non-migratory (resident) species, and (c) a clade-specific analysis exclusively focusing on the genus of Passeriformes, as this is the only monophyletic clade in our dataset with a sufficiently high number of species for both obligate migrants and non-migratory species, thus allowing for a more fine-tuned assessment on a narrower phylogenetic scale. The clade-specific subset allows us to test if the migratory phenotype might be controlled by a different clade-specific subset of genes. This comparative approach allows us to identify the presence or absence of general patterns, using genetic variation at candidate genes to distinguish between patterns related to phylogenetic relationships and migratory behaviour. The divergence driven by phylogeny (i.e. the gene trees matching the species tree, ‘phylogenetic topology’) was constructed using the total evidence nucleotide species tree (TENT) phylogeny, published by Jarvis et al. (2014). For species not included in the TENT phylogeny we used timetree (Hedges et al. 2015) divergence times to position these species in our phylogeny. The divergence constrained by migration scenario (‘migratory phenotype topology’) was constructed by clustering each phenotype (once exclusively focusing on migratory versus resident species for the restricted dataset; and also for the full dataset including other movement as a third phenotype category) in one separate branch while keeping the evolutionary relationships of the phylogenetic topology within each branch. Random divergence (‘random tree’) was generated shuffling branches randomly from the gene trees obtained, in order to avoid bias regarding the method of random trees generation by TOPD/fmts that only randomises taxa, but not branches for the statistical comparison. Restricting these analyses to exclusively Passerine species allowed us to analyse the effects of the evolutionary patterns of each candidate genes on a smaller scale. An example of the topologies is illustrated for the candidate gene PER3 in Fig. S1.

Comparisons of these three focal topologies were carried out in TOPD/fmtS (Puigbò et al. 2007) using three different approaches: nodal, splits and disagree from the program. In brief, the ‘nodal approach’ counts the number of nodes that separate two taxa in a given topology and calculates the root mean squared deviation (RSMD) between each pair of trees. For identical topologies RMSD results in a value of zero. To calculate the significance of the RMSD obtained, TOPD/fmts calculates the distance between two contrasted tree pairs and 100 random trees obtaining one standard deviation (SD) confidence interval (CI). Compared topologies are characterised as statistically similar, within noise or different, depending on their distance with respect to CI. Specifically distances below CI denote similar topologies; distances above CI indicate statistical difference (distances around CI are within noise). The ‘disagree method’ characterises how many taxa need to be removed from the compared topology in order to end up with the exact same topologies for both trees (assessed as count of taxa/total taxa). Consequently, a value of 0/total indicates identical topologies. The ‘splits method’ evaluates if there are common branches between both trees, the lower the distance the more branches the tree pair shares.

dN/dS analysis

Accounting for the fact that different species might have found different ways to alter similar phenotypes in the same gene (i.e. different changes in sequence), we also analysed synonymous and non-synonymous mutations of all candidate genes across species. In order to pick up on putative selective pressures on candidate genes for migration, a gene-wide dN/dS analysis was carried out on the Datamonkey server (Pond and Frost 2005). We used three different datasets for each candidate gene: one including all the species, one restricted to migratory, and another restricted to non-migratory species.

Gene-wide dN/dS ratios (w) were estimated by maximum likelihood (ML) methods using a different model for each gene. Each model was obtained from the CMS module of the server. Neighbour Joining (NJ) phylogenies obtained for each candidate gene were used as input to assess likelihood of the tree comparing the neutral null model M1 (w < 1) and the model M2 that allows w > 1. Positive selection was assessed if the likelihood shows a p < 0.05.

To evaluate if lineages with migratory species show a signature of selection, a branch-specific analysis of dN/dS was also carried on Datamonkey with the Branch-Site REL program (Pond et al. 2011). The dataset for each candidate included migratory and non-migratory species. Branches under episodic diversifying selection were identified with a Holm–Bonferroni corrected p ≤ 0.05.

Structural features and predictors of migration

We used a linear regression analysis to test for a correlation between the genotype of migratory species at each focal candidate gene and both breeding latitude and migratory distance. Models for both predictors were run separately. The genotype used for each gene was the microsatellite length (as number of bases) of the 3′UTR of ADCYAP1 and CREB1, or poly-Q (as number of predicted glutamine amino acids) on exon 20 of CLOCK and NPAS. For the CLOCK gene we included two separate polymorphic regions with variable poly-Q repeats in our analysis (both variable regions are located in the same exon). The significance of the fit was assessed with a simple linear regression, using a significance threshold p ≤ 0.05.

Within and across population variability in candidate gene sequence

Our comparative analyses focus exclusively on the sequence of one reference genome; inter-individual variation is not taken into account, mostly due to the limitation of available data to examine this level of variation. In order to make an attempt to see if variance within on candidate gene could be higher/lower in a specific migratory phenotype, we focused on CLOCK gene polymorphism, the only candidate gene with a sufficiently high number of individual sequence data available for several species (n = 10), including both migratory (n = 8) and resident (n = 2) species. Here we compare datasets of individually genotyped migratory species: flycatcher Ficedula hypoleuca (Saino et al. 2015; n = 226), willow warbler Phylloscopus trochilus ssp (unpublished data, n = 384), chiffchaff Phylloscopus collybita ssp (unpublished data, n = 56), nightingale Luscinia megarhynchos (Saino et al. 2015; n = 151), tree pipit Anthus trivialis (Saino et al. 2015; n = 144), barn swallow Hirundo rustica (Dor et al. 2011; n = 830), whinchat Saxicola rubetra (Saino et al. 2015, n = 374); and two non-migratory species: blue tit Cyanistes caerulea (Liedvogel et al. 2009; n = 950), great tit Parus major (Liedvogel and Sheldon 2010; n = 804). We compared averages and variances among different pairs of groups or species, employing a Welch t test and F test, respectively. We assume as statistically similar, distributions with a p > 0.001.

Results

Comparing phylogenetic, migratory phenotype and random topologies

We show constructed gene trees for a select number of candidate genes in Fig. 1, right column. In general, candidate gene trees do not separate migratory from non-migratory birds. Instead, they resemble the phylogenetic topologies expected based on Jarvis et al. (2014). Nonetheless, some groups of birds that comprised exclusively of migrants group together in most of the candidate gene trees. For example, Falconifomes, Accipithridae and some Passeriformes show a clustering pattern. Nevertheless, this is most likely due to their higher levels of similarity and common ancestry rather than relationships based on migratory phenotype (also see Fig. S2). Note that in addition to separate analyses for each of the 25 candidate genes, we also concatenated sequence data from all genes by species and constructed a combined candidate genes tree (Fig. 2). In this concatenated tree the clustering patterns persist in the aforementioned lineages.

Fig. 1
figure 1

Clustering patterns of migratory candidate genes. Right panel shows gene phylogeny for each candidate gene. Neighbour Joining (NJ) analyses shown for the four most widely used candidate genes for migration. Coloured dots indicate migratory (blue), non-migratory (red) and partial migrant/dispersive (yellow) taxa: node support is indicated by the size of nodes. Coloured clouds highlighted in red, yellow, green and purple highlight represent Passeriformes (red), Falconiformes (yellow), Accipitridae (green) and Anseriformes (purple), respectively. Left panel shows repeat lengths of focal genetic variants of each candidate gene, exemplarily illustrated for the most widely used candidate genes in the context of migration: ADCYAP1, CLOCK, NPAS and CREB1. Genotype is plotted in relation to migratory distance (open circles, dashed lines) and breeding latitude (filled circles, continuous lines). Dashed and continuous lines indicate the trend for linear regression of repeat variation at the candidate locus versus migratory distance and breeding, respectively. Variation at the CLOCK genes is shown for both variable length repeat regions (also see Fig. S3). R squared values are state for fitted linear models

Full size image
Fig. 2
figure 2

Phylogenetic tree of concatenated candidate gene sequences. Neighbour Joining topology of genomic sequences for all 25 candidate genes. Colour scheme and node support as in Fig. 1

Full size image

Results from nodal, splits and disagree methods for statistically distinguishing between phylogenetic, migratory phenotype and random trees can be found in Table 1 (full dataset). Our results clearly show that the phylogenetic topology tree provides the best fit to most gene trees. The lack of support for the migratory phenotype topology shows no evidence for a monophyletic origin for a migratory phenotype across avian taxa. This pattern did not change when restricting our analysis to only the extreme phenotypes, i.e. exclusively obligate migrants and clear non-migrants. We further zoomed into just one clade (specifically Passerines). Although overall, we see general patterns either more similar to the speciation phylogeny or not showing differences from a comparison with random topologies (see nodal approach in Table 1), the monophyletic clade analysis on Passeriformes indicates that respective gene topologies for some candidates, specifically HRSP12 and HSPA5, are consistently more similar to the migration topology than the speciation topology (Fig. 3). This clustering pattern is also evident in the tree topologies (see Fig. S2) where most of the migrant species within Passerine tend to cluster in one branch. The other approaches did not show evidence for a general trend towards one or a set of candidate genes being recurrently involved in distinguishing migrants from non-migrants.

Table 1 Topology comparison among candidate gene trees and target trees
Full size table
Fig. 3
figure 3

Topologies on a constrained phylogenetic scale. Clade-specific analysis exclusively focusing on all species within the genus Passerines. The candidate gene trees obtained from the Neighbour Joining analysis for HRSP12 and HSPA5, two candidate genes that follow clustering patterns consistently more similar to the migration topology than the speciation topology. Colouring scheme and node support as in Fig. 1

Full size image

Recall that we also concatenated sequence data from all genes by species and constructed a complete gene tree. This concatenated phylogeny shows that most of the taxonomic groups clustering migratory species together are statistically well supported and show larger branch lengths. This could suggest an evolutionary process of acceleration in these lineages for some of the candidate genes.

Selection in candidate genes is not related to migration

To identify selective forces in lineages with migratory species, we performed a gene-wide and branch-specific dN/dS analysis. The gene-wide analysis for selection in candidate genes across all bird species indicates that CLOCK, DRD4, NEK2, HSPA5 and CSNKE1, have been under positive selection at p < 0.05 (Table 2). Nonetheless, and important within the focus of our analysis, this signature does not seem to be linked to the migratory phenotype in any of our comparisons, as independent datasets of separately analysed migratory, partial migratory and non-migratory species generally do not show any signature of selection. The only exception here was the CLOCK gene that showed significant values for selection in partial migrants. Thus, despite the general overall pattern of selection across all avian clade, none of these candidate genes showed lineages under selective pressure.

Table 2 Gene-wide dN/dS analysis of migratory gene candidates
Full size table

Association between structural features and predictors of migration

Alignments of the genomic sequences (Fig. S3) of any of the 25 gene candidates did not separate migratory and non-migratory species. ADCYAP1 has variable sequence lengths of a microsatellite with an AG repeat at the 3′ UTR region ranging from 26 to 56 bp for migratory species, and 10–54 bp for non-migratory species. CREB1 also has a microsatellite with a TG/CG repeat motive at the 3′ UTR region, ranging from 12 to 26 bp in migratory and 14–42 bp in non-migratory species, respectively (the longest length variant was exclusively found in tits). NPAS has one variable region of poly-Q repeats, in migratory species; this has length variation between 6 and 10 amino acids (aa), in non-migratory species the repeat length varied between 5 and 12 aa repeats. The CLOCK gene has 2 poly-Q variable regions, both located on exon 20. The first region (R1) shows length variation between 6 and 13 poly-Q repeats in migrants and 4 to 14 repeats in non-migratory species. Length polymorphism in variable region two (R2) varied between 4 and 9 poly-Q repeats in both phenotypes. Lengths variation in neither of the focal candidates showed significant differences between migratory and non-migratory birds.

The intra-specific analysis focusing on individual variability of CLOCK gene polymorphism shows a tendency for higher variability in migratory species (Fig. S4); however, given that only two non-migratory species were included (namely great tits and blue tits—great tits being mostly monomorphic, blue tits highly variable), this does not allow us to draw any statistically supported conclusions (Table S3). Nonetheless, this individual based analysis of variability suggests that most common allele lengths (i.e. number of poly-Q repeats) varied considerably between species, as did the overall degree of variability within species (Fig. S4, Table S3). For example, migratory species like H. rustica and F. hypoleuca, show constrained levels of variability, being mostly monomorphic, while C. caeruleus, a non-migrant species, shows a degree of variance comparable to those migrant species. Statistical differences were found between the group of migrant and non-migrant species. Nonetheless when the comparisons were performed between individual species, migrant species and non-migrant species showed statistical differences, but also comparisons between a migrant and a non-migrant species did not show statistical differences (see Table S3).

We used linear models to quantify the correlation of structural repeats (poly-Q in CLOCK and NPAS2; and allele lengths in ADCYAP1 and CREB1) and both breeding latitude and migratory distance (Fig. 1, left column). The fit of the data to a linear regression model did not show significant associations between genotype and migratory phenotype (Table 3). The best fit, albeit not significant, that was following a similar pattern as earlier work showing a correlation between candidate gene variation and migratory phenotype results from our comparison between ADCYAP1 genotype variation in relation to migratory distance (r 2 = 0.2120, p = 0.01132). Although not significant, here the regression trend is in line with earlier studies, and might be taken as support for earlier findings reporting genotype variation at this locus with relevance to the migratory phenotype on the within-species level.

Table 3 Predictors of linear models
Full size table

Discussion

We used publicly available archived genomic data from non-migratory and migratory bird species to test for the presence or absence of a general clustering pattern in candidate genes for migration. In a cross-species comparative approach, we characterised sequence features in 25 candidate genes in a dataset including birds across all orders of Aves. Our study thus illustrates the potential of public genomic data to test for general patterns of migratory gene candidates in a cross-species comparative context. Despite this powerful dataset we were not able to identify genetic variation that allowed us to distinguish migratory from non-migratory birds based on sequence difference in any of the candidate genes included here.

Most patterns we recover based on candidate gene sequence variation were driven by species phylogeny, and did not show indications separating migratory and non-migratory birds. This could suggest that candidate genes do not play a general role in controlling migration across the full avian clade, or that patterns cannot be evident in a wide phylogenetic analysis. In addition to analysing each candidate gene individually, we also analysed concatenated sequence data from all candidate genes, which also did not allow clear separation of migratory and non-migratory birds. The latter finding might not be too surprising, as this scenario would assume divergence of candidate genes among lineages in parallel to species divergence, a rather unlikely scenario.

In order to address the fact that we cannot exclude that the relevance and contribution of certain traits may differ between avian lineages in a way that different traits may be controlled by different mutations in the same genes, and/or varying patterns of epistatic interactions with other genes, we also screened for putative features of selection within these candidate genes. However, our dN/dS approach did not allow identifying patterns of selection at either the gene-wide or branch-sites level for the migratory lineages. This result does not provide support for lineage-specific selection on the gene candidates.

Indicative patterns might get lost if not approached at the appropriate phylogenetic scale. Our more focused analyses on a restricted data assembly exclusive to the Passeriformes allow aims to address this issue at least partly allowing us to potentially recover general trends towards some or all candidate genes being recurrently involved in the control of migration. Analyses limited to this family show improved resolution and come close to the margin of randomness in nodal comparisons. However, even on this more focal scale only a few of these comparisons show better values than the phylogeny species tree, suggesting that even on a more constrained phylogenetic scale, the genes keep following the trends of phylogeny or speciation rather than clustering according to migratory phenotype.

Most studies of candidate genes look for correlates of length polymorphism in a limited number of species. Our approach aims to extend this framework to assessing sequence variability in candidate genes to a broader evolutionary scale. However, our data do not uncover statistical differences between species with different phenotypes, neither do they recover correlations with specific migratory features, specifically distance and breeding latitude. ADCYAP1 did show a tendency for positive, though non-significant, correlation with migratory distance in our dataset, supporting earlier findings in species-specific analyses. As the length polymorphism in this gene is located in a 3′ UTR this could suggest that a specific pattern of gene expression is playing a role for this gene, and this feature might be relevant in long-distance migrants.

Our focus is put on general sequence differences in candidate genes for migration across species, but we also note that there is a lot of variation in the level of individual and inter-population diversity across avian species in many of the candidate genes included here, such as CLOCK poly-Q. Within-species variability at polymorphic loci often is high in some, but certainly not all species, irrespective of migratory phenotype. Previous work has suggested that variability in lengths polymorphism of CLOCK related to the timing of migration is enhanced in migratory birds (Saino et al. 2015). Our analyses testing for a role of inter-individual variation in CLOCK lengths polymorphism between migratory and non-migratory species falls short in strengthening this suggestion. We detect similarly high levels of polymorphisms as well as scenarios with almost monomorphic genotypes for various species, irrespective of migratory phenotype. This could mean that the observed variability in CLOCK is more tightly linked to other parameters that shape the variability of the migratory phenotype, such as breeding latitudes, distances and timing regimes, rather than specifically controlling timing traits in a migration context.

But even when leaving inter-individual variability aside and focusing on just one reference sequence, we expect general patterns that allow separating migrants from non-migratory species of suggested candidate genes for migration to show up if they were in fact under parallel selection. Our results do not support this hypothesis and thus highlight limitations and question the usefulness of a candidate gene approach in the context of understanding migratory behaviour. Consequently, we call for caution and highlight the limitation of candidate gene approaches in macro-evolutionary contexts, as in most cases functionality cannot be easily inferred. Also, most candidate gene studies do not allow for a clear distinction from genetic drift, come with uncertainty about the amount of loci involved in the trait and do not allow controlling for possible effect of linkage disequilibrium. Our study further points out how this approach can erroneously simplify a highly complex phenotype that most probably is a multilocus trait.

One further point we feel is especially important to keep in mind when using a candidate gene approach in the context of migration: many candidate genes have been identified in genetic model organisms—none of that migrate or show any correlate to a migratory phenotype, and functional validation or relevance is mostly lacking. In many of these studies, the subset of candidate genes for migration were selected based on their assumed effect on anticipated candidate traits that feed into the complex phenomenon of migration. Importantly, one possibility explaining the lack of pattern might be the limitation of the starter set of ad hoc hypothesised candidate genes we investigate and may not be completely unexpected. Once more genomic data for migratory and non-migratory species will become available, the genomic toolbox allowing to investigate the genetic basis of migratory traits will grow and allow for an increasingly more informed list of de novo identified candidates to be tested in migratory birds.