Environment-Dependent Genotype-Phenotype Associations in Avian Breeding Time

Understanding how genes shape phenotypes is essential to assess the evolutionary potential of a trait. Identifying the genes underlying quantitative behavioral or life-history traits has, however, proven to be a major challenge. The majority of these traits are phenotypically plastic and different parts of the genome can be involved in shaping the trait under different environmental conditions. These variable genotype-phenotype associations could be one explanation for the limited success of genome-wide association studies in such traits. We here use avian seasonal timing of breeding, a trait that is highly plastic in response to spring temperature, to explore effects of such genotype-by-environment interactions in genome-wide association studies. We genotyped 2045 great tit females for 384081 single nucleotide polymorphisms (SNPs) and recorded their egg-laying dates in the wild. When testing for associations between SNPs and egg-laying dates, no SNP reached genome-wide significance. We then explored whether SNP effects were modified by annual spring temperature by formally testing for an interaction between SNP effect and temperature. The models including the SNP^∗temperature interaction performed consistently better although no SNP reached genome-wide significance. Our results suggest that the effects of genes shaping seasonal timing depended on annual spring temperature. Such environment-dependent effects are expected for any phenotypically plastic trait. Taking these effects into account will thus improve the success of detecting genes involved in phenotypically plastic traits, thereby leading to a better understanding of their evolutionary potential.

Introduction

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.

Materials and Methods

Phenotypes and Sample Collection

The timing of egg-laying of great tits (Parus major) has been studied in our long-term study populations on the ‘Veluwe’ area near Arnhem (52° 00′ N, 5° 50′ E, Netherlands) since the 1950s. Great tits are small, hole-breeding passerine birds that readily accept nest boxes, which has made them a ‘model system’ in behavioral and evolutionary ecology. In our study populations nest boxes are supplied in abundance so that almost the entire population breeds in boxes and can be monitored. Nest boxes were checked weekly in spring and the dates when the first egg of a clutch was laid (from now on: egg-laying date) were back-calculated from the number of recorded eggs assuming one egg is laid per day. All nestlings were banded with standard aluminum bird-bands and blood sampled when seven to 10-days-old. Adult great tits were caught when feeding their chicks with spring-traps in the nest box, identified from their bands, and blood-sampled. Whole blood samples were stored in either 1 ml Cell Lysis Solution (Gentra Puregene Kit, Qiagen, United States) or Queens buffer (Seutin et al., 1991). Field work, including capturing and blood-sampling birds, was carried out under a license of the Animal Experimental Committee of the Royal Dutch Academy of Sciences (KNAW) protocol NIOO-10.07.

Genotyping of Birds

DNA was extracted from blood samples by using the FavorPrep 96-Well Genomic DNA Extraction Kit (Favorgen Biotech corp.). DNA quality and DNA concentration were measured on a Nanodrop 2000 (Thermo Scientific). A total of 2311 female great tits were genotyped using a custom made Affymetrix^® great tit 650K SNP chip at Edinburgh Genomics (Edinburgh, United Kingdom). SNP calling was done following the Affymetrix^® best practices workflow by using the Axiom^® Analysis Suite 1.1. Eight individuals with dish quality control value of <0.82 were discarded. Dish quality control is an Affymetrix-specific QC measure and we used the default threshold to exclude individual samples. The recommended SNP group (PolyHighResolution, NoMinorHom, MonoHighResolution, CallRateBelowThreshold, Hemizygous) consisted of 537174 SNPs while 73796 SNPs were discarded. In addition to the SNPs that did not pass the quality control steps, an additional 234 SNPs were removed because they were duplicates or the genomic position was missing (NCBI Parus major genome version 1.1, GCA_001522545.2). Altogether 536940 SNPs passed initial quality control and 2303 individuals were included for downstream analyses.

Quality Control of Genotypes

The GenABEL package v1.8-0 (Aulchenko et al., 2007) implemented in R v3.2.0 (R Core Team, 2015) was used to perform a quality control on the dataset. 15719 SNPs and 24 individuals were discarded, because of being monomorphic and/or a low call rate (<95%). The mean autosomal heterozygosity was 0.342 and nine additional individuals were discarded because of a high heterozygosity (FDR < 1%). The mean identity by state of SNPs (IBS) was 0.696 and 14 individuals were discarded because of high IBS value (≥0.95). We did not exclude SNPs that deviated from Hardy-Weinberg equilibrium (HWE) because egg-laying date in our study populations is under directional selection and causal loci may therefore deviate from HWE. We used a MAF-threshold of 0.1 because we wanted to avoid small numbers of individuals for any genotype-temperature class group. We furthermore excluded all SNPs for that not all three genotypes were present and also SNPs that had less than 20 individuals in any genotype-temperature class group.

