# Divergence in evolutionary potential of life history traits among wild populations is predicted by differences in climatic conditions

Stéphane Chantepie, Anne Charmantier, Boris Delahaie, Frank Adriaensen, Erik Matthysen, Marcel Visser, Elena Álvarez, Emilio Barba, Markku Orell, Ben Sheldon, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Stéphane Chantepie, Anne Charmantier, Boris Delahaie, Frank Adriaensen, Erik Matthysen, et al. "Divergence in evolutionary potential of life history traits among wild populations is predicted by differences in climatic conditions." Evolution Letters, 2024, 8 (1), pp.29-42. 10.1093/evlett/qrad067. hal-04548215

HAL Id: hal-04548215

https://hal.science/hal-04548215

Submitted on 26 Apr 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Divergence in evolutionary potential of life history traits among wild populations is predicted by differences in climatic conditions

Stéphane Chantepie, Anne Charmantier, Boris Delahaie, Frank Adriaensen, Erik Matthysen, Marcel Visser, Elena Álvarez, Emilio Barba, Markku Orell, Ben Sheldon, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Stéphane Chantepie, Anne Charmantier, Boris Delahaie, Frank Adriaensen, Erik Matthysen, et al. Divergence in evolutionary potential of life history traits among wild populations is predicted by differences in climatic conditions. Evolution Letters, 2024, 8 (1), pp.29-42. 10.1093/evlett/qrad067. hal-04548215

HAL Id: hal-04548215

https://hal.science/hal-04548215v1

Submitted on 26 Apr 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Divergence in evolutionary potential of life history traits among wild populations is predicted by differences in climatic conditions

Stéphane Chantepie<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, Anne Charmantier<fnref n="4" />, Boris Delahaie<fnref n="5" />,<fnref n="6" />, Frank Adriaensen<fnref n="7" />, Erik Matthysen<fnref n="7" />, Marcel E. Visser<fnref n="8" />, Elena Álvarez<fnref n="9" />, Emilio Barba<fnref n="9" />, Markku Orell<fnref n="10" />, Ben Sheldon<fnref n="11" />, Elena Ivankina<fnref n="12" />, Anvar Kerimov<fnref n="13" /> ID, Sébastien Lavergne<fnref n="1" />, Céline Teplitsky<fnref n="4" />

Corresponding author: INRAE, ONF, BIOFORA, 2163 Avenue de la Pomme de Pin, CS 40001 Ardon, 45075, Orléans, France. Email: stephane.chantepie@inrae.fr

## Abstract

Short-term adaptive evolution represents one of the primary mechanisms allowing species to persist in the face of global change. Predicting the adaptive response at the species level requires reliable estimates of the evolutionary potential of traits involved in adaptive responses, as well as understanding how evolutionary potential varies across a species' range. Theory suggests that spatial variation in the fitness landscape due to environmental variation will directly impact the evolutionary potential of traits. However, empirical evidence on the link between environmental variation and evolutionary potential across a species range in the wild is lacking. In this study, we estimate multivariate evolutionary potential (via the genetic variance-covariance matrix, or G-matrix) for six morphological and life history traits in 10 wild populations of great tits (*Parus major*) distributed across Europe. The G-matrix significantly varies in size, shape, and orientation across populations for both types of traits. For life history traits, the differences in G-matrix are larger when populations are more distant in their climatic niche. This suggests that local climates contribute to shaping the evolutionary potential of phenotypic traits that are strongly related to fitness. However, we found no difference in the overall evolutionary potential (i.e., G-matrix size) between populations closer to the core or the edge of the distribution area. This large-scale comparison of G-matrices across wild populations emphasizes that integrating variation in multivariate evolutionary potential is important to understand and predict species' adaptive responses to new selective pressures.

**Keywords:** G-matrix, wild populations, climatic niche, great tits

## Lay Summary

Adaptive evolution of phenotypes is a key mechanism for species persistence in the face of global change, and predicting these responses is a major challenge to implementing biodiversity scenarios. The quality of these predictions relies notably on our estimates of evolutionary potential, and several questions need to be answered: How different is evolutionary potential across a species distribution? Will it change under new environmental conditions? Can we use information on evolutionary potential from one population to infer adaptive responses in another population? Here we address these questions in the great tit (*Parus major*) by evaluating how evolutionary potential is influenced by climatic conditions across the species' range. Our study is based on the long-term monitoring of 10 wild populations of great tits distributed across Europe for which information about life history and morphological traits have been collected. We find significant differences in evolutionary potential between populations for both categories of traits, and for life history traits, these differences are linked to population differences in climatic conditions. These results show that evolutionary potential is not a "one-size-fits-all" measure and that its variation needs to be taken into account when predicting adaptive responses to climate change.

<footnote n="1">¹Université Grenoble Alpes, CNRS, Univ. Savoie Mont Blanc, CNRS, LECA, Laboratoire d'Écologie Alpine, Grenoble, France</footnote>

<footnote n="2">²Institut de Systématique, Évolution, Biodiversité, École Pratique des Hautes Études, Paris Sciences et Lettres, Muséum National d'Histoire Naturelle, CNRS, Sorbonne Université, Université des Antilles, Paris, France</footnote>

<footnote n="3">³INRAE, ONF, BIOFORA, 45075, Orléans, France</footnote>

<footnote n="4">⁴CEFE, Univ Montpellier, CNRS, EPHE, IRD, Montpellier, France</footnote>

<footnote n="5">⁵CIRAD, UMR DIADE, Montpellier 34398, France</footnote>

<footnote n="6">⁶UMR DIADE, Université de Montpellier, CIRAD, IRD, Montpellier, France</footnote>

<footnote n="7">⁷Evolutionary Ecology Group, Department of Biology, University of Antwerp, Wilrijk, Belgium</footnote>

<footnote n="8">⁸Department of Animal Ecology, Netherlands Institute of Ecology, Wageningen, The Netherlands</footnote>

<footnote n="9">⁹Cavanilles Institute of Biodiversity and Evolutionary Biology, Department of Microbiology and Ecology, University of Valencia, Paterna, Spain</footnote>

<footnote n="10">¹⁰Ecology and Genetics Research Unit, Faculty of Science, University of Oulu, Oulu, Finland</footnote>

<footnote n="11">¹¹Edward Grey Institute, Department of Biology, University of Oxford, Oxford, United Kingdom</footnote>

<footnote n="12">¹²Zvenigorod Biological Station of Lomonosov Moscow State University, Moscow, Russia</footnote>

<footnote n="13">¹³Department of Vertebrate Zoology, Faculty of Biology, Lomonosov Moscow State University, Moscow, Russia</footnote>
# Divergence in evolutionary potential of life history traits among wild populations is predicted by differences in climatic conditions

Stéphane Chantepie<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, Anne Charmantier<fnref n="4" />, Boris Delahaie<fnref n="5" />,<fnref n="6" />, Frank Adriaensen<fnref n="7" />, Erik Matthysen<fnref n="7" />, Marcel E. Visser<fnref n="8" />, Elena Álvarez<fnref n="9" />, Emilio Barba<fnref n="9" />, Markku Orell<fnref n="10" />, Ben Sheldon<fnref n="11" />, Elena Ivankina<fnref n="12" />, Anvar Kerimov<fnref n="13" /> ID, Sébastien Lavergne<fnref n="1" />, Céline Teplitsky<fnref n="4" />

Corresponding author: INRAE, ONF, BIOFORA, 2163 Avenue de la Pomme de Pin, CS 40001 Ardon, 45075, Orléans, France. Email: stephane.chantepie@inrae.fr

## Abstract

Short-term adaptive evolution represents one of the primary mechanisms allowing species to persist in the face of global change. Predicting the adaptive response at the species level requires reliable estimates of the evolutionary potential of traits involved in adaptive responses, as well as understanding how evolutionary potential varies across a species' range. Theory suggests that spatial variation in the fitness landscape due to environmental variation will directly impact the evolutionary potential of traits. However, empirical evidence on the link between environmental variation and evolutionary potential across a species range in the wild is lacking. In this study, we estimate multivariate evolutionary potential (via the genetic variance-covariance matrix, or G-matrix) for six morphological and life history traits in 10 wild populations of great tits (*Parus major*) distributed across Europe. The G-matrix significantly varies in size, shape, and orientation across populations for both types of traits. For life history traits, the differences in G-matrix are larger when populations are more distant in their climatic niche. This suggests that local climates contribute to shaping the evolutionary potential of phenotypic traits that are strongly related to fitness. However, we found no difference in the overall evolutionary potential (i.e., G-matrix size) between populations closer to the core or the edge of the distribution area. This large-scale comparison of G-matrices across wild populations emphasizes that integrating variation in multivariate evolutionary potential is important to understand and predict species' adaptive responses to new selective pressures.

**Keywords:** G-matrix, wild populations, climatic niche, great tits

## Lay Summary

Adaptive evolution of phenotypes is a key mechanism for species persistence in the face of global change, and predicting these responses is a major challenge to implementing biodiversity scenarios. The quality of these predictions relies notably on our estimates of evolutionary potential, and several questions need to be answered: How different is evolutionary potential across a species distribution? Will it change under new environmental conditions? Can we use information on evolutionary potential from one population to infer adaptive responses in another population? Here we address these questions in the great tit (*Parus major*) by evaluating how evolutionary potential is influenced by climatic conditions across the species' range. Our study is based on the long-term monitoring of 10 wild populations of great tits distributed across Europe for which information about life history and morphological traits have been collected. We find significant differences in evolutionary potential between populations for both categories of traits, and for life history traits, these differences are linked to population differences in climatic conditions. These results show that evolutionary potential is not a "one-size-fits-all" measure and that its variation needs to be taken into account when predicting adaptive responses to climate change.

<footnote n="1">¹Université Grenoble Alpes, CNRS, Univ. Savoie Mont Blanc, CNRS, LECA, Laboratoire d'Écologie Alpine, Grenoble, France</footnote>

<footnote n="2">²Institut de Systématique, Évolution, Biodiversité, École Pratique des Hautes Études, Paris Sciences et Lettres, Muséum National d'Histoire Naturelle, CNRS, Sorbonne Université, Université des Antilles, Paris, France</footnote>

<footnote n="3">³INRAE, ONF, BIOFORA, 45075, Orléans, France</footnote>

<footnote n="4">⁴CEFE, Univ Montpellier, CNRS, EPHE, IRD, Montpellier, France</footnote>

<footnote n="5">⁵CIRAD, UMR DIADE, Montpellier 34398, France</footnote>

<footnote n="6">⁶UMR DIADE, Université de Montpellier, CIRAD, IRD, Montpellier, France</footnote>

<footnote n="7">⁷Evolutionary Ecology Group, Department of Biology, University of Antwerp, Wilrijk, Belgium</footnote>

<footnote n="8">⁸Department of Animal Ecology, Netherlands Institute of Ecology, Wageningen, The Netherlands</footnote>

<footnote n="9">⁹Cavanilles Institute of Biodiversity and Evolutionary Biology, Department of Microbiology and Ecology, University of Valencia, Paterna, Spain</footnote>

<footnote n="10">¹⁰Ecology and Genetics Research Unit, Faculty of Science, University of Oulu, Oulu, Finland</footnote>

<footnote n="11">¹¹Edward Grey Institute, Department of Biology, University of Oxford, Oxford, United Kingdom</footnote>

<footnote n="12">¹²Zvenigorod Biological Station of Lomonosov Moscow State University, Moscow, Russia</footnote>

<footnote n="13">¹³Department of Vertebrate Zoology, Faculty of Biology, Lomonosov Moscow State University, Moscow, Russia</footnote>
Introduction

Along with phenotypic plasticity and range shifts toward suit-
able environments, evolution represents a key mechanism for
natural populations to persist under consistently adverse envi-
ronmental conditions (Chevin et al., 2010). The rate and direc-
tion of evolutionary response critically depend on both natural
selection and available standing genetic variance (Walsh & Lynch,
2018). The joint response to selection of multiple traits can be
facilitated or constrained by the genetic association among traits
(Agrawal & Stinchcombe, 2009; Bolstad et al., 2014; Schluter,
1996; Stinchcombe & Kirkpatrick, 2012; Teplitsky et al., 2014).
Under this framework, the matrix of additive genetic variances
and covariances (G-matrix), which summarizes the multivariate
standing genetic variance, plays a fundamental role in predict-
ing the response to selection as the adaptive response across a
single generation is given by the product of the G-matrix and the
selection gradient (Lande, 1979). The extrapolation of the adap-
tive response to selection over multiple generations or across
space thus largely depends on the constancy of the G-matrix
in time and space. The stability of the G-matrix is a paramount
assumption, prevalent in quantitative genetic models (Duputié et
al., 2012).

The structure of G-matrices results from the interactive effects
of selection, mutation, drift, and migration processes (Steppan et
al., 2002). Since these processes, which jointly affect allele fre-
quencies across populations, can change both in space and time,
theory predicts that G-matrices should also vary (Chantepie &
Chevin, 2020). Spatiotemporal changes in the G-matrix can be
best interpreted in geometrical terms (Jones et al., 2004) such as
their variation in size (overall evolutionary potential), orientation
 dir ect of fastest evolution), and shape (constraint resulting
from genetic correlation). In finite populations, the shape of the
G-matrix is a compromise between the shape and orientation of
the matrix of selective gradients and the matrix of pleiotropic
mutations (Jones et al., 2007, 2014; Revell, 2007). At mutation-
selection equilibrium, when drift is high, the G-matrix tends to
orient toward the direction of the mutational matrix (Chantepie &
Chevin, 2020). Differences in these forces among populations will
contribute to differentiation among G-matrices, while gene flow
is expected to homogenize allele frequencies across populations
hence reducing the spatial variation in G-matrices (Guillaume &
Whitlock, 2007; Lande, 1992). Finally, drift may also lead to ran-
dom divergence of G-matrices among populations (Phillips et
al., 2001). Experimental studies conducted over the last decades
have demonstrated that each of the above-mentioned evolution-
ary forces can significantly affect the geometry of the G-matrix
(examples for drift: Phillips et al., 2001, selection: Cano et al.,
2004; Careau et al., 2015; Doroszuk et al., 2008; Shaw et al., 1995,
dispersal: Nosil et al., 2006). In addition to evolutionary forces,
gene expression depending on the environment (G × E) can gen-
erate variation in the G-matrix independently of changes in allele
frequencies (Wood & Brodie, 2015).

Although experimental studies are valuable for understanding how evolutionary forces can shape the G-matrix and testing theoretical hypotheses, the extrapolation of laboratory results to wild populations living under natural conditions remains challenging (Kruuk et al., 2008). In particular, the difficulty to measure precisely each of these evolutionary forces as well as their interactions makes it impossible to predict theoretically the variation in G-matrices in the wild (Arnold et al., 2008; Roff, 2000). Investigating this variation among populations can provide inference on forces driving spatial differences but also, using space

for time substitution, it can contribute to improve predictions of
temporal changes in G-matrices. Addressing the spatial variation
in the G-matrix requires thus not only to characterize the extent
to which the G-matrix varies among populations but also the
spatial determinants (e.g., geographic distance, environmental
conditions) that have shaped this variation. So far, the literature
relative to G-matrix variation among populations has provided
contrasting results with evidence for both stability (Ashman,
2003; Brodie, 1993; Delahaie et al., 2017; Hangartner et al., 2020;
Henry & Stinchcombe, 2023; Puentes et al., 2016; Sniegula et
al., 2018) and variation (Berger et al., 2013; Eroukhmanoff &
Svensson, 2011; Paccard et al., 2016; Roff et al., 2004; Teplitsky et
al., 2011) of the G-matrix across various spatial gradients in plant
and animal populations. The extent to which differing environ-
mental conditions contribute to explain G-matrix differentiation
is still rarely explicitly investigated, e.g., by correlating G-matrix
geometry to environmental variables. Changes in environmental
conditions are expected to have a greater impact on the selection
regime for life history traits than for morphological traits (e.g.,
Lundblad & Conway, 2020 in birds) because life history traits are
directly affected by the environment. We thus predict that vari-
ation among the G-matrix for life history traits should be more
strongly correlated to environmental variation than variation
among G-matrices for morphological traits (Arnold et al., 2008).

Linking variation in G-matrices among populations with envi-
ronmental variation requires to define relevant environmental
factors describing global differences in selective pressures across
a species distribution range. A meta-analysis across plants and
animals highlighted the relation between climatic conditions
(temperature and precipitation) and spatial variation in selection
coefficients (Siepielski et al., 2017, 2018). The set of climatic condi-
tions where a species' population can have a positive growth rate
 termed climatic niche hereafter) thus appears as an appropriate
integrative description of the environmental space that might
shape the selection regime of traits and the among-population var-
iation in G-matrices. In line with this, evolvability of single traits
in bird populations tends to be higher for intermediate levels of
climatic favorability (Martínez-Padilla et al., 2017). More generally,
the reduction in effective sample size due to harsh environmental
conditions, such as low-climatic favorability experienced by pop-
ulations at the niche border, could result in a global reduction of
additive genetic variance (Polechová & Barton, 2015) and has been
advocated as the primary mechanism for explaining species distri-
bution boundaries (Holt et al., 2008). Thanks to the development of
long-term monitoring of wild populations (Culina et al., 2021) with
replicated populations for some model species in ecology, it is now
possible to compare G-matrices across populations ranging from
the center to the edge of the species' climatic niche.

