On the importance of time scales when studying adaptive evolution

Abstract Long‐term field studies coupled with quantitative genomics offer a powerful means to understand the genetic bases underlying quantitative traits and their evolutionary changes. However, analyzing and interpreting the time scales at which adaptive evolution occurs is challenging. First, while evolution is predictable in the short term, with strikingly rapid phenotypic changes in data series, it remains unpredictable in the long term. Second, while the temporal dynamics of some loci with large effects on phenotypic variation and fitness have been characterized, this task can be complicated in cases of highly polygenic trait architecture implicating numerous small effect size loci, or when statistical tests are sensitive to the heterogeneity of some key characteristics of the genome, like variation in recombination rate along the chromosomes. After introducing these aforementioned challenges, we discuss a recent investigation of the genomic architecture and spatio‐temporal variation in great tit bill length, which was related to the recent use of bird feeders. We discuss how this case study illustrates the importance of considering different temporal scales and evolutionary mechanisms both while analyzing trait temporal trends and when searching for and interpreting the signals of putative genomic footprints of selection. More generally this commentary discusses interesting challenges for unraveling the time scale at which adaptive traits evolve and their genomic bases.


Impact summary
An important goal in evolutionary biology is to understand how individual traits evolve, leading to fascinating variations in time and space. Long-term field studies have been crucial in trying to understand the timing, extent, and ecological determinants of such trait variation in wild populations. In this context, recent genomic tools can be used to look for the genetic bases underlying such trait variation and can provide clues on the nature and timing of their evolution. However, the analysis and the interpretation of the time scales at which evolution occurs remain challenging. First, analyzing long-term data series can be tricky; short-term changes are highly predictable whereas long-term evolution is much less predictable. A second difficult task is to study the architecture of complex quantitative traits and to decipher the timing and roles of the several genomic mechanisms involved in their evolution. This commentary introduces these challenges and discusses a recent investigation of the nature and timing of ecological and genomic factors responsible for variation in great tit bill length. Overall, we raise cautionary warnings regarding several conceptual and technical features and limitations when coupling analyses of long-term and genomic data to study trait evolution in wild populations.
strength and directionality of natural selection underlying such variation. Some of these long-term examinations of key phenotypic traits detected strikingly fast phenotypic change, driven by rapid ecological changes (Grant and Grant 2006). In contrast, many long-term studies failed to reveal micro-evolutionary change and response to selection (Merila et al. 2001;Grant and Grant 2002). While examining longitudinal data looking for both long-term trends as well as short-term fluctuations has the potential to shed light on evolutionary trajectories in natural populations, our ability to understand and more notably predict evolution remains limited outside the laboratory. Quantitative genetic models were initially developed with the aim of predicting evolutionary change, based on estimates of selection and additive genetic variation (Falconer 1960). Their predictive power worked efficiently for the genetic improvement of complex traits in many animal and plant breeding programs. Yet when spatio-temporal ecological heterogeneity is involved, evolution in the wild remains largely unpredictable (Pemberton 2010;Pujol et al. 2018).
Coupling such long-term studies with genomic tools is a powerful way to improve our understanding of the genetic bases underlying evolutionary changes in response to environmental variation. Rapid and recent monogenic adaptations based on de-novo mutations are often used as examples, for instance the rise of the melanic morph of the peppered moth Biston betularia following the industrial revolution (van't Hof et al. 2016). Similarly, there are famous examples of rapid parallel monogenic or oligogenic adaptation based on long lasting standing genetic variation segregating in heterogeneous environments, for example coloration in deer mice Peromyscus and armor plates in stickleback Gasterosteus aculeatus (Barret and Hoekstra 2011;Nelson and Cresko 2018). However, genomic analyses are facing several challenges when it comes to making inferences about polygenic adaptation and quantitative trait evolution. First, loci effect sizes are often small, requiring thousands, if not millions, of both single nucleotide polymorphisms (SNPs) and individuals for genome wide association studies (GWAS) to reveal significant effects. Such an investigation requires technical commitments (Wellenreuther and Hansson 2016;Gienapp et al. 2017a), but also a conceptual shift toward suppressing our desire to discover large effect alleles that are, in theory, rarely responsible for quantitative variation (Rockman 2011). Second, genome characteristics and especially variation in recombination rate along the genome can cause variation in the extent of background selection (defined as the loss of genetic diversity at a neutral locus due to negative selection against linked deleterious alleles) (Charlesworth et al. 1993;Charlesworth 2013;Nordborg et al. 1996). This variation in recombination can confound or bias detections of positive selection and hence our comprehension of the timing and nature of evolutionary trajectories based on genomic data (Roesti et al. 2012;Burri et al. 2015;Berner and Roesti 2017;Burri 2017;Comeron 2017;Delmore et al. 2018). Nevertheless, several studies have begun to decipher the polygenic mechanisms of rapid evolution of quantitative traits.
Avian bill morphology has played a prominent role in empirical studies of evolution and natural selection (Lack 1947;Grant 1999), perhaps because the size and shape of bills show large variations across and within bird species and are shaped by strong selective forces since they directly determine foraging efficiency on various food sources. For instance, the emblematic study of Darwin's finches on the Galapagos island of Daphne Mayor aimed at capturing evolutionary changes in bill size (Boag and Grant 1981;Grant and Grant 1993). From 1977 to 1978, bill size increased markedly after a severe drought in 1977. This analysis clearly demonstrated that extreme climatic events such as an El niño event are strong drivers of bill size evolution in this species. After 30 years of perspective, however, Grant and Grant (2002) concluded that while evolution of bill length was predictable as a rapid response to strong selection, it remained unpredictable on a slightly longer microevolutionary scale. Although the question of predictability of evolution across time scales remains challenging, even in the genomic era (Nosil et al. 2018), genomic tools did provide insights on the evolution of bill size in the context of the rapid diversification of Darwin's finches. A handful of genes were found to be significantly associated with bill size and shape in the medium ground finch Geospiza fortis (Lamichhaney et al. 2016), among which a major locus has been shown to influence bill dimensions in the Darwin finches' entire radiation (Lamichhaney et al. 2015). These genomic analyses hence cracked the genomic architecture of this trait variation at both small and large microevolutionary scale, with the predominant control of a few large effect loci.
In a recent study, Bosse et al. (2017) investigated the genetic architecture of bill length in the Great tit Parus major, using a tremendous amount of data from long-term research programs in Wytham woods in the UK and in the Netherlands (NL). Applying modern analyses of both population and quantitative genomics using 500 k SNP, the authors provide insight into the signatures of divergent selection in the studied populations and the genomic architecture of variation in bill length. In line with the quantitative nature of variation in bill length and with quantitative genetics theory (Lynch and Walsh 1998), the authors showed that the genetic architecture of bill length was highly polygenic. Specifically, the authors showed, using a mixture analysis fitting all the SNPs simultaneously, that 3009 SNPs explained collectively 31% of bill length phenotypic variation. None of these SNPs reached genome wide significance in the GWAS with bill length, revealing small effect sizes of individual variants. In accordance with recent quantitative genomic findings notably for other traits in great tits (e.g., Robinson et al. 2013), the proportion of variance in bill length explained by each chromosome amazingly scaled with its size, demonstrating that the many SNPs additively explaining bill length are distributed throughout the genome. The polygenic analysis also predicted the difference in bill length observed between the UK and the NL, further illustrating the polygenic nature of this trait variation. This evidence for a polygenic control of bill length with no large effect SNPs is hence very different from the previous example on Darwin finches' bill with large effect loci. Nevertheless, Bosse et al. (2017) also showed that variation at a single gene, col4A5, was associated with bill length (although not reaching genome wide significance). Bosse et al. (2017) then discussed whether an extended use of feeders, that are more abundant in the UK, might have driven the evolution of larger bills in the UK compared to the NL. Among the arguments pointing to feeders as drivers of longer bill lengths in UK great tits, Bosse et al. (2017) reported that col4A5 was associated with bill length in the UK but not in the NL, was highly differentiated between the UK and the NL, and was associated with greater reproductive success and higher activity at feeding sites in the UK. In addition, they reported that bills were longer in the UK compared to mainland Europe ( We argue in this comment that the speculation on the role of feeders on bill length evolution is not well supported by these arguments. Based on both phenotypic and genomic data, we propose instead that differences in bill length between the UK and the NL might have been evolving on a longer time scale than the contemporary use of feeders. First, an examination of bill length monitoring shows puzzling temporal and geographic patterns that may not incriminate the use of feeders. Second, the genomic patterns found at the region containing col4A5 might be compatible with the combined effect of variation in recombination rate and background selection, hence questioning the putative recent role of feeders as agents of positive selection. The bill length trend inferred from the long-term monitoring is highly dependent on the time scale considered. Inspired by the readings of articles such as the Grant and Grant (2006) study relating the effects of rare and extreme events on the evolution of phenotypic traits, we carefully inspected the evolution of bill length during the studied period, looking for particularly rapid changes. We reanalyzed the data using a breakpoint computations method (coin R-package) aiming at localizing such striking change. The best cutpoint was found between 1986 and 1987 (maxT = 7.41, P-value < 2.2 × 10 -16 ), which corresponds to a conspicuous change in bill length. Measures from the five years preceding 1987 differed from subsequent years (t-test: t = −7.28, df = 674.53, P-value = 9.1 × 10 -13 ). Although a linear regression through the entire period, from 1982 to 2007, as implemented by Bosse et al., yielded a positive slope, removing the five measures prior to 1987, bill length significantly decreased from 1987 to 2007 (Fig. 1A 8 1982 1986 1990 1994 1998 2002 2006 Bill length (mm)  6 1982 1986 1990 1994 1998 2002 2006 Tarsus length (mm)  59 1982 1986 1990 1994 1998 2002 2006 Year Bill length / Tarsus length tarsus length into account, or not, F = 7.13, P-value = 8.2 × 10 -4 ; F = 13.2, P-value: 2.9 × 10 -4 , respectively). A LOESS model confirms the decrease in bill length during the second part of the record. More generally, slopes of linear models linking bill length to time, using all of the possible combinations of 10-25 consecutive years (simply removing 1-16 years at the beginning, or the end, or both, of the data-series), were often positive when including one or more data points collected between 1982 and 1987, while negative slopes were often observed when excluding these years (Fig. 2), further illustrating that both time periods yield opposite patterns. Furthermore, tarsus length also decreased .9 × 10 -7 , respectively). Correspondingly, the bill length/tarsus length ratio did not change significantly from 1987 to 2007 ( Fig. 1 C, F = 1.51, P-value = 0.28), indicating no change in bill length, when scaled to tarsus length, over this period. The origin of this cutpoint is important to investigate, since it could constitute evidence for a response, either genetic or plastic, of both bill and tarsus lengths to a rapid abiotic or biotic change or to a methodological change. It seems at present difficult to decipher whether bill length or tarsus length or both, if any, were targeted by selection. Indeed, both traits are commonly phenotypically and genetically positively correlated in passerines (Teplitsky et al. 2014;Poissant et al. 2016), hence should often evolve together. Overall, these results very likely rule out the possibility of a contemporary  positive effect of feeders on bill length. One could nevertheless argue that bill length could have increased at a longer time scale than the 1982-2007 period, but still recent enough to incriminate the use of feeders. Figure 3 however shows the absence of a clear pattern of bill length increase in museum specimens collected in the UK. Although more data are needed to confirm the pattern, it suggests that larger bill length in the UK seems to have evolved over a longer time period than the one during which feeders have been used. Moreover, bill length across Europe does not display a clear dichotomy between the UK and mainland Europe but rather smooth spatial variations

. Variation in bill length across European countries for both female (in blue) and male (in green) great tit specimens from museums. Letters A and B illustrate significant differences between sites and "ns" refers to nonsignificant differences (Tukey HSD).
( Fig. 4), with an ANOVA showing a significant effect of country (P = 3.75 × 10 -7 ) but not of sex (P = 0.345), with UK birds having longer bills compared to three other countries (France, Italy, and the NL). This spatial variation also challenges the suggestion by Bosse et al. (2017) of a recent increase in bill length in British great tits caused by feeders.
We then question whether the evolution of the genomic region containing the gene col4A5, which was the corner stone linking bill length, activity at feeders, reproductive success, and divergent selection, could have been influenced by recombination rate variation and background selection, rather than recent positive selection due to feeders. While the decay of LD is very fast in the great tit genome (marginal after 2 kb, as shown in Fig. S1 in Bosse et al. 2017, and in Fig. 5A here (see supplementary material 1 for method details)), col4A5 was found in a large (>1 Mb) genomic region with high long distance (20-200 Kb) LD in the UK (Fig. 5B). The F ST between the UK and the NL was high along this region (Fig. 5D here; Fig. S3A in Bosse et al. 2017). The eigenGWAS in this region was significant using both populations simultaneously for the entire set of SNPs (upper panel Fig. S2 in Bosse et al. 2017) and for chromosome 4A only (Fig. 5E). As argued by Bosse et al. (2017), such results are compatible with recent strong positive selection in the UK over this large stretch of DNA on chromosome 4A. However, this large stretch of high LD around col4A5 was not only found in the UK but also in the NL (Fig. 5C). The eigenGWAS in this region was also significant for both populations separately ( Fig. 5F and G). Furthermore, this large region was previously identified by Laine et al. (2016) as showing a signature of selective sweeps and reduced nucleotide diversity at the scale of the entire species distribution rather than only in the UK (Sixth chromosome in Fig. 2 in Laine et al. 2016). This same region was also identified as showing elevated differentiation in several lineages and populations of Ficedula flycatchers (Ellegren et al. 2012;Fig. 1C in Burri et al. 2015). Therefore, given the existence of this large genomic region with reduced variation, increased LD, and increased divergence at several spatial scales in great tits, and in flycatchers, it seems unlikely that the mechanism shaping this region has been acting only in the UK, recently, and implying mainly positive section. Burri et al. (2015) determined that the high differentiation and high LD at this region, shared across flycatcher lineages, was due to the effect of linked selection combined with low recombination and issued a crucial warning: "scans are likely to identify recombination-mediated elevations of differentiation not necessarily attributable to selective sweeps." Accordingly, we propose that low recombination (potentially reflecting pericentromeric regions) and background selection in the region containing col4A5 in great tits could have resulted in locally reduced genetic diversity, increased differentiation and increased LD that could have altogether mimicked signatures of recent positive selection.
To illustrate the general challenge raised here in differentiating recombination rate variation combined with background selection from positive selection, we present a very simple simulation of a 20 Mb chromosome (hence comparable to the length of chromosome 4A in great tits) containing a region of 1.6 Mb with greatly reduced recombination (comparable to the length of the col4A5 haplotypes), with random occurrence of neutral and deleterious mutations but no beneficial ones (see supplementary material 2 for methods details). We show that for an average LD decay comparable to what was found in great tits (Fig. 6A), long distance LD (Fig. 6B), eigenGWAS (Fig. 6C), integrated extended haplotype homozygosity (Fig. 6D) and nucleotide diversity (Fig. 6E) all show striking deviations in the recombination coldspot compared to the rest of the chromosome. Therefore, this simple simulation illustrates that such a pattern as that observed by Bosse et al. (2017) is compatible with the combined action of large-scale variation in recombination and background selection.
Neglecting the effects of variation in recombination and LD along the genome might have not only resulted in false-positive footprints of selection in coldspots of recombination but also in false negatives in regions with a higher recombination rate (Berner and Roesti 2017). The large sliding window used in Bosse et al. (sliding over 200 kb while the average LD is highly reduced after 2 Kb) potentially worsened this problem by capturing principally elevated differentiation in large regions with increased LD possibly resulting from extended lower recombination due to structural variations. It may have also diluted narrow (i.e., narrower than the sliding window) peaks of elevated differentiation in more common regions with lower LD outside of recombination coldspots. Unfortunately, this probably relatively common caveat holds regardless of the nature of selection (i.e., positive or purifying) and is a purely mathematical problem of mismatch between the sliding window length and the extent of LD variation along the genome. It typically occurs when the unit of the sliding window is the base pair instead of the centimorgan. In fact, such large sliding windows relative to the LD decay should be used as neutral local baselines to ascertain the effects of variation in LD and recombination along the genome, with local outliers detected when comparing local residuals to such baselines (Roesti et al. 2012;Burri 2017).
All these considerations suggest that, although the pattern of increased divergence at col4A5 is apparently compatible with strong positive selection, as suggested by Bosse et al (2017), the combined role of background selection, strong recombination rate variation, and invariably large averaging window while long distance LD is variable, should be more comprehensively tested. Correspondingly, accounting for these factors may also uncover more variants under stronger positive selection located outside of the few low recombination regions. In this context, the causality of the associations between col4A5 and bill length, activity at feeders and reproductive success in the UK, requires further clarifications.
This case study offers exciting avenues of research to unravel the determinants of both recent and long-term as well as spatial variations in quantitative traits in the Great tit but also in other emblematic species displaying quantitative trait variations. Spatial trait variation could be unraveled in multitrait GWAS using polygenic frameworks, inspired by a pioneering study on great tit breeding phenology (Gienapp et al. 2017b). Additionally, a formal selection analysis relating bill length to overwinter survival when birds are most likely to benefit from food provided by people, is required to elucidate the nature of the evolutionary forces behind bill length variation. Finally, including the effect of variation in recombination rate, background selection, and LD along the genome to draw local neutral envelopes of genomic differentiation and modulate the local width of sliding windows, will likely help identify further candidate loci (Roesti et al. 2012(Roesti et al. , 2013Burri 2017;Comeron 2017;Delmore et al. 2018). Genomic analyses performed in other populations across Europe could also help determine the timing and the nature of selection, by taking into account the temporal dynamics of the differentiation landscape (Burri 2017). We wish to conclude by emphasizing the importance of integrating, or at least recognizing, the widespread and sometimes strong variation in recombination rate along genomes, which can in some circumstances distort our understanding of evolutionary processes based on genomic investigations (Roesti et al. 2012;Burri et al. 2015;Berner and Roesti 2017;Burri 2017;Comeron 2017;Delmore et al. 2018).

AUTHOR CONTRIBUTIONS
Conceived the study: CP & AC; Performed the analyses: CP; Wrote the manuscript: CP & AC.

ACKNOWLEDGMENTS
We thank several colleagues from the Institute of Evolutionary Sciences of Montpellier (ISEM) and from the Center for Functional and Evolutionary Ecology (CEFE), especially Paul Jay, Tim Janicke, Nicolas Bierne, Pierre-Alexandre Gagnaire, François Bonhomme, and Patrice David, for discussions regarding the impact of recombination rate variation on the landscape of diversity and differentiation. We thank Caroline Rose for comments on the English syntax. We thank the associate editor, three anonymous reviewers, and JC Senar for constructive comments on previous versions of the manuscript. Both authors are presently funded by a European Research Council Starting grant-Selection in Heterogeneous Environment (ERC-2013-StG-337365-SHE). The authors declare no conflict of interest.

DATA ARCHIVING
Location and description of the genomic and phenotypic data are listed in Bosse et al. 2017.