In PLINK 1.07 (Purcell et al., 2007) we calculated genomic relatedness, based on IBS, for all individuals across all SNPs and a genetic distance from this. This genetic distance matrix was then used, also in PLINK, for multi-dimensional scaling (MDS) to identify whether there was clustering of genetically similar individuals. Twenty-one individuals were identified as outliers and excluded (Supplementary Figure S1). In total 384081 SNPs and 2249 individuals passed the quality control. Given the size of the great tit genome of about 1.02 Gbp (Laine et al., 2015) we have a coverage of roughly 1 SNP per 2700 bp. Of these individuals 2045 had records for at least one egg-laying date and could be included in the analysis.

Statistical Analyses

Individual egg-laying dates were corrected for among-year and among-area variation prior to analyses by centering to the annual mean egg-laying date separately per area. To test whether SNP effects differed depending on spring temperature we divided the years in our data set into three categories according to their spring temperature. For every year we calculated the average daily mean temperature for the period March 11 to April 20. Average temperature during this period correlates best with annual mean egg-laying dates for our study population (Visser et al., 2006). The used daily mean temperatures from the weather station De Bilt were obtained from the website of the Royal Dutch Meteorological Institute (KNMI¹). We then calculated the 33%iles and grouped the first, second and third 33%ile into the ‘cold,’ ‘intermediate’ and ‘warm’ category, respectively. The ‘cold’ category comprised 1074 observed egg-laying date of 941 individuals in 6 years, the ‘intermediate’ category 770 observed egg-laying dates of 651 individuals in 6 years and the ‘warm’ category 1615 observed egg-laying dates of 1310 individuals in 5 years. Note that the total number of individuals in our analysis (2045) is not the sum of the individuals breeding in cold, intermediate and warm years, since the same individual can be included in more than one category. The average (and standard deviation in brackets) of the temperature during the period March 11 to April 20 was 6.2°C (0.96°C) for the ‘cold’ category, 7.8°C (0.27°C) for the ‘intermediate’ category, and 9.4°C (0.66°C) for the ‘warm’ category, respectively. Our model to test for SNP effect and their interaction with temperature was:

${\mathit{y}}_{\mathit{i}\mathit{,}\mathit{j}}\mathrm{=}\mathrm{\mu }\mathrm{+}{\mathrm{a}\mathrm{g}\mathrm{e}}_{\mathit{i}}\mathrm{+}{\mathrm{S}\mathrm{N}\mathrm{P}}_{\mathit{i}}\mathrm{+}{\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathit{i}}\mathrm{:}{\mathrm{S}\mathrm{N}\mathrm{P}}_{\mathit{i}}\mathrm{+}{\mathit{p}\mathit{e}}_{\mathit{i}}\mathrm{+}{\mathit{a}}_{\mathit{i}}\mathrm{+}{\mathit{e}}_{\mathit{i}\mathit{,}\mathit{j}}$

where the phenotype y_i,j of observation j of individual i is a function of the population mean μ, its age (two-level factor: first year breeder or older), its SNP genotype, the interaction between SNP genotype and temperature class temp, its ‘permanent environment’ (i.e., individual) effect pe_i, its breeding values a_i and the residual effect e_i. μ, age, the SNP genotype and temperature class are fixed effects and pe_i and a_i random effects. The variance-covariance matrix for a is derived from the expected covariance between individuals due to their additive genetic effects. Fitting the additive genetic variance-covariance matrix accounts for the fact that individuals in our study population are related.

All GWA analyses were run in ASReml-R 3.0 (Gilmour et al., 2015). Due to computation-time constraints we fitted the (sparse) relatedness matrix calculated from a pedigree instead of the (full) ‘genomic’ relatedness matrix. This should not affect the results as results for a simple model testing only the SNP effect (RepeatABEL cannot fit interactions) run in RepeatABEL (genomic relatedness matrix) and ASReml (pedigree relatedness matrix) correlated highly (Supplementary Figure S2).