Estimation of the G-matrix in wild populations requires datasets
with measures of phenotypic traits and associated pedigree over
several decades. In this study, we took advantage of 10 exceptional
long-term datasets on great tits (Parus major) populations moni-
tored in the wild across Europe from Spain to Finland and Russia
to estimate and compare G-matrices among contrasted environ-
mental conditions. The overall genetic differentiation (Fst) between
European populations of great tits has been shown to be low (Laine
et al., 2016; Lemoine et al., 2016; Spurgin et al., 2019), suggesting a
limited impact of drift and founder effects on the variation among
G-matrices compared to environmental drivers. G-matrices were
inferred separately for morphological and life history traits in each
population using "animal models" under a Bayesian framework.
Based on climatic data available at large spatial scales and high
Introduction

Along with phenotypic plasticity and range shifts toward suit-
able environments, evolution represents a key mechanism for
natural populations to persist under consistently adverse envi-
ronmental conditions (Chevin et al., 2010). The rate and direc-
tion of evolutionary response critically depend on both natural
selection and available standing genetic variance (Walsh & Lynch,
2018). The joint response to selection of multiple traits can be
facilitated or constrained by the genetic association among traits
(Agrawal & Stinchcombe, 2009; Bolstad et al., 2014; Schluter,
1996; Stinchcombe & Kirkpatrick, 2012; Teplitsky et al., 2014).
Under this framework, the matrix of additive genetic variances
and covariances (G-matrix), which summarizes the multivariate
standing genetic variance, plays a fundamental role in predict-
ing the response to selection as the adaptive response across a
single generation is given by the product of the G-matrix and the
selection gradient (Lande, 1979). The extrapolation of the adap-
tive response to selection over multiple generations or across
space thus largely depends on the constancy of the G-matrix
in time and space. The stability of the G-matrix is a paramount
assumption, prevalent in quantitative genetic models (Duputié et
al., 2012).

The structure of G-matrices results from the interactive effects
of selection, mutation, drift, and migration processes (Steppan et
al., 2002). Since these processes, which jointly affect allele fre-
quencies across populations, can change both in space and time,
theory predicts that G-matrices should also vary (Chantepie &
Chevin, 2020). Spatiotemporal changes in the G-matrix can be
best interpreted in geometrical terms (Jones et al., 2004) such as
their variation in size (overall evolutionary potential), orientation
 dir ect of fastest evolution), and shape (constraint resulting
from genetic correlation). In finite populations, the shape of the
G-matrix is a compromise between the shape and orientation of
the matrix of selective gradients and the matrix of pleiotropic
mutations (Jones et al., 2007, 2014; Revell, 2007). At mutation-
selection equilibrium, when drift is high, the G-matrix tends to
orient toward the direction of the mutational matrix (Chantepie &
Chevin, 2020). Differences in these forces among populations will
contribute to differentiation among G-matrices, while gene flow
is expected to homogenize allele frequencies across populations
hence reducing the spatial variation in G-matrices (Guillaume &
Whitlock, 2007; Lande, 1992). Finally, drift may also lead to ran-
dom divergence of G-matrices among populations (Phillips et
al., 2001). Experimental studies conducted over the last decades
have demonstrated that each of the above-mentioned evolution-
ary forces can significantly affect the geometry of the G-matrix
(examples for drift: Phillips et al., 2001, selection: Cano et al.,
2004; Careau et al., 2015; Doroszuk et al., 2008; Shaw et al., 1995,
dispersal: Nosil et al., 2006). In addition to evolutionary forces,
gene expression depending on the environment (G × E) can gen-
erate variation in the G-matrix independently of changes in allele
frequencies (Wood & Brodie, 2015).

Although experimental studies are valuable for understand-
ing how evolutionary forces can shape the G-matrix and testing
theoretical hypotheses, the extrapolation of laboratory results to
wild populations living under natural conditions remains chal-
lenging (Kruuk et al., 2008). In particular, the difficulty to meas-
ure precisely each of these evolutionary forces as well as their
interactions makes it impossible to predict theoretically the var-
iation in G-matrices in the wild (Arnold et al., 2008; Roff, 2000).
Investigating this variation among populations can provide infer-
ence on forces driving spatial differences but also, using space

for time substitution, it can contribute to improve predictions of
temporal changes in G-matrices. Addressing the spatial variation
in the G-matrix requires thus not only to characterize the extent
to which the G-matrix varies among populations but also the
spatial determinants (e.g., geographic distance, environmental
conditions) that have shaped this variation. So far, the literature
relative to G-matrix variation among populations has provided
contrasting results with evidence for both stability (Ashman,
2003; Brodie, 1993; Delahaie et al., 2017; Hangartner et al., 2020;
Henry & Stinchcombe, 2023; Puentes et al., 2016; Sniegula et
al., 2018) and variation (Berger et al., 2013; Eroukhmanoff &
Svensson, 2011; Paccard et al., 2016; Roff et al., 2004; Teplitsky et
al., 2011) of the G-matrix across various spatial gradients in plant
and animal populations. The extent to which differing environ-
mental conditions contribute to explain G-matrix differentiation
is still rarely explicitly investigated, e.g., by correlating G-matrix
geometry to environmental variables. Changes in environmental
conditions are expected to have a greater impact on the selection
regime for life history traits than for morphological traits (e.g.,
Lundblad & Conway, 2020 in birds) because life history traits are
directly affected by the environment. We thus predict that vari-
ation among the G-matrix for life history traits should be more
strongly correlated to environmental variation than variation
among G-matrices for morphological traits (Arnold et al., 2008).

Linking variation in G-matrices among populations with envi-
ronmental variation requires to define relevant environmental
factors describing global differences in selective pressures across
a species distribution range. A meta-analysis across plants and
animals highlighted the relation between climatic conditions
(temperature and precipitation) and spatial variation in selection
coefficients (Siepielski et al., 2017, 2018). The set of climatic condi-
tions where a species' population can have a positive growth rate
 termed climatic niche hereafter) thus appears as an appropriate
integrative description of the environmental space that might
shape the selection regime of traits and the among-population var-
iation in G-matrices. In line with this, evolvability of single traits
in bird populations tends to be higher for intermediate levels of
climatic favorability (Martínez-Padilla et al., 2017). More generally,
the reduction in effective sample size due to harsh environmental
conditions, such as low-climatic favorability experienced by pop-
ulations at the niche border, could result in a global reduction of
additive genetic variance (Polechová & Barton, 2015) and has been
advocated as the primary mechanism for explaining species distri-
bution boundaries (Holt et al., 2008). Thanks to the development of
long-term monitoring of wild populations (Culina et al., 2021) with
replicated populations for some model species in ecology, it is now
possible to compare G-matrices across populations ranging from
the center to the edge of the species' climatic niche.

Estimation of the G-matrix in wild populations requires datasets
with measures of phenotypic traits and associated pedigree over
several decades. In this study, we took advantage of 10 exceptional
long-term datasets on great tits (Parus major) populations moni-
tored in the wild across Europe from Spain to Finland and Russia
to estimate and compare G-matrices among contrasted environ-
mental conditions. The overall genetic differentiation (Fst) between
European populations of great tits has been shown to be low (Laine
et al., 2016; Lemoine et al., 2016; Spurgin et al., 2019), suggesting a
limited impact of drift and founder effects on the variation among
G-matrices compared to environmental drivers. G-matrices were
inferred separately for morphological and life history traits in each
population using "animal models" under a Bayesian framework.
Based on climatic data available at large spatial scales and high
A

<figure><img src="image_3.png" /></figure>

B

<figure><img src="image_4.png" /></figure>

**Figure 1.** Location of populations within space and climatic niche. (A) Geographical distribution of great tits with the location of the 10 focal populations. (B) Climatic niche of great tits. Results of the principal component analysis realized on the climatic variables (array) used at a resolution of 5 km². Points correspond to the predicted position of the 10 populations within the climatic niche. The gray gradient represents the frequency of the climatic conditions across the distribution and ranges from 1 (lightest gray) to 10⁴ cells (darkest gray). The full names of the temperature variables are available in Methods section.

resolution, we constructed the climatic niche of great tits in their distribution range in order to (a) investigate the association between the relative variation in G-matrices among populations and their distances within the climatic niche for both trait categories and (b) test whether populations closer to the niche center showed more genetic variation than populations at the edge of the niche.

# Methods

## Great tit datasets

The meta-dataset used in this study represents a collection of 10 populations of great tits monitored in the wild at different

locations across Europe (Figure 1). These populations have been monitored for at least 20 years with a maximum of 59 years (Supporting Information 1). In all these populations, birds were banded either as chicks or as breeding parents at their nest with a unique identifying metal ring. This identification was used to build a trans-generational social pedigree in each population (Supporting Information 2, Supplementary Tables S8 and S9).

Focal phenotypic traits consisted of both morphological and life history features. Three morphological measurements were taken on both male and female breeders: tarsus length (in mm), flattened wing length (in mm), and body mass (in g). For the two Belgian populations, mean tarsus length was slightly smaller
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Location of populations within space and climatic niche. (A) Geographical distribution of great tits with the location of the 10 focal populations. (B) Climatic niche of great tits. Results of the principal component analysis realized on the climatic variables (array) used at a resolution of 5 km². Points correspond to the predicted position of the 10 populations within the climatic niche. The gray gradient represents the frequency of the climatic conditions across the distribution and ranges from 1 (lightest gray) to 10⁴ cells (darkest gray). The full names of the temperature variables are available in Methods section.</figcaption></figure>

resolution, we constructed the climatic niche of great tits in their distribution range in order to (a) investigate the association between the relative variation in G-matrices among populations and their distances within the climatic niche for both trait categories and (b) test whether populations closer to the niche center showed more genetic variation than populations at the edge of the niche.

# Methods

## Great tit datasets

The meta-dataset used in this study represents a collection of 10 populations of great tits monitored in the wild at different

locations across Europe (Figure 1). These populations have been monitored for at least 20 years with a maximum of 59 years (Supporting Information 1). In all these populations, birds were banded either as chicks or as breeding parents at their nest with a unique identifying metal ring. This identification was used to build a trans-generational social pedigree in each population (Supporting Information 2, Supplementary Tables S8 and S9).

Focal phenotypic traits consisted of both morphological and life history features. Three morphological measurements were taken on both male and female breeders: tarsus length (in mm), flattened wing length (in mm), and body mass (in g). For the two Belgian populations, mean tarsus length was slightly smaller
than in the other populations due to a different measurement method (pers. comm., Matthysen) (Supporting Information 1, Supplementary Table S1). Similarly, for the Dutch populations, the measurement method led to smaller wing lengths than in other populations. However, additive differences in measures do not impact variance estimates (Hansen et al., 2011). Unlike for morphological traits, only females were considered for the three life history traits analyzed, namely laying date (date of first egg laid in the nest, day 1 = Jan 1st), clutch size (number of eggs), and fledging success (number of chicks leaving the nest). Although great tits can produce a second clutch in some years or may replace a first clutch after a breeding failure, only the first clutch of the year was considered for each female. Fledging success was zero truncated because breeding failure was not recorded similarly in all populations (failures to lay any egg vs. failures following the first egg laid) and can result from many different causes (e.g., climatic event, predation, and experimental effect).

## Estimation of G-matrix

G-matrices were estimated for each population by using multivariate animal models (Kruuk, 2004; Walsh & Lynch, 2018). Distinct animal models were fitted for morphological and life history traits. Multivariate animal models represent a particular form of multivariate mixed models that allow the partitioning of phenotypic variance into its genetic and nongenetic components. The G-matrix was estimated by fitting a random effect based on the coefficient of coancestry among pairwise individuals, inferred from pedigree data (Kruuk, 2004). The nongenetic components in our animal models included a permanent environmental effect (PE-matrix) that accounted for the nonindependence of repeated measurements of the same individual, a year effect (Y-matrix) for nonindependence of values with respect to years and a residual variance (R-matrix). We also fitted the observer ID (Obs-matrix) as a random effect on wing and tarsus length when information was available. For all traits, we fitted age as a fixed effect (as a quadratic continuous function or age classes according to populations). For morphological traits, additional fixed effects for sex, sex by age interaction, and date of measurement (as cubic continuous or time periods) were fitted. The fixed and random parameters were adapted to each dataset to account for differences in data collection (Supporting Information 2, Supplementary Table S10). We dealt with missing fixed quantitative parameters with two strategies using (a) the average data augmentation or (b) the missing data deletion. Although both strategies may bias the parameter estimates (Nakagawa, 2015), no difference was found between the G-matrix estimates depending on the strategy used (not shown). Only the results with the average data augmentation are presented here.

We ran the animal mixed models in a Bayesian framework using the MCMCglmm R-package (Hadfield, 2010). The posterior distribution of each parameter estimates corresponded to $10^3$ independent samples along the Monte Carlo chain. For all animal models, we started with a burn-in phase of $2.10^5$ iterations followed by an estimation phase of $10^6$ iterations where parameter estimates were sampled every $10^3$ steps. The convergence of the models was assessed by visually checking the posterior estimates and ensuring that autocorrelation of all parameter estimates was lower than 0.1. When necessary, both the number of iterations and sample steps were increased to meet the autocorrelation threshold requirement. To facilitate convergence, the animal models were run using slightly informative inverse Wishart prior distributions for each component of the variance. For each population, the scale parameter prior was defined as a diagonal

matrix containing the phenotypic variance of traits divided by the number of random factor in the model on the diagonal (Hadfield, 2010). Finally, the degree of belief for priors was set equal to the number of traits.

## G-matrix standardization

We used the same traits in all populations, but these traits differ widely in terms of scales and units (e.g., gram vs. mm), leading to a risk that larger measurements could disproportionately affect potential differences among G-matrices (Hansen et al., 2011). When the values of a trait are several orders of magnitude higher than in other traits, almost all the variation in G-matrix is explained by this trait. We thus standardized the G-matrices posterior estimates before running comparisons among populations.

Standardization is classically done relatively to the trait overall mean or to the phenotypic variance. However, standardization to the phenotypic mean, considered as a natural measure of evolvability, could not be used here because no natural mean exists for variables expressed on an interval scale such as laying date and may not be appropriate for clutch size (Hansen & Houle, 2008; Pélabon et al., 2020). Furthermore, in natural populations, standardization to the phenotypic variance may lead to strongly biased results due to the correlation between additive genetic variance and environmental variance (Hansen et al., 2011).

Variance ratios of similar traits (e.g., phenotypes and breeding values) are unitless, scale invariant, and biologically interpretable (as the percentage of change between the two trait variances) (Le Maître & Mitteroecker, 2019). In this study, G-matrices were standardized to the interpopulation mean of traits additive genetic variance ($Va$). Within traits, additive genetic variances of all populations are divided by the same constant value (mean traits additive genetic variance), so that the relative variation of these genetic variances among populations is not affected by standardization. Similarly, genetic covariances were divided by the square root of the product of mean traits additive genetic variances to preserve the genetic correlation among traits. Overall, this standardization is focusing on preserving the relative variation in variance and covariance of traits among populations (see Supporting Information 3 for details) in order to adequately describe the variations among G-matrices but was not intended to give insights on the difference in absolute evolvability among traits. To estimate $Va$, we first drew a thousand multivariate normal estimates for each posterior sample of all G-matrix, and then calculated the variance of all merged estimates. The interpopulation mean additive genetic variance for each trait can be interpreted as the trait-specific genetic variance of species, assuming that the populations are panmictic. A variance ratio above one thus means that the population harbors more genetic variance than average for this trait. We provided unstandardized G-matrix, P-matrix, and heritability for all traits in Supporting Information 4 (Supplementary Tables S11-S17).

## Null model

To ensure that our results do not arise solely from sampling effects and uncertainty in the estimation of G-matrices, a null model was developed following the general framework of Morrissey et al. (2019). For each population, the breeding values of founders were drawn from a multivariate distribution with zero mean and a variance-covariance matrix equal to the G-matrix estimated for this population with the animal model. These breeding values were then randomly shuffled among the founders from the different populations. This random attribution of breeding values simulates a unique G-matrix for all populations, from which
than in the other populations due to a different measurement method (pers. comm., Matthysen) (Supporting Information 1, Supplementary Table S1). Similarly, for the Dutch populations, the measurement method led to smaller wing lengths than in other populations. However, additive differences in measures do not impact variance estimates (Hansen et al., 2011). Unlike for morphological traits, only females were considered for the three life history traits analyzed, namely laying date (date of first egg laid in the nest, day 1 = Jan 1st), clutch size (number of eggs), and fledging success (number of chicks leaving the nest). Although great tits can produce a second clutch in some years or may replace a first clutch after a breeding failure, only the first clutch of the year was considered for each female. Fledging success was zero truncated because breeding failure was not recorded similarly in all populations (failures to lay any egg vs. failures following the first egg laid) and can result from many different causes (e.g., climatic event, predation, and experimental effect).

