Environmental change, as, e.g., global warming, will inevitably lead to novel selection pressures, especially in phenological traits, as, for example, timing of breeding, migration or hibernation (Gienapp et al., 2014), and in the long run only adaptive microevolution (‘evolutionary rescue’) will enable population persistence (Bürger and Lynch, 1995; Visser, 2008; Gienapp et al., 2013; Gonzalez et al., 2013). Understanding the genetics of traits that are affected by environmental change is important, as this will allow us to assess the adaptive potential of these traits. From quantitative genetic studies we know that many traits, including ecologically relevant phenological, life-history or behavioral traits, are heritable (Mousseau and Roff, 1987), and hence that evolutionary adaptation to climate change would be possible. Heritability is, however, not necessarily constant across environments (Pigliucci, 2001, 2005; Nussey et al., 2007; Husby et al., 2011), which complicates predictions about evolutionary adaptation. This is especially true when environmental conditions will shift beyond the currently observed range due to climate change. Consequently, in the end a better understanding of the genomics of the traits under selection will help us to predict whether they will be able to adapt fast enough to novel environmental conditions.

Due to the advancement in molecular genetics, it is now possible to map the loci that underlie traits of interest in natural populations of ‘ecological model species’ and this has moved the research field of ecological genetics into the genomics era (Slate et al., 2010). So far, however, the number of published mapping studies in natural populations is limited. Over the last decade, over 20 published studies aimed to identify loci underlying almost 40 traits in natural populations (Supplementary Table S1). A number of studies successfully identified some of the loci underlying, for example, wing length in birds, age at maturity in Salmon, or birth weight and recombination rate in mammals (Slate et al., 2002; Tarka et al., 2010; Barson et al., 2015; Johnston et al., 2016). In these studies the identified loci explained moderate to substantial amounts of the genetic variation but in other studies that could identify loci the amount of variation explained by them was (very) small (e.g., Husby et al., 2015; Wenzel et al., 2015). In general, in only about one third of the analyzed traits, one or more loci or genomic regions associated with the trait could be identified. The total number of studies aiming at mapping loci in wild populations that could not identify any locus may even be larger because these studies may simply not have been published due to their non-significant results. There are a number of possible reasons behind this limited success of gene-mapping in wild populations. For example, in some studies the number of individuals was only in the low hundreds. Heritabilities of life-history traits also tend to be low to moderate (Mousseau and Roff, 1987), which would make detecting loci more difficult. Another potential reason for the limited success of these gene-mapping studies in wild populations could be that phenotypic plasticity complicates the detection of relationships between genotype and phenotypes. Many life-history and behavioral traits but also morphological traits are phenotypically plastic (Pigliucci, 2005).

In phenotypically plastic traits the same genotype can produce different phenotypes, and these different phenotypes can hence not be the product of genomic sequence variation alone but must be caused by environmentally induced differential gene expression or post-translational mechanisms (Snell-Rood et al., 2010; Weake and Workman, 2010; Huntzinger and Izaurralde, 2011). This means that the effects of certain genes will be reduced, entirely shut down, or altered under specific environmental conditions. Consider, for example, three genes, A, B and C, affecting a plastic trait: In environment X the expression of gene A may be up-regulated and the expression of gene B down-regulated, while in environment Y the expression of gene A may be down-regulated and the expression of gene B up-regulated. In contrast, the expression of gene C may be neither up- nor downregulated in these environments. A mapping study not taking phenotypic plasticity into account will only identify gene C as associated with the plastic trait because sequence variation in genes A and B is only associated with the trait in part of the data. If post-translational mechanisms are affected, the effect of a gene could not only disappear, due to environmentally dependent downregulation, but its phenotypic effects could even be reversed in one environment compared to the other. Only when modeling an interaction between marker effects and environment the effects of genes A and B will be detected. Such a variable genomic architecture of plastic traits may make detecting loci underlying such traits more difficult if phenotypes are measured in different environment, which is likely to be the case in natural populations, and phenotypic plasticity is not taken into account. Genotype-by-environment interactions have been tested in laboratory or medical studies (Mackay et al., 2009; Thomas, 2010) but not in studies of natural populations where their effects should, however, be even more important.

