# Low but significant evolutionary potential for growth, phenology and reproduction traits in European beech

Marjana Westergren, Juliette Archambeau, Marko Bajc, Rok Damjanić, Adélaïde Theraroz, Hojka Kraigher, Sylvie Oddou-muratorio, Santiago C González-martínez

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

## To cite this version:

Marjana Westergren, Juliette Archambeau, Marko Bajc, Rok Damjanić, Adélaïde Theraroz, et al. Low but significant evolutionary potential for growth, phenology and reproduction traits in European beech. Molecular Ecology, 2023, 10.1111/mec.17196. hal-04598151

HAL Id: hal-04598151

https://hal.science/hal-04598151

Submitted on 5 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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>
# Low but significant evolutionary potential for growth, phenology and reproduction traits in European beech

Marjana Westergren, Juliette Archambeau, Marko Bajc, Rok Damjanić, Adélaïde Theraroz, Hojka Kraigher, Sylvie Oddou-muratorio, Santiago C González-martínez

► To cite this version:

Marjana Westergren, Juliette Archambeau, Marko Bajc, Rok Damjanić, Adélaïde Theraroz, et al. Low but significant evolutionary potential for growth, phenology and reproduction traits in European beech. Molecular Ecology, 2023, 10.1111/mec.17196. hal-04598151

HAL Id: hal-04598151

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

Submitted on 5 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Low but significant evolutionary potential for growth, phenology and reproduction traits in European beech

Marjana Westergren<fnref n="1" /> ID

Juliette Archambeau<fnref n="2" /> ID

Marko Bajc<fnref n="1" />

Rok Damjanić<fnref n="1" />

Adélaïde Theraroz<fnref n="2" />

Hojka Kraigher<fnref n="1" />

Sylvie Oddou-Muratorio<fnref n="3" />,<fnref n="4" />

Santiago C. González-Martínez<fnref n="2" /> ID

### Correspondence

Marjana Westergren, Slovenian Forestry Institute, Večna Pot 2, Ljubljana 1000, Slovenia.

Email: marjana.westergren@gozdis.si

### Funding information

European Cooperation in Science and Technology; European Union H2020, Grant/Award Number: 862221; European Union LIFE Programme, Grant/Award Number: ENV/SI/000148; LabEx COTE; Slovenian Research Agency, Grant/Award Number: P4-0107

Handling Editor: Victoria L. Sork

### Abstract

Local survival of forest tree populations under climate change depends on existing genetic variation and their adaptability to changing environments. Responses to selection were studied in European beech (*Fagus sylvatica*) under field conditions. A total of 1087 adult trees, seeds, 1-year-old seedlings and established multiyear saplings were genotyped with 16 nuSSRs. Adult trees were assessed for phenotypic traits related to growth, phenology and reproduction. Parentage and paternity analyses were used to estimate effective female and male fecundity as a proxy of fitness and showed that few parents contributed to successful regeneration. Selection gradients were estimated from the relationship between traits and fecundity, while heritability and evolvability were estimated using mixed models and the breeder's equation. Larger trees bearing more fruit and early male flowering had higher total fecundity, while trees with longer growth season had lower total fecundity (directional selection). Stabilizing selection on spring phenology was found for female fecundity, highlighting the role of late frosts as a selection driver. Selection gradients for other traits varied between measurement years and the offspring cohort used to estimate parental fecundity. Compared to other studies in natural populations, we found low to moderate heritability and evolvability for most traits. Response to selection was higher for growth than for budburst, leaf senescence or reproduction traits, reflecting more consistent selection gradients across years and sex functions, and higher phenotypic variability in the population. Our study provides empirical evidence suggesting that populations of long-lived organisms such as forest trees can adapt locally, even at short-time scales.

### KEYWORDS

climate change, *Fagus sylvatica*, heritability, in situ adaptation, response to selection, selection gradients

<footnote n="1">¹Slovenian Forestry Institute, Ljubljana, Slovenia</footnote>

<footnote n="2">²INRAE, Univ. Bordeaux, BIOGECO, Cestas, France</footnote>

<footnote n="3">³INRAE, URFM, Avignon, France</footnote>

<footnote n="4">⁴INRAE, Univ. de Pau et des Pays de l'Adour, E2S UPPA, ECOBIOP, Saint-Pée-sur-Nivelle, France</footnote>
<figure><img src="image_2.png" /></figure>

# Low but significant evolutionary potential for growth, phenology and reproduction traits in European beech

Marjana Westergren<fnref n="1" /> ID

Juliette Archambeau<fnref n="2" /> ID

Marko Bajc<fnref n="1" />

Rok Damjanić<fnref n="1" />

Adélaïde Theraroz<fnref n="2" />

Hojka Kraigher<fnref n="1" />

Sylvie Oddou-Muratorio<fnref n="3" />,<fnref n="4" />

Santiago C. González-Martínez<fnref n="2" /> ID

### Correspondence

Marjana Westergren, Slovenian Forestry Institute, Večna Pot 2, Ljubljana 1000, Slovenia.

Email: marjana.westergren@gozdis.si

### Funding information

European Cooperation in Science and Technology; European Union H2020, Grant/Award Number: 862221; European Union LIFE Programme, Grant/Award Number: ENV/SI/000148; LabEx COTE; Slovenian Research Agency, Grant/Award Number: P4-0107

Handling Editor: Victoria L. Sork

### Abstract

Local survival of forest tree populations under climate change depends on existing genetic variation and their adaptability to changing environments. Responses to selection were studied in European beech (Fagus sylvatica) under field conditions. A total of 1087 adult trees, seeds, 1-year-old seedlings and established multiyear saplings were genotyped with 16 nuSSRs. Adult trees were assessed for phenotypic traits related to growth, phenology and reproduction. Parentage and paternity analyses were used to estimate effective female and male fecundity as a proxy of fitness and showed that few parents contributed to successful regeneration. Selection gradients were estimated from the relationship between traits and fecundity, while heritability and evolvability were estimated using mixed models and the breeder's equation. Larger trees bearing more fruit and early male flowering had higher total fecundity, while trees with longer growth season had lower total fecundity (directional selection). Stabilizing selection on spring phenology was found for female fecundity, highlighting the role of late frosts as a selection driver. Selection gradients for other traits varied between measurement years and the offspring cohort used to estimate parental fecundity. Compared to other studies in natural populations, we found low to moderate heritability and evolvability for most traits. Response to selection was higher for growth than for budburst, leaf senescence or reproduction traits, reflecting more consistent selection gradients across years and sex functions, and higher phenotypic variability in the population. Our study provides empirical evidence suggesting that populations of long-lived organisms such as forest trees can adapt locally, even at short-time scales.

### KEYWORDS

climate change, *Fagus sylvatica*, heritability, in situ adaptation, response to selection, selection gradients

<footnote n="1">¹Slovenian Forestry Institute, Ljubljana, Slovenia</footnote>

<footnote n="2">²INRAE, Univ. Bordeaux, BIOGECO, Cestas, France</footnote>

<footnote n="3">³INRAE, URFM, Avignon, France</footnote>

<footnote n="4">⁴INRAE, Univ. de Pau et des Pays de l'Adour, E2S UPPA, ECOBIOP, Saint-Pée-sur-Nivelle, France</footnote>
# 1 | INTRODUCTION

Forest trees, the long-lived foundation species of forest ecosystems, modulate important ecosystem services, such as habitat provision, carbon storage, biomass and oxygen production, and support forest biodiversity (Ellison, 2019). Different sources of evidence suggest that tree populations will likely not be able to track their ecological optima given the projected rate and magnitude of climate change (Petit et al., 2008; Petit-Cailleux et al., 2021). Instead, they will need to adapt in situ to novel climatic conditions and disrupted species interactions (e.g. new pests and diseases; Ramsfield et al., 2016), or go extinct. Insights on the short-term adaptive potential of forest trees is thus needed to predict the evolutionary response of entire forest ecosystems to climate change (Hoffmann & Sgrò, 2011) and the utility of assisted gene flow (Aitken & Whitlock, 2013).

Evolutionary adaptive change can only take place when two conditions are fulfilled; firstly, there should be enough heritable phenotypic variation in a population. Secondly, this variation must be correlated to fitness, that is, the individual's ability to survive and reproduce. Under particular environmental change drivers, for example, strong selection pressures associated with stochastic extreme events, adaptation may then take place. In long-lived species such as forest trees adaptation may be slowed due to interannual climate variability, while forest management would often favour progeny from a single or a few reproduction events. Still, trees may adapt faster than previously thought, as some limited experimental evidence suggests (e.g. Petit et al., 2004). For example, in six common gardens where haphazard management resulted in spreading of naturalized populations of non-native loblolly pine (Pinus taeda L.) during the last 50 years, significant shifts in allele frequencies in 25 genes were observed along the invasion front, the gene-level patterns of evolution being population-specific (Zenni & Hoban, 2015). In the same system, Zenni et al. (2016) found evidence of rapid evolution in growth rates and a lack of trade-offs between growth and defence traits. In addition, recent theoretical models based on polygenic trait architectures, which seem common in forest trees (e.g. de Miguel et al., 2022), also predict rapid adaptation resulting from small allele frequency shifts in multiple loci (Remington, 2015).

Ideally, evolutionary adaptive change should be estimated by studying the same population(s) over time, either through molecular markers or phenotypic traits (Hansen et al., 2012). But for long-lived temperate forest trees, long-term monitoring studies are difficult because generations change over centuries rather than years or a few decades, even when different developmental stages may occur simultaneously at the same site (adult trees, saplings, seedlings and seeds). Currently, monitoring programmes are being established (e.g. Bajc et al., 2021), but until fully implemented, a valid approach is still to combine pedigree-tracing using molecular markers and standard computation of quantitative genetic parameters using the animal model (e.g. heritability), with estimation of selection gradients, that is, regression coefficients of relative fitness on a trait (Lande

& Arnold, 1983). Subsequently, the predicted response to selection can be calculated using the breeder's equation (Hansen et al., 2012). An advantage of using selection gradients is that they measure both the strength and direction of selection and, when regressing multiple traits, are (at least theoretically) able to distinguish direct selection on the trait from indirect selection resulting from correlation with other traits (Lande & Arnold, 1983). They are also standardized, either by the standard deviation or the phenotypic mean, which provides a dimensionless measure of selection strength that can be compared across traits, populations or species (Matsumura et al., 2012).

Because natural selection occurs in nature, key evolutionary parameters should be studied in situ and under multiple environments. However, cumulative lifetime fitness of iteroparous trees is extremely difficult to assess, especially in the field (Barringer et al., 2013), and even measuring fitness components (e.g. survival to a given life stage, mating success, fecundity under a certain period) is challenging. Seed set, abundant during mast seeding, is easily observable, but only a tenth of this seed will germinate and 1 in 1000 germinated seeds will result in a sapling after a few years (Davis et al., 2005). A more robust approach is to use effective fecundity as fitness proxy, which can be estimated through parentage analysis of established saplings (Geber & Griffen, 2003; Orr, 2009), that is, the saplings that have already survived the initial selection phases and have a high probability to form the core of the next generation. In addition, measuring phenotypic traits and controlling for confounding environmental variation is more difficult in nature than in common gardens or greenhouses, so such studies are rare for long-lived tree species and often contradict results from controlled environments. In maritime pine, for example, Alía et al. (2014) found significant differences between selection gradients and expected responses to selection examined for the same traits in growth chambers and near-natural conditions in the field. Estimates of narrow-sense heritability ($h^2$) also tend to differ under controlled and natural conditions and do not follow a clear pattern. In the above experiment, heritability estimates were up to ten times lower under near-natural conditions than in growth chambers for traits related to seedling mass, but not for height or survival, where heritability was higher under near-natural conditions. Heritability for serotiny, a fire-related reproduction trait, was twofold lower in Aleppo and maritime pines in nature than in a common garden (Castellanos et al., 2015), while heritability of budburst for sessile and pedunculate oaks and beech was higher in nature than in common gardens (Alexandre, Truffaut, Ducousso, et al., 2020; Bontemps et al., 2016; Gauzere et al., 2016).

This study aims at estimating the evolutionary potential of a European beech (Fagus sylvatica L.) population at the core of its range. This major broadleaved tree species, which is monoecious and protogynous (Merzeau et al., 1994) and does not maintain a seed bank in the soil (Packham et al., 2012), dominates many lowland and mountain forests across Europe. It requires a humid atmosphere with precipitation throughout the year and tolerates winter cold but is sensitive to late spring frosts (Houston Durrant et al., 2016). Moreover, extensive genetic knowledge based on previous studies
# 1 | INTRODUCTION

Forest trees, the long-lived foundation species of forest ecosystems, modulate important ecosystem services, such as habitat provision, carbon storage, biomass and oxygen production, and support forest biodiversity (Ellison, 2019). Different sources of evidence suggest that tree populations will likely not be able to track their ecological optima given the projected rate and magnitude of climate change (Petit et al., 2008; Petit-Cailleux et al., 2021). Instead, they will need to adapt in situ to novel climatic conditions and disrupted species interactions (e.g. new pests and diseases; Ramsfield et al., 2016), or go extinct. Insights on the short-term adaptive potential of forest trees is thus needed to predict the evolutionary response of entire forest ecosystems to climate change (Hoffmann & Sgrò, 2011) and the utility of assisted gene flow (Aitken & Whitlock, 2013).

Evolutionary adaptive change can only take place when two conditions are fulfilled; firstly, there should be enough heritable phenotypic variation in a population. Secondly, this variation must be correlated to fitness, that is, the individual's ability to survive and reproduce. Under particular environmental change drivers, for example, strong selection pressures associated with stochastic extreme events, adaptation may then take place. In long-lived species such as forest trees adaptation may be slowed due to interannual climate variability, while forest management would often favour progeny from a single or a few reproduction events. Still, trees may adapt faster than previously thought, as some limited experimental evidence suggests (e.g. Petit et al., 2004). For example, in six common gardens where haphazard management resulted in spreading of naturalized populations of non-native loblolly pine (Pinus taeda L.) during the last 50 years, significant shifts in allele frequencies in 25 genes were observed along the invasion front, the gene-level patterns of evolution being population-specific (Zenni & Hoban, 2015). In the same system, Zenni et al. (2016) found evidence of rapid evolution in growth rates and a lack of trade-offs between growth and defence traits. In addition, recent theoretical models based on polygenic trait architectures, which seem common in forest trees (e.g. de Miguel et al., 2022), also predict rapid adaptation resulting from small allele frequency shifts in multiple loci (Remington, 2015).

Ideally, evolutionary adaptive change should be estimated by studying the same population(s) over time, either through molecular markers or phenotypic traits (Hansen et al., 2012). But for long-lived temperate forest trees, long-term monitoring studies are difficult because generations change over centuries rather than years or a few decades, even when different developmental stages may occur simultaneously at the same site (adult trees, saplings, seedlings and seeds). Currently, monitoring programmes are being established (e.g. Bajc et al., 2021), but until fully implemented, a valid approach is still to combine pedigree-tracing using molecular markers and standard computation of quantitative genetic parameters using the animal model (e.g. heritability), with estimation of selection gradients, that is, regression coefficients of relative fitness on a trait (Lande

& Arnold, 1983). Subsequently, the predicted response to selection can be calculated using the breeder's equation (Hansen et al., 2012). An advantage of using selection gradients is that they measure both the strength and direction of selection and, when regressing multiple traits, are (at least theoretically) able to distinguish direct selection on the trait from indirect selection resulting from correlation with other traits (Lande & Arnold, 1983). They are also standardized, either by the standard deviation or the phenotypic mean, which provides a dimensionless measure of selection strength that can be compared across traits, populations or species (Matsumura et al., 2012).

Because natural selection occurs in nature, key evolutionary parameters should be studied in situ and under multiple environments. However, cumulative lifetime fitness of iteroparous trees is extremely difficult to assess, especially in the field (Barringer et al., 2013), and even measuring fitness components (e.g. survival to a given life stage, mating success, fecundity under a certain period) is challenging. Seed set, abundant during mast seeding, is easily observable, but only a tenth of this seed will germinate and 1 in 1000 germinated seeds will result in a sapling after a few years (Davis et al., 2005). A more robust approach is to use effective fecundity as fitness proxy, which can be estimated through parentage analysis of established saplings (Geber & Griffen, 2003; Orr, 2009), that is, the saplings that have already survived the initial selection phases and have a high probability to form the core of the next generation. In addition, measuring phenotypic traits and controlling for confounding environmental variation is more difficult in nature than in common gardens or greenhouses, so such studies are rare for long-lived tree species and often contradict results from controlled environments. In maritime pine, for example, Alía et al. (2014) found significant differences between selection gradients and expected responses to selection examined for the same traits in growth chambers and near-natural conditions in the field. Estimates of narrow-sense heritability ($h^2$) also tend to differ under controlled and natural conditions and do not follow a clear pattern. In the above experiment, heritability estimates were up to ten times lower under near-natural conditions than in growth chambers for traits related to seedling mass, but not for height or survival, where heritability was higher under near-natural conditions. Heritability for serotiny, a fire-related reproduction trait, was twofold lower in Aleppo and maritime pines in nature than in a common garden (Castellanos et al., 2015), while heritability of budburst for sessile and pedunculate oaks and beech was higher in nature than in common gardens (Alexandre, Truffaut, Ducousso, et al., 2020; Bontemps et al., 2016; Gauzere et al., 2016).

This study aims at estimating the evolutionary potential of a European beech (Fagus sylvatica L.) population at the core of its range. This major broadleaved tree species, which is monoecious and protogynous (Merzeau et al., 1994) and does not maintain a seed bank in the soil (Packham et al., 2012), dominates many lowland and mountain forests across Europe. It requires a humid atmosphere with precipitation throughout the year and tolerates winter cold but is sensitive to late spring frosts (Houston Durrant et al., 2016). Moreover, extensive genetic knowledge based on previous studies
in common gardens allows to formulate hypotheses on the phenotypic traits contributing to beech adaptation. For example, phenological traits (e.g. budburst and leaf senescence dates) show genetic differentiation along latitudinal or elevational clines at various spatial scales, with populations of higher elevation or latitude flushing earlier than populations from low elevation or latitude in common garden conditions (Gauzere et al., 2020; Gömöry & Paule, 2011; Vitasse et al., 2009). These studies thus suggest that beech populations evolved different phenological trait values across the species range since the course of post-glacial recolonization. Based on current ecological and genetic knowledge for European beech, we expect significant selection gradients and responses to selection for life history/phenology traits, including reproduction (e.g. flowering time and abundance, and seed production), growth phenology (budburst or leaf unfolding, and leaf senescence) and growth traits, as detailed in Figure 1 and below.

Firstly, larger individuals are expected to have higher reproductive success (Klinkhamer et al., 1997) because they have access to more resources and can produce more gametes (budget effects). In wind-pollinated plants, large size may also benefit males by allowing the liberation of pollen from a greater height, fostering its dispersal and increasing mating success (Klinkhamer et al., 1997; Tonnabel

et al., 2019). Consistently, many studies reported that fecundity increased with tree size (e.g. Oddou-Muratorio et al., 2018, 2021; Wickert et al., 2017). Therefore, we expect directional selection for tree size and amount of fructification.

Secondly, early budburst is expected to increase the length of the vegetative period and spring photosynthesis (Keenan et al., 2014; Richardson et al., 2006) and thereby the resources available for growth and reproduction of the current year. However, early budburst also increases the risk of late frost damage (Augspurger, 2013; Bigler & Bugmann, 2018), which can result in reduced growth and increased risk of mortality (Kreyling et al., 2014; Menzel et al., 2015). We expect directional selection for early budburst, modulated by late frost risk (higher risk would result in selection against early budburst) leading to perhaps stabilizing selection.

Thirdly, early flowering plants are generally favoured in short-lived herbaceous plants (Munguía-Rosas et al., 2011) but not in *Quercus lobata*, where trees that flowered early or late set fewer acorns than those that flowered at the population's peak, when pollen is most abundant (Koenig et al., 2012). Maximizing the overlap in flowering with potential mates, such as during mast seeding years (Hacket-Pain & Bogdziewicz, 2021; Pearse et al., 2016), is likely to increase both male and female reproductive success, particularly if

<figure><img src="image_4.png" /><figcaption>**FIGURE 1** Expectations of selection gradients for life history/phenology traits in protogynous European beech. Pane a shows expectation of selection gradients throughout the trees' lifespan and pane b selection gradients modulated by temporal dimension within the growing season. Grey is the current situation and black the expected situation after evolutionary change; the arrow indicates the direction of selection. For tree size, seed amount, senescence and budburst we expect directional selection, with budburst also modulated by late spring frosts. For flowering, we expect stabilizing selection around the population average with a component of directional selection for earlier male flowering. On the graphs, the X axis indicates the phenotypic trait value, Y axis its frequency in a population and $\omega$ is the fitness.</figcaption></figure>
in common gardens allows to formulate hypotheses on the phenotypic traits contributing to beech adaptation. For example, phenological traits (e.g. budburst and leaf senescence dates) show genetic differentiation along latitudinal or elevational clines at various spatial scales, with populations of higher elevation or latitude flushing earlier than populations from low elevation or latitude in common garden conditions (Gauzere et al., 2020; Gömöry & Paule, 2011; Vitasse et al., 2009). These studies thus suggest that beech populations evolved different phenological trait values across the species range since the course of post-glacial recolonization. Based on current ecological and genetic knowledge for European beech, we expect significant selection gradients and responses to selection for life history/phenology traits, including reproduction (e.g. flowering time and abundance, and seed production), growth phenology (budburst or leaf unfolding, and leaf senescence) and growth traits, as detailed in Figure 1 and below.

Firstly, larger individuals are expected to have higher reproductive success (Klinkhamer et al., 1997) because they have access to more resources and can produce more gametes (budget effects). In wind-pollinated plants, large size may also benefit males by allowing the liberation of pollen from a greater height, fostering its dispersal and increasing mating success (Klinkhamer et al., 1997; Tonnabel

et al., 2019). Consistently, many studies reported that fecundity increased with tree size (e.g. Oddou-Muratorio et al., 2018, 2021; Wickert et al., 2017). Therefore, we expect directional selection for tree size and amount of fructification.

Secondly, early budburst is expected to increase the length of the vegetative period and spring photosynthesis (Keenan et al., 2014; Richardson et al., 2006) and thereby the resources available for growth and reproduction of the current year. However, early budburst also increases the risk of late frost damage (Augspurger, 2013; Bigler & Bugmann, 2018), which can result in reduced growth and increased risk of mortality (Kreyling et al., 2014; Menzel et al., 2015). We expect directional selection for early budburst, modulated by late frost risk (higher risk would result in selection against early budburst) leading to perhaps stabilizing selection.

Thirdly, early flowering plants are generally favoured in short-lived herbaceous plants (Munguía-Rosas et al., 2011) but not in *Quercus lobata*, where trees that flowered early or late set fewer acorns than those that flowered at the population's peak, when pollen is most abundant (Koenig et al., 2012). Maximizing the overlap in flowering with potential mates, such as during mast seeding years (Hacket-Pain & Bogdziewicz, 2021; Pearse et al., 2016), is likely to increase both male and female reproductive success, particularly if

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