## Estimation of G-matrix

G-matrices were estimated for each population by using multivariate animal models (Kruuk, 2004; Walsh & Lynch, 2018). Distinct animal models were fitted for morphological and life history traits. Multivariate animal models represent a particular form of multivariate mixed models that allow the partitioning of phenotypic variance into its genetic and nongenetic components. The G-matrix was estimated by fitting a random effect based on the coefficient of coancestry among pairwise individuals, inferred from pedigree data (Kruuk, 2004). The nongenetic components in our animal models included a permanent environmental effect (PE-matrix) that accounted for the nonindependence of repeated measurements of the same individual, a year effect (Y-matrix) for nonindependence of values with respect to years and a residual variance (R-matrix). We also fitted the observer ID (Obs-matrix) as a random effect on wing and tarsus length when information was available. For all traits, we fitted age as a fixed effect (as a quadratic continuous function or age classes according to populations). For morphological traits, additional fixed effects for sex, sex by age interaction, and date of measurement (as cubic continuous or time periods) were fitted. The fixed and random parameters were adapted to each dataset to account for differences in data collection (Supporting Information 2, Supplementary Table S10). We dealt with missing fixed quantitative parameters with two strategies using (a) the average data augmentation or (b) the missing data deletion. Although both strategies may bias the parameter estimates (Nakagawa, 2015), no difference was found between the G-matrix estimates depending on the strategy used (not shown). Only the results with the average data augmentation are presented here.

We ran the animal mixed models in a Bayesian framework using the MCMCglmm R-package (Hadfield, 2010). The posterior distribution of each parameter estimates corresponded to $10^3$ independent samples along the Monte Carlo chain. For all animal models, we started with a burn-in phase of $2.10^5$ iterations followed by an estimation phase of $10^6$ iterations where parameter estimates were sampled every $10^3$ steps. The convergence of the models was assessed by visually checking the posterior estimates and ensuring that autocorrelation of all parameter estimates was lower than 0.1. When necessary, both the number of iterations and sample steps were increased to meet the autocorrelation threshold requirement. To facilitate convergence, the animal models were run using slightly informative inverse Wishart prior distributions for each component of the variance. For each population, the scale parameter prior was defined as a diagonal

matrix containing the phenotypic variance of traits divided by the number of random factor in the model on the diagonal (Hadfield, 2010). Finally, the degree of belief for priors was set equal to the number of traits.

## G-matrix standardization

We used the same traits in all populations, but these traits differ widely in terms of scales and units (e.g., gram vs. mm), leading to a risk that larger measurements could disproportionately affect potential differences among G-matrices (Hansen et al., 2011). When the values of a trait are several orders of magnitude higher than in other traits, almost all the variation in G-matrix is explained by this trait. We thus standardized the G-matrices posterior estimates before running comparisons among populations.

Standardization is classically done relatively to the trait overall mean or to the phenotypic variance. However, standardization to the phenotypic mean, considered as a natural measure of evolvability, could not be used here because no natural mean exists for variables expressed on an interval scale such as laying date and may not be appropriate for clutch size (Hansen & Houle, 2008; Pélabon et al., 2020). Furthermore, in natural populations, standardization to the phenotypic variance may lead to strongly biased results due to the correlation between additive genetic variance and environmental variance (Hansen et al., 2011).

Variance ratios of similar traits (e.g., phenotypes and breeding values) are unitless, scale invariant, and biologically interpretable (as the percentage of change between the two trait variances) (Le Maître & Mitteroecker, 2019). In this study, G-matrices were standardized to the interpopulation mean of traits additive genetic variance ($Va$). Within traits, additive genetic variances of all populations are divided by the same constant value (mean traits additive genetic variance), so that the relative variation of these genetic variances among populations is not affected by standardization. Similarly, genetic covariances were divided by the square root of the product of mean traits additive genetic variances to preserve the genetic correlation among traits. Overall, this standardization is focusing on preserving the relative variation in variance and covariance of traits among populations (see Supporting Information 3 for details) in order to adequately describe the variations among G-matrices but was not intended to give insights on the difference in absolute evolvability among traits. To estimate $Va$, we first drew a thousand multivariate normal estimates for each posterior sample of all G-matrix, and then calculated the variance of all merged estimates. The interpopulation mean additive genetic variance for each trait can be interpreted as the trait-specific genetic variance of species, assuming that the populations are panmictic. A variance ratio above one thus means that the population harbors more genetic variance than average for this trait. We provided unstandardized G-matrix, P-matrix, and heritability for all traits in Supporting Information 4 (Supplementary Tables S11-S17).

## Null model

To ensure that our results do not arise solely from sampling effects and uncertainty in the estimation of G-matrices, a null model was developed following the general framework of Morrissey et al. (2019). For each population, the breeding values of founders were drawn from a multivariate distribution with zero mean and a variance-covariance matrix equal to the G-matrix estimated for this population with the animal model. These breeding values were then randomly shuffled among the founders from the different populations. This random attribution of breeding values simulates a unique G-matrix for all populations, from which
population genetic values are sampled. Based on these random breeding values for founders and each population's pedigree, we estimated the breeding values of all remaining individuals following the infinitesimal assumption. The breeding values of an individual were drawn from a multivariate normal distribution with a vector of means equal to the mean breeding values of its parents and a deviation equal to the segregation variance (Walsh & Lynch, 2018). For an individual, the segregation variance was estimated as half the product of the randomized G-matrix of their population (estimated on the new population's founders) and one minus the average inbreeding coefficient of its parents (Walsh & Lynch, 2018). For each population, the phenotypes were then recomposed by adding to the simulated breeding values the random and fixed effects estimated by the animal model run on the wild dataset. The fixed effects were estimated as the product of the fixed-effect design matrix and the vector of fixed effects drawn from the posterior sample of the animal model for each population. Random effects were drawn from a multivariate distribution with a mean of zero and their respective variance-covariance matrix resulting from the animal models (i.e., PE-matrix, Y-matrix, Obs-matrix, and R-matrix). Original populations estimates for random effects other than additive genetic variance were kept because the amount of heritability can affect the power of the analysis (Morrissey et al., 2007) and effects such as year effect or permanent environment are environmental and thus population specific.

This sequence of generating phenotypes was repeated for each posterior sample of animal models and resulted in 1,000 null datasets for each population. The G-matrix of null datasets was estimated with the same animal models used for the corresponding wild population dataset. To reduce computational burden, we sampled 200 posterior samples for each null animal model (Morrissey et al., 2019) to compute a mean null distribution G-matrix for each population.

## Climatic data and niche modeling

The climatic niche of great tits was estimated using known species occurrences and climatic data available at a large geographical scale. We used the great tit distribution map from Birdlife International and Nature Serve (Birdlife International and Nature Serve, 2018), which includes all subspecies. In this studies, we only considered European subspecies of the great tit and restricted geographic distribution to the West Palearctic region (Kvist et al., 1999). We extracted temperature data from WorldClim2 Global Climate GIS data at a 2.5 min resolution (~5 km) over the species distribution (Fick & Hijmans, 2017). The WorldClim2 contained climate data that were monthly averaged for the 1970-2000 period, a scale relevant for our study as we were interested in large-scale climate variation (rather than weather anomalies) and since founders of the populations (i.e., individuals whose parent is unknown) could enter the datasets any year during the monitoring (either as immigrants or chicks born outside of nest boxes). The temperature variables included to build the climatic niche were: the annual mean temperature, mean diurnal temperature range, temperature seasonality, maximum temperature of warmest month, minimum temperature of coldest month, temperature annual range, mean temperature of warmest quarter, and mean temperature of coldest quarter (all temperatures in degrees Celsius).

To estimate the niche space, a principal component analysis (PCA) was performed on these eight climatic variables over all cells included in the distribution of great tits (Figure 1A). We assessed the position of each population within the bioclimatic

niche by (a) extracting the local WorldClim2 climatic conditions at their respective geographical coordinates and (b) projecting these local conditions on the two first PCA axes (Figure 1B). The niche distances among populations were estimated as the Euclidean distance among the projected population positions within the first two axes PCA space. As the niche distance and geographic distance were highly correlated, only the former was used for analyses. We considered the origin of the PCA space as the center of the bioclimatic niche and estimated the Euclidean distance between the center of the niche and the position of the projected populations.

## Statistical differences in G-matrices among populations

### Genetic covariance tensor analysis

Differences in genetic covariance structure across populations were tested using genetic covariance tensors (Hine et al., 2009). Instead of performing pairwise comparisons between matrices, the tensor method investigates the changes in the structure of G-matrices along a gradient. The tensor approach describes the variation among G-matrices based on successive linear transformations and identifies independent combinations of traits displaying a change in genetic variance. Detailed descriptions of the method can be found in Aguirre et al. (2014), Walter et al. (2018), and Supporting Information 5.

Succinctly, a genetic covariance tensor was built from the variances and covariances between the elements of the G-matrices. An eigen-decomposition of this covariance tensor was performed to assess the independent axes of variation among the G-matrices. The eigenvectors of this decomposition, called eigenvectors, represent the directions which maximize the variance among the original G-matrices. The magnitude of the variation among the G-matrices along a particular eigentensor was determined by the corresponding eigenvalues. The coordinates of each G-matrix in the space defined by the eigentensors describe how much the populations differ along a given eigentensor. In other words, when two G-matrices are close with respect to a particular tensor, then these matrices are described in the same way by the eigentensor. The eigentensor itself can be eigen-decomposed to eigenvectors that assess which independent linear combinations of traits have contributed to the change in G-matrix structure along this tensor. The corresponding eigenvalue of the eigenvectors describes the magnitude of the contribution of trait in these changes. We applied this method to the posterior estimates of the animal models (Aguirre et al., 2014) and compared the results to our null model.

The distances among populations were estimated as the Euclidean distance among coordinates in the space defined by the two first tensors (hereafter "tensor distance"). To evaluate the extent to which climatic conditions are driving the differentiation among G-matrices, first, we tested the relation between the tensor distances estimated with posterior mode of coordinates and niche distances using a Mantel test. Second, to account for the uncertainty in coordinates estimates on the tensors, we also estimated 1,000 tensor distances matrix using the posterior distribution of coordinates. To remove the dependence among the tensor distances and allow for estimating correlation with niche distance, we shuffled each element among the 1,000 distance matrices. Finally, we performed linear regressions with randomized tensor distance as the dependent variable and the niche distance as the explanatory variable for all the tensor distance matrices and tested significance of the slope using the 95% credible interval (CI).
population genetic values are sampled. Based on these random breeding values for founders and each population's pedigree, we estimated the breeding values of all remaining individuals following the infinitesimal assumption. The breeding values of an individual were drawn from a multivariate normal distribution with a vector of means equal to the mean breeding values of its parents and a deviation equal to the segregation variance (Walsh & Lynch, 2018). For an individual, the segregation variance was estimated as half the product of the randomized G-matrix of their population (estimated on the new population's founders) and one minus the average inbreeding coefficient of its parents (Walsh & Lynch, 2018). For each population, the phenotypes were then recomposed by adding to the simulated breeding values the random and fixed effects estimated by the animal model run on the wild dataset. The fixed effects were estimated as the product of the fixed-effect design matrix and the vector of fixed effects drawn from the posterior sample of the animal model for each population. Random effects were drawn from a multivariate distribution with a mean of zero and their respective variance-covariance matrix resulting from the animal models (i.e., PE-matrix, Y-matrix, Obs-matrix, and R-matrix). Original populations estimates for random effects other than additive genetic variance were kept because the amount of heritability can affect the power of the analysis (Morrissey et al., 2007) and effects such as year effect or permanent environment are environmental and thus population specific.

This sequence of generating phenotypes was repeated for each posterior sample of animal models and resulted in 1,000 null datasets for each population. The G-matrix of null datasets was estimated with the same animal models used for the corresponding wild population dataset. To reduce computational burden, we sampled 200 posterior samples for each null animal model (Morrissey et al., 2019) to compute a mean null distribution G-matrix for each population.

## Climatic data and niche modeling

The climatic niche of great tits was estimated using known species occurrences and climatic data available at a large geographical scale. We used the great tit distribution map from Birdlife International and Nature Serve (Birdlife International and Nature Serve, 2018), which includes all subspecies. In this studies, we only considered European subspecies of the great tit and restricted geographic distribution to the West Palearctic region (Kvist et al., 1999). We extracted temperature data from WorldClim2 Global Climate GIS data at a 2.5 min resolution (~5 km) over the species distribution (Fick & Hijmans, 2017). The WorldClim2 contained climate data that were monthly averaged for the 1970-2000 period, a scale relevant for our study as we were interested in large-scale climate variation (rather than weather anomalies) and since founders of the populations (i.e., individuals whose parent is unknown) could enter the datasets any year during the monitoring (either as immigrants or chicks born outside of nest boxes). The temperature variables included to build the climatic niche were: the annual mean temperature, mean diurnal temperature range, temperature seasonality, maximum temperature of warmest month, minimum temperature of coldest month, temperature annual range, mean temperature of warmest quarter, and mean temperature of coldest quarter (all temperatures in degrees Celsius).

To estimate the niche space, a principal component analysis (PCA) was performed on these eight climatic variables over all cells included in the distribution of great tits (Figure 1A). We assessed the position of each population within the bioclimatic

niche by (a) extracting the local WorldClim2 climatic conditions at their respective geographical coordinates and (b) projecting these local conditions on the two first PCA axes (Figure 1B). The niche distances among populations were estimated as the Euclidean distance among the projected population positions within the first two axes PCA space. As the niche distance and geographic distance were highly correlated, only the former was used for analyses. We considered the origin of the PCA space as the center of the bioclimatic niche and estimated the Euclidean distance between the center of the niche and the position of the projected populations.

## Statistical differences in G-matrices among populations

### Genetic covariance tensor analysis

Differences in genetic covariance structure across populations were tested using genetic covariance tensors (Hine et al., 2009). Instead of performing pairwise comparisons between matrices, the tensor method investigates the changes in the structure of G-matrices along a gradient. The tensor approach describes the variation among G-matrices based on successive linear transformations and identifies independent combinations of traits displaying a change in genetic variance. Detailed descriptions of the method can be found in Aguirre et al. (2014), Walter et al. (2018), and Supporting Information 5.

Succinctly, a genetic covariance tensor was built from the variances and covariances between the elements of the G-matrices. An eigen-decomposition of this covariance tensor was performed to assess the independent axes of variation among the G-matrices. The eigenvectors of this decomposition, called eigenvectors, represent the directions which maximize the variance among the original G-matrices. The magnitude of the variation among the G-matrices along a particular eigentensor was determined by the corresponding eigenvalues. The coordinates of each G-matrix in the space defined by the eigentensors describe how much the populations differ along a given eigentensor. In other words, when two G-matrices are close with respect to a particular tensor, then these matrices are described in the same way by the eigentensor. The eigentensor itself can be eigen-decomposed to eigenvectors that assess which independent linear combinations of traits have contributed to the change in G-matrix structure along this tensor. The corresponding eigenvalue of the eigenvectors describes the magnitude of the contribution of trait in these changes. We applied this method to the posterior estimates of the animal models (Aguirre et al., 2014) and compared the results to our null model.

The distances among populations were estimated as the Euclidean distance among coordinates in the space defined by the two first tensors (hereafter "tensor distance"). To evaluate the extent to which climatic conditions are driving the differentiation among G-matrices, first, we tested the relation between the tensor distances estimated with posterior mode of coordinates and niche distances using a Mantel test. Second, to account for the uncertainty in coordinates estimates on the tensors, we also estimated 1,000 tensor distances matrix using the posterior distribution of coordinates. To remove the dependence among the tensor distances and allow for estimating correlation with niche distance, we shuffled each element among the 1,000 distance matrices. Finally, we performed linear regressions with randomized tensor distance as the dependent variable and the niche distance as the explanatory variable for all the tensor distance matrices and tested significance of the slope using the 95% credible interval (CI).
## Orientation, shape, and size of G-matrices

We assessed the overall structure of each G-matrix through eigenanalyses. The major eigenvector of a G-matrix called $g_{max}$ is oriented in the direction of maximum genetic variance in the population in the standardized space of traits. The ratio between the eigenvalue of $g_{max}$ and the sum of all the eigenvalues provides information about the eccentricity of the G-matrix, or how much the response to selection is favored in the direction of $g_{max}$.

The angle between two populations refers to the variation in their direction of fastest evolution in the standardized space of traits. We estimated the critical angles between $g_{max}$ of two G-matrices using the method of Krzanowski (1979). We examined the statistical support for two populations having different $g_{max}$ by comparing the intrapopulation variation of $g_{max}$ (drawn from the posterior distribution) with the interpopulation variation $g_{max}$ (Ovaskainen et al., 2008) (Supporting Information 6). To evaluate whether the direction of fastest evolution diverged according to the difference in climatic conditions, the angles between $g_{max}$ were compared against niche distance using Mantel test and linear regressions (see details above).

