# Evolutionary rescue in a mixed beech-fir forest: insights from a quantitative-genetics approach in a process-based model

Louis Devresse, Freya Way, Tanguy Postic, François de Coligny, Xavier Morin

► To cite this version:

Louis Devresse, Freya Way, Tanguy Postic, François de Coligny, Xavier Morin. Evolutionary rescue in a mixed beech-fir forest: insights from a quantitative-genetics approach in a process-based model. Oikos, In press, 10.1002/oik.11661. hal-05150321

HAL Id: hal-05150321

https://hal.inrae.fr/hal-05150321v1

Submitted on 8 Jul 2025

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_2.png" /></figure>

OIKOS

Research article

# Evolutionary rescue in a mixed beech-fir forest: insights from a quantitative-genetics approach in a process-based model

Louis Devresse ID,✉<fnref n="1" />, Freya Way ID<fnref n="2" />, Tanguy Postic ID<fnref n="1" />, François de Coligny ID<fnref n="3" /> and Xavier Morin ID,✉<fnref n="1" />

**Correspondence:** Xavier Morin (xavier.morin@cefe.cnrs.fr)

Oikos

2025: e11661

doi: 10.1002/oik.11661

Subject Editor: Dries Bonte

Editor-in-Chief: Pedro Peres-Neto

Accepted 17 May 2025

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

Questions have been raised about the ability of long-lived organisms, such as trees, to adapt to rapid climate change, and to what extent forest management actions influence the evolutionary responses of tree species. Given the life history of trees and the time scales involved, these questions are often addressed through modeling approaches. Yet, most of these studies focus on single-species case studies. The main objective and originality of our work is to explore the evolutionary responses of tree species to climate change using a process-based model, in a multi-specific context. This approach allows us to investigate the conditions necessary for evolutionary rescue in a mixed beech-fir forest. Furthermore, we explored how climate change adaptation and mitigation solutions, such as assisted gene flow and assisted migration, affect the conditions for evolutionary rescue in this forest type. To achieve these objectives, we integrated a quantitative genetic module into a process-based forest gap model, enabling species-specific parameters to evolve as quantitative traits under selective pressure and drift. Our results show that increased trait variability and heritability reduce the loss of forest cover following climatic warming in the short term (over a century). We also found that assisted gene flow had the expected effect of aiding species adapt to climate change. Finally, our study suggests that introducing new pre-adapted species into the forest could improve recovery after climate change but could also hinder the evolutionary rescue of local species. We conclude that integrating evolutionary dynamics into process-based models significantly enhances their predictive power by incorporating genetic adaptation scenarios that would otherwise be overlooked. This approach also allows us to test eco-evolutionary hypotheses and better understand the potential consequences of adaptation measures to climate change for tree species.

Keywords: climate change, eco-evolution, ecology, evolution, forest ecology, modeling, process-based modeling, quantitative genetic

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

<footnote n="1">¹CEFE, University of Montpellier, CNRS, EPHE, IRD, Montpellier, France</footnote>

<footnote n="2">²University of Edinburgh, Edinburgh, UK</footnote>

<footnote n="3">³botAnique et Modélisation de l'Architecture des Plantes et des végétations AMAP, Montferrier-sur-Lez, France</footnote>
Introduction

There is growing concern about the viability of temper-
ate forests in the coming years given the speed at which
climatic conditions are expected to change (Dyderski et al.
2018, Wessely et al. 2024). At the global scale, a large num-
ber of forests have already experienced increased mortality
induced by heat waves and extreme drought events related to
climate change (Hartmann et al. 2022). Among the emerg-
ing adaptation solutions to climate change, the promotion
of species diversity and forest structure complexity is gaining
ground, supported by a growing body of empirical evidence
(Hisano et al. 2018, Depauw et al. 2024). More precisely,
increasing evidence suggests that tree species diversity can
stabilize or even promote forest growth in a changing cli-
mate – according to observations (Pardos et al. 2021) and
modeling work (Morin et al. 2018) – and also buffer the
impact of disturbances in communities (Silva Pedro et al.
2015). Higher diversity may also raise the probability that
a forest will contain species pre-adapted to future environ-
mental conditions through a sampling effect (Loreau and
Hector 2001). However, forecasting the long-term trajectory
of forest ecosystems in a changing environment requires con-
sideration of the evolutionary dynamics of tree species. In
fact, the ability of long-lived species such as trees to adapt
to rapid climate change, or to simply avoid local extinctions
(Alberto et al. 2013) is often questioned. The evolutionary
rescue is defined by a population's ability to avoid extinction
(persistence) and recover through natural selection acting on
heritable variation (Bell 2017). Recovery refers here to the
population's capacity to regain a viable size following a distur-
bance. Indeed, the prospect of evolutionary rescue of tree spe-
cies in the face of climate change remains largely unknown, as
well as the effect of management actions on the evolutionary
responses of tree species (Lefèvre et al. 2014). In addition,
evolution and species diversity can interact in complex ways,
stressing the importance of taking eco-evolutionary feedback
loops into account when studying the fate of forests. For
example, in species-rich communities, certain species could
be prevented from adapting to a new environment because
of competitive exclusion from other pre-adapted species
(Johansson 2008, Price and Kirkpatrick 2009, Thompson
and Fronhofer 2019). This effect is even more important for
rare species, for which a reduction in population size due
to competition can have particularly harmful consequences
on their adaptive capacity, by increasing genetic drift (De
Mazancourt et al. 2008, Eldijk et al. 2020). Lastly, the selec-
tion pressures exerted by interactions between species can act
in the same or in the opposite directions to those imposed by
environmental changes, thereby promoting or hindering spe-
cies adaptation (Osmond and de Mazancourt 2013).

Given the difficulty of experimenting on mature forests
and integrating the effects of future climate shifts, these
questions are frequently investigated through modeling
approaches. Among these, process-based models (PBMs)
have been increasingly promoted in recent years. We define
PBMs as simplified representations of a system, characterized

by functional components that interact with each other and
with the environment through explicitly simulated physical
and mechanistic processes over time (Mäkelä et al. 2000).
In comparison with statistical or phenomenological models,
PBMs should be more robust in predicting outcomes under
unprecedented conditions, as they do not need to extrapo-
late trends based on observed statistical relationships between
variables (Cuddington et al. 2013, Van der Meersch et al.
2025). For instance, PBMs have already been used to sim-
ulate species range shifts (Morin et al. 2008, Keenan et al.
2011), changes in ecosystem productivity (Guillemot et al.
2017, Morin et al. 2018) and responses to drought events
(Ruffault et al. 2022). In addition, the development of PBMs
often involves the addition of new functional components,
with their associated parameters and equations. The devel-
opment of PBMs has therefore been limited by the gaps in
ecological knowledge to implement and parameterize new
functions, as well as (although to a lesser extent) comput-
ing power and code architecture. This is even more true for
finer-scale models such as the individual-based models often
used to simulate forest dynamics, because such models are
often parameterized for the same handful of dominant spe-
cies and focus mainly on monospecific and even-aged forests
(discussed by Morin et al. 2021). For many of the same rea-
sons, evolutionary processes are still rarely implemented in
PBMs, as they would require the simulation of new processes,
such as the transmission of heritable traits or alleles from par-
ents to offspring. This involves parameters that are costly and
difficult to estimate in field studies such as the heritabilities
and intraspecific variances of traits. As a result, most predic-
tions of the future of tree species have been made using fixed
responses, without incorporating evolutionary adaptation.
However, authors who have incorporated evolution into
forest PBMs have shown, for example, that species-specific
parameters that are subjected to evolution can change sig-
nificantly under the pressure of natural selection, after only
two to six generations (Oddou-Muratorio and Davi 2014,
Kramer et al. 2015). Furthermore, modeling studies suggest
that forest management practices such as early thinning can
hinder adaptation by substituting the self-thinning process,
which naturally selects better-adapted allele combinations
over time (Fririon et al. 2024). Consequently, ignoring evo-
lution can lead to inappropriate climate change adaptation
and mitigation measures, and can result in an overestimation
of the negative consequences of climate change on forests
(Gloy et al. 2023).

However, studies that have integrated evolution into PBMs
so far have limited their scope to monospecific and often
even-aged forests. For this reason we propose here a novel
approach, that involves considering evolutionary dynamics in
a multi-specific context in order to improve the robustness of
PBM predictions over long time scales. To this end, we first
developed an evolutionary module, capable of simulating the
evolution of species-specific parameters of a PBM through
natural selection and drift in diverse and uneven-aged forests.
This module was added to the ForCEEPS (Forest Community
Ecology and Ecosystem Processes) PBM (see Material and
methods for a model description and Morin et al. 2020, 2021). Among the parameters describing species in the model, we focused on the three most relevant for our objectives – intrinsic growth rate, tolerance to shade and tolerance to drought – as described below. In addition, our approach differs from previously mentioned studies (Oddou-Muratorio and Davi 2014, Kramer et al. 2015, Fririon et al. 2024) in that we considered overlapping generations in order to better represent unmanaged forests. These studies also used quantitative trait loci (QTL) to simulate phenotypic variation and heritability, whereas we have opted for the infinitesimal model (Fisher 1919, Barton et al. 2017). To draw a parallel between these two genetic models, the QTL model assumes that traits are governed by a finite number of loci of quantitative effect and their inheritance is explicitly modeled when reproduction occurs. In the infinitesimal model, however, the phenotypes of the offspring are drawn in a Gaussian distribution around the average phenotype of two parents and this distribution is assumed to result from the contribution of a large number of genes with small and additive effects. Furthermore, contrary to QTL approaches, the infinitesimal model does not explicitly represent changes in allele frequencies in the population and includes simpler (and therefore faster) calculations, while retaining a certain genetic realism. This model is recommended by some authors for more ecologically focused models that aim to take account of species evolution (Romero-Mujalli et al. 2019).