**FIGURE 1** Expectations of selection gradients for life history/phenology traits in protogynous European beech. Pane a shows expectation of selection gradients throughout the trees' lifespan and pane b selection gradients modulated by temporal dimension within the growing season. Grey is the current situation and black the expected situation after evolutionary change; the arrow indicates the direction of selection. For tree size, seed amount, senescence and budburst we expect directional selection, with budburst also modulated by late spring frosts. For flowering, we expect stabilizing selection around the population average with a component of directional selection for earlier male flowering. On the graphs, the X axis indicates the phenotypic trait value, Y axis its frequency in a population and $\omega$ is the fitness.
the availability of ovules/pollen is limited (e.g. in isolated populations or under extreme climatic events). In protogynous beech, we expect a combination of stabilizing selection for female/male flowering date around the population average, and overall population directional selection for early flowering, in particular for males.

Lastly, the impact of autumn phenology (leaf senescence) on fitness is poorly understood. However, late leaf senescence would increase autumn photosynthesis (Toomey et al., 2015), and thereby the level of resources available for trees to initiate next-year reproductive buds (Hoch et al., 2013). We, therefore, expect directional selection for late senescence, in order to increase available resources for the following growth and reproductive season.

In this study, we tested these expectations by computing selection gradients based on effective fecundity for a wide variety of reproduction, phenology and growth traits. Moreover, the evolutionary relevance of significant selection gradients was evaluated by estimating trait heritability and the expected responses to selection in natural field conditions.

## 2 | MATERIALS AND METHODS

### 2.1 | Study site, plant material and sampling

The study site Pri Studencu on eutric cambisols is located in the pre-Dinaric Alps of Slovenia (45.7750N, 14.9488 E), 520 m above the sea level (Figure 2), on a slight slope where the highest elevation distance between any two trees was 18.08 m. It is part of a larger natural beech forest, where beech accounts for 68% of the tree biomass, followed by *Picea abies* (20.4%), and with the occasional presence of *Abies alba*, *Acer pseudoplatanus* and *Prunus avium*. However, on the 100 m by 100 m study plot itself, beech accounts for 97.7% of the tree biomass. The beech populations in the vicinity have survived the last glaciation in situ (Brus, 2010; Šercelj, 1996) and are suggested to be the source of genetic variation for other European beech populations, after the ice retreated (Magri, 2008; Magri et al., 2006). With an average yearly temperature of 8.8°C and rainfall of 1449 mm (Slovenian National

<figure><img src="image_5.png" /><figcaption>FIGURE 2 Study plot and sampling (a), and budburst (b) and leaf senescence (c) phenology across years (2017–2019); min, max (first and last date when a tree achieved a certain stage, dotted line) and standard deviation (SD, full line). Letters denote statistical significance after a Kruskal–Wallis test at $\alpha=0.05$ (b and c).</figcaption></figure>
the availability of ovules/pollen is limited (e.g. in isolated populations or under extreme climatic events). In protogynous beech, we expect a combination of stabilizing selection for female/male flowering date around the population average, and overall population directional selection for early flowering, in particular for males.

Lastly, the impact of autumn phenology (leaf senescence) on fitness is poorly understood. However, late leaf senescence would increase autumn photosynthesis (Toomey et al., 2015), and thereby the level of resources available for trees to initiate next-year reproductive buds (Hoch et al., 2013). We, therefore, expect directional selection for late senescence, in order to increase available resources for the following growth and reproductive season.

In this study, we tested these expectations by computing selection gradients based on effective fecundity for a wide variety of reproduction, phenology and growth traits. Moreover, the evolutionary relevance of significant selection gradients was evaluated by estimating trait heritability and the expected responses to selection in natural field conditions.

## 2 | MATERIALS AND METHODS

### 2.1 | Study site, plant material and sampling

The study site Pri Studencu on eutric cambisols is located in the pre-Dinaric Alps of Slovenia (45.7750N, 14.9488 E), 520 m above the sea level (Figure 2), on a slight slope where the highest elevation distance between any two trees was 18.08 m. It is part of a larger natural beech forest, where beech accounts for 68% of the tree biomass, followed by *Picea abies* (20.4%), and with the occasional presence of *Abies alba*, *Acer pseudoplatanus* and *Prunus avium*. However, on the 100 m by 100 m study plot itself, beech accounts for 97.7% of the tree biomass. The beech populations in the vicinity have survived the last glaciation in situ (Brus, 2010; Šercelj, 1996) and are suggested to be the source of genetic variation for other European beech populations, after the ice retreated (Magri, 2008; Magri et al., 2006). With an average yearly temperature of 8.8°C and rainfall of 1449 mm (Slovenian National

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

**FIGURE 2** Study plot and sampling (a), and budburst (b) and leaf senescence (c) phenology across years (2017–2019); min, max (first and last date when a tree achieved a certain stage, dotted line) and standard deviation (SD, full line). Letters denote statistical significance after a Kruskal–Wallis test at $\alpha=0.05$ (b and c).
Meteorological Service), the environmental conditions are optimal for beech growth. Nevertheless, since the beginning of monitoring in 2016, late frost events of varying severity were recorded in 2016, 2017, 2019 and 2020 (see Figure S1).

All adult trees in a one-ha plot were marked, mapped and sampled for genetic analysis in 2016. A tree was considered as adult (i.e. potentially reproductive) if its diameter at 1.3m (DBH) was larger than 15 cm. On the site, there were 251 such trees, estimated to be between 70 and 90 years old (even-aged stand), based on the information from the forestry archives. Five trees had DBH less than 16 cm and did not bear fruit in the 2018 monitoring, but trees with DBH of 16 cm or higher flowered abundantly. In addition to the adult trees, 200 saplings at four randomly chosen locations (50 saplings per location; Figure 2, red circles) were also marked, mapped and sampled for genetic analysis. At the time of sampling in 2016, the saplings, although of similar size, were 5–9 years old (estimated by counting growth rings under microscope on a sample of 15 saplings), most likely belonging to two different regeneration events, a full mast in 2011 and a poorer seed crop in 2007 (information obtained from forestry archives). In addition, 436 seeds from the fructification event in 2016 were sampled from 22 reproducing trees, 20 seeds per tree (except one tree with 16 seeds only). These trees were randomly chosen among the dominant trees in the study plot, that is, trees with their crown at the top of the canopy layer. Finally, in 2019, 200 one-year-old seedlings, a direct product of the 2018 fructification event, located in proximity of the 2016 sapling locations, were also sampled for genetic analysis (Figure 2, blue circles).

To assess the effect of tree competition, the distance-dependent Neighbour-based Competition Index (NCI) was calculated for each of the 251 adult trees. The index assumes that the net effect of neighbouring trees varies as a function of the size of the closest five neighbours (conspecific or not) with DBH larger than 15 cm and as an inverse function of their distance to the target tree (Canham et al., 2004), that is:

$$NCI = \sum_{i=1}^{5} \left( \frac{DBH_{i}}{distance_{i}} \right)$$

For the trees on the edge of the sampling plot, trees outside of the plot were included in NCI calculations.

## 2.2 | Phenotypic traits

Tree size was evaluated by its DBH in 2018. For this study, we regularly monitored all 251 adult trees for phenological traits in 2018 and 2019 using binoculars and a DSLR Camera with a telephoto lens. Observations of budburst phenology were carried out every 6 to 8 days between 20 March and 24 April in 2018 and 27 February and 30 April in 2019. We used scores corresponding to five phenological stages (Table S1; Westergren et al., 2020), from dormant buds (stage 1) to fully unfolded leaves that are smooth and bright (stage 5). We derived the length of budburst (LBB) and the bud phenological

stage at Julian day 100 (BB-JD100) from these observations. We also computed a spring phenology score sum (SPSS) for each tree by summing daily scores during the observation period, following Bontemps et al. (2016). The higher the SPSS, the earlier and quicker the process of budburst (Bontemps et al., 2016). Like budburst, leaf senescence was monitored in weekly intervals, between 23 August and 7 November in 2018 and 28 August and 6 November in 2019. We used a four-stage scale to record the phenology of senescence, from trees with fully green leaves (stage 1) to trees with brown or shed leaves (stage 4) (Table S1; Westergren et al., 2020). We derived the length of leaf senescence (LLS) and the leaf senescence stage at Julian day 287 in 2018 and Julian day 267 in 2019 (LS-JD287 and LS-JD267 respectively) from these observations. The Julian days at which phenological events were considered for further analysis (i.e. BB-JD100, LS-JD287 and LS-JD267) were selected in a way to maximize phenotypic variation. For both budburst and leaf senescence, a tree was assigned to a given phenological stage only when more than two thirds of its crown had reached it. In addition, we recorded reproduction traits in 2018 (there was no flowering and consequently no fructification in 2019). We assessed male flowering on 24 April (FLW) using a three-stage scale and fruit abundance on 18 April (FT-JD108) and 26 July (FT-JD217) using a four-stage scale (Table S1). Finally, the growing season length (GL) was computed as the difference between initiation of budburst (stage 2) and leaf senescence stage 3, when more than two thirds of the leaves on a tree have turned yellow and stopped photosynthesizing. Evaluated traits are summarized in Table 1.

## 2.3 | Molecular markers

Total DNA for seeds, seedlings, saplings and adult trees (N=1051) was isolated from lyophilized leaves and embryos using the DNeasy 96 plant kit (QIAGEN), as per the manufacturer's specifications. All individuals were genotyped at 16 highly polymorphic SSR loci using primers and PCR protocols in Lefevre et al. (2012), but with modified premixes; primer kit (premix) 1 in Lefevre et al. (2012) was divided into two separate premixes in order to avoid potential overlap of allelic ranges of markers sfc_0036 and csolfagus_29; primer concentrations in all three premixes were also optimized. The sizing of the PCR products was performed on the Applied Biosystems 3500 Genetic Analyser with accompanying GeneMapper 5.0 software. One locus, DZ447_A_0, was excluded from further analysis because of presence of null alleles with frequency >0.05 in the 2016 dataset (detected by Genepop v4.6; Rousset, 2008).

## 2.4 | Effective male and female fecundity

We estimated effective male and female fecundity, as a proxy of fitness, from the genotypes and spatial locations of adults and offspring, considering seeds, seedlings and saplings in separate analyses. To that aim, we used the Mixed Effects Mating Model (MEMM)
Meteorological Service), the environmental conditions are optimal for beech growth. Nevertheless, since the beginning of monitoring in 2016, late frost events of varying severity were recorded in 2016, 2017, 2019 and 2020 (see Figure S1).

All adult trees in a one-ha plot were marked, mapped and sampled for genetic analysis in 2016. A tree was considered as adult (i.e. potentially reproductive) if its diameter at 1.3m (DBH) was larger than 15cm. On the site, there were 251 such trees, estimated to be between 70 and 90 years old (even-aged stand), based on the information from the forestry archives. Five trees had DBH less than 16cm and did not bear fruit in the 2018 monitoring, but trees with DBH of 16cm or higher flowered abundantly. In addition to the adult trees, 200 saplings at four randomly chosen locations (50 saplings per location; Figure 2, red circles) were also marked, mapped and sampled for genetic analysis. At the time of sampling in 2016, the saplings, although of similar size, were 5–9 years old (estimated by counting growth rings under microscope on a sample of 15 saplings), most likely belonging to two different regeneration events, a full mast in 2011 and a poorer seed crop in 2007 (information obtained from forestry archives). In addition, 436 seeds from the fructification event in 2016 were sampled from 22 reproducing trees, 20 seeds per tree (except one tree with 16 seeds only). These trees were randomly chosen among the dominant trees in the study plot, that is, trees with their crown at the top of the canopy layer. Finally, in 2019, 200 one-year-old seedlings, a direct product of the 2018 fructification event, located in proximity of the 2016 sapling locations, were also sampled for genetic analysis (Figure 2, blue circles).

To assess the effect of tree competition, the distance-dependent Neighbour-based Competition Index (NCI) was calculated for each of the 251 adult trees. The index assumes that the net effect of neighbouring trees varies as a function of the size of the closest five neighbours (conspecific or not) with DBH larger than 15cm and as an inverse function of their distance to the target tree (Canham et al., 2004), that is:

$$NCI = \sum_{i=1}^{5} \left( \frac{DBH_{i}}{distance_{i}} \right)$$

For the trees on the edge of the sampling plot, trees outside of the plot were included in NCI calculations.

## 2.2 | Phenotypic traits

Tree size was evaluated by its DBH in 2018. For this study, we regularly monitored all 251 adult trees for phenological traits in 2018 and 2019 using binoculars and a DSLR Camera with a telephoto lens. Observations of budburst phenology were carried out every 6 to 8 days between 20 March and 24 April in 2018 and 27 February and 30 April in 2019. We used scores corresponding to five phenological stages (Table S1; Westergren et al., 2020), from dormant buds (stage 1) to fully unfolded leaves that are smooth and bright (stage 5). We derived the length of budburst (LBB) and the bud phenological

stage at Julian day 100 (BB-JD100) from these observations. We also computed a spring phenology score sum (SPSS) for each tree by summing daily scores during the observation period, following Bontemps et al. (2016). The higher the SPSS, the earlier and quicker the process of budburst (Bontemps et al., 2016). Like budburst, leaf senescence was monitored in weekly intervals, between 23 August and 7 November in 2018 and 28 August and 6 November in 2019. We used a four-stage scale to record the phenology of senescence, from trees with fully green leaves (stage 1) to trees with brown or shed leaves (stage 4) (Table S1; Westergren et al., 2020). We derived the length of leaf senescence (LLS) and the leaf senescence stage at Julian day 287 in 2018 and Julian day 267 in 2019 (LS-JD287 and LS-JD267 respectively) from these observations. The Julian days at which phenological events were considered for further analysis (i.e. BB-JD100, LS-JD287 and LS-JD267) were selected in a way to maximize phenotypic variation. For both budburst and leaf senescence, a tree was assigned to a given phenological stage only when more than two thirds of its crown had reached it. In addition, we recorded reproduction traits in 2018 (there was no flowering and consequently no fructification in 2019). We assessed male flowering on 24 April (FLW) using a three-stage scale and fruit abundance on 18 April (FT-JD108) and 26 July (FT-JD217) using a four-stage scale (Table S1). Finally, the growing season length (GL) was computed as the difference between initiation of budburst (stage 2) and leaf senescence stage 3, when more than two thirds of the leaves on a tree have turned yellow and stopped photosynthesizing. Evaluated traits are summarized in Table 1.

## 2.3 | Molecular markers

Total DNA for seeds, seedlings, saplings and adult trees (N=1051) was isolated from lyophilized leaves and embryos using the DNeasy 96 plant kit (QIAGEN), as per the manufacturer's specifications. All individuals were genotyped at 16 highly polymorphic SSR loci using primers and PCR protocols in Lefevre et al. (2012), but with modified premixes; primer kit (premix) 1 in Lefevre et al. (2012) was divided into two separate premixes in order to avoid potential overlap of allelic ranges of markers sfc_0036 and csolfagus_29; primer concentrations in all three premixes were also optimized. The sizing of the PCR products was performed on the Applied Biosystems 3500 Genetic Analyser with accompanying GeneMapper 5.0 software. One locus, DZ447_A_0, was excluded from further analysis because of presence of null alleles with frequency >0.05 in the 2016 dataset (detected by Genepop v4.6; Rousset, 2008).

## 2.4 | Effective male and female fecundity

We estimated effective male and female fecundity, as a proxy of fitness, from the genotypes and spatial locations of adults and offspring, considering seeds, seedlings and saplings in separate analyses. To that aim, we used the Mixed Effects Mating Model (MEMM)
<table><caption>TABLE 1 Loadings of the first three Principal Components (PCs), explaining 72.2% of the variance of phenotypic traits measured in 2018 and 74.2% of the variance of those measured in 2019, after a VARIMAX rotation.</caption><thead><tr><th rowspan="2">Trait</th><th colspan="3">2018</th><th colspan="3">2019</th><th rowspan="2">Description (units)</th></tr><tr><th>PC1 (34.8)</th><th>PC2 (22.2)</th><th>PC3 (15.2)</th><th>PC1 (29.1)</th><th>PC2 (23.5)</th><th>PC3 (21.6)</th></tr></thead><tbody><tr><td><strong>Growth</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>DBH</td><td><strong>0.70</strong></td><td>-0.08</td><td>0.25</td><td>-0.26</td><td><strong>0.69</strong></td><td>-0.12</td><td>Diameter at 1.30 m (cm)</td></tr><tr><td>GL</td><td>-0.45</td><td><strong>0.64</strong></td><td>-0.19</td><td>0.33</td><td>-0.38</td><td><strong>0.76</strong></td><td>Length of the growing season (days)</td></tr><tr><td><strong>Budburst</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>LBB</td><td>0.15</td><td><strong>-0.85</strong></td><td>0.06</td><td>-0.28</td><td>0.15</td><td><strong>-0.62</strong></td><td>Length of budburst (GDD5)</td></tr><tr><td>BB-JD100</td><td><strong>0.61</strong></td><td>0.50</td><td>0.08</td><td><strong>0.93</strong></td><td>-0.11</td><td>0.19</td><td>Bud phenological stage at Julian day 100 (1-5 scale)</td></tr><tr><td>SPSS</td><td>0.33</td><td><strong>0.90</strong></td><td>0.02</td><td><strong>0.94</strong></td><td>-0.17</td><td>0.13</td><td>Spring phenology score sum (see text)</td></tr><tr><td><strong>Reproduction</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FLW</td><td><strong>0.88</strong></td><td>-0.06</td><td>0.11</td><td>NA</td><td>NA</td><td>NA</td><td>Male flowering stage at Julian day 114 (1-3 scale)</td></tr><tr><td>FT-JD108</td><td><strong>0.89</strong></td><td>-0.06</td><td>0.07</td><td>NA</td><td>NA</td><td>NA</td><td>Fruit abundance at Julian day 108 (1-4 scale)</td></tr><tr><td>FT-JD217</td><td><strong>0.73</strong></td><td>0.07</td><td>0.23</td><td>NA</td><td>NA</td><td>NA</td><td>Fruit abundance at Julian day 217 (1-4 scale)</td></tr><tr><td><strong>Leaf senescence</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>LLS</td><td>0.07</td><td>-0.03</td><td><strong>0.91</strong></td><td>-0.21</td><td>0.48</td><td><strong>0.69</strong></td><td>Length of leaf senescence (days)</td></tr><tr><td>LS-JD287/JD267</td><td>0.44</td><td>-0.10</td><td><strong>0.72</strong></td><td>0.00</td><td><strong>0.86</strong></td><td>-0.06</td><td>Leaf senescence stage at Julian day 287 in 2018 or 267 in 2019 (1-4 scale)</td></tr></tbody></table>

Note: For each PC, the % of explained variance is given in brackets. In bold, loadings >0.60. See interpretation of the PCs in the main text. Variable description and measurement units are given for each phenotypic trait. GDD5 stands for Growing Degree Days with base temperature of 5°C.

which jointly infers mating and dispersal parameters for a hermaphroditic plant population in a Bayesian framework (Oddou-Muratorio et al., 2018; code available from the authors). First, we used MEMM to estimate the individual male fecundities jointly with the pollen dispersal kernel, the selfing rate and the pollen migration rate, from the seeds collected on mother-trees in 2016. Second, we used the MEMM seedlings extension of MEMM to estimate jointly male and female effective fecundities with the pollen and seed dispersal kernels, the selfing rate and the pollen and seed migration rates, from established multiyear saplings (collected in 2016) and one-year-old seedlings (collected in 2019) separately.

These statistical models account for the effect of the relative position of putative parents and offspring. At the same basic fecundity, putative parents closer to an offspring would have a higher parentage probability in uncorrected models. Hence, by using MEMM, estimates of fecundity are not sensitive to spatial biases due to sampling design or edge effects. Moreover, MEMM estimates male fecundity as the relative amount of pollen achieving successful pollination, and female fecundity as the relative number of seeds achieving successful germination and establishment in the population. Hence, MEMM estimates of fecundities are 'effective', in the sense that they account for the individual effects (maternal or genetic) that modify the success of mating, including differences in seed maturation,

germination or early survival during the post-dispersal processes. Effective fecundity provides more realistic estimates of individual plant contribution to the next generation than simpler estimates, such as fruit or seed set (Conner et al., 1996; Elle & Meagher, 2000; Oddou-Muratorio et al., 2018). Also, as MEMM estimates of fecundity are relative; they do not convey any information on the total number of pollen grain or seeds produced by a given individual, and they are not affected by demographic processes driving seeds/ seedlings survival independently of their genotypes (e.g. density-dependent mortality). Finally, they take into account uncertainty in parentage reconstruction. Indeed, MEMM does not categorically assign parents to offspring, but rather considers the likelihood of all adults to be the parent of each offspring, accounting for genotyping error.

MEMM models were run for 50,000 iterations of the MCMC with a thinning step of 10 (i.e. 5000 steps recorded) and a burn-in of 500. Genotyping error was set to typical values for nuclear SSRs, that is, 0.01 for $p_{err1}$ (probability of the allele read to differ by only one motif repeat from the true allele) and 0.001 for $p_{err2}$ (probability of the allele read to differ by any number of motif repeats from the true allele), with one mismatch allowed. Runs that allowed a maximum of zero (no genotyping error) to three mismatches produced similar results (not shown).
<table><caption>TABLE 1 Loadings of the first three Principal Components (PCs), explaining 72.2% of the variance of phenotypic traits measured in 2018 and 74.2% of the variance of those measured in 2019, after a VARIMAX rotation.</caption><thead><tr><th rowspan="2">Trait</th><th colspan="3">2018</th><th colspan="3">2019</th><th rowspan="2">Description (units)</th></tr><tr><th>PC1 (34.8)</th><th>PC2 (22.2)</th><th>PC3 (15.2)</th><th>PC1 (29.1)</th><th>PC2 (23.5)</th><th>PC3 (21.6)</th></tr></thead><tbody><tr><td><strong>Growth</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>DBH</td><td><strong>0.70</strong></td><td>-0.08</td><td>0.25</td><td>-0.26</td><td><strong>0.69</strong></td><td>-0.12</td><td>Diameter at 1.30 m (cm)</td></tr><tr><td>GL</td><td>-0.45</td><td><strong>0.64</strong></td><td>-0.19</td><td>0.33</td><td>-0.38</td><td><strong>0.76</strong></td><td>Length of the growing season (days)</td></tr><tr><td><strong>Budburst</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>LBB</td><td>0.15</td><td><strong>-0.85</strong></td><td>0.06</td><td>-0.28</td><td>0.15</td><td><strong>-0.62</strong></td><td>Length of budburst (GDD5)</td></tr><tr><td>BB-JD100</td><td><strong>0.61</strong></td><td>0.50</td><td>0.08</td><td><strong>0.93</strong></td><td>-0.11</td><td>0.19</td><td>Bud phenological stage at Julian day 100 (1-5 scale)</td></tr><tr><td>SPSS</td><td>0.33</td><td><strong>0.90</strong></td><td>0.02</td><td><strong>0.94</strong></td><td>-0.17</td><td>0.13</td><td>Spring phenology score sum (see text)</td></tr><tr><td><strong>Reproduction</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FLW</td><td><strong>0.88</strong></td><td>-0.06</td><td>0.11</td><td>NA</td><td>NA</td><td>NA</td><td>Male flowering stage at Julian day 114 (1-3 scale)</td></tr><tr><td>FT-JD108</td><td><strong>0.89</strong></td><td>-0.06</td><td>0.07</td><td>NA</td><td>NA</td><td>NA</td><td>Fruit abundance at Julian day 108 (1-4 scale)</td></tr><tr><td>FT-JD217</td><td><strong>0.73</strong></td><td>0.07</td><td>0.23</td><td>NA</td><td>NA</td><td>NA</td><td>Fruit abundance at Julian day 217 (1-4 scale)</td></tr><tr><td><strong>Leaf senescence</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>LLS</td><td>0.07</td><td>-0.03</td><td><strong>0.91</strong></td><td>-0.21</td><td>0.48</td><td><strong>0.69</strong></td><td>Length of leaf senescence (days)</td></tr><tr><td>LS-JD287/JD267</td><td>0.44</td><td>-0.10</td><td><strong>0.72</strong></td><td>0.00</td><td><strong>0.86</strong></td><td>-0.06</td><td>Leaf senescence stage at Julian day 287 in 2018 or 267 in 2019 (1-4 scale)</td></tr></tbody></table>