To test for the hypothesis of a decrease in the total relative genetic variance or total relative size from the center to the edge of the niche, we used the volume of each G-matrix as a proxy of the total genetic variance. For each G-matrix, we estimated the posterior distribution of volume as the sum of the eigenvalues for all independent posterior samples of the G-matrices. Then, we performed a linear regression modeling the posterior volume distribution as a function of the distance of the population to the center of the niche (significance of the slope tested using the 95% CI).

## Results

### Climatic niche and phenotype divergence across populations

The 10 populations of great tits used in this study were distributed over a wide geographic area (distance between populations ranged from 2 to 3,300 km) and were found in contrasting environments (Figure 1B). The first two axes of the PCA accounted for 61% and 33% of the variance in bioclimatic data, respectively. The first axis of the PCA explained broad variations in temperature within the distribution range of great tits and the second axis mostly described environmental variability (Figure 1B).

Populations living in the most contrasted temperature environments showed significant variation in phenotypic traits. For example, individuals from the Valencia population (Spain) fledged on average 3.5 less offspring than those from the Zvenigorod population (Russia) (Supporting Information 1). The mean laying date in the Oulu population (Finland) occurred 1 month later than in the Mediterranean populations. Regarding body mass, the only morphological trait which is not influenced by differences in measurement methods (see Methods for details), we found that Spanish individuals were more than 1.7 g lighter than their Russian counterparts. All statistical tests on mean differences are provided in Supporting Information 1 (Supplementary Tables S2–S7).

### Variation in life history G-matrix linked with climatic conditions

Given our standardization strategy, results highlighted the relative change in G-matrix across the 10 populations studied. For life history traits, the distance among standardized G-matrices

in tensors space can be summarized by the two first eigentensors accounting altogether for 64% of variances among G-matrices (different from the null model at *p*-value of .03 and .09, respectively, Supporting Information 5, Supplementary Figure S5). The distances among G-matrices in tensors space is positively linked with the distances among populations within the climatic niche (Figure 2C, slope\_data[95% CI] = 0.13[0.06:0.21]; Mantel test: *p*-value= .001, *r* = 0.69). This differentiation among G-matrices according to climatic differences was not detected under the null model (slope\_null[95% CI] = 0.02[−0.02:0.09]; Mantel test: *p*-value = .38, *r* = 0.05, slope\_data is statistically higher than slope\_null at a *p*-value of .016). This result was robust to the removal of the two most extreme populations (Supporting Information 5, Supplementary Figures S8–S10). The first tensor accounted for 48% [26%:71%] of the variation among G-matrices and revealed that the G-matrices of populations from Southern (Spain, south of France) and Northern Europe (Finland, Russia) tended to be different from the other populations. Along this tensor, the major difference in G-matrices were led by a coordinated change in genetic variance for clutch size and fledging success (the first eigenvector of the first tensor explained 97% of the variation with loadings of −0.74 and −0.64 for clutch size and fledging success, respectively). Significant differences in coordinates on the first tensor were found between the Oulu population (Finland) and the populations located around 50°N (Belgium, Netherlands, England) (Figure 2A, Supporting Information 5, Supplementary Figure S7) reflecting the significantly lower additive genetic variance in clutch size and lower genetic correlation between clutch size and fledging success in the Finnish population (Figure 3A). Overall, the genetic correlation between clutch size and fledging success was lower for populations experiencing higher environmental variability (PC2 in Figures 1 and 3B).

The second tensor accounted for 16% [5%:41%] of the variation among standardized G-matrices and mostly represented a change in additive genetic variance of laying date (the first eigenvector of the second tensor explained 81% of the variation with loadings of 0.93 for the laying date). This tensor highlighted differences between Valencia (Spain) and Zvenigorod (Russia): the additive genetic variance of laying date was significantly higher in the former (Figures 2A and 3A, Supplementary Figure S7). This difference resulted in a significant relation between the genetic variance of laying date and the distance to the niche center (Figure 3A; slope[95% CI] = 0.29[0.05:0.56]) mostly driven by the temperature gradient (Figure 3B; slope[95% CI] = 0.11[0.02:0.23]).

The angles among populations' $g_{max}$ also increased with the niche distance between populations (Figure 4A) (slope\_data[95% CI] = 8.79[5.40:10.71]; Mantel *p*-value = .01, *r* = 0.70; Supporting Information 6). These differences arose because, $g_{max}$ was oriented in the direction of laying date for both the Spanish and Finnish populations, of fledging success for the Russian population, and of both clutch size and fledging success in the other populations. As the standardized genetic variances are of the same magnitude for both clutch size and fledging success of these populations, the direction of $g_{max}$ is mostly defined by the genetic correlation between these two traits (Figure 4B, bidimensional ellipses are provided in Supporting Information 7). The variances explained by $g_{max}$ were similar among all populations (between 0.55 and 0.68, Supporting Information 7, Supplementary Figures S14 and S15). The relation between angles among populations' $g_{max}$ and the niche distance estimated on the null model was not significant (slope\_null[95% CI] = 1.05[−1.52:4.3]; Mantel test: *p*-value = .10, *r* = 0.38, slope\_data is statistically higher than slope\_null at a *p*-value of .001).
## Orientation, shape, and size of G-matrices

We assessed the overall structure of each G-matrix through eigenanalyses. The major eigenvector of a G-matrix called $g_{max}$ is oriented in the direction of maximum genetic variance in the population in the standardized space of traits. The ratio between the eigenvalue of $g_{max}$ and the sum of all the eigenvalues provides information about the eccentricity of the G-matrix, or how much the response to selection is favored in the direction of $g_{max}$.

The angle between two populations refers to the variation in their direction of fastest evolution in the standardized space of traits. We estimated the critical angles between $g_{max}$ of two G-matrices using the method of Krzanowski (1979). We examined the statistical support for two populations having different $g_{max}$ by comparing the intrapopulation variation of $g_{max}$ (drawn from the posterior distribution) with the interpopulation variation $g_{max}$ (Ovaskainen et al., 2008) (Supporting Information 6). To evaluate whether the direction of fastest evolution diverged according to the difference in climatic conditions, the angles between $g_{max}$ were compared against niche distance using Mantel test and linear regressions (see details above).

To test for the hypothesis of a decrease in the total relative genetic variance or total relative size from the center to the edge of the niche, we used the volume of each G-matrix as a proxy of the total genetic variance. For each G-matrix, we estimated the posterior distribution of volume as the sum of the eigenvalues for all independent posterior samples of the G-matrices. Then, we performed a linear regression modeling the posterior volume distribution as a function of the distance of the population to the center of the niche (significance of the slope tested using the 95% CI).

## Results

### Climatic niche and phenotype divergence across populations

The 10 populations of great tits used in this study were distributed over a wide geographic area (distance between populations ranged from 2 to 3,300 km) and were found in contrasting environments (Figure 1B). The first two axes of the PCA accounted for 61% and 33% of the variance in bioclimatic data, respectively. The first axis of the PCA explained broad variations in temperature within the distribution range of great tits and the second axis mostly described environmental variability (Figure 1B).

Populations living in the most contrasted temperature environments showed significant variation in phenotypic traits. For example, individuals from the Valencia population (Spain) fledged on average 3.5 less offspring than those from the Zvenigorod population (Russia) (Supporting Information 1). The mean laying date in the Oulu population (Finland) occurred 1 month later than in the Mediterranean populations. Regarding body mass, the only morphological trait which is not influenced by differences in measurement methods (see Methods for details), we found that Spanish individuals were more than 1.7 g lighter than their Russian counterparts. All statistical tests on mean differences are provided in Supporting Information 1 (Supplementary Tables S2–S7).

### Variation in life history G-matrix linked with climatic conditions

Given our standardization strategy, results highlighted the relative change in G-matrix across the 10 populations studied. For life history traits, the distance among standardized G-matrices

in tensors space can be summarized by the two first eigentensors accounting altogether for 64% of variances among G-matrices (different from the null model at *p*-value of .03 and .09, respectively, Supporting Information 5, Supplementary Figure S5). The distances among G-matrices in tensors space is positively linked with the distances among populations within the climatic niche (Figure 2C, slope\_data[95% CI] = 0.13[0.06:0.21]; Mantel test: *p*-value= .001, *r* = 0.69). This differentiation among G-matrices according to climatic differences was not detected under the null model (slope\_null[95% CI] = 0.02[-0.02:0.09]; Mantel test: *p*-value = .38, *r* = 0.05, slope\_data is statistically higher than slope\_null at a *p*-value of .016). This result was robust to the removal of the two most extreme populations (Supporting Information 5, Supplementary Figures S8–S10). The first tensor accounted for 48% [26%:71%] of the variation among G-matrices and revealed that the G-matrices of populations from Southern (Spain, south of France) and Northern Europe (Finland, Russia) tended to be different from the other populations. Along this tensor, the major difference in G-matrices were led by a coordinated change in genetic variance for clutch size and fledging success (the first eigenvector of the first tensor explained 97% of the variation with loadings of -0.74 and -0.64 for clutch size and fledging success, respectively). Significant differences in coordinates on the first tensor were found between the Oulu population (Finland) and the populations located around 50°N (Belgium, Netherlands, England) (Figure 2A, Supporting Information 5, Supplementary Figure S7) reflecting the significantly lower additive genetic variance in clutch size and lower genetic correlation between clutch size and fledging success in the Finnish population (Figure 3A). Overall, the genetic correlation between clutch size and fledging success was lower for populations experiencing higher environmental variability (PC2 in Figures 1 and 3B).

The second tensor accounted for 16% [5%:41%] of the variation among standardized G-matrices and mostly represented a change in additive genetic variance of laying date (the first eigenvector of the second tensor explained 81% of the variation with loadings of 0.93 for the laying date). This tensor highlighted differences between Valencia (Spain) and Zvenigorod (Russia): the additive genetic variance of laying date was significantly higher in the former (Figures 2A and 3A, Supplementary Figure S7). This difference resulted in a significant relation between the genetic variance of laying date and the distance to the niche center (Figure 3A; slope[95% CI] = 0.29[0.05:0.56]) mostly driven by the temperature gradient (Figure 3B; slope[95% CI] = 0.11[0.02:0.23]).

The angles among populations' $g_{max}$ also increased with the niche distance between populations (Figure 4A) (slope\_data[95% CI] = 8.79[5.40:10.71]; Mantel *p*-value = .01, *r* = 0.70; Supporting Information 6). These differences arose because, $g_{max}$ was oriented in the direction of laying date for both the Spanish and Finnish populations, of fledging success for the Russian population, and of both clutch size and fledging success in the other populations. As the standardized genetic variances are of the same magnitude for both clutch size and fledging success of these populations, the direction of $g_{max}$ is mostly defined by the genetic correlation between these two traits (Figure 4B, bidimensional ellipses are provided in Supporting Information 7). The variances explained by $g_{max}$ were similar among all populations (between 0.55 and 0.68, Supporting Information 7, Supplementary Figures S14 and S15). The relation between angles among populations' $g_{max}$ and the niche distance estimated on the null model was not significant (slope\_null[95% CI] = 1.05[-1.52:4.3]; Mantel test: *p*-value = .10, *r* = 0.38, slope\_data is statistically higher than slope\_null at a *p*-value of .001).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Genetic covariance tensor results. Points in the upper panels represent the posterior modes of the coordinates of each matrix along the two first tensors for (A) life history G-matrices and (B) morphological G-matrices. Colors correspond to the different populations (see Figure 1 for details) and gray bars to the 95%, 75%, and 50% credible intervals (CI). Relation between the Euclidean distance among populations in the two first tensor spaces and the distance of population within their climatic niche are represented for (C) life history and (D) morphological G-matrices. Points represent the distances between pairs of populations (posterior mode for the tensor distances). Colored hemispheres correspond to the population considered for pairwise comparison. The light blue and gray envelops represent the linear regressions realized on posterior distributions of distance estimated using data and null model, respectively. Only regressions with a slope within 95% CI are shown.</figcaption></figure>

Although major changes occurred among G-matrices of life history traits, no decrease in their total volume was found from the center to the edge of the niche (slope[95% CI] = 0.31[-0.12:0.28], Figure 3C).

## Variation of morphological G-matrix

For morphological traits, no significant linear relation was found between the distance among standardized G-matrices in the space defined by the two first tensors and the distance among populations within the climatic niche (Figure 2D, slope[95% CI] = -0.008[-0.07:0.04], not significantly different from the null model). Distance among G-matrices in tensors space was based on the two first eigentensors accounting for 61% of variation among G-matrices (different from the null model at *p*-value of .19 and .04, respectively, Supporting Information 5, Supplementary Figure S6).

The first tensor accounted for 37% [20%:65%] of the variation among G-matrices. This axis segregated Peerdsbos (Belgium) and Zvenigorod (Russia) from Oosterhout (Netherland) populations (Figure 2B). The differences among the G-matrices of these populations were related to a coordinated change in the genetic variance of wing length, tarsus length, and body mass (the first eigenvector of the first tensor explained 98% of the variation

with loadings of both traits, respectively, -0.66, -0.41, and -0.63). The genetic correlation between wing length and body mass was significantly higher in Peerdsbos population than in Oulu, Oosterhout, and Boshoeck populations (Figure 5).

The second tensor (24% [11%:47%]) mostly represented a difference in G-matrices between Peerdsbos (Belgium) and several other populations (Valencia, Rouvière, Oosterhout, Zvenigorod, and Hoge Veluwe, Supplementary Figure S7) due to a low additive genetic variance of tarsus length in this population (76% of the variation along the first eigenvector of the second eigentensor; loadings of tarsus length 0.91; Figure 5).

No relation between the angles among populations' $g_{max}$ and the niche distance was found (slope[95% CI] = -0.24[-2.70:1.19], Supplementary Figure S13) as well as between the total size and the distance to the niche center (slope[95% CI] = 0.02[-0.28:0.26], Supplementary Figure S16).

## Discussion

No general consensus can be reached regarding the spatial variation of G-matrices in the wild without evaluating how this variation is environmentally driven. In particular, an important question that needs to be addressed is whether a G-matrix
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Genetic covariance tensor results. Points in the upper panels represent the posterior modes of the coordinates of each matrix along the two first tensors for (A) life history G-matrices and (B) morphological G-matrices. Colors correspond to the different populations (see Figure 1 for details) and gray bars to the 95%, 75%, and 50% credible intervals (CI). Relation between the Euclidean distance among populations in the two first tensor spaces and the distance of population within their climatic niche are represented for (C) life history and (D) morphological G-matrices. Points represent the distances between pairs of populations (posterior mode for the tensor distances). Colored hemispheres correspond to the population considered for pairwise comparison. The light blue and gray envelops represent the linear regressions realized on posterior distributions of distance estimated using data and null model, respectively. Only regressions with a slope within 95% CI are shown.</figcaption></figure>

Although major changes occurred among G-matrices of life history traits, no decrease in their total volume was found from the center to the edge of the niche (slope[95% CI] = 0.31[-0.12:0.28], Figure 3C).

## Variation of morphological G-matrix

For morphological traits, no significant linear relation was found between the distance among standardized G-matrices in the space defined by the two first tensors and the distance among populations within the climatic niche (Figure 2D, slope[95% CI] = -0.008[-0.07:0.04], not significantly different from the null model). Distance among G-matrices in tensors space was based on the two first eigentensors accounting for 61% of variation among G-matrices (different from the null model at *p*-value of .19 and .04, respectively, Supporting Information 5, Supplementary Figure S6).

The first tensor accounted for 37% [20%:65%] of the variation among G-matrices. This axis segregated Peerdsbos (Belgium) and Zvenigorod (Russia) from Oosterhout (Netherland) populations (Figure 2B). The differences among the G-matrices of these populations were related to a coordinated change in the genetic variance of wing length, tarsus length, and body mass (the first eigenvector of the first tensor explained 98% of the variation

with loadings of both traits, respectively, -0.66, -0.41, and -0.63). The genetic correlation between wing length and body mass was significantly higher in Peerdsbos population than in Oulu, Oosterhout, and Boshoeck populations (Figure 5).

The second tensor (24% [11%:47%]) mostly represented a difference in G-matrices between Peerdsbos (Belgium) and several other populations (Valencia, Rouvière, Oosterhout, Zvenigorod, and Hoge Veluwe, Supplementary Figure S7) due to a low additive genetic variance of tarsus length in this population (76% of the variation along the first eigenvector of the second eigentensor; loadings of tarsus length 0.91; Figure 5).

No relation between the angles among populations' $g_{max}$ and the niche distance was found (slope[95% CI] = -0.24[-2.70:1.19], Supplementary Figure S13) as well as between the total size and the distance to the niche center (slope[95% CI] = 0.02[-0.28:0.26], Supplementary Figure S16).

## Discussion

No general consensus can be reached regarding the spatial variation of G-matrices in the wild without evaluating how this variation is environmentally driven. In particular, an important question that needs to be addressed is whether a G-matrix
<figure><img src="image_6.png" /><figcaption>**Figure 3.** (A) Standardized additive genetic variances and genetic correlations estimated for life history traits. (B) Genetic variance of the laying date and genetic correlation between clutch size and fledging success related with the principal components of the environmental niche (see Supplementary Information 4 for all traits). (C) Volume of the life history G-matrix organized along the distance to the niche center. Points represent the posterior mode estimates filled with colors corresponding to their respective populations and gray bars the 95%, 75%, and 50% credible intervals (CI) (see Figures 1 and 2 for details). On the laying date panel, the gray envelope represents the 95% CI of the significant linear relationship between genetic variance and distance to niche center. Significant differences at a *p*-value of .05 (**) are displayed on top of figures for panel (A).</figcaption></figure>