As a proof of concept for this new evolutionary module, we focused on the evolution of a beech-fir forest (*Fagus sylvatica* and *Abies alba*) based on a specific site in south-eastern France, which we simulated with a climate change scenario corresponding to the representative concentration pathways (RCP) 4.5 (with an average warming of 3°C and a reduction in annual rainfall of 120 mm over the century). We chose the beech-fir forest because it is a common species assemblage in the European montane ecosystem (Toïgo et al. 2015). Moreover, the chosen location is in the southern part of the range of both species, where they are particularly vulnerable to global warming. Indeed, in this specific context, the ForCEEPS model predicts the extinction of beech and fir

in the absence of species evolution, making it an ideal setup for assessing the conditions for which evolutionary rescue could emerge, and the effect of management practices on these conditions. We had three main objectives: 1) identifying the conditions for evolutionary rescue of beech and fir, in monospecific (hereafter 'monospecific treatment') and mixed forests (hereafter beech-fir 'mixture treatment'); 2) assessing the impact of assisted gene flow on the prospects for evolutionary rescue and on evolutionary trajectories of beech and fir populations, by importing beech and fir individuals pre-adapted to more arid conditions (hereafter 'assisted gene flow treatment') (Hlásny et al. 2014, Fady et al. 2016); and 3) evaluating the consequences of introducing new species in the forest that are more tolerant of water stress (Table 1) on the evolutionary rescue conditions and evolutionary trajectories of beech and fir (hereafter 'assisted migration treatment') (Hlásny et al. 2014, Hisano et al. 2018). We considered the two treatments 'assisted gene flow' and 'assisted migration' to be human-driven movement of reproductive material for climate change adaptation. Therefore, we refer to these treatments as 'climate change adaptation measures' or 'management practices'. However, these treatments could also be interpreted as the natural migration of individuals from an external pre-adapted population into a population maladapted to future conditions, such as the natural migration of drought-tolerant species into the beech-fir forest of our focus site.

## Material and methods

We integrated our evolutionary module into the PBM ForCEEPS, a model from the gap model family developed to simulate both managed and freely evolving forests (Morin et al. 2020, 2021) (Fig. 1, blue, green and yellow boxes). Gap models simulate the growth, mortality and regeneration of trees in a collection of patches of land. ForCEEPS is individual based, and relies on calculating an optimal growth rate for each tree using empirical equations from Bugmann (1996). This optimal growth rate is subsequently reduced

<table><caption>Table 1. Initial species-specific parameter values for both local and imported species at the Drôme site (26 French department). Initial values of the parameters related to the three candidate traits for evolution, i.e. tolerance to drought (*DrTol*) or shade (*ShTol*) and its counterpart, the seedlings' tolerance to shade (*ShTolSI*), and intrinsic growth rate (*g*). Additionally, *AM*<sub>s</sub>, the age at maturity of the species is not subject to evolution but is used for the computation of evolutionary rates. For the parameter values of other species parameterized in ForCEEPS, see the Supporting information (for the other species-specific parameters values and description we refer to Morin et al. 2021).</caption><thead><tr><th>Species name</th><th>g<sub>s</sub></th><th>DrTol<sub>s</sub></th><th>ShTol<sub>s</sub></th><th>ShTolSI<sub>s</sub></th><th>AM<sub>s</sub></th></tr></thead><tbody><tr><td>Abies alba</td><td>350</td><td>0.23</td><td>1</td><td>0.05</td><td>40</td></tr><tr><td>Fagus sylvatica</td><td>260</td><td>0.25</td><td>1</td><td>0.05</td><td>45</td></tr><tr><td><strong>Imported species</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Quercus ilex</td><td>79</td><td>0.45</td><td>5</td><td>0.1</td><td>30</td></tr><tr><td>Q. petraea</td><td>246</td><td>0.33</td><td>7</td><td>0.2</td><td>30</td></tr><tr><td>Q. pubescens</td><td>146</td><td>0.33</td><td>7</td><td>0.3</td><td>30</td></tr><tr><td>Q. robur</td><td>249</td><td>0.3</td><td>9</td><td>0.3</td><td>40</td></tr><tr><td>Pinus cembra</td><td>115</td><td>0.3</td><td>5</td><td>0.2</td><td>30</td></tr><tr><td>P. montana</td><td>138</td><td>0.37</td><td>9</td><td>0.4</td><td>30</td></tr><tr><td>P. pinaster</td><td>350</td><td>0.4</td><td>7</td><td>0.25</td><td>24.5</td></tr><tr><td>P. sylvestris</td><td>150</td><td>0.37</td><td>9</td><td>0.3</td><td>20</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Figure 1. Overview of the ForCEEPS (Forest Community Ecology and Ecosystem Processes) model and recent developments discussed in this article. (for further details, refer to Morin et al. 2021). Teal box illustrates the feedback between the forest species composition and regeneration: more abundant species contribute more to the regeneration. Furthermore, like in most classic gap models, an option is available to activate a uniform seed rain for all species, irrespective of abundance, and is referred to here as 'migration'. Purple box: The evolution module allows the evolution of a set of traits (species parameters) under the assumptions of the infinitesimal quantitative genetic model. Each individual trait is comprised of a genetic ($G_{i,s}$) and an environmental ($E_i$) contribution. The module accommodates simultaneous evolution of covariant traits along established trade-offs. The initial state of any evolution simulation is a forest inventory generated from</figcaption></figure>
Figure 1. Continued

a 'baseline simulation' (see methods). (A) For a locally generated offspring $o$, its phenotype is determined by the average genotype of two randomly selected local parents ($G_{m,s}$ and $G_{f,j}$) to which we add within-family variation ($\gamma_o$) and environmental noise ($E_o$). An individual can be randomly drawn as a parent if it is older than the species age of maturity parameter (Table 1) and if it is growing above a certain threshold (not in a 'slow-growth' state). Larger reproducers are more likely to be parents. (B) The genetic value of initial individuals and migrant offspring are drawn from predefined distributions.

according to several environmental stressors, such as lack of light, water, heat or nitrogen. Tree mortality depends on two components: background mortality and growth-related mortality. The former depends on the species' maximum longevity and a set background death probability; the latter has a chance of occuring if a tree experiences three consecutive years with growth below a certain threshold (hereafter 'slow growth state'). Furthermore, we modified the regeneration process in ForCEEPS to account for interspecific competition. Our approach to modeling regeneration is more precisely described later in this section. In gap models, individuals are generally not spatialized within a patch to reduce simulation time. Trees in the same patch compete for light and share the same environmental conditions, and so the carrying capacity of a patch is an emergent property. The sum of several patches simulates a forest or a representative sample of a given forest within a forest matrix (Shugart et al. 2018). The main process in the model is competition for light within these patches, considered to be the main driver of temperate forest dynamics (Kohyama and Takada 2012). However, tree to tree competition for light is indirectly influenced by the effect of environmental conditions, especially climate, on growth (Morin et al. 2020). The ForCEEPS model operates on an annual time step, but takes into account monthly climatic variations (temperature and precipitation) to calculate the water balance of the patches. Species in ForCEEPS are defined by a specific set of parameters (Table 1). To simplify the interpretation of simulation results, we have limited our scope to the evolution of three of these traits, most related to our research questions. We focused on the adaptive dynamics of drought tolerance (hereafter *DrTol*) as it is strongly selected by environmental change, and therefore presumably also under selection from global warming. Negative trade-offs were imposed between *ShTol* and *DrTol*, and between *g* and *DrTol*, to account for the relationships observed between drought tolerance and competitive abilities in tree species (Niinemets and Valladares 2006, Kattge et al. 2020). The other species-specific parameters do not undergo evolution, and a single value per parameter is assigned to all individuals of the same species.

## Regeneration and species composition feedback module

In individual-based forest PBMs, two main approaches are used to simulate regeneration. The first approach explicitly links seedling establishment to the fecundity of mature trees: the number of new seedlings is determined by the number of seeds produced within the simulation, followed by the application of a mortality filter. This method is commonly used in models that incorporate evolutionary processes. The second approach, more typical of gap models such as ForCEEPS,