Breeding time, i.e., the timing of egg-laying, is a highly plastic trait in birds. Depending on spring temperature breeding time can vary by weeks from 1 year to the next (Gienapp et al., 2005, 2010; Charmantier et al., 2008). Besides temperature, photoperiod is known to affect breeding time in birds (Lambrechts et al., 1997; Lambrechts and Perret, 2000) and there is also evidence that the roles of temperature and photoperiod differ in early and late spring (Gienapp et al., 2005, 2010). As pointed out above, phenological traits are expected to come under selection from climate change and increasing spring temperatures have led to mismatches between breeding time and food supply in birds and other species (Visser et al., 1998; Visser and Both, 2005). Understanding the genomics of avian breeding time is hence interesting in the context of climate change, which makes it a very suitable candidate to test whether a variable genomic architecture can obscure genotype-phenotype links. We here use a unique data set consisting of wild 2045 female great tits with recorded egg-laying dates that also have been genotyped on a recently developed 650k SNP chip to test whether the association between SNPs and egg-laying date interacts with spring temperature.

The heritability of egg-laying date varied among temperature classes. In line with earlier results (Husby et al., 2011) it was lowest under cold spring temperatures (0.14 ± 0.05, estimate and SE) and increased with temperature to 0.38 ± 0.06 and 0.41 ± 0.06 in the medium and warm temperature classes, respectively. The additive genetic variance was statistically significant for all three temperature classes (LRT with 1 df: χ² = 9.39, p = 0.002; χ² = 41.9, p < 0.001; χ² = 57.5, p < 0.001, Table 1).

When all years were analyzed together, no SNP was associated with year- and area-centered egg-laying date at the genome-wide significance level (Figure 1). When we tested whether SNP effects differed between cold, intermediate and warm springs, without fitting heterogeneous residuals, two SNPs reached genome-wide significance and one SNP was close to genome-wide significance (Supplementary Table S2). In this model p-values were, however, substantially inflated (Supplementary Figure S3A). Fitting heterogeneous residuals in the SNP^∗temperature interaction model removed this inflation (Supplementary Figure S3B) but also meant that no SNP reached genome-wide significance anymore (Figure 2 and Table 2). The two SNPs that reached genome-wide significance and the SNP that almost reached genome-wide significance were, however, still among the 10 most significant SNPs when heterogeneous residuals were fitted (Supplementary Table S2). Three SNPs that were among the 10 most significant SNPs were located within the thyroglobulin (TG) gene, the heparan sulfate-glucosamine 3-sulfotransferase 5 (HS3ST5) gene and the teneurin transmembrane protein 4 (TENM4) gene, respectively (Table 2). For illustration, we plotted mean egg-laying dates against genotypes and spring temperature for these three SNPs (Figure 3). In cold springs the birds laid about 10 days later than in warm springs. The maximum difference in egg-laying dates between SNP genotypes was 0.94 days in cold springs, and 0.77 days in warm springs but 2.46 days in intermediate springs. In all three SNPs no genotype bred consistently earlier under all spring temperatures.

To compare model fit of the models with and without SNP^∗temperature interaction we calculated the marginal r² following Nakagawa and Schielzeth (2013), i.e., the variance explained by the fixed effects alone. Even after adjusting for the higher model complexity the model including the SNP^∗temperature interaction performed consistently better than the model fitting only the SNP effect (Figure 4).

Using data from 2045 wild great tits that were genotyped for more than 500k SNPs we could not identify any SNPs that were associated with egg-laying date. However, including an interaction between SNP effects and spring temperature improved model fit (Figure 4) although no SNP reached significance after multiple testing correction (Figure 2). This result shows that including the effects that the environment has on genotype-phenotype relationship can improve our understanding of this relationship. This has been realized previously in some research fields, as, e.g., animal and plant breeding or psychiatrical genetics (Kraft et al., 2007; e.g., Lindström et al., 2009; Korte and Farlow, 2013; Silva et al., 2014) but not in evolutionary genomics of wild populations.

Fitting an interaction between loci and an environmental variable can only identify loci that are responsible for trait variation in a given environment but not the loci that are responsible for phenotypic plasticity itself. Mapping the genes that are responsible for phenotypic plasticity by altering gene expression or post-translational mechanisms is possible if genetic variation in reaction norm slopes exist and reaction norms of individuals or genotypes, e.g., using clones, can be reliably quantified (e.g., Schadt et al., 2003; Tétard-Jones et al., 2011).