estimated in one population can be extrapolated to other populations or environments (current or future). Using long-term phenotypic data associated with pedigrees in 10 wild populations of great tits, we explored intraspecific variation in G-matrices across the species range. We provide compelling evidence for differences in G-matrices among populations of great tits across Europe for both life history and morphological traits. Interestingly, variation in G-matrices for life history traits (but not morphological traits) could be related to the population's position along a climatic niche gradient. The major divergence in genetic architecture for life history traits was measured between the populations located close to the 50°N latitude (Belgium, Netherlands, and England) and the other populations located in the Mediterranean and Subarctic climates.

## Variation in the life history G-matrix linked with the bioclimatic niche

### G-matrix variation among populations

We detected a significant link between the distance between two populations within the species' climatic niche and their difference in G-matrices for life history traits.

Among all potential mechanisms that may have shaped this relative variation in G-matrices among populations in relation with climate, a response to a change in selection regimes due to environmental variation is a possible explanation. Specifically, we

showed that the variation in the G-matrix observed among these 10 wild populations was correlated with the difference in environment for the life history traits, traits most likely to be affected by the change in selection regime due to environmental variation (e.g., Lundblad & Conway, 2020). The large impact of a change in selection regime on the shape of the G-matrix was previously shown in controlled experiments performed on a wide range species (e.g., Cano et al., 2004; Wood & Brodie, 2015). Both the strength and direction of selection vary in space among wild populations of plants and animals especially among contrasted environmental conditions (Siepielski et al., 2017): beyond the question of how G-matrices vary spatially, the open challenge is to identify the nature and strength of environmental variations that alter G-matrix structure such as volume and orientation. In this regard, nonsignificant variation in G-matrices between contrasted environments found in some empirical studies (e.g., reviewed in Arnold et al., 2008) might reveal a lack of contrast between environmental conditions. For example, the lack of variation among life history and morphological G-matrices estimated in island and mainland Mediterranean populations of blue tits (Delahaie et al., 2017) may simply reflect the limited climatic differences that exist between these four Mediterranean populations. On the contrary, significant variation in G-matrices was found between a Mediterranean (Spain) and a Continental (Denmark) populations of barn swallows (Hirundo rustica) (Teplitsky et al., 2011).
<figure><img src="image_6.png" /><figcaption>**Figure 3.** (A) Standardized additive genetic variances and genetic correlations estimated for life history traits. (B) Genetic variance of the laying date and genetic correlation between clutch size and fledging success related with the principal components of the environmental niche (see Supplementary Information 4 for all traits). (C) Volume of the life history G-matrix organized along the distance to the niche center. Points represent the posterior mode estimates filled with colors corresponding to their respective populations and gray bars the 95%, 75%, and 50% credible intervals (CI) (see Figures 1 and 2 for details). On the laying date panel, the gray envelope represents the 95% CI of the significant linear relationship between genetic variance and distance to niche center. Significant differences at a *p*-value of .05 (**) are displayed on top of figures for panel (A).</figcaption></figure>

estimated in one population can be extrapolated to other populations or environments (current or future). Using long-term phenotypic data associated with pedigrees in 10 wild populations of great tits, we explored intraspecific variation in G-matrices across the species range. We provide compelling evidence for differences in G-matrices among populations of great tits across Europe for both life history and morphological traits. Interestingly, variation in G-matrices for life history traits (but not morphological traits) could be related to the population's position along a climatic niche gradient. The major divergence in genetic architecture for life history traits was measured between the populations located close to the 50°N latitude (Belgium, Netherlands, and England) and the other populations located in the Mediterranean and Subarctic climates.

## Variation in the life history G-matrix linked with the bioclimatic niche

### G-matrix variation among populations

We detected a significant link between the distance between two populations within the species' climatic niche and their difference in G-matrices for life history traits.

Among all potential mechanisms that may have shaped this relative variation in G-matrices among populations in relation with climate, a response to a change in selection regimes due to environmental variation is a possible explanation. Specifically, we

showed that the variation in the G-matrix observed among these 10 wild populations was correlated with the difference in environment for the life history traits, traits most likely to be affected by the change in selection regime due to environmental variation (e.g., Lundblad & Conway, 2020). The large impact of a change in selection regime on the shape of the G-matrix was previously shown in controlled experiments performed on a wide range species (e.g., Cano et al., 2004; Wood & Brodie, 2015). Both the strength and direction of selection vary in space among wild populations of plants and animals especially among contrasted environmental conditions (Siepielski et al., 2017): beyond the question of how G-matrices vary spatially, the open challenge is to identify the nature and strength of environmental variations that alter G-matrix structure such as volume and orientation. In this regard, nonsignificant variation in G-matrices between contrasted environments found in some empirical studies (e.g., reviewed in Arnold et al., 2008) might reveal a lack of contrast between environmental conditions. For example, the lack of variation among life history and morphological G-matrices estimated in island and mainland Mediterranean populations of blue tits (Delahaie et al., 2017) may simply reflect the limited climatic differences that exist between these four Mediterranean populations. On the contrary, significant variation in G-matrices was found between a Mediterranean (Spain) and a Continental (Denmark) populations of barn swallows (Hirundo rustica) (Teplitsky et al., 2011).
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Difference in orientation between G-matrices of life history traits. (A) Relationship between the difference in $g_{max}$ direction for each pair of populations and their respective niche distance. Points represent the distances between pairs of populations (posterior mode for the angle distances). Colored hemispheres refer to the population considered for pairwise comparison. The light blue and gray envelops represent the linear regressions realized on posterior distributions of distance estimated using data and null model, respectively. Only regressions with a slope within 95% CI are shown. (B) Projection of the posterior mode of $g_{max}$ axes for each population into the additive genetic variance trait space with colors corresponding to their respective populations (see Figure 1 for details).</figcaption></figure>

Altogether and despite a lack of generality, these studies and our results tend to show that contrasted climatic conditions, or more generally broad environmental condition, represent a significant factor that partly shape the G-matrix of life history traits. Hence, while G-matrices are not stable, using environment for time substitution could allow to predict more accurate responses to selection, for example, in future warmer climate.

While the spatial variation of G-matrices can partly result from differential selection due to particular local climatic conditions, several alternative scenarios may also explain our results. First, genotype-environment interactions (G × E) can produce correlation patterns between climatic conditions and G-matrix variation. Recent studies that compared DNA markers among populations of great tits distributed across Europe highlighted low overall genetic differentiation among European populations of great tits except for populations in small Mediterranean islands (i.e., Corsica, Sardinia, and Crete) (Laine et al., 2016; Lemoine et al., 2016; Spurgin et al., 2019). In this context, it is likely that a G × E interaction may contribute to the variation of G-matrices and produce a pattern that may be similar to a change in allelic

frequencies between populations. Disentangling both processes requires common garden experiments, a difficult task for many species, but future studies based on gene expression may help solve this issue. Whether G-matrices differences result from G × E and/or selection processes, this study highlights the importance of climatic conditions on the structure of the G-matrix and consequences for predicting responses to selection in these environments.

Second, genetic drift and founder effects may lead to large G-matrix variation among wild populations (e.g., Paccard et al., 2016). The low overall genetic differentiation between populations at the European scale and generally large effective population sizes suggest that drift may not be a determining factor for most of the differences observed between our G-matrix (Kvist et al., 1999; Laine et al., 2016). Nevertheless, significant genetic differentiation has been reported between Spanish populations and other populations, suggesting that this population may be more isolated and experienced higher selective pressure or stronger drift (Lemoine et al., 2016). Third, indirect genetic effects (i.e., parental care, social behavior) may alter the distribution of
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Difference in orientation between G-matrices of life history traits. (A) Relationship between the difference in $g_{max}$ direction for each pair of populations and their respective niche distance. Points represent the distances between pairs of populations (posterior mode for the angle distances). Colored hemispheres refer to the population considered for pairwise comparison. The light blue and gray envelops represent the linear regressions realized on posterior distributions of distance estimated using data and null model, respectively. Only regressions with a slope within 95% CI are shown. (B) Projection of the posterior mode of $g_{max}$ axes for each population into the additive genetic variance trait space with colors corresponding to their respective populations (see Figure 1 for details).</figcaption></figure>

Altogether and despite a lack of generality, these studies and our results tend to show that contrasted climatic conditions, or more generally broad environmental condition, represent a significant factor that partly shape the G-matrix of life history traits. Hence, while G-matrices are not stable, using environment for time substitution could allow to predict more accurate responses to selection, for example, in future warmer climate.

While the spatial variation of G-matrices can partly result from differential selection due to particular local climatic conditions, several alternative scenarios may also explain our results. First, genotype-environment interactions (G × E) can produce correlation patterns between climatic conditions and G-matrix variation. Recent studies that compared DNA markers among populations of great tits distributed across Europe highlighted low overall genetic differentiation among European populations of great tits except for populations in small Mediterranean islands (i.e., Corsica, Sardinia, and Crete) (Laine et al., 2016; Lemoine et al., 2016; Spurgin et al., 2019). In this context, it is likely that a G × E interaction may contribute to the variation of G-matrices and produce a pattern that may be similar to a change in allelic

frequencies between populations. Disentangling both processes requires common garden experiments, a difficult task for many species, but future studies based on gene expression may help solve this issue. Whether G-matrices differences result from G × E and/or selection processes, this study highlights the importance of climatic conditions on the structure of the G-matrix and consequences for predicting responses to selection in these environments.

Second, genetic drift and founder effects may lead to large G-matrix variation among wild populations (e.g., Paccard et al., 2016). The low overall genetic differentiation between populations at the European scale and generally large effective population sizes suggest that drift may not be a determining factor for most of the differences observed between our G-matrix (Kvist et al., 1999; Laine et al., 2016). Nevertheless, significant genetic differentiation has been reported between Spanish populations and other populations, suggesting that this population may be more isolated and experienced higher selective pressure or stronger drift (Lemoine et al., 2016). Third, indirect genetic effects (i.e., parental care, social behavior) may alter the distribution of
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Standardized additive genetic variances and genetic correlations estimated for morphological traits. Details of the legend are provided in Figure 3.</figcaption></figure>

phenotypic variance within and among families and affect the estimation of additive genetic variance when using animal models (Kruuk & Hadfield, 2007). Consequently, our results might also partly reflect a potential variation across Europe of some indirect genetic effects in great tits such as parental effect (Garant et al., 2008; Postma, 2005). Unfortunately, a precise partitioning of the additive genetic and indirect genetic effects was not possible for all the populations due to lack of sufficient data.

## G-matrix variation along niche gradient

Contrary to expectations, no coordinated decrease in genetic variance was found toward the edge of the environmental niche, suggesting a limited impact of genetic drift in the edge populations for this species. Edge effects on evolutionary potential have been previously found in European wild birds when comparing trait evolvabilities along climatic gradients for 12 species (Martínez-Padilla et al., 2017). Our result is however in line with the overall low genetic differentiation among European populations of great tits populations in this study suggesting low-drift effect and that none of the populations here is truly at the edge of the climatic niche. Actually, we even found a significant increase in the additive genetic variance of laying date along a temperature gradient, with the highest genetic variance in the Spanish population, closest to the edge of the great tits niche. This pattern of increase in genetic variance from the center to the niche margin has been predicted to arise when gene flow between populations living in contrasted environments with different phenotypic optima (Polechová, 2018; Polechová & Barton, 2015). It could be the case here as laying date is generally considered to display an optimum in great tits (Gamelon et al., 2018; De Villemereuil et al., 2020). Under this theoretical model, genetic variance is expected to decrease at the extreme margin of the niche only, where the effects of drift and demographic stochasticity become predominant.

We did not detect a consistent change of all (co)-variance components along the niche gradient, and this should not be expected unless selection vary in a coordinated way across the climatic

range for all traits. Moreover, the interpretation of variations of each element of the G-matrix may be misleading when the genetic correlation between traits varies greatly from one population to another (Blows, 2007). Here, we found strong genetic correlations between clutch size and fledging success only in populations located around 50°N latitude, likely defined by an environmental variability gradient. Thus, a finer understanding of change in genetic variance and covariance of the clutch size and fledging success will require to compare the phenotypic optima and selection strength, including correlational selection, across populations.

## Variation in the G-matrix for morphological traits

While the G-matrix for morphological traits differed among populations, no link with a climatic gradient was found, as opposed to our findings for the G-matrix for life history traits. According to the tensor analysis, the standardized morphological G-matrix was most different between the Russian populations and one of the Belgian and the Dutch populations (Figure 2B) due to both a change in genetic variances and genetic correlations between traits. These results suggested that factors other than climatic conditions were responsible for this variation in morphological genetic (co)variances. Overall, the stability of the genetic architecture was higher for morphological traits than life history traits (Supplementary Figures S14 and S15).

Several nonmutually exclusive explanations may be involved in the relative stability of the morphological G-matrix compared to the life history G-matrix. First, morphological traits have often been shown to be under stabilizing selection in birds (e.g., mass: Covas et al., 2002; tarsus length: Alatalo & Lundberg, 1986) and the constraint imposed by flight traits is high and possibly shared across populations independently of weather conditions. The relative stability of the G-matrix of morphological traits in our study was in line with predictions based on the expected stability of selection pressures on these traits (Arnold et al., 2001). Second, the existence of unmeasured traits strongly correlated with focal morphological traits can constrain and even prevent
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Standardized additive genetic variances and genetic correlations estimated for morphological traits. Details of the legend are provided in Figure 3.</figcaption></figure>

phenotypic variance within and among families and affect the estimation of additive genetic variance when using animal models (Kruuk & Hadfield, 2007). Consequently, our results might also partly reflect a potential variation across Europe of some indirect genetic effects in great tits such as parental effect (Garant et al., 2008; Postma, 2005). Unfortunately, a precise partitioning of the additive genetic and indirect genetic effects was not possible for all the populations due to lack of sufficient data.

## G-matrix variation along niche gradient

Contrary to expectations, no coordinated decrease in genetic variance was found toward the edge of the environmental niche, suggesting a limited impact of genetic drift in the edge populations for this species. Edge effects on evolutionary potential have been previously found in European wild birds when comparing trait evolvabilities along climatic gradients for 12 species (Martínez-Padilla et al., 2017). Our result is however in line with the overall low genetic differentiation among European populations of great tits populations in this study suggesting low-drift effect and that none of the populations here is truly at the edge of the climatic niche. Actually, we even found a significant increase in the additive genetic variance of laying date along a temperature gradient, with the highest genetic variance in the Spanish population, closest to the edge of the great tits niche. This pattern of increase in genetic variance from the center to the niche margin has been predicted to arise when gene flow between populations living in contrasted environments with different phenotypic optima (Polechová, 2018; Polechová & Barton, 2015). It could be the case here as laying date is generally considered to display an optimum in great tits (Gamelon et al., 2018; De Villemereuil et al., 2020). Under this theoretical model, genetic variance is expected to decrease at the extreme margin of the niche only, where the effects of drift and demographic stochasticity become predominant.

We did not detect a consistent change of all (co)-variance components along the niche gradient, and this should not be expected unless selection vary in a coordinated way across the climatic

range for all traits. Moreover, the interpretation of variations of each element of the G-matrix may be misleading when the genetic correlation between traits varies greatly from one population to another (Blows, 2007). Here, we found strong genetic correlations between clutch size and fledging success only in populations located around 50°N latitude, likely defined by an environmental variability gradient. Thus, a finer understanding of change in genetic variance and covariance of the clutch size and fledging success will require to compare the phenotypic optima and selection strength, including correlational selection, across populations.

## Variation in the G-matrix for morphological traits

While the G-matrix for morphological traits differed among populations, no link with a climatic gradient was found, as opposed to our findings for the G-matrix for life history traits. According to the tensor analysis, the standardized morphological G-matrix was most different between the Russian populations and one of the Belgian and the Dutch populations (Figure 2B) due to both a change in genetic variances and genetic correlations between traits. These results suggested that factors other than climatic conditions were responsible for this variation in morphological genetic (co)variances. Overall, the stability of the genetic architecture was higher for morphological traits than life history traits (Supplementary Figures S14 and S15).