assumes a fixed seedling density per unit area, calibrated from inventory data, regardless of the actual fecundity of mature trees. Here too, a mortality filter is applied (Bugmann 1996, Didion et al. 2009, Morin et al. 2021). With this approach, regeneration therefore consists in a constant annual 'seed rain' including all species considered to be able to reproduce in the simulation, independently of their current abundance and density in the forest. This second approach strongly constrains evolution, first because it eliminates fecundity selection – the fact that some adults have higher fecundity than others due to differences in functional trait values. Moreover, it prevents interspecific competition from impacting regeneration. To overcome these limitations, in this study we separated regeneration into two distinct processes: 'local regeneration' and 'migration regeneration'. On the one hand, for local regeneration, and building on earlier work (Vallet, unpubl.) we multiplied the original target number of seedlings predicted by ForCEEPS for each species by its relative abundance (Fig. 1 teal box, arrow A). The relative abundance of species is assessed based on their basal area (BA) ($m^2 ha^{-1}$) across all patches, i.e. the sum of trunk areas at a height of 1.30 $m ha^{-1}$. Moreover, individual reproductive success explicitly depends on individual BA and thus on heritable trait values, which allows accounting for fecundity selection. On the other hand, migration regeneration simply allows the introduction of individuals of different species or with different trait values into the forest during the simulation (Fig. 1 teal box, arrow B). The migration regeneration is identical to the classic 'seed rain' gap model scheme and the original ForCEEPS regeneration, but weighted by a migration parameter ($M$). The parameter $M$ is a continuous index ranging from 0 to 1, where 0 indicates that only local regeneration occurs, and 1 indicates that the complete seed rain is added on top of the local regeneration, effectively doubling the target number of seedlings per unit area. This last point signifies that migration can only raise the seeding density. Therefore, $M$ must be maintained at a low value ($< 0.2$) to prevent the risk of the model operating outside its validity range. Regarding the seedling-related parameters such as $ShTolSl_s$ (Table 1), it is worth noting that in the ForCEEPS model, these parameters only play a role during the regeneration process on the mortality filter. Once a tree has passed this filter, it is established in the forest with a height of 1.3 $m$, and its adult parameters such as $ShTol_s$ start affecting its growth. Additionally, in our framework, the 'age at maturity' parameter ($AM_s$) (Table 1) is only used to determine the age at which a tree can participate in reproduction.

## Evolution module: evolvable traits

The core of this work was the incorporation of a module simulating evolutionary dynamics into the ForCEEPS model
(Fig. 1 purple box). This module required species traits to be both variable and heritable, with natural selection and drift impacting the distribution of these traits through time. The validations of drift and natural selection processes in simple contexts are available in the Supporting information. As stated above, for this study, we specifically selected three species-specific parameters as focal evolving traits, chosen for their relevance in responding to selection pressures imposed by climate change and competition:

* The individual's drought tolerance (DrTol<sub>i,t</sub>): a continuous index ranging from 0 (where drought stress prevents all growth and regeneration) to 1 (where drought stress has no impact on the individual's growth or regeneration).

* The individual's shade tolerance (ShTol<sub>i,t</sub>) and its counterpart, the seedlings' shade tolerance (ShTol<sub>sl,t</sub>). ShTol<sub>i,t</sub> is a continuous relative index ranging from 1 for the most shade tolerant species to 9 for the most shade intolerant species. ShTol<sub>i,t</sub> determines the impact light deprivation has on growth. Similarly, ShTol<sub>sl,t</sub> is a relative index ranging from 0 (tolerant) to 1 (intolerant) that determines the light availability required for seedling establishment.

* The individual's maximal growth rate ($g_{i,t}$), unitless. The fastest growing species in our species pool is Abies alba with an initial $g$ of 350 and the slowest growing species is Quercus ilex with an initial $g$ of 79.

The initial trait values before evolution: DrTol<sub>i,t</sub>, ShTol<sub>i,t</sub> and ShTol<sub>sl,t</sub> (Table 1) come from a compilation of previous studies (Ellenberg and Mueller-Dombois 1965, Bugmann 1996, Niinemets and Valladares 2006, Rameau et al. 2008) and calibrations using climatic and growth data from forest inventories (Morin et al. 2021). Initial values for $g$ are inferred from actual measurements from the French national forest inventory (Morin et al. 2021).

## Quantitative genetic model: initialization of traits values

The three traits described above are treated as quantitative and are assumed to be influenced by a large number of genes with small and additive effects, each of which contributes infinitesimally to an individual's phenotype (infinitesimal model sensu Fisher 1919) (Barton et al. 2017). Consequently, traits initially follow a normal distribution in populations. Each initial state of an evolution simulation must consist of a forest inventory, derived from field observations or previous simulations. From there, a phenotype is assigned to each individual in the form of an array of values for each evolving trait by randomly drawing from multi-normal distributions (Eq. 1):

$$ \bar{P}_{i,t} = \bar{G}_{i,t} + \bar{E}_i \qquad (1) $$

Where $\bar{G}_{i,t} \sim \mathcal{N}(\bar{\mu}_s, \Sigma_G)$ and, $\bar{E}_i \sim \mathcal{N}(\bar{0}, \Sigma_E)$

The phenotype array $P_{i,t}$ of individual $i$ of species $s$ is composed of a genetic and an environmental component, respectively

$G_{i,t}$ and $E_{i,t}$, $\mu_s$ is an array containing the initial values of the ForCEEPS species-specific parameters (DrTol<sub>i,t</sub>, ShTol<sub>i,t</sub> and $g$). Consequently, $G_{i,t}$ and $P_{i,t}$ are centered around the initial species-specific parameters described above, whereas $E_{i,t}$ is a zero-centered random non heritable noise. $\Sigma_G$ and $\Sigma_E$ are the genetic and environmental variance-covariance matrices as described below (Eq. 2):

$$ \Sigma_G = h^2 \times \Sigma_P \quad \text{and} \quad \Sigma_E = (1-h^2) \times \Sigma_P \qquad (2) $$

where $h^2$ is the heritability and $\Sigma_P$ is the phenotypic variance-covariance matrix. Equation 2 indicates that $\Sigma_G$ and $\Sigma_E$ are proportional to $\Sigma_P$, which is a major simplification of reality. This choice was necessary to reduce the number of input parameters for our evolutionary module, especially since it is extremely difficult to obtain estimates for the genetic and environmental covariances of traits. $\Sigma_P$ values depend on the initial intraspecific trait variability ratio (VR) which is the intraspecific SD of a given trait divided by its interspecific SD (SD<sub>intra</sub> and SD<sub>inter</sub>, respectively) (see Eq. 3 and the Supporting information for the species pool used to compute SD<sub>inter</sub>). As a result, increasing VR mechanically increases both the genetic and environmental trait variation (Eq. 3):

$$ \Sigma_P = \operatorname{diag}(\overline{\mathrm{SD}_{\mathrm{intra}}}) \times \Sigma_c \times \operatorname{diag}(\overline{\mathrm{SD}_{\mathrm{intra}}}) \\ \text{Where } \overline{\mathrm{SD}_{\mathrm{intra}}} = \mathrm{VR} \times \overline{\mathrm{SD}_{\mathrm{inter}}} \qquad (3) $$

where $\overline{\mathrm{SD}_{\mathrm{intra}}}$ is the intraspecific traits standard deviation array and $\Sigma_c$, the trait correlation matrix. Trait correlation values were obtained from interspecific correlations of proxy traits for DrTol, ShTol and $g$ found in literature and databases (Niinemets and Valladares 2006, Kattge et al. 2020). We used proxies to expand the scope so as to include the largest possible set of species. By applying these interspecific correlations to the intraspecific scale, we assumed that the trade-offs shaping trait relationships at the macroevolutionary scale also apply at the microevolutionary scale. We set those correlations to 0.36 between DrTol and ShTol and 0.38 between DrTol and $g$ (see the Supporting information for more details).

Seedling shade tolerance (ShTol<sub>sl</sub>) evolves only through changes in ShTol, making it an extension of this trait. The seedling shade tolerance for an individual, ShTol<sub>i,t</sub>, was linearly derived from the difference between the individual's shade tolerance for growth, ShTol<sub>i,t</sub>, and the species' initial shade tolerance value (ShTol). A 1/8 ratio was applied to convert the difference ShTol<sub>i,t</sub> - ShTol<sub>s</sub> on a scale of 1 to 9 into the difference ShTol<sub>sl,t</sub> - ShTol<sub>s</sub> on a scale of 0 to 1 (Eq. 4):

$$ \mathrm{ShTol}_{\mathrm{sl},t} = \mathrm{ShTol}_{\mathrm{sl},s} + \frac{1}{8}(\mathrm{ShTol}_{\mathrm{sl},t} - \mathrm{ShTol}_{\mathrm{sl},s}) \qquad (4) $$
## Quantitative genetic model: trait inheritance

After the initialization, traits are inherited through the process described below and illustrated in Fig. 1 (purple box, arrow A) during the simulation runtime. For locally regenerated individuals, new seedlings of species s are attributed to two parents randomly chosen among a pool of reproducer trees. A tree is considered as reproducer if it meets two criteria: 1) the tree is older than its species-specific maturity age (Table 1) and 2) it is not in state of ‘slow growth’. An individual tree is considered in a slow growth state if it has been growing less than a specific threshold for three consecutive years. This phenomenon occurs when environmental factors such as temperature or light levels are too restrictive for the tree and simulates the process of carbon starvation (Morin et al. 2021). We assume that, in such conditions, trees do not invest in reproduction. The probability of each tree being drawn as a parent is weighted by its relative size among the reproducers (assessed by trunk diameter), making the assumption that larger trees produce more offsprings. Once a seedling has been attributed two parents, each of its non-evolving traits is set to the species’ base trait value (Morin et al. 2021) while the phenotypic values of each of the three evolving traits is defined as the mean of the genetic trait values of both parents, to which a genetic ($\gamma_o$) and an environmental noise are added ($e_o$) (Eq. 5):