An assumption underlying linear regression and ANOVA is that the variance of the dependent variable is homoscedastic with respect to the independent variable, i.e., that its variance does not change with the covariate in linear regression or that its variance is roughly equal in the different factor groups in ANOVA. If this assumption is violated, the resulting heteroscedasticity leads to biased results. This problem is exacerbated by fitting SNP^∗environment interactions and can lead to inflation of p-values as has been recognized in the field of psychiatric genetics (Voorman et al., 2011; Almli et al., 2014). To account for this problem we allowed the residual variance to differ among temperature classes in the models fitted with ASReml. Residual variances did indeed differ among temperature classes and fitting them avoided inflating p-values. See Supplementary Materials for further details.

Because the standard Bonferroni correction is generally considered overly conservative, an adjusted Bonferroni correction taking linkage disequilibrium (LD) between all SNP loci into account has been suggested (Moskvina and Schmidt, 2008). We calculated this correction with a sliding window of 20 SNPs using the software KEFF VSEP 2007 (Moskvina and Schmidt, 2008). However, because of low LD between SNPs in our great tit population, the amount of effective tests, on which the Bonferroni correction depends, dropped by just 4% and we hence did not apply this correction.

Although previous work in of the study populations analyzed here showed that egg-laying dates are heritable in general (Gienapp et al., 2006; Husby et al., 2010) and also under different spring temperatures (Husby et al., 2011), we here analyzed whether egg-laying date would be heritable in each of the three temperature classes. To increase the power of our analysis, we not only included the genotyped individuals but all females with known identity in our analysis. This data set consisted of 4624 observations of 4032 females, 3737 observations of 3019 females and 4147 of 3532 females for the cold, medium and warm temperature class, respectively. We hence had to include a pedigree-based relatedness matrix instead of a genomic relatedness matrix in our animal model. Besides the additive genetic (random) effect, we fitted female identity as random effect to account for repeated observations as well as age (two level factor as above), year (as factor), area, and the interaction year^∗area as fixed effects. These analyses were also run in ASReml-R 3.0.

Results

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).

TABLE 1

TABLE 1. Variance components (and SE) from quantitative genetic analysis of egg-laying date, separately by temperature class.

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.

FIGURE 1

FIGURE 1. Significances (shown as –log10) of all SNPs included in the genome-wide association analyses of egg-laying date assuming temperature-independent SNP effects. The black line indicates the genome-wide significance level equivalent to P = 0.05 after applying a Bonferroni correction. Points are color-coded by chromosome.

FIGURE 2

FIGURE 2. Significances (shown as –log10) of all SNPs included in the genome-wide association analyses of egg-laying date testing for an interaction between temperature and SNP effects on laying date. The black line indicates the genome-wide significance level equivalent to P = 0.05 after applying a Bonferroni correction. Points are color-coded by chromosome.

TABLE 2

TABLE 2. The 10 most significant SNPs for genome-wide association analyses testing for an interaction between SNP effects and temperature, ordered by significance.

FIGURE 3

FIGURE 3. Variation in egg-laying dates of SNP genotypes with temperature for the three SNPs with the most significant interaction with temperature (Table 2), AX-100216683 (A), AX-100724221 (B) and AX-100140360 (C). Shown are means and standard errors, indicated by whiskers, for homozygote 1 (light gray), heterozygote (medium gray) and homozygote 2 (black). Plotting the area- and year-centered egg-laying dates that were used in the analysis would give a horizontal mean reaction norm. To show the real differences in egg-laying dates with temperature, we here plotted ‘raw’ egg-laying dates that were not corrected for year- and area effects. Note that therefore part of the variation within SNP- and temperature-classes is due to variation among areas. Sample sizes per temperature class-genotype group ranged between 77 and 757 individuals and were on average 381.2.

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).

FIGURE 4

FIGURE 4. Adjusted marginal r² of models with and without SNP^∗temperature interaction. The marginal r² measures the variation explained by fixed effects only and is hence appropriate here because the random effects structure of the models is identical. The r² is only plotted for models in which either the SNP effect or the interaction SNP^∗temperature was significant at the 0.05 level.

Discussion

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.

Author Contributions

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.

Funding

Part of this work was funded by an ERC Advanced Grant (339092 - E-Response) to MV.

Conflict of Interest Statement

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.

Acknowledgments