The phenotypic variance in egg-laying dates varied with temperature and to account for this heteroscedasticity we fitted heterogeneous residuals in more SNP^∗temperature interaction model, i.e., allowed the residual variance to differ among temperature classes. Not doing so would have let to inflated p-values because an important assumption of linear models was violated. This is a potential problem in many statistical analysis and may be more common but the p-value inflation often goes undetected and becomes only apparent when a (very) large number of tests are performed, as for example in GWAS. The inflation of p-values due to heteroscedasticity is likely more problematic when SNP^∗environment interactions are fitted but could theoretically also occur when only SNP effects are tested. It is normally assumed that variants at a locus affect the mean but they may also affect the variance of the trait (e.g., Ansel et al., 2008; Shen et al., 2012) and thereby lead to heteroscedasticity among genotypes. An observed inflation in p-values in GWAS may hence not only be due to population structure, which can be dealt with by applying ‘genomic control’ but also due to heteroscedasticity. In the latter case, however, applying ‘genomic control’ will not be sufficient and fitting heterogeneous residuals, as we did here, or other approaches (Voorman et al., 2011) to correct for this problem are needed.

Among the 10 most significant SNPs whose effects on breeding time differed depending on spring temperature in the year the phenotype was expressed, three were associated with genes of known function. One SNP is located within the thyroglobulin (TG) gene. Thyroglobulin is a precursor for thyroid hormones which form an integral part of the Hypothalamus–pituitary–thyroid axis, which has an important role in controlling general metabolism (Dietrich et al., 2012). In bulls the genetic variance in TG has been linked to age of puberty (Fernández et al., 2014) and thyroid hormones affect reproduction in both sexes (Flood et al., 2013). Another SNP was located within the heparan sulfate-glucosamine 3-sulfotransferase 5 (HS3ST5) gene, which is involved in heparan sulfate biosynthesis and glycosaminoglycan metabolism. Heparan sulfate is involved in cell-signaling and has various essential functions in development and homeostasis (Li and Kusche-Gullberg, 2016). Like for HS3ST5, Glycosaminoglycans are involved in cell-signaling and have structural functions in connective tissue, bone and blood vessels (Esko et al., 2009). The third SNP was located within the teneurin transmembrane protein 4 (TENM4). Teneurins are a highly conserved gene family that function as cell surface signal molecules and transcriptional regulators (Tucker and Chiquet-Ehrismann, 2006). TENM4 has been shown to play a role in neural development in chicken (Tucker et al., 2000).

Genome-wide association studies have worked well in certain fields, e.g., in case-control studies of human diseases (Visscher et al., 2012) but identifying genes underlying quantitative traits, especially in natural populations has proven to be much more difficult. Besides other reasons, as, e.g., unknown environmental effects on phenotypes that could not be accounted for, limited power or overly conservative multiple testing corrections, one potential reason for this is variation in phenotypes due to phenotypic plasticity. This variation can be accounted for if the environmental driver of the plasticity or the temporal or spatial scale on that it varies is known. For example, egg-laying dates vary strongly from year-to-year driven by spring temperatures. We can therefore account for this plasticity by correcting phenotypes for annual means; even if we had no idea which environmental variable was driving the annual variation. Phenotypic plasticity can, however, also ‘obscure’ genotype-phenotype associations because different loci are associated with the trait in different environments due to differential gene expression levels or post-translational mechanisms. An important aspect when testing whether genotype-phenotype associations vary across environments is obviously to identify the correct environmental variable that drives plastic variation in the trait. We here used spring temperature from March 11 to April 20, which has been shown to predict egg-laying dates fairly well in our study population (Husby et al., 2011), to classify our environment. However, using two temperature classes instead of three led to different results, in terms of the SNPs included in the best models. This highlights the problem of identifying a meaningful environmental variable when fitting SNP^∗environment interactions.

Most traits are phenotypically plastic (Pigliucci, 2001) and the lack of awareness of the variable genetic architecture underlying plastic traits may explain why identifying genes underlying quantitative traits has proven to be a major challenge. We here found indication for such variable SNP effects in plastic traits. This demonstrates that explicitly modeling phenotypic plasticity can be crucial for genome-wide association studies, especially – but not only – in studies of wild populations, which would thereby contribute to our still limited knowledge of the genome-phenome link in ecologically relevant traits, which is especially important as species need to adapt to their changing world.

PG and MV compiled the phenotypic data. VL compiled the genomic data. PG and VL conducted all analyses. AM organized the DNA samples and extracted the DNA. PG, VL, KvO, and MV developed this research, PG wrote the manuscript with contributions from VL and all co-authors commented on it.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.