$$ \bar{P}_{o,s} = \frac{\overline{G_{m,s}} + \overline{G_{f,s}}}{\underbrace{2}_{\overline{G_{o,s}}}} + \bar{\gamma}_o + \bar{E}_o \quad (5) $$

where $\bar{\gamma}_o \sim N(\bar{0}, \Sigma_G \times \frac{1}{2})$ and $\bar{E}_o \sim N(\bar{0}, \Sigma_E)$

with $P_{o,s}$ the phenotype of offspring o, for a given species s and $G_{m,s}$ and $G_{f,s}$ the respective genotypes of the mother and father. The genetic noise $\gamma_o$ originates from the within-family genetic variance, also called the segregation variance as it results from the process of Mendelian segregation of alleles. Under the infinitesimal model, assuming sufficiently large and stable ancestral populations and ignoring inbreeding, this genetic noise variance is constant across generation and normally distributed, regardless of the distributions of the parents’ genetic values (Bulmer 1971). Without selection, setting the within-family genetic variance to half the initial genetic variance keeps the overall population’s genetic variance constant from one generation to the next, as computing the mean breeding value of two parents to obtain that of the offspring will result in reducing population’s genetic variance by half each generation. In contrast, for seedlings originating from migration (Fig. 1 purple box arrow B), phenotypes are drawn in multi-normal distributions as described in Eq. 1, with the difference that the means can differ from those of local seedlings. These distributions can either be set at the beginning of the simulation, or can be updated during the on-going simulation to best represent the genetic

composition of replicate populations evolving simultaneously. In the present case, migrants’ trait means and variances are initially set depending on the considered treatment and evolutionary scenario. Finally, it is worth noting that while heritability $h^2$ and VR are prescribed as input parameters to initialize trait values and introduce genetic noise in offspring (through $\sum_G$, $\sum_F$ and $\sum_p$), the realized genetic and phenotypic variances (and accordingly, heritability) change dynamically throughout the simulation as selection and genetic drift act on the reproducing individuals.

## Site and climate data

Our study is based on two sites: a focal site in south-east France (hereafter ‘Drôme’), where beech and fir forests are common ecosystems (44°55'01.2"N, 05°17'27.6" E: 1146 mm of yearly precipitation and mean temperatures of 5.9°C), and another site with drier conditions (hereafter ‘Montagne Noire’), at the southern limit where beech-fir forests can be found (43°27'55.4"N, 02°16'43.2"E: 983 mm of yearly precipitation and mean temperatures of 9.3°C). The model’s site parameters (i.e. other than climate conditions) for the Drôme site were derived from field measurements, and those for the Montagne Noire site were arbitrarily set to be identical to ensure comparability between the two sites (Supporting information). Climate data consisted of monthly temperature and precipitation data available from the CHELSA dataset (CMIP5, BCG model) (Karger et al. 2020), which were split into three periods:

* Historical period (years 1000–2000 in the simulation): climate data from 1969 to 2000 (from the historical CHELSA dataset) were repeatedly used to represent 1000 years of past climate conditions.

* Warming period (years 2000–2100 in the simulation): climate data from 2000 to 2100, based on the RCP 4.5 scenario, were applied to simulate a century of climate change.

* Stabilization period (years 2100–3000 in the simulation): climate conditions from 2090 to 2100 were cyclically repeated for 1900 years to represent a stabilized climate.

Simulations begin in the year 1000, starting with 1000 years of historical climates, followed by climate change between the years 2000 and 2100, and concluded by a stabilization period from 2100 to 3000, i.e. a total of 2000 years of evolution.

## Simulation plan

### Does mixing fir and beech dampen evolutionary rescue?

The evolutionary rescue conditions of each species and the evolutionary trajectories of traits during and after the three climatic periods were compared between the two treatments: beech-fir mixtures and beech-fir monospecific forests. We ran 50 replicates of each of the three compositions for each evolutionary scenario at the Drôme site, with associated climate projections (a total of 50 × 3 × 16 = 2400 simulations). Methods for quantitatively describing evolutionary rescue
and evolutionary trajectories are described in the ‘Analyses of model output and implementation’ section.

**Does assisted gene flow facilitate evolutionary rescue?**

Second, we compared evolutionary rescue condition and evolutionary trajectories of traits between the beech-fir mixture treatment and an 'assisted gene flow' treatment. This latter was identical to the beech-fir mixture treatment, with the difference that we added immigration of pre-adapted beech and fir individuals during the 100-year warming period. The migration parameter *M* was arbitrarily set to 0.1 to allow new genotypes to colonize with a marginal demographic impact (based on trial simulations). The phenotypic distributions of migrant traits (Supporting information) are centered around the average phenotypes across 25 replicates of beech-fir mixed forests that have evolved over 1000 years in the historical climates of the Montagne Noire (with values of *h*² and VR corresponding to the evolutionary scenario considered). We carried out only 25 replicates of this treatment, as computation time was significantly higher (25 × 16=400 simulations).

**Does assisted migration prevent local species' evolutionary rescue?**

Last, the evolutionary rescue conditions and evolutionary trajectories of traits were confronted between the beech-fir mixtures treatment and an 'assisted migration' treatment. This treatment was similar to the assisted gene flow treatment, but instead of adding beech and fir individuals during the warming period, individuals from a set of eight other species were introduced through migration (*M*=0.1). The criteria for choosing imported species are 1) that they are parameterized in ForCEEPS; and 2) that they have a higher drought tolerance than the two local species (Table 1). The phenotypic values of imported species were drawn from multivariate normal distributions, with the means for the evolving traits corresponding to the initial values parameterized for these species (Table 1). Imported species were also subject to evolution as soon as mature individuals established themselves in the simulated forests. For reasons of computational efficiency, we also carried out only 25 replicates of this treatment, starting exclusively from mixed beech-fir forests (25 × 16=400 simulations).

**Simulating initial conditions**

To obtain mature, age-structured forests in a pseudo-equilibrium state to serve as the initial state for all our later simulations, we carried ‘baseline simulations’ without evolution (years 0–1000 in the simulation). These simulations consist of 1000 years of ForCEEPS simulations from bare soil under historical climates without the evolution module. For the beech-fir mixture treatment, the baseline simulations included 50 patches of 1000 m² each of mixed beech-fir forest, with beech being the most abundant species. These simulations produced inventories containing on average 720 beech trees and 90 fir trees, representing a respective BA of 11.60 m²

ha⁻¹ and 1.91 m² ha⁻¹ (Supporting information). To isolate the effect of interspecific interactions between species from the effect of the maximum population size, we constrained the number of individuals to the same order of magnitude between beech-fir mixtures and monospecific baseline simulations. In the Drôme site and climate, the ForCEEPS model predicts that a mixed beech-fir forest will be strongly dominated by beech; therefore, the final size of the fir population will vary greatly between monospecific or mixed forest simulations, complicating the comparison of evolutionary trajectories, as the effects of interspecific interactions become indistinguishable from those caused by genetic drift. To remedy this, we reduced the number of patches in monospecific baseline simulations by running 42 patches of monospecific beech and 11 patches of monospecific fir. This yielded an average of 693 beech trees with a BA density of 14.13 m² ha⁻¹, and 76 fir trees with a density of 8.10 m² ha⁻¹. Then, to address our three research questions, we implemented the simulation framework depicted in the Supporting information, with the evolution module activated on the basis of the inventories resulting from the baseline simulations. As *DrTol* and *ShTol* cannot be measured directly on individuals in field studies, it is impossible to rely on empirical estimates of their heritability or variability in real populations. Consequently, to identify the conditions for evolutionary rescue, each question was explored through four scenarios of *h*² (0.1, 0.2, 0.3, 0.4) and four scenarios of VR (0.1, 0.2, 0.3, 0.4), amounting to 16 evolutionary scenarios. These scenarios cover cases ranging from the systematic extinction of both species in all replicates to a recovery above or near the initial BAs for both species after 1000 years of climate stabilization. These evolutionary scenarios also cover realistic ranges of *h*² (Alberto et al. 2013) and VR (Rosas et al. 2019, Fririon et al. 2024) for traits related to growth and hydraulic stress response.

**Analyses of model output and implementation**