Note: For each PC, the % of explained variance is given in brackets. In bold, loadings >0.60. See interpretation of the PCs in the main text. Variable description and measurement units are given for each phenotypic trait. GDD5 stands for Growing Degree Days with base temperature of 5°C.

which jointly infers mating and dispersal parameters for a hermaphroditic plant population in a Bayesian framework (Oddou-Muratorio et al., 2018; code available from the authors). First, we used MEMM to estimate the individual male fecundities jointly with the pollen dispersal kernel, the selfing rate and the pollen migration rate, from the seeds collected on mother-trees in 2016. Second, we used the MEMM seedlings extension of MEMM to estimate jointly male and female effective fecundities with the pollen and seed dispersal kernels, the selfing rate and the pollen and seed migration rates, from established multiyear saplings (collected in 2016) and one-year-old seedlings (collected in 2019) separately.

These statistical models account for the effect of the relative position of putative parents and offspring. At the same basic fecundity, putative parents closer to an offspring would have a higher parentage probability in uncorrected models. Hence, by using MEMM, estimates of fecundity are not sensitive to spatial biases due to sampling design or edge effects. Moreover, MEMM estimates male fecundity as the relative amount of pollen achieving successful pollination, and female fecundity as the relative number of seeds achieving successful germination and establishment in the population. Hence, MEMM estimates of fecundities are 'effective', in the sense that they account for the individual effects (maternal or genetic) that modify the success of mating, including differences in seed maturation,

germination or early survival during the post-dispersal processes. Effective fecundity provides more realistic estimates of individual plant contribution to the next generation than simpler estimates, such as fruit or seed set (Conner et al., 1996; Elle & Meagher, 2000; Oddou-Muratorio et al., 2018). Also, as MEMM estimates of fecundity are relative; they do not convey any information on the total number of pollen grain or seeds produced by a given individual, and they are not affected by demographic processes driving seeds/ seedlings survival independently of their genotypes (e.g. density-dependent mortality). Finally, they take into account uncertainty in parentage reconstruction. Indeed, MEMM does not categorically assign parents to offspring, but rather considers the likelihood of all adults to be the parent of each offspring, accounting for genotyping error.

MEMM models were run for 50,000 iterations of the MCMC with a thinning step of 10 (i.e. 5000 steps recorded) and a burn-in of 500. Genotyping error was set to typical values for nuclear SSRs, that is, 0.01 for $p_{err1}$ (probability of the allele read to differ by only one motif repeat from the true allele) and 0.001 for $p_{err2}$ (probability of the allele read to differ by any number of motif repeats from the true allele), with one mismatch allowed. Runs that allowed a maximum of zero (no genotyping error) to three mismatches produced similar results (not shown).
## 2.5 | Single-trait and multivariate selection gradients

We computed linear ($\beta_p$) and quadratic ($\gamma_q$) SD-standardized selection gradients (with 'SD' standing for 'standard deviation', noted here also as $\sigma_p$; see below) using effective fecundities estimated from seeds collected in 2016 from mother trees and 1-year-old seedlings collected in 2019 in the field (which resulted from seeds sired in 2018) and standardized trait values obtained for the 2018 growing season, following Lande and Arnold (1983). Similarly, SD-standardized selection gradients were computed using effective fecundities estimated from a multiyear sapling cohort collected in 2016 and phenotypic data obtained for both 2018 and 2019 growing seasons. While selection gradients based on seeds and 1-year-old seedlings should mostly reflect the processes taking place at early stages of selection and be strongly related to phenotypic traits measured the previous year, those based on multiyear saplings should provide a better view of the established cohort that will form the next generation and be affected by mother phenotypes associated with multiple regeneration events (i.e. years).

First, univariate SD-standardized selection gradients were computed for the male function using individual effective male fecundity estimates from the MEMM model (seeds), and for both the male and female functions using effective fecundity estimates from the extended MEMM seedlings model (seedlings and saplings), and considering both linear and quadratic components. Effective fecundity estimates were log-transformed prior to selection gradient analysis as in Oddou-Muratorio et al. (2018). Models with non-significant quadratic components were rerun considering only linear selection gradients. As beech has both male and female flowers on the same individual, an individual selection gradient can be computed that considers both female and male functions. Hence, univariate SD-standardized selection gradients for total (i.e. individual) effective fecundity were computed averaging the male and female fecundity obtained from the extended MEMM seedlings model based on saplings, and used for estimation of the overall expected selection responses in the field (see below).

Second, multivariate SD-standardized selection gradients, which better account for correlations between traits, were computed for each of the study years separately (2018, 2019) based on standardized trait values resulting from a Principal Component Analysis (PCA) after VARIMAX rotation for the traits measured each year. The first three Principal Components (PCs), explaining 72.2% of the variance in 2018 and 74.2% of the variance in 2019, were selected based on their eigenvalues (i.e. eigenvalue > 1; Figure S4A,B). Then, the following multivariate equation was adjusted (Alía et al., 2014; Lande & Arnold, 1983) using the STATS package in R (R Core Team, 2020):

$$w = \alpha + \sum_{i=1}^{3} \beta_{\sigma_i} z_i + \frac{1}{2} \sum_{i=1}^{3} \gamma_{\sigma_i} z_i^2 + \varepsilon,$$

where w is the relative fitness estimated as the log-transformed effective fecundity, $\alpha$ is the population mean, $i=1$ to 3 are the first three

PCs, and $\beta_{\sigma_i}$ and $\gamma_{\sigma_i}$ are the SD-standardized linear and quadratic selection gradients respectively.

## 2.6 | Heritability, evolvability and expected selection response

To obtain a pedigree for calculation of the additive genetic variance (see below), the genealogical relationships among adult trees were reconstructed following Fernández and Toro (2006), as implemented in MOL_COANC vs. 2 (available from the authors). This method incorporates selfing (although selfing was not significantly different from zero in our population, see Table S2) and uses simulated annealing to identify the pedigree that best matches the coancestry matrix estimated with molecular markers, in our case 15 nuSSRs. Runs assuming 1, 3, 5, 8 and 20 previous generations (i.e. pedigree depths) were performed with a maximum number of allowed step ('temperatures') set at 400, the number of solutions tested in each step fixed to 5000 and a rate of decrease in temperature in each step of 0.9. The run with the highest correlation between the observed (based on molecular markers) and the estimated (based on genealogical relationships) coancestry matrix was kept for further analyses (i.e. the one for three previous generations, with correlation of 0.769; called below 'reconstructed pedigree').

Quantitative genetic parameters were estimated by fitting animal models in MCMCglmm R package (de Villemereuil et al., 2012; Hadfield, 2010). For Gaussian traits, the model was built as follows:

$$x = \mu + Y_1c + Y_2e + Za + \varepsilon,$$

where $x$ is the vector of observed phenotypic data (after normalization using the 'qqnorm' function in R), $\mu$ is the population mean phenotypic value, $c$ is the vector of the fixed effect of competition (NCI), $e$ is the vector of the fixed effect of tree size (DBH), $a$ is the vector of the random genetic (additive) effect of $x$, $\varepsilon$ is the vector of the random effect of the residuals, and $Y_1$, $Y_2$ and $Z$ are the index matrices related to each effect. The main parameter we want to estimate with these models is the additive genetic variance $\sigma_A^2$,

$$a \sim N(0; A\sigma_A^2),$$

where $A$ is the additive genetic relationship matrix computed from the reconstructed pedigree. For non-Gaussian traits (i.e. those based on an interval scale, Table 1), the model was extended by using the 'ordinal' distribution family with a probit link (de Villemereuil et al., 2012; Hadfield, 2010).

Burn-in was set to 150,000 and chain length to one to four million iterations, depending on the model. NCI and DBH were used as covariates to reduce environmental noise, which is expected to be relatively high in field conditions (compared to common gardens). Models with two, one or no covariates were run and those with smaller Deviance Information Criterion (DIC) were selected. Chain convergence was evaluated using the Heidelberg stationarity test and visual observation of the trace plots as in de Villemereuil
## 2.5 | Single-trait and multivariate selection gradients

We computed linear ($\beta_p$) and quadratic ($\gamma_q$) SD-standardized selection gradients (with 'SD' standing for 'standard deviation', noted here also as $\sigma_p$; see below) using effective fecundities estimated from seeds collected in 2016 from mother trees and 1-year-old seedlings collected in 2019 in the field (which resulted from seeds sired in 2018) and standardized trait values obtained for the 2018 growing season, following Lande and Arnold (1983). Similarly, SD-standardized selection gradients were computed using effective fecundities estimated from a multiyear sapling cohort collected in 2016 and phenotypic data obtained for both 2018 and 2019 growing seasons. While selection gradients based on seeds and 1-year-old seedlings should mostly reflect the processes taking place at early stages of selection and be strongly related to phenotypic traits measured the previous year, those based on multiyear saplings should provide a better view of the established cohort that will form the next generation and be affected by mother phenotypes associated with multiple regeneration events (i.e. years).

First, univariate SD-standardized selection gradients were computed for the male function using individual effective male fecundity estimates from the MEMM model (seeds), and for both the male and female functions using effective fecundity estimates from the extended MEMM seedlings model (seedlings and saplings), and considering both linear and quadratic components. Effective fecundity estimates were log-transformed prior to selection gradient analysis as in Oddou-Muratorio et al. (2018). Models with non-significant quadratic components were rerun considering only linear selection gradients. As beech has both male and female flowers on the same individual, an individual selection gradient can be computed that considers both female and male functions. Hence, univariate SD-standardized selection gradients for total (i.e. individual) effective fecundity were computed averaging the male and female fecundity obtained from the extended MEMM seedlings model based on saplings, and used for estimation of the overall expected selection responses in the field (see below).

Second, multivariate SD-standardized selection gradients, which better account for correlations between traits, were computed for each of the study years separately (2018, 2019) based on standardized trait values resulting from a Principal Component Analysis (PCA) after VARIMAX rotation for the traits measured each year. The first three Principal Components (PCs), explaining 72.2% of the variance in 2018 and 74.2% of the variance in 2019, were selected based on their eigenvalues (i.e. eigenvalue > 1; Figure S4A,B). Then, the following multivariate equation was adjusted (Alía et al., 2014; Lande & Arnold, 1983) using the STATS package in R (R Core Team, 2020):

$$w = \alpha + \sum_{i=1}^{3} \beta_{\sigma_i} z_i + \frac{1}{2} \sum_{i=1}^{3} \gamma_{\sigma_i} z_i^2 + \varepsilon,$$

where w is the relative fitness estimated as the log-transformed effective fecundity, $\alpha$ is the population mean, $i=1$ to 3 are the first three

PCs, and $\beta_{\sigma_i}$ and $\gamma_{\sigma_i}$ are the SD-standardized linear and quadratic selection gradients respectively.

## 2.6 | Heritability, evolvability and expected selection response

To obtain a pedigree for calculation of the additive genetic variance (see below), the genealogical relationships among adult trees were reconstructed following Fernández and Toro (2006), as implemented in MOL_COANC vs. 2 (available from the authors). This method incorporates selfing (although selfing was not significantly different from zero in our population, see Table S2) and uses simulated annealing to identify the pedigree that best matches the coancestry matrix estimated with molecular markers, in our case 15 nuSSRs. Runs assuming 1, 3, 5, 8 and 20 previous generations (i.e. pedigree depths) were performed with a maximum number of allowed step ('temperatures') set at 400, the number of solutions tested in each step fixed to 5000 and a rate of decrease in temperature in each step of 0.9. The run with the highest correlation between the observed (based on molecular markers) and the estimated (based on genealogical relationships) coancestry matrix was kept for further analyses (i.e. the one for three previous generations, with correlation of 0.769; called below 'reconstructed pedigree').

Quantitative genetic parameters were estimated by fitting animal models in MCMCglmm R package (de Villemereuil et al., 2012; Hadfield, 2010). For Gaussian traits, the model was built as follows:

$$x = \mu + Y_1c + Y_2e + Za + \varepsilon,$$

where $x$ is the vector of observed phenotypic data (after normalization using the 'qqnorm' function in R), $\mu$ is the population mean phenotypic value, $c$ is the vector of the fixed effect of competition (NCI), $e$ is the vector of the fixed effect of tree size (DBH), $a$ is the vector of the random genetic (additive) effect of $x$, $\varepsilon$ is the vector of the random effect of the residuals, and $Y_1$, $Y_2$ and $Z$ are the index matrices related to each effect. The main parameter we want to estimate with these models is the additive genetic variance $\sigma_A^2$,

$$a \sim N(0; A\sigma_A^2),$$

where $A$ is the additive genetic relationship matrix computed from the reconstructed pedigree. For non-Gaussian traits (i.e. those based on an interval scale, Table 1), the model was extended by using the 'ordinal' distribution family with a probit link (de Villemereuil et al., 2012; Hadfield, 2010).

Burn-in was set to 150,000 and chain length to one to four million iterations, depending on the model. NCI and DBH were used as covariates to reduce environmental noise, which is expected to be relatively high in field conditions (compared to common gardens). Models with two, one or no covariates were run and those with smaller Deviance Information Criterion (DIC) were selected. Chain convergence was evaluated using the Heidelberg stationarity test and visual observation of the trace plots as in de Villemereuil
et al. (2012). MCMC diagnostics (i.e. autocorrelation and effective sample size) were assessed to ensure algorithm efficiency.

Estimates of variance components for total effective fecundity (i.e. our proxy for fitness) and each trait (mode), obtained from best MCMCglmm runs, were used to compute narrow-sense heritability ($h^2$), as follows:

$$h^2 = \frac{\sigma_A^2}{\sigma_A^2 + \sigma_R^2}$$

where $\sigma_A^2$ is the additive genetic variance and $\sigma_R^2$ the residual variance. For non-Gaussian traits (Table 1), it is necessary to account for an additional source of variance coming from the probit link (Nakagawa & Schielzeth, 2010), resulting in:

$$h^2 = \frac{\sigma_A^2}{\sigma_A^2 + \sigma_R^2 + 1}$$

Evolvability ($I_A$) is often preferable to heritability as a measure of adaptive capacity, because it allows a more direct interpretation as the expected percent change in a trait under a unit strength of selection and it is not affected by microenvironmental variation within populations (Hansen et al., 2011). Evolvability is meaningless for traits measured in an interval scale and thus it was only computed for Gaussian traits, as follows:

$$I_A = \left( \frac{\sigma_A}{\mu} \right)^2$$

Finally, the expected selection response ($R$), that is, the expected genetic change in a phenotypic trait value between generations (in the same unit as the measured trait), was computed as:

$$R = h^2 \beta_a \sigma_p,$$

where $h^2$ is the narrow-sense heritability (as computed above), $\beta_a$ is the overall SD-standardized selection gradient based on established multiyear saplings for total effective fecundity, and $\sigma_p$ the phenotypic standard deviation.

# 3 | RESULTS

## 3.1 | Male and female effective fecundity based on seedlings/saplings

Effective fecundity of individual trees varied among years and sexes (Figure S2) and was highly skewed: many trees had a negligible female/male fecundity while some trees had a large fecundity. This skewness was more marked in 2016 when estimated from established saplings (5–9 years old) than in 2019 when estimated from one-year-old seedlings. Female fecundity was higher than one in 41 trees (~16%) in 2016 and 44 (~18%) in 2019, while male fecundity was higher than one in 36 trees (~14%) in 2016 and 53 trees (~21%) in 2019.

## 3.2 | Phenotypic trait variation and correlation

Average trait values and their standard deviation (SD) for reproductive and growth phenology and growth (size) measured in 2018 and 2019 are given in Figure 2 and Table 3.