Several nonmutually exclusive explanations may be involved in the relative stability of the morphological G-matrix compared to the life history G-matrix. First, morphological traits have often been shown to be under stabilizing selection in birds (e.g., mass: Covas et al., 2002; tarsus length: Alatalo & Lundberg, 1986) and the constraint imposed by flight traits is high and possibly shared across populations independently of weather conditions. The relative stability of the G-matrix of morphological traits in our study was in line with predictions based on the expected stability of selection pressures on these traits (Arnold et al., 2001). Second, the existence of unmeasured traits strongly correlated with focal morphological traits can constrain and even prevent
structural changes in G-matrix. Indeed, strong pleiotropic genes or selective correlation between traits could constrain and stabilize the G-matrix (Chantepie & Chevin, 2020; Jones et al., 2003). A further understanding of variation in the G-matrix of morphological traits would require the comparison of populations experimenting substantially different selection regimes on morphology. Ecological specialization has, for example, been shown to shape the G-matrix for morphological traits across seven Anolis lizards species (McGlothlin et al., 2018). A recent study showed that great tits that live in cities are smaller than forest birds (Caizergues et al., 2021), which could be due to a change in the selection regime and lead to a variation of the G-matrix. Selection due to migration might also act on morphological traits of passerines (Vágási et al., 2016). Although Northern populations have been shown to be partial residents migrating irregularly (Kvist et al., 1999), the lack of longitudinal data with pedigree and migration information for these populations preclude actual testing of the impact of migration on the G-matrix of morphological traits. However, undertaking comparisons of G-matrices in these ecological specialization contexts would help to test our hypothesis that G-matrix of morphological traits is only relatively insensitive to variation in climatic conditions but not to all environmental variations.

## Conclusion

The increasing number of long-term datasets collected including both phenotypic and pedigree information in wild populations opens new perspectives to assess the link between evolutionary processes and environmental conditions. This study shows that assuming a constant G-matrix among populations may not be a valid approximation for life history traits and to a lesser extent for morphological traits. Furthermore, the heterogeneity of climatic conditions across the spatial distribution of great tits partially explained the differences in G-matrices among populations for life history traits which open new perspectives in understanding and predicting the adaptive response of life history traits to global changes at the species level.

## Supplementary material

Supplementary material is available online at Evolution Letters.

## Data and code availability

Data are available upon request at https://spibirds.org. This hub acts as a central contact point for accessing the data, ensures quality checks on the data to correct for any mistakes in the datasets and warns users about potential issues when using the data. This hub is maintained by the Netherlands Institute of Ecology. Full output of animal models, unstandardized G-matrix, standardized G-matrix, unstandardized G-matrix estimated with the null models, and standardized G-matrix estimated with the null models are available at https://doi.org/10.6084/m9.fishshare.24783900.v1. Code for null model simulation is available at https://github.com/schantepie/NullGmatrix/.

## Author contributions

S.C. and C.T. designed the study. S.C. performed the analyses and wrote the manuscript with major contributions from C.T., S.L., and A.C. B.D., F.A., E.M., E.Á., E.B., M.O., B.S., E.I., and A.K. contributed

to the interpretation of the results, provided critical feedback, and significantly edited and approved the final manuscript.

## Funding

A.C. and C.T. thank the long-standing funding support of Observatoire des Sciences de l'Univers OREME (Observatoire de REcherche Montpelliérain de l'Environnement) since 2009 for the monitoring in Rouvière. E.B. and E.A. thank the financial support by the Spanish Ministry of Science and Innovation (PID2021-122171NB-I00/AEI/10.13039/501100011033/FEDER, UE). C.T. is funded by the CNRS and by a grant from the BNP-Paribas Foundation ("Climate and Biodiversity initiative," project Mommy Knows Best).

*Conflict of interest:* The authors declare no conflict of interest. Editorial processing of the manuscript was done independently of A.C. and C.T., who is an Associate and Guest Editor, respectively, of Evolution Letters.

## Acknowledgments

We wish to thank Thomas Hansen, David Houle, and Christophe Pélabon for discussions on standardization and preliminary results. We would like to thank the two anonymous reviewers for their detailed and meaningful comments, which greatly contributed to the improvement of the manuscript.

## References

Agrawal, A. F., & Stinchcombe, J. R. (2009). How much do genetic covariances alter the rate of adaptation? *Proceedings of the Royal Society B: Biological Sciences*, 276(1659), 1183–1191. https://doi.org/10.1098/rspb.2008.1671

Aguirre, J. D., Hine, E., McGuigan, K., & Blows, M. W. (2014). Comparing G: Multivariate analysis of genetic variation in multiple populations. *Heredity*, 112(1), 21–29. https://doi.org/10.1038/hdy.2013.12

Alatalo, R. V., & Lundberg, A. (1986). Heritability and selection on taxes length in the pied flycatcher (*Ficedula hypoleuca*). *Evolution*, 40(3), 574–583. https://doi.org/10.1111/j.1558-5646.1986.tb00508.x

Arnold, S. J., Bürger, R., Hohenlohe, P. A., Ajie, B. C., & Jones, A. G. (2008). Understanding the evolution and stability of the G-matrix. *Evolution*, 62(10), 2451–2461. https://doi.org/10.1111/j.1558-5646.2008.00472.x

Arnold, S. J., Pfrender, M. E., & Jones, A. G. (2001). The adaptive landscape as a conceptual bridge between micro- and macroevolution. *Genetica*, 112–113, 9–32.

Ashman, T. -L. (2003). Constraints on the evolution of males and sexual dimorphism: Field estimates of genetic architecture of reproductive traits in three populations of Gynodioecious *Fragaria virginiana*. *Evolution*, 57(9), 2012–2025. https://doi.org/10.1111/j.0014-3820.2003.tb00381.x

Berger, D., Postma, E., Blanckenhorn, W. U., & Walters, R. J. (2013). Quantitative genetic divergence and standing genetic (co)variability in thermal reaction norms along latitude. *Evolution*, 67(8), 2385–2399. https://doi.org/10.1111/evo.12138

Birdlife International and Nature Serve. (2018). Bird species distribution maps of the world [WWW Document]. Version 7.0. http://datazone.birdlife.org/species/requestdis [accessed 2018].

Blows, M. W. (2007). A tale of two matrices: Multivariate approaches in evolutionary biology. *Journal of Evolutionary Biology*, 20(1), 1–8. https://doi.org/10.1111/j.1420-9101.2006.01164.x
structural changes in G-matrix. Indeed, strong pleiotropic genes or selective correlation between traits could constrain and stabilize the G-matrix (Chantepie & Chevin, 2020; Jones et al., 2003). A further understanding of variation in the G-matrix of morphological traits would require the comparison of populations experimenting substantially different selection regimes on morphology. Ecological specialization has, for example, been shown to shape the G-matrix for morphological traits across seven Anolis lizards species (McGlothlin et al., 2018). A recent study showed that great tits that live in cities are smaller than forest birds (Caizergues et al., 2021), which could be due to a change in the selection regime and lead to a variation of the G-matrix. Selection due to migration might also act on morphological traits of passerines (Vágási et al., 2016). Although Northern populations have been shown to be partial residents migrating irregularly (Kvist et al., 1999), the lack of longitudinal data with pedigree and migration information for these populations preclude actual testing of the impact of migration on the G-matrix of morphological traits. However, undertaking comparisons of G-matrices in these ecological specialization contexts would help to test our hypothesis that G-matrix of morphological traits is only relatively insensitive to variation in climatic conditions but not to all environmental variations.

## Conclusion

The increasing number of long-term datasets collected including both phenotypic and pedigree information in wild populations opens new perspectives to assess the link between evolutionary processes and environmental conditions. This study shows that assuming a constant G-matrix among populations may not be a valid approximation for life history traits and to a lesser extent for morphological traits. Furthermore, the heterogeneity of climatic conditions across the spatial distribution of great tits partially explained the differences in G-matrices among populations for life history traits which open new perspectives in understanding and predicting the adaptive response of life history traits to global changes at the species level.

## Supplementary material

Supplementary material is available online at Evolution Letters.

## Data and code availability

Data are available upon request at https://spibirds.org. This hub acts as a central contact point for accessing the data, ensures quality checks on the data to correct for any mistakes in the datasets and warns users about potential issues when using the data. This hub is maintained by the Netherlands Institute of Ecology. Full output of animal models, unstandardized G-matrix, standardized G-matrix, unstandardized G-matrix estimated with the null models, and standardized G-matrix estimated with the null models are available at https://doi.org/10.6084/m9.fishshare.24783900.v1. Code for null model simulation is available at https://github.com/schantepie/NullGmatrix/.

## Author contributions

S.C. and C.T. designed the study. S.C. performed the analyses and wrote the manuscript with major contributions from C.T., S.L., and A.C. B.D., F.A., E.M., E.Á., E.B., M.O., B.S., E.I., and A.K. contributed

to the interpretation of the results, provided critical feedback, and significantly edited and approved the final manuscript.

## Funding

A.C. and C.T. thank the long-standing funding support of Observatoire des Sciences de l'Univers OREME (Observatoire de REcherche Montpelliérain de l'Environnement) since 2009 for the monitoring in Rouvière. E.B. and E.A. thank the financial support by the Spanish Ministry of Science and Innovation (PID2021-122171NB-I00/AEI/10.13039/501100011033/FEDER, UE). C.T. is funded by the CNRS and by a grant from the BNP-Paribas Foundation ("Climate and Biodiversity initiative," project Mommy Knows Best).

*Conflict of interest:* The authors declare no conflict of interest. Editorial processing of the manuscript was done independently of A.C. and C.T., who is an Associate and Guest Editor, respectively, of Evolution Letters.

## Acknowledgments

We wish to thank Thomas Hansen, David Houle, and Christophe Pélabon for discussions on standardization and preliminary results. We would like to thank the two anonymous reviewers for their detailed and meaningful comments, which greatly contributed to the improvement of the manuscript.

## References

Agrawal, A. F., & Stinchcombe, J. R. (2009). How much do genetic covariances alter the rate of adaptation? *Proceedings of the Royal Society B: Biological Sciences*, 276(1659), 1183–1191. https://doi.org/10.1098/rspb.2008.1671

Aguirre, J. D., Hine, E., McGuigan, K., & Blows, M. W. (2014). Comparing G: Multivariate analysis of genetic variation in multiple populations. *Heredity*, 112(1), 21–29. https://doi.org/10.1038/hdy.2013.12

Alatalo, R. V., & Lundberg, A. (1986). Heritability and selection on taxes length in the pied flycatcher (*Ficedula hypoleuca*). *Evolution*, 40(3), 574–583. https://doi.org/10.1111/j.1558-5646.1986.tb00508.x

Arnold, S. J., Bürger, R., Hohenlohe, P. A., Ajie, B. C., & Jones, A. G. (2008). Understanding the evolution and stability of the G-matrix. *Evolution*, 62(10), 2451–2461. https://doi.org/10.1111/j.1558-5646.2008.00472.x

Arnold, S. J., Pfrender, M. E., & Jones, A. G. (2001). The adaptive landscape as a conceptual bridge between micro- and macroevolution. *Genetica*, 112–113, 9–32.

Ashman, T. -L. (2003). Constraints on the evolution of males and sexual dimorphism: Field estimates of genetic architecture of reproductive traits in three populations of Gynodioecious *Fragaria virginiana*. *Evolution*, 57(9), 2012–2025. https://doi.org/10.1111/j.0014-3820.2003.tb00381.x

Berger, D., Postma, E., Blanckenhorn, W. U., & Walters, R. J. (2013). Quantitative genetic divergence and standing genetic (co)variability in thermal reaction norms along latitude. *Evolution*, 67(8), 2385–2399. https://doi.org/10.1111/evo.12138

Birdlife International and Nature Serve. (2018). Bird species distribution maps of the world [WWW Document]. Version 7.0. http://datazone.birdlife.org/species/requestdis [accessed 2018].

Blows, M. W. (2007). A tale of two matrices: Multivariate approaches in evolutionary biology. *Journal of Evolutionary Biology*, 20(1), 1–8. https://doi.org/10.1111/j.1420-9101.2006.01164.x
Bolstad, G. H., Hansen, T. F., Pélabon, C., Falahati-Anbaran, M., Pérez-Barrales, R., & Armbruster, W. S. (2014). Genetic constraints predict evolutionary divergence in Dalechampia blossoms. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 369(1649), 20130255.

Brodie, E. D. (1993). Homogeneity of the genetic variance-covariance matrix for antipredator traits in two natural populations of the garter snake *Thamnophis ordinoides*. *Evolution*, 47(3), 844-854.
https://doi.org/10.1111/j.1558-5646.1993.tb01238.x

Caizergues, A. E., Charmantier, A., Lambrechts, M. M., Perret, S., Demeyrier, V., Lucas, A., & Grégoire, A. (2021). An avian urban morphotype: How the city environment shapes great tit morphology at different life stages. *Urban Ecosystems*, 24(5), 929-941.
https://doi.org/10.1007/s11252-020-01077-0

Cano, J. M., Laurila, A., Palo, J., & Merilä, J. (2004). Population differentiation in G matrix structure due to natural selection in *Rana temporaria*. *Evolution*, 58(9), 2013–2020. https://doi.org/10.1111/j.0014-3820.2004.tb00486.x

Careau, V., Wolak, M. E., Carter, P. A., & Garland, T. (2015). Evolution of the additive genetic variance-covariance matrix under continuous directional selection on a complex behavioural phenotype. *Proceedings of the Royal Society B: Biological Sciences*, 282(1819), 20151119.

Chantepie, S., & Chevin, L. (2020). How does the strength of selection influence genetic correlations? *Evolution Letters*, 4(6), 468–478.
https://doi.org/10.1002/evl3.201

Chevin, L. -M., Lande, R., & Mace, G. M. (2010). Adaptation, plasticity, and extinction in a changing environment: Towards a predictive theory. *PLoS Biology*, 8(4), e1000357. https://doi.org/10.1371/journal.pbio.1000357

Covas, R., Brown, C. R., Anderson, M. D., & Brown, M. B. (2002). Stabilizing selection on body mass in the sociable weaver *Philetairus socius*. *Proceedings of the Royal Society B: Biological Sciences*, 269(1503), 1905–1909. https://doi.org/10.1098/rspb.2002.2106

Culina, A., Adriaensen, F., Bailey, L. D., Burgess, M. D., Charmantier, A., Cole, E. F., Eeva, T., Matthysen, E., Nater, C. R., Sheldon, B. C., Sæther, B. -E., Vriend, S. J. G., Zajkova, Z., Adamík, P., Aplin, L. M., Angulo, E., Artemyev, A., Barba, E., Barišić, S., ... Visser, M. E. (2021). Connecting the data landscape of long-term ecological studies: The SPI-Birds data hub. *Journal of Animal Ecology*, 90(9), 2147–2160. https://doi.org/10.1111/1365-2656.13388

Delahaie, B., Charmantier, A., Chantepie, S., Garant, D., Porlier, M., & Teplitsky, C. (2017). Conserved G-matrices of morphological and life-history traits among continental and island blue tit populations. *Heredity*, 119(2), 76–87. https://doi.org/10.1038/hdy.2017.15

Doroszuk, A., Wojewodzic, M. W., Gort, G., & Kammenga, J. E. (2008). Rapid divergence of genetic variance-covariance matrix within a natural population. *The American Naturalist*, 171(3), 291–304.
https://doi.org/10.1086/527478

Duputié, A., Massol, F., Chuine, I., Kirkpatrick, M., & Ronce, O. (2012). How do genetic correlations affect species range shifts in a changing environment? *Ecology Letters*, 15(3), 251–259. https://doi.org/10.1111/j.1461-0248.2011.01734.x

Eroukhmanoff, F., & Svensson, E. I. (2011). Evolution and stability of the G-matrix during the colonization of a novel environment. *Journal of Evolutionary Biology*, 24(6), 1363–1373. https://doi.org/10.1111/j.1420-9101.2011.02270.x

Fick, S. E., & Hijmans, R. J. (2017). WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. *International Journal of Climatology*, 37(12), 4302–4315.

Gamelon, M., Tufto, J., Nilsson, A. L. K., Jerstad, K., Røstad, O. W., Stenseth, N. C., & Sæther, B. -E. (2018). Environmental drivers of varying selective optima in a small passerine: A multivariate,

multiepisodic approach. *Evolution*, 72(11), 2325-2342. https://doi.org/10.1111/evo.13610

Garant, D., Hadfield, J. D., Kruuk, L. E. B., & Sheldon, B. C. (2008). Stability of genetic variance and covariance for reproductive characters in the face of climate change in a wild bird population. *Molecular Ecology*, 17(1), 179-188. https://doi.org/10.1111/j.1365-294X.2007.03436.x

Guillaume, F., & Whitlock, M. C. (2007). Effects of migration on the genetic covariance matrix. *Evolution*, 61(10), 2398-2409. https://doi.org/10.1111/j.1558-5646.2007.00193.x

Hadfield, J. D. (2010). MCMC methods for multi-response generalized linear mixed models: The MCMCglmm R package. *Journal of Statistical Software*, 33(2), 1-22.

Hangartner, S., Lasne, C., Sgrò, C. M., Connallon, T., & Monro, K. (2020). Genetic covariances promote climatic adaptation in Australian Drosophila*. *Evolution*, 74(2), 326-337. https://doi.org/10.1111/evo.13831

Hansen, T. F., & Houle, D. (2008). Measuring and comparing evolvability and constraint in multivariate characters. *Journal of Evolutionary Biology*, 21(5), 1201-1219. https://doi.org/10.1111/j.1420-9101.2008.01573.x

Hansen, T. F., Pélabon, C., & Houle, D. (2011). Heritability is not evolvability. *Evolutionary Biology*, 38(3), 258-277. https://doi.org/10.1007/s11692-011-9127-6