Over the years many people were involved with the extensive field work for data collection. The National Park ‘Hoge Veluwe’, Staatsbosbeheer, Geldersch Landschap & Kasteelen, and the van Boetzelaer family kindly gave us permission to work on their properties. Marieke Beltman helped with data handling. Jun-Mo Kim and Jon Slate were involved in developing the SNP chip we used for our analysis. Mirte Bosse, Martien Groenen, Ilse van Grevenhof and two reviewers commented on the manuscript and gave insightful advice.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fgene.2017.00102/full#supplementary-material

Footnotes

1. ^https://data.knmi.nl/datasets

References

Almli, L. M., Duncan, R., Feng, H., Ghosh, D., Binder, E. B., Bradley, B., et al. (2014). Correcting systematic inflation in genetic association tests that consider interaction effects. JAMA Psychiatry 71, 1392–1399. doi: 10.1001/jamapsychiatry.2014.1339

PubMed Abstract | CrossRef Full Text | Google Scholar

Ansel, J., Bottin, H., Rodriguez-Beltran, C., Damon, C., Nagarajan, M., Fehrmann, S., et al. (2008). Cell-to-cell stochastic variation in gene expression is a complex genetic trait. PLoS Genet. 4:e1000049. doi: 10.1371/journal.pgen.1000049

PubMed Abstract | CrossRef Full Text | Google Scholar

Aulchenko, Y. S., Ripke, S., Isaacs, A., and van Duijn, C. M. (2007). GenABEL: an R library for genome-wide association analysis. Bioinformatics 23, 1294–1296. doi: 10.1093/bioinformatics/btm108

PubMed Abstract | CrossRef Full Text | Google Scholar

Barson, N. J., Aykanat, T., Hindar, K., Baranski, M., Bolstad, G. H., Fiske, P., et al. (2015). Sex-dependent dominance at a single locus maintains variation in age at maturity in salmon. Nature 528, 405–408. doi: 10.1038/nature16062

PubMed Abstract | CrossRef Full Text | Google Scholar

Bürger, R., and Lynch, M. (1995). Evolution and extinction in a changing environment: a quantitative-genetic analysis. Evolution 49, 151–163. doi: 10.1111/j.1558-5646.1995.tb05967.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Charmantier, A., McCleery, R. H., Cole, L. R., Perrins, C., Kruuk, L. E. B., and Sheldon, B. C. (2008). Adaptive phenotypic plasticity in response to climate change in a wild bird population. Science 320, 800–803. doi: 10.1126/science.1157174

PubMed Abstract | CrossRef Full Text | Google Scholar

Dietrich, J. W., Landgrafe, G., and Fotiadou, E. H. (2012). TSH and thyrotropic agonists: key actors in thyroid homeostasis. J. Thyroid Res. 2012:351864. doi: 10.1155/2012/351864

PubMed Abstract | CrossRef Full Text | Google Scholar

Esko, J. D., Kimata, K., and Lindahl, U. (2009). “Proteoglycans and sulfated glycosaminoglycans,” in Essentials of Glycobiology, eds A. Varki, R. D. Cummings, J. D. Esko, H. H. Freeze, P. Stanley, C. R. Bertozzi, et al. (Cold Spring Harbor, NY: Cold Spring Harbor Laboratory Press).

Google Scholar

Fernández, M. E., Goszczynski, D. E., Prando, A. J., Peral-García, P., Baldo, A., Giovambattista, G., et al. (2014). Assessing the association of single nucleotide polymorphisms in thyroglobulin gene with age of puberty in bulls. J. Anim. Sci. Technol. 56:17. doi: 10.1186/2055-0391-56-17

PubMed Abstract | CrossRef Full Text | Google Scholar

Flood, D. E. K., Fernandino, J. I., and Langlois, V. S. (2013). Thyroid hormones in male reproductive development: evidence for direct crosstalk between the androgen and thyroid hormone axes. Gen. Comp. Endocrinol. 192, 2–14. doi: 10.1016/j.ygcen.2013.02.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Gienapp, P., Hemerik, L., and Visser, M. E. (2005). A new statistical tool to predict phenology under climate change scenarios. Glob. Change Biol. 11, 600–606. doi: 10.1111/j.1365-2486.2005.00925.x

CrossRef Full Text | Google Scholar

Gienapp, P., Lof, M., Reed, T. E., McNamara, J., Verhulst, S., and Visser, M. E. (2013). Predicting demographically sustainable rates of adaptation: Can great tit breeding time keep pace with climate change? Philos. Trans. R. Soc. Lond. B Biol. Sci. 368:20120289. doi: 10.1098/rstb.2012.0289