For further analysis, the demographic state of a forest was assessed using the BA (m² ha⁻¹) of each species. BA is commonly used in silviculture and forest ecology as a proxy of forests health status. We consider that the BA can be used as a proxy for population size, since in a population where all individuals have similar age and dimensions, a reduction in BA directly reflects a decline in population size, and this relationship can be extended to populations with mixed age structures at equilibrium (Supporting information). Incidentally, large mortality events can also lead to a slight rejuvenation of the populations (Supporting information). To facilitate comparison between the treatments and species studied, BA can be expressed in relative terms (RBA) (%). In this case, the BA in a given year is expressed as a percentage relative to a reference historical mean BA between the years 1500 and 2000. We are interested either in the RBA of the year in which the BA is lowest (minimal RBA) or in the RBA of the last year of simulation (final RBA). To represent in a single figure the results of all our simulations, we compute the median values of these quantities for all replicates in (Eq. 6):
$$
\text{Median minimal RBA} = \text{Median} \left( \frac{\text{Min}(\text{BA}_{1000:3000,r})}{\text{Mean}(\text{BA}_{1500:2000,r})} \right) \times 100
$$
$$
\text{Median final RBA} = \text{Median} \left( \frac{\text{BA}_{3000,r}}{\text{Mean}(\text{BA}_{1500:2000,r})} \right) \times 100 \quad (6)
$$
$$
\text{for } r = 1, 2, \infty \dots N_{\text{rep}}
$$

Next, we consider a population to have undergone evolution-
ary rescue if two criteria are met. First, its BA and therefore its
population size never goes down to 0 (persistence criterion);
and second, its growth rate between the date when the BA is
minimal and the last date of the simulation is significantly
different from 0 with a confidence threshold of 5% (recovery
criterion). The growth rate is assessed by the slope of a linear
regression between the BA and the time during the recovery
period.

To describe trait evolution, trait values per species *s* were
averaged by weighting them by the BA of each individual.
This approach ensures the values are representative of the
population of a given replicate *r* in a given year *y* and prevents
over-representation of a large number of unfit seedlings that
will not survive several consecutive years. We refer to those
weighted arithmetic averages as *Trait*<sub>s,r,y</sub> for the three respec-
tive traits *DrTol*, *ShTol* and *g*. For each evolutionary scenario,
average trait evolutionary rates across replicates (*Nrep*) were
also measured in Haldanes (*H*<sub>0</sub>) between the years 2000 (*y*1)
and 2100 (*y*2), depending on species *s* maturation age (*AM*<sub>s</sub>)
(see Eq. 7 and *Gingerich* 1983). Evolutionary rates values
are provided, along with their 95% confidence intervals
across replicates for the *h*<sup>2</sup> and VR=0.3 evolutionary sce-
narios, as this setting led to the most realistic evolutionary
rates according to field studies and experiments (Bone and
Farres 2001). We also computed weighted average *DrTol* dif-
ferences between the end of the simulations (*DrTol*<sub>s,r,3000</sub>)
and the mean over a reference historical period (*DrTol*<sub>s,r,ref</sub>)
i.e. Δ*DrTol*<sub>s,r</sub>. The chosen reference period spans from 1500
to 2000. We did not include years from 1000 to 1500 in
the reference period to allow time for ‘local adaptation’.
Evolutionary rates in *H*<sub>0</sub> enable comparison with other stud-
ies in evolutionary ecology, while Δ*DrTol*<sub>s,r</sub> reflects the mag-
nitude of evolutionary change over the entire duration of the
simulations in our particular framework, enabling us to com-
pare the scenarios with each other (Eq. 7):

$$
\overline{H_{0,\text{Trait},s,y1,y2}} = \frac{1}{N_{\text{rep}}} \sum_{r=1}^{N_{\text{rep}}} \left( \frac{\overline{\text{Trait}_{s,r,y2}} - \overline{\text{Trait}_{s,r,y1}}}{\text{SD}(\overline{\text{Trait}_{r,y1,y2}})} \right) \times \frac{\overline{AM_s}}{y2 - y1} \quad (7)
$$

Climate, site and trait values were analyzed and trans-
formed into the desired formats and figures using R ver.e
4.3.0 (www.r-project.org). The ForCEEPS model, as well
as its evolution module, is available via Capsis, a simulation

platform for forest growth and dynamics models (Dufour-Kowalski et al. 2012).

**Results and discussion**

Does mixing fir and beech dampen evolutionary rescue?

The median minimal relative basal area (RBA) (Fig. 2) can be
interpreted as the population’s safety margin for persistence.
VR and h² had a strong impact on minimal RBA (Fig. 2B)
and population size (Supporting information), highlighting
the relevance of the evolutionary processes even in the short-
term, as minimal RBA typically occurs soon after the end of
the warming phase (year 2150–2250) (Fig. 2C, 3, Supporting
information). We use the terms ‘slower evolutionary scenar-
ios’ for cases with lower h² and VR values, and ‘faster evo-
lutionary scenarios’ for those with higher h² and VR values,
for which we expect respectively lower and higher evolution
rates. For fir, evolutionary rescue conditions were met for
minimum values of h² and VR of 0.2 and 0.3, respectively, in
both monospecific and beech–fir mixture treatment popula-
tions. In the monospecific beech treatment, however, evolu-
tionary rescue conditions were met for ‘slower’ evolutionary
scenarios of h² of 0.1 and VR of 0.3 (Fig. 2B). The most
pre-adapted species encountered the lowest risks of extinc-
tion, i.e. beech with ΔDrTol<sub>s</sub>=0.25 as opposed to fir with
DrTol<sub>s</sub>=0.23. The lack of differences between the mono-
specific and mixture treatments in terms of safety margin
was due to a relaxation of competition for light constraints,
resulting from the high mortality in the warming period
(Fig. 3, Supporting information). Therefore mixing did not
impact the chances of persistence for both species. We used
the final RBA to assess the population’s ability to recover its
pre-warming BA by the end of the 1000-year stabilization
period in a warmed climate. In contrast to the short-term
response, the final RBA of fir populations showed differences
between the monospecific and mixture treatments. Fir recov-
ery was generally lower in the beech–fir mixture treatment
than in the monospecific treatment for the fastest evolution-
ary scenario (h² and VR=0.4, Fig. 2D), with a median final
RBA of 77% and 118%, respectively.

An analysis of trait evolutionary trajectories is necessary
to understand why evolution tends to favor beech over fir
populations. First, *DrTol* evolution towards higher values
was critical for population recovery, with populations that
reached higher Δ*DrTol* showing the strongest recoveries. As
expected, higher Δ*DrTol* values were observed for higher *h*²
and VR scenarios (Fig. 4). However, interspecific interactions
between beech and fir had no impact on the evolutionary
trajectories of the three traits in the short term, as there were
no differences in fir and beech evolutionary rates between the
monospecific and the beech–fir mixture treatment (Fig. 5A),
regardless of the evolutionary scenario (Supporting informa-
tion). Fir *DrTol* evolutionary rates were 0.23 *H*<sub>0</sub> [0.22; 0.25]
in the monospecific treatment compared with 0.23 *H*<sub>0</sub> [0.20;
0.25] in the mixture. For beech, the rates were 0.35 *H*<sub>0</sub> [0.34;
0.36] for monospecific and beech–fir mixture treatments. The
<figure><img src="image_6.png" /><figcaption>(B) Median minimal relative basal area (RBA) (%)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(C) Median basal area change over time
$h^2 = 0.3 \ \& \ VR = 0.3$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(D) Median final relative basal area (RBA) (%)</figcaption></figure>

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

Figure 2. (A) Illustration for an example simulation of the reference, minimum, and final basal area (BA) used to calculate minimal and final relative basal area (RBA) (see methods). (B) Median across replicates of the minimal RBA in a forest of Drôme site for the various scenarios of heritability ($h^2$) and intra- to interspecific trait standard deviation ratio (VR). (C). Median across replicates of the BA change over time, for the $h^2$ and VR = 0.3 evolutionary scenario in beech-fir mixture and in the assisted migration treatment. (D) Median across replicates of the final RBA, in a forest of Drôme site for the various scenarios of heritability ($h^2$) and intra- to interspecific trait standard deviation ratio (VR). For (B) and (D), the number of replicates is shown under the treatment type. In addition, the numerical values of the minimum and final RBA are shown in white for the evolutionary scenario in which the RBA is greatest, as well as for the most parsimonious evolutionary scenarios for which more than half of the simulations experienced an evolutionary rescue. The parsimony criterion for an evolutionary scenario *k* is $\text{Min}(\text{Max}(h_k^2, VR_k))$. Parsimony is evaluated by minimizing this criterion; in the case of a tie, the scenario where $h^2$ or VR is the lowest is retained.

same observations can be made for the two other traits (*ShTol* and *g*) (Fig. 5A). This result can be explained by the alle- viation of competition following the massive mortality epi- sodes that accompanied climate change. More surprisingly,

long-term evolutionary trajectories were also unaffected by competitive interactions between beech and fir individuals. This is evidenced by the overlap of evolutionary trajectories and destinations of monospecific and the beech-fir mixture
<figure><img src="image_10.png" /><figcaption>Figure 3. Median across 25 to 50 replicates of the mean relative light availability perceived by individuals of beech and fir (model output) in the four treatments for the evolutionary scenario with *h*² = 0.3 and VR = 0.3. The envelopes show the quantiles 0.025 and 0.975 over all repetitions of light availability. Bottom stripe color indicates different climatic periods. Results for other evolutionary scenarios are presented in the Supporting information.</figcaption></figure>