Henry, G. A., & Stinchcombe, J. R. (2023). G-matrix stability in clinically diverging populations of an annual weed. *Evolution*, 77(1), 49-62. https://doi.org/10.1093/evolut/qpac005

Hine, E., Chenoweth, S. F., Rundle, H. D., & Blows, M. W. (2009). Characterizing the evolution of genetic variance using genetic covariance tensors. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 364(1523), 1567-1578. https://doi.org/10.1098/rstb.2008.0313

Holt, A. R., Davies, Z. G., Tyler, C., & Staddon, S. (2008). Meta-analysis of the effects of predation on animal prey abundance: Evidence from UK vertebrates. *PLoS One*, 3(6), e2400. https://doi.org/10.1371/journal.pone.0002400

Jones, A. G., Arnold, S. J., & Bürger, R. (2003). Stability of the G-matrix in a population experiencing pleiotropic mutation, stabilizing selection, and genetic drift. *Evolution*, 57(8), 1747-1760. https://doi.org/10.1111/j.0014-3820.2003.tb00583.x

Jones, A. G., Arnold, S. J., & Bürger, R. (2004). Evolution and stability of the G-matrix on a landscape with a moving optimum. *Evolution*, 58(8), 1639-1654. https://doi.org/10.1111/j.0014-3820.2004.tb00450.x

Jones, A. G., Arnold, S. J., & Bürger, R. (2007). The mutation matrix and the evolution of evolvability. *Evolution*, 61(4), 727-745. https://doi.org/10.1111/j.1558-5646.2007.00071.x

Jones, A. G., Bürger, R., & Arnold, S. J. (2014). Epistasis and natural selection shape the mutational architecture of complex traits. *Nature Communications*, 5, 3709. https://doi.org/10.1038/ncomms4709

Kruuk, L. E. B. (2004). Estimating genetic parameters in natural populations using the “Animal Model”. *Proceedings of the Royal Society B: Biological Sciences*, 359, 873-890.

Kruuk, L. E. B., & Hadfield, J. D. (2007). How to separate genetic and environmental causes of similarity between relatives. *Journal of Evolutionary Biology*, 20(5), 1890-1903. https://doi.org/10.1111/j.1420-9101.2007.01377.x

Kruuk, L. E. B., Slate, J., & Wilson, A. J. (2008). New answers for old questions: The evolutionary quantitative genetics of wild animal populations. *Annual Review of Ecology, Evolution, and Systematics*, 39, 525-548.

Krzanowski, W. J. (1979). Between-groups comparison of principal components. *Journal of the American Statistical Association*, 74, 703-707.
Bolstad, G. H., Hansen, T. F., Pélabon, C., Falahati-Anbaran, M., Pérez-Barrales, R., & Armbruster, W. S. (2014). Genetic constraints predict evolutionary divergence in Dalechampia blossoms. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 369(1649), 20130255.

Brodie, E. D. (1993). Homogeneity of the genetic variance-covariance matrix for antipredator traits in two natural populations of the garter snake *Thamnophis ordinoides*. *Evolution*, 47(3), 844–854.
https://doi.org/10.1111/j.1558-5646.1993.tb01238.x

Caizergues, A. E., Charmantier, A., Lambrechts, M. M., Perret, S., Demeyrier, V., Lucas, A., & Grégoire, A. (2021). An avian urban morphotype: How the city environment shapes great tit morphology at different life stages. *Urban Ecosystems*, 24(5), 929–941.
https://doi.org/10.1007/s11252-020-01077-0

Cano, J. M., Laurila, A., Palo, J., & Merilä, J. (2004). Population differentiation in G matrix structure due to natural selection in *Rana temporaria*. *Evolution*, 58(9), 2013–2020. https://doi.org/10.1111/j.0014-3820.2004.tb00486.x

Careau, V., Wolak, M. E., Carter, P. A., & Garland, T. (2015). Evolution of the additive genetic variance-covariance matrix under continuous directional selection on a complex behavioural phenotype. *Proceedings of the Royal Society B: Biological Sciences*, 282(1819), 20151119.

Chantepie, S., & Chevin, L. (2020). How does the strength of selection influence genetic correlations? *Evolution Letters*, 4(6), 468–478.
https://doi.org/10.1002/evl3.201

Chevin, L. -M., Lande, R., & Mace, G. M. (2010). Adaptation, plasticity, and extinction in a changing environment: Towards a predictive theory. *PLoS Biology*, 8(4), e1000357. https://doi.org/10.1371/journal.pbio.1000357

Covas, R., Brown, C. R., Anderson, M. D., & Brown, M. B. (2002). Stabilizing selection on body mass in the sociable weaver *Philetairus socius*. *Proceedings of the Royal Society B: Biological Sciences*, 269(1503), 1905–1909. https://doi.org/10.1098/rspb.2002.2106

Culina, A., Adriaensen, F., Bailey, L. D., Burgess, M. D., Charmantier, A., Cole, E. F., Eeva, T., Matthysen, E., Nater, C. R., Sheldon, B. C., Sæther, B. -E., Vriend, S. J. G., Zajkova, Z., Adamík, P., Aplin, L. M., Angulo, E., Artemyev, A., Barba, E., Barišić, S., ... Visser, M. E. (2021). Connecting the data landscape of long-term ecological studies: The SPI-Birds data hub. *Journal of Animal Ecology*, 90(9), 2147–2160. https://doi.org/10.1111/1365-2656.13388

Delahaie, B., Charmantier, A., Chantepie, S., Garant, D., Porlier, M., & Teplitsky, C. (2017). Conserved G-matrices of morphological and life-history traits among continental and island blue tit populations. *Heredity*, 119(2), 76–87. https://doi.org/10.1038/hdy.2017.15

Doroszuk, A., Wojewodzic, M. W., Gort, G., & Kammenga, J. E. (2008). Rapid divergence of genetic variance-covariance matrix within a natural population. *The American Naturalist*, 171(3), 291–304.
https://doi.org/10.1086/527478

Duputié, A., Massol, F., Chuine, I., Kirkpatrick, M., & Ronce, O. (2012). How do genetic correlations affect species range shifts in a changing environment? *Ecology Letters*, 15(3), 251–259. https://doi.org/10.1111/j.1461-0248.2011.01734.x

Eroukhmanoff, F., & Svensson, E. I. (2011). Evolution and stability of the G-matrix during the colonization of a novel environment. *Journal of Evolutionary Biology*, 24(6), 1363–1373. https://doi.org/10.1111/j.1420-9101.2011.02270.x

Fick, S. E., & Hijmans, R. J. (2017). WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. *International Journal of Climatology*, 37(12), 4302–4315.

Gamelon, M., Tufto, J., Nilsson, A. L. K., Jerstad, K., Røstad, O. W., Stenseth, N. C., & Sæther, B. -E. (2018). Environmental drivers of varying selective optima in a small passerine: A multivariate,

multiepisodic approach. *Evolution*, 72(11), 2325–2342. https://doi.org/10.1111/evo.13610

Garant, D., Hadfield, J. D., Kruuk, L. E. B., & Sheldon, B. C. (2008). Stability of genetic variance and covariance for reproductive characters in the face of climate change in a wild bird population. *Molecular Ecology*, 17(1), 179–188. https://doi.org/10.1111/j.1365-294X.2007.03436.x

Guillaume, F., & Whitlock, M. C. (2007). Effects of migration on the genetic covariance matrix. *Evolution*, 61(10), 2398–2409. https://doi.org/10.1111/j.1558-5646.2007.00193.x

Hadfield, J. D. (2010). MCMC methods for multi-response generalized linear mixed models: The MCMCglmm R package. *Journal of Statistical Software*, 33(2), 1–22.

Hangartner, S., Lasne, C., Sgrò, C. M., Connallon, T., & Monro, K. (2020). Genetic covariances promote climatic adaptation in Australian Drosophila*. *Evolution*, 74(2), 326–337. https://doi.org/10.1111/evo.13831

Hansen, T. F., & Houle, D. (2008). Measuring and comparing evolvability and constraint in multivariate characters. *Journal of Evolutionary Biology*, 21(5), 1201–1219. https://doi.org/10.1111/j.1420-9101.2008.01573.x

Hansen, T. F., Pélabon, C., & Houle, D. (2011). Heritability is not evolvability. *Evolutionary Biology*, 38(3), 258–277. https://doi.org/10.1007/s11692-011-9127-6

Henry, G. A., & Stinchcombe, J. R. (2023). G-matrix stability in clinically diverging populations of an annual weed. *Evolution*, 77(1), 49–62. https://doi.org/10.1093/evolut/qpac005

Hine, E., Chenoweth, S. F., Rundle, H. D., & Blows, M. W. (2009). Characterizing the evolution of genetic variance using genetic covariance tensors. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 364(1523), 1567–1578. https://doi.org/10.1098/rstb.2008.0313

Holt, A. R., Davies, Z. G., Tyler, C., & Staddon, S. (2008). Meta-analysis of the effects of predation on animal prey abundance: Evidence from UK vertebrates. *PLoS One*, 3(6), e2400. https://doi.org/10.1371/journal.pone.0002400

Jones, A. G., Arnold, S. J., & Bürger, R. (2003). Stability of the G-matrix in a population experiencing pleiotropic mutation, stabilizing selection, and genetic drift. *Evolution*, 57(8), 1747–1760. https://doi.org/10.1111/j.0014-3820.2003.tb00583.x

Jones, A. G., Arnold, S. J., & Bürger, R. (2004). Evolution and stability of the G-matrix on a landscape with a moving optimum. *Evolution*, 58(8), 1639–1654. https://doi.org/10.1111/j.0014-3820.2004.tb00450.x

Jones, A. G., Arnold, S. J., & Bürger, R. (2007). The mutation matrix and the evolution of evolvability. *Evolution*, 61(4), 727–745. https://doi.org/10.1111/j.1558-5646.2007.00071.x

Jones, A. G., Bürger, R., & Arnold, S. J. (2014). Epistasis and natural selection shape the mutational architecture of complex traits. *Nature Communications*, 5, 3709. https://doi.org/10.1038/ncomms4709

Kruuk, L. E. B. (2004). Estimating genetic parameters in natural populations using the “Animal Model”. *Proceedings of the Royal Society B: Biological Sciences*, 359, 873–890.

Kruuk, L. E. B., & Hadfield, J. D. (2007). How to separate genetic and environmental causes of similarity between relatives. *Journal of Evolutionary Biology*, 20(5), 1890–1903. https://doi.org/10.1111/j.1420-9101.2007.01377.x

Kruuk, L. E. B., Slate, J., & Wilson, A. J. (2008). New answers for old questions: The evolutionary quantitative genetics of wild animal populations. *Annual Review of Ecology, Evolution, and Systematics*, 39, 525–548.

Krzanowski, W. J. (1979). Between-groups comparison of principal components. *Journal of the American Statistical Association*, 74, 703–707.
Kvist, L., Ruokonen, M., Lumme, J., & Orell, M. (1999). The colonization history and present-day population structure of the European great tit (Parus major major). *Heredity*, *82*(5), 495–502.
https://doi.org/10.1038/sj.hdy.6885130

Laine, V. N., Gossmann, T. I., Schachtschneider, K. M., Garroway, C. J., Madsen, O., Verhoeven, K. J. F., de Jager, V., Megens, H. -J., Warren, W. C., Minx, P., Crooijmans, R. P. M. A., Corcoran, P., Sheldon, B. C., Slate, J., Zeng, K., van Oers, K., Visser, M. E., & Groenen, M. A. M.; Great Tit HapMap Consortium. (2016). Evolutionary signals of selection on cognition from the great tit genome and methylome. *Nature Communications*, *7*, 10474. https://doi.org/10.1038/ncomms10474

Lande, R. (1979). Quantitative genetic analysis of multivariate evolution, applied to brain: Body size allometry. *Evolution*, *33*(1Part2), 402–416. https://doi.org/10.1111/j.1558-5646.1979.tb04694.x

Lande, R. (1992). Neutral theory of quantitative genetic variance in an island model with local extinction and colonization. *Evolution*, *46*(2), 381–389. https://doi.org/10.1111/j.1558-5646.1992.tb02046.x

Le Maître, A., & Mitteroecker, P. (2019). Multivariate comparison of variance in R. *Methods in Ecology and Evolution*, *10*(9), 1380–1392.
https://doi.org/10.1111/2041-210x.13253

Lemoine, M., Lucek, K., Perrier, C., Saladin, V., Adriaensen, F., Barba, E., Belda, E. J., Charmantier, A., Cichoń, M., Eeva, T., Grégoire, A., Hinde, C. A., Johnsen, A., Komdeur, J., Mänd, R., Matthysen, E., Norte, A. C., Pitala, N., Sheldon, B. C., ... Richner, H. (2016). Low but contrasting neutral genetic differentiation shaped by winter temperature in European great tits. *Biological Journal of the Linnean Society*, *118*(3), 668–685. https://doi.org/10.1111/bij.12745

Lundblad, C. G., & Conway, C. J. (2020). Variation in selective regimes drives intraspecific variation in life-history traits and migratory behaviour along an elevational gradient. *The Journal of Animal Ecology*, *89*(2), 397–411. https://doi.org/10.1111/1365-2656.13134

Martinez-Padilla, J., Estrada, A., Early, R., & Garcia-Gonzalez, F. (2017). Evolvability meets biogeography: Evolutionary potential decreases at high and low environmental favourability. *Proceedings of the Royal Society B: Biological Sciences*, *284*(1856), 20170516. https://doi.org/10.1098/rspb.2017.0516

McGlothlin, J. W., Kobiela, M. E., Wright, H. V., Mahler, D. L., Kolbe, J. J., Losos, J. B., & Brodie, E. D. (2018). Adaptive radiation along a deeply conserved genetic line of least resistance in Anolis lizards. *Evolution Letters*, *2*(4), 310–322. https://doi.org/10.1002/evl3.72

Morrissey, M. B., Hangartner, S., & Monro, K. (2019). A note on simulating null distributions for G matrix comparisons. *Evolution*, *73*(12), 2512–2517. https://doi.org/10.1111/evo.13842

Morrissey, M. B., Wilson, J., Pemberton, J. M., & Ferguson, M. M. (2007). A framework for power and sensitivity analyses for quantitative genetic studies of natural populations, and case studies in Soay sheep (*Ovis aries*). *Journal of Evolutionary Biology*, *20*(6), 2309–2321.
https://doi.org/10.1111/j.1420-9101.2007.01412.x

Nakagawa, S. (2015). Missing data: Mechanisms, methods, and messages. In G. A. Fox, S. Negrete-Yankelevich, & V. J. Sosa (Eds.), *Ecological statistics: Contemporary theory and application*. Oxford University Press.

Nosil, P., Crespi, B. J., Sandoval, C. P., & Kirkpatrick, M. (2006). Migration and the genetic covariance between habitat preference and performance. *The American Naturalist*, *167*(3), E66–E78.
https://doi.org/10.1086/499383

Ovaskainen, O., Cano, J. M., & Merilä, J. (2008). A Bayesian framework for comparative quantitative genetics. *Proceedings of the Royal Society B: Biological Sciences*, *275*(1635), 669–678. https://doi.org/10.1098/rspb.2007.0949

Paccard, A., Van Buskirk, J., & Willi, Y. (2016). Quantitative genetic architecture at latitudinal range boundaries: Reduced variation but higher trait independence. *American Naturalist*, *187*(5), 667–677. https://doi.org/10.1086/685643

Pélabon, C., Hilde, C. H., Einum, S., & Gamelon, M. (2020). On the use of the coefficient of variation to quantify and compare trait variation. *Evolution Letters*, *4*(3), 180–188. https://doi.org/10.1002/evl3.171

Phillips, P. C., Whitlock, M. C., & Fowler, K. (2001). Inbreeding changes the shape of the genetic covariance matrix in *Drosophila melanogaster*. *Genetics*, *158*(3), 1137–1145. https://doi.org/10.1093/genetics/158.3.1137

Polechová, J. (2018). Is the sky the limit? On the expansion threshold of a species' range. *PLoS Biology*, *16*(6), e2005372. https://doi.org/10.1371/journal.pbio.2005372

Polechová, J., & Barton, N. H. (2015). Limits to adaptation along environmental gradients. *Proceedings of the National Academy of Sciences of the United States of America*, *112*(20), 6401–6406. https://doi.org/10.1073/pnas.1421515112

Postma, E. (2005). Evolutionary genetics of life-history traits in a structured environment. Understanding variation in clutch size and laying date in great tits (*Parus major*) [PhD thesis], Utrecht University.

Puentes, A., Granath, G., & Ågren, J. (2016). Similarity in G matrix structure among natural populations of *Arabidopsis lyrata*. *Evolution*, *70*(10), 2370–2386. https://doi.org/10.1111/evo.13034

Revell, L. J. (2007). The G matrix under fluctuating correlational mutation and selection. *Evolution*, *61*(8), 1857–1872. https://doi.org/10.1111/j.1558-5646.2007.00161.x

Roff, D. (2000). The evolution of the G matrix: Selection or drift? *Heredity*, *84*(Pt2), 135–142. https://doi.org/10.1046/j.1365-2540.2000.00695.x