Concerning traits measured in 2018, larger trees had earlier reproductive phenology and, to a lesser extent earlier growth phenology (budburst and leaf senescence), based on pairwise trait correlation (see Pearson's correlation plot in Figure S3A). Large trees with early budburst (high BB-JD100) produced a higher number of flowers and fruits and had slightly shorter length of the growing season. These correlations are mostly captured by PC1$_{2018}$ (see PC loadings in Table 1; see also Figure S4C–F for the PCA based on 2018 traits). In spring, trees with early budburst also had a faster budburst (low LBB). These two phenology variables were well-captured by SPSS, which has a heavy loading on PC2$_{2018}$. In contrast, early leaf senescence (high LS-JD287) was associated with longer duration of leaf senescence, as suggested by Pearson correlations (Figure S3A) and reflected in PC3$_{2018}$.

Correlations between budburst traits were weak across years and only significant for SPSS (see Figure S3C and Figure 2b). Correlations across years were higher for autumn leaf senescence traits, although interannual differences were still significant (see Figure S3C and Figure 2c). Given this variation in growth phenology across years, larger trees had later budburst in 2019 than in 2018 but maintained early leaf senescence. As a result, they had a shorter growing season in 2019 than in 2018 (see Figure S3B). Accordingly, PC1$_{2019}$ is mostly explained by budburst traits but not by tree size, while PC2$_{2019}$ is positively related to large trees with early autumn leaf senescence (see loading of PCs in Table 1; see also Figure S4G–J). For both 2018 and 2019 measurements, the longer the duration of the budburst and the earlier the leaf senescence, the shorter the duration of the growing season. The length of the growing season was mostly captured by PC2$_{2018}$ and PC3$_{2019}$ (Table 1; Figure S4).

## 3.3 | Single-trait and multivariate selection gradients

Overall, single-trait and multivariate SD-standardized selection gradients showed similar patterns. However, selection gradients based on seed and one-year-old seedling models (Table 2) often showed notable differences from the models for established saplings, which represent a multiyear cohort, despite some common trends (e.g. for male flowering phenology; Table 3). Trees with early budburst (high BB-JD100) sired more seeds and fathered more one-year-old seedlings the year after the reproduction event took place (as shown by a significant positive directional selection gradient; Table 2) but did not contribute a higher number of established saplings overall (non-significant selection gradient in Table 3). In addition, although early leaf senescence was associated with a higher male fecundity based on models for established multiyear saplings (see below; Table 3 and Figure 3d), it had no effect on the number of sired seeds or
et al. (2012). MCMC diagnostics (i.e. autocorrelation and effective sample size) were assessed to ensure algorithm efficiency.

Estimates of variance components for total effective fecundity (i.e. our proxy for fitness) and each trait (mode), obtained from best MCMCglmm runs, were used to compute narrow-sense heritability ($h^2$), as follows:

$$h^2 = \frac{\sigma_A^2}{\sigma_A^2 + \sigma_R^2}$$

where $\sigma_A^2$ is the additive genetic variance and $\sigma_R^2$ the residual variance. For non-Gaussian traits (Table 1), it is necessary to account for an additional source of variance coming from the probit link (Nakagawa & Schielzeth, 2010), resulting in:

$$h^2 = \frac{\sigma_A^2}{\sigma_A^2 + \sigma_R^2 + 1}$$

Evolvability ($I_A$) is often preferable to heritability as a measure of adaptive capacity, because it allows a more direct interpretation as the expected percent change in a trait under a unit strength of selection and it is not affected by microenvironmental variation within populations (Hansen et al., 2011). Evolvability is meaningless for traits measured in an interval scale and thus it was only computed for Gaussian traits, as follows:

$$I_A = \left( \frac{\sigma_A}{\mu} \right)^2$$

Finally, the expected selection response ($R$), that is, the expected genetic change in a phenotypic trait value between generations (in the same unit as the measured trait), was computed as:

$$R = h^2 \beta_a \sigma_p,$$

where $h^2$ is the narrow-sense heritability (as computed above), $\beta_a$ is the overall SD-standardized selection gradient based on established multiyear saplings for total effective fecundity, and $\sigma_p$ the phenotypic standard deviation.

# 3 | RESULTS

## 3.1 | Male and female effective fecundity based on seedlings/saplings

Effective fecundity of individual trees varied among years and sexes (Figure S2) and was highly skewed: many trees had a negligible female/male fecundity while some trees had a large fecundity. This skewness was more marked in 2016 when estimated from established saplings (5–9 years old) than in 2019 when estimated from one-year-old seedlings. Female fecundity was higher than one in 41 trees (~16%) in 2016 and 44 (~18%) in 2019, while male fecundity was higher than one in 36 trees (~14%) in 2016 and 53 trees (~21%) in 2019.

## 3.2 | Phenotypic trait variation and correlation

Average trait values and their standard deviation (SD) for reproductive and growth phenology and growth (size) measured in 2018 and 2019 are given in Figure 2 and Table 3.

Concerning traits measured in 2018, larger trees had earlier reproductive phenology and, to a lesser extent earlier growth phenology (budburst and leaf senescence), based on pairwise trait correlation (see Pearson's correlation plot in Figure S3A). Large trees with early budburst (high BB-JD100) produced a higher number of flowers and fruits and had slightly shorter length of the growing season. These correlations are mostly captured by PC1$_{2018}$ (see PC loadings in Table 1; see also Figure S4C–F for the PCA based on 2018 traits). In spring, trees with early budburst also had a faster budburst (low LBB). These two phenology variables were well-captured by SPSS, which has a heavy loading on PC2$_{2018}$. In contrast, early leaf senescence (high LS-JD287) was associated with longer duration of leaf senescence, as suggested by Pearson correlations (Figure S3A) and reflected in PC3$_{2018}$.

Correlations between budburst traits were weak across years and only significant for SPSS (see Figure S3C and Figure 2b). Correlations across years were higher for autumn leaf senescence traits, although interannual differences were still significant (see Figure S3C and Figure 2c). Given this variation in growth phenology across years, larger trees had later budburst in 2019 than in 2018 but maintained early leaf senescence. As a result, they had a shorter growing season in 2019 than in 2018 (see Figure S3B). Accordingly, PC1$_{2019}$ is mostly explained by budburst traits but not by tree size, while PC2$_{2019}$ is positively related to large trees with early autumn leaf senescence (see loading of PCs in Table 1; see also Figure S4G–J). For both 2018 and 2019 measurements, the longer the duration of the budburst and the earlier the leaf senescence, the shorter the duration of the growing season. The length of the growing season was mostly captured by PC2$_{2018}$ and PC3$_{2019}$ (Table 1; Figure S4).

## 3.3 | Single-trait and multivariate selection gradients

Overall, single-trait and multivariate SD-standardized selection gradients showed similar patterns. However, selection gradients based on seed and one-year-old seedling models (Table 2) often showed notable differences from the models for established saplings, which represent a multiyear cohort, despite some common trends (e.g. for male flowering phenology; Table 3). Trees with early budburst (high BB-JD100) sired more seeds and fathered more one-year-old seedlings the year after the reproduction event took place (as shown by a significant positive directional selection gradient; Table 2) but did not contribute a higher number of established saplings overall (non-significant selection gradient in Table 3). In addition, although early leaf senescence was associated with a higher male fecundity based on models for established multiyear saplings (see below; Table 3 and Figure 3d), it had no effect on the number of sired seeds or
<table><caption>TABLE 2 Univariate SD-standardized selection gradients (linear models) based on standardized trait values measured in 2018 and fitness estimates based on seeds (only male selection gradients) and 1-year-old seedlings sired in 2018 (both male and female selection gradients).</caption><thead><tr><th rowspan="3">Trait</th><th rowspan="3">N</th><th colspan="4">Female selection gradients (β<sub>fα</sub>)</th><th colspan="2">Male selection gradients (β<sub>mα</sub>)</th></tr><tr><th colspan="2">Seedling model</th><th colspan="2">Seed model</th><th colspan="2">Seedling model</th></tr><tr><th>Linear</th><th>Quadratic</th><th>Linear</th><th>Quadratic</th><th>Linear</th><th>Quadratic</th></tr></thead><tbody><tr><td>Growth</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>DBH</td><td>251</td><td><strong>0.088***</strong></td><td>-0.014<sup>ns</sup></td><td><strong>0.179***</strong></td><td>0.066<sup>ns</sup></td><td><strong>0.069***</strong></td><td>0.000<sup>ns</sup></td></tr><tr><td>GL</td><td>251</td><td>-0.018<sup>ns</sup></td><td>0.016<sup>ns</sup></td><td>-0.037<sup>ns</sup></td><td>-0.038<sup>ns</sup></td><td>-0.003<sup>ns</sup></td><td>0.002<sup>ns</sup></td></tr><tr><td>Budburst</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>LBB</td><td>251</td><td>0.013<sup>ns</sup></td><td>-0.140<sup>ns</sup></td><td><strong>0.057*</strong></td><td>0.016<sup>ns</sup></td><td>0.002<sup>ns</sup></td><td>-0.178<sup>ns</sup></td></tr><tr><td>BB-JD100</td><td>251</td><td>0.040<sup>ns</sup></td><td>0.002<sup>ns</sup></td><td><strong>0.063*</strong></td><td>-0.074<sup>ns</sup></td><td><strong>0.043**</strong></td><td><strong>-0.046*</strong></td></tr><tr><td>SPSS</td><td>251</td><td>0.019<sup>ns</sup></td><td>0.006<sup>ns</sup></td><td>0.000<sup>ns</sup></td><td>0.006<sup>ns</sup></td><td>0.024<sup>ns</sup></td><td>-0.006<sup>ns</sup></td></tr><tr><td>Reproduction</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FLW</td><td>251</td><td><strong>0.074**</strong></td><td>0.020<sup>ns</sup></td><td><strong>0.105***</strong></td><td>-0.084<sup>ns</sup></td><td><strong>0.047***</strong></td><td>-0.066<sup>ns</sup></td></tr><tr><td>FT-JD108</td><td>251</td><td>-0.011<sup>ns</sup></td><td><strong>-0.188*</strong></td><td>-0.009<sup>ns</sup></td><td><strong>-0.264**</strong></td><td><strong>0.046***</strong></td><td>-0.064<sup>ns</sup></td></tr><tr><td>FT-JD217</td><td>251</td><td>0.042<sup>ns</sup></td><td><strong>-0.106*</strong></td><td>0.007<sup>ns</sup></td><td><strong>-0.180**</strong></td><td><strong>0.046***</strong></td><td>-0.026<sup>ns</sup></td></tr><tr><td>Leaf senescence</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>LLS</td><td>239</td><td>0.013<sup>ns</sup></td><td>-0.010<sup>ns</sup></td><td>0.013<sup>ns</sup></td><td>0.050<sup>ns</sup></td><td>0.008<sup>ns</sup></td><td>-0.026<sup>ns</sup></td></tr><tr><td>LS-JD287</td><td>251</td><td>0.028<sup>ns</sup></td><td>-0.044<sup>ns</sup></td><td>0.020<sup>ns</sup></td><td>-0.094<sup>ns</sup></td><td>0.016<sup>ns</sup></td><td>-0.028<sup>ns</sup></td></tr></tbody></table>

Note: Fitness estimates based on fecundity were log-transformed following Oddou-Muratorio et al. (2018). Models with non-significant quadratic component were rerun considering only linear selection gradients. GDD5 stands for Growing Degree Days with base temperature of 5°C. Significant univariate selection gradients are in bold. *α < 0.05; **α < 0.01; ***α < 0.001; ns, not significant.

Abbreviations: BB-JD100, Bud phenological stage at Julian day 100; DBH, Diameter at 1.30 m; FLW, Male flowering stage at Julian day 114; FT-JD108/JD217, Fruit abundance at Julian day 108/217; GL, Length of the growing season; LBB, Length of budburst; LLS, Length of leaf senescence; LS-JD287, Leaf senescence stage at Julian day 287; SPSS, Spring phenology score sum.

one-year-old seedlings (Table 2). Differences in selection gradients for seeds/one-year-old seedlings versus multiyear saplings suggest that post-dispersal processes and mortality during early establishment are important factors filtering the genotypes that will form the next adult cohort in the studied population. They also hint that different selective processes may affect single reproductive events. Since only seeds and seedlings that become part of the next generation of breeders (i.e. reproductive trees) are relevant for evolution, we will henceforth focus on selection gradients based on estimates of effective fecundity obtained for established multiyear saplings. These saplings, which have survived the initial steps of selection, will form the core of the next generation of the forest and, therefore, the next generation of reproductive trees. We acknowledge that selection gradients were estimated based on 2 years of phenotypic measurements, which are not the years during which the multiyear saplings were produced. However, using established saplings seems more appropriate for investigating which phenotypic traits in parent trees are associated with a greater number of offspring among future reproductive trees.

For traits measured in 2018, larger trees with early male flowering and more abundant fructification had higher effective fecundity based on established multiyear saplings, consistent with a significant directional selection gradient (Table 3; Figure 3a; Table 5). Interestingly, for both female and male functions, there was also a

significant negative quadratic selection gradient for fruit abundance at the end of summer (FT-JD217) indicating stabilizing selection (Table 3a; Figure S5A,B). Unfortunately, because of masting, the studied population did not produce fruit in 2019 and we could not validate this result in a second year.

Based on models for established multiyear saplings, selection gradients of phenological traits were slightly different across years, in particular for budburst and female function. Considering 2018's growth phenology, female effective fecundity increased with earlier budburst, whereas male effective fecundity increased with early leaf senescence, as shown by the significant positive female directional selection gradients for BB-JD100 and SPSS, and male directional selection gradients for LLS and LS-JD287 (Table 3a; also partially reflected in Figure 3a based on PCs and total effective fecundity). For 2019's growth phenology, however, not only mothers with late budburst, as in 2018, but also those with very early budburst were selected against, resulting in overall negative directional and quadratic selection gradients for BB-JD100 that indicate stabilizing selection (Table 3b, Figure S5C,D). Multivariate selection gradients also supported stabilizing selection for female budburst (Figure 3c; Table 5). For both 2018's and 2019's autumn phenology, early leaf senescence was associated with higher male fecundity (Figure 3a,d; Table 5). Lastly, both single-trait (for traits measured in 2018 and 2019) and multivariate selection gradients based on total effective
<table><caption>TABLE 2 Univariate SD-standardized selection gradients (linear models) based on standardized trait values measured in 2018 and fitness estimates based on seeds (only male selection gradients) and 1-year-old seedlings sired in 2018 (both male and female selection gradients).</caption><thead><tr><th rowspan="3">Trait</th><th rowspan="3">N</th><th colspan="4">Female selection gradients (β<sub>fα</sub>)</th><th colspan="2">Male selection gradients (β<sub>me</sub>)</th></tr><tr><th colspan="2">Seedling model</th><th colspan="2">Seed model</th><th colspan="2">Seedling model</th></tr><tr><th>Linear</th><th>Quadratic</th><th>Linear</th><th>Quadratic</th><th>Linear</th><th>Quadratic</th></tr></thead><tbody><tr><td>Growth</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>DBH</td><td>251</td><td><strong>0.088</strong>***</td><td>-0.014<sup>ns</sup></td><td><strong>0.179</strong>***</td><td>0.066<sup>ns</sup></td><td><strong>0.069</strong>***</td><td>0.000<sup>ns</sup></td></tr><tr><td>GL</td><td>251</td><td>-0.018<sup>ns</sup></td><td>0.016<sup>ns</sup></td><td>-0.037<sup>ns</sup></td><td>-0.038<sup>ns</sup></td><td>-0.003<sup>ns</sup></td><td>0.002<sup>ns</sup></td></tr><tr><td>Budburst</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>LBB</td><td>251</td><td>0.013<sup>ns</sup></td><td>-0.140<sup>ns</sup></td><td><strong>0.057</strong>*</td><td>0.016<sup>ns</sup></td><td>0.002<sup>ns</sup></td><td>-0.178<sup>ns</sup></td></tr><tr><td>BB-JD100</td><td>251</td><td>0.040<sup>ns</sup></td><td>0.002<sup>ns</sup></td><td><strong>0.063</strong>*</td><td>-0.074<sup>ns</sup></td><td><strong>0.043</strong>**</td><td><strong>-0.046</strong>*</td></tr><tr><td>SPSS</td><td>251</td><td>0.019<sup>ns</sup></td><td>0.006<sup>ns</sup></td><td>0.000<sup>ns</sup></td><td>0.006<sup>ns</sup></td><td>0.024<sup>ns</sup></td><td>-0.006<sup>ns</sup></td></tr><tr><td>Reproduction</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>FLW</td><td>251</td><td><strong>0.074</strong>**</td><td>0.020<sup>ns</sup></td><td><strong>0.105</strong>***</td><td>-0.084<sup>ns</sup></td><td><strong>0.047</strong>***</td><td>-0.066<sup>ns</sup></td></tr><tr><td>FT-JD108</td><td>251</td><td>-0.011<sup>ns</sup></td><td><strong>-0.188</strong>*</td><td>-0.009<sup>ns</sup></td><td><strong>-0.264</strong>**</td><td><strong>0.046</strong>***</td><td>-0.064<sup>ns</sup></td></tr><tr><td>FT-JD217</td><td>251</td><td>0.042<sup>ns</sup></td><td><strong>-0.106</strong>*</td><td>0.007<sup>ns</sup></td><td><strong>-0.180</strong>**</td><td><strong>0.046</strong>***</td><td>-0.026<sup>ns</sup></td></tr><tr><td>Leaf senescence</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>LLS</td><td>239</td><td>0.013<sup>ns</sup></td><td>-0.010<sup>ns</sup></td><td>0.013<sup>ns</sup></td><td>0.050<sup>ns</sup></td><td>0.008<sup>ns</sup></td><td>-0.026<sup>ns</sup></td></tr><tr><td>LS-JD287</td><td>251</td><td>0.028<sup>ns</sup></td><td>-0.044<sup>ns</sup></td><td>0.020<sup>ns</sup></td><td>-0.094<sup>ns</sup></td><td>0.016<sup>ns</sup></td><td>-0.028<sup>ns</sup></td></tr></tbody></table>

Note: Fitness estimates based on fecundity were log-transformed following Oddou-Muratorio et al. (2018). Models with non-significant quadratic component were rerun considering only linear selection gradients. GDD5 stands for Growing Degree Days with base temperature of 5°C. Significant univariate selection gradients are in bold. *α < 0.05; **α < 0.01; ***α < 0.001; ns, not significant.

Abbreviations: BB-JD100, Bud phenological stage at Julian day 100; DBH, Diameter at 1.30 m; FLW, Male flowering stage at Julian day 114; FT-JD108/JD217, Fruit abundance at Julian day 108/217; GL, Length of the growing season; LBB, Length of budburst; LLS, Length of leaf senescence; LS-JD287, Leaf senescence stage at Julian day 287; SPSS, Spring phenology score sum.

one-year-old seedlings (Table 2). Differences in selection gradients for seeds/one-year-old seedlings versus multiyear saplings suggest that post-dispersal processes and mortality during early establishment are important factors filtering the genotypes that will form the next adult cohort in the studied population. They also hint that different selective processes may affect single reproductive events. Since only seeds and seedlings that become part of the next generation of breeders (i.e. reproductive trees) are relevant for evolution, we will henceforth focus on selection gradients based on estimates of effective fecundity obtained for established multiyear saplings. These saplings, which have survived the initial steps of selection, will form the core of the next generation of the forest and, therefore, the next generation of reproductive trees. We acknowledge that selection gradients were estimated based on 2 years of phenotypic measurements, which are not the years during which the multiyear saplings were produced. However, using established saplings seems more appropriate for investigating which phenotypic traits in parent trees are associated with a greater number of offspring among future reproductive trees.

For traits measured in 2018, larger trees with early male flowering and more abundant fructification had higher effective fecundity based on established multiyear saplings, consistent with a significant directional selection gradient (Table 3; Figure 3a; Table 5). Interestingly, for both female and male functions, there was also a

significant negative quadratic selection gradient for fruit abundance at the end of summer (FT-JD217) indicating stabilizing selection (Table 3a; Figure S5A,B). Unfortunately, because of masting, the studied population did not produce fruit in 2019 and we could not validate this result in a second year.

Based on models for established multiyear saplings, selection gradients of phenological traits were slightly different across years, in particular for budburst and female function. Considering 2018's growth phenology, female effective fecundity increased with earlier budburst, whereas male effective fecundity increased with early leaf senescence, as shown by the significant positive female directional selection gradients for BB-JD100 and SPSS, and male directional selection gradients for LLS and LS-JD287 (Table 3a; also partially reflected in Figure 3a based on PCs and total effective fecundity). For 2019's growth phenology, however, not only mothers with late budburst, as in 2018, but also those with very early budburst were selected against, resulting in overall negative directional and quadratic selection gradients for BB-JD100 that indicate stabilizing selection (Table 3b, Figure S5C,D). Multivariate selection gradients also supported stabilizing selection for female budburst (Figure 3c; Table 5). For both 2018's and 2019's autumn phenology, early leaf senescence was associated with higher male fecundity (Figure 3a,d; Table 5). Lastly, both single-trait (for traits measured in 2018 and 2019) and multivariate selection gradients based on total effective
<table><caption>TABLE 3 Trait variability and univariate SD-standardized selection gradients (linear models) based on standardized trait values measured in 2018 (a) and 2019 (b), and fitness estimates based on established multiyear saplings.</caption><thead><tr><th rowspan="2">Trait</th><th rowspan="2">N</th><th colspan="2">Trait values</th><th colspan="2">Female selection gradients (β<sub>fσ</sub>)</th><th colspan="2">Male selection gradients (β<sub>mσ</sub>)</th></tr><tr><th>Mean (μ<sub>p</sub>)</th><th>SD (σ<sub>p</sub>)</th><th>Linear</th><th>Quadratic</th><th>Linear</th><th>Quadratic</th></tr></thead><tbody><tr><td colspan="8"><strong>(a) Traits measured in year 2018</strong></td></tr><tr><td colspan="8"><strong>Growth</strong></td></tr><tr><td>DBH</td><td>251</td><td>34.761</td><td>10.667</td><td><strong>0.074**</strong></td><td>-0.012<sup>ns</sup></td><td><strong>0.070***</strong></td><td>0.022<sup>ns</sup></td></tr><tr><td>GL</td><td>251</td><td>188.414</td><td>7.595</td><td>0.013<sup>ns</sup></td><td>0.022<sup>ns</sup></td><td><strong>-0.038*</strong></td><td>0.024<sup>ns</sup></td></tr><tr><td colspan="8"><strong>Budburst</strong></td></tr><tr><td>LBB</td><td>251</td><td>106.141</td><td>37.888</td><td>-0.014<sup>ns</sup></td><td>-0.304<sup>ns</sup></td><td>0.029<sup>ns</sup></td><td>0.016<sup>ns</sup></td></tr><tr><td>BB-JD100</td><td>251</td><td>2.765</td><td>0.555</td><td><strong>0.093***</strong></td><td>-0.052<sup>ns</sup></td><td>0.020<sup>ns</sup></td><td>-0.048<sup>ns</sup></td></tr><tr><td>SPSS</td><td>251</td><td>16.263</td><td>0.821</td><td><strong>0.060*</strong></td><td>-0.006<sup>ns</sup></td><td>-0.007<sup>ns</sup></td><td>-0.012<sup>ns</sup></td></tr><tr><td colspan="8"><strong>Reproduction</strong></td></tr><tr><td>FLW</td><td>251</td><td>1.578</td><td>0.902</td><td><strong>0.060*</strong></td><td>0.016<sup>ns</sup></td><td><strong>0.049**</strong></td><td>0.038<sup>ns</sup></td></tr><tr><td>FT-JD108</td><td>251</td><td>1.948</td><td>1.227</td><td>0.050<sup>ns</sup></td><td>-0.100<sup>ns</sup></td><td><strong>0.042**</strong></td><td>-0.026<sup>ns</sup></td></tr><tr><td>FT-JD217</td><td>251</td><td>1.267</td><td>0.970</td><td><strong>0.069*</strong></td><td><strong>-0.142*</strong></td><td><strong>0.037*</strong></td><td><strong>-0.092*</strong></td></tr><tr><td colspan="8"><strong>Leaf senescence</strong></td></tr><tr><td>LLS</td><td>239</td><td>24.494</td><td>10.855</td><td>0.036<sup>ns</sup></td><td>-0.048<sup>ns</sup></td><td><strong>0.061***</strong></td><td>0.010<sup>ns</sup></td></tr><tr><td>LS-JD287</td><td>251</td><td>2.263</td><td>0.761</td><td>0.004<sup>ns</sup></td><td>-0.064<sup>ns</sup></td><td><strong>0.042**</strong></td><td>0.012<sup>ns</sup></td></tr><tr><td colspan="8"><strong>(b) Traits measured in year 2019</strong></td></tr><tr><td colspan="8"><strong>Growth</strong></td></tr><tr><td>GL</td><td>250</td><td>174.996</td><td>7.008</td><td><strong>-0.086**</strong></td><td>-0.008<sup>ns</sup></td><td><strong>-0.052**</strong></td><td>-0.006<sup>ns</sup></td></tr><tr><td colspan="8"><strong>Budburst</strong></td></tr><tr><td>LBB</td><td>251</td><td>54.853</td><td>18.965</td><td><strong>0.081**</strong></td><td>-0.006<sup>ns</sup></td><td>0.019<sup>ns</sup></td><td>-0.000<sup>ns</sup></td></tr><tr><td>BB-JD100</td><td>251</td><td>3.135</td><td>1.015</td><td><strong>-0.084**</strong></td><td><strong>-0.220**</strong></td><td>-0.025<sup>ns</sup></td><td>-0.070<sup>ns</sup></td></tr><tr><td>SPSS</td><td>251</td><td>16.837</td><td>2.231</td><td>-0.009<sup>ns</sup></td><td>-0.038<sup>ns</sup></td><td>-0.023<sup>ns</sup></td><td>-0.022<sup>ns</sup></td></tr><tr><td colspan="8"><strong>Leaf senescence</strong></td></tr><tr><td>LLS</td><td>247</td><td>28.109</td><td>8.956</td><td>-0.036<sup>ns</sup></td><td>-0.070<sup>ns</sup></td><td>0.017<sup>ns</sup></td><td><strong>-0.054*</strong></td></tr><tr><td>LS-JD267</td><td>251</td><td>2.171</td><td>0.709</td><td>0.043<sup>ns</sup></td><td>0.050<sup>ns</sup></td><td><strong>0.040*</strong></td><td>-0.004<sup>ns</sup></td></tr></tbody></table>

Note: Fitness estimates based on fecundity were log-transformed following Oddou-Muratorio et al. (2018). Models with non-significant quadratic component were rerun considering only linear selection gradients. GDD5 stands for Growing Degree Days with base temperature of 5°C. Significant univariate selection gradients are in bold. *α < 0.05; **α < 0.01; ***α < 0.001; ns, not significant.

Abbreviations: BB-JD100, Bud phenological stage at Julian day 100; DBH, Diameter at 1.30 m; FLW, Male flowering stage at Julian day 114; FT-JD108/JD217, Fruit abundance at Julian day 108/217; GL, Length of the growing season; LBB, Length of budburst; LLS, Length of leaf senescence; LS-JD287/JD267, Leaf senescence stage at Julian day 287/267; SPSS, Spring phenology score sum.

fecundity (only traits measured in 2019) showed negative directional selection for the length of the growing season, GL (Tables 3 and 5; Figure 3b).

## 3.4 | Heritability, evolvability and expected selection response

A prerequisite for a selection response (R) is that fitness is heritable. This condition was met since total effective fecundity had a narrow-sense heritability ($h^2$) of 0.158 (95% CIs: 0.049, 0.460; using DBH as covariate, see Material and Methods). Overall, estimates of narrow-sense heritability were low to moderate (albeit having a high

probability of being non-zero) across traits and years (mean 0.165), with only three phenotypic traits having $h^2 > 0.15$, all of which were related to growth (size) and spring growth phenology (GL, LBB and SPSS; Table 4). Evolvability was consistent and relatively low across years, with higher values for tree size ($l_A$ for DBH of ~0.01) and some phenological traits (LBB ~0.02 and LSS ~0.02 for 2018's measurements, and LBB ~0.04 and LSS ~0.01 for 2019's measurements). The highest expected responses to selection R, which were also consistent across years, were observed for growth traits (DBH, GL). Accordingly, natural selection in the studied stand is expected to increase tree size (DBH) and reduce GL, but the magnitude of the change in phenotypic values would be low (Table 4). With respect to budburst, the highest R was found for LBB, while the other traits are
<table><caption>TABLE 3 Trait variability and univariate SD-standardized selection gradients (linear models) based on standardized trait values measured in 2018 (a) and 2019 (b), and fitness estimates based on established multiyear saplings.</caption><thead><tr><th rowspan="2">Trait</th><th rowspan="2">N</th><th colspan="2">Trait values</th><th colspan="2">Female selection gradients (β<sub>fσ</sub>)</th><th colspan="2">Male selection gradients (β<sub>mσ</sub>)</th></tr><tr><th>Mean (μ<sub>p</sub>)</th><th>SD (σ<sub>p</sub>)</th><th>Linear</th><th>Quadratic</th><th>Linear</th><th>Quadratic</th></tr></thead><tbody><tr><td colspan="8"><strong>(a) Traits measured in year 2018</strong></td></tr><tr><td colspan="8"><strong>Growth</strong></td></tr><tr><td>DBH</td><td>251</td><td>34.761</td><td>10.667</td><td><strong>0.074**</strong></td><td>-0.012<sup>ns</sup></td><td><strong>0.070***</strong></td><td>0.022<sup>ns</sup></td></tr><tr><td>GL</td><td>251</td><td>188.414</td><td>7.595</td><td>0.013<sup>ns</sup></td><td>0.022<sup>ns</sup></td><td><strong>-0.038*</strong></td><td>0.024<sup>ns</sup></td></tr><tr><td colspan="8"><strong>Budburst</strong></td></tr><tr><td>LBB</td><td>251</td><td>106.141</td><td>37.888</td><td>-0.014<sup>ns</sup></td><td>-0.304<sup>ns</sup></td><td>0.029<sup>ns</sup></td><td>0.016<sup>ns</sup></td></tr><tr><td>BB-JD100</td><td>251</td><td>2.765</td><td>0.555</td><td><strong>0.093***</strong></td><td>-0.052<sup>ns</sup></td><td>0.020<sup>ns</sup></td><td>-0.048<sup>ns</sup></td></tr><tr><td>SPSS</td><td>251</td><td>16.263</td><td>0.821</td><td><strong>0.060*</strong></td><td>-0.006<sup>ns</sup></td><td>-0.007<sup>ns</sup></td><td>-0.012<sup>ns</sup></td></tr><tr><td colspan="8"><strong>Reproduction</strong></td></tr><tr><td>FLW</td><td>251</td><td>1.578</td><td>0.902</td><td><strong>0.060*</strong></td><td>0.016<sup>ns</sup></td><td><strong>0.049**</strong></td><td>0.038<sup>ns</sup></td></tr><tr><td>FT-JD108</td><td>251</td><td>1.948</td><td>1.227</td><td>0.050<sup>ns</sup></td><td>-0.100<sup>ns</sup></td><td><strong>0.042**</strong></td><td>-0.026<sup>ns</sup></td></tr><tr><td>FT-JD217</td><td>251</td><td>1.267</td><td>0.970</td><td><strong>0.069*</strong></td><td><strong>-0.142*</strong></td><td><strong>0.037*</strong></td><td><strong>-0.092*</strong></td></tr><tr><td colspan="8"><strong>Leaf senescence</strong></td></tr><tr><td>LLS</td><td>239</td><td>24.494</td><td>10.855</td><td>0.036<sup>ns</sup></td><td>-0.048<sup>ns</sup></td><td><strong>0.061***</strong></td><td>0.010<sup>ns</sup></td></tr><tr><td>LS-JD287</td><td>251</td><td>2.263</td><td>0.761</td><td>0.004<sup>ns</sup></td><td>-0.064<sup>ns</sup></td><td><strong>0.042**</strong></td><td>0.012<sup>ns</sup></td></tr><tr><td colspan="8"><strong>(b) Traits measured in year 2019</strong></td></tr><tr><td colspan="8"><strong>Growth</strong></td></tr><tr><td>GL</td><td>250</td><td>174.996</td><td>7.008</td><td><strong>-0.086**</strong></td><td>-0.008<sup>ns</sup></td><td><strong>-0.052**</strong></td><td>-0.006<sup>ns</sup></td></tr><tr><td colspan="8"><strong>Budburst</strong></td></tr><tr><td>LBB</td><td>251</td><td>54.853</td><td>18.965</td><td><strong>0.081**</strong></td><td>-0.006<sup>ns</sup></td><td>0.019<sup>ns</sup></td><td>-0.000<sup>ns</sup></td></tr><tr><td>BB-JD100</td><td>251</td><td>3.135</td><td>1.015</td><td><strong>-0.084**</strong></td><td><strong>-0.220**</strong></td><td>-0.025<sup>ns</sup></td><td>-0.070<sup>ns</sup></td></tr><tr><td>SPSS</td><td>251</td><td>16.837</td><td>2.231</td><td>-0.009<sup>ns</sup></td><td>-0.038<sup>ns</sup></td><td>-0.023<sup>ns</sup></td><td>-0.022<sup>ns</sup></td></tr><tr><td colspan="8"><strong>Leaf senescence</strong></td></tr><tr><td>LLS</td><td>247</td><td>28.109</td><td>8.956</td><td>-0.036<sup>ns</sup></td><td>-0.070<sup>ns</sup></td><td>0.017<sup>ns</sup></td><td><strong>-0.054*</strong></td></tr><tr><td>LS-JD267</td><td>251</td><td>2.171</td><td>0.709</td><td>0.043<sup>ns</sup></td><td>0.050<sup>ns</sup></td><td><strong>0.040*</strong></td><td>-0.004<sup>ns</sup></td></tr></tbody></table>

Note: Fitness estimates based on fecundity were log-transformed following Oddou-Muratorio et al. (2018). Models with non-significant quadratic component were rerun considering only linear selection gradients. GDD5 stands for Growing Degree Days with base temperature of 5°C. Significant univariate selection gradients are in bold. *α < 0.05; **α < 0.01; ***α < 0.001; ns, not significant.

Abbreviations: BB-JD100, Bud phenological stage at Julian day 100; DBH, Diameter at 1.30 m; FLW, Male flowering stage at Julian day 114; FT-JD108/JD217, Fruit abundance at Julian day 108/217; GL, Length of the growing season; LBB, Length of budburst; LLS, Length of leaf senescence; LS-JD287/JD267, Leaf senescence stage at Julian day 287/267; SPSS, Spring phenology score sum.

fecundity (only traits measured in 2019) showed negative directional selection for the length of the growing season, GL (Tables 3 and 5; Figure 3b).

## 3.4 | Heritability, evolvability and expected selection response

A prerequisite for a selection response (R) is that fitness is heritable. This condition was met since total effective fecundity had a narrow-sense heritability ($h^2$) of 0.158 (95% CIs: 0.049, 0.460; using DBH as covariate, see Material and Methods). Overall, estimates of narrow-sense heritability were low to moderate (albeit having a high

probability of being non-zero) across traits and years (mean 0.165), with only three phenotypic traits having $h^2 > 0.15$, all of which were related to growth (size) and spring growth phenology (GL, LBB and SPSS; Table 4). Evolvability was consistent and relatively low across years, with higher values for tree size ($l_A$ for DBH of ~0.01) and some phenological traits (LBB ~0.02 and LSS ~0.02 for 2018's measurements, and LBB ~0.04 and LSS ~0.01 for 2019's measurements). The highest expected responses to selection R, which were also consistent across years, were observed for growth traits (DBH, GL). Accordingly, natural selection in the studied stand is expected to increase tree size (DBH) and reduce GL, but the magnitude of the change in phenotypic values would be low (Table 4). With respect to budburst, the highest R was found for LBB, while the other traits are
<figure><img src="image_6.png" /></figure>

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

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

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

FIGURE 3 Multivariate SD-standardized selection gradients for total (a and b), female (c) and male (d) relative fitness estimated as effective fecundity using the model for established multiyear saplings and phenotypic traits measured in 2018 and 2019 after a Principal Component Analysis with VARIMAX rotation (3 PCs with eigenvalues >1, 72% and 74% of variance explained in 2018 and 2019 respectively; see main text and Supplemental Information); log-transformed relative fecundity >0 (i.e. relative fecundity >1) indicates higher than average contribution to reproductive output. ns, not significant.

not expected to respond to selection despite relatively high heritability for the SPSS synthetic index. This is mainly due to inconsistent selection gradients across male and female functions, and years, sometimes even showing opposite trends (e.g. for SPSS measured in 2018 and 2019). Finally, reproduction traits are expected to show a very low response to selection, due to both low phenotypic variance for fructification in the studied stand and low heritability (average $h^2$ of 0.099 for reproduction traits).

4 | DISCUSSION

In the present study, we have investigated the short-term evolutionary potential of a beech population in the core area of the species' distribution, that is, under (current) optimal growing conditions. The particular interest of this work lies in estimating key parameters of the expected selection responses (heritability, evolvability, selection gradients) for phenology, growth (size) and reproduction traits directly in the field. The vast majority of studies examining the evolutionary potential of forest tree populations are based on common garden experiments (but see Alexandre, Truffaut, Ducousso, et al., 2020; Alexandre, Truffaut, Klein, et al., 2020; Andrew et al., 2005; Bontemps et al., 2016; Castellanos et al., 2015) while in situ studies can provide a more realistic picture of the selection pressures that trees currently face and the genetic variation that is available for adaptation (Kruuk et al., 2008).

4.1 | Selection gradients for growth, phenology and reproduction traits

Growth, phenology, and reproduction traits are all potential targets of direct or indirect selection in plants, as suggested by their variation along environmental gradients (see Caignard et al., 2017; Vitasse et al., 2009, for forest trees). However, in their review of 871 plant traits, Geber and Griffen (2003) showed that direct selection accounts for only a small proportion of total selection and reported an absolute median value of 0.118 for 653 linear selection gradients and 0.029 for 214 quadratic selection gradients. In our study, the absolute median value was 0.061 for the 28 significant linear selection gradients and 0.142 for the seven significant quadratic selection gradients (all models confounded). Thus, our results seem to be typical for plant species. Evaluating multiple traits simultaneously is necessary as selection acts on integrated phenotypes (Schlichting, 1989), where (physiological) processes are balanced across many levels of organization (Bloom et al., 1985), leading to positive or negative correlations between traits (Arntz & Delph, 2001). Two traits can be correlated because of resource limitation (inducing trade-offs), shared genetic (i.e. pleiotropy) or developmental basis, shared functional role or a common plastic response to microenvironmental variation (Arntz & Delph, 2001; Denney et al., 2020; Schlichting, 1989).

In our study, larger trees bore more fruit and had higher fitness, as shown by the phenotypic correlation between tree size and fruit abundance, and by the significant directional selection gradients for
<figure><img src="image_5.png" /></figure>

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

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

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

FIGURE 3 Multivariate SD-standardized selection gradients for total (a and b), female (c) and male (d) relative fitness estimated as effective fecundity using the model for established multiyear saplings and phenotypic traits measured in 2018 and 2019 after a Principal Component Analysis with VARIMAX rotation (3 PCs with eigenvalues >1, 72% and 74% of variance explained in 2018 and 2019 respectively; see main text and Supplemental Information); log-transformed relative fecundity >0 (i.e. relative fecundity >1) indicates higher than average contribution to reproductive output. ns, not significant.

not expected to respond to selection despite relatively high heritability for the SPSS synthetic index. This is mainly due to inconsistent selection gradients across male and female functions, and years, sometimes even showing opposite trends (e.g. for SPSS measured in 2018 and 2019). Finally, reproduction traits are expected to show a very low response to selection, due to both low phenotypic variance for fructification in the studied stand and low heritability (average $h^2$ of 0.099 for reproduction traits).

4 | DISCUSSION

In the present study, we have investigated the short-term evolutionary potential of a beech population in the core area of the species' distribution, that is, under (current) optimal growing conditions. The particular interest of this work lies in estimating key parameters of the expected selection responses (heritability, evolvability, selection gradients) for phenology, growth (size) and reproduction traits directly in the field. The vast majority of studies examining the evolutionary potential of forest tree populations are based on common garden experiments (but see Alexandre, Truffaut, Ducousso, et al., 2020; Alexandre, Truffaut, Klein, et al., 2020; Andrew et al., 2005; Bontemps et al., 2016; Castellanos et al., 2015) while in situ studies can provide a more realistic picture of the selection pressures that trees currently face and the genetic variation that is available for adaptation (Kruuk et al., 2008).

4.1 | Selection gradients for growth, phenology and reproduction traits

Growth, phenology, and reproduction traits are all potential targets of direct or indirect selection in plants, as suggested by their variation along environmental gradients (see Caignard et al., 2017; Vitasse et al., 2009, for forest trees). However, in their review of 871 plant traits, Geber and Griffen (2003) showed that direct selection accounts for only a small proportion of total selection and reported an absolute median value of 0.118 for 653 linear selection gradients and 0.029 for 214 quadratic selection gradients. In our study, the absolute median value was 0.061 for the 28 significant linear selection gradients and 0.142 for the seven significant quadratic selection gradients (all models confounded). Thus, our results seem to be typical for plant species. Evaluating multiple traits simultaneously is necessary as selection acts on integrated phenotypes (Schlichting, 1989), where (physiological) processes are balanced across many levels of organization (Bloom et al., 1985), leading to positive or negative correlations between traits (Arntz & Delph, 2001). Two traits can be correlated because of resource limitation (inducing trade-offs), shared genetic (i.e. pleiotropy) or developmental basis, shared functional role or a common plastic response to microenvironmental variation (Arntz & Delph, 2001; Denney et al., 2020; Schlichting, 1989).

In our study, larger trees bore more fruit and had higher fitness, as shown by the phenotypic correlation between tree size and fruit abundance, and by the significant directional selection gradients for
<table><caption>TABLE 4 Narrow-sense heritability ($h^2$), evolvability ($I_n$) and expected response to selection ($R$) for phenotypic traits measured in 2018 (a) and 2019 (b).</caption><thead><tr><th colspan="7">(a) Traits measured in year 2018</th></tr><tr><th rowspan="2">Trait</th><th rowspan="2">Model (covariates)<sup>a</sup></th><th colspan="2">Variances (mode)</th><th rowspan="2">h<sup>2</sup> (95% CIs)</th><th rowspan="2">Evolvability (I<sub>n</sub>)</th><th rowspan="2">Selection response (R)</th></tr><tr><th>Additive</th><th>Residual</th></tr></thead><tbody><tr><td colspan="7"><strong>Growth</strong></td></tr><tr><td>DBH</td><td>Gaussian (NCI)</td><td>0.125</td><td>0.803</td><td><strong>0.103</strong> (0.042, 0.377)</td><td>0.0097</td><td>0.0945</td></tr><tr><td>GL</td><td>Gaussian (DBH)</td><td>0.531</td><td>0.338</td><td><strong>0.529</strong> (0.235, 0.904)</td><td>0.0009</td><td>-0.1012</td></tr><tr><td colspan="7"><strong>Budburst</strong></td></tr><tr><td>LBB</td><td>Gaussian (NCI)</td><td>0.177</td><td>0.829</td><td><strong>0.175</strong> (0.049, 0.445)</td><td>0.0223</td><td>0.1161</td></tr><tr><td>BB-JD100</td><td>Probit (DBH)</td><td>0.144</td><td>1</td><td><strong>0.069</strong> (0.022, 0.622)</td><td>NA</td><td>0.0021</td></tr><tr><td>SPSS</td><td>Gaussian (DBH+NCI)</td><td>0.184</td><td>0.771</td><td><strong>0.191</strong> (0.049, 0.588)</td><td>0.0005</td><td>0.0034</td></tr><tr><td colspan="7"><strong>Reproduction</strong></td></tr><tr><td>FLW</td><td>Probit (NCI)</td><td>0.185</td><td>1</td><td><strong>0.113</strong> (0.023, 0.481)</td><td>NA</td><td>0.0067</td></tr><tr><td>FT-JD108</td><td>Probit</td><td>0.177</td><td>1</td><td><strong>0.112</strong> (0.023, 0.543)</td><td>NA</td><td>0.0080</td></tr><tr><td>FT-JD217</td><td>Probit (NCI)</td><td>0.165</td><td>1</td><td><strong>0.072</strong> (0.024, 0.767)</td><td>NA</td><td>0.0043</td></tr><tr><td colspan="7"><strong>Leaf senescence</strong></td></tr><tr><td>LLS</td><td>Gaussian (DBH)</td><td>0.110</td><td>0.807</td><td><strong>0.109</strong> (0.040, 0.340)</td><td>0.0214</td><td>0.0758</td></tr><tr><td>LS-JD287</td><td>Probit (NCI)</td><td>0.156</td><td>1</td><td><strong>0.139</strong> (0.024, 0.680)</td><td>NA</td><td>0.0031</td></tr><tr><th colspan="7">(b) Traits measured in year 2019</th></tr><tr><th rowspan="2">Trait</th><th rowspan="2">Model (covariates)<sup>b</sup></th><th colspan="2">Variances (mode)</th><th rowspan="2">h<sup>2</sup> (95% CIs)</th><th rowspan="2">Evolvability (I<sub>n</sub>)</th><th rowspan="2">Selection response (R)</th></tr><tr><th>Additive</th><th>Residual</th></tr><tr><td colspan="7"><strong>Growth</strong></td></tr><tr><td>GL</td><td>Gaussian (DBH+NCI)</td><td>0.126</td><td>0.727</td><td><strong>0.145</strong> (0.046, 0.512)</td><td>0.0002</td><td>-0.0723</td></tr><tr><td colspan="7"><strong>Budburst</strong></td></tr><tr><td>LBB</td><td>Gaussian (NCI)</td><td>0.340</td><td>0.672</td><td><strong>0.368</strong> (0.096, 0.717)</td><td>0.0440</td><td>0.3035</td></tr><tr><td>BB-JD100</td><td>Probit (DBH)</td><td>0.121</td><td>1</td><td><strong>0.058</strong> (0.025, 0.357)</td><td>NA</td><td>-0.0024</td></tr><tr><td>SPSS</td><td>Gaussian (DBH+NCI)</td><td>0.186</td><td>0.648</td><td><strong>0.238</strong> (0.057, 0.551)</td><td>0.0042</td><td>-0.0144</td></tr><tr><td colspan="7"><strong>Leaf senescence</strong></td></tr><tr><td>LLS</td><td>Gaussian (DBH)</td><td>0.141</td><td>0.784</td><td><strong>0.122</strong> (0.042, 0.464)</td><td>0.0124</td><td>-0.0013</td></tr><tr><td>LS-JD267</td><td>Probit (NCI)</td><td>0.233</td><td>1</td><td><strong>0.091</strong> (0.023, 0.672)</td><td>NA</td><td>0.0029</td></tr></tbody></table>

Note: Both naïve models (no covariates) and models including diameter (DBH), a proxy for age, and/or the Neighbour-based Competition Index (NCI) as fixed factors were tested. Only the best models (i.e. the models with lowest DIC and all fixed factors significant) are shown. Heritability estimates are considered different from zero when 95% CIs do not overlap 0.010 (notice that heritability values are higher than zero by construction). Significant heritability estimates are in bold. Expected selection responses were computed considering univariate directional selection gradients for total effective fecundity from the models based on established multiyear saplings (see text). NA: not applicable for traits based on an interval scale, for example, phenology stages (Hansen et al., 2011).

Abbreviations: BB-JD100, Bud phenological stage at Julian day 100; DBH, Diameter at 1.30 m; FLW, Male flowering stage at Julian day 114; FT-JD108/JD217, Fruit abundance at Julian day 108/217; GL, Length of the growing season; LBB, Length of budburst; LLS, Length of leaf senescence; LS-JD287/JD267, Leaf senescence stage at Julian day 287/267; SPSS, Spring phenology score sum.

DBH in all models. This was expected as large trees generally allocate more resources to reproduction (Klinkhamer et al., 1997; Oddou-Muratorio et al., 2018), allowing them to produce considerable amounts of seeds and thus dominate recruitment at local and landscape scales (Lindenmayer & Laurance, 2017; Thomas, 2011). This result is also consistent with a study of natural oak populations that found significant linear selection gradients for radial growth-related

traits (Alexandre, Truffaut, Ducoussou, et al., 2020). Larger trees also had an earlier male flowering phenology, resulting in a consistent signal of directional selection for this trait. Selection for early flowering is common in short-lived herbaceous plants (Munguía-Rosas et al., 2011) and in the case of beech, it can probably also be explained by protogyny, as early male flowers are more likely to be in phase with female flowers.

<footnote><sup>a</sup>Models FLW, FT-JD217, LS-JD287~DBH/DBH+NCI, FT-JD108~DBH/NCI/DBH+NCI did not converge.</footnote>

<footnote><sup>b</sup>Model LS-JD267~DBH+NCI did not converge.</footnote>
<table><caption>TABLE 4 Narrow-sense heritability ($h^2$), evolvability ($I_A$) and expected response to selection ($R$) for phenotypic traits measured in 2018 (a) and 2019 (b).</caption><thead><tr><th colspan="7">(a) Traits measured in year 2018</th></tr><tr><th rowspan="2">Trait</th><th rowspan="2">Model (covariates)<sup>a</sup></th><th colspan="2">Variances (mode)</th><th rowspan="2">h<sup>2</sup> (95% CIs)</th><th rowspan="2">Evolvability (I<sub>A</sub>)</th><th rowspan="2">Selection response (R)</th></tr><tr><th>Additive</th><th>Residual</th></tr></thead><tbody><tr><td colspan="7"><strong>Growth</strong></td></tr><tr><td>DBH</td><td>Gaussian (NCI)</td><td>0.125</td><td>0.803</td><td><strong>0.103</strong> (0.042, 0.377)</td><td>0.0097</td><td>0.0945</td></tr><tr><td>GL</td><td>Gaussian (DBH)</td><td>0.531</td><td>0.338</td><td><strong>0.529</strong> (0.235, 0.904)</td><td>0.0009</td><td>-0.1012</td></tr><tr><td colspan="7"><strong>Budburst</strong></td></tr><tr><td>LBB</td><td>Gaussian (NCI)</td><td>0.177</td><td>0.829</td><td><strong>0.175</strong> (0.049, 0.445)</td><td>0.0223</td><td>0.1161</td></tr><tr><td>BB-JD100</td><td>Probit (DBH)</td><td>0.144</td><td>1</td><td><strong>0.069</strong> (0.022, 0.622)</td><td>NA</td><td>0.0021</td></tr><tr><td>SPSS</td><td>Gaussian (DBH+NCI)</td><td>0.184</td><td>0.771</td><td><strong>0.191</strong> (0.049, 0.588)</td><td>0.0005</td><td>0.0034</td></tr><tr><td colspan="7"><strong>Reproduction</strong></td></tr><tr><td>FLW</td><td>Probit (NCI)</td><td>0.185</td><td>1</td><td><strong>0.113</strong> (0.023, 0.481)</td><td>NA</td><td>0.0067</td></tr><tr><td>FT-JD108</td><td>Probit</td><td>0.177</td><td>1</td><td><strong>0.112</strong> (0.023, 0.543)</td><td>NA</td><td>0.0080</td></tr><tr><td>FT-JD217</td><td>Probit (NCI)</td><td>0.165</td><td>1</td><td><strong>0.072</strong> (0.024, 0.767)</td><td>NA</td><td>0.0043</td></tr><tr><td colspan="7"><strong>Leaf senescence</strong></td></tr><tr><td>LLS</td><td>Gaussian (DBH)</td><td>0.110</td><td>0.807</td><td><strong>0.109</strong> (0.040, 0.340)</td><td>0.0214</td><td>0.0758</td></tr><tr><td>LS-JD287</td><td>Probit (NCI)</td><td>0.156</td><td>1</td><td><strong>0.139</strong> (0.024, 0.680)</td><td>NA</td><td>0.0031</td></tr><tr><th colspan="7">(b) Traits measured in year 2019</th></tr><tr><th rowspan="2">Trait</th><th rowspan="2">Model (covariates)<sup>b</sup></th><th colspan="2">Variances (mode)</th><th rowspan="2">h<sup>2</sup> (95% CIs)</th><th rowspan="2">Evolvability (I<sub>A</sub>)</th><th rowspan="2">Selection response (R)</th></tr><tr><th>Additive</th><th>Residual</th></tr><tr><td colspan="7"><strong>Growth</strong></td></tr><tr><td>GL</td><td>Gaussian (DBH+NCI)</td><td>0.126</td><td>0.727</td><td><strong>0.145</strong> (0.046, 0.512)</td><td>0.0002</td><td>-0.0723</td></tr><tr><td colspan="7"><strong>Budburst</strong></td></tr><tr><td>LBB</td><td>Gaussian (NCI)</td><td>0.340</td><td>0.672</td><td><strong>0.368</strong> (0.096, 0.717)</td><td>0.0440</td><td>0.3035</td></tr><tr><td>BB-JD100</td><td>Probit (DBH)</td><td>0.121</td><td>1</td><td><strong>0.058</strong> (0.025, 0.357)</td><td>NA</td><td>-0.0024</td></tr><tr><td>SPSS</td><td>Gaussian (DBH+NCI)</td><td>0.186</td><td>0.648</td><td><strong>0.238</strong> (0.057, 0.551)</td><td>0.0042</td><td>-0.0144</td></tr><tr><td colspan="7"><strong>Leaf senescence</strong></td></tr><tr><td>LLS</td><td>Gaussian (DBH)</td><td>0.141</td><td>0.784</td><td><strong>0.122</strong> (0.042, 0.464)</td><td>0.0124</td><td>-0.0013</td></tr><tr><td>LS-JD267</td><td>Probit (NCI)</td><td>0.233</td><td>1</td><td><strong>0.091</strong> (0.023, 0.672)</td><td>NA</td><td>0.0029</td></tr></tbody></table>

Note: Both naïve models (no covariates) and models including diameter (DBH), a proxy for age, and/or the Neighbour-based Competition Index (NCI) as fixed factors were tested. Only the best models (i.e. the models with lowest DIC and all fixed factors significant) are shown. Heritability estimates are considered different from zero when 95% CIs do not overlap 0.010 (notice that heritability values are higher than zero by construction). Significant heritability estimates are in bold. Expected selection responses were computed considering univariate directional selection gradients for total effective fecundity from the models based on established multiyear saplings (see text). NA: not applicable for traits based on an interval scale, for example, phenology stages (Hansen et al., 2011).

Abbreviations: BB-JD100, Bud phenological stage at Julian day 100; DBH, Diameter at 1.30 m; FLW, Male flowering stage at Julian day 114; FT-JD108/JD217, Fruit abundance at Julian day 108/217; GL, Length of the growing season; LBB, Length of budburst; LLS, Length of leaf senescence; LS-JD287/JD267, Leaf senescence stage at Julian day 287/267; SPSS, Spring phenology score sum.

DBH in all models. This was expected as large trees generally allocate more resources to reproduction (Klinkhamer et al., 1997; Oddou-Muratorio et al., 2018), allowing them to produce considerable amounts of seeds and thus dominate recruitment at local and landscape scales (Lindenmayer & Laurance, 2017; Thomas, 2011). This result is also consistent with a study of natural oak populations that found significant linear selection gradients for radial growth-related

traits (Alexandre, Truffaut, Ducousso, et al., 2020). Larger trees also had an earlier male flowering phenology, resulting in a consistent signal of directional selection for this trait. Selection for early flowering is common in short-lived herbaceous plants (Munguía-Rosas et al., 2011) and in the case of beech, it can probably also be explained by protogyny, as early male flowers are more likely to be in phase with female flowers.

<footnote><sup>a</sup>Models FLW, FT-JD217, LS-JD287~DBH/DBH+NCI, FT-JD108~DBH/NCI/DBH+NCI did not converge.</footnote>

<footnote><sup>b</sup>Model LS-JD267~DBH+NCI did not converge.</footnote>
<table><caption>TABLE 5 Multivariate SD-standardized selection gradients for total, female and male relative fitness estimated as effective fecundity using the model for established multiyear saplings and phenotypic traits measured in 2018 (a) and 2019 (b) after a Principal Component Analysis with *VARIMAX* rotation (3 PCs with eigenvalues >1, 72% and 74% of variance explained for traits measured in 2018 and 2019 respectively).</caption><thead><tr><th rowspan="2">Trait</th><th colspan="2">Female selection gradients (β<sub>fα</sub>)</th><th colspan="2">Male selection gradients (β<sub>mα</sub>)</th><th colspan="2">Total fecundity selection gradients (β<sub>tα</sub>)</th></tr><tr><th>Linear</th><th>Quadratic</th><th>Linear</th><th>Quadratic</th><th>Linear</th><th>Quadratic</th></tr></thead><tbody><tr><td colspan="7">(a) Traits measured in year 2018</td></tr><tr><td>PC1</td><td><strong>0.074**</strong></td><td>NA</td><td><strong>0.046**</strong></td><td>NA</td><td><strong>0.070***</strong></td><td>NA</td></tr><tr><td>PC2</td><td><strong>0.051*</strong></td><td>NA</td><td>-0.022<sup>ns</sup></td><td>NA</td><td>0.004<sup>ns</sup></td><td>NA</td></tr><tr><td>PC3</td><td>0.014<sup>ns</sup></td><td>NA</td><td><strong>0.052**</strong></td><td>NA</td><td><strong>0.044*</strong></td><td>NA</td></tr><tr><td colspan="7">(b) Traits measured in year 2019</td></tr><tr><td>PC1</td><td>-0.020<sup>ns</sup></td><td><strong>-0.092*</strong></td><td>-0.022<sup>ns</sup></td><td>-0.008<sup>ns</sup></td><td>-0.020<sup>ns</sup></td><td>NA</td></tr><tr><td>PC2</td><td><strong>0.056*</strong></td><td>-0.004<sup>ns</sup></td><td><strong>0.076***</strong></td><td><strong>0.052*</strong></td><td><strong>0.069***</strong></td><td>NA</td></tr><tr><td>PC3</td><td>-0.095**</td><td>-0.018<sup>ns</sup></td><td>-0.022<sup>ns</sup></td><td>-0.040<sup>ns</sup></td><td><strong>-0.052**</strong></td><td>NA</td></tr></tbody></table>

Note: Models with quadratic components are shown only when one of the quadratic components was significant. Significant univariate selection gradients are in bold. α < 0.1; *α < 0.05; **α < 0.01; ***α < 0.001; ns, not significant.

Counterintuitively, large trees also had an earlier leaf senescence and a shorter growing season (GL), a pattern also found in Gömöry et al. (2007), which in turn was associated with higher male and female (only for GL measured in 2019) effective fecundity in the studied population. Early leaf senescence could be an advantage for trees that efficiently transfer the nutrients needed for growth and fruitification to the next year and do not waste energy protecting leaves during a dry and cold fall/winter (Netzer et al., 2018; Staaf & Stjernquist, 1986). However, late leaf senescence would also increase photosynthesis in the fall (Toomey et al., 2015), increasing the amount of resources available for leaf and flower formation in the next growing season (Hoch et al., 2013). Noticeably, xylem increment of beech in Slovenia generally extends from April to early August, with 75% of growth being achieved by the beginning of July, while cells require approximately a month longer to mature (Prislan et al., 2018). When most growth is finalized early, the length of the growing season may be less critical than, for example, spring phenology traits controlling the start of the growing season (see below). Early senescence, and other growth phenology traits, may also be indirectly selected because of pleiotropy involving flowering date, as shown in poplar (Böhlenius et al., 2006) and Norway spruce (Gyllenstrand et al., 2007).

## 4.2 | Fluctuating selection across ontogenetic stages and temporal environmental conditions

Only tree size and male flowering phenology, a trait correlated to tree size (Figure S3A), showed consistent selection gradients across seed, one-year-old seedling and established-multiyear sapling (5–9 years old) models. A possible explanation for inconsistent selection gradients across cohorts is the existence of genotype-by-age interactions, as natural selection can act in opposite directions at different ontogenetic stages (see reviews in Donohue et al., 2010; Schulter et al., 1991). For example, earlier leaf emergence of trees in the understory is explained by their earlier ontogenetic stage

(Vitasse, 2013). However, interestingly, no significant selection gradients for one-year-old seedlings or seeds were observed for several traits with significant selection gradients operating on the established multiyear sapling cohort (e.g. growing season length, SPSS and leaf senescence) highlighting the importance of the strong selection phase that marks the transition from seedlings to saplings. The role of early natural selection modifying initial seed rain patterns has been highlighted in other forest trees (e.g. González-Martínez et al., 2006; Warwell & Shaw, 2019). In addition, Hufford and Hamrick (2003) showed that changes in offspring genetic composition were greatest between the mature seeds and the established saplings that will form the core of the next generation. Strong selection on seedlings could be caused by a variety of factors, both biotic (e.g. the absence of mycorrhiza or predation) and abiotic (e.g. soil infertility, lack of shelter by shadow or low humidity).

Another possible explanation for the varying selection gradients found in our study is the environmental variation associated with different regeneration events. Selection pressures may vary from year to year depending on year-specific environmental conditions (especially climatic conditions). This has already been reported for the growth rhythm of whitebark pine (Warwell & Shaw, 2018), and various reviews (albeit mostly based on vertebrates) suggest that changes in selection direction across years are common (Siepielski et al., 2009; but see Morrissey & Hadfield, 2012). Moreover, genotype-by-environment interaction (G×E) is substantial in forest trees (see Table 1 in Li et al., 2017) and thus plastic responses of parent trees may also underlie differences in selection gradients across years. Interestingly, Lee et al. (2016) found that variation in climate across years generated more G×E than soil water availability across experimental sites in maize. In our study, for example, budburst showed an inconsistent signal of selection across years, with female effective fecundity being higher for trees with faster and earlier budburst based on traits measured in 2018 while it was higher for trees with longer duration of budburst, stabilizing around a relatively early budburst date (Julian day 100), based on traits measured in 2019.
<table><caption>TABLE 5 Multivariate SD-standardized selection gradients for total, female and male relative fitness estimated as effective fecundity using the model for established multiyear saplings and phenotypic traits measured in 2018 (a) and 2019 (b) after a Principal Component Analysis with *VARIMAX* rotation (3 PCs with eigenvalues >1, 72% and 74% of variance explained for traits measured in 2018 and 2019 respectively).</caption><thead><tr><th rowspan="2">Trait</th><th colspan="2">Female selection gradients (β<sub>fα</sub>)</th><th colspan="2">Male selection gradients (β<sub>mα</sub>)</th><th colspan="2">Total fecundity selection gradients (β<sub>tα</sub>)</th></tr><tr><th>Linear</th><th>Quadratic</th><th>Linear</th><th>Quadratic</th><th>Linear</th><th>Quadratic</th></tr></thead><tbody><tr><td colspan="7">(a) Traits measured in year 2018</td></tr><tr><td>PC1</td><td><strong>0.074**</strong></td><td>NA</td><td><strong>0.046**</strong></td><td>NA</td><td><strong>0.070***</strong></td><td>NA</td></tr><tr><td>PC2</td><td><strong>0.051*</strong></td><td>NA</td><td>-0.022<sup>ns</sup></td><td>NA</td><td>0.004<sup>ns</sup></td><td>NA</td></tr><tr><td>PC3</td><td>0.014<sup>ns</sup></td><td>NA</td><td><strong>0.052**</strong></td><td>NA</td><td><strong>0.044*</strong></td><td>NA</td></tr><tr><td colspan="7">(b) Traits measured in year 2019</td></tr><tr><td>PC1</td><td>-0.020<sup>ns</sup></td><td><strong>-0.092*</strong></td><td>-0.022<sup>ns</sup></td><td>-0.008<sup>ns</sup></td><td>-0.020<sup>ns</sup></td><td>NA</td></tr><tr><td>PC2</td><td><strong>0.056*</strong></td><td>-0.004<sup>ns</sup></td><td><strong>0.076***</strong></td><td><strong>0.052*</strong></td><td><strong>0.069***</strong></td><td>NA</td></tr><tr><td>PC3</td><td>-0.095**</td><td>-0.018<sup>ns</sup></td><td>-0.022<sup>ns</sup></td><td>-0.040<sup>ns</sup></td><td><strong>-0.052**</strong></td><td>NA</td></tr></tbody></table>

Note: Models with quadratic components are shown only when one of the quadratic components was significant. Significant univariate selection gradients are in bold. α < 0.1; *α < 0.05; **α < 0.01; ***α < 0.001; ns, not significant.

Counterintuitively, large trees also had an earlier leaf senescence and a shorter growing season (GL), a pattern also found in Gömöry et al. (2007), which in turn was associated with higher male and female (only for GL measured in 2019) effective fecundity in the studied population. Early leaf senescence could be an advantage for trees that efficiently transfer the nutrients needed for growth and fruitification to the next year and do not waste energy protecting leaves during a dry and cold fall/winter (Netzer et al., 2018; Staaf & Stjernquist, 1986). However, late leaf senescence would also increase photosynthesis in the fall (Toomey et al., 2015), increasing the amount of resources available for leaf and flower formation in the next growing season (Hoch et al., 2013). Noticeably, xylem increment of beech in Slovenia generally extends from April to early August, with 75% of growth being achieved by the beginning of July, while cells require approximately a month longer to mature (Prislan et al., 2018). When most growth is finalized early, the length of the growing season may be less critical than, for example, spring phenology traits controlling the start of the growing season (see below). Early senescence, and other growth phenology traits, may also be indirectly selected because of pleiotropy involving flowering date, as shown in poplar (Böhlenius et al., 2006) and Norway spruce (Gyllenstrand et al., 2007).

## 4.2 | Fluctuating selection across ontogenetic stages and temporal environmental conditions

Only tree size and male flowering phenology, a trait correlated to tree size (Figure S3A), showed consistent selection gradients across seed, one-year-old seedling and established-multiyear sapling (5–9 years old) models. A possible explanation for inconsistent selection gradients across cohorts is the existence of genotype-by-age interactions, as natural selection can act in opposite directions at different ontogenetic stages (see reviews in Donohue et al., 2010; Schulter et al., 1991). For example, earlier leaf emergence of trees in the understory is explained by their earlier ontogenetic stage

(Vitasse, 2013). However, interestingly, no significant selection gradients for one-year-old seedlings or seeds were observed for several traits with significant selection gradients operating on the established multiyear sapling cohort (e.g. growing season length, SPSS and leaf senescence) highlighting the importance of the strong selection phase that marks the transition from seedlings to saplings. The role of early natural selection modifying initial seed rain patterns has been highlighted in other forest trees (e.g. González-Martínez et al., 2006; Warwell & Shaw, 2019). In addition, Hufford and Hamrick (2003) showed that changes in offspring genetic composition were greatest between the mature seeds and the established saplings that will form the core of the next generation. Strong selection on seedlings could be caused by a variety of factors, both biotic (e.g. the absence of mycorrhiza or predation) and abiotic (e.g. soil infertility, lack of shelter by shadow or low humidity).

Another possible explanation for the varying selection gradients found in our study is the environmental variation associated with different regeneration events. Selection pressures may vary from year to year depending on year-specific environmental conditions (especially climatic conditions). This has already been reported for the growth rhythm of whitebark pine (Warwell & Shaw, 2018), and various reviews (albeit mostly based on vertebrates) suggest that changes in selection direction across years are common (Siepielski et al., 2009; but see Morrissey & Hadfield, 2012). Moreover, genotype-by-environment interaction (GxE) is substantial in forest trees (see Table 1 in Li et al., 2017) and thus plastic responses of parent trees may also underlie differences in selection gradients across years. Interestingly, Lee et al. (2016) found that variation in climate across years generated more GxE than soil water availability across experimental sites in maize. In our study, for example, budburst showed an inconsistent signal of selection across years, with female effective fecundity being higher for trees with faster and earlier budburst based on traits measured in 2018 while it was higher for trees with longer duration of budburst, stabilizing around a relatively early budburst date (Julian day 100), based on traits measured in 2019.
This may be the result of the notable differences in temperature sum accumulation in spring 2018 and 2019, suggesting substantial G×E (see Figure S3D). This interpretation is also consistent with the weak trait correlations found across years in our study, especially for bud-burst, and if true, would likely lead to increased genetic variation within established saplings, even if responses to selection are limited as predicted in this study.

Alternatively to directional selection for earlier phenology, the budburst patterns in 2018 can be explained based on allocation theory. Indeed, in accordance with the resource-switching hypothesis (see Pearse et al., 2016), faster and earlier budburst in 2018 could be related to the need to acquire more resources and allocate them to costly flowering and fruit production in mast years (see Hoch et al., 2013). Other studies in beech have shown that trees producing more fruit also invest less in leaf production (Innes, 1994) and growth (Hackett-Pain et al., 2018), and that vegetative growth and reproduction compete for carbon and nitrogen resources in mast years (Han et al., 2011). In the case of the piñon pine, it has been shown that this species prefers to allocate resources for reproduction rather than for defence against stress factors during mast years (Redmond et al., 2019). Reproduction is fundamental to the adaptive capacity of plants (Bonser, 2013), and therefore patterns of resource allocation for reproduction may alter selection gradients between mast and non-mast years. The patterns of stabilizing selection suggested by selection gradients for budburst measured in 2019 could reflect some adaptation to late frosts in non-mast years in the study area (where late frost is relatively frequent; see Figure S1).

Finally, the stable selection gradients across models for tree size and male flowering phenology can be explained either by consistent selection pressure across ontogenetic stages, by strong selection gradients for early regeneration events (seeds and one-year-old seedlings) that cannot be erased by subsequent processes, or by same-direction selection across years, independent of environmental conditions. Because our study covered only a short timescale in a long-lived tree species, even the consistent selection observed for these traits over the study period could be offset by strong and punctual opposing pressures due to stochastic extreme events, for example, intense droughts that have a greater impact on taller trees (Bennett et al., 2015) or late frosts that may affect early flowering males (Inouye, 2000).

## 4.3 | Low to moderate heritability and evolvability

Narrow-sense heritability, *h*², and evolvability, *I*<sub>A</sub>, are key components of the response to selection (Visscher et al., 2008). Assessing *h*² and *I*<sub>A</sub> of populations in situ has only recently become possible thanks to advances in genomic quantitative genetics (Kruuk et al., 2008), which in particular allows the estimation of relatedness using markers (Gienapp et al., 2017). For this reason, heritability in forest trees has been estimated almost exclusively in common garden experiments. In their review, Lind et al. (2018) found higher heritability for reproduction and phenology traits (*h*²~0.50) than for

growth traits (*h*²~0.25), a pattern also found in Caignard et al. (2019) and Scotti-Saintagne et al. (2004). This was not the case in our study, where the highest heritability was found for length of the growing season (GL) measured in 2018 (*h*²~0.53), while tree size (DBH) showed similar heritability to phenology and reproduction traits (*h*²~0.13). More specifically, the mean heritability across years was 0.26 for growth traits, 0.18 for budburst, 0.10 for reproduction traits and 0.12 for leaf senescence. By comparison, in forest tree common gardens, the timing of budburst generally has heritability above 0.80 (Alberto et al., 2011; Baliuckas & Pliura, 2003; Christophe & Birot, 1979; Rehfeldt, 1983), which is generally slightly higher than that of budset (Howe et al., 2003). Like other growth traits, DBH has low to moderate heritability, for example, a median value of 0.19 in a review of 67 studies (Cornelius, 1994), between 0.15 and 0.25 in pedunculate oak (Caignard et al., 2019), 0.21 in stone pine (Mutke et al., 2005) or between 0.09 and 0.43 in Scots pine (Bilir et al., 2006). Estimates of heritability for reproduction traits are rarer in forest trees, mainly because they require work on old, mature trees, and are highly age-dependent (Lind et al., 2018). For cone or seed production, *h*² was found to vary between 0.39 and 0.64 in lodgepole pine (Hannerz et al., 2001) and between 0.46 and 0.69 in pedunculate oak (Caignard et al., 2019).

Overall, our estimates were therefore lower than those previously reported based on common gardens, in particular for growth phenology (budburst and leaf senescence) and reproduction traits, and had wide credibility intervals. This is to be expected for heritability estimated directly in the field (Geber & Griffen, 2003), as there are fewer individuals per family to estimate the relatedness matrix accurately (Alexandre, Truffaut, Ducousso, et al., 2020) and larger microenvironmental variation than in controlled experiments (Hermida et al., 2013; Sgrò & Hoffmann, 1998). Moreover, if relatives tend to share the same environment, which is expected in natural populations of forest trees due to restricted seed dispersal (Asuka et al., 2005; Chybicki & Burczyk, 2010), this can inflate the estimated genetic variance (Rausher, 1992). Genetic variance may also be overestimated by maternal effects (Kruuk et al., 2008), which are difficult to account for in the field. To statistically control for such potential biases, we included a Neighbour-based Competition Index (NCI, a proxy for competition) and tree size (as estimated by DBH, a proxy for age and soil variation) as covariates in the models. Moreover, maternal effects in beech adult trees are expected to be negligible (see Gauzere et al., 2016) and the low to moderate heritability estimates for most traits in this study suggest that genetic variance was not overestimated.

Studies estimating heritability and evolvability in situ remain rare in forest trees (but see (Alexandre, Truffaut, Ducousso, et al., 2020; Andrew et al., 2005; Bontemps et al., 2016; Castellanos et al., 2015), and even in plants in general (but see Kulbaba et al., 2019; Sedlacek et al., 2016). Bontemps et al. (2016) and Alexandre, Truffaut, Ducousso, et al. (2020) evaluated comparable traits to our study, but still obtained much higher heritability estimates, even for the same tree species. For example, in beech, Bontemps et al. (2016) estimated narrow-sense heritability (*h*²) of 0.84–0.92 for SPSS (0.19–0.24 in
This may be the result of the notable differences in temperature sum accumulation in spring 2018 and 2019, suggesting substantial G×E (see Figure S3D). This interpretation is also consistent with the weak trait correlations found across years in our study, especially for bud-burst, and if true, would likely lead to increased genetic variation within established saplings, even if responses to selection are limited as predicted in this study.

Alternatively to directional selection for earlier phenology, the budburst patterns in 2018 can be explained based on allocation theory. Indeed, in accordance with the resource-switching hypothesis (see Pearse et al., 2016), faster and earlier budburst in 2018 could be related to the need to acquire more resources and allocate them to costly flowering and fruit production in mast years (see Hoch et al., 2013). Other studies in beech have shown that trees producing more fruit also invest less in leaf production (Innes, 1994) and growth (Hackett-Pain et al., 2018), and that vegetative growth and reproduction compete for carbon and nitrogen resources in mast years (Han et al., 2011). In the case of the piñon pine, it has been shown that this species prefers to allocate resources for reproduction rather than for defence against stress factors during mast years (Redmond et al., 2019). Reproduction is fundamental to the adaptive capacity of plants (Bonser, 2013), and therefore patterns of resource allocation for reproduction may alter selection gradients between mast and non-mast years. The patterns of stabilizing selection suggested by selection gradients for budburst measured in 2019 could reflect some adaptation to late frosts in non-mast years in the study area (where late frost is relatively frequent; see Figure S1).

Finally, the stable selection gradients across models for tree size and male flowering phenology can be explained either by consistent selection pressure across ontogenetic stages, by strong selection gradients for early regeneration events (seeds and one-year-old seedlings) that cannot be erased by subsequent processes, or by same-direction selection across years, independent of environmental conditions. Because our study covered only a short timescale in a long-lived tree species, even the consistent selection observed for these traits over the study period could be offset by strong and punctual opposing pressures due to stochastic extreme events, for example, intense droughts that have a greater impact on taller trees (Bennett et al., 2015) or late frosts that may affect early flowering males (Inouye, 2000).

## 4.3 | Low to moderate heritability and evolvability

Narrow-sense heritability, *h*², and evolvability, *I*<sub>A</sub>, are key components of the response to selection (Visscher et al., 2008). Assessing *h*² and *I*<sub>A</sub> of populations in situ has only recently become possible thanks to advances in genomic quantitative genetics (Kruuk et al., 2008), which in particular allows the estimation of relatedness using markers (Gienapp et al., 2017). For this reason, heritability in forest trees has been estimated almost exclusively in common garden experiments. In their review, Lind et al. (2018) found higher heritability for reproduction and phenology traits (*h*²~0.50) than for

growth traits (*h*²~0.25), a pattern also found in Caignard et al. (2019) and Scotti-Saintagne et al. (2004). This was not the case in our study, where the highest heritability was found for length of the growing season (GL) measured in 2018 (*h*²~0.53), while tree size (DBH) showed similar heritability to phenology and reproduction traits (*h*²~0.13). More specifically, the mean heritability across years was 0.26 for growth traits, 0.18 for budburst, 0.10 for reproduction traits and 0.12 for leaf senescence. By comparison, in forest tree common gardens, the timing of budburst generally has heritability above 0.80 (Alberto et al., 2011; Baliuckas & Pliura, 2003; Christophe & Birot, 1979; Rehfeldt, 1983), which is generally slightly higher than that of budset (Howe et al., 2003). Like other growth traits, DBH has low to moderate heritability, for example, a median value of 0.19 in a review of 67 studies (Cornelius, 1994), between 0.15 and 0.25 in pedunculate oak (Caignard et al., 2019), 0.21 in stone pine (Mutke et al., 2005) or between 0.09 and 0.43 in Scots pine (Bilir et al., 2006). Estimates of heritability for reproduction traits are rarer in forest trees, mainly because they require work on old, mature trees, and are highly age-dependent (Lind et al., 2018). For cone or seed production, *h*² was found to vary between 0.39 and 0.64 in lodgepole pine (Hannerz et al., 2001) and between 0.46 and 0.69 in pedunculate oak (Caignard et al., 2019).

Overall, our estimates were therefore lower than those previously reported based on common gardens, in particular for growth phenology (budburst and leaf senescence) and reproduction traits, and had wide credibility intervals. This is to be expected for heritability estimated directly in the field (Geber & Griffen, 2003), as there are fewer individuals per family to estimate the relatedness matrix accurately (Alexandre, Truffaut, Ducousso, et al., 2020) and larger microenvironmental variation than in controlled experiments (Hermida et al., 2013; Sgrò & Hoffmann, 1998). Moreover, if relatives tend to share the same environment, which is expected in natural populations of forest trees due to restricted seed dispersal (Asuka et al., 2005; Chybicki & Burczyk, 2010), this can inflate the estimated genetic variance (Rausher, 1992). Genetic variance may also be overestimated by maternal effects (Kruuk et al., 2008), which are difficult to account for in the field. To statistically control for such potential biases, we included a Neighbour-based Competition Index (NCI, a proxy for competition) and tree size (as estimated by DBH, a proxy for age and soil variation) as covariates in the models. Moreover, maternal effects in beech adult trees are expected to be negligible (see Gauzere et al., 2016) and the low to moderate heritability estimates for most traits in this study suggest that genetic variance was not overestimated.

Studies estimating heritability and evolvability in situ remain rare in forest trees (but see (Alexandre, Truffaut, Ducousso, et al., 2020; Andrew et al., 2005; Bontemps et al., 2016; Castellanos et al., 2015), and even in plants in general (but see Kulbaba et al., 2019; Sedlacek et al., 2016). Bontemps et al. (2016) and Alexandre, Truffaut, Ducousso, et al. (2020) evaluated comparable traits to our study, but still obtained much higher heritability estimates, even for the same tree species. For example, in beech, Bontemps et al. (2016) estimated narrow-sense heritability (*h*²) of 0.84–0.92 for SPSS (0.19–0.24 in
our study). In sessile and pedunculate oaks, budburst (leaf unfolding) had h²~0.85 and 0.81, and tree size (evaluated as the circumference at 1.30m) had h²~0.77 and 0.12 respectively (Alexandre, Truffaut, Ducousso, et al., 2020). The generally low heritability estimates in this study may be a result of increased residual variance, which encompasses (micro)environmental variance, G×E interactions and dominance, rather than reduced additive genetic variance, as previously shown in animals (Kruuk et al., 2000; Hansen et al., 2011; Wheelwright et al., 2014; see Table 4). Fine-grained spatial heterogeneity may generate various plastic responses across individuals, thus increasing the phenotypic variance in the population. Even when we controlled for microenvironmental variation with composite indices such as NCI and DBH, other covariates such as soil moisture content, which in beech is a key factor for development (Houston Durrant et al., 2016), if available, could have explained a larger proportion of the microenvironmental variance, leading to higher estimates of heritability.

To reduce the noise in heritability caused by environmental variance, Hansen et al. (2011) suggest calculating evolvability. Unlike heritability, evolvability provides a measure of standardized additive genetic variance and is, therefore, not confounded by environmental variation, allowing for direct comparison between traits, years and species (Hansen et al., 2011). The estimates of evolvability in our study were all below the median value ($I_A$=0.26) obtained in a review on plants and animals (Hansen et al., 2011), suggesting low to moderate genetic variation in the studied beech population (median over all traits of ~0.01). For example, we found $I_A$ ~0.01 for tree size (DBH) in our study, whereas $I_A$ was ~0.09 for size traits in Hansen et al. (2011) and ~0.02 for DBH in common gardens of pedunculate oak (Caignard et al., 2019). Interestingly, Alexandre, Truffaut, Ducousso, et al. (2020) reported $I_A$ ~0.02 for tree size (circumference at 1.3m) in the field in sessile oak, but ten-fold lower evolvability ($I_A$ ~0.0025) in co-occurring pedunculate oaks, which also showed lower genetic variation in other traits. The authors suggested that selection pressures may have depleted genetic variation of pedunculate oak in this forest, which is supported by its decreasing demographic trend. Thus, rather than increased residual variance, past selection leading to reduced additive genetic variance may also underlie the low to moderate heritability and evolvability for beech in our case study. Indeed, despite some conflicting evidence (e.g. in vertebrates, Wheelwright et al., 2014), quantitative genetic variation of fitness traits is generally negatively correlated with strength of selection (Kruuk et al., 2000; Mousseau & Roff, 1987; Wheelwright et al., 2014). The study area has been inhabited by beech since the last ice age (Brus, 2010; Šercelj, 1996). A favourable environment for beech growth may thus have led to strong competition and selection among individuals growing at high density, decreasing the genetic variance. For example, high reproductive output is likely to be selected for in such a mature forest as the study population, based on allocation theory (Bonser, 2013). Considering limited recruitment opportunities, this would have exacerbated competition at early offspring establishment. This is corroborated by the fact that even relatively small trees (i.e. DBH of 16cm) bore fruit in our study population, but further selection studies in the field are needed to support this hypothesis.

## 4.4 | Weak predicted selection responses

Using the breeder's equation, we were able to predict selection responses (R), by disentangling their ecological ($\beta_n$ and $\sigma_p$) and genetic (h²) components, assuming that the phenotypic value of the focal trait is the only target of selection (Walsh & Lynch, 2018). We found that natural selection primarily increased the length of budburst (average R across measurements of 0.21; due to relatively high heritability and selection gradients), decreased the length of the growing season (R=-0.085; due to relatively high heritability) and increased tree size (R=0.09; due to a relatively high selection gradient). However, responses to selection for most phenology and reproduction traits were predicted to be low, which can be explained by a combination of low heritability estimates, weak to moderate selection gradients and low phenotypic variance. Low phenotypic variance was particularly relevant for phenology traits, whose timing was synchronized within measurement years, resulting in low within-year variance (Table 3). Furthermore, some studies reported that phenological variation in beech may follow a counter-gradient (Gauzere et al., 2020; Vitasse et al., 2009), with phenotypic and genetic clines having different signs. In this case, the breeder's equation fails to correctly predict the expected change in trait mean (Morrissey et al., 2010; Pemberton, 2010). However, non-significant selection responses for phenological traits have also been reported in natural oak populations (Alexandre, Truffaut, Klein, et al., 2020), despite substantial genetic variation and co-gradient genetic clines for these traits (Alberto et al., 2013; Bontemps et al., 2016; Vitasse et al., 2009). The authors argued that this can be explained either by a low response to selection during the study period (e.g. no strong late frosts) or, more likely, by the absence of any selection, with co-gradients being explained by assortative mating (Alexandre, Truffaut, Klein, et al., 2020). Interestingly, compared to our study, Alexandre, Truffaut, Klein, et al. (2020) also reported a much higher predicted selection response for tree size (R=11.70 for circumference at 1.3m) in the expanding natural population of sessile oak, but not in the co-occurring and declining population of pedunculate oak. Biological factors related to specific population characteristics also seem a likely explanation of the contrasted selection response for tree size between sessile oak in Alexandre, Truffaut, Klein, et al. (2020) and our study.

Importantly, predicted selection responses based on phenotypic traits measured in a single or few years may not be representative of lifetime selection responses in long-lived forest trees such as our study population. As mentioned earlier, selection gradients in trees often fluctuate over years or across ontogenetic stages, especially for labile traits such as reproduction and phenological traits. Therefore, traits that appear to be under directional selection on short time scales may be subject to stabilizing selection on longer time scales, resulting in no response to selection in the long term. This is one of the main explanations for the so-called 'missing response to selection in the wild', that is, the absence of an observed response to selection despite evidence of genetic variation and selection pressures on a trait (Pujol et al., 2018). Other factors relevant
our study). In sessile and pedunculate oaks, budburst (leaf unfolding) had h²~0.85 and 0.81, and tree size (evaluated as the circumference at 1.30m) had h²~0.77 and 0.12 respectively (Alexandre, Truffaut, Ducousso, et al., 2020). The generally low heritability estimates in this study may be a result of increased residual variance, which encompasses (micro)environmental variance, G×E interactions and dominance, rather than reduced additive genetic variance, as previously shown in animals (Kruuk et al., 2000; Hansen et al., 2011; Wheelwright et al., 2014; see Table 4). Fine-grained spatial heterogeneity may generate various plastic responses across individuals, thus increasing the phenotypic variance in the population. Even when we controlled for microenvironmental variation with composite indices such as NCI and DBH, other covariates such as soil moisture content, which in beech is a key factor for development (Houston Durrant et al., 2016), if available, could have explained a larger proportion of the microenvironmental variance, leading to higher estimates of heritability.

To reduce the noise in heritability caused by environmental variance, Hansen et al. (2011) suggest calculating evolvability. Unlike heritability, evolvability provides a measure of standardized additive genetic variance and is, therefore, not confounded by environmental variation, allowing for direct comparison between traits, years and species (Hansen et al., 2011). The estimates of evolvability in our study were all below the median value ($I_A$=0.26) obtained in a review on plants and animals (Hansen et al., 2011), suggesting low to moderate genetic variation in the studied beech population (median over all traits of ~0.01). For example, we found $I_A$ ~0.01 for tree size (DBH) in our study, whereas $I_A$ was ~0.09 for size traits in Hansen et al. (2011) and ~0.02 for DBH in common gardens of pedunculate oak (Caignard et al., 2019). Interestingly, Alexandre, Truffaut, Ducousso, et al. (2020) reported $I_A$ ~0.02 for tree size (circumference at 1.3m) in the field in sessile oak, but ten-fold lower evolvability ($I_A$ ~0.0025) in co-occurring pedunculate oaks, which also showed lower genetic variation in other traits. The authors suggested that selection pressures may have depleted genetic variation of pedunculate oak in this forest, which is supported by its decreasing demographic trend. Thus, rather than increased residual variance, past selection leading to reduced additive genetic variance may also underlie the low to moderate heritability and evolvability for beech in our case study. Indeed, despite some conflicting evidence (e.g. in vertebrates, Wheelwright et al., 2014), quantitative genetic variation of fitness traits is generally negatively correlated with strength of selection (Kruuk et al., 2000; Mousseau & Roff, 1987; Wheelwright et al., 2014). The study area has been inhabited by beech since the last ice age (Brus, 2010; Šercelj, 1996). A favourable environment for beech growth may thus have led to strong competition and selection among individuals growing at high density, decreasing the genetic variance. For example, high reproductive output is likely to be selected for in such a mature forest as the study population, based on allocation theory (Bonser, 2013). Considering limited recruitment opportunities, this would have exacerbated competition at early offspring establishment. This is corroborated by the fact that even relatively small trees (i.e. DBH of 16cm) bore fruit in our study population, but further selection studies in the field are needed to support this hypothesis.

## 4.4 | Weak predicted selection responses

Using the breeder's equation, we were able to predict selection responses (R), by disentangling their ecological ($\beta_n$ and $\sigma_p$) and genetic (h²) components, assuming that the phenotypic value of the focal trait is the only target of selection (Walsh & Lynch, 2018). We found that natural selection primarily increased the length of budburst (average R across measurements of 0.21; due to relatively high heritability and selection gradients), decreased the length of the growing season (R=-0.085; due to relatively high heritability) and increased tree size (R=0.09; due to a relatively high selection gradient). However, responses to selection for most phenology and reproduction traits were predicted to be low, which can be explained by a combination of low heritability estimates, weak to moderate selection gradients and low phenotypic variance. Low phenotypic variance was particularly relevant for phenology traits, whose timing was synchronized within measurement years, resulting in low within-year variance (Table 3). Furthermore, some studies reported that phenological variation in beech may follow a counter-gradient (Gauzere et al., 2020; Vitasse et al., 2009), with phenotypic and genetic clines having different signs. In this case, the breeder's equation fails to correctly predict the expected change in trait mean (Morrissey et al., 2010; Pemberton, 2010). However, non-significant selection responses for phenological traits have also been reported in natural oak populations (Alexandre, Truffaut, Klein, et al., 2020), despite substantial genetic variation and co-gradient genetic clines for these traits (Alberto et al., 2013; Bontemps et al., 2016; Vitasse et al., 2009). The authors argued that this can be explained either by a low response to selection during the study period (e.g. no strong late frosts) or, more likely, by the absence of any selection, with co-gradients being explained by assortative mating (Alexandre, Truffaut, Klein, et al., 2020). Interestingly, compared to our study, Alexandre, Truffaut, Klein, et al. (2020) also reported a much higher predicted selection response for tree size (R=11.70 for circumference at 1.3m) in the expanding natural population of sessile oak, but not in the co-occurring and declining population of pedunculate oak. Biological factors related to specific population characteristics also seem a likely explanation of the contrasted selection response for tree size between sessile oak in Alexandre, Truffaut, Klein, et al. (2020) and our study.

Importantly, predicted selection responses based on phenotypic traits measured in a single or few years may not be representative of lifetime selection responses in long-lived forest trees such as our study population. As mentioned earlier, selection gradients in trees often fluctuate over years or across ontogenetic stages, especially for labile traits such as reproduction and phenological traits. Therefore, traits that appear to be under directional selection on short time scales may be subject to stabilizing selection on longer time scales, resulting in no response to selection in the long term. This is one of the main explanations for the so-called 'missing response to selection in the wild', that is, the absence of an observed response to selection despite evidence of genetic variation and selection pressures on a trait (Pujol et al., 2018). Other factors relevant
to this process arise both from statistical factors related to the measurement uncertainty of genetic parameters, especially additive genetic variance (including the uncertainty associated to pedigree reconstruction in natural populations), and from biological mechanisms (e.g. phenotypic plasticity, genetic correlations) that can modify selection pressures and thus prevent the loss of genetic variation induced by selection (Pujol et al., 2018).

## 4.5 | Conclusion

Our study shows that key parameters related to the potential of populations to respond to selection can be estimated in situ in forest trees. Importantly, our estimates are probably more realistic than those currently obtained in common gardens because we directly assessed selection pressures faced by trees in the wild. We were able to estimate a proxy of fitness (i.e. male and female effective fecundity) that is more accurate than survival or reproduction, and quantitative trait analyses were based on adult trees, which is often not possible in common gardens. Our study highlighted the fluctuating selection pressures forest trees face in the population studied and the resulting need to conduct long-term studies (e.g. Grant & Grant, 2014) that also take into account microenvironment variation (Denney et al., 2020). Forest trees are already facing reduced productivity, increased mortality and range shifts (Allen et al., 2010; McDowell et al., 2020). Therefore, long-term in situ studies are undoubtedly needed to anticipate the potential of natural populations across species' distribution ranges to respond to new selection pressures resulting from climate change. The generalization of genomic tools offers the possibility of applying similar methods to those used in this study to a growing number of forest trees, which is particularly valuable for planning appropriate conservation and management strategies.

## AUTHOR CONTRIBUTIONS

MW and SCG-M designed the research, MB and MW carried out laboratory analysis, MW, RD and SCG-M carried out formal analysis, MW, SO-M, JA and SCG-M wrote the manuscript while all authors critically reviewed it.

## ACKNOWLEDGEMENTS

This research was supported by LIFEGENMON (LIFE13 ENV/ SI/000148; sampling, trait measurements, genotyping), the Slovenian Research Agency (research core funding no. P4-0107; manuscript preparation), a Labex COTE mobility grant (stay of MW at BIOGECO for data analysis) and a G-BIKE COST Action STSM (stay of SCGM at SFI for data analysis and manuscript preparation). This study was also funded by the European Union's Horizon 2020 research and innovation programme under grant agreement No 862221 (FORGENIUS). We thank Antoine Kremer and Hermine Alexandre from INRAE and Ricardo Alía from INIA-CSIC for discussions that greatly assisted the research, Barbara Štupar for DNA extraction and Boris Rantaša from SFI for his help with figure preparation. The drawings of beech

phenology, seed and habitus in Figures 1 and 2 are reproduced from Westergren et al. (2020) with permission of the publisher.

## CONFLICT OF INTEREST STATEMENT

The authors have no conflicts of interest.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are openly available in DataDryad with DOI: http://doi.org/10.5061/dryad.x0k6djhrn (Westergren et al., 2023).

## ORCID

Marjana Westergren https://orcid.org/0000-0002-4204-0161

Juliette Archambeau https://orcid.org/0000-0003-0155-9538

Santiago C. González-Martínez https://orcid.org/0000-0002-4534-3766

## REFERENCES

Aitken, S., & Whitlock, M. C. (2013). Assisted gene flow to facilitate local adaptation to climate change. *Annual Review of Ecology, Evolution, and Systematics*, *44*, 367–388. https://doi.org/10.1146/annurev-ecolsys-110512-135747

Alberto, F., Bouffier, L., Louvet, J.-M., Lamy, J.-B., Delzon, S., & Kremer, A. (2011). Adaptive responses for seed and leaf phenology in natural populations of sessile oak along an altitudinal gradient: Variation of phenological traits in *Q. petraea*. *Journal of Evolutionary Biology*, *24*(7), 1442–1454. https://doi.org/10.1111/j.1420-9101.2011.02277.x

Alberto, F. J., Aitken, S. N., Alía, R., González-Martínez, S. C., Hänninen, H., Kremer, A., Lefèvre, F., Lenormand, T., Yeaman, S., Whetten, R., & Savolainen, O. (2013). Potential for evolutionary responses to climate change—Evidence from tree populations. *Global Change Biology*, *19*(6), 1645–1661. https://doi.org/10.1111/gcb.12181

Alexandre, H., Truffaut, L., Ducousso, A., Louvet, J.-M., Nepveu, G., Torres-Ruiz, J. M., Lagane, F., Firmat, C., Musch, B., Delzon, S., & Kremer, A. (2020). In situ estimation of genetic variation of functional and ecological traits in *Quercus petraea* and *Q. robur*. *Tree Genetics & Genomes*, *16*(2), 32. https://doi.org/10.1007/s11295-019-1407-9

Alexandre, H., Truffaut, L., Klein, E., Ducousso, A., Chancerel, E., Lesur, I., Dencausse, B., Louvet, J., Nepveu, G., Torres-Ruiz, J. M., Lagane, F., Musch, B., Delzon, S., & Kremer, A. (2020). How does contemporary selection shape oak phenotypes? *Evolutionary Applications*, *13*(10), 2772–2790. https://doi.org/10.1111/eva.13082

Alía, R., Chambel, R., Notivol, E., Climent, J., & González-Martínez, S. C. (2014). Environment-dependent microevolution in a Mediterranean pine (*Pinus pinaster* Aiton). *BMC Evolutionary Biology*, *14*(1), 200. https://doi.org/10.1186/s12862-014-0200-5

Allen, C. D., Macalady, A. K., Chenchouni, H., Bachelet, D., McDowell, N., Vennetier, M., Kitzberger, T., Rigling, A., Breshears, D. D., Hogg, E. H. (T.), Gonzalez, P., Fensham, R., Zhang, Z., Castro, J., Demidova, N., Lim, J.-H., Allard, G., Running, S. W., Semerci, A., & Cobb, N. (2010). A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. *Forest Ecology and Management*, *259*(4), 660–684. https://doi.org/10.1016/j.foreco.2009.09.001

Andrew, R. L., Peakall, R., Wallis, I. R., Wood, J. T., Knight, E. J., & Foley, W. J. (2005). Marker-based quantitative genetics in the wild?: The heritability and genetic correlation of chemical defenses in
to this process arise both from statistical factors related to the measurement uncertainty of genetic parameters, especially additive genetic variance (including the uncertainty associated to pedigree reconstruction in natural populations), and from biological mechanisms (e.g. phenotypic plasticity, genetic correlations) that can modify selection pressures and thus prevent the loss of genetic variation induced by selection (Pujol et al., 2018).

## 4.5 | Conclusion

Our study shows that key parameters related to the potential of populations to respond to selection can be estimated in situ in forest trees. Importantly, our estimates are probably more realistic than those currently obtained in common gardens because we directly assessed selection pressures faced by trees in the wild. We were able to estimate a proxy of fitness (i.e. male and female effective fecundity) that is more accurate than survival or reproduction, and quantitative trait analyses were based on adult trees, which is often not possible in common gardens. Our study highlighted the fluctuating selection pressures forest trees face in the population studied and the resulting need to conduct long-term studies (e.g. Grant & Grant, 2014) that also take into account microenvironment variation (Denney et al., 2020). Forest trees are already facing reduced productivity, increased mortality and range shifts (Allen et al., 2010; McDowell et al., 2020). Therefore, long-term in situ studies are undoubtedly needed to anticipate the potential of natural populations across species' distribution ranges to respond to new selection pressures resulting from climate change. The generalization of genomic tools offers the possibility of applying similar methods to those used in this study to a growing number of forest trees, which is particularly valuable for planning appropriate conservation and management strategies.

## AUTHOR CONTRIBUTIONS

MW and SCG-M designed the research, MB and MW carried out laboratory analysis, MW, RD and SCG-M carried out formal analysis, MW, SO-M, JA and SCG-M wrote the manuscript while all authors critically reviewed it.

## ACKNOWLEDGEMENTS

This research was supported by LIFEGENMON (LIFE13 ENV/ SI/000148; sampling, trait measurements, genotyping), the Slovenian Research Agency (research core funding no. P4-0107; manuscript preparation), a Labex COTE mobility grant (stay of MW at BIOGECO for data analysis) and a G-BIKE COST Action STSM (stay of SCGM at SFI for data analysis and manuscript preparation). This study was also funded by the European Union's Horizon 2020 research and innovation programme under grant agreement No 862221 (FORGENIUS). We thank Antoine Kremer and Hermine Alexandre from INRAE and Ricardo Alía from INIA-CSIC for discussions that greatly assisted the research, Barbara Štupar for DNA extraction and Boris Rantaša from SFI for his help with figure preparation. The drawings of beech

phenology, seed and habitus in Figures 1 and 2 are reproduced from Westergren et al. (2020) with permission of the publisher.

## CONFLICT OF INTEREST STATEMENT

The authors have no conflicts of interest.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are openly available in DataDryad with DOI: http://doi.org/10.5061/dryad.x0k6djhrn (Westergren et al., 2023).

## ORCID

Marjana Westergren https://orcid.org/0000-0002-4204-0161

Juliette Archambeau https://orcid.org/0000-0003-0155-9538

Santiago C. González-Martínez https://orcid.org/0000-0002-4534-3766

## REFERENCES

Aitken, S., & Whitlock, M. C. (2013). Assisted gene flow to facilitate local adaptation to climate change. *Annual Review of Ecology, Evolution, and Systematics*, *44*, 367–388. https://doi.org/10.1146/annurev-ecolsys-110512-135747

Alberto, F., Bouffier, L., Louvet, J.-M., Lamy, J.-B., Delzon, S., & Kremer, A. (2011). Adaptive responses for seed and leaf phenology in natural populations of sessile oak along an altitudinal gradient: Variation of phenological traits in *Q. petraea*. *Journal of Evolutionary Biology*, *24*(7), 1442–1454. https://doi.org/10.1111/j.1420-9101.2011.02277.x

Alberto, F. J., Aitken, S. N., Alía, R., González-Martínez, S. C., Hänninen, H., Kremer, A., Lefèvre, F., Lenormand, T., Yeaman, S., Whetten, R., & Savolainen, O. (2013). Potential for evolutionary responses to climate change—Evidence from tree populations. *Global Change Biology*, *19*(6), 1645–1661. https://doi.org/10.1111/gcb.12181

Alexandre, H., Truffaut, L., Ducousso, A., Louvet, J.-M., Nepveu, G., Torres-Ruiz, J. M., Lagane, F., Firmat, C., Musch, B., Delzon, S., & Kremer, A. (2020). In situ estimation of genetic variation of functional and ecological traits in *Quercus petraea* and *Q. robur*. *Tree Genetics & Genomes*, *16*(2), 32. https://doi.org/10.1007/s11295-019-1407-9

Alexandre, H., Truffaut, L., Klein, E., Ducousso, A., Chancerel, E., Lesur, I., Dencausse, B., Louvet, J., Nepveu, G., Torres-Ruiz, J. M., Lagane, F., Musch, B., Delzon, S., & Kremer, A. (2020). How does contemporary selection shape oak phenotypes? *Evolutionary Applications*, *13*(10), 2772–2790. https://doi.org/10.1111/eva.13082

Alía, R., Chambel, R., Notivol, E., Climent, J., & González-Martínez, S. C. (2014). Environment-dependent microevolution in a Mediterranean pine (*Pinus pinaster* Aiton). *BMC Evolutionary Biology*, *14*(1), 200. https://doi.org/10.1186/s12862-014-0200-5

Allen, C. D., Macalady, A. K., Chenchouni, H., Bachelet, D., McDowell, N., Vennetier, M., Kitzberger, T., Rigling, A., Breshears, D. D., Hogg, E. H. (T.), Gonzalez, P., Fensham, R., Zhang, Z., Castro, J., Demidova, N., Lim, J.-H., Allard, G., Running, S. W., Semerci, A., & Cobb, N. (2010). A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. *Forest Ecology and Management*, *259*(4), 660–684. https://doi.org/10.1016/j.foreco.2009.09.001

Andrew, R. L., Peakall, R., Wallis, I. R., Wood, J. T., Knight, E. J., & Foley, W. J. (2005). Marker-based quantitative genetics in the wild?: The heritability and genetic correlation of chemical defenses in
Eucalyptus. *Genetics*, 171(4), 1989–1998. https://doi.org/10.1534/genetics.105.042952

Arntz, M. A., & Delph, L. F. (2001). Pattern and process: Evidence for the evolution of photosynthetic traits in natural populations. *Oecologia*, 127(4), 455–467. https://doi.org/10.1007/s004420100650

Asuka, Y., Tomaru, N., Munehara, Y., Tani, N., Tsumura, Y., & Yamamoto, S. (2005). Half-sib family structure of *Fagus crenata* saplings in an old-growth beech-dwarf bamboo forest. *Molecular Ecology*, 14(8), 2565–2575. https://doi.org/10.1111/j.1365-294x.2005.02585.x

Augspurger, C. K. (2013). Reconstructing patterns of temperature, phe-nology, and frost damage over 124 years: Spring damage risk is in-creasing. *Ecology*, 94(1), 41–50. https://doi.org/10.1890/12-0200.1

Bajc, M., Aravanopoulos, F. A., Westergren, M., Fussi, B., Kavaliauskas, D., Alizoti, P., Kiourtsis, F., & Kraigher, H. (2021). *Manual for forest genetic monitoring* (p. 326). Slovenian Forestry Institute. https://doi.org/10.20315/SFS.167

Baliuckas, V., & Pliura, A. (2003). Genetic variation and phenotypic plasticity of *Quercus robur* populations and open-pollinated families in Lithuania. *Scandinavian Journal of Forest Research*, 18(4), 305–319. https://doi.org/10.1080/02827580310005153

Barringer, B. C., Koenig, W. D., & Knops, J. M. H. (2013). Interrelationships among life-history traits in three California oaks. *Oecologia*, 171(1), 129–139. https://doi.org/10.1007/s00442-012-2386-9

Bennett, A. C., McDowell, N. G., Allen, C. D., & Anderson-Teixeira, K. J. (2015). Larger trees suffer most during drought in forests world-wide. *Nature Plants*, 1(10), 15139. https://doi.org/10.1038/nplants.2015.139

Bigler, C., & Bugmann, H. (2018). Climate-induced shifts in leaf unfolding and frost risk of European trees and shrubs. *Scientific Reports*, 8(1), 9865. https://doi.org/10.1038/s41598-018-27893-1

Bilir, N., Prescher, F., Ayan, S., & Lindgren, D. (2006). Growth charac-ters and number of strobili in clonal seed orchards of *Pinus syl-vestris*. *Euphytica*, 152(2), 1–9. https://doi.org/10.1007/s10681-006-9216-2

Bloom, A. J., Chapin, F. S., & Mooney, H. A. (1985). Resource limitation in plants—An economic analogy. *Annual Review of Ecology and Systematics*, 92(16), 363.

Böhlenius, H., Huang, T., Charbonnel-Campaa, L., Brunner, A. M., Jansson, S., Strauss, S. H., & Nilsson, O. (2006). CO/FT regulatory module controls timing of flowering and seasonal growth cessation in trees. *Science*, 312(5776), 1040–1043. https://doi.org/10.1126/science.1126038

Bonser, S. P. (2013). High reproductive efficiency as an adaptive strategy in competitive environments. *Functional Ecology*, 27(4), 876–885. https://doi.org/10.1111/1365-2435.12064

Bontemps, A., Lefèvre, F., Davi, H., & Oddou-Muratorio, S. (2016). In situ marker-based assessment of leaf trait evolutionary potential in a marginal European beech population. *Journal of Evolutionary Biology*, 29(3), 514–527. https://doi.org/10.1111/jeb.12801

Brus, R. (2010). Growing evidence for the existence of glacial refugia of European beech (*Fagus sylvatica* L.) in the south-eastern Alps and north-western Dinaric Alps. *Periodicum Biologorum*, 112(3), 239–246.

Caignard, T., Delzon, S., Bodénès, C., Dencausse, B., & Kremer, A. (2019). Heritability and genetic architecture of reproduction-related traits in a temperate oak species. *Tree Genetics & Genomes*, 15(1), 1. https://doi.org/10.1007/s11295-018-1309-2

Caignard, T., Kremer, A., Firmat, C., Nicolas, M., Venner, S., & Delzon, S. (2017). Increasing spring temperatures favor oak seed production in temperate areas. *Scientific Reports*, 7(1), 8555. https://doi.org/10.1038/s41598-017-09172-7

Canham, C. D., LePage, P. T., & Coates, K. D. (2004). A neighborhood analysis of canopy tree competition: Effects of shading versus crowding. *Canadian Journal of Forest Research*, 34(4), 778–787. https://doi.org/10.1139/x03-232

Castellanos, M. C., González-Martínez, S. C., & Pausas, J. G. (2015). Field heritability of a plant adaptation to fire in heterogeneous land-scapes. *Molecular Ecology*, 24(22), 5633–5642. https://doi.org/10.1111/mec.13421

Christophe, C., & Birot, Y. (1979). Genetic variation within and between populations of douglas-fir. *Silvae Genetica*, 28(5–8), 197–206.

Chybicki, I. J., & Burczyk, J. (2010). Realized gene flow within mixed stands of *Quercus robur* L. and *Q. petraea* (Matt.) L. revealed at the stage of naturally established seedling. *Molecular Ecology*, 19(10), 2137–2151. https://doi.org/10.1111/j.1365-294X.2010.04632.x

Conner, J. K., Rush, S., & Jennetten, P. (1996). Measurements of natural selection on floral traits in wild radish (*Raphanus raphanistrum*). I. Selection through lifetime female fitness. *Evolution*, 50(3), 1127–1136. https://doi.org/10.1111/j.1558-5646.1996.tb02353.x

Cornelius, J. (1994). Heritabilities and additive genetic coefficients of variation in forest trees. *Canadian Journal of Forest Research*, 24(2), 372–379. https://doi.org/10.1139/x94-050

Davis, M. B., Shaw, R. G., & Ettersp, J. R. (2005). Evolutionary responses to changing climate. *Ecology*, 86(7), 1704–1714. https://doi.org/10.1890/03-0788

de Miguel, M., Rodríguez-Quilón, I., Heuertz, M., Hurel, A., Grivet, D., Jaramillo-Correa, J. P., Vendramin, G. G., Plomion, C., Majada, J., Alía, R., Eckert, A. J., & González-Martínez, S. C. (2022). Polygenic adaptation and negative selection across traits, years and environ-ments in a long-lived plant species (*Pinus pinaster* Ait., Pinaceae). *Molecular Ecology*, 31(7), 2089–2105. https://doi.org/10.1111/mec.16367

de Villemereuil, P., Wells, J. A., Edwards, R. D., & Blomberg, S. P. (2012). Bayesian models for comparative analysis integrating phylogenetic uncertainty. *BMC Evolutionary Biology*, 12(1), 102. https://doi.org/10.1186/1471-2148-12-102

Denney, D. A., Jameel, M. I., Bemmels, J. B., Rochford, M. E., & Anderson, J. T. (2020). Small spaces, big impacts: Contributions of micro-environmental variation to population persistence under climate change. *AoB Plants*, 12(2), plaa005. https://doi.org/10.1093/aobpla/pla0005

Donohue, K., Rubio de Casas, R., Burghardt, L., Kovach, K., & Willis, C. G. (2010). Germination, postgermination adaptation, and species eco-logical ranges. *Annual Review of Ecology, Evolution, and Systematics*, 41(1), 293–319. https://doi.org/10.1146/annurev-ecolsys-102209-9-144715

Elle, E., & Meagher, T. R. (2000). Sex allocation and reproductive suc-cess in the andromonoecious perennial *Solanum carolinense* (Solanaceae). II. Paternity and functional gender. *The American Naturalist*, 156(6), 622–636. https://doi.org/10.1086/316997

Ellison, A. M. (2019). Foundation species, non-trophic interactions, and the value of being common. *IScience*, 13, 254–268. https://doi.org/10.1016/j.isci.2019.02.020

Fernández, J., & Toro, M. A. (2006). A new method to estimate relatedness from molecular markers. *Molecular Ecology*, 15(6), 1657–1667. https://doi.org/10.1111/j.1365-294X.2006.02873.x

Gauzere, J., Klein, E. K., Brendel, O., Davi, H., & Oddou-Muratorio, S. (2016). Using partial genotyping to estimate the genetic and ma-terial determinants of adaptive traits in a progeny trial of *Fagus sylvatica*. *Tree Genetics & Genomes*, 12(6), 115. https://doi.org/10.1007/s11295-016-1062-3

Gauzere, J., Teuf, B., Davi, H., Chevin, L., Caignard, T., Leys, B., Delzon, S., Ronce, O., & Chuine, I. (2020). Where is the optimum? Predicting the variation of selection along climatic gradients and the adap-tive value of plasticity. A case study on tree phenology. *Evolution Letters*, 4(2), 109–123. https://doi.org/10.1002/evl3.160

Geber, M. A., & Griffen, L. R. (2003). Inheritance and natural selection on functional traits. *International Journal of Plant Sciences*, 164(S3), S21–S42. https://doi.org/10.1086/368233
Eucalyptus. *Genetics*, 171(4), 1989–1998. https://doi.org/10.1534/genetics.105.042952

Arntz, M. A., & Delph, L. F. (2001). Pattern and process: Evidence for the evolution of photosynthetic traits in natural populations. *Oecologia*, 127(4), 455–467. https://doi.org/10.1007/s004420100650

Asuka, Y., Tomaru, N., Munehara, Y., Tani, N., Tsumura, Y., & Yamamoto, S. (2005). Half-sib family structure of *Fagus crenata* saplings in an old-growth beech-dwarf bamboo forest. *Molecular Ecology*, 14(8), 2565–2575. https://doi.org/10.1111/j.1365-294x.2005.02585.x

Augspurger, C. K. (2013). Reconstructing patterns of temperature, phe-nology, and frost damage over 124 years: Spring damage risk is in-creasing. *Ecology*, 94(1), 41–50. https://doi.org/10.1890/12-0200.1

Bajc, M., Aravanopoulos, F. A., Westergren, M., Fussi, B., Kavaliauskas, D., Alizoti, P., Kiourtsis, F., & Kraigher, H. (2021). *Manual for forest genetic monitoring* (p. 326). Slovenian Forestry Institute. https://doi.org/10.20315/SFS.167

Baliuckas, V., & Pliura, A. (2003). Genetic variation and phenotypic plas-ticity of *Quercus robur* populations and open-pollinated families in Lithuania. *Scandinavian Journal of Forest Research*, 18(4), 305–319. https://doi.org/10.1080/02827580310005153

Barringer, B. C., Koenig, W. D., & Knops, J. M. H. (2013). Interrelationships among life-history traits in three California oaks. *Oecologia*, 171(1), 129–139. https://doi.org/10.1007/s00442-012-2386-9

Bennett, A. C., McDowell, N. G., Allen, C. D., & Anderson-Teixeira, K. J. (2015). Larger trees suffer most during drought in forests world-wide. *Nature Plants*, 1(10), 15139. https://doi.org/10.1038/nplants.2015.139

Bigler, C., & Bugmann, H. (2018). Climate-induced shifts in leaf unfolding and frost risk of European trees and shrubs. *Scientific Reports*, 8(1), 9865. https://doi.org/10.1038/s41598-018-27893-1

Bilir, N., Prescher, F., Ayan, S., & Lindgren, D. (2006). Growth charac-ters and number of strobili in clonal seed orchards of *Pinus syl-vestris*. *Euphytica*, 152(2), 1–9. https://doi.org/10.1007/s10681-006-9216-2

Bloom, A. J., Chapin, F. S., & Mooney, H. A. (1985). Resource limitation in plants—An economic analogy. *Annual Review of Ecology and Systematics*, 92(16), 363.

Böhlenius, H., Huang, T., Charbonnel-Campaa, L., Brunner, A. M., Jansson, S., Strauss, S. H., & Nilsson, O. (2006). CO/FT regulatory module controls timing of flowering and seasonal growth cessation in trees. *Science*, 312(5776), 1040–1043. https://doi.org/10.1126/science.1126038

Bonser, S. P. (2013). High reproductive efficiency as an adaptive strategy in competitive environments. *Functional Ecology*, 27(4), 876–885. https://doi.org/10.1111/1365-2435.12064

Bontemps, A., Lefèvre, F., Davi, H., & Oddou-Muratorio, S. (2016). In situ marker-based assessment of leaf trait evolutionary potential in a marginal European beech population. *Journal of Evolutionary Biology*, 29(3), 514–527. https://doi.org/10.1111/jeb.12801

Brus, R. (2010). Growing evidence for the existence of glacial refugia of European beech (*Fagus sylvatica* L.) in the south-eastern Alps and north-western Dinaric Alps. *Periodicum Biologorum*, 112(3), 239–246.

Caignard, T., Delzon, S., Bodénès, C., Dencausse, B., & Kremer, A. (2019). Heritability and genetic architecture of reproduction-related traits in a temperate oak species. *Tree Genetics & Genomes*, 15(1), 1. https://doi.org/10.1007/s11295-018-1309-2

Caignard, T., Kremer, A., Firmat, C., Nicolas, M., Venner, S., & Delzon, S. (2017). Increasing spring temperatures favor oak seed production in temperate areas. *Scientific Reports*, 7(1), 8555. https://doi.org/10.1038/s41598-017-09172-7

Canham, C. D., LePage, P. T., & Coates, K. D. (2004). A neighborhood analysis of canopy tree competition: Effects of shading versus crowding. *Canadian Journal of Forest Research*, 34(4), 778–787. https://doi.org/10.1139/x03-232

Castellanos, M. C., González-Martínez, S. C., & Pausas, J. G. (2015). Field heritability of a plant adaptation to fire in heterogeneous land-scapes. *Molecular Ecology*, 24(22), 5633–5642. https://doi.org/10.1111/mec.13421

Christophe, C., & Birot, Y. (1979). Genetic variation within and between populations of douglas-fir. *Silvae Genetica*, 28(5–8), 197–206.

Chybicki, I. J., & Burczyk, J. (2010). Realized gene flow within mixed stands of *Quercus robur* L. and *Q. petraea* (Matt.) L. revealed at the stage of naturally established seedling. *Molecular Ecology*, 19(10), 2137–2151. https://doi.org/10.1111/j.1365-294X.2010.04632.x

Conner, J. K., Rush, S., & Jennetten, P. (1996). Measurements of natural selection on floral traits in wild radish (*Raphanus raphanistrum*). I. Selection through lifetime female fitness. *Evolution*, 50(3), 1127–1136. https://doi.org/10.1111/j.1558-5646.1996.tb02353.x

Cornelius, J. (1994). Heritabilities and additive genetic coefficients of variation in forest trees. *Canadian Journal of Forest Research*, 24(2), 372–379. https://doi.org/10.1139/x94-050

Davis, M. B., Shaw, R. G., & Ettersp, J. R. (2005). Evolutionary responses to changing climate. *Ecology*, 86(7), 1704–1714. https://doi.org/10.1890/03-0788

de Miguel, M., Rodríguez-Quilón, I., Heuertz, M., Hurel, A., Grivet, D., Jaramillo-Correa, J. P., Vendramin, G. G., Plomion, C., Majada, J., Alía, R., Eckert, A. J., & González-Martínez, S. C. (2022). Polygenic adaptation and negative selection across traits, years and environ-ments in a long-lived plant species (*Pinus pinaster* Ait., Pinaceae). *Molecular Ecology*, 31(7), 2089–2105. https://doi.org/10.1111/mec.16367

de Villemereuil, P., Wells, J. A., Edwards, R. D., & Blomberg, S. P. (2012). Bayesian models for comparative analysis integrating phylogenetic uncertainty. *BMC Evolutionary Biology*, 12(1), 102. https://doi.org/10.1186/1471-2148-12-102

Denney, D. A., Jameel, M. I., Bemmels, J. B., Rochford, M. E., & Anderson, J. T. (2020). Small spaces, big impacts: Contributions of micro-environmental variation to population persistence under climate change. *AoB Plants*, 12(2), plaa005. https://doi.org/10.1093/aobpla/plaa005

Donohue, K., Rubio de Casas, R., Burghardt, L., Kovach, K., & Willis, C. G. (2010). Germination, postgermination adaptation, and species eco-logical ranges. *Annual Review of Ecology, Evolution, and Systematics*, 41(1), 293–319. https://doi.org/10.1146/annurev-ecolsys-102209-9-144715

Elle, E., & Meagher, T. R. (2000). Sex allocation and reproductive suc-cess in the andromonoecious perennial *Solanum carolinense* (Solanaceae). II. Paternity and functional gender. *The American Naturalist*, 156(6), 622–636. https://doi.org/10.1086/316997

Ellison, A. M. (2019). Foundation species, non-trophic interactions, and the value of being common. *IScience*, 13, 254–268. https://doi.org/10.1016/j.isci.2019.02.020

Fernández, J., & Toro, M. A. (2006). A new method to estimate relatedness from molecular markers. *Molecular Ecology*, 15(6), 1657–1667. https://doi.org/10.1111/j.1365-294X.2006.02873.x

Gauzere, J., Klein, E. K., Brendel, O., Davi, H., & Oddou-Muratorio, S. (2016). Using partial genotyping to estimate the genetic and ma-terial determinants of adaptive traits in a progeny trial of *Fagus sylvatica*. *Tree Genetics & Genomes*, 12(6), 115. https://doi.org/10.1007/s11295-016-1062-3

Gauzere, J., Teuf, B., Davi, H., Chevin, L., Caignard, T., Leys, B., Delzon, S., Ronce, O., & Chuine, I. (2020). Where is the optimum? Predicting the variation of selection along climatic gradients and the adap-tive value of plasticity. A case study on tree phenology. *Evolution Letters*, 4(2), 109–123. https://doi.org/10.1002/evl3.160

Geber, M. A., & Griffen, L. R. (2003). Inheritance and natural selection on functional traits. *International Journal of Plant Sciences*, 164(S3), S21–S42. https://doi.org/10.1086/368233