treatments for both species (Fig. 5B). Populations of both species showed an increase in *DrTol* throughout the simulations, accompanied by a decrease in *ShTol* and *g* over time due to the negative trade-off between *DrTol* and the two competitive traits (Fig. 5B, Supporting information). Lower recoveries were observed in fir populations for the beech-fir mixture treatment compared to the monospecific treatment for equivalent Δ*DrTol* evolution (see the curvature of the final RBA to Δ*DrTol* relationship in Fig. 4). We interpret this result as the return of competition for light as the main process shaping the community composition in the beech-fir mixture treatment, preventing fir recovery despite its tolerance to drought. In beech-fir mixtures, evolution favored the dominant species, which is in line with our expectations and theoretical works (De Mazancourt et al. 2008, Thompson and Fronhofer, 2019, Eldijk et al. 2020). We propose two non-exclusive interpretations of this phenomenon: 1) beech was favored because its populations were initially larger, reducing drift and increasing the probability of the population to host more drought-tolerant individuals, explaining the faster evolutionary rates in beech populations overall (Fig. 5A); and/or 2) beech was initially fitter due to its higher *DrTol*, allowing its populations to reach higher final *DrTol* values (Fig. 5B).

**Does assisted gene flow facilitate evolutionary rescue?**

As expected, the assisted gene flow treatment led to the highest *DrTol* evolutionary rates during the warming period, with beech $H_0$=0.50 [0.46; 0.51] and fir $H_0$=0.39 [0.35; 0.43] for the *h*² and VR=0.3 scenario (Fig. 5A). These values are within the largest known $H_0$ for plant adaptation to abrupt changes in natural populations (Bone and Farres 2001). Therefore, assisted gene flow led to a higher safety margin (minimal RBA) for both species and to a larger number of

evolutionary scenarios in which both beech and fir experienced evolutionary rescue (Fig. 2B and D). While assisted gene flow helped both species adapt to climate change in the short term (minimal RBA), it led to more complex relationships between evolutionary scenarios (*h*² and VR) and recovery (final RBA) in the long term. In contrast to the beech-fir mixture treatment, highest final RBA no longer occurred for the fastest evolutionary scenarios (*h*² and VR=0.4) for fir and beech populations individually and at the whole-forest scale (Fig. 2D). In the assisted gene flow treatment, higher evolutionary rates led to even lower values of *ShTol* in beech (Fig. 5B). These reduced *ShTol* values resulted in slower effective growth rates, as shading intensified during forest recovery (Fig. 3). Moreover for fir, lower recoveries occurred for equivalent Δ*DrTol* in the assisted gene flow treatment compared to the simple mixture (flatter Δ*DrTol* to final RBA recovery curve) (Fig. 4). This result can be explained by the fact that assisted gene flow may have asymmetrically favored the dominant species of the mixture (beech). Altogether, this treatment suggests that beyond a certain point, stronger evolutionary responses regarding the Δ*DrTol* no longer guarantee better recoveries. In this context, the trade-off between *ShTol* and Δ*DrTol*, along with interspecific competition, played a more decisive role in shaping the eco-evolutionary dynamics of forests.

**Does assisted migration prevent local species' evolutionary rescue?**

In the assisted migration treatment, fir and beech became extinct in the majority of evolutionary scenarios, including those where both species experienced evolutionary rescue in the beech-fir treatment (Fig. 2B and D). This result is consistent with our expectations and the literature on the subject, as the arrival of newly pre-adapted species can endanger local evolving species (De Mazancourt et al. 2008). We
<figure><img src="image_11.png" /><figcaption>Figure 4. Final relative basal area (RBA) (see Fig. 2A) as a function of evolutionary change in tolerance to drought i.e. ΔDrTol for fir and beech species, and for the whole forest (all species confounded). ΔDrTol is computed as the difference between the arithmetic weighted mean trait value in the year 3000 and the weighted mean value of the trait averaged over a reference period (1500–2000 years). The weights correspond to the basal area (BA) of the individuals in relation to the BA of their species or that of the forest as a whole (all species combined). Each point represents a population of a given replicate, species (or community) and treatment.</figcaption></figure>

observed a turnover between beech and sessile oak *Quercus petraea*, a more thermophilous species (with higher *DrTol*, see Table 1). Contrary to the previous treatments, in all evolutionary scenarios, fir populations decreased relatively less than beech populations (Fig. 2B and D). It is important to note that changes in BA are expressed in relative terms, so the fact that beech is no longer the dominant species largely explains this result. Second, if both beech and fir populations are subjected to a similar shading regime (Fig. 3), beech is expected to grow more slowly due to its lower g. (Table 1, Fig. 5B). Furthermore, fir and beech populations recoveries (final RBA) were lower in the assisted migration treatment for equivalent ΔDrTol in comparison with the two-species mixture treatment. The difference between the two treatments was greater for beech populations for which +0.07 in ΔDrTol allowed values to reach over 100% of final RBA

in the beech-fir mixtures against 25% or less in the assisted migration treatment (Fig. 4). *DrTol* evolutionary rates were unaffected by species introduction in the short term (during the warming period) with *H*₀=0.26 [0.23; 0.29] in fir and *H*₀=0.37 [0.36; 0.38] in beech for the assisted migration treatment against 0.23 [0.20; 0.25] in fir and 0.35 [0.34; 0.36] in beech for the beech-fir mixture treatment (Fig. 5A). Assisted migration tended to slow down *DrTol* evolution only in beech populations over the long term (Fig. 5B, Supporting information). Slower evolution in beech *DrTol* was associated with retention of a higher *ShTol* compared to the beech-fir mixture (Fig. 5B, Supporting information). This result is in line with the view of a niche conservatism induced by interspecific interactions (Peterson et al. 1999, Wiens et al. 2010). It is worth noting that forests with assisted migration
<figure><img src="image_12.png" /><figcaption>(B) Long-term evolutionary trajectories ($h^2$ = 0.3 & VR = 0.3)</figcaption></figure>

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

Figure 5. (A) Evolutionary rates in Haldanes for the three evolving traits before and after the warming period, across the four treatments for beech and fir. The grey dashed line marks the 0 evolutionary rate threshold. Only the $h^2$ and VR = 0.3 evolutionary scenario is presented here; other scenarios are available in the supporting information. (B) Median across replicates of the weighted average trait values in two dimensions over time. Line transparency reflects time progression, also highlighted by date labels. Framed numbers represent the number of replicates in which populations survived until the last simulation year (3000), with a maximum of 50 for the monospecific and beech-fir mixture treatments, and 25 for the assisted gene flow and assisted migration treatments. Vertical and horizontal error bars show the 0.025 and 0.975 quantiles around the median weighted average trait values across replicates at the simulation's end for each treatment. For individual population evolutionary trajectories, please see the Supporting information.
had the highest leaf cover and, consequently, were the forests where competition for light was the most important (Fig. 3, Supporting information). Altogether, these results indicate that assisted migration induces a faster return to a light-competition-driven regime in forests as new species colonize. At the whole-forest scale (sum of all species BA), assisted migration allowed for the most significant recoveries and the highest minimal RBA, with all final recoveries over 130% and all minimal RBA above 75% (Fig. 2B and D). Therefore, our case study tends to confirm the view that species range shifts could mitigate the impact of climate change on forest ecosystems in some scenarios, but at the risk of losing local species, provided migration rates are sufficiently high (Hisano et al. 2018).

## Model and simulation plan limitations

The initial-state forest densities (resulting from baseline simulations in historical climates) at the Drôme site (Supporting information) can be considered low compared with the forests that actually grow in this region (IGN|Institut national de l'information géographique et forestière 2024). This is because our study sites are located in the southern part of the range of the two focal species and that our model was calibrated for a single value of *DrTol*, per species for their range throughout Europe. The *DrTol* of the Drôme populations has therefore most certainly been underestimated, resulting in smaller populations and consequently an overestimation of the drift in our simulations. It is also important to mention that, during the course of the simulation, population sizes can fall to extremely low levels, especially for the most parsimonious evolutionary rescue scenarios (Supporting information). If we had taken into account large-scale disturbances such as fires and winds, it is highly likely that these populations, constrained to small sizes for long periods of time, would simply have become extinct. To clarify, we do not, in practice, predict evolutionary rescues for the most parsimonious scenarios of $b^2$ and VR, but rather this approach allows us to comment on the lower part of the graphs and the minimum conditions for observing evolutionary rescues in our framework between our different treatments.

On a different topic, the arbitrary choice of stabilizing the climate after 100 years of warming (at the end of the predictions available in climate projections) may limit our conclusions as this scenario is both arbitrary and likely unrealistic. Different post-warming trajectories (e.g. continued warming or cooling) could have led to different outcomes, highlighting the importance of halting climate change to allow species to adapt. In particular, the long-term effects observed in our study are influenced by the assumption of climate stabilization after 2100, as we noticed that in all treatments but the assisted migration one, populations started to recover only after the climate had stabilized (Fig. 2C, 3). This result is consistent with evolutionary theory. When a population undergoes a change of environment, we expect the evolution of its traits to always lag behind a moving optimal value (Bürger and Lynch 1995). This principle remained applicable even in the fastest evolutionary scenarios we explored. This shows that to simulate long-term impacts of climate change

on species evolution, not only do the speed and the absolute values of changing climate matter, but also the duration of the simulated change.

Another consequence of a rapid climate change is a significant decrease in population sizes, thus increasing the likelihood of kinship between individuals. However, inbreeding has the effect of eroding the genetic variance between siblings or the within-family variance or $\gamma$ noise in our case Barton et al. (2017). As a result, the overall genetic variance, and therefore the evolutionary speed, may have been overestimated as a consequence of considering the within-family variance to be constant. Finally, only two selective forces operate in our simulations: 1) climate change, which will positively select *DrTol*; and 2) competition between individuals for light, which will select for higher *ShTol* and *g*. Given that we have imposed negative trade-offs between *DrTol* and the traits related to competition ability (*ShTol* and *g*), intraspecific interactions could only hinder the evolution of *DrTol*, consequently preventing evolutionary rescue. Introducing pre-adapted species to future environments in this context leads to a competitive hierarchy for light, ultimately resulting in the dominance of imported species over local ones. This has imposed two strong and purely contradictory selective pressures on local species, creating particularly stringent conditions for evolutionary rescue.