Roff, D. A., Mousseau, T., Møller, A. P., Lope, F. de, & Saino, N. (2004). Geographic variation in the G matrices of wild populations of the barn swallow. *Heredity*, *93*, 8–14.

Schluter, D. (1996). Adaptive radiation along genetic lines of least resistance. *Evolution*, *50*(5), 1766–1774. https://doi.org/10.1111/j.1558-5646.1996.tb03563.x

Shaw, F. H., Shaw, R. G., Wilkinson, G. S., & Turelli, M. (1995). Changes in genetic variances and covariances: G Whiz! *Evolution*, *49*(6), 1260–1267. https://doi.org/10.1111/j.1558-5646.1995.tb04452.x

Siepielski, A. M., Morrissey, M. B., Buoro, M., Carlson, S. M., Caruso, C. M., Clegg, S. M., Coulson, T., DiBattista, J., Gotanda, K. M., Francis, C. D., Hereford, J., Kingsolver, J. G., Augustine, K. E., Kruuk, L. E. B., Martin, R. A., Sheldon, B. C., Sletvold, N., Svensson, E. I., Wade, M. J., & MacColl, A. D. C. (2017). Precipitation drives global variation in natural selection. *Science*, *355*(6328), 959–962. https://doi.org/10.1126/science.aag2773

Siepielski, A. M., Morrissey, M. B., Buoro, M., Carlson, S. M., Caruso, C. M., Clegg, S. M., Coulson, T., DiBattista, J., Gotanda, K. M., Francis, C. D., Hereford, J., Kingsolver, J. G., Augustine, K. E., Kruuk, L. E. B., Martin, R. A., Sheldon, B. C., Sletvold, N., Svensson, E. I., Wade, M. J., & MacColl, A. D. C. (2018). Response to comment on "Precipitation drives global variation in natural selection". *Science*, *359*(6374), eea5760. https://doi.org/10.1126/science.aan5760

Sniegula, S., Golab, M. J., Drobniak, S. M., & Johansson, F. (2018). The genetic variance but not the genetic covariance of life-history traits changes towards the north in a time-constrained insect. *Journal of Evolutionary Biology*, *31*(6), 853–865. https://doi.org/10.1111/jeb.13269

Spurgin, L. G., Bosse, M., Adriaensen, F., Albayrak, T., Barboutis, C., Belda, E., Bushuev, A., Cecere, J. G., Charmantier, A., Cichoń, M., Dingemanse, N. J., Doligez, B., Eeva, T., Erikstad, K. E., Fedorov,
Kvist, L., Ruokonen, M., Lumme, J., & Orell, M. (1999). The colonization history and present-day population structure of the European great tit (Parus major major). *Heredity*, *82*(5), 495–502.
https://doi.org/10.1038/sj.hdy.6885130

Laine, V. N., Gossmann, T. I., Schachtschneider, K. M., Garroway, C. J., Madsen, O., Verhoeven, K. J. F., de Jager, V., Megens, H. -J., Warren, W. C., Minx, P., Crooijmans, R. P. M. A., Corcoran, P., Sheldon, B. C., Slate, J., Zeng, K., van Oers, K., Visser, M. E., & Groenen, M. A. M.; Great Tit HapMap Consortium. (2016). Evolutionary signals of selection on cognition from the great tit genome and methylome. *Nature Communications*, *7*, 10474. https://doi.org/10.1038/ncomms10474

Lande, R. (1979). Quantitative genetic analysis of multivariate evolution, applied to brain: Body size allometry. *Evolution*, *33*(1Part2), 402–416. https://doi.org/10.1111/j.1558-5646.1979.tb04694.x

Lande, R. (1992). Neutral theory of quantitative genetic variance in an island model with local extinction and colonization. *Evolution*, *46*(2), 381–389. https://doi.org/10.1111/j.1558-5646.1992.tb02046.x

Le Maître, A., & Mitteroecker, P. (2019). Multivariate comparison of variance in R. *Methods in Ecology and Evolution*, *10*(9), 1380–1392.
https://doi.org/10.1111/2041-210x.13253

Lemoine, M., Lucek, K., Perrier, C., Saladin, V., Adriaensen, F., Barba, E., Belda, E. J., Charmantier, A., Cichoń, M., Eeva, T., Grégoire, A., Hinde, C. A., Johnsen, A., Komdeur, J., Mänd, R., Matthysen, E., Norte, A. C., Pitala, N., Sheldon, B. C., ... Richner, H. (2016). Low but contrasting neutral genetic differentiation shaped by winter temperature in European great tits. *Biological Journal of the Linnean Society*, *118*(3), 668–685. https://doi.org/10.1111/bij.12745

Lundblad, C. G., & Conway, C. J. (2020). Variation in selective regimes drives intraspecific variation in life-history traits and migratory behaviour along an elevational gradient. *The Journal of Animal Ecology*, *89*(2), 397–411. https://doi.org/10.1111/1365-2656.13134

Martínez-Padilla, J., Estrada, A., Early, R., & Garcia-Gonzalez, F. (2017). Evolvability meets biogeography: Evolutionary potential decreases at high and low environmental favourability. *Proceedings of the Royal Society B: Biological Sciences*, *284*(1856), 20170516. https://doi.org/10.1098/rspb.2017.0516

McGlothlin, J. W., Kobiela, M. E., Wright, H. V., Mahler, D. L., Kolbe, J. J., Losos, J. B., & Brodie, E. D. (2018). Adaptive radiation along a deeply conserved genetic line of least resistance in Anolis lizards. *Evolution Letters*, *2*(4), 310–322. https://doi.org/10.1002/evl3.72

Morrissey, M. B., Hangartner, S., & Monro, K. (2019). A note on simulating null distributions for G matrix comparisons. *Evolution*, *73*(12), 2512–2517. https://doi.org/10.1111/evo.13842

Morrissey, M. B., Wilson, J., Pemberton, J. M., & Ferguson, M. M. (2007). A framework for power and sensitivity analyses for quantitative genetic studies of natural populations, and case studies in Soay sheep (*Ovis aries*). *Journal of Evolutionary Biology*, *20*(6), 2309–2321.
https://doi.org/10.1111/j.1420-9101.2007.01412.x

Nakagawa, S. (2015). Missing data: Mechanisms, methods, and messages. In G. A. Fox, S. Negrete-Yankelevich, & V. J. Sosa (Eds.), *Ecological statistics: Contemporary theory and application*. Oxford University Press.

Nosil, P., Crespi, B. J., Sandoval, C. P., & Kirkpatrick, M. (2006). Migration and the genetic covariance between habitat preference and performance. *The American Naturalist*, *167*(3), E66–E78.
https://doi.org/10.1086/499383

Ovaskainen, O., Cano, J. M., & Merilä, J. (2008). A Bayesian framework for comparative quantitative genetics. *Proceedings of the Royal Society B: Biological Sciences*, *275*(1635), 669–678. https://doi.org/10.1098/rspb.2007.0949

Paccard, A., Van Buskirk, J., & Willi, Y. (2016). Quantitative genetic architecture at latitudinal range boundaries: Reduced variation but higher trait independence. *American Naturalist*, *187*(5), 667–677. https://doi.org/10.1086/685643

Pélabon, C., Hilde, C. H., Einum, S., & Gamelon, M. (2020). On the use of the coefficient of variation to quantify and compare trait variation. *Evolution Letters*, *4*(3), 180–188. https://doi.org/10.1002/evl3.171

Phillips, P. C., Whitlock, M. C., & Fowler, K. (2001). Inbreeding changes the shape of the genetic covariance matrix in *Drosophila melanogaster*. *Genetics*, *158*(3), 1137–1145. https://doi.org/10.1093/genetics/158.3.1137

Polechová, J. (2018). Is the sky the limit? On the expansion threshold of a species' range. *PLoS Biology*, *16*(6), e2005372. https://doi.org/10.1371/journal.pbio.2005372

Polechová, J., & Barton, N. H. (2015). Limits to adaptation along environmental gradients. *Proceedings of the National Academy of Sciences of the United States of America*, *112*(20), 6401–6406. https://doi.org/10.1073/pnas.1421515112

Postma, E. (2005). Evolutionary genetics of life-history traits in a structured environment. Understanding variation in clutch size and laying date in great tits (*Parus major*) [PhD thesis], Utrecht University.

Puentes, A., Granath, G., & Ågren, J. (2016). Similarity in G matrix structure among natural populations of *Arabidopsis lyrata*. *Evolution*, *70*(10), 2370–2386. https://doi.org/10.1111/evo.13034

Revell, L. J. (2007). The G matrix under fluctuating correlational mutation and selection. *Evolution*, *61*(8), 1857–1872. https://doi.org/10.1111/j.1558-5646.2007.00161.x

Roff, D. (2000). The evolution of the G matrix: Selection or drift? *Heredity*, *84*(Pt2), 135–142. https://doi.org/10.1046/j.1365-2540.2000.00695.x

Roff, D. A., Mousseau, T., Møller, A. P., Lope, F. de, & Saino, N. (2004). Geographic variation in the G matrices of wild populations of the barn swallow. *Heredity*, *93*, 8–14.

Schluter, D. (1996). Adaptive radiation along genetic lines of least resistance. *Evolution*, *50*(5), 1766–1774. https://doi.org/10.1111/j.1558-5646.1996.tb03563.x

Shaw, F. H., Shaw, R. G., Wilkinson, G. S., & Turelli, M. (1995). Changes in genetic variances and covariances: G Whiz! *Evolution*, *49*(6), 1260–1267. https://doi.org/10.1111/j.1558-5646.1995.tb04452.x

Siepielski, A. M., Morrissey, M. B., Buoro, M., Carlson, S. M., Caruso, C. M., Clegg, S. M., Coulson, T., DiBattista, J., Gotanda, K. M., Francis, C. D., Hereford, J., Kingsolver, J. G., Augustine, K. E., Kruuk, L. E. B., Martin, R. A., Sheldon, B. C., Sletvold, N., Svensson, E. I., Wade, M. J., & MacColl, A. D. C. (2017). Precipitation drives global variation in natural selection. *Science*, *355*(6328), 959–962. https://doi.org/10.1126/science.aag2773

Siepielski, A. M., Morrissey, M. B., Buoro, M., Carlson, S. M., Caruso, C. M., Clegg, S. M., Coulson, T., DiBattista, J., Gotanda, K. M., Francis, C. D., Hereford, J., Kingsolver, J. G., Augustine, K. E., Kruuk, L. E. B., Martin, R. A., Sheldon, B. C., Sletvold, N., Svensson, E. I., Wade, M. J., & MacColl, A. D. C. (2018). Response to comment on "Precipitation drives global variation in natural selection". *Science*, *359*(6374), eaan5760. https://doi.org/10.1126/science.aan5760

Sniegula, S., Golab, M. J., Drobniak, S. M., & Johansson, F. (2018). The genetic variance but not the genetic covariance of life-history traits changes towards the north in a time-constrained insect. *Journal of Evolutionary Biology*, *31*(6), 853–865. https://doi.org/10.1111/jeb.13269

Spurgin, L. G., Bosse, M., Adriaensen, F., Albayrak, T., Barboutis, C., Belda, E., Bushuev, A., Cecere, J. G., Charmantier, A., Cichoń, M., Dingemanse, N. J., Doligez, B., Eeva, T., Erikstad, K. E., Fedorov,
V., Griggio, M., Heylen, D., Hille, S., Hinde, C. A., ... Slate, J. (2019). The great tit HapMap project: A continental-scale analysis of genomic variation in a songbird. *bioRxiv* 561399. https://doi.org/10.1101/561399

Steppan, S.J., Phillips, P.C., & Houle, D. (2002). Comparative quantitative genetics: Evolution of the G matrix. *Trends in Ecology & Evolution*, 17(7), 320-327. https://doi.org/10.1016/s0169-5347(02)02505-3

Stinchcombe, J. R., & Kirkpatrick, M. (2012). Genetics and evolution of function-valued traits: Understanding environmentally responsive phenotypes. *Trends in Ecology & Evolution*, 27, 637-647.

Teplitsky, C., Mouawad, N. G., Balbontin, J., De Lope, F., & Møller, A. P. (2011). Quantitative genetics of migration syndromes: A study of two barn swallow populations. *Journal of Evolutionary Biology*, 24(9), 2025-2039. https://doi.org/10.1111/j.1420-9101.2011.02342.x

Teplitsky, C., Tarka, M., Møller, A. P., Nakagawa, S., Balbontín, J., Burke, T. A., Doutrelant, C., Gregoire, A., Hansson, B., Hasselquist, D., Gustafsson, L., de Lope, F., Marzal, A., Mills, J. A., Wheelwright, N. T., Yarrall, J. W., & Charmantier, A. (2014). Assessing multivariate constraints to evolution across ten long-term avian studies. *PLoS One*, 9(3), e90444. https://doi.org/10.1371/journal.pone.0090444

Vágási, C., Pap, P., Vincze, O., Osváth, G., Erritzøe, J., & Moller, A. (2016). Morphological adaptations to migration in birds. *Evolutionary Biology*, 43, 48-59.

De Villemereuil, P., Charmantier, A., Arlt, D., Bize, P., Brekke, P., Brouwer, L., Cockburn, A., Côté, S. D., Dobson, F. S., Evans, S. R., Festa-Bianchet, M., Gamelon, M., Hamel, S., Hegelbach, J., Jerstad, K., Kempenaers, B., Kruuk, L. E. B., Kumpula, J., Kvalnes, T., ... Chevin, L. M. (2020). Fluctuating optimum and temporally variable selection on breeding date in birds and mammals. *Proceedings of the National Academy of Sciences of the United States of America*, 117, 31969-31978.

Walsh, B., & Lynch, M. (2018). *Evolution and selection of quantitative traits*. Oxford University Press.

Walter, G. M., Aguirre, J. D., Blows, M. W., & Ortiz-Barrientos, D. (2018). Evolution of genetic variance during adaptive radiation. *The American Naturalist*, 191(4), E108-E128. https://doi.org/10.1086/696123

Wood, C. W., & Brodie, E. D. (2015). Environmental effects on the structure of the G-matrix. *Evolution*, 69(11), 2927-2940. https://doi.org/10.1111/evo.12795
V., Griggio, M., Heylen, D., Hille, S., Hinde, C. A., ... Slate, J. (2019). The great tit HapMap project: A continental-scale analysis of genomic variation in a songbird. *bioRxiv* 561399. https://doi.org/10.1101/561399

Steppan, S.J., Phillips, P.C., & Houle, D. (2002). Comparative quantitative genetics: Evolution of the G matrix. *Trends in Ecology & Evolution*, 17(7), 320-327. https://doi.org/10.1016/s0169-5347(02)02505-3

Stinchcombe, J. R., & Kirkpatrick, M. (2012). Genetics and evolution of function-valued traits: Understanding environmentally responsive phenotypes. *Trends in Ecology & Evolution*, 27, 637-647.

Teplitsky, C., Mouawad, N. G., Balbontin, J., De Lope, F., & Møller, A. P. (2011). Quantitative genetics of migration syndromes: A study of two barn swallow populations. *Journal of Evolutionary Biology*, 24(9), 2025-2039. https://doi.org/10.1111/j.1420-9101.2011.02342.x

Teplitsky, C., Tarka, M., Møller, A. P., Nakagawa, S., Balbontín, J., Burke, T. A., Doutrelant, C., Gregoire, A., Hansson, B., Hasselquist, D., Gustafsson, L., de Lope, F., Marzal, A., Mills, J. A., Wheelwright, N. T., Yarrall, J. W., & Charmantier, A. (2014). Assessing multivariate constraints to evolution across ten long-term avian studies. *PLoS One*, 9(3), e90444. https://doi.org/10.1371/journal.pone.0090444

Vágási, C., Pap, P., Vincze, O., Osváth, G., Erritzøe, J., & Moller, A. (2016). Morphological adaptations to migration in birds. *Evolutionary Biology*, 43, 48-59.

De Villemereuil, P., Charmantier, A., Arlt, D., Bize, P., Brekke, P., Brouwer, L., Cockburn, A., Côté, S. D., Dobson, F. S., Evans, S. R., Festa-Bianchet, M., Gamelon, M., Hamel, S., Hegelbach, J., Jerstad, K., Kempenaers, B., Kruuk, L. E. B., Kumpula, J., Kvalnes, T., ... Chevin, L. M. (2020). Fluctuating optimum and temporally variable selection on breeding date in birds and mammals. *Proceedings of the National Academy of Sciences of the United States of America*, 117, 31969-31978.

Walsh, B., & Lynch, M. (2018). *Evolution and selection of quantitative traits*. Oxford University Press.

Walter, G. M., Aguirre, J. D., Blows, M. W., & Ortiz-Barrientos, D. (2018). Evolution of genetic variance during adaptive radiation. *The American Naturalist*, 191(4), E108-E128. https://doi.org/10.1086/696123

Wood, C. W., & Brodie, E. D. (2015). Environmental effects on the structure of the G-matrix. *Evolution*, 69(11), 2927-2940. https://doi.org/10.1111/evo.12795