PubMed Abstract | CrossRef Full Text | Google Scholar

Gienapp, P., Postma, E., and Visser, M. E. (2006). Why breeding time has not responded to selection for earlier breeding in a songbird population. Evolution 60, 2381–2388. doi: 10.1111/j.0014-3820.2006.tb01872.x

CrossRef Full Text | Google Scholar

Gienapp, P., Reed, T. E., and Visser, M. E. (2014). Why climate change will invariably alter selection pressures on phenology. Proc. Biol. Sci. 281:20141611. doi: 10.1098/rspb.2014.1611

PubMed Abstract | CrossRef Full Text | Google Scholar

Gienapp, P., Väisänen, R. A., and Brommer, J. E. (2010). Latitudinal variation in phenotypic plasticity of avian breeding time. J. Anim. Ecol. 79, 836–842.

Gilmour, A. R., Gogel, B. J., Cullis, B. R., Welham, S. J., and Thompson, R. (2015). ASReml User Guide Release 4.1 Functional Specification. Hemel Hempstead: VSN International Ltd.

Google Scholar

Gonzalez, A., Ronce, O., Ferriere, R., and Hochberg, M. E. (2013). Evolutionary rescue: an emerging focus at the intersection between ecology and evolution. Philos. Trans. R. Soc. Lond. B Biol. Sci. 268:20120404. doi: 10.1098/rstb.2012.0404

PubMed Abstract | CrossRef Full Text | Google Scholar

Huntzinger, E., and Izaurralde, E. (2011). Gene silencing by microRNAs: contributions of translational repression and mRNA decay. Nat. Rev. Genet. 12, 99–110. doi: 10.1038/nrg2936

PubMed Abstract | CrossRef Full Text | Google Scholar

Husby, A., Kawakami, T., Rönnegård, L., Smeds, L., Ellegren, H., and Qvarnström, A. (2015). Genome-wide association mapping in a wild avian population identifies a link between genetic and phenotypic variation in a life-history trait. Proc. Biol. Sci. 282:20150156. doi: 10.1098/rspb.2015.0156

PubMed Abstract | CrossRef Full Text | Google Scholar

Husby, A., Nussey, D. H., Visser, M. E., Wilson, A. J., Sheldon, B. C., and Kruuk, L. E. B. (2010). Contrasting patterns pf phenotypic plasticity in reproductive traits in two great tit (Parus major) populations. Evolution 64, 2221–2237. doi: 10.1111/j.1558-5646.2010.00991.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Husby, A., Visser, M. E., and Kruuk, L. E. B. (2011). Speeding up microevolution: the effects of increasing temperature on selection and genetic variance in a wild bird population. PLoS Biol. 9:e1000585. doi: 10.1371/journal.pbio.1000585

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnston, S. E., Bérénos, C., Slate, J., and Pemberton, J. M. (2016). Conserved genetic architecture underlying individual recombination rate variation in a wild population of Soay sheep (Ovis aries). Genetics 203, 583–598. doi: 10.1534/genetics.115.185553

PubMed Abstract | CrossRef Full Text | Google Scholar

Korte, A., and Farlow, A. (2013). The advantages and limitations of trait analysis with GWAS: a review. Plant Methods 9:29. doi: 10.1186/1746-4811-9-29

PubMed Abstract | CrossRef Full Text | Google Scholar

Kraft, P., Yen, Y. C., Stram, D. O., Morrison, J., and Gauderman, W. J. (2007). Exploiting gene-environment interaction to detect genetic associations. Hum. Hered. 63, 111–119. doi: 10.1159/000099183

PubMed Abstract | CrossRef Full Text | Google Scholar

Laine, V. N., Gossmann, T. I., Schachtschneider, K. M., Garroway, C. J., Madsen, O., Verhoeven, K. J. F., et al. (2015). Evolutionary signals of selection on cognition from the great tit genome and methylome. Nat. Commun. 7:10474. doi: 10.1038/ncomms10474

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambrechts, M. M., Blondel, J., Maistre, M., and Perret, P. (1997). A single response mechanism is responsible for evolutionary adaptive variation in a bird’s laying date. Proc. Natl. Acad. Sci. U.S.A. 94, 5153–5155. doi: 10.1073/pnas.94.10.5153