## Perspectives

The erosion of genetic variance discussed earlier could be accounted for in the future by tracking the pedigree of individuals (Barton et al. 2017). That being said, this effect is counterbalanced by gene flow, which generally keeps inbreeding at low levels in natural forest ecosystems (Petit and Hampe 2006) except in the case of isolated populations (Mimura and Aitken 2007). Tracking pedigrees between reproducers and reducing the within-family variance according to the kinship of the two parents in the present context would, however, result in assuming that the simulated forest is a completely closed system. The best solution we can think of would be to simulate realistic migration processes and assume that immigrants share no ancestry with the local individuals. Genetic variation would then be reduced by inbreeding and replenished by migration, with the balance between these processes depending on effective population size. However, tracking pedigree would also be computationally expensive for large populations. Implementing this solution could, however, be relevant for simulations of small populations where higher inbreeding is expected and tracking pedigree is less costly. Another reason to improve the integration of migration into models is to explore its demographic impact on populations. In our case, we have deliberately limited these impacts by keeping the migration parameter low ($M=0.1$) to highlight only the genetic consequences of the arrival of pre-adapted individuals. Nevertheless, assisted gene flow is also advocated for the demographic support they offer (Aitken and Whitlock 2013), and therefore, higher migration rates should ideally be tested. Studying the consequences of different migration intensities on evolutionary rescue remains a crucial step towards a better understanding of the consequences of
assisted gene flow. Beyond improving a model's realism and predictive capabilities, our priority is to better understand how species' adaptation to a new environment is impacted by ecological interactions. We therefore think it is urgent to explore cases where intra- and interspecific interactions are at the origin of selection pressures that accompany environmental change, and not only cases where they act in opposite directions. With moderate competitive selection pressure, a positive effect of competition on evolutionary rescue could be observed (Osmond and de Mazancourt 2013). The conditions for this to occur in the context of global warming could be explored within our framework – for instance if competition for water between individuals were explicitly considered in ForCEEPS. In such a case, increased competition between individuals as well as climate change may select for more drought-tolerant individuals.

## Conclusions

To our knowledge, this work constitutes the first attempt to simulate the evolution of multiple interacting species using a functional and mechanistic framework (PBM). Most importantly, we showed that both trait variability and heritability had a major impact on the model's predictions, even over a time scale of one or two centuries. Our results support the view that increasing genetic variability within stands is a key lever to enhance forests' capacity to cope with climate change. More broadly, they highlight the necessity of accounting for the evolutionary responses of tree species when projecting the impacts of climate change on forests. In addition, 'assisted gene flow' enabled species to keep pace with climate change without any apparent drawbacks emerging in our framework. We also showed that abiotic factors, such as climate, are not the only dynamic forces influencing species' evolution. Population dynamics also shape species' evolutionary trajectories through competitive interactions. However, interspecific competition did not lower the chances of evolutionary rescue in the short term (century) in the beech-fir mixture as the two focal species were primarily affected by climate change. The evolutionary trajectories of these two species were only affected by competition in the long term (millennium), particularly in the assisted migration treatment when more arid-tolerant species colonized the forests. In this treatment, more intense competition reduced the likelihood of evolutionary rescue for beech and fir. Nevertheless, assisted migration allowed for the fastest forest recoveries. Although it is a proof of concept, our approach represents an important step towards understanding the interactions between competition and adaptive constraints, as well as obtaining more robust predictions of the impacts of climate change on forests while accounting for evolutionary dynamics.

**Acknowledgements** – We are particularly grateful to the various developers of the Capsis platform and to Nogoye Diouf, who worked on the trait data from the TRY database.

## Author contributions

**Louis Devresse:** Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (equal); Validation (lead); Visualization (lead); Writing – original draft (lead); **Freya Way:** Conceptualization (supporting); Methodology (supporting); Software (supporting); Writing – review and editing (supporting). **Tanguy Postic:** Conceptualization (supporting); Writing – review and editing (supporting); **François de Coligny:** Conceptualization (equal); Methodology (equal); Software (equal). **Xavier Morin:** Conceptualization (supporting); Formal analysis (supporting); Methodology (supporting); Supervision (lead); Writing – review and editing (supporting).

## Data availability statement

The version of ForCEEPS used to produce the data presented here, and the R scripts for data analysis and figure production, as well as the Supporting information, are available online (https://zenodo.org/records/15172843; Devresse et al. 2024).

## Supporting information

The Supporting information associated with this article is available with the online version.

## References

Aitken, S. N. and Whitlock, M. C. 2013. Assisted gene flow to facilitate local adaptation to climate change. – Annu. Rev. Ecol. Evol. Syst. 44: 367–388.

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. and Savolainen, O. 2013. Potential for evolutionary responses to climate change – evidence from tree populations. – Global Change Biol. 19: 1645–1661.

Barton, N. H., Etheridge, A. M. and Véber, A. 2017. The infinitesimal model: definition, derivation and implications. – Theor. Popul. Biol. 118: 50–73.

Bell, G. 2017. Evolutionary rescue. – Annu. Rev. Ecol. Evol. Syst. 48: 605–627.

Bone, E. and Farres, A. 2001. Trends and rates of microevolution in plants. – *Genetica* 112–113: 165–182.

Bugmann, H. K. M. 1996. A simplified forest model to study species composition along climate gradients. – *Ecology* 77: 2055–2074.

Bulmer, M. G. 1971. The effect of selection on genetic variability. – Am. Nat. 105: 201–211.

Bürger, R. and Lynch, M. 1995. Evolution and extinction in a changing environment: a quantitative-genetic analysis. – *Evolution* 49: 151–163.

Cuddington, K., Fortin, M. J., Gerber, L. R., Hastings, A., Liebhold, A., O'Connor, M. and Ray, C. 2013. Process-based models are required to manage ecological systems in a changing world. – *Ecosphere* 4: art20.

De Mazancourt, C., Johnson, E. and Barraclough, T. G. 2008. Biodiversity inhibits species' evolutionary responses to changing environments. – *Ecol. Lett.* 11: 380–388.
Depauw, L. et al. 2024. Enhancing tree performance through species mixing: review of a quarter-century of TreeDivNet experiments reveals research gaps and practical insights. – Curr. For. Rep. 10: 1–20.

Devresse, L., Coligny, F. and Morin, X. 2024. Evolutionary rescue in a mixed beech-fir forest: insights from a quantitative-genetics approach in a process-based model – model, data, scripts and supporting Information. – https://doi.org/10.5281/zenodo.14197821.

Didion, M., Kupferschmid, A. D. and Bugmann, H. 2009. Long-term effects of ungulate browsing on forest composition and structure. – For. Ecol. Manage. 258: S44–S55.

Dufour-Kowalski, S., Courbaud, B., Dreyfus, P., Meredieu, C. and de Coligny, F. 2012. Capsis: an open software framework and community for forest growth modelling. – Ann. For. Sci. 69: 221–233.

Dyderski, M. K., Paź, S., Frelich, L. E. and Jagodziński, A. M. 2018. How much does climate change threaten European forest tree species distributions? – Global Change Biol. 24: 1150–1163. van Eldijk, T. J. B., Bisschop, K. and Etienne, R. S. 2020. Uniting community ecology and evolutionary rescue theory: community-wide rescue leads to a rapid loss of rare species. – Front. Ecol. Evol. 8: 552268

Ellenberg, H. and Mueller-Dombois, D. 1965. Tentative physiognomic-ecological classification of plant formations of the earth. – Ber. Geobot. Inst. Rubel. 37: 21–55.

Fady, B. et al. 2016. Evolution-based approach needed for the conservation and silviculture of peripheral forest tree populations. – For. Ecol. Manage. 375: 66–75.

Fisher, R. A. 1919. XV.—The correlation between relatives on the supposition of Mendelian inheritance. – Trans. R. Soc. Edinb. 52: 399–433.

Fririon, V., Davi, H., Oddou-Muratorio, S., Ligot, G. and Lefèvre, F. 2024. Can thinning foster forest genetic adaptation to drought? A demo-genetic modelling approach with disturbance regimes. – Evol. Appl. 17: e70051.

Gingerich, P. D. 1983. Rates of evolution: effects of time and temporal scaling. – Science 222: 159–161.

Gloy, J., Herzschuh, U. and Kruse, S. 2023. Evolutionary adaptation of trees and modelled future larch forest extent in Siberia. – Ecol. Modell. 478: 110278.

Guillemot, J., Francois, C., Hmimina, G., Dufrêne, E., Martin-StPaul, N. K., Soudani, K., Marie, G., Ourcival, J. M. and Delpierre, N. 2017. Environmental control of carbon allocation matters for modelling forest growth. – New Phytol. 214: 180–193.

Hartmann, H., Bastos, A., Das, A. J., Esquivel-Muelbert, A., Hammond, W. M., Martínez-Vilalta, J., McDowell, N. G., Powers, J. S., Pugh, T. A. M., Ruthrof, K. X. and Allen, C. D. 2022. Climate change risks to global forest health: emergence of unexpected events of elevated tree mortality worldwide. – Annu. Rev. Plant Biol. 73: 673–702.

Hisano, M., Searle, E. B. and Chen, H. Y. H. 2018. Biodiversity as a solution to mitigate climate change impacts on the functioning of forest ecosystems. – Biol. Rev. 93: 439–456.

Hlásny, T., Mátyás, C., Seidl, R., Kulla, L., Merganičová, K., Trombik, J., Dobor, L., Barcza, Z. and Konôpka, B. 2014. Climate change increases the drought risk in central European forests: what are the options for adaptation? – For. J. 60: 5–18.

Johansson, J. 2008. Evolutionary responses to environmental changes: how does competition affect adaptation? – Evolution 62: 421–435.

Karger, D. N., Schmatz, D. R., Dettling, G. and Zimmermann, N. E. 2020. High-resolution monthly precipitation and temperature time series from 2006 to 2100. – Sci. Data 7: 248.

Kattge, J. et al. 2020. TRY plant trait database - enhanced coverage and open access. – Global Change Biol. 26: 119–188.

Keenan, T., Maria Serra, J., Lloret, F., Ninyerola, M. and Sabate, S. 2011. Predicting the future of forests in the Mediterranean under climate change, with niche- and process-based models: CO₂ matters! – Global Change Biol. 17: 565–579.

Kohyama, T. S. and Takada, T. 2012. One-sided competition for light promotes coexistence of forest trees that share the same adult height. – J. Ecol. 100: 1501–1511.

Kramer, K., van der Werf, B. and Schelhaas, M.-J. 2015. Bring in the genes: genetic-ecophysiological modeling of the adaptive response of trees to environmental change. With application to the annual cycle. – Front. Plant Sci. 5: 742.

Lefèvre, F., Boivin, T., Bontemps, A., Courbet, F., Davi, H., Durand-Gillmann, M., Fady, B., Gauzere, J., Gidoin, C., Karam, M. J., Lalagüe, H., Oddou-Muratorio, S. and Pichot, C. 2014. Considering evolutionary processes in adaptive forestry. – Ann. For. Sci. 71: 723–739.

Loreau, M. and Hector, A. 2001. Partitioning selection and complementarity in biodiversity experiments. – Nature 412: 72–76.

Mäkelä, A., Landsberg, J., Ek, A. R., Burk, T. E., Ter-Mikaelian, M., Ågren, G. I., Oliver, C. D. and Puttonen, P. 2000. Process-based models for forest ecosystem management: current state of the art and challenges for practical implementation. – Tree Physiol. 20: 289–298.

Mimura, M. and Aitken, S. N. 2007. Increased selfing and decreased effective pollen donor number in peripheral relative to central populations in *Picea sitchensis* (Pinaceae). – Am. J. Bot. 94: 991–998.

Morin, X., Viner, D. and Chuine, I. 2008. Tree species range shifts at a continental scale: new predictive insights from a process-based model. – J. Ecol. 96: 784–794.

Morin, X., Fahse, L., Jactel, H., Scherer-Lorenzen, M., García-Valdés, R. and Bugmann, H. 2018. Long-term response of forest productivity to climate change is mostly driven by change in tree species composition. – Sci. Rep. 8: 5627.

Morin, X., Damestoy, T., Toigo, M., Castagneyrol, B., Jactel, H., de Coligny, F. and Meredieu, C. 2020. Using forest gap models and experimental data to explore long-term effects of tree diversity on the productivity of mixed planted forests. – Ann. For. Sci. 77: 1–19.

Morin, X., Bugmann, H., de Coligny, F., Martin-StPaul, N., Cailleret, M., Limousin, J. M., Ourcival, J. M., Prevosto, B., Simioni, G., Toigo, M., Vennetier, M., Catteau, E. and Guillemot, J. 2021. Beyond forest succession: a gap model to study ecosystem functioning and tree community composition under climate change. – Funct. Ecol. 35: 955–975.

Niinemets, Ü. and Valladares, F. 2006. Tolerance to shade, drought, and waterlogging of temperate Northern Hemisphere trees and shrubs. – Ecol. Monogr. 76: 521–547.

Oddou-Muratorio, S. and Davi, H. 2014. Simulating local adaptation to climate of forest trees with a physio-demo-genetics model. – Evol. Appl. 7: 453–467.

Osmond, M. M. and de Mazancourt, Cd. 2013. How competition affects evolutionary rescue. – Philos. Trans. R. Soc. B 368: 20120085.

Pardos, M. et al. 2021. The greater resilience of mixed forests to drought mainly depends on their composition: analysis along a climate gradient across Europe. – For. Ecol. Manage. 481: 118687.
Peterson, A. T., Soberón, J. and Sánchez-Cordero, V. 1999. Conservation of ecological niches in evolutionary time. – Science 285: 1265–1267.

Petit, R. J. and Hampe, A. 2006. Some evolutionary consequences of being a tree. – Annu. Rev. Ecol. Evol. Syst. 37: 187–214.

Price, T. D. and Kirkpatrick, M. 2009. Evolutionarily stable range limits set by interspecific competition. – Proc. R. Soc. B 276: 1429–1434.

Rameau, J. C., Mansion, D., Dumé, G., Timbal, J., Lecointe, A., Dupont, P. and Keller, R. 1989. Flore forestière française (guide écologique illustré), tome 1: Plaines et collines, tome 2: Montagnes and tome 3: Région méditerranéenne - IDF.

Romero-Mujalli, D., Jeltsch, F. and Tiedemann, R. 2019. Individual-based modeling of eco-evolutionary dynamics: state of the art and future directions. – Reg. Environ. Change 19: 1–12.

Rosas, T., Mencuccini, M., Barba, J., Cochard, H., Saura-Mas, S. and Martínez-Vilalta, J. 2019. Adjustments and coordination of hydraulic, leaf and stem traits along a water availability gradient. – New Phytol. 223: 632–646.

Ruffault, J., Pimont, F., Cochard, H., Dupuy, J. L. and Martin-StPaul, N. 2022. SurEau-Ecos v2.0: a trait-based plant hydraulics model for simulations of plant water status and drought-induced mortality at the ecosystem level. – Geosci. Model Dev. 15: 5593–5626.

Shugart, H. H., Wang, B., Fischer, R., Ma, J., Fang, J., Yan, X., Huth, A. and Armstrong, A. H. 2018. Gap models and their

individual-based relatives in the assessment of the consequences of global change. – Environ. Res. Lett. 13: 033001.

Silva Pedro, M., Rammer, W. and Seidl, R. 2015. Tree species diversity mitigates disturbance impacts on the forest carbon cycle. – Oecologia 177: 619–630.

Thompson, P. L. and Fronhofer, E. A. 2019. The conflict between adaptation and dispersal for maintaining biodiversity in changing environments. – Proc. Natl Acad. Sci. USA 116: 21061–21067.

Toïgo, M., Vallet, P., Perot, T., Bontemps, J. D., Piedallu, C. and Courbaud, B. 2015. Overyielding in mixed forests decreases with site productivity. – J. Ecol. 103: 502–512.

Van der Meersch, V., Armstrong, E., Mouillot, F., Duputié, A., Davi, H., Saltré, F. and Chuine, I. 2025. Paleorecords reveal biological mechanisms crucial for reliable species range shift projections amid rapid climate change. – Ecol. Lett. 28: e70080.

Wessely, J., Essl, F., Fiedler, K., Gattringer, A., Hülber, B., Ignateva, O., Moser, D., Rammer, W., Dullinger, S. and Seidl, R. 2024. A climate-induced tree species bottleneck for forest management in Europe. – Nat. Ecol. Evol. 8: 1109–1117.

Wiens, J. J., Ackerly, D. D., Allen, A. P., Anacker, B. L., Buckley, L. B., Cornell, H. V., Damschen, E. I., Jonathan Davies, T., Grytnes, J. A., Harrison, S. P., Hawkins, B. A., Holt, R. D., McCain, C. M. and Stephens, P. R. 2010. Niche conservatism as an emerging principle in ecology and conservation biology. – Ecol. Lett. 13: 1310–1324.