CrossRef Full Text | Google Scholar

Lambrechts, M. M., and Perret, P. (2000). A long photoperiod overrides non-photoperiodic factors in blue tits’ timing of reproduction. Proc. R. Soc. Lond. B Biol. Sci. 267, 585–588. doi: 10.1098/rspb.2000.1041

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J. P., and Kusche-Gullberg, M. (2016). Heparan Sulfate: biosynthesis, structure, and function. Int. Rev. Cell Mol. Biol. 325, 215–273. doi: 10.1016/bs.ircmb.2016.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindström, S., Yen, Y.-C., Spiegelman, D., and Kraft, P. (2009). The impact of gene-environment dependence and misclassification in genetic association studies incorporating gene-environment interactions. Hum. Hered. 68, 171–181. doi: 10.1159/000224637

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackay, T. F., Stone, E. A., and Ayroles, J. F. (2009). The genetics of quantitative traits: challenges and prospects. Nat. Rev. Genet. 10, 565–577. doi: 10.1038/nrg2612

PubMed Abstract | CrossRef Full Text | Google Scholar

Moskvina, V., and Schmidt, K. M. (2008). On multiple-testing correction in genome-wide association studies. Genet. Epidemiol. 32, 567–573. doi: 10.1002/gepi.20331

PubMed Abstract | CrossRef Full Text | Google Scholar

Mousseau, T. A., and Roff, D. A. (1987). Natural selection and the heritability of fitness components. Heredity 59, 181–198. doi: 10.1038/hdy.1987.113

CrossRef Full Text | Google Scholar

Nakagawa, S., and Schielzeth, H. (2013). A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods Ecol. Evol. 4, 133–142. doi: 10.1111/j.2041-210x.2012.00261.x

CrossRef Full Text | Google Scholar

Nussey, D. H., Wilson, A. J., and Brommer, J. E. (2007). The evolutionary ecology of individual phenotypic plasticity in wild populations. J. Evol. Biol. 20, 831–844. doi: 10.1111/j.1420-9101.2007.01300.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pigliucci, M. (2001). Phenotypic Plasticity. Baltimore, MD: John Hopkins University Press.

Google Scholar

Pigliucci, M. (2005). Evolution of phenotypic plasticity: where are we going now? Trends Ecol. Evol. 20, 481–486.

PubMed Abstract | Google Scholar

Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A. R., Bender, D., et al. (2007). PLINK: a toolset for whole-genome association and population-based linkage analysis. Am. J. Hum. Genet. 81, 559–575. doi: 10.1086/519795

PubMed Abstract | CrossRef Full Text | Google Scholar

R Core Team (2015). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.

Google Scholar

Schadt, E. E., Monks, S. A., Drake, T. A., Lusis, A. J., Che, N., Colinayo, V., et al. (2003). Genetics of gene expression surveyed in maize, mouse and man. Nature 422, 297–302. doi: 10.1038/nature01434

PubMed Abstract | CrossRef Full Text | Google Scholar

Seutin, G., White, B. N., and Boag, P. T. (1991). Preservation of avian blood and tissue samples for DNA analyses. Can. J. Zool. 69, 82–90. doi: 10.1139/z91-013

CrossRef Full Text

Shen, X., Petterson, M., Rönnegård, L., and Carlborg, Ö. (2012). Inheritance beyond plain heritability: variance-controlling genes in Arabidopsis thaliana. PLoS Genet. 8:e1002839. doi: 10.1371/journal.pgen.1002839

PubMed Abstract | CrossRef Full Text | Google Scholar

Silva, F. F., Mulder, H. A., Knol, E. F., Lopes, M. S., Guimarães, S. E. F., Lopes, P. S., et al. (2014). Sire evaluation for total number born in pigs using a genomic reaction norms approach. J. Anim. Sci. 92, 3825–3834. doi: 10.2527/jas.2013-6486

PubMed Abstract | CrossRef Full Text | Google Scholar

Slate, J., Visscher, P. M., MacGregor, S., Stevens, D., Tate, M. L., and Pemberton, J. M. (2002). A genome scan for quantitative trait loci in a wild population of red deer (Cervus elaphus). Genetics 162, 1863–1873.

PubMed Abstract | Google Scholar

Slate, L., Santure, A. E., Feulner, P. G. D., Brown, E. A., Ball, A. D., Johnston, S. E., et al. (2010). Genome mapping in intensively studied wild vertebrate populations. Trends Genet. 26, 275–284. doi: 10.1016/j.tig.2010.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Snell-Rood, E. C., Van Dyken, J. D., Cruikshank, T., Wade, M. J., and Moczek, A. P. (2010). Toward a population genetic framework of developmental evolution: the costs, limits, and consequences of phenotypic plasticity. Bioessays 32, 71–81. doi: 10.1002/bies.200900132

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarka, M., Åkesson, M., Beraldi, D., Hernández-Sánchez, J., Hasselquist, D., Bensch, S., et al. (2010). A strong quantitative trait locus for wing length on chromosome 2 in a wild population of great reed warblers. Proc. Biol. Sci. 277, 2361–2369. doi: 10.1098/rspb.2010.0033

PubMed Abstract | CrossRef Full Text | Google Scholar

Tétard-Jones, C., Kertesz, M. A., and Preziosi, R. F. (2011). Quantitative trait loci mapping of phenotypic plasticity and genotype-environment interactions in plant and insect performance. Philos. Trans. R. Soc. Lond. B Biol. Sci. 366, 1368–1379. doi: 10.1098/rstb.2010.0356

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomas, D. (2010). Gene–environment-wide association studies: emerging approaches. Nat. Rev. Genet. 11, 259–272. doi: 10.1038/nrg2764

PubMed Abstract | CrossRef Full Text | Google Scholar

Tucker, R. P., and Chiquet-Ehrismann, R. (2006). Teneurins: a conserved family of transmembrane proteins involved in intercellular signaling during development. Dev. Biol. 290, 237–245. doi: 10.1016/j.ydbio.2005.11.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Tucker, R. P., Martin, D., Kos, R., and Chiquet-Ehrismann, R. (2000). The expression of teneurin-4 in the avian embryo. Mech. Dev. 98, 187–191. doi: 10.1016/S0925-4773(00)00444-5

CrossRef Full Text | Google Scholar

Visscher, P. M., Brown, M. A., McCarthy, M. I., and Yang, J. (2012). Five years of GWAS discovery. Am. J. Hum. Genet. 90, 7–24. doi: 10.1016/j.ajhg.2011.11.029

PubMed Abstract | CrossRef Full Text | Google Scholar

Visser, M. E. (2008). Keeping up with a warming world; assessing the rate of adaptation to climate change. Proc. Biol. Sci. 275, 649–659. doi: 10.1098/rspb.2007.0997

PubMed Abstract | CrossRef Full Text | Google Scholar

Visser, M. E., and Both, C. (2005). Shifts in phenology due to global climate change: the need for a yardstick. Proc. R. Soc. Lond. B Biol. Sci. 272, 2561–2569. doi: 10.1098/rspb.2005.3356

PubMed Abstract | CrossRef Full Text | Google Scholar

Visser, M. E., Holleman, L. J. M., and Gienapp, P. (2006). Shifts in caterpillar biomass phenology due to climate change and its impact on the breeding biology of an insectivorous bird. Oecologia 147, 164–172. doi: 10.1007/s00442-005-0299-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Visser, M. E., van Noordwijk, A. J., Tinbergen, J. M., and Lessells, C. M. (1998). Warmer springs lead to mistimed reproduction in great tits (Parus major). Proc. R. Soc. Lond. B Biol. Sci. 265, 1867–1870. doi: 10.1098/rspb.1998.0514

CrossRef Full Text | Google Scholar

Voorman, A., Lumley, T., McKnight, B., and Rice, K. (2011). Behavior of QQ-plots and genomic control in studies of gene-environment interaction. PLoS ONE 6:e19416. doi: 10.1371/journal.pone.0019416

PubMed Abstract | CrossRef Full Text | Google Scholar

Weake, V. M., and Workman, J. L. (2010). Inducible gene expression: diverse regulatory mechanisms. Nat. Rev. Genet. 11, 429–437. doi: 10.1038/nrg2781

PubMed Abstract | CrossRef Full Text | Google Scholar

Wenzel, M. A., James, M. C., Douglas, A., and Piertney, S. B. (2015). Genome-wide association and genome partitioning reveal novel genomic regions underlying variation in gastrointestinal nematode burden in a wild bird. Mol. Ecol. 24, 4175–4192. doi: 10.1111/mec.13313

PubMed Abstract | CrossRef Full Text | Google Scholar